Nektar++
Loading...
Searching...
No Matches
Static Public Member Functions | List of all members
Nektar::FullMatrixFuncs Struct Reference

#include <MatrixFuncs.h>

Static Public Member Functions

static unsigned int GetRequiredStorageSize (unsigned int rows, unsigned int columns)
 
static unsigned int CalculateIndex (unsigned int totalRows, unsigned int totalColumns, unsigned int curRow, unsigned int curColumn)
 
static std::tuple< unsigned int, unsigned int > Advance (const unsigned int totalRows, const unsigned int totalColumns, const unsigned int curRow, const unsigned int curColumn)
 
template<typename DataType >
static void Invert (unsigned int rows, unsigned int columns, Array< OneD, DataType > &data, const char transpose)
 
template<typename DataType >
static void PseudoInverse (unsigned int rows, unsigned int columns, Array< OneD, DataType > &data)
 
template<typename DataType >
static void EigenSolve (unsigned int n, const Array< OneD, const DataType > &A, Array< OneD, DataType > &EigValReal, Array< OneD, DataType > &EigValImag, Array< OneD, DataType > &EigVecs)
 

Detailed Description

Definition at line 80 of file MatrixFuncs.h.

Member Function Documentation

◆ Advance()

std::tuple< unsigned int, unsigned int > Nektar::FullMatrixFuncs::Advance ( const unsigned int  totalRows,
const unsigned int  totalColumns,
const unsigned int  curRow,
const unsigned int  curColumn 
)
static

Definition at line 121 of file MatrixFuncs.cpp.

124{
125 unsigned int nextRow = curRow;
126 unsigned int nextColumn = curColumn;
127
128 if (nextRow < totalRows)
129 {
130 ++nextRow;
131 }
132
133 if (nextRow >= totalRows)
134 {
135 nextRow = 0;
136 ++nextColumn;
137 }
138
139 if (nextColumn >= totalColumns)
140 {
141 nextRow = std::numeric_limits<unsigned int>::max();
142 nextColumn = std::numeric_limits<unsigned int>::max();
143 }
144
145 return std::tuple<unsigned int, unsigned int>(nextRow, nextColumn);
146}

Referenced by Nektar::NekMatrix< DataType, StandardMatrixTag >::Advance(), and Nektar::FullMatrixStoragePolicyUnitTests::BOOST_AUTO_TEST_CASE().

◆ CalculateIndex()

unsigned int Nektar::FullMatrixFuncs::CalculateIndex ( unsigned int  totalRows,
unsigned int  totalColumns,
unsigned int  curRow,
unsigned int  curColumn 
)
static

Definition at line 114 of file MatrixFuncs.cpp.

117{
118 return curColumn * totalRows + curRow;
119}

Referenced by Nektar::FullMatrixStoragePolicyUnitTests::BOOST_AUTO_TEST_CASE().

◆ EigenSolve()

template<typename DataType >
static void Nektar::FullMatrixFuncs::EigenSolve ( unsigned int  n,
const Array< OneD, const DataType > &  A,
Array< OneD, DataType > &  EigValReal,
Array< OneD, DataType > &  EigValImag,
Array< OneD, DataType > &  EigVecs 
)
inlinestatic

Definition at line 298 of file MatrixFuncs.h.

