From a737dd308f154c2f9dcfc4ab39f0b4ab5473b44a Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Wed, 21 Jul 2010 23:35:42 +0000 Subject: [PATCH] Doxygen update to DenseMatrix --- Cantera/src/numerics/DenseMatrix.cpp | 279 +++++++++++++++------------ Cantera/src/numerics/DenseMatrix.h | 129 +++++++------ 2 files changed, 222 insertions(+), 186 deletions(-) diff --git a/Cantera/src/numerics/DenseMatrix.cpp b/Cantera/src/numerics/DenseMatrix.cpp index 762e4e393..4cc2c46cf 100644 --- a/Cantera/src/numerics/DenseMatrix.cpp +++ b/Cantera/src/numerics/DenseMatrix.cpp @@ -16,134 +16,165 @@ #include "stringUtils.h" namespace Cantera { - - /// assignment. - DenseMatrix& DenseMatrix::operator=(const DenseMatrix& y) { - if (&y == this) return *this; - Array2D::operator=(y); - m_ipiv = y.ipiv(); - return *this; - } - - void DenseMatrix::resize(int n, int m, doublereal v) { - Array2D::resize(n,m,v); - m_ipiv.resize( max(n,m) ); - } - - void DenseMatrix::mult(const double* b, double* prod) const { - ct_dgemv(ctlapack::ColMajor, ctlapack::NoTranspose, - static_cast(nRows()), - static_cast(nRows()), 1.0, ptrColumn(0), //begin(), - static_cast(nRows()), b, 1, 0.0, prod, 1); - } - - void DenseMatrix::leftMult(const double* b, double* prod) const { - int nc = static_cast(nColumns()); - int nr = static_cast(nRows()); - int n, i; - double sum = 0.0; - for (n = 0; n < nc; n++) { - sum = 0.0; - for (i = 0; i < nr; i++) { - sum += value(i,n)*b[i]; - } - prod[n] = sum; - } - } - - int solve(DenseMatrix& A, double* b) { - 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)); - ct_dgetrs(ctlapack::NoTranspose, - static_cast(A.nRows()), 1, A.ptrColumn(0), //begin(), - static_cast(A.nRows()), - &A.ipiv()[0], b, - static_cast(A.nColumns()), info); - if (info != 0) - throw CanteraError("DenseMatrix::solve", - "DGETRS returned INFO = "+int2str(info)); - return 0; + //==================================================================================================================== + // Default Constructor + DenseMatrix::DenseMatrix() + { + } + //==================================================================================================================== + /* + * Constructor. Create an \c n by \c m matrix, and initialize + * all elements to \c v. + */ + DenseMatrix::DenseMatrix(int n, int m, doublereal v) : + Array2D(n, m, v) + { + m_ipiv.resize(max(n, m)); + } + //==================================================================================================================== + // copy constructor + /* + * @param y Object to be copied + */ + DenseMatrix::DenseMatrix(const DenseMatrix& y) : + Array2D(y) + { + m_ipiv = y.ipiv(); + } + //==================================================================================================================== + // assignment + DenseMatrix& DenseMatrix::operator=(const DenseMatrix& y) { + if (&y == this) return *this; + Array2D::operator=(y); + m_ipiv = y.ipiv(); + return *this; + } + //==================================================================================================================== + // Destructor. Does nothing. + DenseMatrix::~DenseMatrix() + { + } + //==================================================================================================================== + void DenseMatrix::resize(int n, int m, doublereal v) { + Array2D::resize(n,m,v); + m_ipiv.resize( max(n,m) ); + } + //==================================================================================================================== + void DenseMatrix::mult(const double* b, double* prod) const { + ct_dgemv(ctlapack::ColMajor, ctlapack::NoTranspose, + static_cast(nRows()), + static_cast(nRows()), 1.0, ptrColumn(0), //begin(), + static_cast(nRows()), b, 1, 0.0, prod, 1); + } + //==================================================================================================================== + void DenseMatrix::leftMult(const double* b, double* prod) const { + int nc = static_cast(nColumns()); + int nr = static_cast(nRows()); + int n, i; + double sum = 0.0; + for (n = 0; n < nc; n++) { + sum = 0.0; + for (i = 0; i < nr; i++) { + sum += value(i,n)*b[i]; + } + prod[n] = sum; } - - int solve(DenseMatrix& A, DenseMatrix& b) { - int info=0; - ct_dgetrf(static_cast(A.nRows()), - static_cast(A.nColumns()), A.ptrColumn(0), - static_cast(A.nRows()), &A.ipiv()[0], info); - if (info != 0) - throw CanteraError("DenseMatrix::solve", - "DGETRF returned INFO = "+int2str(info)); - ct_dgetrs(ctlapack::NoTranspose, static_cast(A.nRows()), - static_cast(b.nColumns()), - A.ptrColumn(0), static_cast(A.nRows()), - &A.ipiv()[0], b.ptrColumn(0), - static_cast(b.nRows()), info); - if (info != 0) - throw CanteraError("DenseMatrix::solve", - "DGETRS returned INFO = "+int2str(info)); - return 0; - } - + } + //==================================================================================================================== + int solve(DenseMatrix& A, double* b) { + 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)); + ct_dgetrs(ctlapack::NoTranspose, + static_cast(A.nRows()), 1, A.ptrColumn(0), //begin(), + static_cast(A.nRows()), + &A.ipiv()[0], b, + static_cast(A.nColumns()), info); + if (info != 0) + throw CanteraError("DenseMatrix::solve", + "DGETRS returned INFO = "+int2str(info)); + return 0; + } + //==================================================================================================================== + int solve(DenseMatrix& A, DenseMatrix& b) { + int info=0; + ct_dgetrf(static_cast(A.nRows()), + static_cast(A.nColumns()), A.ptrColumn(0), + static_cast(A.nRows()), &A.ipiv()[0], info); + if (info != 0) + throw CanteraError("DenseMatrix::solve", + "DGETRF returned INFO = "+int2str(info)); + ct_dgetrs(ctlapack::NoTranspose, static_cast(A.nRows()), + static_cast(b.nColumns()), + A.ptrColumn(0), static_cast(A.nRows()), + &A.ipiv()[0], b.ptrColumn(0), + static_cast(b.nRows()), info); + if (info != 0) + throw CanteraError("DenseMatrix::solve", + "DGETRS returned INFO = "+int2str(info)); + return 0; + } + //==================================================================================================================== #ifdef INCL_LEAST_SQUARES - /** @todo fix lwork */ - int leastSquares(DenseMatrix& A, double* b) { - int info = 0; - int rank = 0; - double rcond = -1.0; - // fix this! - int lwork = 6000; // 2*(3*min(m,n) + max(2*min(m,n), max(m,n))); - vector_fp work(lwork); - vector_fp s(min(static_cast(A.nRows()), - static_cast(A.nColumns()))); - ct_dgelss(static_cast(A.nRows()), - static_cast(A.nColumns()), 1, A.ptrColumn(0), - static_cast(A.nRows()), b, - static_cast(A.nColumns()), &s[0], //.begin(), - rcond, rank, &work[0], work.size(), info); - if (info != 0) - throw CanteraError("DenseMatrix::leaseSquares", - "DGELSS returned INFO = "+int2str(info)); - return 0; - } + /** @todo fix lwork */ + int leastSquares(DenseMatrix& A, double* b) { + int info = 0; + int rank = 0; + double rcond = -1.0; + // fix this! + int lwork = 6000; // 2*(3*min(m,n) + max(2*min(m,n), max(m,n))); + vector_fp work(lwork); + vector_fp s(min(static_cast(A.nRows()), + static_cast(A.nColumns()))); + ct_dgelss(static_cast(A.nRows()), + static_cast(A.nColumns()), 1, A.ptrColumn(0), + static_cast(A.nRows()), b, + static_cast(A.nColumns()), &s[0], //.begin(), + rcond, rank, &work[0], work.size(), info); + if (info != 0) + throw CanteraError("DenseMatrix::leaseSquares", + "DGELSS returned INFO = "+int2str(info)); + return 0; + } #endif + //==================================================================================================================== + void multiply(const DenseMatrix& A, const double* b, double* 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) { + 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); + } + //==================================================================================================================== + int invert(DenseMatrix& A, int nn) { + integer n = (nn > 0 ? nn : static_cast(A.nRows())); + int info=0; + ct_dgetrf(n, n, A.ptrColumn(0), static_cast(A.nRows()), + &A.ipiv()[0], info); + if (info != 0) + throw CanteraError("invert", + "DGETRF returned INFO="+int2str(info)); - void multiply(const DenseMatrix& A, const double* b, double* 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) { - 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); - } - - int invert(DenseMatrix& A, int nn) { - integer n = (nn > 0 ? nn : static_cast(A.nRows())); - int info=0; - ct_dgetrf(n, n, A.ptrColumn(0), static_cast(A.nRows()), - &A.ipiv()[0], info); - if (info != 0) - throw CanteraError("invert", - "DGETRF returned INFO="+int2str(info)); - - vector_fp work(n); - integer lwork = static_cast(work.size()); - ct_dgetri(n, A.ptrColumn(0), static_cast(A.nRows()), - &A.ipiv()[0], - &work[0], lwork, info); - if (info != 0) - throw CanteraError("invert", - "DGETRI returned INFO="+int2str(info)); - return 0; - } + vector_fp work(n); + integer lwork = static_cast(work.size()); + ct_dgetri(n, A.ptrColumn(0), static_cast(A.nRows()), + &A.ipiv()[0], + &work[0], lwork, info); + if (info != 0) + throw CanteraError("invert", + "DGETRI returned INFO="+int2str(info)); + return 0; + } + //==================================================================================================================== } diff --git a/Cantera/src/numerics/DenseMatrix.h b/Cantera/src/numerics/DenseMatrix.h index 613bb948b..0686eca25 100644 --- a/Cantera/src/numerics/DenseMatrix.h +++ b/Cantera/src/numerics/DenseMatrix.h @@ -22,83 +22,88 @@ namespace Cantera { - /** - * A class for full (non-sparse) matrices with Fortran-compatible - * data storage. Adds matrix operations to class Array2D. + /** + * A class for full (non-sparse) matrices with Fortran-compatible + * data storage. Adds matrix operations to class Array2D. + */ + class DenseMatrix : public Array2D { + + public: + + //! Default Constructor + DenseMatrix(); + + /** + * Constructor. Create an \c n by \c m matrix, and initialize + * all elements to \c v. */ - class DenseMatrix : public Array2D { + DenseMatrix(int n, int m, doublereal v = 0.0); - public: + //! copy constructor + /*! + * @param y Object to be copied + */ + DenseMatrix(const DenseMatrix& y); - DenseMatrix(){} + //! assignment operator + /*! + * @param y Object to be copied + */ + DenseMatrix& operator=(const DenseMatrix& y); - /** - * Constructor. Create an \c n by \c m matrix, and initialize - * all elements to \c v. - */ - DenseMatrix(int n, int m, doublereal v = 0.0) : Array2D(n,m,v) { - m_ipiv.resize( max(n, m) ); - } + //! Destructor. Does nothing. + virtual ~DenseMatrix(); + - /// copy constructor - DenseMatrix(const DenseMatrix& y) : Array2D(y) { - m_ipiv = y.ipiv(); - } - - /// assignment. - DenseMatrix& operator=(const DenseMatrix& y); - - void resize(int n, int m, doublereal v = 0.0); - - /// Destructor. Does nothing. - virtual ~DenseMatrix(){} - - - /** - * Multiply A*b and write result to \c prod. - */ - virtual void mult(const double* b, double* prod) const; - - /** - * Left-multiply the matrix by transpose(b), and write the - * result to prod. - */ - virtual void leftMult(const double* b, double* prod) const; - - vector_int& ipiv() { return m_ipiv; } - const vector_int& ipiv() const { return m_ipiv; } - - protected: - - vector_int m_ipiv; - }; + void resize(int n, int m, doublereal v = 0.0); + /** - * Solve Ax = b. Array b is overwritten on exit with x. + * Multiply A*b and write result to \c prod. */ - int solve(DenseMatrix& A, double* b); + virtual void mult(const double* b, double* prod) const; - /** Solve Ax = b for multiple right-hand-side vectors. */ - int solve(DenseMatrix& A, DenseMatrix& b); + /** + * Left-multiply the matrix by transpose(b), and write the + * result to prod. + */ + virtual void leftMult(const double* b, double* prod) const; + + vector_int& ipiv() { return m_ipiv; } + const vector_int& ipiv() const { return m_ipiv; } + + protected: + + vector_int m_ipiv; + }; + + + /** + * Solve Ax = b. Array b is overwritten on exit with x. + */ + int solve(DenseMatrix& A, double* b); + + /** Solve Ax = b for multiple right-hand-side vectors. */ + int solve(DenseMatrix& A, DenseMatrix& b); #ifdef INCL_LEAST_SQUARES - /** @todo fix lwork */ - int leastSquares(DenseMatrix& A, double* b); + /** @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. + */ + void multiply(const DenseMatrix& A, const double* b, double* prod); - void increment(const DenseMatrix& A, - const double* b, double* prod); + void increment(const DenseMatrix& A, + const double* b, double* prod); - /** - * invert A. A is overwritten with A^-1. - */ - int invert(DenseMatrix& A, int nn=-1); + /** + * invert A. A is overwritten with A^-1. + */ + int invert(DenseMatrix& A, int nn=-1); }