Doxygen updates
-> Also moved implementation functions for BandMatrix into the cpp file
This commit is contained in:
Harry Moffat
parent
6cf9c46630
commit
e8d6649c7b
4 changed files with 497 additions and 282 deletions
|
|
@ -22,160 +22,243 @@ using namespace std;
|
|||
|
||||
namespace Cantera {
|
||||
|
||||
//====================================================================================================================
|
||||
BandMatrix::BandMatrix() :
|
||||
m_factored(false),
|
||||
m_n(0),
|
||||
m_kl(0),
|
||||
m_ku(0),
|
||||
m_zero(0.0)
|
||||
{
|
||||
data.clear();
|
||||
ludata.clear();
|
||||
}
|
||||
//====================================================================================================================
|
||||
BandMatrix::BandMatrix(int n, int kl, int ku, doublereal v) :
|
||||
m_factored(false),
|
||||
m_n(n),
|
||||
m_kl(kl),
|
||||
m_ku(ku),
|
||||
m_zero(0.0)
|
||||
{
|
||||
data.resize(n*(2*kl + ku + 1));
|
||||
ludata.resize(n*(2*kl + ku + 1));
|
||||
fill(data.begin(), data.end(), v);
|
||||
fill(ludata.begin(), ludata.end(), 0.0);
|
||||
m_ipiv.resize(m_n);
|
||||
}
|
||||
//====================================================================================================================
|
||||
BandMatrix::BandMatrix(const BandMatrix& y) :
|
||||
m_factored(false),
|
||||
m_n(0),
|
||||
m_kl(0),
|
||||
m_ku(0),
|
||||
m_zero(0.0)
|
||||
{
|
||||
m_n = y.m_n;
|
||||
m_kl = y.m_kl;
|
||||
m_ku = y.m_ku;
|
||||
data = y.data;
|
||||
ludata = y.ludata;
|
||||
m_factored = y.m_factored;
|
||||
m_ipiv = y.m_ipiv;
|
||||
}
|
||||
//====================================================================================================================
|
||||
BandMatrix::~BandMatrix() {
|
||||
|
||||
/// Default constructor.
|
||||
BandMatrix::BandMatrix() : m_factored(false), m_n(0),
|
||||
m_kl(0), m_ku(0), m_zero(0.0) {
|
||||
data.clear(); ludata.clear();
|
||||
}
|
||||
//====================================================================================================================
|
||||
BandMatrix& BandMatrix::operator=(const BandMatrix & y) {
|
||||
if (&y == this) return *this;
|
||||
m_n = y.m_n;
|
||||
m_kl = y.m_kl;
|
||||
m_ku = y.m_ku;
|
||||
m_ipiv = y.m_ipiv;
|
||||
data = y.data;
|
||||
ludata = y.ludata;
|
||||
m_factored = y.m_factored;
|
||||
return *this;
|
||||
}
|
||||
//====================================================================================================================
|
||||
void BandMatrix::resize(int n, int kl, int ku, doublereal v) {
|
||||
m_n = n;
|
||||
m_kl = kl;
|
||||
m_ku = ku;
|
||||
data.resize(n*(2*kl + ku + 1));
|
||||
ludata.resize(n*(2*kl + ku + 1));
|
||||
m_ipiv.resize(m_n);
|
||||
fill(data.begin(), data.end(), v);
|
||||
m_factored = false;
|
||||
}
|
||||
//====================================================================================================================
|
||||
void BandMatrix::bfill(doublereal v) {
|
||||
std::fill(data.begin(), data.end(), v);
|
||||
m_factored = false;
|
||||
}
|
||||
//====================================================================================================================
|
||||
doublereal& BandMatrix::operator()(int i, int j) {
|
||||
return value(i,j);
|
||||
}
|
||||
//====================================================================================================================
|
||||
doublereal BandMatrix::operator()(int i, int j) const {
|
||||
return value(i,j);
|
||||
}
|
||||
//====================================================================================================================
|
||||
doublereal& BandMatrix::value(int i, int j) {
|
||||
m_factored = false;
|
||||
if (i < j - m_ku || i > j + m_kl) {
|
||||
return m_zero;
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Constructor. Create an n by n banded matrix.
|
||||
* @param n number of rows and columns
|
||||
* @param kl number of subdiagonals
|
||||
* @param ku number of superdiagonals
|
||||
* @param v initial value (default = 0.0)
|
||||
*/
|
||||
BandMatrix::BandMatrix(int n, int kl, int ku, doublereal v)
|
||||
: m_factored(false), m_n(n), m_kl(kl), m_ku(ku) {
|
||||
data.resize(n*(2*kl + ku + 1));
|
||||
ludata.resize(n*(2*kl + ku + 1));
|
||||
fill(data.begin(), data.end(), v);
|
||||
fill(ludata.begin(), ludata.end(), 0.0);
|
||||
m_ipiv.resize(m_n);
|
||||
return data[index(i,j)];
|
||||
}
|
||||
//====================================================================================================================
|
||||
doublereal BandMatrix::value( int i, int j) const {
|
||||
if (i < j - m_ku || i > j + m_kl) return 0.0;
|
||||
return data[index(i,j)];
|
||||
}
|
||||
//====================================================================================================================
|
||||
int BandMatrix::index(int i, int j) const {
|
||||
int rw = m_kl + m_ku + i - j;
|
||||
return (2*m_kl + m_ku + 1)*j + rw;
|
||||
}
|
||||
//====================================================================================================================
|
||||
doublereal BandMatrix::_value(int i, int j) const {
|
||||
return data[index(i,j)];
|
||||
}
|
||||
//====================================================================================================================
|
||||
// Number of rows
|
||||
int BandMatrix::nRows() const {
|
||||
return m_n;
|
||||
}
|
||||
//====================================================================================================================
|
||||
// Number of columns
|
||||
int BandMatrix::nColumns() const {
|
||||
return m_n;
|
||||
}
|
||||
//====================================================================================================================
|
||||
// Number of subdiagonals
|
||||
int BandMatrix::nSubDiagonals() const {
|
||||
return m_kl;
|
||||
}
|
||||
//====================================================================================================================
|
||||
// Number of superdiagonals
|
||||
int BandMatrix::nSuperDiagonals() const {
|
||||
return m_ku;
|
||||
}
|
||||
//====================================================================================================================
|
||||
int BandMatrix::ldim() const {
|
||||
return 2*m_kl + m_ku + 1;
|
||||
}
|
||||
//====================================================================================================================
|
||||
vector_int & BandMatrix::ipiv() {
|
||||
return m_ipiv;
|
||||
}
|
||||
//====================================================================================================================
|
||||
/*
|
||||
* Multiply A*b and write result to \c prod.
|
||||
*/
|
||||
void BandMatrix::mult(const doublereal * const b, doublereal * const prod) const {
|
||||
int nr = nRows();
|
||||
doublereal sum = 0.0;
|
||||
for (int m = 0; m < nr; m++) {
|
||||
sum = 0.0;
|
||||
for (int j = m - m_kl; j <= m + m_ku; j++) {
|
||||
if (j >= 0 && j < m_n)
|
||||
sum += _value(m,j) * b[j];
|
||||
}
|
||||
prod[m] = sum;
|
||||
}
|
||||
|
||||
/// copy constructor
|
||||
BandMatrix::BandMatrix(const BandMatrix& y) {
|
||||
m_n = y.m_n;
|
||||
m_kl = y.m_kl;
|
||||
m_ku = y.m_ku;
|
||||
data = y.data;
|
||||
ludata = y.ludata;
|
||||
m_factored = y.m_factored;
|
||||
m_ipiv = y.m_ipiv;
|
||||
}
|
||||
//====================================================================================================================
|
||||
/*
|
||||
* Multiply b*A and write result to \c prod.
|
||||
*/
|
||||
void BandMatrix::leftMult(const doublereal * const b, doublereal * const prod) const {
|
||||
int nc = nColumns();
|
||||
doublereal sum = 0.0;
|
||||
for (int n = 0; n < nc; n++) {
|
||||
sum = 0.0;
|
||||
for (int i = n - m_ku; i <= n + m_kl; i++) {
|
||||
if (i >= 0 && i < m_n)
|
||||
sum += _value(i,n) * b[i];
|
||||
}
|
||||
prod[n] = sum;
|
||||
}
|
||||
}
|
||||
//====================================================================================================================
|
||||
/*
|
||||
* Perform an LU decomposition. LAPACK routine DGBTRF is used.
|
||||
* The factorization is saved in ludata.
|
||||
*/
|
||||
int BandMatrix::factor() {
|
||||
int info=0;
|
||||
copy(data.begin(), data.end(), ludata.begin());
|
||||
ct_dgbtrf(nRows(), nColumns(), nSubDiagonals(), nSuperDiagonals(),
|
||||
DATA_PTR(ludata), ldim(), DATA_PTR(ipiv()), info);
|
||||
|
||||
BandMatrix& BandMatrix::operator=(const BandMatrix& y) {
|
||||
if (&y == this) return *this;
|
||||
m_n = y.m_n;
|
||||
m_kl = y.m_kl;
|
||||
m_ku = y.m_ku;
|
||||
m_ipiv = y.m_ipiv;
|
||||
data = y.data;
|
||||
ludata = y.ludata;
|
||||
m_factored = y.m_factored;
|
||||
return *this;
|
||||
// if info = 0, LU decomp succeeded.
|
||||
if (info == 0) {
|
||||
m_factored = true;
|
||||
} else {
|
||||
m_factored = false;
|
||||
ofstream fout("bandmatrix.csv");
|
||||
fout << *this << endl;
|
||||
fout.close();
|
||||
}
|
||||
return info;
|
||||
}
|
||||
//====================================================================================================================
|
||||
int BandMatrix::solve(int n, const doublereal * const b, doublereal * const x) {
|
||||
copy(b, b+n, x);
|
||||
return solve(n, x);
|
||||
}
|
||||
//====================================================================================================================
|
||||
int BandMatrix::solve(int n, doublereal* b) {
|
||||
int info = 0;
|
||||
if (!m_factored) info = factor();
|
||||
if (info == 0)
|
||||
ct_dgbtrs(ctlapack::NoTranspose, nColumns(), nSubDiagonals(),
|
||||
nSuperDiagonals(), 1, DATA_PTR(ludata), ldim(),
|
||||
DATA_PTR(ipiv()), b, nColumns(), info);
|
||||
|
||||
void BandMatrix::resize(int n, int kl, int ku, doublereal v) {
|
||||
m_n = n;
|
||||
m_kl = kl;
|
||||
m_ku = ku;
|
||||
data.resize(n*(2*kl + ku + 1));
|
||||
ludata.resize(n*(2*kl + ku + 1));
|
||||
m_ipiv.resize(m_n);
|
||||
fill(data.begin(), data.end(), v);
|
||||
fill(data.begin(), data.end(), 0.0);
|
||||
m_factored = false;
|
||||
// error handling
|
||||
if (info != 0) {
|
||||
ofstream fout("bandmatrix.csv");
|
||||
fout << *this << endl;
|
||||
fout.close();
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Multiply A*b and write result to \c prod.
|
||||
*/
|
||||
void BandMatrix::mult(const double* b, double* prod) const {
|
||||
int nr = rows();
|
||||
int m, j;
|
||||
double sum = 0.0;
|
||||
for (m = 0; m < nr; m++) {
|
||||
sum = 0.0;
|
||||
for (j = m - m_kl; j <= m + m_ku; j++) {
|
||||
if (j >= 0 && j < m_n)
|
||||
sum += _value(m,j)*b[j];
|
||||
}
|
||||
prod[m] = sum;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Multiply b*A and write result to \c prod.
|
||||
*/
|
||||
void BandMatrix::leftMult(const double* b, double* prod) const {
|
||||
int nc = columns();
|
||||
int n, i;
|
||||
double sum = 0.0;
|
||||
for (n = 0; n < nc; n++) {
|
||||
sum = 0.0;
|
||||
for (i = n - m_ku; i <= n + m_kl; i++) {
|
||||
if (i >= 0 && i < m_n)
|
||||
sum += _value(i,n)*b[i];
|
||||
}
|
||||
prod[n] = sum;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Perform an LU decomposition. LAPACK routine DGBTRF is used.
|
||||
* The factorization is saved in ludata.
|
||||
*/
|
||||
int BandMatrix::factor() {
|
||||
int info=0;
|
||||
copy(data.begin(), data.end(), ludata.begin());
|
||||
ct_dgbtrf(rows(), columns(), nSubDiagonals(), nSuperDiagonals(),
|
||||
DATA_PTR(ludata), ldim(), DATA_PTR(ipiv()), info);
|
||||
|
||||
// if info = 0, LU decomp succeeded.
|
||||
if (info == 0) {
|
||||
m_factored = true;
|
||||
}
|
||||
else {
|
||||
m_factored = false;
|
||||
ofstream fout("bandmatrix.csv");
|
||||
fout << *this << endl;
|
||||
fout.close();
|
||||
}
|
||||
return info;
|
||||
}
|
||||
|
||||
|
||||
|
||||
int BandMatrix::solve(int n, const doublereal* b, doublereal* x) {
|
||||
copy(b, b+n, x);
|
||||
return solve(n, x);
|
||||
}
|
||||
|
||||
int BandMatrix::solve(int n, doublereal* b) {
|
||||
int info = 0;
|
||||
if (!m_factored) info = factor();
|
||||
if (info == 0)
|
||||
ct_dgbtrs(ctlapack::NoTranspose, columns(), nSubDiagonals(),
|
||||
nSuperDiagonals(), 1, DATA_PTR(ludata), ldim(),
|
||||
DATA_PTR(ipiv()), b, columns(), info);
|
||||
|
||||
// error handling
|
||||
if (info != 0) {
|
||||
ofstream fout("bandmatrix.csv");
|
||||
fout << *this << endl;
|
||||
fout.close();
|
||||
}
|
||||
return info;
|
||||
}
|
||||
|
||||
ostream& operator<<(ostream& s, const BandMatrix& m) {
|
||||
int nr = m.rows();
|
||||
int nc = m.columns();
|
||||
int i,j;
|
||||
for (i = 0; i < nr; i++) {
|
||||
for (j = 0; j < nc; j++) {
|
||||
s << m(i,j) << ", ";
|
||||
}
|
||||
s << endl;
|
||||
}
|
||||
return s;
|
||||
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) {
|
||||
int nr = m.nRows();
|
||||
int nc = m.nColumns();
|
||||
for (int i = 0; i < nr; i++) {
|
||||
for (int j = 0; j < nc; j++) {
|
||||
s << m(i,j) << ", ";
|
||||
}
|
||||
s << endl;
|
||||
}
|
||||
return s;
|
||||
}
|
||||
//====================================================================================================================
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -23,130 +23,258 @@
|
|||
|
||||
namespace Cantera {
|
||||
|
||||
/**
|
||||
* A class for banded matrices.
|
||||
/**
|
||||
* A class for banded matrices. This class has matrix inversion processes.
|
||||
* The class is based upon the LAPACK banded storage matrix format.
|
||||
*/
|
||||
class BandMatrix {
|
||||
|
||||
public:
|
||||
|
||||
//! Base Constructor
|
||||
/*!
|
||||
* * Create an \c 0 by \c 0 matrix, and initialize all elements to \c 0.
|
||||
*/
|
||||
class BandMatrix {
|
||||
BandMatrix();
|
||||
|
||||
public:
|
||||
//! Creates a banded matrix and sets all elements to zero
|
||||
/*!
|
||||
* Create an \c n by \c n banded matrix, and initialize all elements to \c v.
|
||||
*
|
||||
* @param n size of the square matrix
|
||||
* @param kl band size on the lower portion of the matrix
|
||||
* @param ku band size on the upper portion of the matrix
|
||||
* @param v intial value of all matrix components.
|
||||
*/
|
||||
BandMatrix(int n, int kl, int ku, doublereal v = 0.0);
|
||||
|
||||
BandMatrix();
|
||||
BandMatrix(int n, int kl, int ku, doublereal v = 0.0);
|
||||
//! Copy constructor
|
||||
/*!
|
||||
* @param y Matrix to be copied
|
||||
*/
|
||||
BandMatrix(const BandMatrix& y);
|
||||
|
||||
/// copy constructor
|
||||
BandMatrix(const BandMatrix& y);
|
||||
//! Destructor. Does nothing.
|
||||
virtual ~BandMatrix();
|
||||
|
||||
/// Destructor. Does nothing.
|
||||
virtual ~BandMatrix(){}
|
||||
//! assignment operator
|
||||
/*!
|
||||
* @param y reference to the matrix to be copied
|
||||
*/
|
||||
BandMatrix& operator=(const BandMatrix& y);
|
||||
|
||||
/// assignment.
|
||||
BandMatrix& operator=(const BandMatrix& y);
|
||||
//! Resize the matrix problem
|
||||
/*!
|
||||
* All data is lost
|
||||
*
|
||||
* @param n size of the square matrix
|
||||
* @param kl band size on the lower portion of the matrix
|
||||
* @param ku band size on the upper portion of the matrix
|
||||
* @param v intial value of all matrix components.
|
||||
*/
|
||||
void resize(int n, int kl, int ku, doublereal v = 0.0);
|
||||
|
||||
void resize(int n, int kl, int ku, doublereal v = 0.0);
|
||||
//! Fill or zero the matrix
|
||||
/*!
|
||||
* @param v Fill value, defaults to zero.
|
||||
*/
|
||||
void bfill(doublereal v = 0.0);
|
||||
|
||||
void bfill(doublereal v) {
|
||||
std::fill(data.begin(), data.end(), v);
|
||||
m_factored = false;
|
||||
}
|
||||
|
||||
doublereal& operator()( int i, int j) {
|
||||
return value(i,j);
|
||||
}
|
||||
|
||||
doublereal operator() ( int i, int j) const {
|
||||
return value(i,j);
|
||||
}
|
||||
|
||||
/// Return a reference to element (i,j). Since this method may
|
||||
/// alter the element value, it may need to be refactored, so
|
||||
/// the flag m_factored is set to false.
|
||||
doublereal& value( int i, int j) {
|
||||
m_factored = false;
|
||||
if (i < j - m_ku || i > j + m_kl) {
|
||||
m_zero = 0.0;
|
||||
return m_zero;
|
||||
}
|
||||
return data[index(i,j)];
|
||||
}
|
||||
|
||||
/// Return the value of element (i,j). This method does not
|
||||
/// alter the array.
|
||||
doublereal value( int i, int j) const {
|
||||
if (i < j - m_ku || i > j + m_kl) return 0.0;
|
||||
return data[index(i,j)];
|
||||
}
|
||||
|
||||
/// Return the location in the internal 1D array corresponding to
|
||||
/// the (i,j) element in the banded array.
|
||||
int index(int i, int j) const {
|
||||
int rw = m_kl + m_ku + i - j;
|
||||
return (2*m_kl + m_ku + 1)*j + rw;
|
||||
}
|
||||
|
||||
/// Return the value of the (i,j) element for (i,j) within the
|
||||
/// bandwidth. For efficiency, this method does not check that
|
||||
/// (i,j) are within the bandwidth; it is up to the calling
|
||||
/// program to insure that this is true.
|
||||
doublereal _value(int i, int j) const {
|
||||
return data[index(i,j)];
|
||||
}
|
||||
|
||||
/// Number of rows
|
||||
int nRows() const { return m_n; }
|
||||
/// @deprecated Redundant.
|
||||
int rows() const { return m_n; }
|
||||
|
||||
/// Number of columns
|
||||
int nColumns() const { return m_n; }
|
||||
/// @deprecated Redundant.
|
||||
int columns() const { return m_n; }
|
||||
|
||||
/// Number of subdiagonals
|
||||
int nSubDiagonals() const { return m_kl; }
|
||||
|
||||
/// Number of superdiagonals
|
||||
int nSuperDiagonals() const { return m_ku; }
|
||||
|
||||
int ldim() const { return 2*m_kl + m_ku + 1; }
|
||||
vector_int& ipiv() { return m_ipiv; }
|
||||
|
||||
/// Multiply A*b and write result to prod.
|
||||
void mult(const double* b, double* prod) const;
|
||||
|
||||
/// Multiply b*A and write result to prod.
|
||||
void leftMult(const double* b, double* prod) const;
|
||||
|
||||
int factor();
|
||||
|
||||
//void solve(const vector_fp& b, vector_fp& x);
|
||||
|
||||
int solve(int n, const doublereal* b, doublereal* x);
|
||||
int solve(int n, doublereal* b);
|
||||
|
||||
vector_fp::iterator begin() {
|
||||
m_factored = false;
|
||||
return data.begin();
|
||||
}
|
||||
vector_fp::iterator end() {
|
||||
m_factored = false;
|
||||
return data.end();
|
||||
}
|
||||
vector_fp::const_iterator begin() const { return data.begin(); }
|
||||
vector_fp::const_iterator end() const { return data.end(); }
|
||||
|
||||
protected:
|
||||
vector_fp data;
|
||||
vector_fp ludata;
|
||||
bool m_factored;
|
||||
//! Index into the (i,j) element
|
||||
/*!
|
||||
* @param i row
|
||||
* @param j column
|
||||
*
|
||||
* Returns a changeable reference to the matrix entry
|
||||
*/
|
||||
doublereal& operator()(int i, int j);
|
||||
|
||||
|
||||
int m_n, m_kl, m_ku;
|
||||
doublereal m_zero;
|
||||
vector_int m_ipiv;
|
||||
//! Constant Index into the (i,j) element
|
||||
/*!
|
||||
* @param i row
|
||||
* @param j column
|
||||
*
|
||||
* Returns an unchangeable reference to the matrix entry
|
||||
*/
|
||||
doublereal operator() (int i, int j) const;
|
||||
|
||||
};
|
||||
//! Return a changeable reference to element (i,j).
|
||||
/*!
|
||||
* Since this method may alter the element value, it may need to be refactored, so
|
||||
* the flag m_factored is set to false.
|
||||
*
|
||||
* @param i row
|
||||
* @param j column
|
||||
*
|
||||
* @return Returns a reference to the value of the matrix entry
|
||||
*/
|
||||
doublereal& value( int i, int j);
|
||||
|
||||
std::ostream& operator<<(std::ostream& s, const BandMatrix& m);
|
||||
|
||||
//! Return the value of element (i,j).
|
||||
/*!
|
||||
* This method does not alter the array.
|
||||
* @param i row
|
||||
* @param j column
|
||||
*
|
||||
* @return Returns the value of the matrix entry
|
||||
*/
|
||||
doublereal value( int i, int j) const;
|
||||
|
||||
//! Returns the location in the internal 1D array corresponding to the (i,j) element in the banded array
|
||||
/*!
|
||||
* @param i row
|
||||
* @param j column
|
||||
*
|
||||
* @return Returns the index of the matrix entry
|
||||
*/
|
||||
int index(int i, int j) const;
|
||||
|
||||
//! Return the value of the (i,j) element for (i,j) within the bandwidth.
|
||||
/*!
|
||||
* For efficiency, this method does not check that (i,j) are within the bandwidth; it is up to the calling
|
||||
* program to insure that this is true.
|
||||
*
|
||||
* @param i row
|
||||
* @param j column
|
||||
*
|
||||
* @return Returns the value of the matrix entry
|
||||
*/
|
||||
doublereal _value(int i, int j) const;
|
||||
|
||||
//! Returns the number of rows
|
||||
int nRows() const;
|
||||
|
||||
//! Number of columns
|
||||
int nColumns() const;
|
||||
|
||||
//! Number of subdiagonals
|
||||
int nSubDiagonals() const;
|
||||
|
||||
//! Number of superdiagonals
|
||||
int nSuperDiagonals() const;
|
||||
|
||||
//! Return the number of rows of storage needed for the band storage
|
||||
int 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 multiplcation
|
||||
* @param prod OUTPUT vector to receive the result
|
||||
*/
|
||||
void mult(const doublereal * const b, doublereal * const prod) const;
|
||||
|
||||
//! Multiply b*A and write result to prod.
|
||||
/*!
|
||||
* @param b Vector to do the lh multiplcation
|
||||
* @param prod OUTPUT vector to receive the result
|
||||
*/
|
||||
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.
|
||||
* 0 indicates a success
|
||||
* ~0 Some error occurred, see the LAPACK documentation
|
||||
*/
|
||||
int factor();
|
||||
|
||||
|
||||
//! Solve the matrix problem Ax = b
|
||||
/*!
|
||||
* @param n size of the matrix
|
||||
* @param b INPUT rhs of the problem
|
||||
* @param x OUTPUT solution to the problem
|
||||
*
|
||||
* @return Return a success flag
|
||||
* 0 indicates a success
|
||||
* ~0 Some error occurred, see the LAPACK documentation
|
||||
*/
|
||||
int solve(int n, const doublereal * const b, doublereal * const x);
|
||||
|
||||
//! Solve the matrix problem Ax = b
|
||||
/*!
|
||||
* @param n size of the matrix
|
||||
* @param b INPUT rhs of the problem
|
||||
* OUTPUT solution to the problem
|
||||
*
|
||||
* @return Return a success flag
|
||||
* 0 indicates a success
|
||||
* ~0 Some error occurred, see the LAPACK documentation
|
||||
*/
|
||||
int solve(int n, doublereal * const b);
|
||||
|
||||
|
||||
//! Returns an iterator for the start of the band storage data
|
||||
/*!
|
||||
* Iterator points to the beginning of the data, and it is changeable.
|
||||
*/
|
||||
vector_fp::iterator begin();
|
||||
|
||||
//! Returns an iterator for the end of the band storage data
|
||||
/*!
|
||||
* Iterator points to the end of the data, and it is changeable.
|
||||
*/
|
||||
vector_fp::iterator end();
|
||||
|
||||
//! Returns a const iterator for the start of the band storage data
|
||||
/*!
|
||||
* Iterator points to the beginning of the data, and it is not changeable.
|
||||
*/
|
||||
vector_fp::const_iterator begin() const;
|
||||
|
||||
//! Returns a const iterator for the end of the band storage data
|
||||
/*!
|
||||
* Iterator points to the end of the data, and it is not changeable.
|
||||
*/
|
||||
vector_fp::const_iterator end() const;
|
||||
|
||||
protected:
|
||||
|
||||
//! Matrix data
|
||||
vector_fp data;
|
||||
|
||||
//! Factorized data
|
||||
vector_fp ludata;
|
||||
|
||||
//! Boolean indicating whether the matrix is factored
|
||||
bool m_factored;
|
||||
|
||||
//! Number of rows and columns of the matrix
|
||||
int m_n;
|
||||
|
||||
//! Number of subdiagonals of the matrix
|
||||
int m_kl;
|
||||
|
||||
//! Number of super diagonals of the matrix
|
||||
int m_ku;
|
||||
|
||||
//! value of zero
|
||||
doublereal m_zero;
|
||||
|
||||
//! Pivot vector
|
||||
vector_int m_ipiv;
|
||||
|
||||
};
|
||||
|
||||
//! Utility routine to print out the matrix
|
||||
/*!
|
||||
* @param s ostream to print the matrix out to
|
||||
* @param m Matrix to be printed
|
||||
*
|
||||
* @return Returns a reference to the ostream
|
||||
*/
|
||||
std::ostream& operator<<(std::ostream& s, const BandMatrix& m);
|
||||
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -716,8 +716,8 @@ namespace Cantera {
|
|||
* @param ydot_comm Current value of the time derivative of the solution vector
|
||||
* @param time_curr current value of the time
|
||||
*/
|
||||
void NonlinearSolver::scaleMatrix(SquareMatrix& jac, double* y_comm, double* ydot_comm, doublereal time_curr,
|
||||
int num_newt_its)
|
||||
void NonlinearSolver::scaleMatrix(SquareMatrix& jac, doublereal * const y_comm, doublereal * const ydot_comm,
|
||||
doublereal time_curr, int num_newt_its)
|
||||
{
|
||||
int irow, jcol;
|
||||
/*
|
||||
|
|
@ -2319,10 +2319,10 @@ namespace Cantera {
|
|||
* NSOLN_RETN_FAIL_DAMPSTEP
|
||||
* Unsuccessful step. We can not find a damping factor that is suitable.
|
||||
*/
|
||||
int NonlinearSolver::dampStep(const doublereal time_curr, const doublereal * y_n_curr,
|
||||
const doublereal * ydot_n_curr, doublereal * const step_1,
|
||||
doublereal * const y_n_1, double* const ydot_n_1, doublereal * const step_2,
|
||||
double& stepNorm_2, SquareMatrix& jac, bool writetitle, int& num_backtracks)
|
||||
int NonlinearSolver::dampStep(const doublereal time_curr, const doublereal * const y_n_curr,
|
||||
const doublereal * const ydot_n_curr, doublereal * const step_1,
|
||||
doublereal * const y_n_1, doublereal * const ydot_n_1, doublereal * const step_2,
|
||||
doublereal & stepNorm_2, SquareMatrix& jac, bool writetitle, int& num_backtracks)
|
||||
{
|
||||
int j, m;
|
||||
int info = 0;
|
||||
|
|
@ -2729,7 +2729,7 @@ namespace Cantera {
|
|||
*/
|
||||
int NonlinearSolver::decideStep(const doublereal time_curr, int leg, doublereal alpha, const doublereal * const y_n_curr,
|
||||
const doublereal * const ydot_n_curr, const std::vector<doublereal> & step_1,
|
||||
const doublereal * const y_n_1, const double* const ydot_n_1, doublereal trustDeltaOld)
|
||||
const doublereal * const y_n_1, const doublereal * const ydot_n_1, doublereal trustDeltaOld)
|
||||
{
|
||||
int retn = 2;
|
||||
bool goodStep = false;
|
||||
|
|
@ -2906,8 +2906,8 @@ namespace Cantera {
|
|||
* @return A positive value indicates a successful convergence
|
||||
* -1 Failed convergence
|
||||
*/
|
||||
int NonlinearSolver::solve_nonlinear_problem(int SolnType, double* y_comm, double* ydot_comm, doublereal CJ,
|
||||
doublereal time_curr, SquareMatrix& jac,
|
||||
int NonlinearSolver::solve_nonlinear_problem(int SolnType, doublereal * const y_comm, doublereal * const ydot_comm,
|
||||
doublereal CJ, doublereal time_curr, SquareMatrix& jac,
|
||||
int &num_newt_its, int &num_linear_solves,
|
||||
int &num_backtracks, int loglevelInput)
|
||||
{
|
||||
|
|
|
|||
|
|
@ -311,6 +311,7 @@ namespace Cantera {
|
|||
*/
|
||||
void setDefaultDeltaBoundsMagnitudes();
|
||||
|
||||
//! Adjust the step minimums
|
||||
void adjustUpStepMinimums();
|
||||
|
||||
//! Set the delta Bounds magnitudes by hand
|
||||
|
|
@ -506,10 +507,10 @@ namespace Cantera {
|
|||
*
|
||||
* @return returns an integer indicating what happened.
|
||||
*/
|
||||
int dampStep(const doublereal time_curr, const double* y_n_curr,
|
||||
const doublereal *ydot_n_curr, double * const step_1,
|
||||
double* const y_n_1, double* const ydot_n_1, double* step_2,
|
||||
double& stepNorm_2, SquareMatrix& jac, bool writetitle,
|
||||
int dampStep(const doublereal time_curr, const doublereal * const y_n_curr,
|
||||
const doublereal * const ydot_n_curr, doublereal * const step_1,
|
||||
doublereal * const y_n_1, doublereal * const ydot_n_1, doublereal * step_2,
|
||||
doublereal & stepNorm_2, SquareMatrix& jac, bool writetitle,
|
||||
int& num_backtracks);
|
||||
|
||||
//! Find the solution to F(X) = 0 by damped Newton iteration.
|
||||
|
|
@ -539,7 +540,7 @@ namespace Cantera {
|
|||
* @return A positive value indicates a successful convergence
|
||||
* -1 Failed convergence
|
||||
*/
|
||||
int solve_nonlinear_problem(int SolnType, double* y_comm,double* ydot_comm, doublereal CJ,
|
||||
int solve_nonlinear_problem(int SolnType, doublereal * const y_comm, doublereal * const ydot_comm, doublereal CJ,
|
||||
doublereal time_curr, SquareMatrix& jac,int &num_newt_its,
|
||||
int &num_linear_solves, int &num_backtracks, int loglevelInput);
|
||||
|
||||
|
|
@ -586,8 +587,10 @@ namespace Cantera {
|
|||
* @param y_comm Current value of the solution vector
|
||||
* @param ydot_comm Current value of the time derivative of the solution vector
|
||||
* @param time_curr current value of the time
|
||||
* @param num_newt_its Current value of the number of newt its
|
||||
*/
|
||||
void scaleMatrix(SquareMatrix& jac, double* y_comm, double* ydot_comm, doublereal time_curr, int num_newt_its);
|
||||
void scaleMatrix(SquareMatrix& jac, doublereal * const y_comm, doublereal * const ydot_comm,
|
||||
doublereal time_curr, int num_newt_its);
|
||||
|
||||
//! Print solution norm contribution
|
||||
/*!
|
||||
|
|
@ -802,11 +805,11 @@ namespace Cantera {
|
|||
* @param time_curr INPUT Current value of the time
|
||||
* @param leg INPUT Leg of the dogleg that we are on
|
||||
* @param alpha INPUT Distance down that leg that we are on
|
||||
* @param y0 INPUT Current value of the solution vector
|
||||
* @param ydot0 INPUT Current value of the derivative of the solution vector
|
||||
* @param step0 INPUT Trial step
|
||||
* @param y1 OUTPUT Solution values at the conditions which are evalulated for success
|
||||
* @param ydot1 OUTPUT Time derivates of solution at the conditions which are evalulated for success
|
||||
* @param y_n_curr INPUT Current value of the solution vector
|
||||
* @param ydot_n_curr INPUT Current value of the derivative of the solution vector
|
||||
* @param step_1 INPUT Trial step
|
||||
* @param y_n_1 OUTPUT Solution values at the conditions which are evalulated for success
|
||||
* @param ydot_n_1 OUTPUT Time derivates of solution at the conditions which are evalulated for success
|
||||
* @param trustDeltaOld INPUT Value of the trust length at the old conditions
|
||||
*
|
||||
*
|
||||
|
|
@ -819,10 +822,10 @@ namespace Cantera {
|
|||
* -2 Current value of the solution vector caused a residual error in its evaluation.
|
||||
* Step is a failure, and the step size must be reduced in order to proceed further.
|
||||
*/
|
||||
int decideStep(const doublereal time_curr, int leg, doublereal alpha, const doublereal * const y0,
|
||||
const doublereal * const ydot0,
|
||||
const std::vector<doublereal> & step0,
|
||||
const doublereal * const y1, const doublereal * const ydot1, doublereal trustDeltaOld);
|
||||
int decideStep(const doublereal time_curr, int leg, doublereal alpha, const doublereal * const y_n_curr,
|
||||
const doublereal * const ydot_n_curr,
|
||||
const std::vector<doublereal> & step_1,
|
||||
const doublereal * const y_n_1, const doublereal * const ydot_n_1, doublereal trustDeltaOld);
|
||||
|
||||
//! Calculated the expected residual along the double dogleg curve.
|
||||
/*!
|
||||
|
|
@ -923,6 +926,7 @@ namespace Cantera {
|
|||
//! Value of the solution time derivative at the new point that is to be considered
|
||||
std::vector<doublereal> m_ydot_n_1;
|
||||
|
||||
//! Value of the step to be taken in the solution
|
||||
std::vector<doublereal> m_step_1;
|
||||
|
||||
//! Vector of column scaling factors
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue