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.");
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);
139 int lda = n,info = 0;
148 Lapack::Dgeev(uplo,lrev,lda, A.get(),lda,
153 &work[0],lwork,info);
160 Lapack::Dgeev(uplo,lrev,lda,
165 &work[0],lwork,info);
167 ASSERTL0(info == 0,
"Info is not zero");
174 static unsigned int GetRequiredStorageSize(
unsigned int rows,
unsigned int columns);
179 static unsigned int CalculateIndex(
unsigned int curRow,
unsigned int curColumn);
181 static boost::tuples::tuple<unsigned int, unsigned int>
182 Advance(
const unsigned int totalRows,
const unsigned int totalColumns,
183 const unsigned int curRow,
const unsigned int curColumn);
189 static unsigned int CalculateIndex(
unsigned int totalColumns,
unsigned int curRow,
unsigned int curColumn);
191 static boost::tuples::tuple<unsigned int, unsigned int>
192 Advance(
const unsigned int totalRows,
const unsigned int totalColumns,
193 const unsigned int curRow,
const unsigned int curColumn,
194 char transpose =
'N');
203 static unsigned int CalculateIndex(
unsigned int curRow,
unsigned int curColumn);
205 template<
typename DataType>
206 static void Invert(
unsigned int rows,
unsigned int columns,
209 #ifdef NEKTAR_USING_BLAS
210 ASSERTL0(rows==columns,
"Only square matrices can be inverted.");
217 Lapack::Dsptrf(
'U', n, data.
get(), ipivot.get(), info);
221 std::string message =
"ERROR: The " + boost::lexical_cast<std::string>(-info) +
"th parameter had an illegal parameter for dsptrf";
226 std::string message =
"ERROR: Element u_" + boost::lexical_cast<std::string>(info) + boost::lexical_cast<std::string>(info) +
" is 0 from dsptrf";
230 Lapack::Dsptri(
'U', n, data.
get(), ipivot.get(),
235 std::string message =
"ERROR: The " + boost::lexical_cast<std::string>(-info) +
"th parameter had an illegal parameter for dsptri";
240 std::string message =
"ERROR: Element u_" + boost::lexical_cast<std::string>(info) + boost::lexical_cast<std::string>(info) +
" is 0 from dsptri";
246 BOOST_STATIC_ASSERT(
sizeof(DataType) == 0);
250 static boost::tuples::tuple<unsigned int, unsigned int>
251 Advance(
const unsigned int totalRows,
const unsigned int totalColumns,
252 const unsigned int curRow,
const unsigned int curColumn);
257 static boost::tuples::tuple<unsigned int, unsigned int>
258 Advance(
const unsigned int totalRows,
const unsigned int totalColumns,
259 const unsigned int curRow,
const unsigned int curColumn);
261 template<
typename DataType>
262 static void Invert(
unsigned int rows,
unsigned int columns,
265 ASSERTL0(rows==columns,
"Only square matrices can be inverted.");
266 for(
unsigned int i = 0; i < rows; ++i)
268 data[i] = 1.0/data[i];
272 static unsigned int GetRequiredStorageSize(
unsigned int rows,
unsigned int columns);
274 static unsigned int CalculateIndex(
unsigned int row,
unsigned int col);
280 static unsigned int GetRequiredStorageSize(
unsigned int rows,
unsigned int columns,
281 unsigned int nSubSuperDiags);
298 static unsigned int CalculateIndex(
unsigned int curRow,
unsigned int curColumn,
299 unsigned int nSuperDiags);
304 #endif //NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_MATRIX_FUNCS_H
#define ASSERTL0(condition, msg)
static Array< OneD, NekDouble > NullNekDouble1DArray
static void Invert(unsigned int rows, unsigned int columns, Array< OneD, DataType > &data, const char transpose)
static unsigned int GetRequiredStorageSize(unsigned int rows, unsigned int columns)
#define LIB_UTILITIES_EXPORT
static unsigned int GetRequiredStorageSize(unsigned int rows, unsigned int columns, unsigned int nSubSuperDiags)
static void EigenSolve(unsigned int n, const Array< OneD, const double > &A, Array< OneD, NekDouble > &EigValReal, Array< OneD, NekDouble > &EigValImag, Array< OneD, NekDouble > &EigVecs=NullNekDouble1DArray)
static void Invert(unsigned int rows, unsigned int columns, Array< OneD, DataType > &data)
static void Invert(unsigned int rows, unsigned int columns, Array< OneD, DataType > &data)