Basic upgrade to the DenseMatrix Class.
We've added a capability to change the handling of singular matrix
conditions.
Options exist to either throw an exception or return an error code.
Options exist to print out a warning or to not print out anything.
This commit is contained in:
Harry Moffat
parent
2e50b8a929
commit
44e03fcf79
5 changed files with 208 additions and 56 deletions
|
|
@ -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<int>(A.nRows()),
|
||||
static_cast<int>(A.nColumns()), A.ptrColumn(0), //begin(),
|
||||
static_cast<int>(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<int>(A.nRows()), 1, A.ptrColumn(0), //begin(),
|
||||
static_cast<int>(A.nRows()),
|
||||
&A.ipiv()[0], b,
|
||||
static_cast<int>(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<int>(A.nRows()),
|
||||
static_cast<int>(A.nColumns()), A.ptrColumn(0),
|
||||
static_cast<int>(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<int>(A.nRows()),
|
||||
static_cast<int>(b.nColumns()),
|
||||
A.ptrColumn(0), static_cast<int>(A.nRows()),
|
||||
&A.ipiv()[0], b.ptrColumn(0),
|
||||
static_cast<int>(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<int>(A.nRows()), b,
|
||||
static_cast<int>(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<int>(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<int>(work.size());
|
||||
ct_dgetri(n, A.ptrColumn(0), static_cast<int>(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;
|
||||
}
|
||||
//====================================================================================================================
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
};
|
||||
|
||||
//==================================================================================================================
|
||||
|
|
|
|||
|
|
@ -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<doublereal> y_curr(neq_, 0.0);
|
||||
std::vector<doublereal> 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;
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -17,6 +17,7 @@
|
|||
#include "stringUtils.h"
|
||||
#include "ctlapack.h"
|
||||
#include "SquareMatrix.h"
|
||||
#include "global.h"
|
||||
|
||||
#include <iostream>
|
||||
#include <vector>
|
||||
|
|
@ -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<int>(nRows()),
|
||||
1, &(*(begin())), static_cast<int>(nRows()),
|
||||
DATA_PTR(ipiv()), b, static_cast<int>(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<int>(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
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue