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