Cantera  2.1.2
SquareMatrix.cpp
Go to the documentation of this file.
1 /**
2  * @file SquareMatrix.cpp
3  */
4
5 /*
6  * Copyright 2004 Sandia Corporation. Under the terms of Contract
7  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
8  * retains certain rights in this software.
9  * See file License.txt for licensing information.
10  */
11
15 #include "cantera/base/global.h"
16
17 using namespace std;
18
19 namespace Cantera
20 {
21
22 SquareMatrix::SquareMatrix() :
23  DenseMatrix(),
24  GeneralMatrix(0),
25  m_factored(0),
26  a1norm_(0.0),
27  useQR_(0)
28 {
29 }
30
31 SquareMatrix::SquareMatrix(size_t n, doublereal v) :
32  DenseMatrix(n, n, v),
33  GeneralMatrix(0),
34  m_factored(0),
35  a1norm_(0.0),
36  useQR_(0)
37
38 {
39 }
40
42  DenseMatrix(y),
43  GeneralMatrix(0),
44  m_factored(y.m_factored),
45  a1norm_(y.a1norm_),
46  useQR_(y.useQR_)
47 {
48 }
49
51 {
52  if (&y == this) {
53  return *this;
54  }
58  a1norm_ = y.a1norm_;
59  useQR_ = y.useQR_;
60  return *this;
61 }
62
63 int SquareMatrix::solve(doublereal* b)
64 {
65  if (useQR_) {
66  return solveQR(b);
67  }
68  int info=0;
69  /*
70  * Check to see whether the matrix has been factored.
71  */
72  if (!m_factored) {
73  int retn = factor();
74  if (retn) {
75  return retn;
76  }
77  }
78  /*
79  * Solve the factored system
80  */
81  ct_dgetrs(ctlapack::NoTranspose, static_cast<int>(nRows()),
82  1, &(*(begin())), static_cast<int>(nRows()),
83  DATA_PTR(ipiv()), b, static_cast<int>(nColumns()), info);
84  if (info != 0) {
85  if (m_printLevel) {
86  writelogf("SquareMatrix::solve(): DGETRS returned INFO = %d\n", info);
87  }
88  if (! m_useReturnErrorCode) {
89  throw CELapackError("SquareMatrix::solve()", "DGETRS returned INFO = " + int2str(info));
90  }
91  }
92  return info;
93 }
94
96 {
97  size_t n = nRows();
98  if (n > 0) {
99  size_t nn = n * n;
100  double* sm = &m_data[0];
101  /*
102  * Using memset is the fastest way to zero a contiguous
103  * section of memory.
104  */
105  (void) memset((void*) sm, 0, nn * sizeof(double));
106  }
107 }
108
109 void SquareMatrix::resize(size_t n, size_t m, doublereal v)
110 {
111  DenseMatrix::resize(n, m, v);
112 }
113
114 void SquareMatrix::mult(const doublereal* b, doublereal* prod) const
115 {
116  DenseMatrix::mult(b, prod);
117 }
118
119 void SquareMatrix::mult(const DenseMatrix& b, DenseMatrix& prod) const
120 {
121  DenseMatrix::mult(b, prod);
122 }
123
124 void SquareMatrix::leftMult(const doublereal* const b, doublereal* const prod) const
125 {
126  DenseMatrix::leftMult(b, prod);
127 }
128
130 {
131  if (useQR_) {
132  return factorQR();
133  }
134  a1norm_ = ct_dlange('1', m_nrows, m_nrows, &(*(begin())), m_nrows, 0);
135  integer n = static_cast<int>(nRows());
136  int info=0;
137  m_factored = 1;
138  ct_dgetrf(n, n, &(*(begin())), static_cast<int>(nRows()), DATA_PTR(ipiv()), info);
139  if (info != 0) {
140  if (m_printLevel) {
141  writelogf("SquareMatrix::factor(): DGETRS returned INFO = %d\n", info);
142  }
143  if (! m_useReturnErrorCode) {
144  throw CELapackError("SquareMatrix::factor()", "DGETRS returned INFO = "+int2str(info));
145  }
146  }
147  return info;
148 }
149
151 {
152  m_factored = 0;
153 }
154
156 {
157  m_factored = 1;
158 }
159
161 {
162  if (tau.size() < m_nrows) {
163  tau.resize(m_nrows, 0.0);
164  work.resize(8 * m_nrows, 0.0);
165  }
166  a1norm_ = ct_dlange('1', m_nrows, m_nrows, &(*(begin())), m_nrows, DATA_PTR(work));
167  int info = 0;
168  m_factored = 2;
169  size_t lwork = work.size();
170  ct_dgeqrf(m_nrows, m_nrows, &(*(begin())), m_nrows, DATA_PTR(tau), DATA_PTR(work), lwork, info);
171  if (info != 0) {
172  if (m_printLevel) {
173  writelogf("SquareMatrix::factorQR(): DGEQRF returned INFO = %d\n", info);
174  }
175  if (! m_useReturnErrorCode) {
176  throw CELapackError("SquareMatrix::factorQR()", "DGEQRF returned INFO = " + int2str(info));
177  }
178  }
179  size_t lworkOpt = static_cast<size_t>(work[0]);
180  if (lworkOpt > lwork) {
181  work.resize(lworkOpt);
182  }
183
184
185  return info;
186 }
187
188 int SquareMatrix::solveQR(doublereal* b)
189 {
190  int info=0;
191  /*
192  * Check to see whether the matrix has been factored.
193  */
194  if (!m_factored) {
195  int retn = factorQR();
196  if (retn) {
197  return retn;
198  }
199  }
200
201  size_t lwork = work.size();
202  if (lwork < m_nrows) {
203  work.resize(8 * m_nrows, 0.0);
204  lwork = 8 * m_nrows;
205  }
206
207  /*
208  * Solve the factored system
209  */
210  ct_dormqr(ctlapack::Left, ctlapack::Transpose, m_nrows, 1, m_nrows, &(*(begin())), m_nrows, DATA_PTR(tau), b, m_nrows,
211  DATA_PTR(work), lwork, info);
212  if (info != 0) {
213  if (m_printLevel) {
214  writelogf("SquareMatrix::solveQR(): DORMQR returned INFO = %d\n", info);
215  }
216  if (! m_useReturnErrorCode) {
217  throw CELapackError("SquareMatrix::solveQR()", "DORMQR returned INFO = " + int2str(info));
218  }
219  }
220  size_t lworkOpt = static_cast<size_t>(work[0]);
221  if (lworkOpt > lwork) {
222  work.resize(lworkOpt);
223  }
224
225  char dd = 'N';
226
227  ct_dtrtrs(ctlapack::UpperTriangular, ctlapack::NoTranspose, &dd, m_nrows, 1, &(*(begin())), m_nrows, b,
228  m_nrows, info);
229  if (info != 0) {
230  if (m_printLevel) {
231  writelogf("SquareMatrix::solveQR(): DTRTRS returned INFO = %d\n", info);
232  }
233  if (! m_useReturnErrorCode) {
234  throw CELapackError("SquareMatrix::solveQR()", "DTRTRS returned INFO = " + int2str(info));
235  }
236  }
237
238  return info;
239 }
240
241 doublereal SquareMatrix::rcond(doublereal anorm)
242 {
243
244  if (iwork_.size() < m_nrows) {
245  iwork_.resize(m_nrows);
246  }
247  if (work.size() <4 * m_nrows) {
248  work.resize(4 * m_nrows);
249  }
250  doublereal rcond = 0.0;
251  if (m_factored != 1) {
252  throw CELapackError("SquareMatrix::rcond()", "matrix isn't factored correctly");
253  }
254
255  // doublereal anorm = ct_dlange('1', m_nrows, m_nrows, &(*(begin())), m_nrows, DATA_PTR(work));
256
257
258  int rinfo = 0;
259  rcond = ct_dgecon('1', m_nrows, &(*(begin())), m_nrows, anorm, DATA_PTR(work),
260  DATA_PTR(iwork_), rinfo);
261  if (rinfo != 0) {
262  if (m_printLevel) {
263  writelogf("SquareMatrix::rcond(): DGECON returned INFO = %d\n", rinfo);
264  }
265  if (! m_useReturnErrorCode) {
266  throw CELapackError("SquareMatrix::rcond()", "DGECON returned INFO = " + int2str(rinfo));
267  }
268  }
269  return rcond;
270 }
271
272 doublereal SquareMatrix::oneNorm() const
273 {
274  return a1norm_;
275 }
276
278 {
279
280  if (iwork_.size() < m_nrows) {
281  iwork_.resize(m_nrows);
282  }
283  if (work.size() <3 * m_nrows) {
284  work.resize(3 * m_nrows);
285  }
286  doublereal rcond = 0.0;
287  if (m_factored != 2) {
288  throw CELapackError("SquareMatrix::rcondQR()", "matrix isn't factored correctly");
289  }
290
291  int rinfo = 0;
292  rcond = ct_dtrcon(0, ctlapack::UpperTriangular, 0, m_nrows, &(*(begin())), m_nrows, DATA_PTR(work),
293  DATA_PTR(iwork_), rinfo);
294  if (rinfo != 0) {
295  if (m_printLevel) {
296  writelogf("SquareMatrix::rcondQR(): DTRCON returned INFO = %d\n", rinfo);
297  }
298  if (! m_useReturnErrorCode) {
299  throw CELapackError("SquareMatrix::rcondQR()", "DTRCON returned INFO = " + int2str(rinfo));
300  }
301  }
302  return rcond;
303 }
304
306 {
307  useQR_ = fAlgorithm;
308 }
309
311 {
312  return (int) useQR_;
313 }
314
316 {
317  return (m_factored != 0);
318 }
319
320 doublereal* SquareMatrix::ptrColumn(size_t j)
321 {
322  return Array2D::ptrColumn(j);
323 }
324
326 {
327  const SquareMatrix* yy_ptr = dynamic_cast<const SquareMatrix*>(& y);
328  Array2D::copyData(*yy_ptr);
329 }
330
331 size_t SquareMatrix::nRows() const
332 {
333  return m_nrows;
334 }
335
336 size_t SquareMatrix::nRowsAndStruct(size_t* const iStruct) const
337 {
338  return m_nrows;
339 }
340
342 {
343  return new SquareMatrix(*this);
344 }
345
346 vector_fp::iterator SquareMatrix::begin()
347 {
348  return m_data.begin();
349 }
350
351 vector_fp::const_iterator SquareMatrix::begin() const
352 {
353  return m_data.begin();
354 }
355
356 doublereal* const* SquareMatrix::colPts()
357 {
358  return DenseMatrix::colPts();
359 }
360
361 size_t SquareMatrix::checkRows(doublereal& valueSmall) const
362 {
363  valueSmall = 1.0E300;
364  size_t iSmall = npos;
365  for (size_t i = 0; i < m_nrows; i++) {
366  double valueS = 0.0;
367  for (size_t j = 0; j < m_nrows; j++) {
368  if (fabs(value(i,j)) > valueS) {
369  valueS = fabs(value(i,j));
370  }
371  }
372  if (valueS < valueSmall) {
373  iSmall = i;
374  valueSmall = valueS;
375  }
376  }
377  return iSmall;
378 }
379
380 size_t SquareMatrix::checkColumns(doublereal& valueSmall) const
381 {
382  valueSmall = 1.0E300;
383  size_t jSmall = npos;
384  for (size_t j = 0; j < m_nrows; j++) {
385  double valueS = 0.0;
386  for (size_t i = 0; i < m_nrows; i++) {
387  if (fabs(value(i,j)) > valueS) {
388  valueS = fabs(value(i,j));
389  }
390  }
391  if (valueS < valueSmall) {
392  jSmall = j;
393  valueSmall = valueS;
394  }
395  }
396  return jSmall;
397 }
398
399 }
DenseMatrix & operator=(const DenseMatrix &y)
Assignment operator.
Definition: DenseMatrix.cpp:55
virtual bool factored() const
true if the current factorization is up to date with the matrix
Dense, Square (not sparse) matrices.
std::string int2str(const int n, const std::string &fmt)
Convert an int to a string using a format converter.
Definition: stringUtils.cpp:40
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:381
vector_fp work
Work vector for QR algorithm.
Definition: SquareMatrix.h:133
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:173
Generic matrix.
Definition: GeneralMatrix.h:23
void copyData(const Array2D &y)
Copy the data from one array into another without doing any checking.
Definition: Array.h:139
int m_printLevel
Print Level.
Definition: DenseMatrix.h:175
doublereal * ptrColumn(size_t j)
Return a pointer to the top of column j, columns are contiguous in memory.
Definition: Array.h:363
This file contains definitions for utility functions and text for modules, inputfiles, logs, textlogs, HTML_logs (see Input File Handling, Diagnostic Output, Writing messages to the screen and Writing HTML Logfiles).
SquareMatrix & operator=(const SquareMatrix &right)
Assignment operator.
int m_factored
the factor flag
Definition: SquareMatrix.h:126
void resize(size_t n, size_t m, doublereal v=0.0)
Resize the matrix.
Definition: DenseMatrix.cpp:71
void writelogf(const char *fmt,...)
Write a formatted message to the screen.
Definition: global.cpp:48
virtual GeneralMatrix * duplMyselfAsGeneralMatrix() const
Duplicator member function.
int solve(doublereal *b)
Solves the Ax = b system returning x in the b spot.
virtual size_t checkColumns(doublereal &valueSmall) const
Check to see if we have any zero columns in the jacobian.
vector_fp tau
Work vector for QR algorithm.
Definition: SquareMatrix.h:130
virtual vector_fp::iterator begin()
Return an iterator pointing to the first element.
virtual void leftMult(const doublereal *const b, doublereal *const prod) const
Left-multiply the matrix by transpose(b), and write the result to prod.
void zero()
Zero the matrix.
size_t nColumns() const
Number of columns.
Definition: Array.h:322
virtual doublereal oneNorm() const
Calculate the one norm of the matrix.
A class for full (non-sparse) matrices with Fortran-compatible data storage.
Definition: SquareMatrix.h:26
int factor()
Factors the A matrix, overwriting A.
size_t m_nrows
Number of rows.
Definition: Array.h:384
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:442
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.
doublereal & value(size_t i, size_t j)
Returns a changeable reference to position in the matrix.
Definition: Array.h:300
virtual void leftMult(const double *const b, double *const prod) const
Left-multiply the matrix by transpose(b), and write the result to prod.
virtual size_t checkRows(doublereal &valueSmall) const
Check to see if we have any zero rows in the jacobian.
virtual void copyData(const GeneralMatrix &y)
Copy the data from one array into another without doing any checking.
virtual size_t nRows() const
Return the number of rows in the matrix.
Contains declarations for string manipulation functions within Cantera.
#define DATA_PTR(vec)
Creates a pointer to the start of the raw data for a vector.
Definition: ct_defs.h:36
virtual int factorAlgorithm() const
Returns the factor algorithm used.
virtual doublereal * ptrColumn(size_t j)
Return a pointer to the top of column j, columns are assumed to be contiguous in memory.
virtual doublereal rcondQR()
Returns an estimate of the inverse of the condition number for the matrix.
std::vector< int > iwork_
Integer work vector for QR algorithms.
Definition: SquareMatrix.h:136
Exception thrown when an LAPACK error is encountered associated with inverting or solving a matrix...
Definition: DenseMatrix.h:34
GeneralMatrix & operator=(const GeneralMatrix &right)
Assignment operator.
int useQR_
Use the QR algorithm to factor and invert the matrix.
Definition: SquareMatrix.h:142
size_t nRowsAndStruct(size_t *const iStruct=0) const
Return the size and structure of the matrix.
int m_useReturnErrorCode
Error Handling Flag.
Definition: DenseMatrix.h:167
void setFactorFlag()
set the factored flag
doublereal a1norm_
1-norm of the matrix. This is determined immediately before every factorization
Definition: SquareMatrix.h:139
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:70
virtual int factorQR()
Factors the A matrix using the QR algorithm, overwriting A.
virtual void clearFactorFlag()
clear the factored flag