302 {
303 int lda = n, info = 0;
304 DataType dum;
305 char uplo = 'N';
306
307 if (EigVecs.size() != 0) // calculate Right Eigen Vectors
308 {
309 int lwork = 4 * lda;
310 Array<OneD, DataType> work(4 * lda);
311 char lrev = 'V';
312 Lapack::DoSgeev(uplo, lrev, lda, A.data(), lda, EigValReal.data(),
313 EigValImag.data(), &dum, 1, EigVecs.data(), lda,
314 &work[0], lwork, info);
315 }
316 else
317 {
318 int lwork = 3 * lda;
319 Array<OneD, DataType> work(3 * lda);
320 char lrev = 'N';
321 Lapack::DoSgeev(uplo, lrev, lda, A.data(), lda, EigValReal.data(),
322 EigValImag.data(), &dum, 1, &dum, 1, &work[0],
323 lwork, info);
324 }
325 ASSERTL0(info == 0, "Info is not zero");
326 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition ErrorUtil.hpp:216
Lapack::DoSgeev
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:335

References ASSERTL0, and Lapack::DoSgeev().

Referenced by Nektar::NekMatrix< DataType, StandardMatrixTag >::EigenSolve().

◆ GetRequiredStorageSize()

unsigned int Nektar::FullMatrixFuncs::GetRequiredStorageSize ( unsigned int  rows,
unsigned int  columns 
)
static

Definition at line 108 of file MatrixFuncs.cpp.

110{
111 return rows * columns;
112}

Referenced by Nektar::ConstMatrix< DataType >::GetRequiredStorageSize().

◆ Invert()

template<typename DataType >
static void Nektar::FullMatrixFuncs::Invert ( unsigned int  rows,
unsigned int  columns,
Array< OneD, DataType > &  data,
const char  transpose 
)
inlinestatic

Definition at line 95 of file MatrixFuncs.h.

97 {
98 ASSERTL0(rows == columns, "Only square matrices can be inverted.");
99 ASSERTL0(transpose == 'N',
100 "Only untransposed matrices may be inverted.");
101
102 int m = rows;
103 int n = columns;
104 int info = 0;
105 Array<OneD, int> ipivot(n);
106 Array<OneD, DataType> work(n);
107
108 if (std::is_floating_point_v<DataType>)
109 {
110 switch (sizeof(DataType))
111 {
112 case sizeof(NekDouble):
113 break;
114 case sizeof(NekSingle):
115 break;
116 default:
117 ASSERTL0(
118 false,
119 "Invert DataType is neither NekDouble nor NekSingle");
120 break;
121 }
122 }
123 else
124 {
125 ASSERTL0(false,
126 "FullMatrixFuncs::Invert DataType is not floating point");
127 }
128
129 Lapack::DoSgetrf(m, n, data.data(), m, ipivot.data(), info);
130
131 if (info < 0)
132 {
133 std::string message =
134 "ERROR: The " + std::to_string(-info) +
135 "th parameter had an illegal parameter for dgetrf";
136 ASSERTL0(false, message.c_str());
137 }
138 else if (info > 0)
139 {
140 std::string message = "ERROR: Element u_" + std::to_string(info) +
141 std::to_string(info) + " is 0 from dgetrf";
142 ASSERTL0(false, message.c_str());
143 }
144
145 Lapack::DoSgetri(n, data.data(), n, ipivot.data(), work.data(), n,
146 info);
147
148 if (info < 0)
149 {
150 std::string message =
151 "ERROR: The " + std::to_string(-info) +
152 "th parameter had an illegal parameter for dgetri";
153 ASSERTL0(false, message.c_str());
154 }
155 else if (info > 0)
156 {
157 std::string message = "ERROR: Element u_" + std::to_string(info) +
158 std::to_string(info) + " is 0 from dgetri";
159 ASSERTL0(false, message.c_str());
160 }
161 }
Lapack::DoSgetrf
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:276
Lapack::DoSgetri
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:305

References ASSERTL0, Lapack::DoSgetrf(), and Lapack::DoSgetri().

Referenced by Nektar::NekMatrix< DataType, StandardMatrixTag >::Invert().

◆ PseudoInverse()

template<typename DataType >
static void Nektar::FullMatrixFuncs::PseudoInverse ( unsigned int  rows,
unsigned int  columns,
Array< OneD, DataType > &  data 
)
inlinestatic

Definition at line 166 of file MatrixFuncs.h.

168 {
169 int m = rows;
170 int n = columns;
171 int lda = m;
172 int ldu = m;
173 int ldvt = n;
174 int info = 0;
175 int min_mn = std::min(m, n);
176 // Copy input data since SVD destroys it
177 Array<OneD, DataType> a(data);
178 // SVD outputs
179 Array<OneD, DataType> s(min_mn, 0.0);
180 Array<OneD, DataType> u(ldu * m, 0.0);
181 Array<OneD, DataType> vt(ldvt * n, 0.0);
182 // Workspace query
183 DataType wkopt = 0.0;
184 int lwork = -1;
185 // SVD job: all singular vectors
186 char jobu = 'A';
187 char jobvt = 'A';
188 //
189 // DGESVD computes the singular value decomposition (SVD) of a real
190 // M-by-N matrix A, optionally computing the left and/or right singular
191 // vectors. The SVD is written
192 //
193 // A = U * SIGMA * transpose(V)
194 //
195 // where SIGMA is an M-by-N matrix which is zero except for its
196 // min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
197 // V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA
198 // are the singular values of A; they are real and non-negative, and
199 // are returned in descending order. The first min(m,n) columns of
200 // U and V are the left and right singular vectors of A.
201 //
202 // Note that the routine returns V**T, not V.
203 //
204 if (sizeof(DataType) == sizeof(NekDouble))
205 {
206 Lapack::Dgesvd(jobu, jobvt, m, n, (double *)a.data(), lda,
207 (double *)s.data(), (double *)u.data(), ldu,
208 (double *)vt.data(), ldvt, (double *)&wkopt, lwork,
209 info);
210 }
211 else if (sizeof(DataType) == sizeof(NekSingle))
212 {
213 Lapack::Sgesvd(jobu, jobvt, m, n, (float *)a.data(), lda,
214 (float *)s.data(), (float *)u.data(), ldu,
215 (float *)vt.data(), ldvt, (float *)&wkopt, lwork,
216 info);
217 }
218 else
219 {
220 ASSERTL0(
221 false,
222 "PseudoInverse DataType is neither NekDouble nor NekSingle");
223 }
224 lwork = static_cast<int>(wkopt);
225 Array<OneD, DataType> work(lwork);
226 if (sizeof(DataType) == sizeof(NekDouble))
227 {
228 Lapack::Dgesvd(jobu, jobvt, m, n, (double *)a.data(), lda,
229 (double *)s.data(), (double *)u.data(), ldu,
230 (double *)vt.data(), ldvt, (double *)work.data(),
231 lwork, info);
232 }
233 else if (sizeof(DataType) == sizeof(NekSingle))
234 {
235 Lapack::Sgesvd(jobu, jobvt, m, n, (float *)a.data(), lda,
236 (float *)s.data(), (float *)u.data(), ldu,
237 (float *)vt.data(), ldvt, (float *)work.data(),
238 lwork, info);
239 }
240 if (info < 0)
241 {
242 std::string message = "ERROR: The " + std::to_string(-info) +
243 "th parameter had an illegal value for gesvd";
244 ASSERTL0(false, message.c_str());
245 }
246 else if (info > 0)
247 {
248 std::string message = "ERROR: SVD did not converge in gesvd";
249 ASSERTL0(false, message.c_str());
250 }
251 // Compute pseudo-inverse: A^+ = V * S^+ * U^T
252 // S^+ is diagonal with 1/s_i for s_i > tol, 0 otherwise
253 DataType tol = std::numeric_limits<DataType>::epsilon() *
254 std::max(m, n) * std::abs(s[0]);
255 Array<OneD, DataType> sinv(min_mn, 0.0);
256 for (int i = 0; i < min_mn; ++i)
257 {
258 if (std::abs(s[i]) > tol)
259 {
260 sinv[i] = 1.0 / s[i];
261 }
262 else
263 {
264 sinv[i] = 0.0;
265 }
266 }
267 // temp = S^+ * U^T (size min_mn x m)
268 Array<OneD, DataType> temp(min_mn * m, 0.0);
269 // scale each row of U^T by sinv
270 for (int i = 0; i < min_mn; ++i)
271 {
272 for (int j = 0; j < m; ++j)
273 {
274 temp[i + j * min_mn] = sinv[i] * u[j + i * ldu];
275 }
276 }
277 // Use BLAS GEMM: data = V (n x min_mn) * temp (min_mn x m)
278 data = Array<OneD, DataType>(n * m, 0.0);
279 char transa = 'T';
280 char transb = 'N';
281 DataType alpha = 1.0;
282 DataType beta = 0.0;
283 if (sizeof(DataType) == sizeof(NekDouble))
284 {
285 Blas::Gemm(transa, transb, n, m, min_mn, alpha, (double *)vt.data(),
286 ldvt, (double *)temp.data(), min_mn, beta,
287 (double *)data.data(), n);
288 }
289 else if (sizeof(DataType) == sizeof(NekSingle))
290 {
291 Blas::Gemm(transa, transb, n, m, min_mn, (float)alpha,
292 (float *)vt.data(), ldvt, (float *)temp.data(), min_mn,
293 (float)beta, (float *)data.data(), n);
294 }
295 }
Blas::Gemm
static void Gemm(const char &transa, const char &transb, const int &m, const int &n, const int &k, const double &alpha, const double *a, const int &lda, const double *b, const int &ldb, const double &beta, double *c, const int &ldc)
BLAS level 3: Matrix-matrix multiply C = A x B where op(A)[m x k], op(B)[k x n], C[m x n] DGEMM perfo...
Definition Blas.hpp:296
Lapack::Sgesvd
static void Sgesvd(const char &jobu, const char &jobvt, const int &m, const int &n, float *a, const int &lda, float *s, float *u, const int &ldu, float *vt, const int &ldvt, float *work, const int &lwork, int &info)
Singular Value Decomposition (SVD) of a real matrix (single precision)
Definition Lapack.hpp:397
Lapack::Dgesvd
static void Dgesvd(const char &jobu, const char &jobvt, const int &m, const int &n, double *a, const int &lda, double *s, double *u, const int &ldu, double *vt, const int &ldvt, double *work, const int &lwork, int &info)
Singular Value Decomposition (SVD) of a real matrix (double precision)
Definition Lapack.hpp:385
Nektar::LibUtilities::beta
@ beta
Gauss Radau pinned at x=-1,.
Definition PointsType.h:59

References ASSERTL0, Nektar::LibUtilities::beta, Lapack::Dgesvd(), Blas::Gemm(), and Lapack::Sgesvd().

Referenced by Nektar::NekMatrix< DataType, StandardMatrixTag >::Invert().

Generated on Sat Jul 4 2026 01:02:26 for Nektar++ by doxygen 1.9.8