Cantera  2.0
SquareMatrix.h
Go to the documentation of this file.
1 /**
2  * @file SquareMatrix.h
3  * Dense, Square (not sparse) matrices.
4  */
5
6 /*
7  * Copyright 2004 Sandia Corporation. Under the terms of Contract
8  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
9  * retains certain rights in this software.
10  * See file License.txt for licensing information.
11  */
12
13 #ifndef CT_SQUAREMATRIX_H
14 #define CT_SQUAREMATRIX_H
15
16 #include "DenseMatrix.h"
17 #include "GeneralMatrix.h"
18
19 namespace Cantera
20 {
21
22 /**
23  * A class for full (non-sparse) matrices with Fortran-compatible
24  * data storage. Adds matrix inversion operations to this class from DenseMatrix.
25  */
26 class SquareMatrix: public DenseMatrix, public GeneralMatrix
27 {
28 public:
29  //! Base Constructor.
30  /*!
31  * Create an \c 0 by \c 0 matrix, and initialize all elements to \c 0.
32  */
33  SquareMatrix();
34
35  //! Constructor.
36  /*!
37  * Create an \c n by \c n matrix, and initialize all elements to \c v.
38  *
39  * @param n size of the square matrix
40  * @param v initial value of all matrix components.
41  */
42  SquareMatrix(size_t n, doublereal v = 0.0);
43
44  //! Copy Constructor
45  /*!
46  * @param right Object to be copied
47  */
48  SquareMatrix(const SquareMatrix& right);
49
50  //! Assignment operator
51  /*!
52  * @param right Object to be copied
53  */
54  SquareMatrix& operator=(const SquareMatrix& right);
55
56
57  //! Destructor. Does nothing.
58  virtual ~SquareMatrix();
59
60  //! Solves the Ax = b system returning x in the b spot.
61  /*!
62  * @param b Vector for the rhs of the equation system
63  */
64  int solve(doublereal* b);
65
66  //! Resize the matrix
67  /*!
68  * @param n Number of rows
69  * @param m Number of columns
70  * @param v double to fill the new space (defaults to zero)
71  */
72  void resize(size_t n, size_t m, doublereal v = 0.0);
73
74  /**
75  * Zero the matrix
76  */
77  void zero();
78
79  //! Multiply A*b and write result to prod.
80  /*!
81  * @param b Vector to do the rh multiplication
82  * @param prod OUTPUT vector to receive the result
83  */
84  virtual void mult(const doublereal* b, doublereal* prod) const;
85
86  //! Multiply b*A and write result to prod.
87  /*!
88  * @param b Vector to do the lh multiplication
89  * @param prod OUTPUT vector to receive the result
90  */
91  virtual void leftMult(const doublereal* const b, doublereal* const prod) const;
92
93  /**
94  * Factors the A matrix, overwriting A. We flip m_factored
95  * boolean to indicate that the matrix is now A-1.
96  */
97  int factor();
98
99  //! Factors the A matrix using the QR algorithm, overwriting A
100  /*!
101  * we set m_factored to 2 to indicate the matrix is now QR factored
102  *
103  * @return Returns the info variable from lapack
104  */
105  virtual int factorQR();
106
107  //! Returns an estimate of the inverse of the condition number for the matrix
108  /*!
109  * The matrix must have been previously factored using the QR algorithm
110  *
111  * @return returns the inverse of the condition number
112  */
113  virtual doublereal rcondQR();
114
115  //! Returns an estimate of the inverse of the condition number for the matrix
116  /*!
117  * The matrix must have been previously factored using the LU algorithm
118  *
119  * @param a1norm Norm of the matrix
120  *
121  * @return returns the inverse of the condition number
122  */
123  virtual doublereal rcond(doublereal a1norm);
124
125  //! Returns the one norm of the matrix
126  virtual doublereal oneNorm() const;
127
128  //! Solves the linear problem Ax=b using the QR algorithm returning x in the b spot
129  /*!
130  * @param b RHS to be solved.
131  */
132  int solveQR(doublereal* b);
133
134
135  //! clear the factored flag
136  virtual void clearFactorFlag();
137
138  //! set the factored flag
139  void setFactorFlag();
140
141  //! Report whether the current matrix has been factored.
142  virtual bool factored() const;
143
144  //! Change the way the matrix is factored
145  /*!
146  * @param fAlgorithm integer
147  * 0 LU factorization
148  * 1 QR factorization
149  */
150  virtual void useFactorAlgorithm(int fAlgorithm);
151
152  //! Returns the factor algorithm used
153  /*!
154  * 0 LU decomposition
155  * 1 QR decomposition
156  *
157  * This routine will always return 0
158  */
159  virtual int factorAlgorithm() const;
160
161  //! Return a pointer to the top of column j, columns are assumed to be contiguous in memory
162  /*!
163  * @param j Value of the column
164  *
165  * @return Returns a pointer to the top of the column
166  */
167  virtual doublereal* ptrColumn(size_t j);
168
169  //! Index into the (i,j) element
170  /*!
171  * @param i row
172  * @param j column
173  *
174  * (note, tried a using directive here, and it didn't seem to work)
175  *
176  * Returns a changeable reference to the matrix entry
177  */
178  virtual doublereal& operator()(size_t i, size_t j) {
179  return Array2D::operator()(i, j);
180  }
181
182  //! Copy the data from one array into another without doing any checking
183  /*!
184  * This differs from the assignment operator as no resizing is done and memcpy() is used.
185  * @param y Array to be copied
186  */
187  virtual void copyData(const GeneralMatrix& y);
188
189  //! Constant Index into the (i,j) element
190  /*!
191  * @param i row
192  * @param j column
193  *
194  * Returns an unchangeable reference to the matrix entry
195  */
196  virtual doublereal operator()(size_t i, size_t j) const {
197  return Array2D::operator()(i, j);
198  }
199
200  //! Return the number of rows in the matrix
201  virtual size_t nRows() const;
202
203  //! Return the size and structure of the matrix
204  /*!
205  * This is inherited from GeneralMatrix
206  *
207  * @param iStruct OUTPUT Pointer to a vector of ints that describe the structure of the matrix.
208  * not used
209  *
210  * @return returns the number of rows and columns in the matrix.
211  */
212  size_t nRowsAndStruct(size_t* const iStruct = 0) const;
213
214  //! Duplicate this object
215  virtual GeneralMatrix* duplMyselfAsGeneralMatrix() const;
216
217
218  //! Return an iterator pointing to the first element
219  /*!
220  */
221  virtual vector_fp::iterator begin();
222
223
224  //! Return a const iterator pointing to the first element
225  virtual vector_fp::const_iterator begin() const;
226
227
228  //! Return a vector of const pointers to the columns
229  /*!
230  * Note the value of the pointers are protected by their being const.
231  * However, the value of the matrix is open to being changed.
232  *
233  * @return returns a vector of pointers to the top of the columns
234  * of the matrices.
235  */
236  virtual doublereal* const* colPts();
237
238  //! Check to see if we have any zero rows in the jacobian
239  /*!
240  * This utility routine checks to see if any rows are zero.
241  * The smallest row is returned along with the largest coefficient in that row
242  *
243  * @param valueSmall OUTPUT value of the largest coefficient in the smallest row
244  *
245  * @return index of the row that is most nearly zero
246  */
247  virtual size_t checkRows(doublereal& valueSmall) const;
248
249  //! Check to see if we have any zero columns in the jacobian
250  /*!
251  * This utility routine checks to see if any columns are zero.
252  * The smallest column is returned along with the largest coefficient in that column
253  *
254  * @param valueSmall OUTPUT value of the largest coefficient in the smallest column
255  *
256  * @return index of the column that is most nearly zero
257  */
258  virtual size_t checkColumns(doublereal& valueSmall) const;
259
260 protected:
261
262  //! the factor flag
264
265 public:
266  //! Work vector for QR algorithm
268
269  //! Work vector for QR algorithm
271
272  //! Integer work vector for QR algorithms
273  std::vector<int> iwork_;
274 protected:
275  //! 1-norm of the matrix. This is determined immediately before every factorization
276  doublereal a1norm_;
277
278  //! Use the QR algorithm to factor and invert the matrix
279  int useQR_;
280 };
281 }
282
283 #endif
284