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