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

Verified that the QR factorization works on one problem.

Browse source
This commit is contained in:
Harry Moffat 2011-04-07 16:14:00 +00:00
parent 62af5826ea
commit 155d7604b5
2 changed files with 58 additions and 28 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

37
Cantera/src/numerics/SquareMatrix.cpp
View file

@ -32,10 +32,12 @@ namespace Cantera {
//====================================================================================================================
SquareMatrix::SquareMatrix() :
DenseMatrix(),
m_factored(0)
m_factored(0),
useQR_(0)
{
}
//====================================================================================================================
// Constructor.
/*
* Create an \c n by \c n matrix, and initialize
@ -46,21 +48,24 @@ namespace Cantera {
*/
SquareMatrix::SquareMatrix(int n, doublereal v) :
DenseMatrix(n, n, v),
m_factored(0)
m_factored(0),
useQR_(0)
{
}
//====================================================================================================================
/*
*
* copy constructor
*/
SquareMatrix::SquareMatrix(const SquareMatrix& y) :
DenseMatrix(y),
m_factored(y.m_factored)
m_factored(y.m_factored),
useQR_(y.useQR_)
{
}
/**
//====================================================================================================================
/*
* Assignment operator
*/
SquareMatrix& SquareMatrix::operator=(const SquareMatrix& y) {
@ -69,12 +74,18 @@ namespace Cantera {
m_factored = y.m_factored;
return *this;
}
/**
//====================================================================================================================
SquareMatrix::~SquareMatrix() {
}
//====================================================================================================================
/*
* Solve Ax = b. Vector b is overwritten on exit with x.
*/
int SquareMatrix::solve(double* b)
{
if (useQR_) {
return solveQR(b);
}
int info=0;
/*
* Check to see whether the matrix has been factored.
@ -122,10 +133,13 @@ namespace Cantera {
DenseMatrix::resize(n, m, v);
}
//====================================================================================================================
/**
/*
* Factor A. A is overwritten with the LU decomposition of A.
*/
int SquareMatrix::factor() {
if (useQR_) {
return factorQR();
}
integer n = static_cast<int>(nRows());
int info=0;
m_factored = 1;
@ -159,9 +173,10 @@ namespace Cantera {
int SquareMatrix::factorQR() {
if ((int) tau.size() < m_nrows) {
tau.resize(m_nrows, 0.0);
work.resize(9 * m_nrows, 0.0);
work.resize(8 * m_nrows, 0.0);
}
int info;
m_factored = 2;
int lwork = work.size();
ct_dgeqrf(m_nrows, m_nrows, &(*(begin())), m_nrows, DATA_PTR(tau), DATA_PTR(work), lwork, info);
if (info != 0) {
@ -173,7 +188,7 @@ namespace Cantera {
}
}
int lworkOpt = work[0];
if (lworkOpt != lwork) {
if (lworkOpt > lwork) {
work.resize(lworkOpt);
}
return info;
@ -215,7 +230,7 @@ namespace Cantera {
}
}
int lworkOpt = work[0];
if (lworkOpt <= lwork) {
if (lworkOpt > lwork) {
work.resize(lworkOpt);
}

49
Cantera/src/numerics/SquareMatrix.h
View file

@ -45,25 +45,35 @@ namespace Cantera {
*/
SquareMatrix(int n, doublereal v = 0.0);
/**
* Copy Constructor
//! Copy Constructor
/*!
* @param right Object to be copied
*/
SquareMatrix(const SquareMatrix&);
/**
* Assignment operator
SquareMatrix(const SquareMatrix& right);
//! Assignment operator
/*!
* @param right Object to be copied
*/
SquareMatrix& operator=(const SquareMatrix&);
SquareMatrix& operator=(const SquareMatrix& right);
/// Destructor. Does nothing.
virtual ~SquareMatrix(){}
//! Destructor. Does nothing.
virtual ~SquareMatrix();
/**
* Solves the Ax = b system returning x in the b spot.
//! Solves the Ax = b system returning x in the b spot.
/*!
* @param b Vector for the rhs of the equation system
*/
int solve(double *b);
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(int n, int m, doublereal v = 0.0);
@ -86,16 +96,14 @@ namespace Cantera {
*/
int factorQR();
//! Solves the linear problem Ax=b using the QR algorithm
//! returning x in the b spot
//! Solves the linear problem Ax=b using the QR algorithm returning x in the b spot
/*!
* @param b RHS to be solved.
*/
int solveQR(doublereal *b);
/**
* clear the factored flag
*/
//! clear the factored flag
void clearFactorFlag();
/**
* set the factored flag
@ -107,8 +115,15 @@ namespace Cantera {
*/
int m_factored;
//! Work vector for QR algorithm
vector_fp tau;
//! Work vector for QR algorithm
vector_fp work;
public:
//! Use the QR algorithm to factor and invert the matrix
int useQR_;
};
}