38 #ifndef NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_MATRIX_FUNCS_H
39 #define NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_MATRIX_FUNCS_H
41 #include <boost/lexical_cast.hpp>
42 #include <boost/tuple/tuple.hpp>
56 static unsigned int GetRequiredStorageSize(
unsigned int totalRows,
unsigned int totalColumns,
57 unsigned int subDiags,
unsigned int superDiags);
59 static unsigned int CalculateNumberOfDiags(
unsigned int totalRows,
unsigned int diags);
61 static unsigned int CalculateNumberOfRows(
unsigned int totalRows,
unsigned int subDiags,
unsigned int superDiags);
63 static unsigned int CalculateIndex(
unsigned int totalRows,
64 unsigned int totalColumns,
65 unsigned int row,
unsigned int column,
66 unsigned int sub,
unsigned int super);
69 static boost::tuples::tuple<unsigned int, unsigned int>
70 Advance(
const unsigned int totalRows,
const unsigned int totalColumns,
71 const unsigned int curRow,
const unsigned int curColumn);
77 static unsigned int GetRequiredStorageSize(
unsigned int rows,
unsigned int columns);
78 static unsigned int CalculateIndex(
unsigned int totalRows,
unsigned int totalColumns,
unsigned int curRow,
unsigned int curColumn);
81 static boost::tuples::tuple<unsigned int, unsigned int>
82 Advance(
const unsigned int totalRows,
const unsigned int totalColumns,
83 const unsigned int curRow,
const unsigned int curColumn);
85 template<
typename DataType>
86 static void Invert(
unsigned int rows,
unsigned int columns,
90 #ifdef NEKTAR_USING_BLAS
91 ASSERTL0(rows==columns,
"Only square matrices can be inverted.");
92 ASSERTL0(transpose==
'N',
"Only untransposed matrices may be inverted.");
97 Array<OneD, int> ipivot(n);
100 Lapack::Dgetrf(m, n, data.
get(), m, ipivot.get(), info);
104 std::string message =
"ERROR: The " + boost::lexical_cast<std::string>(-info) +
"th parameter had an illegal parameter for dgetrf";
109 std::string message =
"ERROR: Element u_" + boost::lexical_cast<std::string>(info) + boost::lexical_cast<std::string>(info) +
" is 0 from dgetrf";
113 Lapack::Dgetri(n, data.
get(), n, ipivot.get(),
114 work.
get(), n, info);
118 std::string message =
"ERROR: The " + boost::lexical_cast<std::string>(-info) +
"th parameter had an illegal parameter for dgetri";
123 std::string message =
"ERROR: Element u_" + boost::lexical_cast<std::string>(info) + boost::lexical_cast<std::string>(info) +
" is 0 from dgetri";
129 BOOST_STATIC_ASSERT(
sizeof(DataType) == 0);
133 static void EigenSolve(
unsigned int n,
134 const Array<OneD, const double>& A,
135 Array<OneD, NekDouble> &EigValReal,
136 Array<OneD, NekDouble> &EigValImag,
139 int lda = n,info = 0;
142 Array<OneD,NekDouble> work(3*lda);
148 Lapack::Dgeev(uplo,lrev,lda, A.get(),lda,
153 &work[0],lwork,info);
158 Lapack::Dgeev(uplo,lrev,lda,
163 &work[0],lwork,info);
165 ASSERTL0(info == 0,
"Info is not zero");
172 static unsigned int GetRequiredStorageSize(
unsigned int rows,
unsigned int columns);
177 static unsigned int CalculateIndex(
unsigned int curRow,
unsigned int curColumn);
179 static boost::tuples::tuple<unsigned int, unsigned int>
180 Advance(
const unsigned int totalRows,
const unsigned int totalColumns,
181 const unsigned int curRow,
const unsigned int curColumn);
187 static unsigned int CalculateIndex(
unsigned int totalColumns,
unsigned int curRow,
unsigned int curColumn);
189 static boost::tuples::tuple<unsigned int, unsigned int>
190 Advance(
const unsigned int totalRows,
const unsigned int totalColumns,
191 const unsigned int curRow,
const unsigned int curColumn,
192 char transpose =
'N');
201 static unsigned int CalculateIndex(
unsigned int curRow,
unsigned int curColumn);
203 static boost::tuples::tuple<unsigned int, unsigned int>
204 Advance(
const unsigned int totalRows,
const unsigned int totalColumns,
205 const unsigned int curRow,
const unsigned int curColumn);
210 static boost::tuples::tuple<unsigned int, unsigned int>
211 Advance(
const unsigned int totalRows,
const unsigned int totalColumns,
212 const unsigned int curRow,
const unsigned int curColumn);
214 template<
typename DataType>
215 static void Invert(
unsigned int rows,
unsigned int columns,
218 ASSERTL0(rows==columns,
"Only square matrices can be inverted.");
219 for(
unsigned int i = 0; i < rows; ++i)
221 data[i] = 1.0/data[i];
225 static unsigned int GetRequiredStorageSize(
unsigned int rows,
unsigned int columns);
227 static unsigned int CalculateIndex(
unsigned int row,
unsigned int col);
233 static unsigned int GetRequiredStorageSize(
unsigned int rows,
unsigned int columns,
234 unsigned int nSubSuperDiags);
251 static unsigned int CalculateIndex(
unsigned int curRow,
unsigned int curColumn,
252 unsigned int nSuperDiags);
257 #endif //NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_MATRIX_FUNCS_H