Worked on getting the Affine solve capability working.
This is now working, so that we can switch to a Hessian with junk on the diagonal when there is a condition number problem.
This commit is contained in:
Harry Moffat
parent
e110f91a12
commit
6eccfa16ed
5 changed files with 532 additions and 18 deletions
|
|
@ -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<double> & 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;
|
||||
}
|
||||
//====================================================================================================================
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
||||
};
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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<int>(nRows());
|
||||
int info=0;
|
||||
m_factored = 1;
|
||||
ct_dgetrf(n, n, &(*(begin())), static_cast<int>(nRows()),
|
||||
DATA_PTR(ipiv()), info);
|
||||
ct_dgetrf(n, n, &(*(begin())), static_cast<int>(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) {
|
||||
|
|
|
|||
|
|
@ -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<int> 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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue