diff --git a/Cantera/src/numerics/NonlinearSolver.cpp b/Cantera/src/numerics/NonlinearSolver.cpp index 1ace7e722..4d24c605e 100644 --- a/Cantera/src/numerics/NonlinearSolver.cpp +++ b/Cantera/src/numerics/NonlinearSolver.cpp @@ -20,6 +20,7 @@ #include "SquareMatrix.h" #include "NonlinearSolver.h" +#include "ctlapack.h" #include "clockWC.h" #include "vec_functions.h" @@ -79,6 +80,12 @@ namespace Cantera { // Turn off printing of dogleg information bool NonlinearSolver::s_print_DogLeg(false); + + // Turn off solving the system twice and comparing the answer. + /* + * Turn this on if you want to compare the Hessian and Newton solve results. + */ + bool NonlinearSolver::s_doBothSolvesAndCompare(false); //==================================================================================================================== // Default constructor /* @@ -130,6 +137,7 @@ namespace Cantera { m_print_flag(0), m_ScaleSolnNormToResNorm(0.001), jacCopy_(0), + Hessian_(0), deltaX_CP_(0), deltaX_Newton_(0), residNorm2Cauchy_(0.0), @@ -146,7 +154,8 @@ namespace Cantera { JdJd_norm_(0.0), normTrust_Newton_(0.0), normTrust_CP_(0.0), - doDogLeg_(0) + doDogLeg_(0), + doAffineSolve_(0) { neq_ = m_func->nEquations(); @@ -173,12 +182,12 @@ namespace Cantera { m_ewt[i] = atolk_[i]; } -#ifdef DEBUG_DOGLEG + jacCopy_.resize(neq_, neq_, 0.0); deltaX_CP_.resize(neq_, 0.0); Jd_.resize(neq_, 0.0); deltaX_trust_.resize(neq_, 1.0); -#endif + } //==================================================================================================================== @@ -228,6 +237,7 @@ namespace Cantera { m_print_flag(0), m_ScaleSolnNormToResNorm(0.001), jacCopy_(0), + Hessian_(0), deltaX_CP_(0), deltaX_Newton_(0), residNorm2Cauchy_(0.0), @@ -244,7 +254,8 @@ namespace Cantera { JdJd_norm_(0.0), normTrust_Newton_(0.0), normTrust_CP_(0.0), - doDogLeg_(0) + doDogLeg_(0), + doAffineSolve_(0) { *this =operator=(right); } @@ -305,6 +316,7 @@ namespace Cantera { m_ScaleSolnNormToResNorm = right.m_ScaleSolnNormToResNorm; jacCopy_ = right.jacCopy_; + Hessian_ = right.Hessian_; deltaX_CP_ = right.deltaX_CP_; deltaX_Newton_ = right.deltaX_Newton_; RJd_norm_ = right.RJd_norm_; @@ -322,6 +334,7 @@ namespace Cantera { normTrust_Newton_ = right.normTrust_Newton_; normTrust_CP_ = right.normTrust_CP_; doDogLeg_ = right.doDogLeg_; + doAffineSolve_ = right.doAffineSolve_; return *this; } @@ -362,8 +375,17 @@ namespace Cantera { } } //==================================================================================================================== - void NonlinearSolver::setSolverScheme(int doDogLeg) { + void NonlinearSolver::setSolverScheme(int doDogLeg, int doAffineSolve) { doDogLeg_ = doDogLeg; + doAffineSolve_ = doAffineSolve; +#ifdef DEBUG_DOGLEG + +#else + if (doDogLeg_) { + throw CanteraError("NonlinearSolver::setSolverScheme", + "ifdef block not on"); + } +#endif } //==================================================================================================================== std::vector & NonlinearSolver::lowBoundsConstraintVector() { @@ -785,6 +807,305 @@ namespace Cantera { return info; } //==================================================================================================================== + // Compute the newton step, either by direct newton's or by solving a close problem that is represented + // by a Hessian ( + /* + * This is algorith A.6.5.1 in Dennis / Schnabel + * + * Compute the QR decomposition + */ + int NonlinearSolver::doAffineNewtonSolve(const doublereal * const y_curr, const doublereal * const ydot_curr, + doublereal * const delta_y, SquareMatrix& jac) + { + bool newtonGood = true; + int irow; + doublereal *delyNewton = 0; + // We can default to QR here ( or not ) + jac.useQR_ = true; + // multiply the residual by -1 + // Scale the residual if there is row scaling. Note, the matrix has already been scaled + if (m_rowScaling && !m_resid_scaled) { + for (int n = 0; n < neq_; n++) { + delta_y[n] = -m_rowScales[n] * m_resid[n]; + } + m_resid_scaled = true; + } else { + for (int n = 0; n < neq_; n++) { + delta_y[n] = -m_resid[n]; + } + } + + // Factor the matrix + double cond = 1.0E300; + int info = 0; + if (!jac.m_factored) { + if (jac.useQR_) { + info = jac.factorQR(); + } else { + info = jac.factor(); + } + } + /* + * Find the condition number of the matrix + * If we have failed to factor, we will fall back to calculating and factoring a modified Hessian + */ + if (info == 0) { + double rcond = 0.0; + if (jac.useQR_) { + rcond = jac.rcondQR(); + } else { + rcond = jac.rcond(jac.a1norm_); + } + if (rcond > 0.0) { + cond = 1.0 / rcond; + } + } + bool doHessian = false; + if (s_doBothSolvesAndCompare) { + doHessian = true; + } + bool doNewton = false; + if (cond < 1.0E7) { + doNewton = true; + if (m_print_flag >= 3) { + printf("\t\t\tdoAffineNewtonSolve: Condition number = %g during regular solve\n", cond); + } + + /* + * Solve the system -> This also involves inverting the matrix + */ + int info = jac.solve(DATA_PTR(delta_y)); + if (info) { + if (m_print_flag >= 2) { + printf("\t\t\tNonlinearSolver::doAffineSolve QRSolve returned INFO = %d. Switching to Hessian solve\n", info); + } + doHessian = true; + newtonGood = false; + } + /* + * reverse the column scaling if there was any on a successful solve + */ + if (m_colScaling) { + for (irow = 0; irow < neq_; irow++) { + delta_y[irow] = delta_y[irow] * m_colScales[irow]; + } + } + + } else { + doHessian = true; + newtonGood = false; + if (m_print_flag >= 3) { + printf("\t\t\tdoAffineNewtonSolve: Condition number too large, %g. Doing a Hessian solve \n", cond); + } + } + + if (doHessian) { + // Store the old value for later comparison + if (doNewton) { + delyNewton = mdp::mdp_alloc_dbl_1(neq_, MDP_DBL_NOINIT); + for (irow = 0; irow < neq_; irow++) { + delyNewton[irow] = delta_y[irow]; + } + } + // Get memory if not done before + if (Hessian_.nRows() == 0) { + Hessian_.resize(neq_, neq_); + } + + /* + * Calculate the symmetric Hessian + */ + Hessian_.zero(); + if (m_rowScaling) { + for (int i = 0; i < neq_; i++) { + for (int j = i; j < neq_; j++) { + for (int k = 0; k < neq_; k++) { + Hessian_(i,j) += jacCopy_(k,i) * jacCopy_(k,j) * m_rowScales[k] * m_rowScales[k]; + } + Hessian_(j,i) = Hessian_(i,j); + } + } + } else { + for (int i = 0; i < neq_; i++) { + for (int j = i; j < neq_; j++) { + for (int k = 0; k < neq_; k++) { + Hessian_(i,j) += jacCopy_(k,i) * jacCopy_(k,j); + } + Hessian_(j,i) = Hessian_(i,j); + } + } + } + + /* + * Calculate the matrix norm of the Hessian + */ + doublereal hnorm = 0.0; + doublereal hcol = 0.0; + if (m_colScaling) { + for (int i = 0; i < neq_; i++) { + for (int j = i; j < neq_; j++) { + hcol += fabs(Hessian_(j,i)) * m_colScales[j]; + } + for (int j = i+1; j < neq_; j++) { + hcol += fabs(Hessian_(i,j)) * m_colScales[j]; + } + hcol *= m_colScales[i]; + if (hcol > hnorm) { + hnorm = hcol; + } + } + } else { + for (int i = 0; i < neq_; i++) { + for (int j = i; j < neq_; j++) { + hcol += fabs(Hessian_(j,i)); + } + for (int j = i+1; j < neq_; j++) { + hcol += fabs(Hessian_(i,j)); + } + if (hcol > hnorm) { + hnorm = hcol; + } + } + } + /* + * Add junk to the Hessian diagonal + */ + hcol = sqrt(neq_) * 1.0E-7 * hnorm; + if (m_colScaling) { + for (int i = 0; i < neq_; i++) { + Hessian_(i,i) += hcol / (m_colScales[i] * m_colScales[i]); + } + } else { + for (int i = 0; i < neq_; i++) { + Hessian_(i,i) += hcol; + } + } + + /* + * Factor the Hessian + */ + int info; + ct_dpotrf(ctlapack::UpperTriangular, neq_, &(*(Hessian_.begin())), neq_, info); + if (info) { + if (m_print_flag >= 2) { + printf("\t\t\tNonlinearSolver::doAffineSolve DPOTRF returned INFO = %d\n", info); + } + return info; + } + + // doublereal *JTF = delta_y; + doublereal *delyH = mdp::mdp_alloc_dbl_1(neq_, MDP_DBL_NOINIT); + // First recalculate the scaled residual. It got wiped out doing the newton solve + if (m_rowScaling) { + for (int n = 0; n < neq_; n++) { + delyH[n] = -m_rowScales[n] * m_resid[n]; + } + } else { + for (int n = 0; n < neq_; n++) { + delyH[n] = -m_resid[n]; + } + } + + if (m_rowScaling) { + for (int j = 0; j < neq_; j++) { + delta_y[j] = 0.0; + for (int i = 0; i < neq_; i++) { + delta_y[j] += delyH[i] * jacCopy_.value(i,j) * m_rowScales[i]; + } + } + } else { + for (int j = 0; j < neq_; j++) { + delta_y[j] = 0.0; + for (int i = 0; i < neq_; i++) { + delta_y[j] += delyH[i] * jacCopy_.value(i,j); + } + } + } + + + /* + * Solve the factored Hessian System + */ + ct_dpotrs(ctlapack::UpperTriangular, neq_, 1,&(*(Hessian_.begin())), neq_, delta_y, neq_, info); + if (info) { + if (m_print_flag >= 2) { + printf("\t\t\tNonlinearSolver::doAffineSolve DPOTRS returned INFO = %d\n", info); + } + return info; + } + /* + * reverse the column scaling if there was any. + */ + if (m_colScaling) { + for (irow = 0; irow < neq_; irow++) { + delta_y[irow] = delta_y[irow] * m_colScales[irow]; + } + } + + + if (s_print_DogLeg || (doDogLeg_ && m_print_flag > 3)) { + printf("\t\t\t Comparison between Hessian deltaX and newton deltaX\n"); + printf("\t\t\t i Hessian+Junk Newton \n"); + printf("\t\t\t--------------------------------------------------------\n"); + for (int i =0; i < neq_; i++) { + printf("\t\t\t%3d %12.5g %12.5g\n", i, delta_y[i], delyNewton[i]); + } + printf("\t\t\t--------------------------------------------------------\n"); + + } + + if (newtonGood) { + mdp::mdp_copy_dbl_1(DATA_PTR(delta_y), CONSTD_DATA_PTR(delyNewton), neq_); + } + mdp::mdp_safe_free((void **) &delyH); + mdp::mdp_safe_free((void **) &delyNewton); + } + +#ifdef DEBUG_JAC + if (printJacContributions) { + for (int iNum = 0; iNum < numRows; iNum++) { + if (iNum > 0) focusRow++; + doublereal dsum = 0.0; + vector_fp& Jdata = jacBack.data(); + doublereal dRow = Jdata[neq_ * focusRow + focusRow]; + printf("\n Details on delta_Y for row %d \n", focusRow); + printf(" Value before = %15.5e, delta = %15.5e," + "value after = %15.5e\n", y_curr[focusRow], + delta_y[focusRow], y_curr[focusRow] + delta_y[focusRow]); + if (!freshJac) { + printf(" Old Jacobian\n"); + } + printf(" col delta_y aij " + "contrib \n"); + printf("-----------------------------------------------------------------------------------------------\n"); + printf(" Res(%d) %15.5e %15.5e %15.5e (Res = %g)\n", + focusRow, delta_y[focusRow], + dRow, RRow[iNum] / dRow, RRow[iNum]); + dsum += RRow[iNum] / dRow; + for (int ii = 0; ii < neq_; ii++) { + if (ii != focusRow) { + doublereal aij = Jdata[neq_ * ii + focusRow]; + doublereal contrib = aij * delta_y[ii] * (-1.0) / dRow; + dsum += contrib; + if (fabs(contrib) > Pcutoff) { + printf("%6d %15.5e %15.5e %15.5e\n", ii, + delta_y[ii] , aij, contrib); + } + } + } + printf("-----------------------------------------------------------------------------------------------\n"); + printf(" %15.5e %15.5e\n", + delta_y[focusRow], dsum); + } + } + +#endif + + m_numTotalLinearSolves++; + return info; + + } + //==================================================================================================================== // Do a steepest descent calculation /* * This call must be made on the unfactored jacobian! @@ -1794,8 +2115,8 @@ namespace Cantera { m_dampRes = 1.0; int j, m; num_backtracks = 0; - double deltaSolnNorm = solnErrorNorm(DATA_PTR(deltaX_CP_)); - double funcDecreaseSDExp = RJd_norm_ / deltaSolnNorm * lambda_; + //double deltaSolnNorm = solnErrorNorm(DATA_PTR(deltaX_CP_)); + //double funcDecreaseSDExp = RJd_norm_ / deltaSolnNorm * lambda_; double tlen; @@ -2185,7 +2506,12 @@ namespace Cantera { #endif // compute the undamped Newton step - info = doNewtonSolve(time_curr, DATA_PTR(m_y_n), DATA_PTR(ydot_curr), DATA_PTR(deltaX_Newton_), jac, m_print_flag); + if (doAffineSolve_) { + info = doAffineNewtonSolve(DATA_PTR(m_y_n), DATA_PTR(ydot_curr), DATA_PTR(deltaX_Newton_), jac); + } else { + info = doNewtonSolve(time_curr, DATA_PTR(m_y_n), DATA_PTR(ydot_curr), DATA_PTR(deltaX_Newton_), jac, m_print_flag); + } + if (info) { m = -1; goto done; @@ -2201,8 +2527,8 @@ namespace Cantera { if (doDogLeg_) { - double trustD = calcTrustDistance(stp); #ifdef DEBUG_DOGLEG + double trustD = calcTrustDistance(stp); if (s_print_DogLeg || m_print_flag > 3) { if (trustD > trustDelta_) { printf("newton's method trustD, %g, larger than trust region, %g\n", trustD, trustDelta_); @@ -2664,6 +2990,12 @@ namespace Cantera { printf("--------------\n"); } } + /* + * Make a copy of the data. Note, this jacobian copy occurs before any matrix scaling operations. + * It's the raw matrix producted by this routine. + */ + jacCopy_.copyData(J); + return retn; } //==================================================================================================================== diff --git a/Cantera/src/numerics/NonlinearSolver.h b/Cantera/src/numerics/NonlinearSolver.h index e586e8cfe..de8aa9a64 100644 --- a/Cantera/src/numerics/NonlinearSolver.h +++ b/Cantera/src/numerics/NonlinearSolver.h @@ -208,6 +208,40 @@ namespace Cantera { const doublereal * const ydot_curr, doublereal * const delta_y, SquareMatrix& jac, int loglevel); + //! Compute the newton step, either by direct newton's or by solving a close problem that is represented + //! by a Hessian ( + /*! + * This is algorith A.6.5.1 in Dennis / Schnabel + * + * Compute the QR decomposition + * + * Compute the undamped Newton step. The residual function is + * evaluated at the current time, t_n, at the current values of the + * solution vector, m_y_n, and the solution time derivative, m_ydot_n. + * The Jacobian is not recomputed. + * + * A factored jacobian is reused, if available. If a factored jacobian + * is not available, then the jacobian is factored. Before factoring, + * the jacobian is row and column-scaled. Column scaling is not + * recomputed. The row scales are recomputed here, after column + * scaling has been implemented. + * + * @param y_curr Current value of the solution + * @param ydot_curr Current value of the solution derivative. + * @param delta_y return value of the raw change in y + * @param jac Jacobian + * + * Internal input + * --------------- + * internal m_resid Storred residual is used as input + * + * + * @return Returns the result code from lapack. A zero means success. Anything + * else indicates a failure. + */ + int doAffineNewtonSolve(const doublereal * const y_curr, const doublereal * const ydot_curr, + doublereal * const delta_y, SquareMatrix& jac); + //! Calculate the size of the current trust region /*! * We carry out a norm of deltaX_trust_ first. Then, we multiply that value @@ -637,8 +671,11 @@ namespace Cantera { * Default is to always use a damping scheme in the Newton Direction. * When this is nonzero, a model trust region approach is used using a double dog leg * with the steepest descent direction used for small step sizes. + * + * @param doAffineSolve Parameter to turn on or off the solution of the system using a Hessian + * if the matrix has a bad condition number. */ - void setSolverScheme(int doDogLeg = 0); + void setSolverScheme(int doDogLeg, int doAffineSolve); /* * ----------------------------------------------------------------------------------------------------------------- @@ -838,8 +875,14 @@ namespace Cantera { double m_ScaleSolnNormToResNorm; //! Copy of the jacobian that doesn't get overwritten when the inverse is determined + /*! + * The jacobian storred here is the raw matrix, before any row or column scaling is carried out + */ Cantera::SquareMatrix jacCopy_; + //! Hessian + Cantera::SquareMatrix Hessian_; + /********************************************************************************************* * VARIABLES ASSOCIATED WITH STEPS AND ASSOCIATED DOUBLE DOGLEG PARAMETERS *********************************************************************************************/ @@ -891,6 +934,9 @@ namespace Cantera { //! General toggle for turning on dog leg damping. int doDogLeg_; + //! General toggle for turning on Affine solve with Hessian + int doAffineSolve_; + /******************************************************************************************* @@ -910,6 +956,13 @@ namespace Cantera { //! Turn on all printing of dogleg information static bool s_print_DogLeg; + + //! Turn on solving both the Newton and Hessian systems and comparing the results + /*! + * This is off by default + */ + static bool s_doBothSolvesAndCompare; + }; } diff --git a/Cantera/src/numerics/SquareMatrix.cpp b/Cantera/src/numerics/SquareMatrix.cpp index 9d9d56861..1724f8c93 100644 --- a/Cantera/src/numerics/SquareMatrix.cpp +++ b/Cantera/src/numerics/SquareMatrix.cpp @@ -32,7 +32,8 @@ namespace Cantera { //==================================================================================================================== SquareMatrix::SquareMatrix() : DenseMatrix(), - m_factored(0), + m_factored(0), + a1norm_(0.0), useQR_(0) { } @@ -48,8 +49,10 @@ namespace Cantera { */ SquareMatrix::SquareMatrix(int n, doublereal v) : DenseMatrix(n, n, v), - m_factored(0), + m_factored(0), + a1norm_(0.0), useQR_(0) + { } //==================================================================================================================== @@ -60,6 +63,7 @@ namespace Cantera { SquareMatrix::SquareMatrix(const SquareMatrix& y) : DenseMatrix(y), m_factored(y.m_factored), + a1norm_(y.a1norm_), useQR_(y.useQR_) { } @@ -71,7 +75,9 @@ namespace Cantera { SquareMatrix& SquareMatrix::operator=(const SquareMatrix& y) { if (&y == this) return *this; DenseMatrix::operator=(y); - m_factored = y.m_factored; + m_factored = y.m_factored; + a1norm_ = y.a1norm_; + useQR_ = y.useQR_; return *this; } //==================================================================================================================== @@ -140,11 +146,11 @@ namespace Cantera { if (useQR_) { return factorQR(); } + a1norm_ = ct_dlange('1', m_nrows, m_nrows, &(*(begin())), m_nrows, DATA_PTR(work)); integer n = static_cast(nRows()); int info=0; m_factored = 1; - ct_dgetrf(n, n, &(*(begin())), static_cast(nRows()), - DATA_PTR(ipiv()), info); + ct_dgetrf(n, n, &(*(begin())), static_cast(nRows()), DATA_PTR(ipiv()), info); if (info != 0) { if (m_printLevel) { writelogf("SquareMatrix::factor(): DGETRS returned INFO = %d\n", info); @@ -174,7 +180,8 @@ namespace Cantera { if ((int) tau.size() < m_nrows) { tau.resize(m_nrows, 0.0); work.resize(8 * m_nrows, 0.0); - } + } + a1norm_ = ct_dlange('1', m_nrows, m_nrows, &(*(begin())), m_nrows, DATA_PTR(work)); int info; m_factored = 2; int lwork = work.size(); @@ -252,6 +259,36 @@ namespace Cantera { return info; } //===================================================================================================================== + doublereal SquareMatrix::rcond(doublereal anorm) { + + if ((int) iwork_.size() < m_nrows) { + iwork_.resize(m_nrows); + } + if ((int) work.size() <4 * m_nrows) { + work.resize(4 * m_nrows); + } + doublereal rcond = 0.0; + if (m_factored != 1) { + throw CELapackError("SquareMatrix::rcond()", "matrix isn't factored correctly"); + } + + // doublereal anorm = ct_dlange('1', m_nrows, m_nrows, &(*(begin())), m_nrows, DATA_PTR(work)); + + + int rinfo; + rcond = ct_dgecon('1', m_nrows, &(*(begin())), m_nrows, anorm, DATA_PTR(work), + DATA_PTR(iwork_), rinfo); + if (rinfo != 0) { + if (m_printLevel) { + writelogf("SquareMatrix::rcond(): DGECON returned INFO = %d\n", rinfo); + } + if (! m_useReturnErrorCode) { + throw CELapackError("SquareMatrix::rcond()", "DGECON returned INFO = " + int2str(rinfo)); + } + } + return rcond; + } + //===================================================================================================================== doublereal SquareMatrix::rcondQR() { if ((int) iwork_.size() < m_nrows) { diff --git a/Cantera/src/numerics/SquareMatrix.h b/Cantera/src/numerics/SquareMatrix.h index 1445bc60d..b95d3b35b 100644 --- a/Cantera/src/numerics/SquareMatrix.h +++ b/Cantera/src/numerics/SquareMatrix.h @@ -104,6 +104,16 @@ namespace Cantera { */ doublereal rcondQR(); + //! Returns an estimate of the inverse of the condition number for the matrix + /*! + * The matrix must have been previously factored using the LU algorithm + * + * @param a1norm Norm of the matrix + * + * @return returns the inverse of the condition number + */ + doublereal rcond(doublereal a1norm); + //! Solves the linear problem Ax=b using the QR algorithm returning x in the b spot /*! * @param b RHS to be solved. @@ -131,6 +141,9 @@ namespace Cantera { //! Integer work vector for QR algorithms std::vector iwork_; + + //! 1-norm of the matrix. This is determined immediately before every factorization + doublereal a1norm_; public: //! Use the QR algorithm to factor and invert the matrix diff --git a/Cantera/src/numerics/ctlapack.h b/Cantera/src/numerics/ctlapack.h index d79c68391..ea7083c1f 100644 --- a/Cantera/src/numerics/ctlapack.h +++ b/Cantera/src/numerics/ctlapack.h @@ -32,6 +32,8 @@ #define _DGBSV_ dgbsv #define _DGBTRF_ dgbtrf #define _DGBTRS_ dgbtrs +#define _DGECON_ dgecon +#define _DLANGE_ dlange #define _DSCAL_ dscal @@ -52,6 +54,8 @@ #define _DGBSV_ dgbsv_ #define _DGBTRF_ dgbtrf_ #define _DGBTRS_ dgbtrs_ +#define _DGECON_ dgecon_ +#define _DLANGE_ dlange_ #define _DSCAL_ dscal_ @@ -147,7 +151,7 @@ extern "C" { doublereal *b, integer *ldb, integer *info); #endif - int _DSCAL_(integer *n, doublereal *da, doublereal *dx, integer *incx); + int _DSCAL_(integer *n, doublereal *da, doublereal *dx, integer *incx); void cblas_dscal(const int N, const double alpha, double *X, const int incX); @@ -207,6 +211,26 @@ extern "C" { #endif +#ifdef LAPACK_FTN_STRING_LEN_AT_END + int _DGECON_(const char *norm, const integer* n, doublereal* a, const integer* lda, + const doublereal *rnorm, const doublereal *rcond, + doublereal* work, const integer* iwork, integer *info, ftnlen nosize); +#else + int _DGECON_(const char *norm, ftnlen nosize, const integer* n, doublereal* a, const integer* lda, + const doublereal *rnorm, const doublereal *rcond, + doublereal* work, const integer* iwork, integer *info); +#endif + + +#ifdef LAPACK_FTN_STRING_LEN_AT_END + doublereal _DLANGE_(const char *norm, const integer* m, const integer* n, doublereal* a, const integer* lda, + doublereal* work, ftnlen nosize); +#else + doublereal _DLANGE_(const char *norm, ftnlen nosize, const integer* m, const integer* n, doublereal* a, const integer* lda, + doublereal* work); +#endif + + } //#endif @@ -315,7 +339,7 @@ namespace Cantera { } //==================================================================================================================== inline void ct_dgetri(int n, doublereal* a, int lda, integer* ipiv, - doublereal* work, int lwork, int& info) { + doublereal* work, int lwork, int& info) { integer f_n = n, f_lda = lda, f_lwork = lwork, f_info = info; _DGETRI_(&f_n, a, &f_lda, ipiv, work, &f_lwork, &f_info); } @@ -486,9 +510,64 @@ namespace Cantera { info = f_info; return; } + //==================================================================================================================== + //! + /*! + */ + inline doublereal ct_dgecon(const char norm, int n, doublereal* a, int lda, doublereal anorm, + doublereal* work, int *iwork, int &info) { + char cnorm = '1'; + if (norm) { + cnorm = norm; + } + integer f_n = n; + integer f_lda = lda; + integer f_info = info; + doublereal rcond; +#ifdef NO_FTN_STRING_LEN_AT_END + _DGECON_(&cnorm, &f_n a, &f_lda, &anorm, &rcond, work, iwork, &f_info); +#else + ftnlen trsize = 1; +#ifdef LAPACK_FTN_STRING_LEN_AT_END + _DGECON_(&cnorm, &f_n, a, &f_lda, &anorm, &rcond, work, iwork, &f_info, trsize); +#else + _DGECON_(&cnorm, trsize, &f_n, a, &f_lda, &anorm, &rcond, work, iwork, &f_info); +#endif +#endif + info = f_info; + return rcond; + } + //==================================================================================================================== + //! + /*! + */ + inline doublereal ct_dlange(const char norm, int m, int n, doublereal* a, int lda, + doublereal* work) { + char cnorm = '1'; + if (norm) { + cnorm = norm; + } + integer f_m = m; + integer f_n = n; + integer f_lda = lda; + doublereal anorm; + +#ifdef NO_FTN_STRING_LEN_AT_END + anorm = _DLANGE_(&cnorm, &f_m, &f_n a, &f_lda, work); +#else + ftnlen trsize = 1; +#ifdef LAPACK_FTN_STRING_LEN_AT_END + anorm = _DLANGE_(&cnorm, &f_m, &f_n, a, &f_lda, work, trsize); +#else + anorm = _DLANGE_(&cnorm, trsize, &f_m, &f_n, a, &f_lda, work); +#endif +#endif + return anorm; + } + //==================================================================================================================== } #endif