diff --git a/Cantera/src/numerics/SquareMatrix.cpp b/Cantera/src/numerics/SquareMatrix.cpp index c8371d889..9d9d56861 100644 --- a/Cantera/src/numerics/SquareMatrix.cpp +++ b/Cantera/src/numerics/SquareMatrix.cpp @@ -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; + } + //===================================================================================================================== } diff --git a/Cantera/src/numerics/SquareMatrix.h b/Cantera/src/numerics/SquareMatrix.h index 1be99b62e..1445bc60d 100644 --- a/Cantera/src/numerics/SquareMatrix.h +++ b/Cantera/src/numerics/SquareMatrix.h @@ -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 iwork_; + public: //! Use the QR algorithm to factor and invert the matrix int useQR_; diff --git a/Cantera/src/numerics/ctlapack.h b/Cantera/src/numerics/ctlapack.h index 3b160184b..40d8b1754 100644 --- a/Cantera/src/numerics/ctlapack.h +++ b/Cantera/src/numerics/ctlapack.h @@ -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