Cantera  2.4.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 // This file is part of Cantera. See License.txt in the top-level directory or
9 // at http://www.cantera.org/license.txt for license and copyright information.
10 
11 #ifndef CT_BANDMATRIX_H
12 #define CT_BANDMATRIX_H
13 
14 #include "GeneralMatrix.h"
15 
16 namespace Cantera
17 {
18 
19 //! A class for banded matrices, involving matrix inversion processes.
20 //! The class is based upon the LAPACK banded storage matrix format.
21 /*!
22  * An important issue with this class is that it stores both the original data
23  * and the LU factorization of the data. This means that the banded matrix
24  * typically will take up twice the room that it is expected to take.
25  *
26  * QR factorizations of banded matrices are not included in the original LAPACK
27  * work. Add-ons are available. However, they are not included here. Instead we
28  * just use the stock LU decompositions.
29  *
30  * This class is a derived class of the base class GeneralMatrix. However,
31  * within the oneD directory, the class is used as is, without reference to the
32  * GeneralMatrix base type.
33  */
34 class BandMatrix : public GeneralMatrix
35 {
36 public:
37  //! Base Constructor
38  /*!
39  * Create an \c 0 by \c 0 matrix, and initialize all elements to \c 0.
40  */
41  BandMatrix();
42  ~BandMatrix();
43 
44  //! Creates a banded matrix and sets all elements to zero
45  /*!
46  * Create an \c n by \c n banded matrix, and initialize all elements to \c v.
47  *
48  * @param n size of the square matrix
49  * @param kl band size on the lower portion of the matrix
50  * @param ku band size on the upper portion of the matrix
51  * @param v initial value of all matrix components.
52  */
53  BandMatrix(size_t n, size_t kl, size_t ku, doublereal v = 0.0);
54 
55  BandMatrix(const BandMatrix& y);
56  BandMatrix& operator=(const BandMatrix& y);
57 
58  //! Resize the matrix problem
59  /*!
60  * All data is lost
61  *
62  * @param n size of the square matrix
63  * @param kl band size on the lower portion of the matrix
64  * @param ku band size on the upper portion of the matrix
65  * @param v initial value of all matrix components.
66  */
67  void resize(size_t n, size_t kl, size_t ku, doublereal v = 0.0);
68 
69  //! Fill or zero the matrix
70  /*!
71  * @param v Fill value, defaults to zero.
72  */
73  void bfill(doublereal v = 0.0);
74 
75  doublereal& operator()(size_t i, size_t j);
76  doublereal operator()(size_t i, size_t j) const;
77 
78  //! Return a changeable reference to element (i,j).
79  /*!
80  * Since this method may alter the element value, it may need to be
81  * refactored, so the flag m_factored is set to false.
82  *
83  * @param i row
84  * @param j column
85  * @returns a reference to the value of the matrix entry
86  */
87  doublereal& value(size_t i, size_t j);
88 
89  //! Return the value of element (i,j).
90  /*!
91  * This method does not alter the array.
92  * @param i row
93  * @param j column
94  * @returns the value of the matrix entry
95  */
96  doublereal value(size_t i, size_t j) const;
97 
98  //! Returns the location in the internal 1D array corresponding to the (i,j)
99  //! element in the banded array
100  /*!
101  * @param i row
102  * @param j column
103  * @returns the index of the matrix entry
104  */
105  size_t index(size_t i, size_t j) const;
106 
107  //! Return the value of the (i,j) element for (i,j) within the bandwidth.
108  /*!
109  * For efficiency, this method does not check that (i,j) are within the
110  * bandwidth; it is up to the calling program to insure that this is true.
111  *
112  * @param i row
113  * @param j column
114  * @returns the value of the matrix entry
115  */
116  doublereal _value(size_t i, size_t j) const;
117 
118  virtual size_t nRows() const;
119 
120  //! Number of columns
121  size_t nColumns() const;
122 
123  //! Number of subdiagonals
124  size_t nSubDiagonals() const;
125 
126  //! Number of superdiagonals
127  size_t nSuperDiagonals() const;
128 
129  //! Return the number of rows of storage needed for the band storage
130  size_t ldim() const;
131 
132  //! Multiply A*b and write result to \c prod.
133  virtual void mult(const doublereal* b, doublereal* prod) const;
134  virtual void leftMult(const doublereal* const b, doublereal* const prod) const;
135 
136  //! Perform an LU decomposition, the LAPACK routine DGBTRF is used.
137  /*!
138  * The factorization is saved in ludata.
139  *
140  * @returns a success flag. 0 indicates a success; ~0 indicates some
141  * error occurred, see the LAPACK documentation
142  */
143  int factor();
144 
145  //! Solve the matrix problem Ax = b
146  /*!
147  * @param b INPUT RHS of the problem
148  * @param x OUTPUT solution to the problem
149  * @return a success flag. 0 indicates a success; ~0 indicates some error
150  * occurred, see the LAPACK documentation
151  */
152  int solve(const doublereal* const b, doublereal* const x);
153 
154  //! Solve the matrix problem Ax = b
155  /*!
156  * @param b INPUT RHS of the problem
157  * OUTPUT solution to the problem
158  * @param nrhs Number of right hand sides to solve
159  * @param ldb Leading dimension of `b`. Default is nColumns()
160  * @returns a success flag. 0 indicates a success; ~0 indicates some error
161  * occurred, see the LAPACK documentation
162  */
163  int solve(doublereal* b, size_t nrhs=1, size_t ldb=0);
164 
165  //! Returns an iterator for the start of the band storage data
166  /*!
167  * Iterator points to the beginning of the data, and it is changeable.
168  */
169  virtual vector_fp::iterator begin();
170 
171  //! Returns an iterator for the end of the band storage data
172  /*!
173  * Iterator points to the end of the data, and it is changeable.
174  */
175  vector_fp::iterator end();
176 
177  //! Returns a const iterator for the start of the band storage data
178  /*!
179  * Iterator points to the beginning of the data, and it is not changeable.
180  */
181  vector_fp::const_iterator begin() const;
182 
183  //! Returns a const iterator for the end of the band storage data
184  /*!
185  * Iterator points to the end of the data, and it is not changeable.
186  */
187  vector_fp::const_iterator end() const;
188 
189  virtual void zero();
190 
191  //! Returns an estimate of the inverse of the condition number for the matrix
192  /*!
193  * The matrix must have been previously factored using the LU algorithm
194  *
195  * @param a1norm Norm of the matrix
196  * @returns the inverse of the condition number
197  */
198  virtual doublereal rcond(doublereal a1norm);
199 
200  //! Returns the factor algorithm used. This method will always return 0
201  //! (LU) for band matrices.
202  virtual int factorAlgorithm() const;
203 
204  //! Returns the one norm of the matrix
205  virtual doublereal oneNorm() const;
206 
207  //! Return a pointer to the top of column j
208  /*!
209  * Column values are assumed to be contiguous in memory (LAPACK band matrix
210  * structure)
211  *
212  * On entry, the matrix A in band storage, in rows KL+1 to 2*KL+KU+1; rows 1
213  * to KL of the array need not be set. The j-th column of A is stored in the
214  * j-th column of the array AB as follows:
215  *
216  * AB(KL + KU + 1 + i - j,j) = A(i,j) for max(1, j - KU) <= i <= min(m, j + KL)
217  *
218  * This routine returns the position of AB(1,j) (fortran-1 indexing) in the
219  * above format
220  *
221  * So to address the (i,j) position, you use the following indexing:
222  *
223  * double *colP_j = matrix.ptrColumn(j);
224  * double a_i_j = colP_j[kl + ku + i - j];
225  *
226  * @param j Value of the column
227  * @returns a pointer to the top of the column
228  */
229  virtual doublereal* ptrColumn(size_t j);
230 
231  //! Return a vector of const pointers to the columns
232  /*!
233  * Note the value of the pointers are protected by their being const.
234  * However, the value of the matrix is open to being changed.
235  *
236  * @returns a vector of pointers to the top of the columns of the matrices.
237  */
238  virtual doublereal* const* colPts();
239 
240  //! Check to see if we have any zero rows in the Jacobian
241  /*!
242  * This utility routine checks to see if any rows are zero. The smallest row
243  * is returned along with the largest coefficient in that row
244  *
245  * @param valueSmall OUTPUT value of the largest coefficient in the smallest row
246  * @return index of the row that is most nearly zero
247  */
248  virtual size_t checkRows(doublereal& valueSmall) const;
249 
250  //! Check to see if we have any zero columns in the Jacobian
251  /*!
252  * This utility routine checks to see if any columns are zero. The smallest
253  * column is returned along with the largest coefficient in that column
254  *
255  * @param valueSmall OUTPUT value of the largest coefficient in the smallest column
256  * @return index of the column that is most nearly zero
257  */
258  virtual size_t checkColumns(doublereal& valueSmall) const;
259 
260  //! LAPACK "info" flag after last factor/solve operation
261  int info() const { return m_info; };
262 
263 protected:
264  //! Matrix data
265  vector_fp data;
266 
267  //! Factorized data
269 
270  //! Number of rows and columns of the matrix
271  size_t m_n;
272 
273  //! Number of subdiagonals of the matrix
274  size_t m_kl;
275 
276  //! Number of super diagonals of the matrix
277  size_t m_ku;
278 
279  //! value of zero
280  doublereal m_zero;
281 
282  class PivData; // pImpl wrapper class
283 
284  //! Pivot vector
285  std::unique_ptr<PivData> m_ipiv;
286 
287  //! Vector of column pointers
288  std::vector<doublereal*> m_colPtrs;
289  std::vector<double*> m_lu_col_ptrs;
290 
291  //! Extra work array needed - size = n
293 
294  //! Extra dp work array needed - size = 3n
296 
297  int m_info;
298 };
299 
300 //! Utility routine to print out the matrix
301 /*!
302  * @param s ostream to print the matrix out to
303  * @param m Matrix to be printed
304  * @returns a reference to the ostream
305  */
306 std::ostream& operator<<(std::ostream& s, const BandMatrix& m);
307 
308 }
309 
310 #endif
size_t m_kl
Number of subdiagonals of the matrix.
Definition: BandMatrix.h:274
virtual vector_fp::iterator begin()
Returns an iterator for the start of the band storage data.
Definition: BandMatrix.cpp:307
vector_int iwork_
Extra work array needed - size = n.
Definition: BandMatrix.h:292
virtual doublereal * ptrColumn(size_t j)
Return a pointer to the top of column j.
Definition: BandMatrix.cpp:428
vector_fp work_
Extra dp work array needed - size = 3n.
Definition: BandMatrix.h:295
BandMatrix()
Base Constructor.
Definition: BandMatrix.cpp:49
doublereal _value(size_t i, size_t j) const
Return the value of the (i,j) element for (i,j) within the bandwidth.
Definition: BandMatrix.cpp:196
virtual size_t checkColumns(doublereal &valueSmall) const
Check to see if we have any zero columns in the Jacobian.
Definition: BandMatrix.cpp:406
int solve(const doublereal *const b, doublereal *const x)
Solve the matrix problem Ax = b.
Definition: BandMatrix.cpp:273
doublereal & operator()(size_t i, size_t j)
Index into the (i,j) element.
Definition: BandMatrix.cpp:164
int info() const
LAPACK "info" flag after last factor/solve operation.
Definition: BandMatrix.h:261
size_t index(size_t i, size_t j) const
Returns the location in the internal 1D array corresponding to the (i,j) element in the banded array...
Definition: BandMatrix.cpp:191
virtual void leftMult(const doublereal *const b, doublereal *const prod) const
Multiply b*A and write result to prod.
Definition: BandMatrix.cpp:239
Generic matrix.
Definition: GeneralMatrix.h:21
vector_fp data
Matrix data.
Definition: BandMatrix.h:261
size_t m_n
Number of rows and columns of the matrix.
Definition: BandMatrix.h:271
size_t nSubDiagonals() const
Number of subdiagonals.
Definition: BandMatrix.cpp:211
void resize(size_t n, size_t kl, size_t ku, doublereal v=0.0)
Resize the matrix problem.
Definition: BandMatrix.cpp:133
size_t nColumns() const
Number of columns.
Definition: BandMatrix.cpp:206
virtual doublereal *const * colPts()
Return a vector of const pointers to the columns.
Definition: BandMatrix.cpp:433
virtual doublereal oneNorm() const
Returns the one norm of the matrix.
Definition: BandMatrix.cpp:369
virtual int factorAlgorithm() const
Returns the factor algorithm used.
Definition: BandMatrix.cpp:364
std::vector< int > vector_int
Vector of ints.
Definition: ct_defs.h:159
std::unique_ptr< PivData > m_ipiv
Pivot vector.
Definition: BandMatrix.h:282
size_t nSuperDiagonals() const
Number of superdiagonals.
Definition: BandMatrix.cpp:216
virtual doublereal rcond(doublereal a1norm)
Returns an estimate of the inverse of the condition number for the matrix.
Definition: BandMatrix.cpp:341
virtual size_t checkRows(doublereal &valueSmall) const
Check to see if we have any zero rows in the Jacobian.
Definition: BandMatrix.cpp:384
doublereal & value(size_t i, size_t j)
Return a changeable reference to element (i,j).
Definition: BandMatrix.cpp:174
vector_fp ludata
Factorized data.
Definition: BandMatrix.h:268
size_t m_ku
Number of super diagonals of the matrix.
Definition: BandMatrix.h:277
virtual size_t nRows() const
Return the number of rows in the matrix.
Definition: BandMatrix.cpp:201
size_t ldim() const
Return the number of rows of storage needed for the band storage.
Definition: BandMatrix.cpp:221
void bfill(doublereal v=0.0)
Fill or zero the matrix.
Definition: BandMatrix.cpp:152
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:157
doublereal m_zero
value of zero
Definition: BandMatrix.h:280
int factor()
Perform an LU decomposition, the LAPACK routine DGBTRF is used.
Definition: BandMatrix.cpp:252
Declarations for the class GeneralMatrix which is a virtual base class for matrices handled by solver...
virtual void zero()
Zero the matrix elements.
Definition: BandMatrix.cpp:158
std::vector< doublereal * > m_colPtrs
Vector of column pointers.
Definition: BandMatrix.h:288
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:8
virtual void mult(const doublereal *b, doublereal *prod) const
Multiply A*b and write result to prod.
Definition: BandMatrix.cpp:226
vector_fp::iterator end()
Returns an iterator for the end of the band storage data.
Definition: BandMatrix.cpp:313
A class for banded matrices, involving matrix inversion processes.
Definition: BandMatrix.h:34