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

Cleaned up Doxygen Documentation for matrix classes

Browse source 
BandMatrix, DenseMatrix, SquareMatrix, and GeneralMatrix
This commit is contained in:
Ray Speth 2013-05-29 00:11:22 +00:00
parent 7e8276809f
commit e9daef3313
8 changed files with 94 additions and 476 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

56
include/cantera/numerics/BandMatrix.h
View file

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

26
include/cantera/numerics/DenseMatrix.h
View file

@ -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<doublereal*> 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

19
include/cantera/numerics/GeneralMatrix.h
View file

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

144
include/cantera/numerics/SquareMatrix.h
View file

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

140
src/numerics/BandMatrix.cpp
View file

@ -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<int>(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<int>(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<int>(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<GeneralMatrix*>(&y);
(void) memcpy(DATA_PTR(data), yyPtr->ptrColumn(0), n);
}
//====================================================================================================================
/*
* clear the factored flag
*/
void BandMatrix::clearFactorFlag()
{
m_factored = 0;
}
//====================================================================================================================
//====================================================================================================================
}
}

45
src/numerics/DenseMatrix.cpp
View file

@ -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<int>(nRows()), 1.0, ptrColumn(0), //begin(),
static_cast<int>(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<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)
{
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);
}
//====================================================================================================================
int invert(DenseMatrix& A, size_t nn)
{
integer n = static_cast<int>(nn != npos ? nn : A.nRows());
@ -284,6 +274,5 @@ int invert(DenseMatrix& A, size_t nn)
}
return info;
}
//====================================================================================================================
}
}

10
src/numerics/GeneralMatrix.cpp
View file

@ -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()
{
}
//====================================================================================================================
}

130
src/numerics/SquareMatrix.cpp
View file

@ -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<const SquareMatrix*>(& 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;
}
//=====================================================================================================================
}