Cantera  2.0
SquareMatrix.cpp
1 /**
2  * @file DenseMatrix.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 
12 #include "cantera/base/ct_defs.h"
16 #include "cantera/base/global.h"
17 
18 #include <iostream>
19 #include <vector>
20 #include <cstring>
21 
22 using namespace std;
23 
24 namespace Cantera
25 {
26 
27 //====================================================================================================================
28 SquareMatrix::SquareMatrix() :
29  DenseMatrix(),
30  GeneralMatrix(0),
31  m_factored(0),
32  a1norm_(0.0),
33  useQR_(0)
34 {
35 }
36 
37 //====================================================================================================================
38 // Constructor.
39 /*
40  * Create an \c n by \c n matrix, and initialize
41  * all elements to \c v.
42  *
43  * @param n size of the square matrix
44  * @param v initial value of all matrix components.
45  */
46 SquareMatrix::SquareMatrix(size_t n, doublereal v) :
47  DenseMatrix(n, n, v),
48  GeneralMatrix(0),
49  m_factored(0),
50  a1norm_(0.0),
51  useQR_(0)
52 
53 {
54 }
55 //====================================================================================================================
56 /*
57  *
58  * copy constructor
59  */
61  DenseMatrix(y),
62  GeneralMatrix(0),
63  m_factored(y.m_factored),
64  a1norm_(y.a1norm_),
65  useQR_(y.useQR_)
66 {
67 }
68 
69 //====================================================================================================================
70 /*
71  * Assignment operator
72  */
74 {
75  if (&y == this) {
76  return *this;
77  }
81  a1norm_ = y.a1norm_;
82  useQR_ = y.useQR_;
83  return *this;
84 }
85 //====================================================================================================================
87 {
88 }
89 //====================================================================================================================
90 /*
91  * Solve Ax = b. Vector b is overwritten on exit with x.
92  */
93 int SquareMatrix::solve(doublereal* b)
94 {
95  if (useQR_) {
96  return solveQR(b);
97  }
98  int info=0;
99  /*
100  * Check to see whether the matrix has been factored.
101  */
102  if (!m_factored) {
103  int retn = factor();
104  if (retn) {
105  return retn;
106  }
107  }
108  /*
109  * Solve the factored system
110  */
111  ct_dgetrs(ctlapack::NoTranspose, static_cast<int>(nRows()),
112  1, &(*(begin())), static_cast<int>(nRows()),
113  DATA_PTR(ipiv()), b, static_cast<int>(nColumns()), info);
114  if (info != 0) {
115  if (m_printLevel) {
116  writelogf("SquareMatrix::solve(): DGETRS returned INFO = %d\n", info);
117  }
118  if (! m_useReturnErrorCode) {
119  throw CELapackError("SquareMatrix::solve()", "DGETRS returned INFO = " + int2str(info));
120  }
121  }
122  return info;
123 }
124 //====================================================================================================================
125 /*
126  * Set all entries to zero
127  */
129 {
130  size_t n = nRows();
131  if (n > 0) {
132  size_t nn = n * n;
133  double* sm = &m_data[0];
134  /*
135  * Using memset is the fastest way to zero a contiguous
136  * section of memory.
137  */
138  (void) memset((void*) sm, 0, nn * sizeof(double));
139  }
140 }
141 //====================================================================================================================
142 void SquareMatrix::resize(size_t n, size_t m, doublereal v)
143 {
144  DenseMatrix::resize(n, m, v);
145 }
146 
147 //====================================================================================================================
148 // Multiply A*b and write result to prod.
149 /*
150  * @param b Vector to do the rh multiplication
151  * @param prod OUTPUT vector to receive the result
152  */
153 void SquareMatrix::mult(const doublereal* b, doublereal* prod) const
154 {
155  DenseMatrix::mult(b, prod);
156 }
157 //====================================================================================================================
158 // Multiply b*A and write result to prod.
159 /*
160  * @param b Vector to do the lh multiplication
161  * @param prod OUTPUT vector to receive the result
162  */
163 void SquareMatrix::leftMult(const doublereal* const b, doublereal* const prod) const
164 {
165  DenseMatrix::leftMult(b, prod);
166 }
167 //====================================================================================================================
168 /*
169  * Factor A. A is overwritten with the LU decomposition of A.
170  */
172 {
173  if (useQR_) {
174  return factorQR();
175  }
176  a1norm_ = ct_dlange('1', m_nrows, m_nrows, &(*(begin())), m_nrows, DATA_PTR(work));
177  integer n = static_cast<int>(nRows());
178  int info=0;
179  m_factored = 1;
180  ct_dgetrf(n, n, &(*(begin())), static_cast<int>(nRows()), DATA_PTR(ipiv()), info);
181  if (info != 0) {
182  if (m_printLevel) {
183  writelogf("SquareMatrix::factor(): DGETRS returned INFO = %d\n", info);
184  }
185  if (! m_useReturnErrorCode) {
186  throw CELapackError("SquareMatrix::factor()", "DGETRS returned INFO = "+int2str(info));
187  }
188  }
189  return info;
190 }
191 //=====================================================================================================================
192 /*
193  * clear the factored flag
194  */
196 {
197  m_factored = 0;
198 }
199 //=====================================================================================================================
200 /*
201  * set the factored flag
202  */
204 {
205  m_factored = 1;
206 }
207 //=====================================================================================================================
209 {
210  if (tau.size() < m_nrows) {
211  tau.resize(m_nrows, 0.0);
212  work.resize(8 * m_nrows, 0.0);
213  }
214  a1norm_ = ct_dlange('1', m_nrows, m_nrows, &(*(begin())), m_nrows, DATA_PTR(work));
215  int info = 0;
216  m_factored = 2;
217  size_t lwork = work.size();
218  ct_dgeqrf(m_nrows, m_nrows, &(*(begin())), m_nrows, DATA_PTR(tau), DATA_PTR(work), lwork, info);
219  if (info != 0) {
220  if (m_printLevel) {
221  writelogf("SquareMatrix::factorQR(): DGEQRF returned INFO = %d\n", info);
222  }
223  if (! m_useReturnErrorCode) {
224  throw CELapackError("SquareMatrix::factorQR()", "DGEQRF returned INFO = " + int2str(info));
225  }
226  }
227  size_t lworkOpt = static_cast<size_t>(work[0]);
228  if (lworkOpt > lwork) {
229  work.resize(lworkOpt);
230  }
231 
232 
233  return info;
234 }
235 //=====================================================================================================================
236 /*
237  * Solve Ax = b. Vector b is overwritten on exit with x.
238  */
239 int SquareMatrix::solveQR(doublereal* b)
240 {
241  int info=0;
242  /*
243  * Check to see whether the matrix has been factored.
244  */
245  if (!m_factored) {
246  int retn = factorQR();
247  if (retn) {
248  return retn;
249  }
250  }
251 
252  size_t lwork = work.size();
253  if (lwork < m_nrows) {
254  work.resize(8 * m_nrows, 0.0);
255  lwork = 8 * m_nrows;
256  }
257 
258  /*
259  * Solve the factored system
260  */
261  ct_dormqr(ctlapack::Left, ctlapack::Transpose, m_nrows, 1, m_nrows, &(*(begin())), m_nrows, DATA_PTR(tau), b, m_nrows,
262  DATA_PTR(work), lwork, info);
263  if (info != 0) {
264  if (m_printLevel) {
265  writelogf("SquareMatrix::solveQR(): DORMQR returned INFO = %d\n", info);
266  }
267  if (! m_useReturnErrorCode) {
268  throw CELapackError("SquareMatrix::solveQR()", "DORMQR returned INFO = " + int2str(info));
269  }
270  }
271  size_t lworkOpt = static_cast<size_t>(work[0]);
272  if (lworkOpt > lwork) {
273  work.resize(lworkOpt);
274  }
275 
276  char dd = 'N';
277 
278  ct_dtrtrs(ctlapack::UpperTriangular, ctlapack::NoTranspose, &dd, m_nrows, 1, &(*(begin())), m_nrows, b,
279  m_nrows, info);
280  if (info != 0) {
281  if (m_printLevel) {
282  writelogf("SquareMatrix::solveQR(): DTRTRS returned INFO = %d\n", info);
283  }
284  if (! m_useReturnErrorCode) {
285  throw CELapackError("SquareMatrix::solveQR()", "DTRTRS returned INFO = " + int2str(info));
286  }
287  }
288 
289  return info;
290 }
291 //=====================================================================================================================
292 doublereal SquareMatrix::rcond(doublereal anorm)
293 {
294 
295  if (iwork_.size() < m_nrows) {
296  iwork_.resize(m_nrows);
297  }
298  if (work.size() <4 * m_nrows) {
299  work.resize(4 * m_nrows);
300  }
301  doublereal rcond = 0.0;
302  if (m_factored != 1) {
303  throw CELapackError("SquareMatrix::rcond()", "matrix isn't factored correctly");
304  }
305 
306  // doublereal anorm = ct_dlange('1', m_nrows, m_nrows, &(*(begin())), m_nrows, DATA_PTR(work));
307 
308 
309  int rinfo = 0;
310  rcond = ct_dgecon('1', m_nrows, &(*(begin())), m_nrows, anorm, DATA_PTR(work),
311  DATA_PTR(iwork_), rinfo);
312  if (rinfo != 0) {
313  if (m_printLevel) {
314  writelogf("SquareMatrix::rcond(): DGECON returned INFO = %d\n", rinfo);
315  }
316  if (! m_useReturnErrorCode) {
317  throw CELapackError("SquareMatrix::rcond()", "DGECON returned INFO = " + int2str(rinfo));
318  }
319  }
320  return rcond;
321 }
322 //=====================================================================================================================
323 doublereal SquareMatrix::oneNorm() const
324 {
325  return a1norm_;
326 }
327 //=====================================================================================================================
329 {
330 
331  if (iwork_.size() < m_nrows) {
332  iwork_.resize(m_nrows);
333  }
334  if (work.size() <3 * m_nrows) {
335  work.resize(3 * m_nrows);
336  }
337  doublereal rcond = 0.0;
338  if (m_factored != 2) {
339  throw CELapackError("SquareMatrix::rcondQR()", "matrix isn't factored correctly");
340  }
341 
342  int rinfo = 0;
343  rcond = ct_dtrcon(0, ctlapack::UpperTriangular, 0, m_nrows, &(*(begin())), m_nrows, DATA_PTR(work),
344  DATA_PTR(iwork_), rinfo);
345  if (rinfo != 0) {
346  if (m_printLevel) {
347  writelogf("SquareMatrix::rcondQR(): DTRCON returned INFO = %d\n", rinfo);
348  }
349  if (! m_useReturnErrorCode) {
350  throw CELapackError("SquareMatrix::rcondQR()", "DTRCON returned INFO = " + int2str(rinfo));
351  }
352  }
353  return rcond;
354 }
355 //=====================================================================================================================
357 {
358  useQR_ = fAlgorithm;
359 }
360 //=====================================================================================================================
362 {
363  return (int) useQR_;
364 }
365 //=====================================================================================================================
367 {
368  return (m_factored != 0);
369 }
370 //=====================================================================================================================
371 // Return a pointer to the top of column j, columns are contiguous in memory
372 /*
373  * @param j Value of the column
374  *
375  * @return Returns a pointer to the top of the column
376  */
377 doublereal* SquareMatrix::ptrColumn(size_t j)
378 {
379  return Array2D::ptrColumn(j);
380 }
381 //=====================================================================================================================
382 // Copy the data from one array into another without doing any checking
383 /*
384  * This differs from the assignment operator as no resizing is done and memcpy() is used.
385  * @param y Array to be copied
386  */
388 {
389  const SquareMatrix* yy_ptr = dynamic_cast<const SquareMatrix*>(& y);
390  Array2D::copyData(*yy_ptr);
391 }
392 //=====================================================================================================================
393 size_t SquareMatrix::nRows() const
394 {
395  return m_nrows;
396 }
397 //=====================================================================================================================
398 size_t SquareMatrix::nRowsAndStruct(size_t* const iStruct) const
399 {
400  return m_nrows;
401 }
402 //=====================================================================================================================
404 {
405  SquareMatrix* dd = new SquareMatrix(*this);
406  return static_cast<GeneralMatrix*>(dd);
407 }
408 //=====================================================================================================================
409 // Return an iterator pointing to the first element
410 vector_fp::iterator SquareMatrix::begin()
411 {
412  return m_data.begin();
413 }
414 //=====================================================================================================================
415 // Return a const iterator pointing to the first element
416 vector_fp::const_iterator SquareMatrix::begin() const
417 {
418  return m_data.begin();
419 }
420 //=====================================================================================================================
421 // Return a vector of const pointers to the columns
422 /*
423  * Note the value of the pointers are protected by their being const.
424  * However, the value of the matrix is open to being changed.
425  *
426  * @return returns a vector of pointers to the top of the columns
427  * of the matrices.
428  */
429 doublereal* const* SquareMatrix::colPts()
430 {
431  return DenseMatrix::colPts();
432 }
433 //=====================================================================================================================
434 
435 size_t SquareMatrix::checkRows(doublereal& valueSmall) const
436 {
437  valueSmall = 1.0E300;
438  size_t iSmall = npos;
439  for (size_t i = 0; i < m_nrows; i++) {
440  double valueS = 0.0;
441  for (size_t j = 0; j < m_nrows; j++) {
442  if (fabs(value(i,j)) > valueS) {
443  valueS = fabs(value(i,j));
444  }
445  }
446  if (valueS < valueSmall) {
447  iSmall = i;
448  valueSmall = valueS;
449  }
450  }
451  return iSmall;
452 }
453 //=====================================================================================================================
454 size_t SquareMatrix::checkColumns(doublereal& valueSmall) const
455 {
456  valueSmall = 1.0E300;
457  size_t jSmall = npos;
458  for (size_t j = 0; j < m_nrows; j++) {
459  double valueS = 0.0;
460  for (size_t i = 0; i < m_nrows; i++) {
461  if (fabs(value(i,j)) > valueS) {
462  valueS = fabs(value(i,j));
463  }
464  }
465  if (valueS < valueSmall) {
466  jSmall = j;
467  valueSmall = valueS;
468  }
469  }
470  return jSmall;
471 }
472 //=====================================================================================================================
473 
474 
475 }
476