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