Cantera  3.1.0a1
Loading...
Searching...
No Matches
DenseMatrix.cpp
Go to the documentation of this file.
1/**
2 * @file DenseMatrix.cpp
3 */
4
5// This file is part of Cantera. See License.txt in the top-level directory or
6// at https://cantera.org/license.txt for license and copyright information.
7
10#include "cantera/base/global.h"
11#if CT_USE_LAPACK
13#else
14 #include "cantera/numerics/eigen_dense.h"
15#endif
16
17namespace Cantera
18{
19
20DenseMatrix::DenseMatrix(size_t n, size_t m, double v) :
21 Array2D(n, m, v)
22{
23 m_ipiv.resize(std::max(n, m));
24 m_colPts.resize(m);
25 if (!m_data.empty()) {
26 for (size_t j = 0; j < m; j++) {
27 m_colPts[j] = &m_data[m_nrows*j];
28 }
29 }
30}
31
33 Array2D(y)
34{
35 m_ipiv = y.ipiv();
36 m_colPts.resize(m_ncols);
37 if (!m_data.empty()) {
38 for (size_t j = 0; j < m_ncols; j++) {
39 m_colPts[j] = &m_data[m_nrows*j];
40 }
41 }
42}
43
44DenseMatrix& DenseMatrix::operator=(const DenseMatrix& y)
45{
46 if (&y == this) {
47 return *this;
48 }
49 Array2D::operator=(y);
50 m_ipiv = y.ipiv();
51 m_colPts.resize(m_ncols);
52 for (size_t j = 0; j < m_ncols; j++) {
53 m_colPts[j] = &m_data[m_nrows*j];
54 }
55 m_useReturnErrorCode = y.m_useReturnErrorCode;
56 m_printLevel = y.m_printLevel;
57 return *this;
58}
59
60void DenseMatrix::resize(size_t n, size_t m, double v)
61{
62 Array2D::resize(n,m,v);
63 m_ipiv.resize(std::max(n,m));
64 m_colPts.resize(m_ncols);
65 if (!m_data.empty()) {
66 for (size_t j = 0; j < m_ncols; j++) {
67 m_colPts[j] = &m_data[m_nrows*j];
68 }
69 }
70}
71
72double* const* DenseMatrix::colPts()
73{
74 return &m_colPts[0];
75}
76
77const double* const* DenseMatrix::const_colPts() const
78{
79 return &m_colPts[0];
80}
81
82void DenseMatrix::mult(const double* b, double* prod) const
83{
84#if CT_USE_LAPACK
85 ct_dgemv(ctlapack::ColMajor, ctlapack::NoTranspose,
86 static_cast<int>(nRows()),
87 static_cast<int>(nColumns()), 1.0, ptrColumn(0),
88 static_cast<int>(nRows()), b, 1, 0.0, prod, 1);
89#else
90 MappedMatrix mat(const_cast<double*>(m_data.data()), nRows(), nColumns());
91 MappedVector bm(const_cast<double*>(b), nColumns());
92 MappedVector pm(prod, nRows());
93 pm = mat * bm;
94#endif
95}
96
97void DenseMatrix::mult(const DenseMatrix& B, DenseMatrix& prod) const
98{
99 if (nColumns() != B.nRows()) {
100 throw CanteraError("DenseMatrix::mult",
101 "Inner matrix dimensions do not agree: {} != {}", nColumns(), B.nRows());
102 }
103 if (nRows() != prod.nRows() || B.nColumns() != prod.nColumns()) {
104 throw CanteraError("DenseMatrix::mult",
105 "Output matrix has wrong dimensions: {}x{} != {}x{}",
106 prod.nRows(), prod.nColumns(), nRows(), B.nColumns());
107 }
108 const double* const* bcols = B.const_colPts();
109 double* const* prodcols = prod.colPts();
110 for (size_t col=0; col < B.nColumns(); ++col) {
111 // Loop over ncols multiplying A*column of B and storing in
112 // corresponding prod column
113 mult(bcols[col], prodcols[col]);
114 }
115}
116
117void DenseMatrix::leftMult(const double* const b, double* const prod) const
118{
119 for (size_t n = 0; n < nColumns(); n++) {
120 double sum = 0.0;
121 for (size_t i = 0; i < nRows(); i++) {
122 sum += value(i,n)*b[i];
123 }
124 prod[n] = sum;
125 }
126}
127
128vector<int>& DenseMatrix::ipiv()
129{
130 return m_ipiv;
131}
132
133int solve(DenseMatrix& A, double* b, size_t nrhs, size_t ldb)
134{
135 if (A.nColumns() != A.nRows()) {
136 if (A.m_printLevel) {
137 writelogf("solve(DenseMatrix& A, double* b): Can only solve a square matrix\n");
138 }
139 throw CanteraError("solve(DenseMatrix& A, double* b)", "Can only solve a square matrix");
140 }
141
142 int info = 0;
143 if (ldb == 0) {
144 ldb = A.nColumns();
145 }
146 #if CT_USE_LAPACK
147 ct_dgetrf(A.nRows(), A.nColumns(), A.ptrColumn(0),
148 A.nRows(), &A.ipiv()[0], info);
149 if (info > 0) {
150 if (A.m_printLevel) {
151 writelogf("solve(DenseMatrix& A, double* b): DGETRF returned INFO = %d U(i,i) is exactly zero. The factorization has"
152 " been completed, but the factor U is exactly singular, and division by zero will occur if "
153 "it is used to solve a system of equations.\n", info);
154 }
155 if (!A.m_useReturnErrorCode) {
156 throw CanteraError("solve(DenseMatrix& A, double* b)",
157 "DGETRF returned INFO = {}. U(i,i) is exactly zero. The factorization has"
158 " been completed, but the factor U is exactly singular, and division by zero will occur if "
159 "it is used to solve a system of equations.", info);
160 }
161 return info;
162 } else if (info < 0) {
163 if (A.m_printLevel) {
164 writelogf("solve(DenseMatrix& A, double* b): DGETRF returned INFO = %d. The argument i has an illegal value\n", info);
165 }
166
167 throw CanteraError("solve(DenseMatrix& A, double* b)",
168 "DGETRF returned INFO = {}. The argument i has an illegal value", info);
169 }
170
171 ct_dgetrs(ctlapack::NoTranspose, A.nRows(), nrhs, A.ptrColumn(0),
172 A.nRows(), &A.ipiv()[0], b, ldb, info);
173 if (info != 0) {
174 if (A.m_printLevel) {
175 writelog("solve(DenseMatrix& A, double* b): DGETRS returned INFO = {}\n", info);
176 }
177 if (info < 0 || !A.m_useReturnErrorCode) {
178 throw CanteraError("solve(DenseMatrix& A, double* b)",
179 "DGETRS returned INFO = {}", info);
180 }
181 }
182 #else
183 MappedMatrix Am(&A(0,0), A.nRows(), A.nColumns());
184 #ifdef NDEBUG
185 auto lu = Am.partialPivLu();
186 #else
187 auto lu = Am.fullPivLu();
188 if (lu.nonzeroPivots() < static_cast<long int>(A.nColumns())) {
189 throw CanteraError("solve(DenseMatrix& A, double* b)",
190 "Matrix appears to be rank-deficient: non-zero pivots = {}; columns = {}",
191 lu.nonzeroPivots(), A.nColumns());
192 }
193 #endif
194 for (size_t i = 0; i < nrhs; i++) {
195 MappedVector bm(b + ldb*i, A.nColumns());
196 bm = lu.solve(bm);
197 }
198 #endif
199 return info;
200}
201
203{
204 return solve(A, b.ptrColumn(0), b.nColumns(), b.nRows());
205}
206
207void multiply(const DenseMatrix& A, const double* const b, double* const prod)
208{
209 A.mult(b, prod);
210}
211
212void increment(const DenseMatrix& A, const double* b, double* prod)
213{
214 #if CT_USE_LAPACK
215 ct_dgemv(ctlapack::ColMajor, ctlapack::NoTranspose,
216 static_cast<int>(A.nRows()), static_cast<int>(A.nColumns()), 1.0,
217 A.ptrColumn(0), static_cast<int>(A.nRows()), b, 1, 1.0, prod, 1);
218 #else
219 MappedMatrix Am(&const_cast<DenseMatrix&>(A)(0,0), A.nRows(), A.nColumns());
220 MappedVector bm(const_cast<double*>(b), A.nColumns());
221 MappedVector pm(prod, A.nRows());
222 pm += Am * bm;
223 #endif
224}
225
226int invert(DenseMatrix& A, size_t nn)
227{
228 int info=0;
229 #if CT_USE_LAPACK
230 integer n = static_cast<int>(nn != npos ? nn : A.nRows());
231 ct_dgetrf(n, n, A.ptrColumn(0), static_cast<int>(A.nRows()),
232 &A.ipiv()[0], info);
233 if (info != 0) {
234 if (A.m_printLevel) {
235 writelogf("invert(DenseMatrix& A, size_t nn): DGETRS returned INFO = %d\n", info);
236 }
237 if (! A.m_useReturnErrorCode) {
238 throw CanteraError("invert(DenseMatrix& A, size_t nn)", "DGETRS returned INFO = {}", info);
239 }
240 return info;
241 }
242
243 vector<double> work(n);
244 integer lwork = static_cast<int>(work.size());
245 ct_dgetri(n, A.ptrColumn(0), static_cast<int>(A.nRows()),
246 &A.ipiv()[0], &work[0], lwork, info);
247 if (info != 0) {
248 if (A.m_printLevel) {
249 writelogf("invert(DenseMatrix& A, size_t nn): DGETRS returned INFO = %d\n", info);
250 }
251 if (! A.m_useReturnErrorCode) {
252 throw CanteraError("invert(DenseMatrix& A, size_t nn)", "DGETRI returned INFO={}", info);
253 }
254 }
255 #else
256 MappedMatrix Am(&A(0,0), A.nRows(), A.nColumns());
257 if (nn == npos) {
258 Am = Am.inverse();
259 } else {
260 Am.topLeftCorner(nn, nn) = Am.topLeftCorner(nn, nn).inverse();
261 }
262 #endif
263 return info;
264}
265
266}
Headers for the DenseMatrix object, which deals with dense rectangular matrices and description of th...
A class for 2D arrays stored in column-major (Fortran-compatible) form.
Definition Array.h:32
vector< double > m_data
Data stored in a single array.
Definition Array.h:219
size_t m_nrows
Number of rows.
Definition Array.h:222
size_t nRows() const
Number of rows.
Definition Array.h:176
size_t m_ncols
Number of columns.
Definition Array.h:225
size_t nColumns() const
Number of columns.
Definition Array.h:181
double * ptrColumn(size_t j)
Return a pointer to the top of column j, columns are contiguous in memory.
Definition Array.h:203
double & value(size_t i, size_t j)
Returns a changeable reference to position in the matrix.
Definition Array.h:160
virtual void resize(size_t n, size_t m, double v=0.0)
Resize the array, and fill the new entries with 'v'.
Definition Array.cpp:47
Base class for exceptions thrown by Cantera classes.
A class for full (non-sparse) matrices with Fortran-compatible data storage, which adds matrix operat...
Definition DenseMatrix.h:55
virtual void leftMult(const double *const b, double *const prod) const
Left-multiply the matrix by transpose(b), and write the result to prod.
int m_printLevel
Print Level.
int m_useReturnErrorCode
Error Handling Flag.
void resize(size_t n, size_t m, double v=0.0) override
Resize the matrix.
vector< int > & ipiv()
Return a changeable value of the pivot vector.
DenseMatrix()=default
Default Constructor.
const double *const * const_colPts() const
Return a const vector of const pointers to the columns.
vector< double * > m_colPts
Vector of column pointers.
vector< int > m_ipiv
Vector of pivots. Length is equal to the max of m and n.
This file contains definitions for utility functions and text for modules, inputfiles and logging,...
void writelogf(const char *fmt, const Args &... args)
Write a formatted message to the screen.
Definition global.h:195
void writelog(const string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition global.h:175
Namespace for the Cantera kernel.
Definition AnyMap.cpp:564
void increment(const DenseMatrix &A, const double *b, double *prod)
Multiply A*b and add it to the result in prod. Uses BLAS routine DGEMV.
const size_t npos
index returned by functions to indicate "no position"
Definition ct_defs.h:180
void multiply(const DenseMatrix &A, const double *const b, double *const prod)
Multiply A*b and return the result in prod. Uses BLAS routine DGEMV.
int invert(DenseMatrix &A, size_t nn)
invert A. A is overwritten with A^-1.
int solve(DenseMatrix &A, double *b, size_t nrhs, size_t ldb)
Solve Ax = b. Array b is overwritten on exit with x.
Contains declarations for string manipulation functions within Cantera.