Cantera  2.2.1
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
16 using namespace std;
17
18 namespace Cantera
19 {
20
21 SquareMatrix::SquareMatrix() :
22  GeneralMatrix(0),
23  a1norm_(0.0),
24  useQR_(0)
25 {
26 }
27
28 SquareMatrix::SquareMatrix(size_t n, doublereal v) :
29  DenseMatrix(n, n, v),
30  GeneralMatrix(0),
31  a1norm_(0.0),
32  useQR_(0)
33
34 {
35 }
36
38  DenseMatrix(y),
39  GeneralMatrix(0),
40  a1norm_(y.a1norm_),
41  useQR_(y.useQR_)
42 {
43 }
44
46 {
47  if (&y == this) {
48  return *this;
49  }
52  a1norm_ = y.a1norm_;
53  useQR_ = y.useQR_;
54  return *this;
55 }
56
57 int SquareMatrix::solve(doublereal* b, size_t nrhs, size_t ldb)
58 {
59  if (useQR_) {
60  return solveQR(b);
61  }
62  int info=0;
63  /*
64  * Check to see whether the matrix has been factored.
65  */
66  if (!m_factored) {
67  int retn = factor();
68  if (retn) {
69  return retn;
70  }
71  }
72  if (ldb == 0) {
73  ldb = nColumns();
74  }
75  /*
76  * Solve the factored system
77  */
78  ct_dgetrs(ctlapack::NoTranspose, static_cast<int>(nRows()),
79  nrhs, &(*(begin())), static_cast<int>(nRows()),
80  DATA_PTR(ipiv()), b, ldb, info);
81  if (info != 0) {
82  if (m_printLevel) {
83  writelogf("SquareMatrix::solve(): DGETRS returned INFO = %d\n", info);
84  }
85  if (! m_useReturnErrorCode) {
86  throw CELapackError("SquareMatrix::solve()", "DGETRS returned INFO = " + int2str(info));
87  }
88  }
89  return info;
90 }
91
93 {
94  size_t n = nRows();
95  if (n > 0) {
96  size_t nn = n * n;
97  double* sm = &m_data;
98  /*
99  * Using memset is the fastest way to zero a contiguous
100  * section of memory.
101  */
102  (void) memset((void*) sm, 0, nn * sizeof(double));
103  }
104 }
105
106 void SquareMatrix::resize(size_t n, size_t m, doublereal v)
107 {
108  DenseMatrix::resize(n, m, v);
109 }
110
111 void SquareMatrix::mult(const doublereal* b, doublereal* prod) const
112 {
113  DenseMatrix::mult(b, prod);
114 }
115
116 void SquareMatrix::mult(const DenseMatrix& b, DenseMatrix& prod) const
117 {
118  DenseMatrix::mult(b, prod);
119 }
120
121 void SquareMatrix::leftMult(const doublereal* const b, doublereal* const prod) const
122 {
123  DenseMatrix::leftMult(b, prod);
124 }
125
127 {
128  if (useQR_) {
129  return factorQR();
130  }
131  a1norm_ = ct_dlange('1', m_nrows, m_nrows, &(*(begin())), m_nrows, 0);
132  integer n = static_cast<int>(nRows());
133  int info=0;
134  m_factored = 1;
135  ct_dgetrf(n, n, &(*(begin())), static_cast<int>(nRows()), DATA_PTR(ipiv()), info);
136  if (info != 0) {
137  if (m_printLevel) {
138  writelogf("SquareMatrix::factor(): DGETRS returned INFO = %d\n", info);
139  }
140  if (! m_useReturnErrorCode) {
141  throw CELapackError("SquareMatrix::factor()", "DGETRS returned INFO = "+int2str(info));
142  }
143  }
144  return info;
145 }
146
148 {
149  m_factored = 1;
150 }
151
153 {
154  if (tau.size() < m_nrows) {
155  tau.resize(m_nrows, 0.0);
156  work.resize(8 * m_nrows, 0.0);
157  }
158  a1norm_ = ct_dlange('1', m_nrows, m_nrows, &(*(begin())), m_nrows, DATA_PTR(work));
159  int info = 0;
160  m_factored = 2;
161  size_t lwork = work.size();
162  ct_dgeqrf(m_nrows, m_nrows, &(*(begin())), m_nrows, DATA_PTR(tau), DATA_PTR(work), lwork, info);
163  if (info != 0) {
164  if (m_printLevel) {
165  writelogf("SquareMatrix::factorQR(): DGEQRF returned INFO = %d\n", info);
166  }
167  if (! m_useReturnErrorCode) {
168  throw CELapackError("SquareMatrix::factorQR()", "DGEQRF returned INFO = " + int2str(info));
169  }
170  }
171  size_t lworkOpt = static_cast<size_t>(work);
172  if (lworkOpt > lwork) {
173  work.resize(lworkOpt);
174  }
175
176
177  return info;
178 }
179
180 int SquareMatrix::solveQR(doublereal* b)
181 {
182  int info=0;
183  /*
184  * Check to see whether the matrix has been factored.
185  */
186  if (!m_factored) {
187  int retn = factorQR();
188  if (retn) {
189  return retn;
190  }
191  }
192
193  size_t lwork = work.size();
194  if (lwork < m_nrows) {
195  work.resize(8 * m_nrows, 0.0);
196  lwork = 8 * m_nrows;
197  }
198
199  /*
200  * Solve the factored system
201  */
202  ct_dormqr(ctlapack::Left, ctlapack::Transpose, m_nrows, 1, m_nrows, &(*(begin())), m_nrows, DATA_PTR(tau), b, m_nrows,
203  DATA_PTR(work), lwork, info);
204  if (info != 0) {
205  if (m_printLevel) {
206  writelogf("SquareMatrix::solveQR(): DORMQR returned INFO = %d\n", info);
207  }
208  if (! m_useReturnErrorCode) {
209  throw CELapackError("SquareMatrix::solveQR()", "DORMQR returned INFO = " + int2str(info));
210  }
211  }
212  size_t lworkOpt = static_cast<size_t>(work);
213  if (lworkOpt > lwork) {
214  work.resize(lworkOpt);
215  }
216
217  char dd = 'N';
218
219  ct_dtrtrs(ctlapack::UpperTriangular, ctlapack::NoTranspose, &dd, m_nrows, 1, &(*(begin())), m_nrows, b,
220  m_nrows, info);
221  if (info != 0) {
222  if (m_printLevel) {
223  writelogf("SquareMatrix::solveQR(): DTRTRS returned INFO = %d\n", info);
224  }
225  if (! m_useReturnErrorCode) {
226  throw CELapackError("SquareMatrix::solveQR()", "DTRTRS returned INFO = " + int2str(info));
227  }
228  }
229
230  return info;
231 }
232
233 doublereal SquareMatrix::rcond(doublereal anorm)
234 {
235
236  if (iwork_.size() < m_nrows) {
237  iwork_.resize(m_nrows);
238  }
239  if (work.size() <4 * m_nrows) {
240  work.resize(4 * m_nrows);
241  }
242  doublereal rcond = 0.0;
243  if (m_factored != 1) {
244  throw CELapackError("SquareMatrix::rcond()", "matrix isn't factored correctly");
245  }
246
247  int rinfo = 0;
248  rcond = ct_dgecon('1', m_nrows, &(*(begin())), m_nrows, anorm, DATA_PTR(work),
249  DATA_PTR(iwork_), rinfo);
250  if (rinfo != 0) {
251  if (m_printLevel) {
252  writelogf("SquareMatrix::rcond(): DGECON returned INFO = %d\n", rinfo);
253  }
254  if (! m_useReturnErrorCode) {
255  throw CELapackError("SquareMatrix::rcond()", "DGECON returned INFO = " + int2str(rinfo));
256  }
257  }
258  return rcond;
259 }
260
261 doublereal SquareMatrix::oneNorm() const
262 {
263  return a1norm_;
264 }
265
267 {
268
269  if (iwork_.size() < m_nrows) {
270  iwork_.resize(m_nrows);
271  }
272  if (work.size() <3 * m_nrows) {
273  work.resize(3 * m_nrows);
274  }
275  doublereal rcond = 0.0;
276  if (m_factored != 2) {
277  throw CELapackError("SquareMatrix::rcondQR()", "matrix isn't factored correctly");
278  }
279
280  int rinfo = 0;
281  rcond = ct_dtrcon(0, ctlapack::UpperTriangular, 0, m_nrows, &(*(begin())), m_nrows, DATA_PTR(work),
282  DATA_PTR(iwork_), rinfo);
283  if (rinfo != 0) {
284  if (m_printLevel) {
285  writelogf("SquareMatrix::rcondQR(): DTRCON returned INFO = %d\n", rinfo);
286  }
287  if (! m_useReturnErrorCode) {
288  throw CELapackError("SquareMatrix::rcondQR()", "DTRCON returned INFO = " + int2str(rinfo));
289  }
290  }
291  return rcond;
292 }
293
295 {
296  useQR_ = fAlgorithm;
297 }
298
300 {
301  return (int) useQR_;
302 }
303
304 doublereal* SquareMatrix::ptrColumn(size_t j)
305 {
306  return Array2D::ptrColumn(j);
307 }
308
310 {
311  const SquareMatrix* yy_ptr = dynamic_cast<const SquareMatrix*>(& y);
312  Array2D::copyData(*yy_ptr);
313 }
314
315 size_t SquareMatrix::nRows() const
316 {
317  return m_nrows;
318 }
319
320 size_t SquareMatrix::nRowsAndStruct(size_t* const iStruct) const
321 {
322  return m_nrows;
323 }
324
326 {
327  return new SquareMatrix(*this);
328 }
329
330 vector_fp::iterator SquareMatrix::begin()
331 {
332  return m_data.begin();
333 }
334
335 vector_fp::const_iterator SquareMatrix::begin() const
336 {
337  return m_data.begin();
338 }
339
340 doublereal* const* SquareMatrix::colPts()
341 {
342  return DenseMatrix::colPts();
343 }
344
345 size_t SquareMatrix::checkRows(doublereal& valueSmall) const
346 {
347  valueSmall = 1.0E300;
348  size_t iSmall = npos;
349  for (size_t i = 0; i < m_nrows; i++) {
350  double valueS = 0.0;
351  for (size_t j = 0; j < m_nrows; j++) {
352  valueS = std::max(fabs(value(i,j)), valueS);
353  }
354  if (valueS < valueSmall) {
355  iSmall = i;
356  valueSmall = valueS;
357  }
358  }
359  return iSmall;
360 }
361
362 size_t SquareMatrix::checkColumns(doublereal& valueSmall) const
363 {
364  valueSmall = 1.0E300;
365  size_t jSmall = npos;
366  for (size_t j = 0; j < m_nrows; j++) {
367  double valueS = 0.0;
368  for (size_t i = 0; i < m_nrows; i++) {
369  valueS = std::max(fabs(value(i,j)), valueS);
370  }
371  if (valueS < valueSmall) {
372  jSmall = j;
373  valueSmall = valueS;
374  }
375  }
376  return jSmall;
377 }
378
379 }
DenseMatrix & operator=(const DenseMatrix &y)
Assignment operator.
Definition: DenseMatrix.cpp:48
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:39
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:375
vector_fp work
Work vector for QR algorithm.
Definition: SquareMatrix.h:126
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.
Generic matrix.
Definition: GeneralMatrix.h:24
void copyData(const Array2D &y)
Copy the data from one array into another without doing any checking.
Definition: Array.h:133
int m_printLevel
Print Level.
Definition: DenseMatrix.h:176
doublereal * ptrColumn(size_t j)
Return a pointer to the top of column j, columns are contiguous in memory.
Definition: Array.h:358
SquareMatrix & operator=(const SquareMatrix &right)
Assignment operator.
void resize(size_t n, size_t m, doublereal v=0.0)
Resize the matrix.
Definition: DenseMatrix.cpp:64
void writelogf(const char *fmt,...)
Write a formatted message to the screen.
Definition: global.cpp:38
virtual GeneralMatrix * duplMyselfAsGeneralMatrix() const
Duplicator member function.
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:123
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 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:317
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:378
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:441
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:295
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)
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:129
Exception thrown when an LAPACK error is encountered associated with inverting or solving a matrix...
Definition: DenseMatrix.h:35
GeneralMatrix & operator=(const GeneralMatrix &right)
Assignment operator.
int useQR_
Use the QR algorithm to factor and invert the matrix.
Definition: SquareMatrix.h:135
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:168
void setFactorFlag()
set the factored flag
doublereal a1norm_
1-norm of the matrix. This is determined immediately before every factorization
Definition: SquareMatrix.h:132
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:71
virtual int factorQR()
Factors the A matrix using the QR algorithm, overwriting A.