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