Nektar++
MatrixFuncs.h
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: MatrixFuncs.h
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: Matrix functions that depend on storage policy. Putting
32// methods in these separate classes makes it easier to use them from
33// normal, scaled, and block matrices.
34//
35///////////////////////////////////////////////////////////////////////////////
36
37#ifndef NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_MATRIX_FUNCS_H
38#define NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_MATRIX_FUNCS_H
39
43#include <limits>
44#include <tuple>
45
46namespace Nektar
47{
49{
50 /// \brief Calculates and returns the storage size required.
51 ///
52 /// This method assumes that the matrix will be used with LU factorizationa
53 /// and allocates additional storage as appropriate.
54 static unsigned int GetRequiredStorageSize(unsigned int totalRows,
55 unsigned int totalColumns,
56 unsigned int subDiags,
57 unsigned int superDiags);
58
59 static unsigned int CalculateNumberOfDiags(unsigned int totalRows,
60 unsigned int diags);
61
62 static unsigned int CalculateNumberOfRows(unsigned int totalRows,
63 unsigned int subDiags,
64 unsigned int superDiags);
65
66 static unsigned int CalculateIndex(unsigned int totalRows,
67 unsigned int totalColumns,
68 unsigned int row, unsigned int column,
69 unsigned int sub, unsigned int super);
70
71 static std::tuple<unsigned int, unsigned int> Advance(
72 const unsigned int totalRows, const unsigned int totalColumns,
73 const unsigned int curRow, const unsigned int curColumn);
74};
75
77{
78
79 static unsigned int GetRequiredStorageSize(unsigned int rows,
80 unsigned int columns);
81 static unsigned int CalculateIndex(unsigned int totalRows,
82 unsigned int totalColumns,
83 unsigned int curRow,
84 unsigned int curColumn);
85
86 static std::tuple<unsigned int, unsigned int> Advance(
87 const unsigned int totalRows, const unsigned int totalColumns,
88 const unsigned int curRow, const unsigned int curColumn);
89
90 template <typename DataType>
91 static void Invert(unsigned int rows, unsigned int columns,
92 Array<OneD, DataType> &data, const char transpose)
93 {
94 ASSERTL0(rows == columns, "Only square matrices can be inverted.");
95 ASSERTL0(transpose == 'N',
96 "Only untransposed matrices may be inverted.");
97
98 int m = rows;
99 int n = columns;
100 int info = 0;
101 Array<OneD, int> ipivot(n);
102 Array<OneD, DataType> work(n);
103
104 if (std::is_floating_point_v<DataType>)
105 {
106 switch (sizeof(DataType))
107 {
108 case sizeof(NekDouble):
109 break;
110 case sizeof(NekSingle):
111 break;
112 default:
113 ASSERTL0(
114 false,
115 "Invert DataType is neither NekDouble nor NekSingle");
116 break;
117 }
118 }
119 else
120 {
121 ASSERTL0(false,
122 "FullMatrixFuncs::Invert DataType is not floating point");
123 }
124
125 Lapack::DoSgetrf(m, n, data.data(), m, ipivot.data(), info);
126
127 if (info < 0)
128 {
129 std::string message =
130 "ERROR: The " + std::to_string(-info) +
131 "th parameter had an illegal parameter for dgetrf";
132 ASSERTL0(false, message.c_str());
133 }
134 else if (info > 0)
135 {
136 std::string message = "ERROR: Element u_" + std::to_string(info) +
137 std::to_string(info) + " is 0 from dgetrf";
138 ASSERTL0(false, message.c_str());
139 }
140
141 Lapack::DoSgetri(n, data.data(), n, ipivot.data(), work.data(), n,
142 info);
143
144 if (info < 0)
145 {
146 std::string message =
147 "ERROR: The " + std::to_string(-info) +
148 "th parameter had an illegal parameter for dgetri";
149 ASSERTL0(false, message.c_str());
150 }
151 else if (info > 0)
152 {
153 std::string message = "ERROR: Element u_" + std::to_string(info) +
154 std::to_string(info) + " is 0 from dgetri";
155 ASSERTL0(false, message.c_str());
156 }
157 }
158
159 template <typename DataType>
160 static void EigenSolve(unsigned int n, const Array<OneD, const DataType> &A,
161 Array<OneD, DataType> &EigValReal,
162 Array<OneD, DataType> &EigValImag,
163 Array<OneD, DataType> &EigVecs)
164 {
165 int lda = n, info = 0;
166 DataType dum;
167 char uplo = 'N';
168
169 if (EigVecs.size() != 0) // calculate Right Eigen Vectors
170 {
171 int lwork = 4 * lda;
172 Array<OneD, DataType> work(4 * lda);
173 char lrev = 'V';
174 Lapack::DoSgeev(uplo, lrev, lda, A.data(), lda, EigValReal.data(),
175 EigValImag.data(), &dum, 1, EigVecs.data(), lda,
176 &work[0], lwork, info);
177 }
178 else
179 {
180 int lwork = 3 * lda;
181 Array<OneD, DataType> work(3 * lda);
182 char lrev = 'N';
183 Lapack::DoSgeev(uplo, lrev, lda, A.data(), lda, EigValReal.data(),
184 EigValImag.data(), &dum, 1, &dum, 1, &work[0],
185 lwork, info);
186 }
187 ASSERTL0(info == 0, "Info is not zero");
188 }
189};
190
192{
193 static unsigned int GetRequiredStorageSize(unsigned int rows,
194 unsigned int columns);
195};
196
198 : public TriangularMatrixFuncs
199{
200 static unsigned int CalculateIndex(unsigned int curRow,
201 unsigned int curColumn);
202
203 static std::tuple<unsigned int, unsigned int> Advance(
204 const unsigned int totalRows, const unsigned int totalColumns,
205 const unsigned int curRow, const unsigned int curColumn);
206};
207
209 : public TriangularMatrixFuncs
210{
211 static unsigned int CalculateIndex(unsigned int totalColumns,
212 unsigned int curRow,
213 unsigned int curColumn);
214
215 static std::tuple<unsigned int, unsigned int> Advance(
216 const unsigned int totalRows, const unsigned int totalColumns,
217 const unsigned int curRow, const unsigned int curColumn,
218 char transpose = 'N');
219};
220
221/// \internal
222/// Symmetric matrices use upper triangular packed storage.
224{
226
227 static unsigned int CalculateIndex(unsigned int curRow,
228 unsigned int curColumn);
229
230 template <typename DataType>
231 static void Invert(unsigned int rows, unsigned int columns,
233 {
234 ASSERTL0(rows == columns, "Only square matrices can be inverted.");
235
236 int n = columns;
237 int info = 0;
238 Array<OneD, int> ipivot(n);
239 Array<OneD, DataType> work(n);
240
241 Lapack::DoSsptrf('U', n, data.data(), ipivot.data(), info);
242
243 if (info < 0)
244 {
245 std::string message =
246 "ERROR: The " + std::to_string(-info) +
247 "th parameter had an illegal parameter for dsptrf";
248 ASSERTL0(false, message.c_str());
249 }
250 else if (info > 0)
251 {
252 std::string message = "ERROR: Element u_" + std::to_string(info) +
253 std::to_string(info) + " is 0 from dsptrf";
254 ASSERTL0(false, message.c_str());
255 }
256
257 Lapack::DoSsptri('U', n, data.data(), ipivot.data(), work.data(), info);
258
259 if (info < 0)
260 {
261 std::string message =
262 "ERROR: The " + std::to_string(-info) +
263 "th parameter had an illegal parameter for dsptri";
264 ASSERTL0(false, message.c_str());
265 }
266 else if (info > 0)
267 {
268 std::string message = "ERROR: Element u_" + std::to_string(info) +
269 std::to_string(info) + " is 0 from dsptri";
270 ASSERTL0(false, message.c_str());
271 }
272 }
273
274 static std::tuple<unsigned int, unsigned int> Advance(
275 const unsigned int totalRows, const unsigned int totalColumns,
276 const unsigned int curRow, const unsigned int curColumn);
277};
278
280{
281 static std::tuple<unsigned int, unsigned int> Advance(
282 const unsigned int totalRows, const unsigned int totalColumns,
283 const unsigned int curRow, const unsigned int curColumn);
284
285 template <typename DataType>
286 static void Invert(unsigned int rows, unsigned int columns,
288 {
289 ASSERTL0(rows == columns, "Only square matrices can be inverted.");
290 for (unsigned int i = 0; i < rows; ++i)
291 {
292 data[i] = 1.0 / data[i];
293 }
294 }
295
296 static unsigned int GetRequiredStorageSize(unsigned int rows,
297 unsigned int columns);
298
299 static unsigned int CalculateIndex(unsigned int row, unsigned int col);
300};
301
303{
304 static unsigned int GetRequiredStorageSize(unsigned int rows,
305 unsigned int columns,
306 unsigned int nSubSuperDiags);
307};
308
311{
312};
313
316{
317};
318
319/// \internal
320/// Symmetric banded matrices use upper triangular banded packed storage.
323{
325
326 static unsigned int CalculateIndex(unsigned int curRow,
327 unsigned int curColumn,
328 unsigned int nSuperDiags);
329};
330
331} // namespace Nektar
332
333#endif // NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_MATRIX_FUNCS_H
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
#define LIB_UTILITIES_EXPORT
1D Array of constant elements with garbage collection and bounds checking.
size_type size() const
Returns the array's size.
static void DoSsptrf(const char &uplo, const int &n, double *ap, int *ipiv, int &info)
factor a real packed-symmetric matrix using Bunch-Kaufman pivoting.
Definition: Lapack.hpp:112
static void DoSsptri(const char &uplo, const int &n, const double *ap, const int *ipiv, double *work, int &info)
Invert a real packed-symmetric matrix problem.
Definition: Lapack.hpp:162
static void DoSgetrf(const int &m, const int &n, double *a, const int &lda, int *ipiv, int &info)
General matrix LU factorisation.
Definition: Lapack.hpp:248
static void DoSgetri(const int &n, double *a, const int &lda, const int *ipiv, double *wk, const int &lwk, int &info)
General matrix inverse.
Definition: Lapack.hpp:277
static void DoSgeev(const char &uplo, const char &lrev, const int &n, const double *a, const int &lda, double *wr, double *wi, double *rev, const int &ldr, double *lev, const int &ldv, double *work, const int &lwork, int &info)
Solve general real matrix eigenproblem.
Definition: Lapack.hpp:307
double NekDouble
static void Invert(unsigned int rows, unsigned int columns, Array< OneD, DataType > &data)
Definition: MatrixFuncs.h:286
static void EigenSolve(unsigned int n, const Array< OneD, const DataType > &A, Array< OneD, DataType > &EigValReal, Array< OneD, DataType > &EigValImag, Array< OneD, DataType > &EigVecs)
Definition: MatrixFuncs.h:160
static void Invert(unsigned int rows, unsigned int columns, Array< OneD, DataType > &data, const char transpose)
Definition: MatrixFuncs.h:91
static void Invert(unsigned int rows, unsigned int columns, Array< OneD, DataType > &data)
Definition: MatrixFuncs.h:231
static unsigned int GetRequiredStorageSize(unsigned int rows, unsigned int columns, unsigned int nSubSuperDiags)
static unsigned int GetRequiredStorageSize(unsigned int rows, unsigned int columns)