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

Added QR factorization algorithms to Cantera's fortran interface

Browse source 
within the numerics directory.
This commit is contained in:
Harry Moffat 2011-04-06 15:40:04 +00:00
parent 3994bb9438
commit fa30f63056
3 changed files with 153 additions and 15 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -32,7 +32,7 @@ namespace Cantera {
//====================================================================================================================
SquareMatrix::SquareMatrix() :
DenseMatrix(),
m_factored(false)
m_factored(0)
{
}
@ -46,7 +46,7 @@ namespace Cantera {
*/
SquareMatrix::SquareMatrix(int n, doublereal v) :
DenseMatrix(n, n, v),
m_factored(false)
m_factored(0)
{
}
@ -128,7 +128,7 @@ namespace Cantera {
int SquareMatrix::factor() {
integer n = static_cast<int>(nRows());
int info=0;
m_factored = true;
m_factored = 1;
ct_dgetrf(n, n, &(*(begin())), static_cast<int>(nRows()),
DATA_PTR(ipiv()), info);
if (info != 0) {
@ -141,17 +141,47 @@ namespace Cantera {
}
return info;
}
//=====================================================================================================================
/*
* clear the factored flag
*/
void SquareMatrix::clearFactorFlag() {
m_factored = false;
m_factored = 0;
}
/**
//=====================================================================================================================
/*
* set the factored flag
*/
void SquareMatrix::setFactorFlag() {
m_factored = true;
m_factored = 1;
}
//=====================================================================================================================
int SquareMatrix::factorQR() {
if ((int) tau.size() < m_nrows) {
tau.resize(m_nrows, 0.0);
work.resize(9 * m_nrows, 0.0);
}
int info;
int lwork = work.size();
ct_dgeqrf(m_nrows, m_nrows, &(*(begin())), m_nrows, DATA_PTR(tau), DATA_PTR(work), lwork, info);
if (info != 0) {
if (m_printLevel) {
writelogf("SquareMatrix::factorQR(): DGEQRF returned INFO = %d\n", info);
}
if (! m_useReturnErrorCode) {
throw CELapackError("SquareMatrix::factorQR()", "DGEQRF returned INFO = " + int2str(info));
}
}
int lworkOpt = work[0];
if (lworkOpt != lwork) {
work.resize(lworkOpt);
}
return info;
}
//=====================================================================================================================
}

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

@ -77,6 +77,15 @@ namespace Cantera {
* boolean to indicate that the matrix is now A-1.
*/
int factor();
//! Factors the A matrix using the QR algorithm, overwriting A
/*!
* we set m_factored to 2 to indicate the matrix is now QR factored
*
* @return Returns the info variable from lapack
*/
int factorQR();
/**
* clear the factored flag
*/
@ -89,7 +98,10 @@ namespace Cantera {
/*
* the factor flag
*/
bool m_factored;
int m_factored;
vector_fp tau;
vector_fp work;
};
}

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

@ -35,6 +35,9 @@
#define _DSCAL_ dscal
#define _DGEQRF_ dgeqrf
#define _DORMQR_ dormqr
#else
#define _DGEMV_ dgemv_
@ -48,14 +51,21 @@
#define _DSCAL_ dscal_
#define _DGEQRF_ dgreqrf_
#define _DORMQR_ dormqr_
#endif
namespace ctlapack {
typedef enum {Transpose = 1, NoTranspose = 0} transpose_t;
typedef enum {ColMajor = 1, RowMajor = 0} storage_t;
typedef enum {UpperTriangular = 0, LowerTriangular = 1} upperlower_t;
typedef enum {Left = 0, Right = 1} side_t;
}
const char no_yes[2] = {'N', 'T'};
const char upper_lower[2] = {'U', 'L'};
const char left_right[2] = {'L', 'R'};
#ifdef USE_CBLAS
#include <Accelerate.h>
@ -129,7 +139,36 @@ extern "C" {
#endif
int _DSCAL_(integer *n, doublereal *da, doublereal *dx, integer *incx);
void cblas_dscal(const int N, const double alpha, double *X, const int incX);
void cblas_dscal(const int N, const double alpha, double *X, const int incX);
int _DGEQRF_(const integer* m, const integer* n, doublereal* a, const integer* lda,
doublereal* tau, doublereal* work, const integer *lwork, integer *info);
#ifdef LAPACK_FTN_STRING_LEN_AT_END
int _DORMQR_(const char* side, const char* trans, const integer* m, const integer* n,
const integer * k, doublereal* a, const integer* lda,
doublereal* tau, doublereal* c, const integer* ldc,
doublereal* work, const integer *lwork, integer *info, ftnlen sisize, ftnlen trsize);
#else
int _DORMQR_(const char* side, ftnlen sisize, const char* trans, ftnlen trsize, const integer* m,
const integer* n, const integer * k, doublereal* a, const integer* lda,
doublereal* tau,doublereal* c, const integer* ldc,
doublereal* work, const integer *lwork, integer *info);
#endif
#ifdef LAPACK_FTN_STRING_LEN_AT_END
int _DTRTRS_(const char* uplo, const char* trans, const char *diag, const integer* n,
const integer * nrhs, doublereal* a, const integer* lda,
doublereal* b, const integer* ldb, integer *info,
ftnlen upsize, ftnlen trsize, ftnlen disize);
#else
int _DTRTRS_(const char* uplo, ftnlen upsize, const char* trans, ftnlen trsize, const char *diag,
ftnlen disize, const integer* n, const integer * nrhs, doublereal* a, const integer* lda,
doublereal* b, const integer* ldb, integer *info);
#endif
}
//#endif
@ -230,11 +269,9 @@ namespace Cantera {
#else
ftnlen trsize = 1;
#ifdef LAPACK_FTN_STRING_LEN_AT_END
_DGETRS_(&tr, &f_n, &f_nrhs, a, &f_lda, ipiv, b, &f_ldb,
&f_info, trsize);
_DGETRS_(&tr, &f_n, &f_nrhs, a, &f_lda, ipiv, b, &f_ldb, &f_info, trsize);
#else
_DGETRS_(&tr, trsize, &f_n, &f_nrhs, a, &f_lda, ipiv, b, &f_ldb,
&f_info);
_DGETRS_(&tr, trsize, &f_n, &f_nrhs, a, &f_lda, ipiv, b, &f_ldb, &f_info);
#endif
#endif
info = f_info;
@ -260,13 +297,72 @@ namespace Cantera {
rank = f_rank;
}
inline void ct_dscal(int n, doublereal da, doublereal* dx, int incx) {
inline void ct_dscal(int n, doublereal da, doublereal* dx, int incx) {
//integer f_n = n, f_incx = incx;
//_DSCAL_(&f_n, &da, dx, &f_incx);
cblas_dscal(n, da, dx, 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;
integer f_n = n;
integer f_lda = lda;
integer f_lwork = lwork;
integer f_info = info;
_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, const integer *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_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);
#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);
#else
_DORMQR_(&side, trsize, &tr, trsize, &f_m, &f_n, &f_k, a, &f_lda, tau, c, &f_ldc, work, lwork, &f_info);
#endif
#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];
char tr = no_yes[trans];
char dd = 'N';
if (diag) {
dd = diag[0];
}
integer f_n = n;
integer f_nrhs = nrhs;
integer f_lda = lda;
integer f_ldb = ldb;
integer f_info = info;
#ifdef NO_FTN_STRING_LEN_AT_END
_DTRTRS_(&uplo, &tr, &dd, &f_n, &f_nrhs, a, &f_lda, b, &f_ldb, &f_info);
#else
ftnlen trsize = 1;
#ifdef LAPACK_FTN_STRING_LEN_AT_END
_DTRTRS_(&uplo, &tr, &dd, &f_n, &f_nrhs, a, &f_lda, b, &f_ldb, &f_info, trsize, trsize, trsize);
#else
_DTRTRS_(&uplo, trsize, &tr, trsize, &dd, trsize, &f_n, &f_nrhs, a, &f_lda, b, &f_ldb, &f_info);
#endif
#endif
info = f_info;
}
}