diff --git a/Cantera/src/numerics/DenseMatrix.cpp b/Cantera/src/numerics/DenseMatrix.cpp index ae5376a77..89f6d7519 100644 --- a/Cantera/src/numerics/DenseMatrix.cpp +++ b/Cantera/src/numerics/DenseMatrix.cpp @@ -14,11 +14,16 @@ #include "utilities.h" #include "DenseMatrix.h" #include "stringUtils.h" +#include "global.h" namespace Cantera { //==================================================================================================================== // Default Constructor - DenseMatrix::DenseMatrix() + DenseMatrix::DenseMatrix() : + Array2D(0,0,0.0), + m_ipiv(0), + m_useReturnErrorCode(0), + m_printLevel(0) { } //==================================================================================================================== @@ -27,17 +32,23 @@ namespace Cantera { * all elements to \c v. */ DenseMatrix::DenseMatrix(int n, int m, doublereal v) : - Array2D(n, m, v) + Array2D(n, m, v), + m_ipiv(0), + m_useReturnErrorCode(0), + m_printLevel(0) { m_ipiv.resize(max(n, m)); } //==================================================================================================================== - // copy constructor + // Copy constructor /* * @param y Object to be copied */ DenseMatrix::DenseMatrix(const DenseMatrix& y) : - Array2D(y) + Array2D(y), + m_ipiv(0), + m_useReturnErrorCode(0), + m_printLevel(0) { m_ipiv = y.ipiv(); } @@ -47,6 +58,8 @@ namespace Cantera { if (&y == this) return *this; Array2D::operator=(y); m_ipiv = y.ipiv(); + m_useReturnErrorCode = y.m_useReturnErrorCode; + m_printLevel = y.m_printLevel; return *this; } //==================================================================================================================== @@ -86,56 +99,113 @@ namespace Cantera { } //==================================================================================================================== int solve(DenseMatrix& A, double* b) { + int info = 0; if (A.nColumns() != A.nRows()) { - throw CanteraError("DenseMatrix::solve", "Can only solve a square matrix"); + if (A.m_printLevel) { + writelogf("solve(DenseMatrix& A, double* b): Can only solve a square matrix\n"); + } + if (! A.m_useReturnErrorCode) { + throw CELapackError("solve(DenseMatrix& A, double* b)", "Can only solve a square matrix"); + } + return -1; } - int info=0; ct_dgetrf(static_cast(A.nRows()), static_cast(A.nColumns()), A.ptrColumn(0), //begin(), static_cast(A.nRows()), &A.ipiv()[0], info); if (info != 0) { if (info > 0) { - throw CanteraError("DenseMatrix::solve", + if (A.m_printLevel) { + writelogf("solve(DenseMatrix& A, double* b): DGETRF returned INFO = %d U(i,i) is exactly zero. The factorization has" + " been completed, but the factor U is exactly singular, and division by zero will occur if " + "it is used to solve a system of equations.\n", info); + } + if (! A.m_useReturnErrorCode) { + throw CELapackError("solve(DenseMatrix& A, double* b)", "DGETRF returned INFO = "+int2str(info) + ". U(i,i) is exactly zero. The factorization has" " been completed, but the factor U is exactly singular, and division by zero will occur if " "it is used to solve a system of equations."); - + } } else { - throw CanteraError("DenseMatrix::solve", + if (A.m_printLevel) { + writelogf("solve(DenseMatrix& A, double* b): DGETRF returned INFO = %d. The argument i has an illegal value\n", info); + } + if (! A.m_useReturnErrorCode) { + throw CELapackError("solve(DenseMatrix& A, double* b)", "DGETRF returned INFO = "+int2str(info) + ". The argument i has an illegal value"); + } } + return info; } ct_dgetrs(ctlapack::NoTranspose, static_cast(A.nRows()), 1, A.ptrColumn(0), //begin(), static_cast(A.nRows()), &A.ipiv()[0], b, static_cast(A.nColumns()), info); - if (info != 0) - throw CanteraError("DenseMatrix::solve", - "DGETRS returned INFO = "+int2str(info)); - return 0; + if (info != 0) { + if (A.m_printLevel) { + writelogf("solve(DenseMatrix& A, double* b): DGETRS returned INFO = %d\n", info); + } + if (! A.m_useReturnErrorCode) { + throw CELapackError("solve(DenseMatrix& A, double* b)", "DGETRS returned INFO = "+int2str(info)); + } + } + return info; } //==================================================================================================================== int solve(DenseMatrix& A, DenseMatrix& b) { + int info = 0; if (A.nColumns() != A.nRows()) { - throw CanteraError("DenseMatrix::solve", "Can only solve a square matrix"); + if (A.m_printLevel) { + writelogf("solve(DenseMatrix& A, DenseMatrix& b): Can only solve a square matrix\n"); + } + if (! A.m_useReturnErrorCode) { + throw CELapackError("solve(DenseMatrix& A, DenseMatrix& b)", "Can only solve a square matrix"); + } + return -1; } - int info=0; ct_dgetrf(static_cast(A.nRows()), static_cast(A.nColumns()), A.ptrColumn(0), static_cast(A.nRows()), &A.ipiv()[0], info); - if (info != 0) - throw CanteraError("DenseMatrix::solve", - "DGETRF returned INFO = "+int2str(info)); + if (info != 0) { + if (info > 0) { + if (A.m_printLevel) { + writelogf("solve(DenseMatrix& A, DenseMatrix& b): DGETRF returned INFO = %d U(i,i) is exactly zero. The factorization has" + " been completed, but the factor U is exactly singular, and division by zero will occur if " + "it is used to solve a system of equations.\n", info); + } + if (! A.m_useReturnErrorCode) { + throw CELapackError("solve(DenseMatrix& A, DenseMatrix& b)", + "DGETRF returned INFO = "+int2str(info) + ". U(i,i) is exactly zero. The factorization has" + " been completed, but the factor U is exactly singular, and division by zero will occur if " + "it is used to solve a system of equations."); + } + } else { + if (A.m_printLevel) { + writelogf("solve(DenseMatrix& A, DenseMatrix& b): DGETRF returned INFO = %d. The argument i has an illegal value\n", info); + } + if (! A.m_useReturnErrorCode) { + throw CELapackError("solve(DenseMatrix& A, DenseMatrix& b)", + "DGETRF returned INFO = "+int2str(info) + ". The argument i has an illegal value"); + } + } + return info; + } + ct_dgetrs(ctlapack::NoTranspose, static_cast(A.nRows()), static_cast(b.nColumns()), A.ptrColumn(0), static_cast(A.nRows()), &A.ipiv()[0], b.ptrColumn(0), static_cast(b.nRows()), info); - if (info != 0) - throw CanteraError("DenseMatrix::solve", - "DGETRS returned INFO = "+int2str(info)); - return 0; + if (info != 0) { + if (A.m_printLevel) { + writelogf("solve(DenseMatrix& A, DenseMatrix& b): DGETRS returned INFO = %d\n", info); + } + if (! A.m_useReturnErrorCode) { + throw CELapackError("solve(DenseMatrix& A, DenseMatrix& b)", "DGETRS returned INFO = "+int2str(info)); + } + } + + return info; } //==================================================================================================================== @@ -155,10 +225,15 @@ namespace Cantera { static_cast(A.nRows()), b, static_cast(A.nColumns()), &s[0], //.begin(), rcond, rank, &work[0], work.size(), info); - if (info != 0) - throw CanteraError("DenseMatrix::leaseSquares", - "DGELSS returned INFO = "+int2str(info)); - return 0; + if (info != 0) { + if (A.m_printLevel) { + writelogf("leastSquares(): DGELSS returned INFO = %d\n", info); + } + if (! A.m_useReturnErrorCode) { + throw CELapackError("leastSquares()", "DGELSS returned INFO = " + int2str(info)); + } + } + return info; } #endif //==================================================================================================================== @@ -179,19 +254,29 @@ namespace Cantera { int info=0; ct_dgetrf(n, n, A.ptrColumn(0), static_cast(A.nRows()), &A.ipiv()[0], info); - if (info != 0) - throw CanteraError("invert", - "DGETRF returned INFO="+int2str(info)); + if (info != 0) { + if (A.m_printLevel) { + writelogf("invert(DenseMatrix& A, int nn): DGETRS returned INFO = %d\n", info); + } + if (! A.m_useReturnErrorCode) { + throw CELapackError("invert(DenseMatrix& A, int nn)", "DGETRS returned INFO = "+int2str(info)); + } + return info; + } vector_fp work(n); integer lwork = static_cast(work.size()); ct_dgetri(n, A.ptrColumn(0), static_cast(A.nRows()), - &A.ipiv()[0], - &work[0], lwork, info); - if (info != 0) - throw CanteraError("invert", - "DGETRI returned INFO="+int2str(info)); - return 0; + &A.ipiv()[0], &work[0], lwork, info); + if (info != 0) { + if (A.m_printLevel) { + writelogf("invert(DenseMatrix& A, int nn): DGETRS returned INFO = %d\n", info); + } + if (! A.m_useReturnErrorCode) { + throw CELapackError("invert(DenseMatrix& A, int nn)", "DGETRI returned INFO="+int2str(info)); + } + } + return info; } //==================================================================================================================== } diff --git a/Cantera/src/numerics/DenseMatrix.h b/Cantera/src/numerics/DenseMatrix.h index 6d43eaef4..c9d38af47 100644 --- a/Cantera/src/numerics/DenseMatrix.h +++ b/Cantera/src/numerics/DenseMatrix.h @@ -31,6 +31,22 @@ namespace Cantera { * */ + + //! Exception thrown when an LAPACK error is encountered associated with inverting or solving a matrix + class CELapackError : public CanteraError { + public: + + //! Constructor passes through to main Cantera error handler + /*! + * @param routine Name of calling routine + * @param msg Informative message + */ + CELapackError(std::string routine, std::string msg) : + CanteraError(routine + " LAPACK ERROR", msg) + { + } + + }; //! A class for full (non-sparse) matrices with Fortran-compatible //! data storage, which adds matrix operations to class Array2D. @@ -116,6 +132,24 @@ namespace Cantera { //! Vector of pivots. Length is equal to the max of m and n. vector_int m_ipiv; + + public: + //! Error Handling Flag + /*! + * The default is to set this to 0. In this case, if a factorization is requested and can't be achieved, + * a CESingularMatrix exception is triggered. No return code is used, because an exception is thrown. + * If this is set to 1, then an exception is not thrown. Routines return with an error code, that is up + * to the calling routine to handle correctly + */ + int m_useReturnErrorCode; + + //! Print Level + /*! + * Printing is done to the log file using the routine writelogf(). + * + * Level of printing that is carried out. Only error conditions are printed out, if this value is nonzero. + */ + int m_printLevel; }; //================================================================================================================== diff --git a/Cantera/src/numerics/NonlinearSolver.cpp b/Cantera/src/numerics/NonlinearSolver.cpp index 9cc67f1ee..c0e65774c 100644 --- a/Cantera/src/numerics/NonlinearSolver.cpp +++ b/Cantera/src/numerics/NonlinearSolver.cpp @@ -601,8 +601,8 @@ namespace Cantera { * recomputed. The row scales are recomputed here, after column * scaling has been implemented. */ - void NonlinearSolver::doNewtonSolve(const doublereal time_curr, const doublereal * const y_curr, - const doublereal * const ydot_curr, double* const delta_y, SquareMatrix& jac, int loglevel) + int NonlinearSolver::doNewtonSolve(const doublereal time_curr, const doublereal * const y_curr, + const doublereal * const ydot_curr, double* const delta_y, SquareMatrix& jac, int loglevel) { int irow; @@ -625,7 +625,7 @@ namespace Cantera { * Solve the system -> This also involves inverting the * matrix */ - (void) jac.solve(DATA_PTR(delta_y)); + int info = jac.solve(DATA_PTR(delta_y)); /* @@ -681,6 +681,7 @@ namespace Cantera { #endif m_numTotalLinearSolves++; + return info; } //==================================================================================================================== void NonlinearSolver::setDefaultDeltaBoundsMagnitudes() @@ -922,7 +923,7 @@ namespace Cantera { double& s1, SquareMatrix& jac, int& loglevel, bool writetitle, int& num_backtracks) { - + int info = 0; int retnTrial = -2; // Compute the weighted norm of the undamped step size step0 doublereal s0 = solnErrorNorm(step0); @@ -1020,9 +1021,15 @@ namespace Cantera { // Compute the next undamped step, step1[], that would result if y1[] were accepted. // We now have two steps that we have calculated step0[] and step1[] if (solnType_ != NSOLN_TYPE_STEADY_STATE) { - doNewtonSolve(time_curr, y1, ydot1, step1, jac, loglevel); + info = doNewtonSolve(time_curr, y1, ydot1, step1, jac, loglevel); } else { - doNewtonSolve(time_curr, y1, ydot0, step1, jac, loglevel); + info = doNewtonSolve(time_curr, y1, ydot0, step1, jac, loglevel); + } + if (info) { + if (loglevel > 0) { + printf("\t\t\tdampStep: current trial step and damping led to LAPACK ERROR %d. Bailing\n", info); + } + return -1; } // compute the weighted norm of step1 @@ -1133,11 +1140,13 @@ namespace Cantera { clockWC wc; int convRes = 0; solnType_ = SolnType; + int info = 0; bool m_residCurrent = false; int m = 0; bool forceNewJac = false; doublereal s1=1.e30; + // std::vector y_curr(neq_, 0.0); std::vector ydot_curr(neq_, 0.0); @@ -1161,6 +1170,11 @@ namespace Cantera { num_backtracks = 0; int i_backtracks; m_print_flag = loglevelInput; + if (m_print_flag > 1) { + jac.m_printLevel = 1; + } else { + jac.m_printLevel = 0; + } if (m_print_flag == 2 || m_print_flag == 3) { @@ -1249,8 +1263,11 @@ namespace Cantera { // compute the undamped Newton step - doNewtonSolve(time_curr, DATA_PTR(m_y_n), DATA_PTR(ydot_curr), DATA_PTR(stp), jac, m_print_flag); - + info = doNewtonSolve(time_curr, DATA_PTR(m_y_n), DATA_PTR(ydot_curr), DATA_PTR(stp), jac, m_print_flag); + if (info) { + m = -1; + goto done; + } if (m_print_flag > 3) { m_normSolnFRaw = solnErrorNorm(DATA_PTR(stp), "Initial Undamped Step of the iteration", 10); @@ -1262,7 +1279,7 @@ namespace Cantera { /* * Filter out bad directions */ - doublereal normFilter = filterNewStep(time_curr, DATA_PTR(m_y_n), DATA_PTR(stp)); + filterNewStep(time_curr, DATA_PTR(m_y_n), DATA_PTR(stp)); // Damp the Newton step @@ -1407,6 +1424,8 @@ namespace Cantera { printf("\t\tNonlinear problem solved successfully in %d its, time elapsed = %g sec\n", num_newt_its, time_elapsed); } + } else { + printf("\t\tNonlinear problem failed to solve after %d its\n", num_newt_its); } } return m; diff --git a/Cantera/src/numerics/NonlinearSolver.h b/Cantera/src/numerics/NonlinearSolver.h index 65136e4c7..21062cc8e 100644 --- a/Cantera/src/numerics/NonlinearSolver.h +++ b/Cantera/src/numerics/NonlinearSolver.h @@ -166,10 +166,12 @@ namespace Cantera { * @param y_current Current value of the solution * @param ydot_current Current value of the solution derivative. * + * @return returns the result code from lapack. A zero means success. Anything + * else indicates a failure. */ - void doNewtonSolve(const doublereal time_curr, const doublereal * const y_curr, - const doublereal * const ydot_curr, doublereal * const delta_y, - SquareMatrix& jac, int loglevel); + int doNewtonSolve(const doublereal time_curr, const doublereal * const y_curr, + const doublereal * const ydot_curr, doublereal * const delta_y, + SquareMatrix& jac, int loglevel); //! Set default deulta bounds amounts diff --git a/Cantera/src/numerics/SquareMatrix.cpp b/Cantera/src/numerics/SquareMatrix.cpp index 8b34c083f..41c7857cc 100644 --- a/Cantera/src/numerics/SquareMatrix.cpp +++ b/Cantera/src/numerics/SquareMatrix.cpp @@ -17,6 +17,7 @@ #include "stringUtils.h" #include "ctlapack.h" #include "SquareMatrix.h" +#include "global.h" #include #include @@ -56,7 +57,10 @@ namespace Cantera { * Check to see whether the matrix has been factored. */ if (!m_factored) { - factor(); + int retn = factor(); + if (retn) { + return retn; + } } /* * Solve the factored system @@ -64,10 +68,15 @@ namespace Cantera { ct_dgetrs(ctlapack::NoTranspose, static_cast(nRows()), 1, &(*(begin())), static_cast(nRows()), DATA_PTR(ipiv()), b, static_cast(nColumns()), info); - if (info != 0) - throw CanteraError("SquareMatrix::solve", - "DGETRS returned INFO = "+int2str(info)); - return 0; + if (info != 0) { + if (m_printLevel) { + writelogf("SquareMatrix::solve(): DGETRS returned INFO = %d\n", info); + } + if (! m_useReturnErrorCode) { + throw CELapackError("SquareMatrix::solve()", "DGETRS returned INFO = " + int2str(info)); + } + } + return info; } /** @@ -96,11 +105,14 @@ namespace Cantera { ct_dgetrf(n, n, &(*(begin())), static_cast(nRows()), DATA_PTR(ipiv()), info); if (info != 0) { - cout << "Singular matrix, info = " << info << endl; - throw CanteraError("invert", - "DGETRF returned INFO="+int2str(info)); + if (m_printLevel) { + writelogf("SquareMatrix::factor(): DGETRS returned INFO = %d\n", info); + } + if (! m_useReturnErrorCode) { + throw CELapackError("SquareMatrix::factor()", "DGETRS returned INFO = "+int2str(info)); + } } - return 0; + return info; } /* * clear the factored flag