36 #ifndef NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_NEK_LINSYS_HPP
37 #define NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_NEK_LINSYS_HPP
49 #include <boost/shared_ptr.hpp>
50 #include <boost/utility/enable_if.hpp>
51 #include <boost/type_traits.hpp>
63 template<
typename DataType>
66 template<
typename DataType>
67 struct IsSharedPointer<boost::shared_ptr<DataType> > :
public boost::true_type {};
74 template<
typename BVectorType,
typename XVectorType>
78 char m_transposeFlag,
unsigned int m_numberOfSubDiagonals,
79 unsigned int m_numberOfSuperDiagonals)
87 Lapack::Dgetrs(
'N',n,1,A.get(),n,(
int *)m_ipivot.get(),x.GetRawPtr(),n,info);
90 std::string message =
"ERROR: The " + boost::lexical_cast<std::string>(-info) +
"th parameter had an illegal parameter for dgetrs";
97 for(
unsigned int i = 0; i < A.num_elements(); ++i)
106 Lapack::Dtptrs(
'U', m_transposeFlag,
'N', n, 1, A.get(), x.GetRawPtr(), n, info);
110 std::string message =
"ERROR: The " + boost::lexical_cast<std::string>(-info) +
"th parameter had an illegal parameter for dtrtrs";
115 std::string message =
"ERROR: The " + boost::lexical_cast<std::string>(-info) +
"th diagonal element of A is 0 for dtrtrs";
124 Lapack::Dtptrs(
'L', m_transposeFlag,
'N', n, 1, A.get(), x.GetRawPtr(), n, info);
128 std::string message =
"ERROR: The " + boost::lexical_cast<std::string>(-info) +
"th parameter had an illegal parameter for dtrtrs";
133 std::string message =
"ERROR: The " + boost::lexical_cast<std::string>(-info) +
"th diagonal element of A is 0 for dtrtrs";
142 Lapack::Dsptrs(
'U', n, 1, A.get(), m_ipivot.get(), x.GetRawPtr(), x.GetRows(), info);
145 std::string message =
"ERROR: The " + boost::lexical_cast<std::string>(-info) +
"th parameter had an illegal parameter for dsptrs";
154 Lapack::Dpptrs(
'U', n, 1, A.get(), x.GetRawPtr(), x.GetRows(), info);
157 std::string message =
"ERROR: The " + boost::lexical_cast<std::string>(-info) +
"th parameter had an illegal parameter for dpptrs";
165 int KL = m_numberOfSubDiagonals;
166 int KU = m_numberOfSuperDiagonals;
169 Lapack::Dgbtrs(m_transposeFlag, n, KL, KU, 1, A.get(), 2*KL+KU+1, m_ipivot.get(), x.GetRawPtr(), n, info);
173 std::string message =
"ERROR: The " + boost::lexical_cast<std::string>(-info) +
"th parameter had an illegal parameter for dgbtrs";
181 int KU = m_numberOfSuperDiagonals;
184 Lapack::Dpbtrs(
'U', n, KU, 1, A.get(), KU+1, x.GetRawPtr(), n, info);
188 std::string message =
"ERROR: The " + boost::lexical_cast<std::string>(-info) +
"th parameter had an illegal parameter for dpbtrs";
209 template<
typename BVectorType,
typename XVectorType>
213 char m_transposeFlag,
unsigned int m_numberOfSubDiagonals,
214 unsigned int m_numberOfSuperDiagonals)
222 Lapack::Dgetrs(
'T',n,1,A.get(),n,(
int *)m_ipivot.get(),x.GetRawPtr(), n,info);
226 std::string message =
"ERROR: The " + boost::lexical_cast<std::string>(-info) +
"th parameter had an illegal parameter for dgetrs";
233 Solve(b, x, m_matrixType, m_ipivot, n, A, m_transposeFlag, m_numberOfSubDiagonals, m_numberOfSuperDiagonals);
237 char trans = m_transposeFlag;
249 Lapack::Dtptrs(
'U', trans,
'N', n, 1, A.get(), x.GetRawPtr(), n, info);
253 std::string message =
"ERROR: The " + boost::lexical_cast<std::string>(-info) +
"th parameter had an illegal parameter for dtrtrs";
258 std::string message =
"ERROR: The " + boost::lexical_cast<std::string>(-info) +
"th diagonal element of A is 0 for dtrtrs";
266 char trans = m_transposeFlag;
277 Lapack::Dtptrs(
'L', trans,
'N', n, 1, A.get(), x.GetRawPtr(), n, info);
281 std::string message =
"ERROR: The " + boost::lexical_cast<std::string>(-info) +
"th parameter had an illegal parameter for dtrtrs";
286 std::string message =
"ERROR: The " + boost::lexical_cast<std::string>(-info) +
"th diagonal element of A is 0 for dtrtrs";
294 Solve(b, x, m_matrixType, m_ipivot, n, A, m_transposeFlag, m_numberOfSubDiagonals, m_numberOfSuperDiagonals);
299 int KL = m_numberOfSubDiagonals;
300 int KU = m_numberOfSuperDiagonals;
303 Lapack::Dgbtrs(m_transposeFlag, n, KL, KU, 1, A.get(), 2*KL+KU+1, m_ipivot.get(), x.GetRawPtr(), n, info);
307 std::string message =
"ERROR: The " + boost::lexical_cast<std::string>(-info) +
"th parameter had an illegal parameter for dgbtrs";
332 template<
typename MatrixType>
344 ASSERTL0(theA->GetTransposeFlag() ==
'N',
"LinearSystem requires a non-transposed matrix.");
347 if( wrapperType ==
eCopy )
356 template<
typename MatrixType>
368 ASSERTL0(theA.GetTransposeFlag() ==
'N',
"LinearSystem requires a non-transposed matrix.");
371 if( wrapperType ==
eCopy )
402 template<
typename VectorType>
411 template<
typename BType,
typename XType>
412 void Solve(
const BType& b, XType& x)
const
420 template<
typename VectorType>
429 template<
typename BType,
typename XType>
441 template<
typename MatrixType>
448 int m = theA.GetRows();
449 int n = theA.GetColumns();
451 int pivotSize = std::max(1, std::min(m, n));
455 Lapack::Dgetrf(m, n,
A.get(), m,
m_ipivot.get(), info);
459 std::string message =
"ERROR: The " + boost::lexical_cast<std::string>(-info) +
"th parameter had an illegal parameter for dgetrf";
464 std::string message =
"ERROR: Element u_" + boost::lexical_cast<std::string>(info) + boost::lexical_cast<std::string>(info) +
" is 0 from dgetrf";
470 for(
unsigned int i = 0; i < theA.GetColumns(); ++i)
472 A[i] = 1.0/theA(i,i);
481 int pivotSize = theA.GetRows();
484 Lapack::Dsptrf(
'U', theA.GetRows(),
A.get(),
m_ipivot.get(), info);
488 std::string message =
"ERROR: The " + boost::lexical_cast<std::string>(-info) +
"th parameter had an illegal parameter for dsptrf";
493 std::string message =
"ERROR: Element u_" + boost::lexical_cast<std::string>(info) + boost::lexical_cast<std::string>(info) +
" is 0 from dsptrf";
501 Lapack::Dpptrf(
'U', theA.GetRows(),
A.get(), info);
505 std::string message =
"ERROR: The " + boost::lexical_cast<std::string>(-info) +
"th parameter had an illegal parameter for dpptrf";
510 std::string message =
"ERROR: The leading minor of order " + boost::lexical_cast<std::string>(info) +
" is not positive definite from dpptrf";
528 unsigned int rawRows = KL+KU+1;
532 for(
unsigned int i = 0; i < theA.GetColumns(); ++i)
534 std::copy(theA.GetRawPtr() + i*rawRows, theA.GetRawPtr() + (i+1)*rawRows,
535 A.get() + (i+1)*KL + i*rawRows);
539 int pivotSize = theA.GetRows();
542 Lapack::Dgbtrf(M, N, KL, KU,
A.get(), 2*KL+KU+1,
m_ipivot.get(), info);
546 std::string message =
"ERROR: The " + boost::lexical_cast<std::string>(-info) +
"th parameter had an illegal parameter for dgbtrf";
551 std::string message =
"ERROR: Element u_" + boost::lexical_cast<std::string>(info) + boost::lexical_cast<std::string>(info) +
" is 0 from dgbtrf";
559 std::string(
"Number of sub- and superdiagonals should ") +
560 std::string(
"be equal for a symmetric banded matrix"));
564 Lapack::Dpbtrf(
'U', theA.GetRows(), KU,
A.get(), KU+1, info);
568 std::string message =
"ERROR: The " + boost::lexical_cast<std::string>(-info) +
"th parameter had an illegal parameter for dpbtrf";
573 std::string message =
"ERROR: The leading minor of order " + boost::lexical_cast<std::string>(info) +
" is not positive definite from dpbtrf";
614 #endif //NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_NEK_LINSYS_HPP
#define ASSERTL0(condition, msg)
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
MatrixStorage m_matrixType
void swap(LinearSystem &rhs)
LinearSystem(const boost::shared_ptr< MatrixType > &theA, PointerWrapper wrapperType=eCopy)
void Solve(const BType &b, XType &x) const
unsigned int GetColumns() const
LinearSystem & operator=(const LinearSystem &rhs)
void SolveTranspose(const BType &b, XType &x) const
unsigned int m_numberOfSuperDiagonals
RawType< VectorType >::type Solve(const VectorType &b)
unsigned int GetRows() const
void FactorMatrix(const MatrixType &theA)
void CopyArray(const Array< OneD, ConstDataType > &source, Array< OneD, DataType > &dest)
PointerWrapper
Specifies if the pointer passed to a NekMatrix or NekVector is copied into an internal representation...
LinearSystem(const LinearSystem &rhs)
LinearSystem(const MatrixType &theA, PointerWrapper wrapperType=eCopy)
static void SolveTranspose(const BVectorType &b, XVectorType &x, MatrixStorage m_matrixType, const Array< OneD, const int > &m_ipivot, unsigned int n, const Array< OneD, const double > &A, char m_transposeFlag, unsigned int m_numberOfSubDiagonals, unsigned int m_numberOfSuperDiagonals)
static void Solve(const BVectorType &b, XVectorType &x, MatrixStorage m_matrixType, const Array< OneD, const int > &m_ipivot, unsigned int n, const Array< OneD, const double > &A, char m_transposeFlag, unsigned int m_numberOfSubDiagonals, unsigned int m_numberOfSuperDiagonals)
unsigned int m_numberOfSubDiagonals
RawType< VectorType >::type SolveTranspose(const VectorType &b)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
static unsigned int GetRequiredStorageSize(unsigned int totalRows, unsigned int totalColumns, unsigned int subDiags, unsigned int superDiags)
Calculates and returns the storage size required.
Array< OneD, int > m_ipivot