From e9daef3313cb1758ff0220d62b1be3fb2cc705fd Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 29 May 2013 00:11:22 +0000 Subject: [PATCH] Cleaned up Doxygen Documentation for matrix classes BandMatrix, DenseMatrix, SquareMatrix, and GeneralMatrix --- include/cantera/numerics/BandMatrix.h | 56 +-------- include/cantera/numerics/DenseMatrix.h | 26 ---- include/cantera/numerics/GeneralMatrix.h | 19 +-- include/cantera/numerics/SquareMatrix.h | 144 +---------------------- src/numerics/BandMatrix.cpp | 140 ++++++---------------- src/numerics/DenseMatrix.cpp | 45 +++---- src/numerics/GeneralMatrix.cpp | 10 +- src/numerics/SquareMatrix.cpp | 130 ++++---------------- 8 files changed, 94 insertions(+), 476 deletions(-) diff --git a/include/cantera/numerics/BandMatrix.h b/include/cantera/numerics/BandMatrix.h index f639fb6b9..c3a7af074 100644 --- a/include/cantera/numerics/BandMatrix.h +++ b/include/cantera/numerics/BandMatrix.h @@ -7,7 +7,6 @@ // Copyright 2001 California Institute of Technology - #ifndef CT_BANDMATRIX_H #define CT_BANDMATRIX_H @@ -86,23 +85,7 @@ public: */ void bfill(doublereal v = 0.0); - //! Index into the (i,j) element - /*! - * @param i row - * @param j column - * - * Returns a changeable reference to the matrix entry - */ doublereal& operator()(size_t i, size_t j); - - - //! Constant index into the (i,j) element - /*! - * @param i row - * @param j column - * - * Returns an unchangeable reference to the matrix entry - */ doublereal operator()(size_t i, size_t j) const; //! Return a changeable reference to element (i,j). @@ -117,7 +100,6 @@ public: */ doublereal& value(size_t i, size_t j); - //! Return the value of element (i,j). /*! * This method does not alter the array. @@ -149,13 +131,10 @@ public: */ doublereal _value(size_t i, size_t j) const; - //! Returns the number of rows virtual size_t nRows() const; //! Return the size and structure of the matrix /*! - * This is inherited from GeneralMatrix - * * @param iStruct OUTPUT Pointer to a vector of ints that describe the structure of the matrix. * istruct[0] = kl * istruct[1] = ku @@ -177,28 +156,13 @@ public: size_t ldim() const; //! Return a reference to the pivot vector - /*! - * @return return a reference to the pivot vector - */ vector_int& ipiv(); - //! Multiply A*b and write result to prod. - /*! - * @param b Vector to do the rh multiplication - * @param prod OUTPUT vector to receive the result - */ virtual void mult(const doublereal* b, doublereal* prod) const; - - //! Multiply b*A and write result to prod. - /*! - * @param b Vector to do the lh multiplication - * @param prod OUTPUT vector to receive the result - */ virtual void leftMult(const doublereal* const b, doublereal* const prod) const; //! Perform an LU decomposition, the LAPACK routine DGBTRF is used. /*! - * * The factorization is saved in ludata. * * @return Return a success flag. @@ -207,7 +171,6 @@ public: */ int factor(); - //! Solve the matrix problem Ax = b /*! * @param b INPUT rhs of the problem @@ -230,7 +193,6 @@ public: */ int solve(doublereal* b); - //! Returns an iterator for the start of the band storage data /*! * Iterator points to the beginning of the data, and it is changeable. @@ -255,9 +217,6 @@ public: */ vector_fp::const_iterator end() const; - /** - * Zero the matrix - */ virtual void zero(); //! Factors the A matrix using the QR algorithm, overwriting A @@ -306,7 +265,6 @@ public: //! Returns the one norm of the matrix virtual doublereal oneNorm() const; - //! Duplicate this object as a GeneralMatrix pointer virtual GeneralMatrix* duplMyselfAsGeneralMatrix() const; //! Report whether the current matrix has been factored. @@ -316,11 +274,11 @@ public: /*! * The LAPACK bandstructure has column values which are contiguous in memory: * - * On entry, the matrix A in band storage, in rows KL+1 to - * 2*KL+KU+1; rows 1 to KL of the array need not be set. - * The j-th column of A is stored in the j-th column of the - * array AB as follows: - * AB(KL + KU + 1 + i - j,j) = A(i,j) for max(1, j - KU) <= i <= min(m, j + KL) + * On entry, the matrix A in band storage, in rows KL+1 to + * 2*KL+KU+1; rows 1 to KL of the array need not be set. + * The j-th column of A is stored in the j-th column of the + * array AB as follows: + * AB(KL + KU + 1 + i - j,j) = A(i,j) for max(1, j - KU) <= i <= min(m, j + KL) * * This routine returns the position of AB(1,j) (fortran-1 indexing) in the above format * @@ -329,7 +287,6 @@ public: * double *colP_j = matrix.ptrColumn(j); * double a_i_j = colP_j[kl + ku + i - j]; * - * * @param j Value of the column * * @return Returns a pointer to the top of the column @@ -353,7 +310,6 @@ public: */ virtual void copyData(const GeneralMatrix& y); - //! Clear the factored flag virtual void clearFactorFlag(); @@ -421,8 +377,6 @@ private: * @param msg String containing the message. */ void err(const std::string& msg) const; - - }; //! Utility routine to print out the matrix diff --git a/include/cantera/numerics/DenseMatrix.h b/include/cantera/numerics/DenseMatrix.h index dc8204617..786cbae98 100644 --- a/include/cantera/numerics/DenseMatrix.h +++ b/include/cantera/numerics/DenseMatrix.h @@ -69,9 +69,7 @@ public: */ class DenseMatrix : public Array2D { - public: - //! Default Constructor DenseMatrix(); @@ -107,14 +105,6 @@ public: */ void resize(size_t n, size_t m, doublereal v = 0.0); - //! Return a vector of const pointers to the columns - /*! - * Note the value of the pointers are protected by their being const. - * However, the value of the matrix is open to being changed. - * - * @return returns a vector of pointers to the top of the columns - * of the matrices. - */ virtual doublereal* const* colPts(); //! Return a const vector of const pointers to the columns @@ -127,17 +117,10 @@ public: */ const doublereal* const* const_colPts() const; - //! 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; //! Multiply A*B and write result to \c prod. /*! - * * @param b input DenseMatrix B of size NxN * @param prod output output DenseMatrix prod size NxN */ @@ -147,7 +130,6 @@ public: /*! * @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 */ @@ -168,7 +150,6 @@ public: } protected: - //! Vector of pivots. Length is equal to the max of m and n. vector_int m_ipiv; @@ -176,7 +157,6 @@ protected: std::vector m_colPts; public: - //! Error Handling Flag /*! * The default is to set this to 0. In this case, if a factorization is requested and can't be achieved, @@ -194,7 +174,6 @@ public: */ int m_printLevel; - // Listing of friend functions which are defined below friend int solve(DenseMatrix& A, double* b); @@ -202,8 +181,6 @@ public: friend int invert(DenseMatrix& A, int nn); }; -//================================================================================================================== - //! Solve Ax = b. Array b is overwritten on exit with x. /*! @@ -265,6 +242,3 @@ int invert(DenseMatrix& A, size_t nn=npos); } #endif - - - diff --git a/include/cantera/numerics/GeneralMatrix.h b/include/cantera/numerics/GeneralMatrix.h index 20524a479..6e9374580 100644 --- a/include/cantera/numerics/GeneralMatrix.h +++ b/include/cantera/numerics/GeneralMatrix.h @@ -32,15 +32,9 @@ public: GeneralMatrix(int matType); //! Copy Constructor - /*! - * @param right Object to be copied - */ GeneralMatrix(const GeneralMatrix& right); //! Assignment operator - /*! - * @param right Object to be copied - */ GeneralMatrix& operator=(const GeneralMatrix& right); //! Destructor. Does nothing. @@ -48,7 +42,7 @@ public: //! Duplicator member function /*! - * This function will duplicate the matrix given a generic GeneralMatrix pointer + * This function will duplicate the matrix given a generic GeneralMatrix * * @return Returns a pointer to the malloced object */ @@ -112,26 +106,16 @@ public: virtual void useFactorAlgorithm(int fAlgorithm) = 0; //! Return the factor algorithm used - /*! - * - */ virtual int factorAlgorithm() const = 0; //! Calculate the one norm of the matrix - /*! - * Returns the one norm of the matrix - */ virtual doublereal oneNorm() const = 0; - //! Return the number of rows in the matrix virtual size_t nRows() const = 0; - //! Return the size and structure of the matrix /*! - * This is inherited from GeneralMatrix - * * @param iStruct OUTPUT Pointer to a vector of ints that describe the structure of the matrix. * * @return returns the number of rows and columns in the matrix. @@ -167,7 +151,6 @@ public: */ virtual doublereal& operator()(size_t i, size_t j) = 0; - //! Constant Index into the (i,j) element /*! * @param i row diff --git a/include/cantera/numerics/SquareMatrix.h b/include/cantera/numerics/SquareMatrix.h index c1f66e97e..f71a37b1b 100644 --- a/include/cantera/numerics/SquareMatrix.h +++ b/include/cantera/numerics/SquareMatrix.h @@ -27,9 +27,6 @@ class SquareMatrix: public DenseMatrix, public GeneralMatrix { public: //! Base Constructor. - /*! - * Create an \c 0 by \c 0 matrix, and initialize all elements to \c 0. - */ SquareMatrix(); //! Constructor. @@ -42,91 +39,28 @@ public: SquareMatrix(size_t n, doublereal v = 0.0); //! Copy Constructor - /*! - * @param right Object to be copied - */ SquareMatrix(const SquareMatrix& right); //! Assignment operator - /*! - * @param right Object to be copied - */ SquareMatrix& operator=(const SquareMatrix& right); - //! Solves the Ax = b system returning x in the b spot. - /*! - * @param b Vector for the rhs of the equation system - */ int solve(doublereal* b); - //! Resize the matrix - /*! - * @param n Number of rows - * @param m Number of columns - * @param v double to fill the new space (defaults to zero) - */ void resize(size_t n, size_t m, doublereal v = 0.0); - /** - * Zero the matrix - */ + //! Zero the matrix void zero(); - //! Multiply A*b and write result to prod. - /*! - * @param b Vector to do the rh multiplication - * @param prod OUTPUT vector to receive the result - */ virtual void mult(const doublereal* b, doublereal* prod) const; - - //! Multiply A*B and write result to \c prod. - /*! - * - * @param b input DenseMatrix B of size NxN - * @param prod output output DenseMatrix prod size NxN - */ virtual void mult(const DenseMatrix& b, DenseMatrix& prod) const; - - //! Multiply b*A and write result to prod. - /*! - * @param b Vector to do the lh multiplication - * @param prod OUTPUT vector to receive the result - */ virtual void leftMult(const doublereal* const b, doublereal* const prod) const; - /** - * Factors the A matrix, overwriting A. We flip m_factored - * boolean to indicate that the matrix is now A-1. - */ int factor(); - - //! Factors the A matrix using the QR algorithm, overwriting A - /*! - * we set m_factored to 2 to indicate the matrix is now QR factored - * - * @return Returns the info variable from lapack - */ virtual int factorQR(); - //! Returns an estimate of the inverse of the condition number for the matrix - /*! - * The matrix must have been previously factored using the QR algorithm - * - * @return returns the inverse of the condition number - */ virtual doublereal rcondQR(); - - //! Returns an estimate of the inverse of the condition number for the matrix - /*! - * The matrix must have been previously factored using the LU algorithm - * - * @param a1norm Norm of the matrix - * - * @return returns the inverse of the condition number - */ virtual doublereal rcond(doublereal a1norm); - //! Returns the one norm of the matrix virtual doublereal oneNorm() const; //! Solves the linear problem Ax=b using the QR algorithm returning x in the b spot @@ -135,22 +69,12 @@ public: */ int solveQR(doublereal* b); - - //! clear the factored flag virtual void clearFactorFlag(); //! set the factored flag void setFactorFlag(); - //! Report whether the current matrix has been factored. virtual bool factored() const; - - //! Change the way the matrix is factored - /*! - * @param fAlgorithm integer - * 0 LU factorization - * 1 QR factorization - */ virtual void useFactorAlgorithm(int fAlgorithm); //! Returns the factor algorithm used @@ -162,46 +86,18 @@ public: */ virtual int factorAlgorithm() const; - //! Return a pointer to the top of column j, columns are assumed to be contiguous in memory - /*! - * @param j Value of the column - * - * @return Returns a pointer to the top of the column - */ virtual doublereal* ptrColumn(size_t j); - //! Index into the (i,j) element - /*! - * @param i row - * @param j column - * - * (note, tried a using directive here, and it didn't seem to work) - * - * Returns a changeable reference to the matrix entry - */ virtual doublereal& operator()(size_t i, size_t j) { return Array2D::operator()(i, j); } - //! Copy the data from one array into another without doing any checking - /*! - * This differs from the assignment operator as no resizing is done and memcpy() is used. - * @param y Array to be copied - */ virtual void copyData(const GeneralMatrix& y); - //! Constant Index into the (i,j) element - /*! - * @param i row - * @param j column - * - * Returns an unchangeable reference to the matrix entry - */ virtual doublereal operator()(size_t i, size_t j) const { return Array2D::operator()(i, j); } - //! Return the number of rows in the matrix virtual size_t nRows() const; //! Return the size and structure of the matrix @@ -215,54 +111,17 @@ public: */ size_t nRowsAndStruct(size_t* const iStruct = 0) const; - //! Duplicate this object virtual GeneralMatrix* duplMyselfAsGeneralMatrix() const; - - //! Return an iterator pointing to the first element - /*! - */ virtual vector_fp::iterator begin(); - - - //! Return a const iterator pointing to the first element virtual vector_fp::const_iterator begin() const; - - //! Return a vector of const pointers to the columns - /*! - * Note the value of the pointers are protected by their being const. - * However, the value of the matrix is open to being changed. - * - * @return returns a vector of pointers to the top of the columns - * of the matrices. - */ virtual doublereal* const* colPts(); - //! Check to see if we have any zero rows in the jacobian - /*! - * This utility routine checks to see if any rows are zero. - * The smallest row is returned along with the largest coefficient in that row - * - * @param valueSmall OUTPUT value of the largest coefficient in the smallest row - * - * @return index of the row that is most nearly zero - */ virtual size_t checkRows(doublereal& valueSmall) const; - - //! Check to see if we have any zero columns in the jacobian - /*! - * This utility routine checks to see if any columns are zero. - * The smallest column is returned along with the largest coefficient in that column - * - * @param valueSmall OUTPUT value of the largest coefficient in the smallest column - * - * @return index of the column that is most nearly zero - */ virtual size_t checkColumns(doublereal& valueSmall) const; protected: - //! the factor flag int m_factored; @@ -285,4 +144,3 @@ protected: } #endif - diff --git a/src/numerics/BandMatrix.cpp b/src/numerics/BandMatrix.cpp index 664a060ac..a449a9049 100644 --- a/src/numerics/BandMatrix.cpp +++ b/src/numerics/BandMatrix.cpp @@ -21,7 +21,6 @@ using namespace std; namespace Cantera { -//==================================================================================================================== BandMatrix::BandMatrix() : GeneralMatrix(1), m_factored(false), @@ -33,7 +32,7 @@ BandMatrix::BandMatrix() : data.clear(); ludata.clear(); } -//==================================================================================================================== + BandMatrix::BandMatrix(size_t n, size_t kl, size_t ku, doublereal v) : GeneralMatrix(1), m_factored(false), @@ -53,7 +52,7 @@ BandMatrix::BandMatrix(size_t n, size_t kl, size_t ku, doublereal v) : m_colPtrs[j] = &(data[ldab * j]); } } -//==================================================================================================================== + BandMatrix::BandMatrix(const BandMatrix& y) : GeneralMatrix(1), m_factored(false), @@ -75,7 +74,7 @@ BandMatrix::BandMatrix(const BandMatrix& y) : m_colPtrs[j] = &(data[ldab * j]); } } -//==================================================================================================================== + BandMatrix& BandMatrix::operator=(const BandMatrix& y) { if (&y == this) { @@ -96,7 +95,7 @@ BandMatrix& BandMatrix::operator=(const BandMatrix& y) } return *this; } -//==================================================================================================================== + void BandMatrix::resize(size_t n, size_t kl, size_t ku, doublereal v) { m_n = n; @@ -113,29 +112,29 @@ void BandMatrix::resize(size_t n, size_t kl, size_t ku, doublereal v) } m_factored = false; } -//==================================================================================================================== + void BandMatrix::bfill(doublereal v) { std::fill(data.begin(), data.end(), v); m_factored = false; } -//==================================================================================================================== + void BandMatrix::zero() { std::fill(data.begin(), data.end(), 0.0); m_factored = false; } -//==================================================================================================================== + doublereal& BandMatrix::operator()(size_t i, size_t j) { return value(i,j); } -//==================================================================================================================== + doublereal BandMatrix::operator()(size_t i, size_t j) const { return value(i,j); } -//==================================================================================================================== + doublereal& BandMatrix::value(size_t i, size_t j) { m_factored = false; @@ -144,7 +143,7 @@ doublereal& BandMatrix::value(size_t i, size_t j) } return data[index(i,j)]; } -//==================================================================================================================== + doublereal BandMatrix::value(size_t i, size_t j) const { if (i + m_ku < j || i > j + m_kl) { @@ -152,7 +151,7 @@ doublereal BandMatrix::value(size_t i, size_t j) const } return data[index(i,j)]; } -//==================================================================================================================== + size_t BandMatrix::index(size_t i, size_t j) const { int jj = static_cast(j); @@ -160,19 +159,17 @@ size_t BandMatrix::index(size_t i, size_t j) const size_t rw = (int) m_kl + (int) m_ku + (int) ii - jj; return (2*m_kl + m_ku + 1)*j + rw; } -//==================================================================================================================== + doublereal BandMatrix::_value(size_t i, size_t j) const { return data[index(i,j)]; } -//==================================================================================================================== -// Number of rows + size_t BandMatrix::nRows() const { return m_n; } -//==================================================================================================================== -// Number of rows + size_t BandMatrix::nRowsAndStruct(size_t* const iStruct) const { if (iStruct) { @@ -229,10 +226,7 @@ void BandMatrix::mult(const doublereal* b, doublereal* prod) const prod[m] = sum; } } -//==================================================================================================================== -/* - * Multiply b*A and write result to \c prod. - */ + void BandMatrix::leftMult(const doublereal* const b, doublereal* const prod) const { int kl = static_cast(m_kl); @@ -250,11 +244,7 @@ void BandMatrix::leftMult(const doublereal* const b, doublereal* const prod) con prod[n] = sum; } } -//==================================================================================================================== -/* - * Perform an LU decomposition. LAPACK routine DGBTRF is used. - * The factorization is saved in ludata. - */ + int BandMatrix::factor() { int info=0; @@ -273,13 +263,13 @@ int BandMatrix::factor() } return info; } -//==================================================================================================================== + int BandMatrix::solve(const doublereal* const b, doublereal* const x) { copy(b, b + m_n, x); return solve(x); } -//==================================================================================================================== + int BandMatrix::solve(doublereal* b) { int info = 0; @@ -299,29 +289,29 @@ int BandMatrix::solve(doublereal* b) } return info; } -//==================================================================================================================== + vector_fp::iterator BandMatrix::begin() { m_factored = false; return data.begin(); } -//==================================================================================================================== + vector_fp::iterator BandMatrix::end() { m_factored = false; return data.end(); } -//==================================================================================================================== + vector_fp::const_iterator BandMatrix::begin() const { return data.begin(); } -//==================================================================================================================== + vector_fp::const_iterator BandMatrix::end() const { return data.end(); } -//==================================================================================================================== + ostream& operator<<(ostream& s, const BandMatrix& m) { size_t nr = m.nRows(); @@ -334,45 +324,24 @@ ostream& operator<<(ostream& s, const BandMatrix& m) } return s; } -//==================================================================================================================== + void BandMatrix::err(const std::string& msg) const { throw CanteraError("BandMatrix() unimplemented function", msg); } -//==================================================================================================================== -// Factors the A matrix using the QR algorithm, overwriting A -/* - * we set m_factored to 2 to indicate the matrix is now QR factored - * - * @return Returns the info variable from lapack - */ + int BandMatrix::factorQR() { factor(); return 0; } -//==================================================================================================================== -// Factors the A matrix using the QR algorithm, overwriting A -// Returns an estimate of the inverse of the condition number for the matrix -/* - * The matrix must have been previously factored using the QR algorithm - * - * @return returns the inverse of the condition number - */ + doublereal BandMatrix::rcondQR() { double a1norm = oneNorm(); return rcond(a1norm); } -//==================================================================================================================== -// Returns an estimate of the inverse of the condition number for the matrix -/* - * The matrix must have been previously factored using the LU algorithm - * - * @param a1norm Norm of the matrix - * - * @return returns the inverse of the condition number - */ + doublereal BandMatrix::rcond(doublereal a1norm) { int printLevel = 0; @@ -403,24 +372,17 @@ doublereal BandMatrix::rcond(doublereal a1norm) } return rcond; } -//==================================================================================================================== -// Change the way the matrix is factored -/* - * @param fAlgorithm integer - * 0 LU factorization - * 1 QR factorization - */ + void BandMatrix::useFactorAlgorithm(int fAlgorithm) { // QR algorithm isn't implemented for banded matrix. } -//==================================================================================================================== + int BandMatrix::factorAlgorithm() const { return 0; } -//==================================================================================================================== -// Returns the one norm of the matrix + doublereal BandMatrix::oneNorm() const { int ku = static_cast(m_ku); @@ -438,7 +400,7 @@ doublereal BandMatrix::oneNorm() const } return value; } -//==================================================================================================================== + size_t BandMatrix::checkRows(doublereal& valueSmall) const { valueSmall = 1.0E300; @@ -464,7 +426,7 @@ size_t BandMatrix::checkRows(doublereal& valueSmall) const } return iSmall; } -//==================================================================================================================== + size_t BandMatrix::checkColumns(doublereal& valueSmall) const { valueSmall = 1.0E300; @@ -490,46 +452,27 @@ size_t BandMatrix::checkColumns(doublereal& valueSmall) const } return jSmall; } -//==================================================================================================================== + GeneralMatrix* BandMatrix::duplMyselfAsGeneralMatrix() const { return new BandMatrix(*this); } -//==================================================================================================================== + bool BandMatrix::factored() const { return m_factored; } -//==================================================================================================================== -// Return a pointer to the top of column j, columns are assumed to be contiguous in memory -/* - * @param j Value of the column - * - * @return Returns a pointer to the top of the column - */ + doublereal* BandMatrix::ptrColumn(size_t j) { return m_colPtrs[j]; } -//==================================================================================================================== -// Return a vector of const pointers to the columns -/* - * Note the value of the pointers are protected by their being const. - * However, the value of the matrix is open to being changed. - * - * @return returns a vector of pointers to the top of the columns - * of the matrices. - */ + doublereal* const* BandMatrix::colPts() { return &(m_colPtrs[0]); } -//==================================================================================================================== -// Copy the data from one array into another without doing any checking -/* - * This differs from the assignment operator as no resizing is done and memcpy() is used. - * @param y Array to be copied - */ + void BandMatrix::copyData(const GeneralMatrix& y) { m_factored = false; @@ -537,15 +480,10 @@ void BandMatrix::copyData(const GeneralMatrix& y) GeneralMatrix* yyPtr = const_cast(&y); (void) memcpy(DATA_PTR(data), yyPtr->ptrColumn(0), n); } -//==================================================================================================================== -/* - * clear the factored flag - */ + void BandMatrix::clearFactorFlag() { m_factored = 0; } -//==================================================================================================================== -//==================================================================================================================== -} +} diff --git a/src/numerics/DenseMatrix.cpp b/src/numerics/DenseMatrix.cpp index d485df4d1..95462fe91 100644 --- a/src/numerics/DenseMatrix.cpp +++ b/src/numerics/DenseMatrix.cpp @@ -13,8 +13,7 @@ namespace Cantera { -//==================================================================================================================== -// Default Constructor + DenseMatrix::DenseMatrix() : Array2D(0,0,0.0), m_ipiv(0), @@ -22,11 +21,7 @@ DenseMatrix::DenseMatrix() : m_printLevel(0) { } -//==================================================================================================================== -/* - * Constructor. Create an \c n by \c m matrix, and initialize - * all elements to \c v. - */ + DenseMatrix::DenseMatrix(size_t n, size_t m, doublereal v) : Array2D(n, m, v), m_ipiv(0), @@ -41,11 +36,7 @@ DenseMatrix::DenseMatrix(size_t n, size_t m, doublereal v) : } } } -//==================================================================================================================== -// Copy constructor -/* - * @param y Object to be copied - */ + DenseMatrix::DenseMatrix(const DenseMatrix& y) : Array2D(y), m_ipiv(0), @@ -60,8 +51,7 @@ DenseMatrix::DenseMatrix(const DenseMatrix& y) : } } } -//==================================================================================================================== -// assignment + DenseMatrix& DenseMatrix::operator=(const DenseMatrix& y) { if (&y == this) { @@ -77,7 +67,7 @@ DenseMatrix& DenseMatrix::operator=(const DenseMatrix& y) m_printLevel = y.m_printLevel; return *this; } -//==================================================================================================================== + void DenseMatrix::resize(size_t n, size_t m, doublereal v) { Array2D::resize(n,m,v); @@ -89,17 +79,17 @@ void DenseMatrix::resize(size_t n, size_t m, doublereal v) } } } -//==================================================================================================================== + doublereal* const* DenseMatrix::colPts() { return &(m_colPts[0]); } -//==================================================================================================================== + const doublereal* const* DenseMatrix::const_colPts() const { return &(m_colPts[0]); } -//==================================================================================================================== + void DenseMatrix::mult(const double* b, double* prod) const { ct_dgemv(ctlapack::ColMajor, ctlapack::NoTranspose, @@ -107,7 +97,7 @@ void DenseMatrix::mult(const double* b, double* prod) const static_cast(nRows()), 1.0, ptrColumn(0), //begin(), static_cast(nRows()), b, 1, 0.0, prod, 1); } -//==================================================================================================================== + void DenseMatrix::mult(const DenseMatrix& B, DenseMatrix& prod) const { if (m_ncols != B.nColumns() || m_nrows != B.nRows() || m_ncols != m_nrows || m_ncols != prod.nColumns() || m_nrows != prod.nColumns()) { @@ -121,7 +111,7 @@ void DenseMatrix::mult(const DenseMatrix& B, DenseMatrix& prod) const mult(bcols[col], prodcols[col]); } } -//==================================================================================================================== + void DenseMatrix::leftMult(const double* const b, double* const prod) const { size_t nc = nColumns(); @@ -135,12 +125,12 @@ void DenseMatrix::leftMult(const double* const b, double* const prod) const prod[n] = sum; } } -//==================================================================================================================== + vector_int& DenseMatrix::ipiv() { return m_ipiv; } -//==================================================================================================================== + int solve(DenseMatrix& A, double* b) { int info = 0; @@ -186,7 +176,7 @@ int solve(DenseMatrix& A, double* b) } return info; } -//==================================================================================================================== + int solve(DenseMatrix& A, DenseMatrix& b) { int info = 0; @@ -239,21 +229,21 @@ int solve(DenseMatrix& A, DenseMatrix& b) return info; } -//==================================================================================================================== + 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) { 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, size_t nn) { integer n = static_cast(nn != npos ? nn : A.nRows()); @@ -284,6 +274,5 @@ int invert(DenseMatrix& A, size_t nn) } return info; } -//==================================================================================================================== -} +} diff --git a/src/numerics/GeneralMatrix.cpp b/src/numerics/GeneralMatrix.cpp index 0d353c50f..84473017b 100644 --- a/src/numerics/GeneralMatrix.cpp +++ b/src/numerics/GeneralMatrix.cpp @@ -14,17 +14,17 @@ using namespace std; namespace Cantera { -//==================================================================================================================== + GeneralMatrix::GeneralMatrix(int matType) : matrixType_(matType) { } -//==================================================================================================================== + GeneralMatrix::GeneralMatrix(const GeneralMatrix& y) : matrixType_(y.matrixType_) { } -//==================================================================================================================== + GeneralMatrix& GeneralMatrix::operator=(const GeneralMatrix& y) { if (&y == this) { @@ -33,9 +33,9 @@ GeneralMatrix& GeneralMatrix::operator=(const GeneralMatrix& y) matrixType_ = y.matrixType_; return *this; } -//==================================================================================================================== + GeneralMatrix::~GeneralMatrix() { } -//==================================================================================================================== + } diff --git a/src/numerics/SquareMatrix.cpp b/src/numerics/SquareMatrix.cpp index 7e2fc35a6..f521a2384 100644 --- a/src/numerics/SquareMatrix.cpp +++ b/src/numerics/SquareMatrix.cpp @@ -24,7 +24,6 @@ using namespace std; namespace Cantera { -//==================================================================================================================== SquareMatrix::SquareMatrix() : DenseMatrix(), GeneralMatrix(0), @@ -34,15 +33,6 @@ SquareMatrix::SquareMatrix() : { } -//==================================================================================================================== -// Constructor. -/* - * Create an \c n by \c n matrix, and initialize - * all elements to \c v. - * - * @param n size of the square matrix - * @param v initial value of all matrix components. - */ SquareMatrix::SquareMatrix(size_t n, doublereal v) : DenseMatrix(n, n, v), GeneralMatrix(0), @@ -52,11 +42,7 @@ SquareMatrix::SquareMatrix(size_t n, doublereal v) : { } -//==================================================================================================================== -/* - * - * copy constructor - */ + SquareMatrix::SquareMatrix(const SquareMatrix& y) : DenseMatrix(y), GeneralMatrix(0), @@ -66,10 +52,6 @@ SquareMatrix::SquareMatrix(const SquareMatrix& y) : { } -//==================================================================================================================== -/* - * Assignment operator - */ SquareMatrix& SquareMatrix::operator=(const SquareMatrix& y) { if (&y == this) { @@ -82,10 +64,7 @@ SquareMatrix& SquareMatrix::operator=(const SquareMatrix& y) useQR_ = y.useQR_; return *this; } -//==================================================================================================================== -/* - * Solve Ax = b. Vector b is overwritten on exit with x. - */ + int SquareMatrix::solve(doublereal* b) { if (useQR_) { @@ -117,10 +96,7 @@ int SquareMatrix::solve(doublereal* b) } return info; } -//==================================================================================================================== -/* - * Set all entries to zero - */ + void SquareMatrix::zero() { size_t n = nRows(); @@ -134,47 +110,27 @@ void SquareMatrix::zero() (void) memset((void*) sm, 0, nn * sizeof(double)); } } -//==================================================================================================================== + void SquareMatrix::resize(size_t n, size_t m, doublereal v) { DenseMatrix::resize(n, m, v); } -//==================================================================================================================== -// Multiply A*b and write result to prod. -/* - * @param b Vector to do the rh multiplication - * @param prod OUTPUT vector to receive the result - */ void SquareMatrix::mult(const doublereal* b, doublereal* prod) const { DenseMatrix::mult(b, prod); } -//==================================================================================================================== -// Multiply A*B and write result to \c prod. -/* - * - * @param b input DenseMatrix B of size NxN - * @param prod output output DenseMatrix prod size NxN - */ + void SquareMatrix::mult(const DenseMatrix& b, DenseMatrix& prod) const { DenseMatrix::mult(b, prod); } -//==================================================================================================================== -// Multiply b*A and write result to prod. -/* - * @param b Vector to do the lh multiplication - * @param prod OUTPUT vector to receive the result - */ + void SquareMatrix::leftMult(const doublereal* const b, doublereal* const prod) const { DenseMatrix::leftMult(b, prod); } -//==================================================================================================================== -/* - * Factor A. A is overwritten with the LU decomposition of A. - */ + int SquareMatrix::factor() { if (useQR_) { @@ -195,23 +151,17 @@ int SquareMatrix::factor() } return info; } -//===================================================================================================================== -/* - * clear the factored flag - */ + void SquareMatrix::clearFactorFlag() { m_factored = 0; } -//===================================================================================================================== -/* - * set the factored flag - */ + void SquareMatrix::setFactorFlag() { m_factored = 1; } -//===================================================================================================================== + int SquareMatrix::factorQR() { if (tau.size() < m_nrows) { @@ -239,10 +189,7 @@ int SquareMatrix::factorQR() return info; } -//===================================================================================================================== -/* - * Solve Ax = b. Vector b is overwritten on exit with x. - */ + int SquareMatrix::solveQR(doublereal* b) { int info=0; @@ -295,7 +242,7 @@ int SquareMatrix::solveQR(doublereal* b) return info; } -//===================================================================================================================== + doublereal SquareMatrix::rcond(doublereal anorm) { @@ -326,12 +273,12 @@ doublereal SquareMatrix::rcond(doublereal anorm) } return rcond; } -//===================================================================================================================== + doublereal SquareMatrix::oneNorm() const { return a1norm_; } -//===================================================================================================================== + doublereal SquareMatrix::rcondQR() { @@ -359,84 +306,62 @@ doublereal SquareMatrix::rcondQR() } return rcond; } -//===================================================================================================================== + void SquareMatrix::useFactorAlgorithm(int fAlgorithm) { useQR_ = fAlgorithm; } -//===================================================================================================================== + int SquareMatrix::factorAlgorithm() const { return (int) useQR_; } -//===================================================================================================================== + bool SquareMatrix::factored() const { return (m_factored != 0); } -//===================================================================================================================== -// Return a pointer to the top of column j, columns are contiguous in memory -/* - * @param j Value of the column - * - * @return Returns a pointer to the top of the column - */ + doublereal* SquareMatrix::ptrColumn(size_t j) { return Array2D::ptrColumn(j); } -//===================================================================================================================== -// Copy the data from one array into another without doing any checking -/* - * This differs from the assignment operator as no resizing is done and memcpy() is used. - * @param y Array to be copied - */ + void SquareMatrix::copyData(const GeneralMatrix& y) { const SquareMatrix* yy_ptr = dynamic_cast(& y); Array2D::copyData(*yy_ptr); } -//===================================================================================================================== + size_t SquareMatrix::nRows() const { return m_nrows; } -//===================================================================================================================== + size_t SquareMatrix::nRowsAndStruct(size_t* const iStruct) const { return m_nrows; } -//===================================================================================================================== + GeneralMatrix* SquareMatrix::duplMyselfAsGeneralMatrix() const { return new SquareMatrix(*this); } -//===================================================================================================================== -// Return an iterator pointing to the first element + vector_fp::iterator SquareMatrix::begin() { return m_data.begin(); } -//===================================================================================================================== -// Return a const iterator pointing to the first element + vector_fp::const_iterator SquareMatrix::begin() const { return m_data.begin(); } -//===================================================================================================================== -// Return a vector of const pointers to the columns -/* - * Note the value of the pointers are protected by their being const. - * However, the value of the matrix is open to being changed. - * - * @return returns a vector of pointers to the top of the columns - * of the matrices. - */ + doublereal* const* SquareMatrix::colPts() { return DenseMatrix::colPts(); } -//===================================================================================================================== size_t SquareMatrix::checkRows(doublereal& valueSmall) const { @@ -456,7 +381,7 @@ size_t SquareMatrix::checkRows(doublereal& valueSmall) const } return iSmall; } -//===================================================================================================================== + size_t SquareMatrix::checkColumns(doublereal& valueSmall) const { valueSmall = 1.0E300; @@ -475,8 +400,5 @@ size_t SquareMatrix::checkColumns(doublereal& valueSmall) const } return jSmall; } -//===================================================================================================================== - } -