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

Added in the estimation of condition numbers for QR factored matrices.

Browse source
This commit is contained in:
Harry Moffat 2011-04-28 20:56:41 +00:00
parent 24ac9bd133
commit 701dd1e1cd
3 changed files with 95 additions and 8 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -112,8 +112,8 @@ namespace Cantera {
}
return info;
}
/**
//====================================================================================================================
/*
* Set all entries to zero
*/
void SquareMatrix::zero() {
@ -191,6 +191,8 @@ namespace Cantera {
if (lworkOpt > lwork) {
work.resize(lworkOpt);
}
return info;
}
//=====================================================================================================================
@ -250,9 +252,33 @@ namespace Cantera {
return info;
}
//=====================================================================================================================
doublereal SquareMatrix::rcondQR() {
if ((int) iwork_.size() < m_nrows) {
iwork_.resize(m_nrows);
}
if ((int) work.size() <3 * m_nrows) {
work.resize(3 * m_nrows);
}
doublereal rcond = 0.0;
if (m_factored != 2) {
throw CELapackError("SquareMatrix::rcondQR()", "matrix isn't factored correctly");
}
int rinfo;
rcond = ct_dtrcon(0, ctlapack::UpperTriangular, 0, m_nrows, &(*(begin())), m_nrows, DATA_PTR(work),
DATA_PTR(iwork_), rinfo);
if (rinfo != 0) {
if (m_printLevel) {
writelogf("SquareMatrix::rcondQR(): DTRCON returned INFO = %d\n", rinfo);
}
if (! m_useReturnErrorCode) {
throw CELapackError("SquareMatrix::rcondQR()", "DTRCON returned INFO = " + int2str(rinfo));
}
}
return rcond;
}
//=====================================================================================================================
}

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

@ -96,6 +96,14 @@ namespace Cantera {
*/
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
*/
doublereal rcondQR();
//! Solves the linear problem Ax=b using the QR algorithm returning x in the b spot
/*!
* @param b RHS to be solved.
@ -121,6 +129,9 @@ namespace Cantera {
//! Work vector for QR algorithm
vector_fp work;
//! Integer work vector for QR algorithms
std::vector<int> iwork_;
public:
//! Use the QR algorithm to factor and invert the matrix
int useQR_;

56
Cantera/src/numerics/ctlapack.h
View file

@ -38,6 +38,7 @@
#define _DGEQRF_ dgeqrf
#define _DORMQR_ dormqr
#define _DTRTRS_ dtrtrs
#define _DTRCON_ dtrcon
#else
@ -55,6 +56,7 @@
#define _DGEQRF_ dgeqrf_
#define _DORMQR_ dormqr_
#define _DTRTRS_ dtrtrs_
#define _DTRCON_ dtrcon_
#endif
@ -171,6 +173,19 @@ extern "C" {
#endif
#ifdef LAPACK_FTN_STRING_LEN_AT_END
int _DTRCON_(const char* norm, const char* uplo, const char *diag, const integer* n,
doublereal* a, const integer* lda, const doublereal *rcond,
doublereal* work, const integer* iwork, integer *info, ftnlen nosize,
ftnlen upsize, ftnlen disize);
#else
int _DTRCON_(const char* norm, ftnlen nosize, const char* uplo, ftnlen upsize, const char *diag,
ftnlen disize, const integer* n, doublereal* a, const integer* lda, const doublereal *rcond,
doublereal* work, const integer* iwork, integer *info);
#endif
}
//#endif
@ -304,7 +319,7 @@ namespace Cantera {
//_DSCAL_(&f_n, &da, dx, &f_incx);
cblas_dscal(n, da, dx, incx);
}
//====================================================================================================================
inline void ct_dgeqrf(int m, int n, doublereal* a, int lda, doublereal *tau,
doublereal* work, int lwork, int &info) {
integer f_m = m;
@ -315,7 +330,7 @@ namespace Cantera {
_DGEQRF_(&f_m, &f_n, a, &f_lda, tau, work, &f_lwork, &f_info);
info = f_info;
}
//====================================================================================================================
inline void ct_dormqr(ctlapack::side_t rlside, ctlapack::transpose_t trans, int m,
int n, int k, doublereal* a, int lda, doublereal *tau, doublereal *c, int ldc,
doublereal *work, int lwork, int &info) {
@ -340,7 +355,7 @@ namespace Cantera {
#endif
info = f_info;
}
//====================================================================================================================
inline void ct_dtrtrs(ctlapack::upperlower_t uplot, ctlapack::transpose_t trans, const char *diag,
int n, int nrhs, doublereal* a, int lda, doublereal *b, int ldb, int &info) {
char uplo = upper_lower[uplot];
@ -366,7 +381,42 @@ namespace Cantera {
#endif
info = f_info;
}
//====================================================================================================================
//!
/*!
* @param work Must be dimensioned equal to greater than 3N
* @param iwork Must be dimensioned equal to or greater than N
*/
inline doublereal ct_dtrcon(const char *norm, ctlapack::upperlower_t uplot, const char *diag,
int n, doublereal* a, int lda, doublereal *work, int *iwork, int &info) {
char uplo = upper_lower[uplot];
char dd = 'N';
if (diag) {
dd = diag[0];
}
char nn = '1';
if (norm) {
nn = norm[0];
}
integer f_n = n;
integer f_lda = lda;
integer f_info = info;
doublereal rcond;
#ifdef NO_FTN_STRING_LEN_AT_END
_DTRCON_(&nn, &uplo, &dd, &f_n, a, &f_lda, &rcond, work, iwork, &f_info);
#else
ftnlen trsize = 1;
#ifdef LAPACK_FTN_STRING_LEN_AT_END
_DTRCON_(&nn, &uplo, &dd, &f_n, a, &f_lda, &rcond, work, iwork, &f_info, trsize, trsize, trsize);
#else
_DTRCON_(&nn, trsize, &uplo, trsize, &dd, trsize, &f_n, a, &f_lda, &rcond work, iwork, &f_info);
#endif
#endif
info = f_info;
return rcond;
}
//====================================================================================================================
}
#endif