Doxygen update to DenseMatrix
This commit is contained in:
Harry Moffat
parent
ec23984da8
commit
a737dd308f
2 changed files with 222 additions and 186 deletions
|
|
@ -16,134 +16,165 @@
|
|||
#include "stringUtils.h"
|
||||
|
||||
namespace Cantera {
|
||||
|
||||
/// assignment.
|
||||
DenseMatrix& DenseMatrix::operator=(const DenseMatrix& y) {
|
||||
if (&y == this) return *this;
|
||||
Array2D::operator=(y);
|
||||
m_ipiv = y.ipiv();
|
||||
return *this;
|
||||
}
|
||||
|
||||
void DenseMatrix::resize(int n, int m, doublereal v) {
|
||||
Array2D::resize(n,m,v);
|
||||
m_ipiv.resize( max(n,m) );
|
||||
}
|
||||
|
||||
void DenseMatrix::mult(const double* b, double* prod) const {
|
||||
ct_dgemv(ctlapack::ColMajor, ctlapack::NoTranspose,
|
||||
static_cast<int>(nRows()),
|
||||
static_cast<int>(nRows()), 1.0, ptrColumn(0), //begin(),
|
||||
static_cast<int>(nRows()), b, 1, 0.0, prod, 1);
|
||||
}
|
||||
|
||||
void DenseMatrix::leftMult(const double* b, double* prod) const {
|
||||
int nc = static_cast<int>(nColumns());
|
||||
int nr = static_cast<int>(nRows());
|
||||
int n, i;
|
||||
double sum = 0.0;
|
||||
for (n = 0; n < nc; n++) {
|
||||
sum = 0.0;
|
||||
for (i = 0; i < nr; i++) {
|
||||
sum += value(i,n)*b[i];
|
||||
}
|
||||
prod[n] = sum;
|
||||
}
|
||||
}
|
||||
|
||||
int solve(DenseMatrix& A, double* b) {
|
||||
int info=0;
|
||||
ct_dgetrf(static_cast<int>(A.nRows()),
|
||||
static_cast<int>(A.nColumns()), A.ptrColumn(0), //begin(),
|
||||
static_cast<int>(A.nRows()), &A.ipiv()[0], info);
|
||||
if (info != 0)
|
||||
throw CanteraError("DenseMatrix::solve",
|
||||
"DGETRF returned INFO = "+int2str(info));
|
||||
ct_dgetrs(ctlapack::NoTranspose,
|
||||
static_cast<int>(A.nRows()), 1, A.ptrColumn(0), //begin(),
|
||||
static_cast<int>(A.nRows()),
|
||||
&A.ipiv()[0], b,
|
||||
static_cast<int>(A.nColumns()), info);
|
||||
if (info != 0)
|
||||
throw CanteraError("DenseMatrix::solve",
|
||||
"DGETRS returned INFO = "+int2str(info));
|
||||
return 0;
|
||||
//====================================================================================================================
|
||||
// Default Constructor
|
||||
DenseMatrix::DenseMatrix()
|
||||
{
|
||||
}
|
||||
//====================================================================================================================
|
||||
/*
|
||||
* Constructor. Create an \c n by \c m matrix, and initialize
|
||||
* all elements to \c v.
|
||||
*/
|
||||
DenseMatrix::DenseMatrix(int n, int m, doublereal v) :
|
||||
Array2D(n, m, v)
|
||||
{
|
||||
m_ipiv.resize(max(n, m));
|
||||
}
|
||||
//====================================================================================================================
|
||||
// copy constructor
|
||||
/*
|
||||
* @param y Object to be copied
|
||||
*/
|
||||
DenseMatrix::DenseMatrix(const DenseMatrix& y) :
|
||||
Array2D(y)
|
||||
{
|
||||
m_ipiv = y.ipiv();
|
||||
}
|
||||
//====================================================================================================================
|
||||
// assignment
|
||||
DenseMatrix& DenseMatrix::operator=(const DenseMatrix& y) {
|
||||
if (&y == this) return *this;
|
||||
Array2D::operator=(y);
|
||||
m_ipiv = y.ipiv();
|
||||
return *this;
|
||||
}
|
||||
//====================================================================================================================
|
||||
// Destructor. Does nothing.
|
||||
DenseMatrix::~DenseMatrix()
|
||||
{
|
||||
}
|
||||
//====================================================================================================================
|
||||
void DenseMatrix::resize(int n, int m, doublereal v) {
|
||||
Array2D::resize(n,m,v);
|
||||
m_ipiv.resize( max(n,m) );
|
||||
}
|
||||
//====================================================================================================================
|
||||
void DenseMatrix::mult(const double* b, double* prod) const {
|
||||
ct_dgemv(ctlapack::ColMajor, ctlapack::NoTranspose,
|
||||
static_cast<int>(nRows()),
|
||||
static_cast<int>(nRows()), 1.0, ptrColumn(0), //begin(),
|
||||
static_cast<int>(nRows()), b, 1, 0.0, prod, 1);
|
||||
}
|
||||
//====================================================================================================================
|
||||
void DenseMatrix::leftMult(const double* b, double* prod) const {
|
||||
int nc = static_cast<int>(nColumns());
|
||||
int nr = static_cast<int>(nRows());
|
||||
int n, i;
|
||||
double sum = 0.0;
|
||||
for (n = 0; n < nc; n++) {
|
||||
sum = 0.0;
|
||||
for (i = 0; i < nr; i++) {
|
||||
sum += value(i,n)*b[i];
|
||||
}
|
||||
prod[n] = sum;
|
||||
}
|
||||
|
||||
int solve(DenseMatrix& A, DenseMatrix& b) {
|
||||
int info=0;
|
||||
ct_dgetrf(static_cast<int>(A.nRows()),
|
||||
static_cast<int>(A.nColumns()), A.ptrColumn(0),
|
||||
static_cast<int>(A.nRows()), &A.ipiv()[0], info);
|
||||
if (info != 0)
|
||||
throw CanteraError("DenseMatrix::solve",
|
||||
"DGETRF returned INFO = "+int2str(info));
|
||||
ct_dgetrs(ctlapack::NoTranspose, static_cast<int>(A.nRows()),
|
||||
static_cast<int>(b.nColumns()),
|
||||
A.ptrColumn(0), static_cast<int>(A.nRows()),
|
||||
&A.ipiv()[0], b.ptrColumn(0),
|
||||
static_cast<int>(b.nRows()), info);
|
||||
if (info != 0)
|
||||
throw CanteraError("DenseMatrix::solve",
|
||||
"DGETRS returned INFO = "+int2str(info));
|
||||
return 0;
|
||||
}
|
||||
|
||||
}
|
||||
//====================================================================================================================
|
||||
int solve(DenseMatrix& A, double* b) {
|
||||
int info=0;
|
||||
ct_dgetrf(static_cast<int>(A.nRows()),
|
||||
static_cast<int>(A.nColumns()), A.ptrColumn(0), //begin(),
|
||||
static_cast<int>(A.nRows()), &A.ipiv()[0], info);
|
||||
if (info != 0)
|
||||
throw CanteraError("DenseMatrix::solve",
|
||||
"DGETRF returned INFO = "+int2str(info));
|
||||
ct_dgetrs(ctlapack::NoTranspose,
|
||||
static_cast<int>(A.nRows()), 1, A.ptrColumn(0), //begin(),
|
||||
static_cast<int>(A.nRows()),
|
||||
&A.ipiv()[0], b,
|
||||
static_cast<int>(A.nColumns()), info);
|
||||
if (info != 0)
|
||||
throw CanteraError("DenseMatrix::solve",
|
||||
"DGETRS returned INFO = "+int2str(info));
|
||||
return 0;
|
||||
}
|
||||
//====================================================================================================================
|
||||
int solve(DenseMatrix& A, DenseMatrix& b) {
|
||||
int info=0;
|
||||
ct_dgetrf(static_cast<int>(A.nRows()),
|
||||
static_cast<int>(A.nColumns()), A.ptrColumn(0),
|
||||
static_cast<int>(A.nRows()), &A.ipiv()[0], info);
|
||||
if (info != 0)
|
||||
throw CanteraError("DenseMatrix::solve",
|
||||
"DGETRF returned INFO = "+int2str(info));
|
||||
ct_dgetrs(ctlapack::NoTranspose, static_cast<int>(A.nRows()),
|
||||
static_cast<int>(b.nColumns()),
|
||||
A.ptrColumn(0), static_cast<int>(A.nRows()),
|
||||
&A.ipiv()[0], b.ptrColumn(0),
|
||||
static_cast<int>(b.nRows()), info);
|
||||
if (info != 0)
|
||||
throw CanteraError("DenseMatrix::solve",
|
||||
"DGETRS returned INFO = "+int2str(info));
|
||||
return 0;
|
||||
}
|
||||
//====================================================================================================================
|
||||
|
||||
#ifdef INCL_LEAST_SQUARES
|
||||
/** @todo fix lwork */
|
||||
int leastSquares(DenseMatrix& A, double* b) {
|
||||
int info = 0;
|
||||
int rank = 0;
|
||||
double rcond = -1.0;
|
||||
// fix this!
|
||||
int lwork = 6000; // 2*(3*min(m,n) + max(2*min(m,n), max(m,n)));
|
||||
vector_fp work(lwork);
|
||||
vector_fp s(min(static_cast<int>(A.nRows()),
|
||||
static_cast<int>(A.nColumns())));
|
||||
ct_dgelss(static_cast<int>(A.nRows()),
|
||||
static_cast<int>(A.nColumns()), 1, A.ptrColumn(0),
|
||||
static_cast<int>(A.nRows()), b,
|
||||
static_cast<int>(A.nColumns()), &s[0], //.begin(),
|
||||
rcond, rank, &work[0], work.size(), info);
|
||||
if (info != 0)
|
||||
throw CanteraError("DenseMatrix::leaseSquares",
|
||||
"DGELSS returned INFO = "+int2str(info));
|
||||
return 0;
|
||||
}
|
||||
/** @todo fix lwork */
|
||||
int leastSquares(DenseMatrix& A, double* b) {
|
||||
int info = 0;
|
||||
int rank = 0;
|
||||
double rcond = -1.0;
|
||||
// fix this!
|
||||
int lwork = 6000; // 2*(3*min(m,n) + max(2*min(m,n), max(m,n)));
|
||||
vector_fp work(lwork);
|
||||
vector_fp s(min(static_cast<int>(A.nRows()),
|
||||
static_cast<int>(A.nColumns())));
|
||||
ct_dgelss(static_cast<int>(A.nRows()),
|
||||
static_cast<int>(A.nColumns()), 1, A.ptrColumn(0),
|
||||
static_cast<int>(A.nRows()), b,
|
||||
static_cast<int>(A.nColumns()), &s[0], //.begin(),
|
||||
rcond, rank, &work[0], work.size(), info);
|
||||
if (info != 0)
|
||||
throw CanteraError("DenseMatrix::leaseSquares",
|
||||
"DGELSS returned INFO = "+int2str(info));
|
||||
return 0;
|
||||
}
|
||||
#endif
|
||||
//====================================================================================================================
|
||||
void multiply(const DenseMatrix& A, const double* b, double* prod) {
|
||||
ct_dgemv(ctlapack::ColMajor, ctlapack::NoTranspose,
|
||||
static_cast<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, int nn) {
|
||||
integer n = (nn > 0 ? nn : static_cast<int>(A.nRows()));
|
||||
int info=0;
|
||||
ct_dgetrf(n, n, A.ptrColumn(0), static_cast<int>(A.nRows()),
|
||||
&A.ipiv()[0], info);
|
||||
if (info != 0)
|
||||
throw CanteraError("invert",
|
||||
"DGETRF returned INFO="+int2str(info));
|
||||
|
||||
void multiply(const DenseMatrix& A, const double* b, double* prod) {
|
||||
ct_dgemv(ctlapack::ColMajor, ctlapack::NoTranspose,
|
||||
static_cast<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, int nn) {
|
||||
integer n = (nn > 0 ? nn : static_cast<int>(A.nRows()));
|
||||
int info=0;
|
||||
ct_dgetrf(n, n, A.ptrColumn(0), static_cast<int>(A.nRows()),
|
||||
&A.ipiv()[0], info);
|
||||
if (info != 0)
|
||||
throw CanteraError("invert",
|
||||
"DGETRF returned INFO="+int2str(info));
|
||||
|
||||
vector_fp work(n);
|
||||
integer lwork = static_cast<int>(work.size());
|
||||
ct_dgetri(n, A.ptrColumn(0), static_cast<int>(A.nRows()),
|
||||
&A.ipiv()[0],
|
||||
&work[0], lwork, info);
|
||||
if (info != 0)
|
||||
throw CanteraError("invert",
|
||||
"DGETRI returned INFO="+int2str(info));
|
||||
return 0;
|
||||
}
|
||||
vector_fp work(n);
|
||||
integer lwork = static_cast<int>(work.size());
|
||||
ct_dgetri(n, A.ptrColumn(0), static_cast<int>(A.nRows()),
|
||||
&A.ipiv()[0],
|
||||
&work[0], lwork, info);
|
||||
if (info != 0)
|
||||
throw CanteraError("invert",
|
||||
"DGETRI returned INFO="+int2str(info));
|
||||
return 0;
|
||||
}
|
||||
//====================================================================================================================
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -22,83 +22,88 @@
|
|||
|
||||
namespace Cantera {
|
||||
|
||||
/**
|
||||
* A class for full (non-sparse) matrices with Fortran-compatible
|
||||
* data storage. Adds matrix operations to class Array2D.
|
||||
/**
|
||||
* A class for full (non-sparse) matrices with Fortran-compatible
|
||||
* data storage. Adds matrix operations to class Array2D.
|
||||
*/
|
||||
class DenseMatrix : public Array2D {
|
||||
|
||||
public:
|
||||
|
||||
//! Default Constructor
|
||||
DenseMatrix();
|
||||
|
||||
/**
|
||||
* Constructor. Create an \c n by \c m matrix, and initialize
|
||||
* all elements to \c v.
|
||||
*/
|
||||
class DenseMatrix : public Array2D {
|
||||
DenseMatrix(int n, int m, doublereal v = 0.0);
|
||||
|
||||
public:
|
||||
//! copy constructor
|
||||
/*!
|
||||
* @param y Object to be copied
|
||||
*/
|
||||
DenseMatrix(const DenseMatrix& y);
|
||||
|
||||
DenseMatrix(){}
|
||||
//! assignment operator
|
||||
/*!
|
||||
* @param y Object to be copied
|
||||
*/
|
||||
DenseMatrix& operator=(const DenseMatrix& y);
|
||||
|
||||
/**
|
||||
* Constructor. Create an \c n by \c m matrix, and initialize
|
||||
* all elements to \c v.
|
||||
*/
|
||||
DenseMatrix(int n, int m, doublereal v = 0.0) : Array2D(n,m,v) {
|
||||
m_ipiv.resize( max(n, m) );
|
||||
}
|
||||
//! Destructor. Does nothing.
|
||||
virtual ~DenseMatrix();
|
||||
|
||||
|
||||
/// copy constructor
|
||||
DenseMatrix(const DenseMatrix& y) : Array2D(y) {
|
||||
m_ipiv = y.ipiv();
|
||||
}
|
||||
|
||||
/// assignment.
|
||||
DenseMatrix& operator=(const DenseMatrix& y);
|
||||
|
||||
void resize(int n, int m, doublereal v = 0.0);
|
||||
|
||||
/// Destructor. Does nothing.
|
||||
virtual ~DenseMatrix(){}
|
||||
|
||||
|
||||
/**
|
||||
* Multiply A*b and write result to \c prod.
|
||||
*/
|
||||
virtual void mult(const double* b, double* prod) const;
|
||||
|
||||
/**
|
||||
* Left-multiply the matrix by transpose(b), and write the
|
||||
* result to prod.
|
||||
*/
|
||||
virtual void leftMult(const double* b, double* prod) const;
|
||||
|
||||
vector_int& ipiv() { return m_ipiv; }
|
||||
const vector_int& ipiv() const { return m_ipiv; }
|
||||
|
||||
protected:
|
||||
|
||||
vector_int m_ipiv;
|
||||
};
|
||||
void resize(int n, int m, doublereal v = 0.0);
|
||||
|
||||
|
||||
|
||||
/**
|
||||
* Solve Ax = b. Array b is overwritten on exit with x.
|
||||
* Multiply A*b and write result to \c prod.
|
||||
*/
|
||||
int solve(DenseMatrix& A, double* b);
|
||||
virtual void mult(const double* b, double* prod) const;
|
||||
|
||||
/** Solve Ax = b for multiple right-hand-side vectors. */
|
||||
int solve(DenseMatrix& A, DenseMatrix& b);
|
||||
/**
|
||||
* Left-multiply the matrix by transpose(b), and write the
|
||||
* result to prod.
|
||||
*/
|
||||
virtual void leftMult(const double* b, double* prod) const;
|
||||
|
||||
vector_int& ipiv() { return m_ipiv; }
|
||||
const vector_int& ipiv() const { return m_ipiv; }
|
||||
|
||||
protected:
|
||||
|
||||
vector_int m_ipiv;
|
||||
};
|
||||
|
||||
|
||||
/**
|
||||
* Solve Ax = b. Array b is overwritten on exit with x.
|
||||
*/
|
||||
int solve(DenseMatrix& A, double* b);
|
||||
|
||||
/** Solve Ax = b for multiple right-hand-side vectors. */
|
||||
int solve(DenseMatrix& A, DenseMatrix& b);
|
||||
|
||||
#ifdef INCL_LEAST_SQUARES
|
||||
/** @todo fix lwork */
|
||||
int leastSquares(DenseMatrix& A, double* b);
|
||||
/** @todo fix lwork */
|
||||
int leastSquares(DenseMatrix& A, double* b);
|
||||
#endif
|
||||
/**
|
||||
* Multiply \c A*b and return the result in \c prod. Uses BLAS
|
||||
* routine DGEMV.
|
||||
*/
|
||||
void multiply(const DenseMatrix& A, const double* b, double* prod);
|
||||
/**
|
||||
* Multiply \c A*b and return the result in \c prod. Uses BLAS
|
||||
* routine DGEMV.
|
||||
*/
|
||||
void multiply(const DenseMatrix& A, const double* b, double* prod);
|
||||
|
||||
void increment(const DenseMatrix& A,
|
||||
const double* b, double* prod);
|
||||
void increment(const DenseMatrix& A,
|
||||
const double* b, double* prod);
|
||||
|
||||
/**
|
||||
* invert A. A is overwritten with A^-1.
|
||||
*/
|
||||
int invert(DenseMatrix& A, int nn=-1);
|
||||
/**
|
||||
* invert A. A is overwritten with A^-1.
|
||||
*/
|
||||
int invert(DenseMatrix& A, int nn=-1);
|
||||
|
||||
}
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue