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