research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Doxygen update for DenseMatrix

Browse source
This commit is contained in:
Harry Moffat 2010-09-16 20:23:56 +00:00
parent 0950fbe256
commit 4ccbaed862
3 changed files with 142 additions and 46 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

5
Cantera/src/numerics/DAE_Solver.h
View file

@ -27,11 +27,6 @@
namespace Cantera {
#ifdef DAE_DEVEL
/**
* @defgroup numerics Numerical Utilities within Cantera
*
*
*/
class Jacobian {
public:

32
Cantera/src/numerics/DenseMatrix.cpp
View file

@ -67,7 +67,7 @@ namespace Cantera {
static_cast<int>(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<int>(nColumns());
int nr = static_cast<int>(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<int>(A.nRows()),
static_cast<int>(A.nColumns()), A.ptrColumn(0), //begin(),
static_cast<int>(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<int>(A.nRows()), 1, A.ptrColumn(0), //begin(),
static_cast<int>(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<int>(A.nRows()),
static_cast<int>(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<int>(A.nRows()), static_cast<int>(A.nColumns()), 1.0,
A.ptrColumn(0), static_cast<int>(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<int>(A.nRows()), static_cast<int>(A.nRows()), 1.0,
A.ptrColumn(0), static_cast<int>(A.nRows()), b, 1, 1.0, prod, 1);

151
Cantera/src/numerics/DenseMatrix.h
View file

@ -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);