Cantera  3.1.0a1
Loading...
Searching...
No Matches
DenseMatrix.h
Go to the documentation of this file.
1/**
2 * @file DenseMatrix.h
3 * Headers for the DenseMatrix object, which deals with dense rectangular matrices and
4 * description of the numerics groupings of objects
5 * (see @ref numerics and @link Cantera::DenseMatrix DenseMatrix @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_DENSEMATRIX_H
12#define CT_DENSEMATRIX_H
13
16#include "cantera/base/Array.h"
17
18namespace Cantera
19{
20/**
21 * @defgroup numerics Numerical Utilities
22 *
23 * @details %Cantera contains some capabilities for solving nonlinear equations and
24 * integrating both ODE and DAE equation systems in time.
25 */
26
27/**
28 * @defgroup matrices Matrix Handling
29 * Classes and methods implementing matrix operations.
30 * @ingroup numerics
31 */
32
33//! A class for full (non-sparse) matrices with Fortran-compatible data storage,
34//! which adds matrix operations to class Array2D.
35/*!
36 * The dense matrix class adds matrix operations onto the Array2D class. These
37 * matrix operations are carried out by the appropriate BLAS and LAPACK routines
38 *
39 * Error handling from BLAS and LAPACK are handled via the following
40 * formulation. Depending on a variable, a singular matrix or other terminal
41 * error condition from LAPACK is handled by either throwing an exception or
42 * by returning the error code condition to the calling routine.
43 *
44 * The int variable, m_useReturnErrorCode, determines which method is used. The
45 * default value of zero means that an exception is thrown. A value of 1 means
46 * that a return code is used.
47 *
48 * Reporting of these LAPACK error conditions is handled by the class variable
49 * m_printLevel. The default is for no reporting. If m_printLevel is nonzero,
50 * the error condition is reported to Cantera's log file.
51 *
52 * @ingroup matrices
53 */
54class DenseMatrix : public Array2D
55{
56public:
57 //! Default Constructor
58 DenseMatrix() = default;
59
60 //! Constructor.
61 /*!
62 * Create an @c n by @c m matrix, and initialize all elements to @c v.
63 *
64 * @param n New number of rows
65 * @param m New number of columns
66 * @param v Default fill value. defaults to zero.
67 */
68 DenseMatrix(size_t n, size_t m, double v = 0.0);
69
70 DenseMatrix(const DenseMatrix& y);
71 DenseMatrix& operator=(const DenseMatrix& y);
72
73 //! Resize the matrix
74 /*!
75 * Resize the matrix to n rows by m cols.
76 *
77 * @param n New number of rows
78 * @param m New number of columns
79 * @param v Default fill value. defaults to zero.
80 */
81 void resize(size_t n, size_t m, double v=0.0) override;
82
83 virtual double* const* colPts();
84
85 //! Return a const vector of const pointers to the columns
86 /*!
87 * Note, the Jacobian can not be altered by this routine, and therefore the
88 * member function is const.
89 *
90 * @returns a vector of pointers to the top of the columns of the matrices.
91 */
92 const double* const* const_colPts() const;
93
94 virtual void mult(const double* b, double* prod) const;
95
96 //! Multiply A*B and write result to @c prod.
97 /*!
98 * Take this matrix to be of size NxM.
99 * @param[in] b DenseMatrix B of size MxP
100 * @param[out] prod DenseMatrix prod size NxP
101 */
102 virtual void mult(const DenseMatrix& b, DenseMatrix& prod) const;
103
104 //! Left-multiply the matrix by transpose(b), and write the result to prod.
105 /*!
106 * @param b left multiply by this vector. The length must be equal to n
107 * the number of rows in the matrix.
108 * @param prod Resulting vector. This is of length m, the number of columns
109 * in the matrix
110 */
111 virtual void leftMult(const double* const b, double* const prod) const;
112
113 //! Return a changeable value of the pivot vector
114 /*!
115 * @returns a reference to the pivot vector as a vector<int>
116 */
117 vector<int>& ipiv();
118
119 //! Return a changeable value of the pivot vector
120 /*!
121 * @returns a reference to the pivot vector as a vector<int>
122 */
123 const vector<int>& ipiv() const {
124 return m_ipiv;
125 }
126
127protected:
128 //! Vector of pivots. Length is equal to the max of m and n.
129 vector<int> m_ipiv;
130
131 //! Vector of column pointers
132 vector<double*> m_colPts;
133
134public:
135 //! Error Handling Flag
136 /*!
137 * The default is to set this to 0. In this case, if a factorization is
138 * requested and can't be achieved, a CESingularMatrix exception is
139 * triggered. No return code is used, because an exception is thrown. If
140 * this is set to 1, then an exception is not thrown. Routines return with
141 * an error code, that is up to the calling routine to handle correctly.
142 * Negative return codes always throw an exception.
143 */
145
146 //! Print Level
147 /*!
148 * Printing is done to the log file using the routine writelogf().
149 *
150 * Level of printing that is carried out. Only error conditions are printed
151 * out, if this value is nonzero.
152 */
154};
155
156
157//! Solve Ax = b. Array b is overwritten on exit with x.
158/*!
159 * The solve function uses the LAPACK routine dgetrf to invert the m xy n matrix.
160 *
161 * The factorization has the form
162 *
163 * A = P * L * U
164 *
165 * where P is a permutation matrix, L is lower triangular with unit diagonal
166 * elements (lower trapezoidal if m > n), and U is upper triangular (upper
167 * trapezoidal if m < n).
168 *
169 * The system is then solved using the LAPACK routine dgetrs
170 *
171 * @param A Dense matrix to be factored
172 * @param b RHS(s) to be solved.
173 * @param nrhs Number of right hand sides to solve
174 * @param ldb Leading dimension of b, if nrhs > 1
175 */
176int solve(DenseMatrix& A, double* b, size_t nrhs=1, size_t ldb=0);
177
178//! Solve Ax = b for multiple right-hand-side vectors.
179/*!
180 * @param A Dense matrix to be factored
181 * @param b Dense matrix of RHS's. Each column is a RHS
182 */
183int solve(DenseMatrix& A, DenseMatrix& b);
184
185//! Multiply @c A*b and return the result in @c prod. Uses BLAS routine DGEMV.
186/*!
187 * @f[
188 * prod_i = sum^N_{j = 1}{A_{ij} b_j}
189 * @f]
190 *
191 * @param[in] A Dense Matrix A with M rows and N columns
192 * @param[in] b vector b with length N
193 * @param[out] prod vector prod length = M
194 */
195void multiply(const DenseMatrix& A, const double* const b, double* const prod);
196
197//! Multiply @c A*b and add it to the result in @c prod. Uses BLAS routine DGEMV.
198/*!
199 * @f[
200 * prod_i += sum^N_{j = 1}{A_{ij} b_j}
201 * @f]
202 *
203 * @param[in] A Dense Matrix A with M rows and N columns
204 * @param[in] b vector b with length N
205 * @param[out] prod vector prod length = M
206 */
207void increment(const DenseMatrix& A, const double* const b, double* const prod);
208
209//! invert A. A is overwritten with A^-1.
210/*!
211 * @param A Invert the matrix A and store it back in place
212 * @param nn Size of A. This defaults to -1, which means that the number of
213 * rows is used as the default size of n
214 */
215int invert(DenseMatrix& A, size_t nn=npos);
216
217}
218
219#endif
Header file for class Cantera::Array2D.
A class for 2D arrays stored in column-major (Fortran-compatible) form.
Definition Array.h:32
A class for full (non-sparse) matrices with Fortran-compatible data storage, which adds matrix operat...
Definition DenseMatrix.h:55
virtual void leftMult(const double *const b, double *const prod) const
Left-multiply the matrix by transpose(b), and write the result to prod.
int m_printLevel
Print Level.
int m_useReturnErrorCode
Error Handling Flag.
const vector< int > & ipiv() const
Return a changeable value of the pivot vector.
void resize(size_t n, size_t m, double v=0.0) override
Resize the matrix.
vector< int > & ipiv()
Return a changeable value of the pivot vector.
DenseMatrix()=default
Default Constructor.
const double *const * const_colPts() const
Return a const vector of const pointers to the columns.
vector< double * > m_colPts
Vector of column pointers.
vector< int > m_ipiv
Vector of pivots. Length is equal to the max of m and n.
This file contains definitions of constants, types and terms that are used in internal routines and a...
Definitions for the classes that are thrown when Cantera experiences an error condition (also contain...
Namespace for the Cantera kernel.
Definition AnyMap.cpp:564
void increment(const DenseMatrix &A, const double *b, double *prod)
Multiply A*b and add it to the result in prod. Uses BLAS routine DGEMV.
const size_t npos
index returned by functions to indicate "no position"
Definition ct_defs.h:180
void multiply(const DenseMatrix &A, const double *const b, double *const prod)
Multiply A*b and return the result in prod. Uses BLAS routine DGEMV.
int invert(DenseMatrix &A, size_t nn)
invert A. A is overwritten with A^-1.
int solve(DenseMatrix &A, double *b, size_t nrhs, size_t ldb)
Solve Ax = b. Array b is overwritten on exit with x.