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

Another iteration to get the QR solve up to speed

Browse source
This commit is contained in:
Harry Moffat 2011-04-06 16:26:32 +00:00
parent fa30f63056
commit 68768e56e7
3 changed files with 68 additions and 4 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -179,6 +179,62 @@ namespace Cantera {
return info;
}
//=====================================================================================================================
/*
* Solve Ax = b. Vector b is overwritten on exit with x.
*/
int SquareMatrix::solveQR(double* b)
{
int info=0;
/*
* Check to see whether the matrix has been factored.
*/
if (!m_factored) {
int retn = factorQR();
if (retn) {
return retn;
}
}
int lwork = work.size();
if (lwork < m_nrows) {
work.resize(8 * m_nrows, 0.0);
lwork = 8 * m_nrows;
}
/*
* Solve the factored system
*/
ct_dormqr(ctlapack::Left, ctlapack::Transpose, m_nrows, 1, m_nrows, &(*(begin())), m_nrows, DATA_PTR(tau), b, m_nrows,
DATA_PTR(work), lwork, info);
if (info != 0) {
if (m_printLevel) {
writelogf("SquareMatrix::solveQR(): DORMQR returned INFO = %d\n", info);
}
if (! m_useReturnErrorCode) {
throw CELapackError("SquareMatrix::solveQR()", "DORMQR returned INFO = " + int2str(info));
}
}
int lworkOpt = work[0];
if (lworkOpt <= lwork) {
work.resize(lworkOpt);
}
char dd = 'N';
ct_dtrtrs(ctlapack::UpperTriangular, ctlapack::NoTranspose, &dd, m_nrows, 1, &(*(begin())), m_nrows, b,
m_nrows, info);
if (info != 0) {
if (m_printLevel) {
writelogf("SquareMatrix::solveQR(): DTRTRS returned INFO = %d\n", info);
}
if (! m_useReturnErrorCode) {
throw CELapackError("SquareMatrix::solveQR()", "DTRTRS returned INFO = " + int2str(info));
}
}
return info;
}
//=====================================================================================================================

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

@ -86,6 +86,13 @@ namespace Cantera {
*/
int factorQR();
//! 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
*/

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

@ -316,23 +316,24 @@ namespace Cantera {
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, const integer *lwork, int &info) {
doublereal *work, int lwork, int &info) {
char side = left_right[rlside];
char tr = no_yes[trans];
integer f_m = m;
integer f_n = n;
integer f_k = k;
integer f_lwork = lwork;
integer f_lda = lda;
integer f_ldc = ldc;
integer f_info = info;
#ifdef NO_FTN_STRING_LEN_AT_END
_DORMQR_(&side, &tr, &f_m, &f_n, &f_k, a, &f_lda, tau, c, &f_ldc, work, lwork, &f_info);
_DORMQR_(&side, &tr, &f_m, &f_n, &f_k, a, &f_lda, tau, c, &f_ldc, work, &f_lwork, &f_info);
#else
ftnlen trsize = 1;
#ifdef LAPACK_FTN_STRING_LEN_AT_END
_DORMQR_(&side, &tr, &f_m, &f_n, &f_k, a, &f_lda, tau, c, &f_ldc, work, lwork, &f_info, trsize, trsize);
_DORMQR_(&side, &tr, &f_m, &f_n, &f_k, a, &f_lda, tau, c, &f_ldc, work, &f_lwork, &f_info, trsize, trsize);
#else
_DORMQR_(&side, trsize, &tr, trsize, &f_m, &f_n, &f_k, a, &f_lda, tau, c, &f_ldc, work, lwork, &f_info);
_DORMQR_(&side, trsize, &tr, trsize, &f_m, &f_n, &f_k, a, &f_lda, tau, c, &f_ldc, work, &f_lwork, &f_info);
#endif
#endif
info = f_info;