Cantera  2.0
BandMatrix.h
Go to the documentation of this file.
1 /**
2  * @file BandMatrix.h
3  * Declarations for the class BandMatrix
4  * which is a child class of GeneralMatrix for banded matrices handled by solvers
5  * (see class \ref numerics and \link Cantera::BandMatrix BandMatrix\endlink).
6  */
7 
8 // Copyright 2001 California Institute of Technology
9 
10 
11 #ifndef CT_BANDMATRIX_H
12 #define CT_BANDMATRIX_H
13 
14 #include "cantera/base/ct_defs.h"
15 #include "ctlapack.h"
16 #include "cantera/base/utilities.h"
18 #include "GeneralMatrix.h"
19 
20 namespace Cantera
21 {
22 
23 //! A class for banded matrices, involving matrix inversion processes.
24 //! The class is based upon the LAPACK banded storage matrix format.
25 /*!
26  * An important issue with this class is that it stores both the original data
27  * and the LU factorization of the data. This means that the banded matrix typically
28  * will take up twice the room that it is expected to take.
29  *
30  * QR factorizations of banded matrices are not included in the original LAPACK work.
31  * Add-ons are available. However, they are not included here. Instead we just use the
32  * stock LU decompositions.
33  *
34  * This class is a derived class of the base class GeneralMatrix. However, within
35  * the oneD directory, the class is used as is, without reference to the GeneralMatrix
36  * base type.
37  */
38 class BandMatrix : public GeneralMatrix
39 {
40 
41 public:
42 
43  //! Base Constructor
44  /*!
45  * * Create an \c 0 by \c 0 matrix, and initialize all elements to \c 0.
46  */
47  BandMatrix();
48 
49  //! Creates a banded matrix and sets all elements to zero
50  /*!
51  * Create an \c n by \c n banded matrix, and initialize all elements to \c v.
52  *
53  * @param n size of the square matrix
54  * @param kl band size on the lower portion of the matrix
55  * @param ku band size on the upper portion of the matrix
56  * @param v initial value of all matrix components.
57  */
58  BandMatrix(size_t n, size_t kl, size_t ku, doublereal v = 0.0);
59 
60  //! Copy constructor
61  /*!
62  * @param y Matrix to be copied
63  */
64  BandMatrix(const BandMatrix& y);
65 
66  //! Destructor. Does nothing.
67  virtual ~BandMatrix();
68 
69  //! assignment operator
70  /*!
71  * @param y reference to the matrix to be copied
72  */
73  BandMatrix& operator=(const BandMatrix& y);
74 
75  //! Resize the matrix problem
76  /*!
77  * All data is lost
78  *
79  * @param n size of the square matrix
80  * @param kl band size on the lower portion of the matrix
81  * @param ku band size on the upper portion of the matrix
82  * @param v initial value of all matrix components.
83  */
84  void resize(size_t n, size_t kl, size_t ku, doublereal v = 0.0);
85 
86  //! Fill or zero the matrix
87  /*!
88  * @param v Fill value, defaults to zero.
89  */
90  void bfill(doublereal v = 0.0);
91 
92  //! Index into the (i,j) element
93  /*!
94  * @param i row
95  * @param j column
96  *
97  * Returns a changeable reference to the matrix entry
98  */
99  doublereal& operator()(size_t i, size_t j);
100 
101 
102  //! Constant index into the (i,j) element
103  /*!
104  * @param i row
105  * @param j column
106  *
107  * Returns an unchangeable reference to the matrix entry
108  */
109  doublereal operator()(size_t i, size_t j) const;
110 
111  //! Return a changeable reference to element (i,j).
112  /*!
113  * Since this method may alter the element value, it may need to be refactored, so
114  * the flag m_factored is set to false.
115  *
116  * @param i row
117  * @param j column
118  *
119  * @return Returns a reference to the value of the matrix entry
120  */
121  doublereal& value(size_t i, size_t j);
122 
123 
124  //! Return the value of element (i,j).
125  /*!
126  * This method does not alter the array.
127  * @param i row
128  * @param j column
129  *
130  * @return Returns the value of the matrix entry
131  */
132  doublereal value(size_t i, size_t j) const;
133 
134  //! Returns the location in the internal 1D array corresponding to the (i,j) element in the banded array
135  /*!
136  * @param i row
137  * @param j column
138  *
139  * @return Returns the index of the matrix entry
140  */
141  size_t index(size_t i, size_t j) const;
142 
143  //! Return the value of the (i,j) element for (i,j) within the bandwidth.
144  /*!
145  * For efficiency, this method does not check that (i,j) are within the bandwidth; it is up to the calling
146  * program to insure that this is true.
147  *
148  * @param i row
149  * @param j column
150  *
151  * @return Returns the value of the matrix entry
152  */
153  doublereal _value(size_t i, size_t j) const;
154 
155  //! Returns the number of rows
156  virtual size_t nRows() const;
157 
158  //! Return the size and structure of the matrix
159  /*!
160  * This is inherited from GeneralMatrix
161  *
162  * @param iStruct OUTPUT Pointer to a vector of ints that describe the structure of the matrix.
163  * istruct[0] = kl
164  * istruct[1] = ku
165  *
166  * @return returns the number of rows and columns in the matrix.
167  */
168  virtual size_t nRowsAndStruct(size_t* const iStruct = 0) const;
169 
170  //! Number of columns
171  size_t nColumns() const;
172 
173  //! Number of subdiagonals
174  size_t nSubDiagonals() const;
175 
176  //! Number of superdiagonals
177  size_t nSuperDiagonals() const;
178 
179  //! Return the number of rows of storage needed for the band storage
180  size_t ldim() const;
181 
182  //! Return a reference to the pivot vector
183  /*!
184  * @return return a reference to the pivot vector
185  */
186  vector_int& ipiv();
187 
188  //! Multiply A*b and write result to prod.
189  /*!
190  * @param b Vector to do the rh multiplication
191  * @param prod OUTPUT vector to receive the result
192  */
193  virtual void mult(const doublereal* b, doublereal* prod) const;
194 
195  //! Multiply b*A and write result to prod.
196  /*!
197  * @param b Vector to do the lh multiplication
198  * @param prod OUTPUT vector to receive the result
199  */
200  virtual void leftMult(const doublereal* const b, doublereal* const prod) const;
201 
202  //! Perform an LU decomposition, the LAPACK routine DGBTRF is used.
203  /*!
204  *
205  * The factorization is saved in ludata.
206  *
207  * @return Return a success flag.
208  * 0 indicates a success
209  * ~0 Some error occurred, see the LAPACK documentation
210  */
211  int factor();
212 
213 
214  //! Solve the matrix problem Ax = b
215  /*!
216  * @param b INPUT rhs of the problem
217  * @param x OUTPUT solution to the problem
218  *
219  * @return Return a success flag
220  * 0 indicates a success
221  * ~0 Some error occurred, see the LAPACK documentation
222  */
223  int solve(const doublereal* const b, doublereal* const x);
224 
225  //! Solve the matrix problem Ax = b
226  /*!
227  * @param b INPUT rhs of the problem
228  * OUTPUT solution to the problem
229  *
230  * @return Return a success flag
231  * 0 indicates a success
232  * ~0 Some error occurred, see the LAPACK documentation
233  */
234  int solve(doublereal* b);
235 
236 
237  //! Returns an iterator for the start of the band storage data
238  /*!
239  * Iterator points to the beginning of the data, and it is changeable.
240  */
241  virtual vector_fp::iterator begin();
242 
243  //! Returns an iterator for the end of the band storage data
244  /*!
245  * Iterator points to the end of the data, and it is changeable.
246  */
247  vector_fp::iterator end();
248 
249  //! Returns a const iterator for the start of the band storage data
250  /*!
251  * Iterator points to the beginning of the data, and it is not changeable.
252  */
253  vector_fp::const_iterator begin() const;
254 
255  //! Returns a const iterator for the end of the band storage data
256  /*!
257  * Iterator points to the end of the data, and it is not changeable.
258  */
259  vector_fp::const_iterator end() const;
260 
261  /**
262  * Zero the matrix
263  */
264  virtual void zero();
265 
266  //! Factors the A matrix using the QR algorithm, overwriting A
267  /*!
268  * we set m_factored to 2 to indicate the matrix is now QR factored
269  *
270  * @return Returns the info variable from lapack
271  */
272  virtual int factorQR();
273 
274  //! Returns an estimate of the inverse of the condition number for the matrix
275  /*!
276  * The matrix must have been previously factored using the QR algorithm
277  *
278  * @return returns the inverse of the condition number
279  */
280  virtual doublereal rcondQR();
281 
282  //! Returns an estimate of the inverse of the condition number for the matrix
283  /*!
284  * The matrix must have been previously factored using the LU algorithm
285  *
286  * @param a1norm Norm of the matrix
287  *
288  * @return returns the inverse of the condition number
289  */
290  virtual doublereal rcond(doublereal a1norm);
291 
292  //! Change the way the matrix is factored
293  /*!
294  * @param fAlgorithm integer
295  * 0 LU factorization
296  * 1 QR factorization
297  */
298  virtual void useFactorAlgorithm(int fAlgorithm);
299 
300  //! Returns the factor algorithm used
301  /*!
302  * 0 LU decomposition
303  * 1 QR decomposition
304  *
305  * This routine will always return 0
306  */
307  virtual int factorAlgorithm() const;
308 
309  //! Returns the one norm of the matrix
310  virtual doublereal oneNorm() const;
311 
312  //! Duplicate this object as a GeneralMatrix pointer
313  virtual GeneralMatrix* duplMyselfAsGeneralMatrix() const;
314 
315  //! Report whether the current matrix has been factored.
316  virtual bool factored() const;
317 
318  //! Return a pointer to the top of column j, column values are assumed to be contiguous in memory
319  /*!
320  * The LAPACK bandstructure has column values which are contiguous in memory:
321  *
322  * On entry, the matrix A in band storage, in rows KL+1 to
323  * 2*KL+KU+1; rows 1 to KL of the array need not be set.
324  * The j-th column of A is stored in the j-th column of the
325  * array AB as follows:
326  * AB(KL + KU + 1 + i - j,j) = A(i,j) for max(1, j - KU) <= i <= min(m, j + KL)
327  *
328  * This routine returns the position of AB(1,j) (fortran-1 indexing) in the above format
329  *
330  * So to address the (i,j) position, you use the following indexing:
331  *
332  * double *colP_j = matrix.ptrColumn(j);
333  * double a_i_j = colP_j[kl + ku + i - j];
334  *
335  *
336  * @param j Value of the column
337  *
338  * @return Returns a pointer to the top of the column
339  */
340  virtual doublereal* ptrColumn(size_t j);
341 
342  //! Return a vector of const pointers to the columns
343  /*!
344  * Note the value of the pointers are protected by their being const.
345  * However, the value of the matrix is open to being changed.
346  *
347  * @return returns a vector of pointers to the top of the columns
348  * of the matrices.
349  */
350  virtual doublereal* const* colPts();
351 
352  //! Copy the data from one array into another without doing any checking
353  /*!
354  * This differs from the assignment operator as no resizing is done and memcpy() is used.
355  * @param y Array to be copied
356  */
357  virtual void copyData(const GeneralMatrix& y);
358 
359 
360  //! Clear the factored flag
361  virtual void clearFactorFlag();
362 
363  //! Check to see if we have any zero rows in the jacobian
364  /*!
365  * This utility routine checks to see if any rows are zero.
366  * The smallest row is returned along with the largest coefficient in that row
367  *
368  * @param valueSmall OUTPUT value of the largest coefficient in the smallest row
369  *
370  * @return index of the row that is most nearly zero
371  */
372  virtual size_t checkRows(doublereal& valueSmall) const;
373 
374  //! Check to see if we have any zero columns in the jacobian
375  /*!
376  * This utility routine checks to see if any columns are zero.
377  * The smallest column is returned along with the largest coefficient in that column
378  *
379  * @param valueSmall OUTPUT value of the largest coefficient in the smallest column
380  *
381  * @return index of the column that is most nearly zero
382  */
383  virtual size_t checkColumns(doublereal& valueSmall) const;
384 
385 protected:
386 
387  //! Matrix data
389 
390  //! Factorized data
392 
393  //! Boolean indicating whether the matrix is factored
395 
396  //! Number of rows and columns of the matrix
397  size_t m_n;
398 
399  //! Number of subdiagonals of the matrix
400  size_t m_kl;
401 
402  //! Number of super diagonals of the matrix
403  size_t m_ku;
404 
405  //! value of zero
406  doublereal m_zero;
407 
408  //! Pivot vector
410 
411  //! Vector of column pointers
412  std::vector<doublereal*> m_colPtrs;
413 
414  //! Extra work array needed - size = n
416 
417  //! Extra dp work array needed - size = 3n
419 
420 private:
421 
422  //! Error function that gets called for unhandled cases
423  /*!
424  * @param msg String containing the message.
425  */
426  void err(std::string msg) const;
427 
428 
429 };
430 
431 //! Utility routine to print out the matrix
432 /*!
433  * @param s ostream to print the matrix out to
434  * @param m Matrix to be printed
435  *
436  * @return Returns a reference to the ostream
437  */
438 std::ostream& operator<<(std::ostream& s, const BandMatrix& m);
439 
440 }
441 
442 #endif