Cantera  2.3.0
SquareMatrix.cpp
Go to the documentation of this file.
1 //! @file SquareMatrix.cpp
2 
3 // This file is part of Cantera. See License.txt in the top-level directory or
4 // at http://www.cantera.org/license.txt for license and copyright information.
5 
8 #if CT_USE_LAPACK
10 #endif
11 
12 using namespace std;
13 
14 namespace Cantera
15 {
16 
17 SquareMatrix::SquareMatrix() :
18  a1norm_(0.0),
19  useQR_(0)
20 {
21  warn_deprecated("class SquareMatrix",
22  "Use class DenseMatrix instead. To be removed after Cantera 2.3.");
23 }
24 
25 SquareMatrix::SquareMatrix(size_t n, doublereal v) :
26  DenseMatrix(n, n, v),
27  a1norm_(0.0),
28  useQR_(0)
29 {
30  warn_deprecated("class SquareMatrix",
31  "Use class DenseMatrix instead. To be removed after Cantera 2.3.");
32 }
33 
35  DenseMatrix(y),
36  GeneralMatrix(y),
37  a1norm_(y.a1norm_),
38  useQR_(y.useQR_)
39 {
40 }
41 
42 SquareMatrix& SquareMatrix::operator=(const SquareMatrix& y)
43 {
44  if (&y == this) {
45  return *this;
46  }
47  DenseMatrix::operator=(y);
48  GeneralMatrix::operator=(y);
49  a1norm_ = y.a1norm_;
50  useQR_ = y.useQR_;
51  return *this;
52 }
53 
54 int SquareMatrix::solve(doublereal* b, size_t nrhs, size_t ldb)
55 {
56 #if CT_USE_LAPACK
57  if (useQR_) {
58  return solveQR(b);
59  }
60  int info=0;
61 
62  // Check to see whether the matrix has been factored.
63  if (!m_factored) {
64  int retn = factor();
65  if (retn) {
66  return retn;
67  }
68  }
69  if (ldb == 0) {
70  ldb = nColumns();
71  }
72 
73  // Solve the factored system
74  ct_dgetrs(ctlapack::NoTranspose, static_cast<int>(nRows()),
75  nrhs, &*begin(), static_cast<int>(nRows()),
76  ipiv().data(), b, ldb, info);
77  if (info != 0) {
78  if (m_printLevel) {
79  writelogf("SquareMatrix::solve(): DGETRS returned INFO = %d\n", info);
80  }
81  if (! m_useReturnErrorCode) {
82  throw CanteraError("SquareMatrix::solve()", "DGETRS returned INFO = {}", info);
83  }
84  }
85  return info;
86 #else
87  throw CanteraError("SquareMatrix::solve",
88  "Not Implemented when LAPACK is not available");
89 #endif
90 }
91 
93 {
94  m_data.assign(m_data.size(), 0.0);
95 }
96 
97 void SquareMatrix::resize(size_t n, size_t m, doublereal v)
98 {
99  DenseMatrix::resize(n, m, v);
100 }
101 
102 void SquareMatrix::mult(const doublereal* b, doublereal* prod) const
103 {
104  DenseMatrix::mult(b, prod);
105 }
106 
107 void SquareMatrix::mult(const DenseMatrix& b, DenseMatrix& prod) const
108 {
109  DenseMatrix::mult(b, prod);
110 }
111 
112 void SquareMatrix::leftMult(const doublereal* const b, doublereal* const prod) const
113 {
114  DenseMatrix::leftMult(b, prod);
115 }
116 
118 {
119 #if CT_USE_LAPACK
120  if (useQR_) {
121  return factorQR();
122  }
123  a1norm_ = ct_dlange('1', m_nrows, m_nrows, &*begin(), m_nrows, 0);
124  integer n = static_cast<int>(nRows());
125  int info=0;
126  m_factored = 1;
127  ct_dgetrf(n, n, &*begin(), static_cast<int>(nRows()), ipiv().data(), info);
128  if (info != 0) {
129  if (m_printLevel) {
130  writelogf("SquareMatrix::factor(): DGETRS returned INFO = %d\n", info);
131  }
132  if (! m_useReturnErrorCode) {
133  throw CanteraError("SquareMatrix::factor()", "DGETRS returned INFO = {}", info);
134  }
135  }
136  return info;
137 #else
138  throw CanteraError("SquareMatrix::factor",
139  "Not Implemented when LAPACK is not available");
140 #endif
141 }
142 
144 {
145  m_factored = 1;
146 }
147 
149 {
150 #if CT_USE_LAPACK
151  if (tau.size() < m_nrows) {
152  tau.resize(m_nrows, 0.0);
153  work.resize(8 * m_nrows, 0.0);
154  }
155  a1norm_ = ct_dlange('1', m_nrows, m_nrows, &*begin(), m_nrows, work.data());
156  int info = 0;
157  m_factored = 2;
158  size_t lwork = work.size();
159  ct_dgeqrf(m_nrows, m_nrows, &*begin(), m_nrows, tau.data(), work.data(), lwork, info);
160  if (info != 0) {
161  if (m_printLevel) {
162  writelogf("SquareMatrix::factorQR(): DGEQRF returned INFO = %d\n", info);
163  }
164  if (! m_useReturnErrorCode) {
165  throw CanteraError("SquareMatrix::factorQR()", "DGEQRF returned INFO = {}", info);
166  }
167  }
168  size_t lworkOpt = static_cast<size_t>(work[0]);
169  if (lworkOpt > lwork) {
170  work.resize(lworkOpt);
171  }
172  return info;
173 #else
174  throw CanteraError("SquareMatrix::factorQR",
175  "Not Implemented when LAPACK is not available");
176 #endif
177 }
178 
179 int SquareMatrix::solveQR(doublereal* b)
180 {
181 #if CT_USE_LAPACK
182  int info=0;
183 
184  // Check to see whether the matrix has been factored.
185  if (!m_factored) {
186  int retn = factorQR();
187  if (retn) {
188  return retn;
189  }
190  }
191 
192  size_t lwork = work.size();
193  if (lwork < m_nrows) {
194  work.resize(8 * m_nrows, 0.0);
195  lwork = 8 * m_nrows;
196  }
197 
198  // Solve the factored system
199  ct_dormqr(ctlapack::Left, ctlapack::Transpose, m_nrows, 1, m_nrows, &*begin(), m_nrows, tau.data(), b, m_nrows,
200  work.data(), lwork, info);
201  if (info != 0) {
202  if (m_printLevel) {
203  writelogf("SquareMatrix::solveQR(): DORMQR returned INFO = %d\n", info);
204  }
205  if (! m_useReturnErrorCode) {
206  throw CanteraError("SquareMatrix::solveQR()", "DORMQR returned INFO = {}", info);
207  }
208  }
209  size_t lworkOpt = static_cast<size_t>(work[0]);
210  if (lworkOpt > lwork) {
211  work.resize(lworkOpt);
212  }
213 
214  char dd = 'N';
215  ct_dtrtrs(ctlapack::UpperTriangular, ctlapack::NoTranspose, &dd, m_nrows, 1, &*begin(), m_nrows, b,
216  m_nrows, info);
217  if (info != 0) {
218  if (m_printLevel) {
219  writelogf("SquareMatrix::solveQR(): DTRTRS returned INFO = %d\n", info);
220  }
221  if (! m_useReturnErrorCode) {
222  throw CanteraError("SquareMatrix::solveQR()", "DTRTRS returned INFO = {}", info);
223  }
224  }
225  return info;
226 #else
227  throw CanteraError("SquareMatrix::solveQR",
228  "Not Implemented when LAPACK is not available");
229 #endif
230 }
231 
232 doublereal SquareMatrix::rcond(doublereal anorm)
233 {
234 #if CT_USE_LAPACK
235  if (iwork_.size() < m_nrows) {
236  iwork_.resize(m_nrows);
237  }
238  if (work.size() <4 * m_nrows) {
239  work.resize(4 * m_nrows);
240  }
241  doublereal rcond = 0.0;
242  if (m_factored != 1) {
243  throw CanteraError("SquareMatrix::rcond()", "matrix isn't factored correctly");
244  }
245 
246  int rinfo = 0;
247  rcond = ct_dgecon('1', m_nrows, &*begin(), m_nrows, anorm, work.data(),
248  iwork_.data(), rinfo);
249  if (rinfo != 0) {
250  if (m_printLevel) {
251  writelogf("SquareMatrix::rcond(): DGECON returned INFO = %d\n", rinfo);
252  }
253  if (! m_useReturnErrorCode) {
254  throw CanteraError("SquareMatrix::rcond()", "DGECON returned INFO = {}", rinfo);
255  }
256  }
257  return rcond;
258 #else
259  throw CanteraError("SquareMatrix::rcond",
260  "Not Implemented when LAPACK is not available");
261 #endif
262 }
263 
264 doublereal SquareMatrix::oneNorm() const
265 {
266  return a1norm_;
267 }
268 
270 {
271 #if CT_USE_LAPACK
272  if (iwork_.size() < m_nrows) {
273  iwork_.resize(m_nrows);
274  }
275  if (work.size() <3 * m_nrows) {
276  work.resize(3 * m_nrows);
277  }
278  doublereal rcond = 0.0;
279  if (m_factored != 2) {
280  throw CanteraError("SquareMatrix::rcondQR()", "matrix isn't factored correctly");
281  }
282 
283  int rinfo = 0;
284  rcond = ct_dtrcon(0, ctlapack::UpperTriangular, 0, m_nrows, &*begin(), m_nrows, work.data(),
285  iwork_.data(), rinfo);
286  if (rinfo != 0) {
287  if (m_printLevel) {
288  writelogf("SquareMatrix::rcondQR(): DTRCON returned INFO = %d\n", rinfo);
289  }
290  if (! m_useReturnErrorCode) {
291  throw CanteraError("SquareMatrix::rcondQR()", "DTRCON returned INFO = {}", rinfo);
292  }
293  }
294  return rcond;
295 #else
296  throw CanteraError("SquareMatrix::rcondQR",
297  "Not Implemented when LAPACK is not available");
298 #endif
299 }
300 
302 {
303  useQR_ = fAlgorithm;
304 }
305 
307 {
308  return (int) useQR_;
309 }
310 
311 doublereal* SquareMatrix::ptrColumn(size_t j)
312 {
313  return Array2D::ptrColumn(j);
314 }
315 
316 size_t SquareMatrix::nRows() const
317 {
318  return m_nrows;
319 }
320 
321 size_t SquareMatrix::nRowsAndStruct(size_t* const iStruct) const
322 {
323  warn_deprecated("SquareMatrix::nRowsAndStruct",
324  "To be removed after Cantera 2.3.");
325  return m_nrows;
326 }
327 
329 {
330  return new SquareMatrix(*this);
331 }
332 
333 vector_fp::iterator SquareMatrix::begin()
334 {
335  return m_data.begin();
336 }
337 
338 vector_fp::const_iterator SquareMatrix::begin() const
339 {
340  return m_data.begin();
341 }
342 
343 doublereal* const* SquareMatrix::colPts()
344 {
345  return DenseMatrix::colPts();
346 }
347 
348 size_t SquareMatrix::checkRows(doublereal& valueSmall) const
349 {
350  valueSmall = 1.0E300;
351  size_t iSmall = npos;
352  for (size_t i = 0; i < m_nrows; i++) {
353  double valueS = 0.0;
354  for (size_t j = 0; j < m_nrows; j++) {
355  valueS = std::max(fabs(value(i,j)), valueS);
356  }
357  if (valueS < valueSmall) {
358  iSmall = i;
359  valueSmall = valueS;
360  }
361  }
362  return iSmall;
363 }
364 
365 size_t SquareMatrix::checkColumns(doublereal& valueSmall) const
366 {
367  valueSmall = 1.0E300;
368  size_t jSmall = npos;
369  for (size_t j = 0; j < m_nrows; j++) {
370  double valueS = 0.0;
371  for (size_t i = 0; i < m_nrows; i++) {
372  valueS = std::max(fabs(value(i,j)), valueS);
373  }
374  if (valueS < valueSmall) {
375  jSmall = j;
376  valueSmall = valueS;
377  }
378  }
379  return jSmall;
380 }
381 
382 }
Dense, Square (not sparse) matrices.
virtual size_t checkRows(doublereal &valueSmall) const
Check to see if we have any zero rows in the Jacobian.
virtual doublereal rcond(doublereal a1norm)
Returns an estimate of the inverse of the condition number for the matrix.
vector_fp m_data
Data stored in a single array.
Definition: Array.h:330
vector_fp work
Work vector for QR algorithm.
Definition: SquareMatrix.h:116
virtual size_t checkColumns(doublereal &valueSmall) const
Check to see if we have any zero columns in the Jacobian.
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
int m_factored
Indicates whether the matrix is factored.
void warn_deprecated(const std::string &method, const std::string &extra)
Print a warning indicating that method is deprecated.
Definition: global.cpp:54
STL namespace.
Generic matrix.
Definition: GeneralMatrix.h:21
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
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.
size_t nRowsAndStruct(size_t *const iStruct=0) const
Return the size and structure of the matrix.
vector_fp tau
Work vector for QR algorithm.
Definition: SquareMatrix.h:113
virtual vector_fp::iterator begin()
Return an iterator pointing to the first element.
int solve(doublereal *b, size_t nrhs=1, size_t ldb=0)
Solves the Ax = b system returning x in the b spot.
virtual size_t nRows() const
Return the number of rows in the matrix.
void zero()
Zero the matrix.
A class for full (non-sparse) matrices with Fortran-compatible data storage.
Definition: SquareMatrix.h:20
int factor()
Factors the A matrix, overwriting A.
size_t m_nrows
Number of rows.
Definition: Array.h:333
doublereal ct_dtrcon(const char *norm, ctlapack::upperlower_t uplot, const char *diag, size_t n, doublereal *a, size_t lda, doublereal *work, int *iwork, int &info)
Definition: ctlapack.h:407
virtual int factorAlgorithm() const
Returns the factor algorithm used.
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:65
virtual void useFactorAlgorithm(int fAlgorithm)
Change the way the matrix is factored.
int solveQR(doublereal *b)
Solves the linear problem Ax=b using the QR algorithm returning x in the b spot.
vector_fp & data()
Return a reference to the data vector.
Definition: Array.h:299
virtual GeneralMatrix * duplMyselfAsGeneralMatrix() const
Duplicator member function.
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
vector_int iwork_
Integer work vector for QR algorithms.
Definition: SquareMatrix.h:119
virtual doublereal oneNorm() const
Calculate the one norm of the matrix.
Contains declarations for string manipulation functions within Cantera.
virtual doublereal * ptrColumn(size_t j)
Return a pointer to the top of column j, columns are assumed to be contiguous in memory.
virtual void leftMult(const doublereal *const b, doublereal *const prod) const
Left-multiply the matrix by transpose(b), and write the result to prod.
virtual doublereal rcondQR()
Returns an estimate of the inverse of the condition number for the matrix.
int useQR_
Use the QR algorithm to factor and invert the matrix.
Definition: SquareMatrix.h:126
int m_useReturnErrorCode
Error Handling Flag.
Definition: DenseMatrix.h:162
void setFactorFlag()
set the factored flag
Namespace for the Cantera kernel.
Definition: application.cpp:29
doublereal a1norm_
1-norm of the matrix.
Definition: SquareMatrix.h:123
size_t nColumns() const
Number of columns.
Definition: Array.h:274
vector_int & ipiv()
Return a changeable value of the pivot vector.
SquareMatrix()
Base Constructor.
A class for full (non-sparse) matrices with Fortran-compatible data storage, which adds matrix operat...
Definition: DenseMatrix.h:72
virtual int factorQR()
Factors the A matrix using the QR algorithm, overwriting A.