From 4ccbaed86234d261fad3b3c890f46e4b82a4a768 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Thu, 16 Sep 2010 20:23:56 +0000 Subject: [PATCH] Doxygen update for DenseMatrix --- Cantera/src/numerics/DAE_Solver.h | 5 - Cantera/src/numerics/DenseMatrix.cpp | 32 ++++-- Cantera/src/numerics/DenseMatrix.h | 151 +++++++++++++++++++++------ 3 files changed, 142 insertions(+), 46 deletions(-) diff --git a/Cantera/src/numerics/DAE_Solver.h b/Cantera/src/numerics/DAE_Solver.h index 4d570cca7..f57d28bbe 100644 --- a/Cantera/src/numerics/DAE_Solver.h +++ b/Cantera/src/numerics/DAE_Solver.h @@ -27,11 +27,6 @@ namespace Cantera { #ifdef DAE_DEVEL - /** - * @defgroup numerics Numerical Utilities within Cantera - * - * - */ class Jacobian { public: diff --git a/Cantera/src/numerics/DenseMatrix.cpp b/Cantera/src/numerics/DenseMatrix.cpp index 4cc2c46cf..ae5376a77 100644 --- a/Cantera/src/numerics/DenseMatrix.cpp +++ b/Cantera/src/numerics/DenseMatrix.cpp @@ -67,7 +67,7 @@ namespace Cantera { static_cast(nRows()), b, 1, 0.0, prod, 1); } //==================================================================================================================== - void DenseMatrix::leftMult(const double* b, double* prod) const { + void DenseMatrix::leftMult(const double* const b, double* const prod) const { int nc = static_cast(nColumns()); int nr = static_cast(nRows()); int n, i; @@ -81,14 +81,30 @@ namespace Cantera { } } //==================================================================================================================== + vector_int& DenseMatrix::ipiv() { + return m_ipiv; + } + //==================================================================================================================== int solve(DenseMatrix& A, double* b) { + if (A.nColumns() != A.nRows()) { + throw CanteraError("DenseMatrix::solve", "Can only solve a square matrix"); + } int info=0; ct_dgetrf(static_cast(A.nRows()), static_cast(A.nColumns()), A.ptrColumn(0), //begin(), static_cast(A.nRows()), &A.ipiv()[0], info); - if (info != 0) - throw CanteraError("DenseMatrix::solve", - "DGETRF returned INFO = "+int2str(info)); + if (info != 0) { + if (info > 0) { + throw CanteraError("DenseMatrix::solve", + "DGETRF returned INFO = "+int2str(info) + ". U(i,i) is exactly zero. The factorization has" + " been completed, but the factor U is exactly singular, and division by zero will occur if " + "it is used to solve a system of equations."); + + } else { + throw CanteraError("DenseMatrix::solve", + "DGETRF returned INFO = "+int2str(info) + ". The argument i has an illegal value"); + } + } ct_dgetrs(ctlapack::NoTranspose, static_cast(A.nRows()), 1, A.ptrColumn(0), //begin(), static_cast(A.nRows()), @@ -101,6 +117,9 @@ namespace Cantera { } //==================================================================================================================== int solve(DenseMatrix& A, DenseMatrix& b) { + if (A.nColumns() != A.nRows()) { + throw CanteraError("DenseMatrix::solve", "Can only solve a square matrix"); + } int info=0; ct_dgetrf(static_cast(A.nRows()), static_cast(A.nColumns()), A.ptrColumn(0), @@ -143,14 +162,13 @@ namespace Cantera { } #endif //==================================================================================================================== - void multiply(const DenseMatrix& A, const double* b, double* prod) { + void multiply(const DenseMatrix& A, const double * const b, double * const prod) { ct_dgemv(ctlapack::ColMajor, ctlapack::NoTranspose, static_cast(A.nRows()), static_cast(A.nColumns()), 1.0, A.ptrColumn(0), static_cast(A.nRows()), b, 1, 0.0, prod, 1); } //==================================================================================================================== - void increment(const DenseMatrix& A, - const double* b, double* prod) { + void increment(const DenseMatrix& A, const double* b, double* prod) { ct_dgemv(ctlapack::ColMajor, ctlapack::NoTranspose, static_cast(A.nRows()), static_cast(A.nRows()), 1.0, A.ptrColumn(0), static_cast(A.nRows()), b, 1, 1.0, prod, 1); diff --git a/Cantera/src/numerics/DenseMatrix.h b/Cantera/src/numerics/DenseMatrix.h index 0686eca25..6d43eaef4 100644 --- a/Cantera/src/numerics/DenseMatrix.h +++ b/Cantera/src/numerics/DenseMatrix.h @@ -1,7 +1,8 @@ /** * @file DenseMatrix.h - * - * Dense (not sparse) matrices. + * Headers for the %DenseMatrix object, which deals with dense rectangular matrices and + * description of the numerics groupings of objects + * (see \ref numerics and \link Cantera::DenseMatrix DenseMatrix \endlink) . */ /* @@ -21,10 +22,23 @@ #include "Array.h" namespace Cantera { - /** - * A class for full (non-sparse) matrices with Fortran-compatible - * data storage. Adds matrix operations to class Array2D. + * @defgroup numerics Numerical Utilities within Cantera + * + * Cantera contains some capabilities for solving nonlinear equations and + * integrating both ODE and DAE equation systems in time. This section describes these + * capabilities. + * + */ + + + //! A class for full (non-sparse) matrices with Fortran-compatible + //! data storage, which adds matrix operations to class Array2D. + /*! + * The dense matrix class adds matrix operations onto the Array2D class. + * These matrix operations are carried out by the appropriate BLAS and LAPACK routines + * + * @ingroup numerics */ class DenseMatrix : public Array2D { @@ -32,20 +46,24 @@ namespace Cantera { //! Default Constructor DenseMatrix(); - - /** - * Constructor. Create an \c n by \c m matrix, and initialize - * all elements to \c v. + + //! Constructor. + /*! + * Create an \c n by \c m matrix, and initialize all elements to \c v. + * + * @param n New number of rows + * @param m New number of columns + * @param v Default fill value. defaults to zero. */ DenseMatrix(int n, int m, doublereal v = 0.0); - //! copy constructor + //! Copy constructor /*! * @param y Object to be copied */ DenseMatrix(const DenseMatrix& y); - //! assignment operator + //! Assignment operator /*! * @param y Object to be copied */ @@ -53,55 +71,120 @@ namespace Cantera { //! Destructor. Does nothing. virtual ~DenseMatrix(); - + //! Resize the matrix + /*! + * Resize the matrix to n rows by m cols. + * + * @param n New number of rows + * @param m New number of columns + * @param v Default fill value. defaults to zero. + */ void resize(int n, int m, doublereal v = 0.0); - - - /** - * Multiply A*b and write result to \c prod. + //! Multiply A*b and write result to \c prod. + /*! + * + * @param b input vector b with length N + * @param prod output output vector prod length = M */ virtual void mult(const double* b, double* prod) const; - /** - * Left-multiply the matrix by transpose(b), and write the - * result to prod. + //! Left-multiply the matrix by transpose(b), and write the result to prod. + /*! + * @param b left multiply by this vector. The length must be equal to n + * the number of rows in the matrix. + * + * @param prod Resulting vector. This is of length m, the number of columns + * in the matrix */ - virtual void leftMult(const double* b, double* prod) const; + virtual void leftMult(const double * const b, double* const prod) const; - vector_int& ipiv() { return m_ipiv; } + //! Return a changeable value of the pivot vector + /*! + * @return Returns a reference to the pivot vector as a vector_int + */ + vector_int& ipiv(); + + //! Return a changeable value of the pivot vector + /*! + * @return Returns a reference to the pivot vector as a vector_int + */ const vector_int& ipiv() const { return m_ipiv; } protected: + //! Vector of pivots. Length is equal to the max of m and n. vector_int m_ipiv; }; + //================================================================================================================== - /** - * Solve Ax = b. Array b is overwritten on exit with x. + + //! Solve Ax = b. Array b is overwritten on exit with x. + /*! + * The solve class uses the LAPACK routine dgetrf to invert the m xy n matrix. + * + * The factorization has the form + * A = P * L * U + * where P is a permutation matrix, L is lower triangular with unit + * diagonal elements (lower trapezoidal if m > n), and U is upper + * triangular (upper trapezoidal if m < n). + * + * The system is then solved using the LAPACK routine dgetrs + * + * @param A Dense matrix to be factored + * @param b rhs to be solved. */ int solve(DenseMatrix& A, double* b); - /** Solve Ax = b for multiple right-hand-side vectors. */ + //! Solve Ax = b for multiple right-hand-side vectors. + /*! + * @param A Dense matrix to be factored + * @param b Dense matrix of rhs's. Each column is a rhs + */ int solve(DenseMatrix& A, DenseMatrix& b); #ifdef INCL_LEAST_SQUARES - /** @todo fix lwork */ + //! Solve Ax = b in the least squares sense + /*! + * @param A Matrix to be inverted in the least squares sense + * @param b Vector b to be solved for + * @todo fix lwork + */ int leastSquares(DenseMatrix& A, double* b); #endif - /** - * Multiply \c A*b and return the result in \c prod. Uses BLAS - * routine DGEMV. - */ - void multiply(const DenseMatrix& A, const double* b, double* prod); + + //! Multiply \c A*b and return the result in \c prod. Uses BLAS routine DGEMV. + /*! + * \f[ + * prod_i = sum^N_{j = 1}{A_{ij} b_j} + * \f] + * + * @param A input Dense Matrix A with M rows and N columns + * @param b input vector b with length N + * @param prod output output vector prod length = M + */ + void multiply(const DenseMatrix& A, const double * const b, double * const prod); - void increment(const DenseMatrix& A, - const double* b, double* prod); + //! Multiply \c A*b and add it to the result in \c prod. Uses BLAS routine DGEMV. + /*! + * \f[ + * prod_i += sum^N_{j = 1}{A_{ij} b_j} + * \f] + * + * @param A input Dense Matrix A with M rows and N columns + * @param b input vector b with length N + * @param prod output output vector prod length = M + */ + void increment(const DenseMatrix& A, const double * const b, double * const prod); - /** - * invert A. A is overwritten with A^-1. + //! invert A. A is overwritten with A^-1. + /*! + * @param A Invert the matrix A and store it back in place + * + * @param nn Size of A. This defaults to -1, which means that the number + * of rows is used as the default size of n */ int invert(DenseMatrix& A, int nn=-1);