STORMM Source Documentation
Loading...
Searching...
No Matches
matrix_ops.h
1// -*-c++-*-
2#ifndef STORMM_MATRIX_OPS_H
3#define STORMM_MATRIX_OPS_H
4
5#include <algorithm>
6#include <cmath>
7#include <vector>
8#include "copyright.h"
9#include "Constants/behavior.h"
10#include "Constants/scaling.h"
11#include "Accelerator/hybrid.h"
12#include "FileManagement/file_enumerators.h"
13#include "Parsing/parse.h"
14#include "Reporting/error_format.h"
15#include "UnitTesting/file_snapshot.h"
16#include "matrix.h"
17#include "vector_ops.h"
18
19namespace stormm {
20namespace stmath {
21
22using constants::ExceptionResponse;
23using card::Hybrid;
24using card::HybridTargetLevel;
25using diskutil::PrintSituation;
26using parse::NumberFormat;
27using parse::PolyNumeric;
28using parse::realToString;
29using testing::writeSnapshot;
30
33constexpr int maximum_ql_iterations = 30;
34
37template <typename T> struct TwinPointer {
38 T* ptr_a;
39 T* ptr_b;
40};
41
43enum class EigenWork {
44 EIGENVALUES,
45 EIGENVECTORS,
46};
47
49enum class TransposeState {
50 AS_IS,
51 TRANSPOSE
53};
54
78template <typename T>
79void matrixMultiply(const T* a, size_t row_a, size_t col_a, const T* b, size_t row_b, size_t col_b,
80 T* c, T scale_a = 1.0, T scale_b = 1.0, T scale_c = 1.0,
81 TransposeState x_a = TransposeState::AS_IS,
82 TransposeState x_b = TransposeState::AS_IS);
83
84template <typename T>
85void matrixMultiply(const std::vector<T> &a, size_t row_a, size_t col_a, const std::vector<T> &b,
86 size_t row_b, size_t col_b, std::vector<T> *c, T scale_a = 1.0,
87 T scale_b = 1.0, T scale_c = 1.0, TransposeState x_a = TransposeState::AS_IS,
88 TransposeState x_b = TransposeState::AS_IS);
89
90template <typename T>
91void matrixMultiply(const Hybrid<T> &a, size_t row_a, size_t col_a, const Hybrid<T> &b,
92 size_t row_b, size_t col_b, Hybrid<T> *c, T scale_a = 1.0, T scale_b = 1.0,
93 T scale_c = 1.0, TransposeState x_a = TransposeState::AS_IS,
94 TransposeState x_b = TransposeState::AS_IS);
95
96template <typename T>
97void matrixMultiply(const HpcMatrix<T> &a, const HpcMatrix<T> &b, HpcMatrix<T> *c, T scale_a = 1.0,
98 T scale_b = 1.0, T scale_c = 1.0, TransposeState x_a = TransposeState::AS_IS,
99 TransposeState x_b = TransposeState::AS_IS);
101
127template <typename T>
128void matrixVectorMultiply(const T* a, const T* x, T* b, size_t row_a, size_t col_a,
129 T scale_a = 1.0, T scale_x = 1.0, T scale_b = 1.0,
130 TransposeState x_a = TransposeState::AS_IS);
131
132template <typename T>
133void matrixVectorMultiply(const std::vector<T> &a, const std::vector<T> &x, std::vector<T> *b,
134 size_t row_a, size_t col_a, T scale_a = 1.0, T scale_x = 1.0,
135 T scale_b = 1.0, TransposeState x_a = TransposeState::AS_IS);
136
137template <typename T>
138void matrixVectorMultiply(const Hybrid<T> &a, const Hybrid<T> &x, Hybrid<T> *b, size_t row_a,
139 size_t col_a, T scale_a = 1.0, T scale_x = 1.0, T scale_b = 1.0,
140 TransposeState x_a = TransposeState::AS_IS);
141
142template <typename T>
143void matrixVectorMultiply(const HpcMatrix<T> &a, const Hybrid<T> &x, Hybrid<T> *b, T scale_a = 1.0,
144 T scale_x = 1.0, T scale_b = 1.0,
145 TransposeState x_a = TransposeState::AS_IS);
147
158template <typename T>
159void transpose(T* a, size_t rows, size_t columns, size_t n_work);
160
161template <typename T>
162void transpose(const T* a, T* a_transpose, size_t rows, size_t columns);
163
164template <typename T>
165void transpose(std::vector<T> *a, size_t rows, size_t columns, size_t n_work);
166
167template <typename T>
168void transpose(const std::vector<T> &a, std::vector<T> *a_transpose, size_t rows, size_t columns);
169
170template <typename T>
171void transpose(Hybrid<T> *a, size_t rows, size_t columns, size_t n_work);
172
173template <typename T>
174void transpose(const Hybrid<T> &a, Hybrid<T> *a_transpose, size_t rows, size_t columns);
176
190template <typename T> void invertSquareMatrix(const T* matrix, T* inverse, size_t rank);
191
192template <typename T> void invertSquareMatrix(T* matrix, T* inverse, size_t rank);
193
194template <typename T> void invertSquareMatrix(const std::vector<T> &matrix,
195 std::vector<T> *inverse);
196
197template <typename T> void invertSquareMatrix(std::vector<T> *matrix, std::vector<T> *inverse);
198
199template <typename T> void invertSquareMatrix(const Hybrid<T> &matrix, Hybrid<T> *inverse);
200
201template <typename T> void invertSquareMatrix(Hybrid<T> *matrix, Hybrid<T> *inverse);
203
212void checkArraySizeToMatrixRank(size_t actual_rank, size_t rank, size_t s_a, size_t s_v,
213 ExceptionResponse policy);
214
230template <typename T>
231void jacobiEigensolver(T* a, T* v, T* d, size_t rank,
232 ExceptionResponse policy = ExceptionResponse::WARN);
233
234template <typename T>
235void jacobiEigensolver(std::vector<T> *a, std::vector<T> *v, std::vector<T> *d, size_t rank = 0,
236 ExceptionResponse policy = ExceptionResponse::WARN);
237
238template <typename T>
239void jacobiEigensolver(Hybrid<T> *a, Hybrid<T> *v, Hybrid<T> *d, size_t rank = 0,
240 ExceptionResponse policy = ExceptionResponse::WARN);
241
242template <typename T>
243void jacobiEigensolver(HpcMatrix<T> *a, HpcMatrix<T> *v, Hybrid<T> *d,
244 ExceptionResponse policy = ExceptionResponse::WARN);
246
264template <typename T>
265void realSymmEigensolver(T* amat, const int rank, T* diag, T* eigv,
266 EigenWork process = EigenWork::EIGENVECTORS);
267
268template <typename T>
269void realSymmEigensolver(std::vector<T> *amat, std::vector<T> *diag, std::vector<T> *eigv,
270 EigenWork process = EigenWork::EIGENVECTORS);
271
272template <typename T>
273void realSymmEigensolver(const T* amat, T* vmat, const int rank, T* diag, T* eigv,
274 EigenWork process = EigenWork::EIGENVECTORS);
275
276template <typename T>
277void realSymmEigensolver(const std::vector<T> &amat, std::vector<T> *vmat, std::vector<T> *diag,
278 std::vector<T> *eigv, EigenWork process = EigenWork::EIGENVECTORS);
280
293template <typename T>
294void pseudoInverse(T* amat, T* workspace, T* inv_workspace, size_t rows, size_t columns);
295
296template <typename T>
297void pseudoInverse(T* amat, size_t rows, size_t columns);
298
299template <typename T>
300void pseudoInverse(const T* amat, T* result, T* workspace, T* inv_workspace, size_t rows,
301 size_t columns);
302
303template <typename T>
304void pseudoInverse(std::vector<T> *amat, std::vector<T> *workspace, std::vector<T> *inv_workspace,
305 size_t rows, size_t columns);
306
307template <typename T>
308void pseudoInverse(std::vector<T> *amat, size_t rows, size_t columns);
309
310template <typename T>
311void pseudoInverse(const std::vector<T> &amat, std::vector<T> *result, std::vector<T> *workspace,
312 std::vector<T> *inv_workspace, size_t rows, size_t columns);
313
314template <typename T>
315std::vector<T> pseudoInverse(const std::vector<T> &amat, size_t rows, size_t columns);
316
317template <typename T>
318std::vector<T> pseudoInverse(const std::vector<T> &amat, std::vector<T> *workspace,
319 std::vector<T> *inv_workspace, size_t rows, size_t columns);
320
321template <typename T>
322void pseudoInverse(Hybrid<T> *amat, Hybrid<T> *workspace, Hybrid<T> *inv_workspace, size_t rows,
323 size_t columns);
324
325template <typename T>
326void pseudoInverse(Hybrid<T> *amat, size_t rows, size_t columns);
327
328template <typename T>
329void pseudoInverse(const Hybrid<T> &amat, Hybrid<T> *result, Hybrid<T> *workspace,
330 Hybrid<T> *inv_workspace, size_t rows, size_t columns);
332
348template <typename T>
349void qrSolver(T* amat, T* xvec, T* bvec, const size_t a_rows, const size_t a_cols);
350
351template <typename T>
352void qrSolver(std::vector<T> *amat, std::vector<T> *xvec, std::vector<T> *bvec);
353
354template <typename T>
355void qrSolver(const std::vector<T> &amat, std::vector<T> *xvec, const std::vector<T> &bvec);
356
357template <typename T>
358void qrSolver(Hybrid<T> *amat, Hybrid<T> *xvec, Hybrid<T> *bvec);
359
360template <typename T>
361void qrSolver(const Hybrid<T> &amat, Hybrid<T> *xvec, const Hybrid<T> &bvec);
363
374template <typename T>
375double leibnizDeterminant(const T* amat, const size_t rank);
376
377template <typename T>
378double leibnizDeterminant(const std::vector<T> &amat);
379
380template <typename T>
381double leibnizDeterminant(const Hybrid<T> &amat);
383
401template <typename T>
402void computeBoxTransform(T lx, T ly, T lz, T alpha, T beta, T gamma, T* umat, T* invu);
403
404template <typename T> void computeBoxTransform(const T* dims, T* umat, T* invu);
405
406template <typename T> void computeBoxTransform(const std::vector<T> &dims, std::vector<T> *umat,
407 std::vector<T> *invu);
409
419template <typename T>
420void extractBoxDimensions(T *lx, T *ly, T *lz, T *alpha, T *beta, T *gamma, const T* invu);
421
434template <typename T, typename T3>
435void rotationMatrixAboutVector(T v_x, T v_y, T v_z, const T theta, T* m);
436
437template <typename T, typename T3>
438void rotationMatrixAboutVector(const T3 v, const T theta, T* m);
439
440template <typename T, typename T3>
441void rotationMatrixAboutVector(const T3 v, const T theta, std::vector<T> *m);
442
443template <typename T, typename T3>
444std::vector<T> rotationMatrixAboutVector(const T3 v, const T theta);
445
446template <typename T>
447std::vector<T> rotationMatrixAboutVector(const std::vector<T> &v, const T theta);
448
449template <typename T>
450void rotationMatrixAboutVector(const T* v, const T theta, T* m);
451
452template <typename T>
453void rotationMatrixAboutVector(const T* v, const T theta, std::vector<T> *m);
455
471template <typename T> std::vector<T> hessianNormalWidths(const T* invu);
472
473template <typename T> std::vector<T> hessianNormalWidths(const std::vector<T> &invu);
474
475template <typename T> void hessianNormalWidths(const T* invu, T *xw, T *yw, T *zw);
476
477template <typename T> void hessianNormalWidths(const std::vector<T> &invu, T *xw, T *yw, T *zw);
479
495template <typename T>
496void printMatrix(const T* matrix, const int rows, const int cols, const std::string &varname,
497 const std::string &file_name = std::string(""),
498 const PrintSituation expectation = PrintSituation::APPEND);
499
500template <typename T>
501void printMatrix(const std::vector<T> &matrix, const int rows, const int cols,
502 const std::string &varname, const std::string &file_name = std::string(""),
503 const PrintSituation expectation = PrintSituation::APPEND);
504
505template <typename T>
506void printMatrix(const Hybrid<T> &matrix, const int rows, const int cols,
507 const std::string &varname, const std::string &file_name = std::string(""),
508 const PrintSituation expectation = PrintSituation::APPEND);
509
510template <typename T>
511void printMatrix(const HpcMatrix<T> &matrix, const std::string &varname,
512 const std::string &file_name = std::string(""),
513 const PrintSituation expectation = PrintSituation::APPEND);
515
516} // namespace stmath
517} // namespace stormm
518
519#include "matrix_ops.tpp"
520
521#endif
stormm::stmath::TwinPointer
A twin pointer, useful for performing partial-pivoting Gaussian elimination on row-major matrices.
Definition matrix_ops.h:37
stormm::stmath::TwinPointer::ptr_b
T * ptr_b
Pointer to the second stretch of data.
Definition matrix_ops.h:39
stormm::stmath::TwinPointer::ptr_a
T * ptr_a
Pointer to the first stretch of data.
Definition matrix_ops.h:38
Generated by doxygen 1.13.2