research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Made the BandMatrix class and SquareMatrix inherit from a common

Browse source 
parent so that they can be used interchangeably in other software.

Removed PsuedBinaryVPSSTP in favor of MolarityIonicVPSSTP
This commit is contained in:
Harry Moffat 2011-10-19 20:45:06 +00:00
parent 0d74e7ddd0
commit b9c9f81761
40 changed files with 5850 additions and 606 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

4
Cantera/fortran/f77demos/Makefile.in
View file

@ -18,6 +18,10 @@ FORT = @F77@
# Fortran compile flags
FORT_FLAGS = @FFLAGS@
FCLIBS= @FCLIBS@
FLIBS = @FLIBS@
# Fortran libraries
FORT_LIBS = @LCXX_FLIBS@ @LCXX_END_LIBS@ @FLIBS@

16
Cantera/src/numerics/BEulerInt.cpp
View file

@ -421,7 +421,7 @@ namespace Cantera {
* Here a small weighting indicates that the change in solution is
* very sensitive to that equation.
*/
void BEulerInt::computeResidWts(SquareMatrix &jac)
void BEulerInt::computeResidWts(GeneralMatrix &jac)
{
int i, j;
double *data = &(*(jac.begin()));
@ -637,7 +637,7 @@ namespace Cantera {
* not have to be computed again.
*
*/
void BEulerInt::beuler_jac(SquareMatrix &J, double * const f,
void BEulerInt::beuler_jac(GeneralMatrix &J, double * const f,
double time_curr, double CJ,
double * const y,
double * const ydot,
@ -1664,7 +1664,7 @@ namespace Cantera {
*/
void BEulerInt::doNewtonSolve(double time_curr, double *y_curr,
double *ydot_curr, double* delta_y,
SquareMatrix& jac, int loglevel)
GeneralMatrix& jac, int loglevel)
{
int irow, jcol;
@ -1682,7 +1682,7 @@ namespace Cantera {
* by the nominal important change in the solution vector
*/
if (m_colScaling) {
if (!jac.m_factored) {
if (!jac.factored()) {
/*
* Go get new scales
*/
@ -1702,7 +1702,7 @@ namespace Cantera {
}
if (m_matrixConditioning) {
if (jac.m_factored) {
if (jac.factored()) {
m_func->matrixConditioning(0, sz, delta_y);
} else {
double *jptr = &(*(jac.begin()));
@ -1716,7 +1716,7 @@ namespace Cantera {
* nonnegative.
*/
if (m_rowScaling) {
if (! jac.m_factored) {
if (! jac.factored()) {
/*
* Ok, this is ugly. jac.begin() returns an vector<double> iterator
* to the first data location.
@ -1940,7 +1940,7 @@ namespace Cantera {
int BEulerInt::dampStep(double time_curr, const double * y0,
const double *ydot0, const double* step0,
double* y1, double* ydot1, double* step1,
double& s1, SquareMatrix& jac,
double& s1, GeneralMatrix & jac,
int& loglevel, bool writetitle,
int& num_backtracks) {
@ -2107,7 +2107,7 @@ namespace Cantera {
int BEulerInt::solve_nonlinear_problem(double * const y_comm,
double * const ydot_comm, double CJ,
double time_curr,
SquareMatrix& jac,
GeneralMatrix & jac,
int &num_newt_its,
int &num_linear_solves,
int &num_backtracks,

14
Cantera/src/numerics/BEulerInt.h
View file

@ -24,7 +24,7 @@
#include "Integrator.h"
#include "ResidJacEval.h"
#include "SquareMatrix.h"
#include "GeneralMatrix.h"
#include "NonlinearSolver.h"
#include "mdp_allo.h"
@ -119,7 +119,7 @@ namespace Cantera {
bool printLargest = false);
virtual void setInitialTimeStep(double delta_t);
void beuler_jac(SquareMatrix &, double * const,
void beuler_jac(GeneralMatrix &, double * const,
double, double, double * const, double * const, int);
@ -175,7 +175,7 @@ namespace Cantera {
int solve_nonlinear_problem(double * const y_comm,
double * const ydot_comm, double CJ,
double time_curr,
SquareMatrix& jac,
GeneralMatrix& jac,
int &num_newt_its,
int &num_linear_solves,
int &num_backtracks,
@ -186,7 +186,7 @@ namespace Cantera {
* evaluated at x, but the Jacobian is not recomputed.
*/
void doNewtonSolve(double, double *, double*, double *,
SquareMatrix&, int);
GeneralMatrix&, int);
//! Bound the Newton step while relaxing the solution
@ -228,12 +228,12 @@ namespace Cantera {
*/
int dampStep(double, const double*, const double*,
const double *, double*, double*,
double*, double&, SquareMatrix&, int&, bool, int&);
double*, double&, GeneralMatrix&, int&, bool, int&);
/*
* Compute Residual Weights
*/
void computeResidWts(SquareMatrix &jac);
void computeResidWts(GeneralMatrix &jac);
/*
* Filter a new step
@ -421,7 +421,7 @@ namespace Cantera {
* Pointer to the jacobian representing the
* time dependent problem
*/
SquareMatrix *tdjac_ptr;
GeneralMatrix *tdjac_ptr;
/**
* Determines the level of printing for each time
* step.

247
Cantera/src/numerics/BandMatrix.cpp
View file

@ -17,6 +17,7 @@
#include "ctexceptions.h"
#include "stringUtils.h"
#include "global.h"
#include <cstring>
using namespace std;
@ -24,6 +25,7 @@ namespace Cantera {
//====================================================================================================================
BandMatrix::BandMatrix() :
GeneralMatrix(1),
m_factored(false),
m_n(0),
m_kl(0),
@ -35,6 +37,7 @@ namespace Cantera {
}
//====================================================================================================================
BandMatrix::BandMatrix(int n, int kl, int ku, doublereal v) :
GeneralMatrix(1),
m_factored(false),
m_n(n),
m_kl(kl),
@ -46,9 +49,15 @@ namespace Cantera {
fill(data.begin(), data.end(), v);
fill(ludata.begin(), ludata.end(), 0.0);
m_ipiv.resize(m_n);
m_colPtrs.resize(n);
int ldab = (2*kl + ku + 1);
for (int j = 0; j < n; j++) {
m_colPtrs[j] = &(data[ldab * j]);
}
}
//====================================================================================================================
BandMatrix::BandMatrix(const BandMatrix& y) :
GeneralMatrix(1),
m_factored(false),
m_n(0),
m_kl(0),
@ -62,6 +71,11 @@ namespace Cantera {
ludata = y.ludata;
m_factored = y.m_factored;
m_ipiv = y.m_ipiv;
m_colPtrs.resize(m_n);
int ldab = (2 *m_kl + m_ku + 1);
for (int j = 0; j < m_n; j++) {
m_colPtrs[j] = &(data[ldab * j]);
}
}
//====================================================================================================================
BandMatrix::~BandMatrix() {
@ -70,6 +84,7 @@ namespace Cantera {
//====================================================================================================================
BandMatrix& BandMatrix::operator=(const BandMatrix & y) {
if (&y == this) return *this;
GeneralMatrix::operator=(y);
m_n = y.m_n;
m_kl = y.m_kl;
m_ku = y.m_ku;
@ -77,6 +92,11 @@ namespace Cantera {
data = y.data;
ludata = y.ludata;
m_factored = y.m_factored;
m_colPtrs.resize(m_n);
int ldab = (2 * m_kl + m_ku + 1);
for (int j = 0; j < m_n; j++) {
m_colPtrs[j] = &(data[ldab * j]);
}
return *this;
}
//====================================================================================================================
@ -88,6 +108,11 @@ namespace Cantera {
ludata.resize(n*(2*kl + ku + 1));
m_ipiv.resize(m_n);
fill(data.begin(), data.end(), v);
m_colPtrs.resize(m_n);
int ldab = (2 * m_kl + m_ku + 1);
for (int j = 0; j < n; j++) {
m_colPtrs[j] = &(data[ldab * j]);
}
m_factored = false;
}
//====================================================================================================================
@ -96,6 +121,11 @@ namespace Cantera {
m_factored = false;
}
//====================================================================================================================
void BandMatrix::zero() {
std::fill(data.begin(), data.end(), 0.0);
m_factored = false;
}
//====================================================================================================================
doublereal& BandMatrix::operator()(int i, int j) {
return value(i,j);
}
@ -127,7 +157,16 @@ namespace Cantera {
}
//====================================================================================================================
// Number of rows
int BandMatrix::nRows() const {
size_t BandMatrix::nRows() const {
return m_n;
}
//====================================================================================================================
// Number of rows
size_t BandMatrix::nRowsAndStruct(int * const iStruct) const {
if (iStruct) {
iStruct[0] = m_kl;
iStruct[1] = m_ku;
}
return m_n;
}
//====================================================================================================================
@ -208,12 +247,12 @@ namespace Cantera {
return info;
}
//====================================================================================================================
int BandMatrix::solve(int n, const doublereal * const b, doublereal * const x) {
copy(b, b+n, x);
return solve(n, x);
int BandMatrix::solve(const doublereal * const b, doublereal * const x) {
copy(b, b + m_n, x);
return solve(x);
}
//====================================================================================================================
int BandMatrix::solve(int n, doublereal* b) {
int BandMatrix::solve(doublereal* b) {
int info = 0;
if (!m_factored) info = factor();
if (info == 0)
@ -259,6 +298,202 @@ namespace Cantera {
}
return s;
}
//====================================================================================================================
//====================================================================================================================
void BandMatrix::err(std::string msg) const {
throw CanteraError("BandMatrix() unimplemented function", msg);
}
//====================================================================================================================
// Factors the A matrix using the QR algorithm, overwriting A
/*
* we set m_factored to 2 to indicate the matrix is now QR factored
*
* @return Returns the info variable from lapack
*/
int BandMatrix::factorQR() {
factor();
return 0;
}
//====================================================================================================================
// Factors the A matrix using the QR algorithm, overwriting A
// 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 BandMatrix::rcondQR() {
double a1norm = oneNorm();
return rcond(a1norm);
}
//====================================================================================================================
// 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 BandMatrix::rcond(doublereal a1norm) {
int printLevel = 0;
int useReturnErrorCode = 0;
if ((int) iwork_.size() < m_n) {
iwork_.resize(m_n);
}
if ((int) work_.size() < 3 * m_n) {
work_.resize(3 * m_n);
}
doublereal rcond = 0.0;
if (m_factored != 1) {
throw CanteraError("BandMatrix::rcond()", "matrix isn't factored correctly");
}
// doublereal anorm = oneNorm();
int ldab = (2 *m_kl + m_ku + 1);
int rinfo;
rcond = ct_dgbcon('1', m_n, m_kl, m_ku, DATA_PTR(ludata), ldab, DATA_PTR(m_ipiv), a1norm, DATA_PTR(work_),
DATA_PTR(iwork_), rinfo);
if (rinfo != 0) {
if (printLevel) {
writelogf("BandMatrix::rcond(): DGBCON returned INFO = %d\n", rinfo);
}
if (! useReturnErrorCode) {
throw CanteraError("BandMatrix::rcond()", "DGBCON returned INFO = " + int2str(rinfo));
}
}
return rcond;
}
//====================================================================================================================
// Change the way the matrix is factored
/*
* @param fAlgorithm integer
* 0 LU factorization
* 1 QR factorization
*/
void BandMatrix::useFactorAlgorithm(int fAlgorithm) {
// QR algorithm isn't implemented for banded matrix.
}
//====================================================================================================================
int BandMatrix::factorAlgorithm() const {
return 0;
}
//====================================================================================================================
// Returns the one norm of the matrix
doublereal BandMatrix::oneNorm() const {
doublereal value = 0.0;
for (int j = 0; j < m_n; j++) {
doublereal sum = 0.0;
doublereal *colP = m_colPtrs[j];
for (int i = j - m_ku; i <= j + m_kl; i++) {
sum += fabs(colP[m_kl + m_ku + i - j]);
}
if (sum > value) {
value = sum;
}
}
return value;
}
//====================================================================================================================
int BandMatrix::checkRows(doublereal &valueSmall) const {
valueSmall = 1.0E300;
int iSmall = -1;
double vv;
for (int i = 0; i < m_n; i++) {
double valueS = 0.0;
for (int j = i - m_kl; j <= i + m_ku; j++) {
if (j >= 0 && (j < m_n)) {
vv = fabs(value(i,j));
if (vv > valueS) {
valueS = vv;
}
}
}
if (valueS < valueSmall) {
iSmall = i;
valueSmall = valueS;
if (valueSmall == 0.0) {
return iSmall;
}
}
}
return iSmall;
}
//====================================================================================================================
int BandMatrix::checkColumns(doublereal &valueSmall) const {
valueSmall = 1.0E300;
int jSmall = -1;
double vv;
for (int j = 0; j < m_n; j++) {
double valueS = 0.0;
for (int i = j - m_ku; i <= j + m_kl; i++) {
if (i >= 0 && (i < m_n)) {
vv = fabs(value(i,j));
if (vv > valueS) {
valueS = vv;
}
}
}
if (valueS < valueSmall) {
jSmall = j;
valueSmall = valueS;
if (valueSmall == 0.0) {
return jSmall;
}
}
}
return jSmall;
}
//====================================================================================================================
GeneralMatrix * BandMatrix::duplMyselfAsGeneralMatrix() const {
BandMatrix *dd = new BandMatrix(*this);
return static_cast<GeneralMatrix *>(dd);
}
//====================================================================================================================
bool BandMatrix::factored() const {
return m_factored;
}
//====================================================================================================================
// Return a pointer to the top of column j, columns are assumed to be contiguous in memory
/*
* @param j Value of the column
*
* @return Returns a pointer to the top of the column
*/
doublereal * BandMatrix::ptrColumn(int j) {
return m_colPtrs[j];
}
//====================================================================================================================
// Return a vector of const pointers to the columns
/*
* Note the value of the pointers are protected by their being const.
* However, the value of the matrix is open to being changed.
*
* @return returns a vector of pointers to the top of the columns
* of the matrices.
*/
doublereal * const * BandMatrix::colPts() {
return &(m_colPtrs[0]);
}
//====================================================================================================================
// Copy the data from one array into another without doing any checking
/*
* This differs from the assignment operator as no resizing is done and memcpy() is used.
* @param y Array to be copied
*/
void BandMatrix::copyData(const GeneralMatrix& y) {
m_factored = false;
size_t n = sizeof(doublereal) * m_n * (2 *m_kl + m_ku + 1);
GeneralMatrix * yyPtr = const_cast<GeneralMatrix *>(&y);
(void) memcpy(DATA_PTR(data), yyPtr->ptrColumn(0), n);
}
//====================================================================================================================
/*
* clear the factored flag
*/
void BandMatrix::clearFactorFlag() {
m_factored = 0;
}
//====================================================================================================================
//====================================================================================================================
}

195
Cantera/src/numerics/BandMatrix.h
View file

@ -1,7 +1,8 @@
/**
* @file BandMatrix.h
*
* Banded matrices.
* Declarations for the class BandMatrix
* which is a child class of GeneralMatrix for banded matrices handled by solvers
* (see class \ref numerics and \link Cantera::BandMatrix BandMatrix\endlink).
*/
/*
@ -20,14 +21,26 @@
#include "ctlapack.h"
#include "utilities.h"
#include "ctexceptions.h"
#include "GeneralMatrix.h"
namespace Cantera {
/**
* A class for banded matrices. This class has matrix inversion processes.
* The class is based upon the LAPACK banded storage matrix format.
//! A class for banded matrices, involving matrix inversion processes.
//! The class is based upon the LAPACK banded storage matrix format.
/*!
* An important issue with this class is that it stores both the original data
* and the LU factorization of the data. This means that the banded matrix typically
* will take up twice the room that it is expected to take.
*
* QR factorizations of banded matrices are not included in the original LAPACK work.
* Add-ons are available. However, they are not included here. Instead we just use the
* stock LU decompositions.
*
* This class is a derived class of the base class GeneralMatrix. However, withinin
* the oneD directory, the class is used as is, without reference to the GeneralMatrix
* base type.
*/
class BandMatrix {
class BandMatrix : public GeneralMatrix {
public:
@ -90,7 +103,7 @@ namespace Cantera {
doublereal& operator()(int i, int j);
//! Constant Index into the (i,j) element
//! Constant index into the (i,j) element
/*!
* @param i row
* @param j column
@ -144,7 +157,19 @@ namespace Cantera {
doublereal _value(int i, int j) const;
//! Returns the number of rows
int nRows() const;
virtual size_t nRows() const;
//! Return the size and structure of the matrix
/*!
* This is inherited from GeneralMatrix
*
* @param iStruct OUTPUT Pointer to a vector of ints that describe the structure of the matrix.
* istruct[0] = kl
* istruct[1] = ku
*
* @return returns the number of rows and columns in the matrix.
*/
virtual size_t nRowsAndStruct(int * const iStruct = 0) const;
//! Number of columns
int nColumns() const;
@ -169,14 +194,14 @@ namespace Cantera {
* @param b Vector to do the rh multiplcation
* @param prod OUTPUT vector to receive the result
*/
void mult(const doublereal * const b, doublereal * const prod) const;
virtual void mult(const doublereal * const b, doublereal * const prod) const;
//! Multiply b*A and write result to prod.
/*!
* @param b Vector to do the lh multiplcation
* @param prod OUTPUT vector to receive the result
*/
void leftMult(const doublereal * const b, doublereal * const prod) const;
virtual void leftMult(const doublereal * const b, doublereal * const prod) const;
//! Perform an LU decomposition, the LAPACK routine DGBTRF is used.
/*!
@ -192,7 +217,6 @@ namespace Cantera {
//! Solve the matrix problem Ax = b
/*!
* @param n size of the matrix
* @param b INPUT rhs of the problem
* @param x OUTPUT solution to the problem
*
@ -200,11 +224,10 @@ namespace Cantera {
* 0 indicates a success
* ~0 Some error occurred, see the LAPACK documentation
*/
int solve(int n, const doublereal * const b, doublereal * const x);
int solve(const doublereal * const b, doublereal * const x);
//! Solve the matrix problem Ax = b
/*!
* @param n size of the matrix
* @param b INPUT rhs of the problem
* OUTPUT solution to the problem
*
@ -212,14 +235,14 @@ namespace Cantera {
* 0 indicates a success
* ~0 Some error occurred, see the LAPACK documentation
*/
int solve(int n, doublereal * const b);
int solve(doublereal * const b);
//! Returns an iterator for the start of the band storage data
/*!
* Iterator points to the beginning of the data, and it is changeable.
*/
vector_fp::iterator begin();
virtual vector_fp::iterator begin();
//! Returns an iterator for the end of the band storage data
/*!
@ -239,6 +262,130 @@ namespace Cantera {
*/
vector_fp::const_iterator end() const;
/**
* Zero the matrix
*/
virtual void zero();
//! Factors the A matrix using the QR algorithm, overwriting A
/*!
* we set m_factored to 2 to indicate the matrix is now QR factored
*
* @return Returns the info variable from lapack
*/
virtual 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
*/
virtual 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
*/
virtual doublereal rcond(doublereal a1norm);
//! Change the way the matrix is factored
/*!
* @param fAlgorithm integer
* 0 LU factorization
* 1 QR factorization
*/
virtual void useFactorAlgorithm(int fAlgorithm);
//! Returns the factor algorithm used
/*!
* 0 LU decomposition
* 1 QR decomposition
*
* This routine will always return 0
*/
virtual int factorAlgorithm() const;
//! Returns the one norm of the matrix
virtual doublereal oneNorm() const;
//! Duplicate this object as a GeneralMatrix pointer
virtual GeneralMatrix * duplMyselfAsGeneralMatrix() const;
//! Report whether the current matrix has been factored.
virtual bool factored() const;
//! Return a pointer to the top of column j, column values are assumed to be contiguous in memory
/*!
* The LAPACK bandstructure has column values which are contiguous in memory:
*
* On entry, the matrix A in band storage, in rows KL+1 to
* 2*KL+KU+1; rows 1 to KL of the array need not be set.
* The j-th column of A is stored in the j-th column of the
* array AB as follows:
* AB(KL + KU + 1 + i - j,j) = A(i,j) for max(1, j - KU) <= i <= min(m, j + KL)
*
* This routine returns the position of AB(1,j) (fortran-1 indexing) in the above format
*
* So to address the (i,j) position, you use the following indexing:
*
* double *colP_j = matrix.ptrColumn(j);
* double a_i_j = colP_j[kl + ku + i - j];
*
*
* @param j Value of the column
*
* @return Returns a pointer to the top of the column
*/
virtual doublereal * ptrColumn(int j);
//! Return a vector of const pointers to the columns
/*!
* Note the value of the pointers are protected by their being const.
* However, the value of the matrix is open to being changed.
*
* @return returns a vector of pointers to the top of the columns
* of the matrices.
*/
virtual doublereal * const * colPts();
//! Copy the data from one array into another without doing any checking
/*!
* This differs from the assignment operator as no resizing is done and memcpy() is used.
* @param y Array to be copied
*/
virtual void copyData(const GeneralMatrix& y);
//! Clear the factored flag
virtual void clearFactorFlag();
//! Check to see if we have any zero rows in the jacobian
/*!
* This utility routine checks to see if any rows are zero.
* The smallest row is returned along with the largest coefficient in that row
*
* @param valueSmall OUTPUT value of the largest coefficient in the smallest row
*
* @return index of the row that is most nearly zero
*/
virtual int checkRows(doublereal & valueSmall) const;
//! Check to see if we have any zero columns in the jacobian
/*!
* This utility routine checks to see if any columns are zero.
* The smallest column is returned along with the largest coefficient in that column
*
* @param valueSmall OUTPUT value of the largest coefficient in the smallest column
*
* @return index of the column that is most nearly zero
*/
virtual int checkColumns(doublereal & valueSmall) const;
protected:
//! Matrix data
@ -265,6 +412,24 @@ namespace Cantera {
//! Pivot vector
vector_int m_ipiv;
//! Vector of column pointers
std::vector<doublereal *> m_colPtrs;
//! Extra work array needed - size = n
vector_int iwork_;
//! Extra dp work array needed - size = 3n
vector_fp work_;
private:
//! Error function that gets called for unhandled cases
/*!
* @param msg String containing the message.
*/
void err(std::string msg) const;
};
//! Utility routine to print out the matrix

2
Cantera/src/numerics/DenseMatrix.cpp
View file

@ -32,7 +32,7 @@ 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)

3
Cantera/src/numerics/DenseMatrix.h
View file

@ -21,6 +21,7 @@
#include "ct_defs.h"
#include "Array.h"
namespace Cantera {
/**
* @defgroup numerics Numerical Utilities within Cantera
@ -123,7 +124,7 @@ namespace Cantera {
* @return returns a vector of pointers to the top of the columns
* of the matrices.
*/
doublereal * const * colPts();
virtual doublereal * const * colPts();
//! Return a const vector of const pointers to the columns
/*!

42
Cantera/src/numerics/GeneralMatrix.cpp Normal file
View file

@ -0,0 +1,42 @@
/**
* @file GeneralMatrix.cpp
*
*/
/*
* $Revision: 725 $
* $Date: 2011-05-16 18:45:08 -0600 (Mon, 16 May 2011) $
*/
/*
* Copywrite 2004 Sandia Corporation. Under the terms of Contract
* DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
* retains certain rights in this software.
* See file License.txt for licensing information.
*/
#include "GeneralMatrix.h"
using namespace std;
namespace Cantera {
//====================================================================================================================
GeneralMatrix::GeneralMatrix(int matType) :
matrixType_(matType)
{
}
//====================================================================================================================
GeneralMatrix::GeneralMatrix(const GeneralMatrix &y) :
matrixType_(y.matrixType_)
{
}
//====================================================================================================================
GeneralMatrix& GeneralMatrix::operator=(const GeneralMatrix &y)
{
if (&y == this) return *this;
matrixType_ = y.matrixType_;
return *this;
}
//====================================================================================================================
GeneralMatrix::~GeneralMatrix()
{
}
//====================================================================================================================
}

245
Cantera/src/numerics/GeneralMatrix.h Normal file
View file

@ -0,0 +1,245 @@
/**
* @file GeneralMatrix.h
* Declarations for the class GeneralMatrix which is a virtual base class for matrices handled by solvers
* (see class \ref numerics and \link Cantera::GeneralMatrix GeneralMatrix\endlink).
*/
/*
* $Date: 2011-10-13 15:16:06 -0600 (Thu, 13 Oct 2011) $
* $Revision: 776 $
*/
/*
* Copywrite 2004 Sandia Corporation. Under the terms of Contract
* DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
* retains certain rights in this software.
* See file License.txt for licensing information.
*/
#ifndef CT_GENERALMATRIX_H
#define CT_GENERALMATRIX_H
#include "ct_defs.h"
namespace Cantera {
//! Generic matrix
class GeneralMatrix {
public:
//! Base Constructor
/*!
* @param matType Matrix type
* 0 full
* 1 banded
*/
GeneralMatrix(int matType);
//! Copy Constructor
/*!
* @param right Object to be copied
*/
GeneralMatrix(const GeneralMatrix& right);
//! Assignment operator
/*!
* @param right Object to be copied
*/
GeneralMatrix& operator=(const GeneralMatrix& right);
//! Destructor. Does nothing.
virtual ~GeneralMatrix();
//! Duplicator member function
/*!
* This function will duplicate the matrix given a generic GeneralMatrix pointer
*
* @return Returns a pointer to the malloced object
*/
virtual GeneralMatrix * duplMyselfAsGeneralMatrix() const = 0;
//! Zero the matrix elements
virtual void zero() = 0;
//! Multiply A*b and write result to prod.
/*!
* @param b Vector to do the rh multiplcation
* @param prod OUTPUT vector to receive the result
*/
virtual void mult(const doublereal * const b, doublereal * const prod) const = 0;
//! Multiply b*A and write result to prod.
/*!
* @param b Vector to do the lh multiplcation
* @param prod OUTPUT vector to receive the result
*/
virtual void leftMult(const doublereal * const b, doublereal * const prod) const = 0;
//! Factors the A matrix, overwriting A.
/*
* We flip m_factored boolean to indicate that the matrix is now A-1.
*/
virtual int factor() = 0;
//! Factors the A matrix using the QR algorithm, overwriting A
/*!
* we set m_factored to 2 to indicate the matrix is now QR factored
*
* @return Returns the info variable from lapack
*/
virtual int factorQR() = 0;
//! 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
*/
virtual doublereal rcondQR() = 0;
//! 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
*/
virtual doublereal rcond(doublereal a1norm) = 0;
//! Change the way the matrix is factored
/*!
* @param fAlgorithm integer
* 0 LU factorization
* 1 QR factorization
*/
virtual void useFactorAlgorithm(int fAlgorithm) = 0;
//! Return the factor algorithm used
/*!
*
*/
virtual int factorAlgorithm() const = 0;
//! Calculate the one norm of the matrix
/*!
* Returns the one norm of the matrix
*/
virtual doublereal oneNorm() const = 0;
//! Return the number of rows in the matrix
virtual size_t nRows() const = 0;
//! Return the size and structure of the matrix
/*!
* This is inherited from GeneralMatrix
*
* @param iStruct OUTPUT Pointer to a vector of ints that describe the structure of the matrix.
*
* @return returns the number of rows and columns in the matrix.
*/
virtual size_t nRowsAndStruct(int * const iStruct = 0) const = 0;
//! clear the factored flag
virtual void clearFactorFlag() = 0;
//! Solves the Ax = b system returning x in the b spot.
/*!
* @param b Vector for the rhs of the equation system
*/
virtual int solve(doublereal *b) = 0;
//! true if the current factorization is up to date with the matrix
virtual bool factored() const = 0;
//! Return a pointer to the top of column j, columns are assumed to be contiguous in memory
/*!
* @param j Value of the column
*
* @return Returns a pointer to the top of the column
*/
virtual doublereal * ptrColumn(int j) = 0;
//! Index into the (i,j) element
/*!
* @param i row
* @param j column
*
* Returns a changeable reference to the matrix entry
*/
virtual doublereal& operator()(int i, int j) = 0;
//! Constant Index into the (i,j) element
/*!
* @param i row
* @param j column
*
* Returns an unchangeable reference to the matrix entry
*/
virtual doublereal operator() (int i, int j) const = 0;
//! Copy the data from one array into another without doing any checking
/*!
* This differs from the assignment operator as no resizing is done and memcpy() is used.
* @param y Array to be copied
*/
virtual void copyData(const GeneralMatrix& y) = 0;
//! Return an iterator pointing to the first element
/*!
* We might drop this later
*/
virtual vector_fp::iterator begin() = 0;
//! Return a const iterator pointing to the first element
/*!
* We might drop this later
*/
virtual vector_fp::const_iterator begin() const = 0;
//! Return a vector of const pointers to the columns
/*!
* Note the value of the pointers are protected by their being const.
* However, the value of the matrix is open to being changed.
*
* @return returns a vector of pointers to the top of the columns
* of the matrices.
*/
virtual doublereal * const * colPts() = 0;
//! Check to see if we have any zero rows in the jacobian
/*!
* This utility routine checks to see if any rows are zero.
* The smallest row is returned along with the largest coefficient in that row
*
* @param valueSmall OUTPUT value of the largest coefficient in the smallest row
*
* @return index of the row that is most nearly zero
*/
virtual int checkRows (doublereal & valueSmall) const = 0;
//! Check to see if we have any zero columns in the jacobian
/*!
* This utility routine checks to see if any columns are zero.
* The smallest column is returned along with the largest coefficient in that column
*
* @param valueSmall OUTPUT value of the largest coefficient in the smallest column
*
* @return index of the column that is most nearly zero
*/
virtual int checkColumns (doublereal & valueSmall) const = 0;
//! Matrix type
/*!
* 0 Square
* 1 Banded
*/
int matrixType_;
};
}
#endif

4
Cantera/src/numerics/Makefile.in
View file

@ -37,14 +37,14 @@ CXX_FLAGS = @CXXFLAGS@ $(LOCAL_DEFS) $(CXX_OPT) $(PIC_FLAG) $(DEBUG_FLAG)
NUMERICS_OBJ = DenseMatrix.o funcs.o Func1.o \
ODE_integrators.o BandMatrix.o DAE_solvers.o \
funcs.o sort.o SquareMatrix.o ResidJacEval.o NonlinearSolver.o \
solveProb.o BEulerInt.o RootFind.o IDA_Solver.o
solveProb.o BEulerInt.o RootFind.o IDA_Solver.o GeneralMatrix.o
NUMERICS_H = ArrayViewer.h DenseMatrix.h \
funcs.h ctlapack.h Func1.h FuncEval.h \
polyfit.h\
BandMatrix.h Integrator.h DAE_Solver.h ResidEval.h sort.h \
SquareMatrix.h ResidJacEval.h NonlinearSolver.h \
solveProb.h BEulerInt.h RootFind.h IDA_Solver.h
solveProb.h BEulerInt.h RootFind.h IDA_Solver.h GeneralMatrix.h
ifeq ($(use_sundials), 1)
ODEPACKAGE_H = CVodesIntegrator.h

519
Cantera/src/numerics/NonlinearSolver.cpp
View file

@ -19,6 +19,7 @@
#include <limits>
#include "SquareMatrix.h"
#include "GeneralMatrix.h"
#include "NonlinearSolver.h"
#include "ctlapack.h"
@ -147,8 +148,8 @@ namespace Cantera {
atolk_(0),
m_print_flag(0),
m_ScaleSolnNormToResNorm(0.001),
jacCopy_(0),
Hessian_(0),
jacCopyPtr_(0),
HessianPtr_(0),
deltaX_CP_(0),
deltaX_Newton_(0),
residNorm2Cauchy_(0.0),
@ -211,7 +212,7 @@ namespace Cantera {
}
jacCopy_.resize(neq_, neq_, 0.0);
// jacCopyPtr_->resize(neq_, 0.0);
deltaX_CP_.resize(neq_, 0.0);
Jd_.resize(neq_, 0.0);
deltaX_trust_.resize(neq_, 1.0);
@ -268,8 +269,8 @@ namespace Cantera {
atolk_(0),
m_print_flag(0),
m_ScaleSolnNormToResNorm(0.001),
jacCopy_(0),
Hessian_(0),
jacCopyPtr_(0),
HessianPtr_(0),
deltaX_CP_(0),
deltaX_Newton_(0),
residNorm2Cauchy_(0.0),
@ -306,6 +307,12 @@ namespace Cantera {
//====================================================================================================================
NonlinearSolver::~NonlinearSolver() {
if (jacCopyPtr_) {
delete jacCopyPtr_;
}
if (HessianPtr_) {
delete HessianPtr_;
}
}
//====================================================================================================================
NonlinearSolver& NonlinearSolver::operator=(const NonlinearSolver &right) {
@ -364,8 +371,15 @@ namespace Cantera {
m_print_flag = right.m_print_flag;
m_ScaleSolnNormToResNorm = right.m_ScaleSolnNormToResNorm;
jacCopy_ = right.jacCopy_;
Hessian_ = right.Hessian_;
if (jacCopyPtr_) {
delete (jacCopyPtr_);
}
jacCopyPtr_ = (right.jacCopyPtr_)->duplMyselfAsGeneralMatrix();
if (HessianPtr_) {
delete (HessianPtr_);
}
HessianPtr_ = (right.HessianPtr_)->duplMyselfAsGeneralMatrix();
deltaX_CP_ = right.deltaX_CP_;
deltaX_Newton_ = right.deltaX_Newton_;
residNorm2Cauchy_ = right.residNorm2Cauchy_;
@ -716,39 +730,57 @@ namespace Cantera {
* @param ydot_comm Current value of the time derivative of the solution vector
* @param time_curr current value of the time
*/
void NonlinearSolver::scaleMatrix(SquareMatrix& jac, doublereal * const y_comm, doublereal * const ydot_comm,
void NonlinearSolver::scaleMatrix(GeneralMatrix& jac, doublereal * const y_comm, doublereal * const ydot_comm,
doublereal time_curr, int num_newt_its)
{
int irow, jcol;
int ku, kl;
int ivec[2];
int n = jac.nRowsAndStruct(ivec);
double *colP_j;
/*
* Column scaling -> We scale the columns of the Jacobian
* by the nominal important change in the solution vector
*/
if (m_colScaling) {
if (!jac.m_factored) {
/*
* Go get new scales -> Took this out of this inner loop.
* Needs to be done at a larger scale.
*/
// setColumnScales();
if (!jac.factored()) {
if (jac.matrixType_ == 0) {
/*
* Go get new scales -> Took this out of this inner loop.
* Needs to be done at a larger scale.
*/
// setColumnScales();
/*
* Scale the new Jacobian
*/
doublereal *jptr = &(*(jac.begin()));
for (jcol = 0; jcol < neq_; jcol++) {
for (irow = 0; irow < neq_; irow++) {
*jptr *= m_colScales[jcol];
jptr++;
/*
* Scale the new Jacobian
*/
doublereal *jptr = &(*(jac.begin()));
for (jcol = 0; jcol < neq_; jcol++) {
for (irow = 0; irow < neq_; irow++) {
*jptr *= m_colScales[jcol];
jptr++;
}
}
} else if (jac.matrixType_ == 1) {
kl = ivec[0];
ku = ivec[1];
for (jcol = 0; jcol < neq_; jcol++) {
colP_j = (doublereal *) jac.ptrColumn(jcol);
for (irow = jcol - ku; irow <= jcol + kl; irow++) {
if (irow >= 0 && irow < neq_) {
colP_j[kl + ku + irow - jcol] *= m_colScales[jcol];
}
}
}
}
}
}
}
/*
* row sum scaling -> Note, this is an unequivical success
* at keeping the small numbers well balanced and nonnegative.
*/
if (! jac.m_factored) {
if (! jac.factored()) {
/*
* Ok, this is ugly. jac.begin() returns an vector<double> iterator
* to the first data location.
@ -759,39 +791,75 @@ namespace Cantera {
m_rowScales[irow] = 0.0;
m_rowWtScales[irow] = 0.0;
}
for (jcol = 0; jcol < neq_; jcol++) {
for (irow = 0; irow < neq_; irow++) {
if (m_rowScaling) {
m_rowScales[irow] += fabs(*jptr);
}
if (m_colScaling) {
// This is needed in order to mitgate the change in J_ij carried out just above this loop.
// Alternatively, we could move this loop up to the top
m_rowWtScales[irow] += fabs(*jptr) * m_ewt[jcol] / m_colScales[jcol];
} else {
m_rowWtScales[irow] += fabs(*jptr) * m_ewt[jcol];
if (jac.matrixType_ == 0) {
for (jcol = 0; jcol < neq_; jcol++) {
for (irow = 0; irow < neq_; irow++) {
if (m_rowScaling) {
m_rowScales[irow] += fabs(*jptr);
}
if (m_colScaling) {
// This is needed in order to mitgate the change in J_ij carried out just above this loop.
// Alternatively, we could move this loop up to the top
m_rowWtScales[irow] += fabs(*jptr) * m_ewt[jcol] / m_colScales[jcol];
} else {
m_rowWtScales[irow] += fabs(*jptr) * m_ewt[jcol];
}
jptr++;
}
}
} else if (jac.matrixType_ == 1) {
kl = ivec[0];
ku = ivec[1];
for (jcol = 0; jcol < neq_; jcol++) {
colP_j = (doublereal *) jac.ptrColumn(jcol);
for (irow = jcol - ku; irow <= jcol + kl; irow++) {
if (irow >= 0 && irow < neq_) {
double vv = fabs(colP_j[kl + ku + irow - jcol]);
if (m_rowScaling) {
m_rowScales[irow] += vv;
}
if (m_colScaling) {
// This is needed in order to mitgate the change in J_ij carried out just above this loop.
// Alternatively, we could move this loop up to the top
m_rowWtScales[irow] += vv * m_ewt[jcol] / m_colScales[jcol];
} else {
m_rowWtScales[irow] += vv * m_ewt[jcol];
}
}
}
jptr++;
}
}
if (m_rowScaling) {
for (irow = 0; irow < neq_; irow++) {
m_rowScales[irow] = 1.0/m_rowScales[irow];
}
for (irow = 0; irow < neq_; irow++) {
m_rowScales[irow] = 1.0/m_rowScales[irow];
}
} else {
for (irow = 0; irow < neq_; irow++) {
m_rowScales[irow] = 1.0;
}
for (irow = 0; irow < neq_; irow++) {
m_rowScales[irow] = 1.0;
}
}
// What we have defined is a maximum value that the residual can be and still pass.
// This isn't sufficient.
if (m_rowScaling) {
jptr = &(*(jac.begin()));
for (jcol = 0; jcol < neq_; jcol++) {
for (irow = 0; irow < neq_; irow++) {
*jptr *= m_rowScales[irow];
jptr++;
if (jac.matrixType_ == 0) {
jptr = &(*(jac.begin()));
for (jcol = 0; jcol < neq_; jcol++) {
for (irow = 0; irow < neq_; irow++) {
*jptr *= m_rowScales[irow];
jptr++;
}
}
} else if (jac.matrixType_ == 1) {
kl = ivec[0];
ku = ivec[1];
for (jcol = 0; jcol < neq_; jcol++) {
colP_j = (doublereal *) jac.ptrColumn(jcol);
for (irow = jcol - ku; irow <= jcol + kl; irow++) {
if (irow >= 0 && irow < neq_) {
colP_j[kl + ku + irow - jcol] *= m_rowScales[irow];
}
}
}
}
}
@ -812,7 +880,7 @@ namespace Cantera {
*/
void NonlinearSolver::calcSolnToResNormVector()
{
if (! jacCopy_.m_factored) {
if (! jacCopyPtr_->factored()) {
doublereal sum = 0.0;
@ -829,7 +897,7 @@ namespace Cantera {
for (int irow = 0; irow < neq_; irow++) {
m_wksp[irow] = 0.0;
}
doublereal *jptr = &(*(jacCopy_.begin()));
doublereal *jptr = &(jacCopyPtr_->operator()(0,0));
for (int jcol = 0; jcol < neq_; jcol++) {
for (int irow = 0; irow < neq_; irow++) {
m_wksp[irow] += (*jptr) * m_ewt[jcol];
@ -873,7 +941,7 @@ namespace Cantera {
*/
int NonlinearSolver::doNewtonSolve(const doublereal time_curr, const doublereal * const y_curr,
const doublereal * const ydot_curr, doublereal * const delta_y,
SquareMatrix& jac)
GeneralMatrix& jac)
{
int irow;
@ -961,16 +1029,22 @@ namespace Cantera {
* This is algorith A.6.5.1 in Dennis / Schnabel
*
* Compute the QR decomposition
*
* Notes on banded Hessian solve:
* The matrix for jT j has a larger band width. Both the top and bottom band widths
* are doubled, going from KU to KU+KL and KL to KU+KL in size. This is not an impossible increase in cost, but
* has to be considered.
*/
int NonlinearSolver::doAffineNewtonSolve(const doublereal * const y_curr, const doublereal * const ydot_curr,
doublereal * const delta_y, SquareMatrix& jac)
doublereal * const delta_y, GeneralMatrix& 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
jac.useFactorAlgorithm(1);
int useQR = jac.factorAlgorithm();
// multiplyl 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++) {
@ -986,8 +1060,8 @@ namespace Cantera {
// Factor the matrix using a standard Newton solve
m_conditionNumber = 1.0E300;
int info = 0;
if (!jac.m_factored) {
if (jac.useQR_) {
if (!jac.factored()) {
if (useQR) {
info = jac.factorQR();
} else {
info = jac.factor();
@ -999,14 +1073,27 @@ namespace Cantera {
*/
if (info == 0) {
doublereal rcond = 0.0;
if (jac.useQR_) {
if (useQR) {
rcond = jac.rcondQR();
} else {
rcond = jac.rcond(jac.a1norm_);
doublereal a1norm = jac.oneNorm();
rcond = jac.rcond(a1norm);
}
if (rcond > 0.0) {
m_conditionNumber = 1.0 / rcond;
}
} else {
m_conditionNumber = 1.0E300;
newtonGood = false;
if (m_print_flag >= 1) {
printf("\t\t doAffineNewtonSolve: ");
if (useQR) {
printf("factorQR()");
} else {
printf("factor()");
}
printf(" returned with info = %d, indicating a zero row or column\n", info);
}
}
bool doHessian = false;
if (s_doBothSolvesAndCompare) {
@ -1040,10 +1127,19 @@ namespace Cantera {
}
} else {
doHessian = true;
newtonGood = false;
if (m_print_flag >= 3) {
printf("\t\t doAffineNewtonSolve() WARNING: Condition number too large, %g. Doing a Hessian solve \n", m_conditionNumber);
if (jac.matrixType_ == 1) {
useNewton = true;
newtonGood = true;
if (m_print_flag >= 3) {
printf("\t\t doAffineNewtonSolve() WARNING: Condition number too large, %g, But Banded Hessian solve "
"not implemented yet \n", m_conditionNumber);
}
} else {
doHessian = true;
newtonGood = false;
if (m_print_flag >= 3) {
printf("\t\t doAffineNewtonSolve() WARNING: Condition number too large, %g. Doing a Hessian solve \n", m_conditionNumber);
}
}
}
@ -1056,30 +1152,32 @@ namespace Cantera {
}
// Get memory if not done before
if (Hessian_.nRows() == 0) {
Hessian_.resize(neq_, neq_);
if (HessianPtr_ == 0) {
HessianPtr_ = jac.duplMyselfAsGeneralMatrix();
}
/*
* Calculate the symmetric Hessian
*/
Hessian_.zero();
GeneralMatrix &hessian = *HessianPtr_;
GeneralMatrix &jacCopy = *jacCopyPtr_;
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(i,j) += jacCopy(k,i) * jacCopy(k,j) * m_rowScales[k] * m_rowScales[k];
}
Hessian_(j,i) = Hessian_(i,j);
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(i,j) += jacCopy(k,i) * jacCopy(k,j);
}
Hessian_(j,i) = Hessian_(i,j);
hessian(j,i) = hessian(i,j);
}
}
}
@ -1092,10 +1190,10 @@ namespace Cantera {
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];
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 += fabs(hessian(i,j)) * m_colScales[j];
}
hcol *= m_colScales[i];
if (hcol > hnorm) {
@ -1105,10 +1203,10 @@ namespace Cantera {
} else {
for (int i = 0; i < neq_; i++) {
for (int j = i; j < neq_; j++) {
hcol += fabs(Hessian_(j,i));
hcol += fabs(hessian(j,i));
}
for (int j = i+1; j < neq_; j++) {
hcol += fabs(Hessian_(i,j));
hcol += fabs(hessian(i,j));
}
if (hcol > hnorm) {
hnorm = hcol;
@ -1127,11 +1225,11 @@ namespace Cantera {
#endif
if (m_colScaling) {
for (int i = 0; i < neq_; i++) {
Hessian_(i,i) += hcol / (m_colScales[i] * m_colScales[i]);
hessian(i,i) += hcol / (m_colScales[i] * m_colScales[i]);
}
} else {
for (int i = 0; i < neq_; i++) {
Hessian_(i,i) += hcol;
hessian(i,i) += hcol;
}
}
@ -1139,7 +1237,7 @@ namespace Cantera {
* Factor the Hessian
*/
int info;
ct_dpotrf(ctlapack::UpperTriangular, neq_, &(*(Hessian_.begin())), neq_, info);
ct_dpotrf(ctlapack::UpperTriangular, neq_, &(*(HessianPtr_->begin())), neq_, info);
if (info) {
if (m_print_flag >= 2) {
printf("\t\t doAffineNewtonSolve() ERROR: Hessian isn't positive definate DPOTRF returned INFO = %d\n", info);
@ -1164,14 +1262,14 @@ namespace Cantera {
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];
delta_y[j] += delyH[i] * jacCopy(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);
delta_y[j] += delyH[i] * jacCopy(i,j);
}
}
}
@ -1180,7 +1278,7 @@ namespace Cantera {
/*
* Solve the factored Hessian System
*/
ct_dpotrs(ctlapack::UpperTriangular, neq_, 1,&(*(Hessian_.begin())), neq_, delta_y, neq_, info);
ct_dpotrs(ctlapack::UpperTriangular, neq_, 1,&(*(hessian.begin())), neq_, delta_y, neq_, info);
if (info) {
if (m_print_flag >= 2) {
printf("\t\t NonlinearSolver::doAffineNewtonSolve() ERROR: DPOTRS returned INFO = %d\n", info);
@ -1287,7 +1385,7 @@ namespace Cantera {
/*
* This call must be made on the unfactored jacobian!
*/
doublereal NonlinearSolver::doCauchyPointSolve(SquareMatrix& jac)
doublereal NonlinearSolver::doCauchyPointSolve(GeneralMatrix& jac)
{
doublereal rowFac = 1.0;
doublereal colFac = 1.0;
@ -1311,7 +1409,7 @@ namespace Cantera {
if (m_rowScaling) {
rowFac = 1.0 / m_rowScales[i];
}
deltaX_CP_[j] -= m_resid[i] * jac.value(i,j) * colFac * rowFac * m_ewt[j] * m_ewt[j]
deltaX_CP_[j] -= m_resid[i] * jac(i,j) * colFac * rowFac * m_ewt[j] * m_ewt[j]
/ (m_residWts[i] * m_residWts[i]);
#ifdef DEBUG_MODE
mdp::checkFinite(deltaX_CP_[j]);
@ -1333,7 +1431,7 @@ namespace Cantera {
if (m_colScaling) {
colFac = 1.0 / m_colScales[j];
}
Jd_[i] += deltaX_CP_[j] * jac.value(i,j) * rowFac * colFac / m_residWts[i];
Jd_[i] += deltaX_CP_[j] * jac(i,j) * rowFac * colFac / m_residWts[i];
}
}
@ -2322,7 +2420,7 @@ namespace Cantera {
int NonlinearSolver::dampStep(const doublereal time_curr, const doublereal * const y_n_curr,
const doublereal * const ydot_n_curr, doublereal * const step_1,
doublereal * const y_n_1, doublereal * const ydot_n_1, doublereal * const step_2,
doublereal & stepNorm_2, SquareMatrix& jac, bool writetitle, int& num_backtracks)
doublereal & stepNorm_2, GeneralMatrix& jac, bool writetitle, int& num_backtracks)
{
int j, m;
int info = 0;
@ -2552,7 +2650,7 @@ namespace Cantera {
int NonlinearSolver::dampDogLeg(const doublereal time_curr, const doublereal* y_n_curr,
const doublereal *ydot_n_curr, std::vector<doublereal> & step_1,
doublereal* const y_n_1, doublereal* const ydot_n_1,
doublereal& stepNorm_1, doublereal& stepNorm_2, SquareMatrix& jac, int& numTrials)
doublereal& stepNorm_1, doublereal& stepNorm_2, GeneralMatrix& jac, int& numTrials)
{
doublereal lambda;
int info;
@ -2907,7 +3005,7 @@ namespace Cantera {
* -1 Failed convergence
*/
int NonlinearSolver::solve_nonlinear_problem(int SolnType, doublereal * const y_comm, doublereal * const ydot_comm,
doublereal CJ, doublereal time_curr, SquareMatrix& jac,
doublereal CJ, doublereal time_curr, GeneralMatrix& jac,
int &num_newt_its, int &num_linear_solves,
int &num_backtracks, int loglevelInput)
{
@ -2921,6 +3019,11 @@ namespace Cantera {
int retnDamp = 0;
int retnCode = 0;
bool forceNewJac = false;
if (jacCopyPtr_) {
delete jacCopyPtr_;
}
jacCopyPtr_ = jac.duplMyselfAsGeneralMatrix();
doublereal stepNorm_1;
doublereal stepNorm_2;
@ -2945,11 +3048,7 @@ namespace Cantera {
num_backtracks = 0;
int i_numTrials;
m_print_flag = loglevelInput;
if (m_print_flag > 1) {
jac.m_printLevel = 1;
} else {
jac.m_printLevel = 0;
}
if (trustRegionInitializationMethod_ == 0) {
trInit = true;
} else if (trustRegionInitializationMethod_ == 1) {
@ -3563,10 +3662,9 @@ namespace Cantera {
* 1 Means a successful operation
* 0 Means an unsuccessful operation
*/
int NonlinearSolver::beuler_jac(SquareMatrix &J, doublereal * const f,
int NonlinearSolver::beuler_jac(GeneralMatrix &J, doublereal * const f,
doublereal time_curr, doublereal CJ,
doublereal * const y,
doublereal * const ydot,
doublereal * const y, doublereal * const ydot,
int num_newt_its)
{
int i, j;
@ -3590,101 +3688,190 @@ namespace Cantera {
return info;
}
} else {
/*******************************************************************
* Generic algorithm to calculate a numerical Jacobian
*/
/*
* Calculate the current value of the rhs given the
* current conditions.
*/
if (J.matrixType_ == 0) {
/*******************************************************************
* Generic algorithm to calculate a numerical Jacobian
*/
/*
* Calculate the current value of the rhs given the
* current conditions.
*/
info = m_func->evalResidNJ(time_curr, delta_t_n, y, ydot, f, JacBase_ResidEval);
m_nfe++;
if (info != 1) {
return info;
}
m_nJacEval++;
info = m_func->evalResidNJ(time_curr, delta_t_n, y, ydot, f, JacBase_ResidEval);
m_nfe++;
if (info != 1) {
return info;
}
m_nJacEval++;
/*
* Malloc a vector and call the function object to return a set of
* deltaY's that are appropriate for calculating the numerical
* derivative.
*/
doublereal *dyVector = mdp::mdp_alloc_dbl_1(neq_, MDP_DBL_NOINIT);
retn = m_func->calcDeltaSolnVariables(time_curr, y, ydot, dyVector, DATA_PTR(m_ewt));
/*
* Malloc a vector and call the function object to return a set of
* deltaY's that are appropriate for calculating the numerical
* derivative.
*/
doublereal *dyVector = mdp::mdp_alloc_dbl_1(neq_, MDP_DBL_NOINIT);
retn = m_func->calcDeltaSolnVariables(time_curr, y, ydot, dyVector, DATA_PTR(m_ewt));
if (s_print_NumJac) {
if (m_print_flag >= 7) {
if (neq_ < 20) {
printf("\t\tUnk m_ewt y dyVector ResN\n");
for (int iii = 0; iii < neq_; iii++){
printf("\t\t %4d %16.8e %16.8e %16.8e %16.8e \n",
iii, m_ewt[iii], y[iii], dyVector[iii], f[iii]);
if (s_print_NumJac) {
if (m_print_flag >= 7) {
if (neq_ < 20) {
printf("\t\tUnk m_ewt y dyVector ResN\n");
for (int iii = 0; iii < neq_; iii++){
printf("\t\t %4d %16.8e %16.8e %16.8e %16.8e \n",
iii, m_ewt[iii], y[iii], dyVector[iii], f[iii]);
}
}
}
}
}
/*
* Loop over the variables, formulating a numerical derivative
* of the dense matrix.
* For the delta in the variable, we will use a variety of approaches
* The original approach was to use the error tolerance amount.
* This may not be the best approach, as it could be overly large in
* some instances and overly small in others.
* We will first protect from being overly small, by using the usual
* sqrt of machine precision approach, i.e., 1.0E-7,
* to bound the lower limit of the delta.
*/
for (j = 0; j < neq_; j++) {
/*
* Loop over the variables, formulating a numerical derivative
* of the dense matrix.
* For the delta in the variable, we will use a variety of approaches
* The original approach was to use the error tolerance amount.
* This may not be the best approach, as it could be overly large in
* some instances and overly small in others.
* We will first protect from being overly small, by using the usual
* sqrt of machine precision approach, i.e., 1.0E-7,
* to bound the lower limit of the delta.
*/
for (j = 0; j < neq_; j++) {
/*
* Get a pointer into the column of the matrix
*/
/*
* Get a pointer into the column of the matrix
*/
col_j = (doublereal *) J.ptrColumn(j);
ysave = y[j];
dy = dyVector[j];
//dy = fmaxx(1.0E-6 * m_ewt[j], fabs(ysave)*1.0E-7);
col_j = (doublereal *) J.ptrColumn(j);
ysave = y[j];
dy = dyVector[j];
//dy = fmaxx(1.0E-6 * m_ewt[j], fabs(ysave)*1.0E-7);
y[j] = ysave + dy;
dy = y[j] - ysave;
if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
ydotsave = ydot[j];
ydot[j] += dy * CJ;
}
/*
* Call the function
*/
info = m_func->evalResidNJ(time_curr, delta_t_n, y, ydot, DATA_PTR(m_wksp),
JacDelta_ResidEval, j, dy);
m_nfe++;
if (info != 1) {
mdp::mdp_safe_free((void **) &dyVector);
return info;
}
doublereal diff;
for (i = 0; i < neq_; i++) {
diff = subtractRD(m_wksp[i], f[i]);
col_j[i] = diff / dy;
}
y[j] = ysave;
if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
ydot[j] = ydotsave;
}
y[j] = ysave + dy;
dy = y[j] - ysave;
if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
ydotsave = ydot[j];
ydot[j] += dy * CJ;
}
/*
* Call the function
*/
/*
* Release memory
*/
mdp::mdp_safe_free((void **) &dyVector);
} else if (J.matrixType_ == 1) {
int ku, kl;
int ivec[2];
int n = J.nRowsAndStruct(ivec);
kl = ivec[0];
ku = ivec[1];
if (n != neq_) {
printf("we have probs\n"); exit(-1);
}
info = m_func->evalResidNJ(time_curr, delta_t_n, y, ydot, DATA_PTR(m_wksp),
JacDelta_ResidEval, j, dy);
m_nfe++;
// --------------------------------- BANDED MATRIX BRAIN DEAD ---------------------------------------------------
info = m_func->evalResidNJ(time_curr, delta_t_n, y, ydot, f, JacBase_ResidEval);
m_nfe++;
if (info != 1) {
mdp::mdp_safe_free((void **) &dyVector);
return info;
}
m_nJacEval++;
doublereal diff;
for (i = 0; i < neq_; i++) {
diff = subtractRD(m_wksp[i], f[i]);
col_j[i] = diff / dy;
}
y[j] = ysave;
if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
ydot[j] = ydotsave;
doublereal *dyVector = mdp::mdp_alloc_dbl_1(neq_, MDP_DBL_NOINIT);
retn = m_func->calcDeltaSolnVariables(time_curr, y, ydot, dyVector, DATA_PTR(m_ewt));
if (s_print_NumJac) {
if (m_print_flag >= 7) {
if (neq_ < 20) {
printf("\t\tUnk m_ewt y dyVector ResN\n");
for (int iii = 0; iii < neq_; iii++){
printf("\t\t %4d %16.8e %16.8e %16.8e %16.8e \n",
iii, m_ewt[iii], y[iii], dyVector[iii], f[iii]);
}
}
}
}
for (j = 0; j < neq_; j++) {
col_j = (doublereal *) J.ptrColumn(j);
ysave = y[j];
dy = dyVector[j];
y[j] = ysave + dy;
dy = y[j] - ysave;
if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
ydotsave = ydot[j];
ydot[j] += dy * CJ;
}
info = m_func->evalResidNJ(time_curr, delta_t_n, y, ydot, DATA_PTR(m_wksp), JacDelta_ResidEval, j, dy);
m_nfe++;
if (info != 1) {
mdp::mdp_safe_free((void **) &dyVector);
return info;
}
doublereal diff;
for (int i = j - ku; i <= j + kl; i++) {
if (i >= 0 && i < neq_) {
diff = subtractRD(m_wksp[i], f[i]);
col_j[kl + ku + i - j] = diff / dy;
}
}
y[j] = ysave;
if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
ydot[j] = ydotsave;
}
}
mdp::mdp_safe_free((void **) &dyVector);
double vSmall;
int ismall = J.checkRows(vSmall);
if (vSmall < 1.0E-100) {
printf("WE have a zero row, %d\n", ismall);
exit(-1);
}
ismall = J.checkColumns(vSmall);
if (vSmall < 1.0E-100) {
printf("WE have a zero column, %d\n", ismall);
exit(-1);
}
// ---------------------BANDED MATRIX BRAIN DEAD -----------------------
}
/*
* Release memory
*/
mdp::mdp_safe_free((void **) &dyVector);
}
if (m_print_flag >= 7 && s_print_NumJac) {
@ -3721,7 +3908,7 @@ namespace Cantera {
* 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);
jacCopyPtr_->copyData(J);
return retn;
}

20
Cantera/src/numerics/NonlinearSolver.h
View file

@ -259,7 +259,7 @@ namespace Cantera {
*/
int doNewtonSolve(const doublereal time_curr, const doublereal * const y_curr,
const doublereal * const ydot_curr, doublereal * const delta_y,
SquareMatrix& jac);
GeneralMatrix& jac);
//! Compute the newton step, either by direct newton's or by solving a close problem that is represented
//! by a Hessian (
@ -293,7 +293,7 @@ namespace Cantera {
* else indicates a failure.
*/
int doAffineNewtonSolve(const doublereal * const y_curr, const doublereal * const ydot_curr,
doublereal * const delta_y, SquareMatrix& jac);
doublereal * const delta_y, GeneralMatrix& jac);
//! Calculate the length of the current trust region in terms of the solution error norm
/*!
@ -438,7 +438,7 @@ namespace Cantera {
* 1 Means a successful operation
* 0 Means an unsuccessful operation
*/
int beuler_jac(SquareMatrix &J, doublereal * const f,
int beuler_jac(GeneralMatrix &J, doublereal * const f,
doublereal time_curr, doublereal CJ, doublereal * const y,
doublereal * const ydot, int num_newt_its);
@ -510,7 +510,7 @@ namespace Cantera {
int dampStep(const doublereal time_curr, const doublereal * const y_n_curr,
const doublereal * const ydot_n_curr, doublereal * const step_1,
doublereal * const y_n_1, doublereal * const ydot_n_1, doublereal * step_2,
doublereal & stepNorm_2, SquareMatrix& jac, bool writetitle,
doublereal & stepNorm_2, GeneralMatrix& jac, bool writetitle,
int& num_backtracks);
//! Find the solution to F(X) = 0 by damped Newton iteration.
@ -541,7 +541,7 @@ namespace Cantera {
* -1 Failed convergence
*/
int solve_nonlinear_problem(int SolnType, doublereal * const y_comm, doublereal * const ydot_comm, doublereal CJ,
doublereal time_curr, SquareMatrix& jac,int &num_newt_its,
doublereal time_curr, GeneralMatrix & jac, int &num_newt_its,
int &num_linear_solves, int &num_backtracks, int loglevelInput);
private:
@ -589,7 +589,7 @@ namespace Cantera {
* @param time_curr current value of the time
* @param num_newt_its Current value of the number of newt its
*/
void scaleMatrix(SquareMatrix& jac, doublereal * const y_comm, doublereal * const ydot_comm,
void scaleMatrix(GeneralMatrix& jac, doublereal * const y_comm, doublereal * const ydot_comm,
doublereal time_curr, int num_newt_its);
//! Print solution norm contribution
@ -689,7 +689,7 @@ namespace Cantera {
*
* @return Returns the norm of the solution update
*/
doublereal doCauchyPointSolve(SquareMatrix& jac);
doublereal doCauchyPointSolve(GeneralMatrix& jac);
//! This is a utility routine that can be used to print out the rates of the initial residual decline
/*!
@ -793,7 +793,7 @@ namespace Cantera {
int dampDogLeg(const doublereal time_curr, const doublereal* y_n_curr,
const doublereal *ydot_n_curr, std::vector<doublereal> & step_1,
doublereal* const y_n_1, doublereal* const ydot_n_1,
doublereal& stepNorm_1, doublereal& stepNorm_2, SquareMatrix& jac, int& num_backtracks);
doublereal& stepNorm_1, doublereal& stepNorm_2, GeneralMatrix& jac, int& num_backtracks);
//! Decide whether the current step is acceptable and adjust the trust region size
/*!
@ -1103,10 +1103,10 @@ namespace Cantera {
/*!
* The jacobian storred here is the raw matrix, before any row or column scaling is carried out
*/
Cantera::SquareMatrix jacCopy_;
Cantera::GeneralMatrix * jacCopyPtr_;
//! Hessian
Cantera::SquareMatrix Hessian_;
Cantera::GeneralMatrix * HessianPtr_;
/*********************************************************************************************
* VARIABLES ASSOCIATED WITH STEPS AND ASSOCIATED DOUBLE DOGLEG PARAMETERS

4
Cantera/src/numerics/ResidJacEval.cpp
View file

@ -333,7 +333,7 @@ namespace Cantera {
evalJacobian(const doublereal t, const doublereal delta_t, doublereal cj,
const doublereal * const y,
const doublereal * const ydot,
SquareMatrix &J,
GeneralMatrix &J,
doublereal * const resid)
{
doublereal * const * jac_colPts = J.colPts();
@ -349,7 +349,7 @@ namespace Cantera {
* @param c_j The current value of the coefficient of the time derivative
* @param y Solution vector (input, do not modify)
* @param ydot Rate of change of solution vector. (input, do not modify)
* @param jac_colPts Reference to the SquareMatrix object to be calculated (output)
* @param jac_colPts Reference to the SquareMatrix object to be calculated (output)
* @param resid Value of the residual that is computed (output)
*/
int ResidJacEval::

9
Cantera/src/numerics/ResidJacEval.h
View file

@ -21,7 +21,7 @@
#define CT_RESIDJACEVAL_H
#include "ResidEval.h"
#include "SquareMatrix.h"
#include "GeneralMatrix.h"
namespace Cantera {
@ -312,7 +312,7 @@ namespace Cantera {
virtual int matrixConditioning(doublereal * const matrix, const int nrows,
doublereal * const rhs);
//! Calculate an analytical jacobian and the residual at the current time and values.
//! Calculate an analytical jacobian and the residual at the current time and values.
/*!
* Only called if the jacFormation method is set to analytical
*
@ -329,8 +329,9 @@ namespace Cantera {
* -0 or neg value Means an unsuccessful operation
*/
virtual int evalJacobian(const doublereal t, const doublereal delta_t, doublereal cj,
const doublereal* const y, const doublereal* const ydot,
SquareMatrix &J, doublereal * const resid);
const doublereal* const y, const doublereal* const ydot,
GeneralMatrix &J, doublereal * const resid);
//! Calculate an analytical jacobian and the residual at the current time and values.
/*!

142
Cantera/src/numerics/SquareMatrix.cpp
View file

@ -32,6 +32,7 @@ namespace Cantera {
//====================================================================================================================
SquareMatrix::SquareMatrix() :
DenseMatrix(),
GeneralMatrix(0),
m_factored(0),
a1norm_(0.0),
useQR_(0)
@ -48,7 +49,8 @@ namespace Cantera {
* @param v intial value of all matrix components.
*/
SquareMatrix::SquareMatrix(int n, doublereal v) :
DenseMatrix(n, n, v),
DenseMatrix(n, n, v),
GeneralMatrix(0),
m_factored(0),
a1norm_(0.0),
useQR_(0)
@ -61,7 +63,8 @@ namespace Cantera {
* copy constructor
*/
SquareMatrix::SquareMatrix(const SquareMatrix& y) :
DenseMatrix(y),
DenseMatrix(y),
GeneralMatrix(0),
m_factored(y.m_factored),
a1norm_(y.a1norm_),
useQR_(y.useQR_)
@ -75,6 +78,7 @@ namespace Cantera {
SquareMatrix& SquareMatrix::operator=(const SquareMatrix& y) {
if (&y == this) return *this;
DenseMatrix::operator=(y);
GeneralMatrix::operator=(y);
m_factored = y.m_factored;
a1norm_ = y.a1norm_;
useQR_ = y.useQR_;
@ -87,7 +91,7 @@ namespace Cantera {
/*
* Solve Ax = b. Vector b is overwritten on exit with x.
*/
int SquareMatrix::solve(double* b)
int SquareMatrix::solve(doublereal * b)
{
if (useQR_) {
return solveQR(b);
@ -138,6 +142,25 @@ namespace Cantera {
void SquareMatrix::resize(int n, int m, doublereal v) {
DenseMatrix::resize(n, m, v);
}
//====================================================================================================================
// Multiply A*b and write result to prod.
/*
* @param b Vector to do the rh multiplcation
* @param prod OUTPUT vector to receive the result
*/
void SquareMatrix::mult(const doublereal * const b, doublereal * const prod) const {
DenseMatrix::mult(b, prod);
}
//====================================================================================================================
// Multiply b*A and write result to prod.
/*
* @param b Vector to do the lh multiplcation
* @param prod OUTPUT vector to receive the result
*/
void SquareMatrix::leftMult(const doublereal * const b, doublereal * const prod) const {
DenseMatrix::leftMult(b, prod);
}
//====================================================================================================================
/*
* Factor A. A is overwritten with the LU decomposition of A.
@ -206,7 +229,7 @@ namespace Cantera {
/*
* Solve Ax = b. Vector b is overwritten on exit with x.
*/
int SquareMatrix::solveQR(double* b)
int SquareMatrix::solveQR(doublereal * b)
{
int info=0;
/*
@ -288,7 +311,11 @@ namespace Cantera {
}
return rcond;
}
//=====================================================================================================================
//=====================================================================================================================
doublereal SquareMatrix::oneNorm() const {
return a1norm_;
}
//=====================================================================================================================
doublereal SquareMatrix::rcondQR() {
if ((int) iwork_.size() < m_nrows) {
@ -316,6 +343,111 @@ namespace Cantera {
return rcond;
}
//=====================================================================================================================
void SquareMatrix::useFactorAlgorithm(int fAlgorithm) {
useQR_ = fAlgorithm;
}
//=====================================================================================================================
int SquareMatrix::factorAlgorithm() const {
return (int) useQR_;
}
//=====================================================================================================================
bool SquareMatrix::factored() const {
return m_factored;
}
//=====================================================================================================================
// Return a pointer to the top of column j, columns are contiguous in memory
/*
* @param j Value of the column
*
* @return Returns a pointer to the top of the column
*/
doublereal * SquareMatrix::ptrColumn(int j) {
return Array2D::ptrColumn(j);
}
//=====================================================================================================================
// Copy the data from one array into another without doing any checking
/*
* This differs from the assignment operator as no resizing is done and memcpy() is used.
* @param y Array to be copied
*/
void SquareMatrix::copyData(const GeneralMatrix& y) {
const SquareMatrix *yy_ptr = dynamic_cast<const SquareMatrix *>(& y);
Array2D::copyData(*yy_ptr);
}
//=====================================================================================================================
size_t SquareMatrix::nRows() const {
return m_nrows;
}
//=====================================================================================================================
size_t SquareMatrix::nRowsAndStruct(int * const iStruct) const {
return m_nrows;
}
//=====================================================================================================================
GeneralMatrix * SquareMatrix::duplMyselfAsGeneralMatrix() const {
SquareMatrix *dd = new SquareMatrix(*this);
return static_cast<GeneralMatrix *>(dd);
}
//=====================================================================================================================
// Return an iterator pointing to the first element
vector_fp::iterator SquareMatrix::begin() {
return m_data.begin();
}
//=====================================================================================================================
// Return a const iterator pointing to the first element
vector_fp::const_iterator SquareMatrix::begin() const {
return m_data.begin();
}
//=====================================================================================================================
// Return a vector of const pointers to the columns
/*
* Note the value of the pointers are protected by their being const.
* However, the value of the matrix is open to being changed.
*
* @return returns a vector of pointers to the top of the columns
* of the matrices.
*/
doublereal * const * SquareMatrix::colPts() {
return DenseMatrix::colPts();
}
//=====================================================================================================================
int SquareMatrix::checkRows(doublereal &valueSmall) const {
valueSmall = 1.0E300;
int iSmall = -1;
for (int i = 0; i < m_nrows; i++) {
double valueS = 0.0;
for (int j = 0; j < m_nrows; j++) {
if (fabs(value(i,j)) > valueS) {
valueS = fabs(value(i,j));
}
}
if (valueS < valueSmall) {
iSmall = i;
valueSmall = valueS;
}
}
return iSmall;
}
//=====================================================================================================================
int SquareMatrix::checkColumns(doublereal &valueSmall) const {
valueSmall = 1.0E300;
int jSmall = -1;
for (int j = 0; j < m_nrows; j++) {
double valueS = 0.0;
for (int i = 0; i < m_nrows; i++) {
if (fabs(value(i,j)) > valueS) {
valueS = fabs(value(i,j));
}
}
if (valueS < valueSmall) {
jSmall = j;
valueSmall = valueS;
}
}
return jSmall;
}
//=====================================================================================================================
}

168
Cantera/src/numerics/SquareMatrix.h
View file

@ -18,6 +18,7 @@
#define CT_SQUAREMATRIX_H
#include "DenseMatrix.h"
#include "GeneralMatrix.h"
namespace Cantera {
@ -25,7 +26,7 @@ namespace Cantera {
* A class for full (non-sparse) matrices with Fortran-compatible
* data storage. Adds matrix inversion operations to this class from DenseMatrix.
*/
class SquareMatrix: public DenseMatrix {
class SquareMatrix: public DenseMatrix, public GeneralMatrix {
public:
@ -61,8 +62,7 @@ namespace Cantera {
//! Destructor. Does nothing.
virtual ~SquareMatrix();
//! Solves the Ax = b system returning x in the b spot.
//! Solves the Ax = b system returning x in the b spot.
/*!
* @param b Vector for the rhs of the equation system
*/
@ -76,12 +76,25 @@ namespace Cantera {
*/
void resize(int n, int m, doublereal v = 0.0);
/**
* Zero the matrix
*/
void zero();
//! Multiply A*b and write result to prod.
/*!
* @param b Vector to do the rh multiplcation
* @param prod OUTPUT vector to receive the result
*/
virtual void mult(const doublereal * const b, doublereal * const prod) const;
//! Multiply b*A and write result to prod.
/*!
* @param b Vector to do the lh multiplcation
* @param prod OUTPUT vector to receive the result
*/
virtual void leftMult(const doublereal * const b, doublereal * const prod) const;
/**
* Factors the A matrix, overwriting A. We flip m_factored
* boolean to indicate that the matrix is now A-1.
@ -94,7 +107,7 @@ namespace Cantera {
*
* @return Returns the info variable from lapack
*/
int factorQR();
virtual int factorQR();
//! Returns an estimate of the inverse of the condition number for the matrix
/*!
@ -102,7 +115,7 @@ namespace Cantera {
*
* @return returns the inverse of the condition number
*/
doublereal rcondQR();
virtual doublereal rcondQR();
//! Returns an estimate of the inverse of the condition number for the matrix
/*!
@ -112,7 +125,10 @@ namespace Cantera {
*
* @return returns the inverse of the condition number
*/
doublereal rcond(doublereal a1norm);
virtual doublereal rcond(doublereal a1norm);
//! Returns the one norm of the matrix
virtual doublereal oneNorm() const;
//! Solves the linear problem Ax=b using the QR algorithm returning x in the b spot
/*!
@ -122,17 +138,136 @@ namespace Cantera {
//! clear the factored flag
void clearFactorFlag();
/**
* set the factored flag
*/
virtual void clearFactorFlag();
//! set the factored flag
void setFactorFlag();
/*
* the factor flag
//! Report whether the current matrix has been factored.
virtual bool factored() const;
//! Change the way the matrix is factored
/*!
* @param fAlgorithm integer
* 0 LU factorization
* 1 QR factorization
*/
virtual void useFactorAlgorithm(int fAlgorithm);
//! Returns the factor algorithm used
/*!
* 0 LU decomposition
* 1 QR decomposition
*
* This routine will always return 0
*/
virtual int factorAlgorithm() const;
//! Return a pointer to the top of column j, columns are assumed to be contiguous in memory
/*!
* @param j Value of the column
*
* @return Returns a pointer to the top of the column
*/
virtual doublereal * ptrColumn(int j);
//! Index into the (i,j) element
/*!
* @param i row
* @param j column
*
* (note, tried a using directive here, and it didn't seem to work)
*
* Returns a changeable reference to the matrix entry
*/
virtual doublereal& operator()(int i, int j) {
return Array2D::operator()(i, j);
}
//! Copy the data from one array into another without doing any checking
/*!
* This differs from the assignment operator as no resizing is done and memcpy() is used.
* @param y Array to be copied
*/
virtual void copyData(const GeneralMatrix& y);
//! Constant Index into the (i,j) element
/*!
* @param i row
* @param j column
*
* Returns an unchangeable reference to the matrix entry
*/
virtual doublereal operator() (int i, int j) const {
return Array2D::operator()(i, j);
}
//! Return the number of rows in the matrix
virtual size_t nRows() const;
//! Return the size and structure of the matrix
/*!
* This is inherited from GeneralMatrix
*
* @param iStruct OUTPUT Pointer to a vector of ints that describe the structure of the matrix.
* not used
*
* @return returns the number of rows and columns in the matrix.
*/
size_t nRowsAndStruct(int * const iStruct = 0) const;
//! Duplicate this object
virtual GeneralMatrix * duplMyselfAsGeneralMatrix() const;
//! Return an iterator pointing to the first element
/*!
*/
virtual vector_fp::iterator begin();
//! Return a const iterator pointing to the first element
virtual vector_fp::const_iterator begin() const;
//! Return a vector of const pointers to the columns
/*!
* Note the value of the pointers are protected by their being const.
* However, the value of the matrix is open to being changed.
*
* @return returns a vector of pointers to the top of the columns
* of the matrices.
*/
virtual doublereal * const * colPts();
//! Check to see if we have any zero rows in the jacobian
/*!
* This utility routine checks to see if any rows are zero.
* The smallest row is returned along with the largest coefficient in that row
*
* @param valueSmall OUTPUT value of the largest coefficient in the smallest row
*
* @return index of the row that is most nearly zero
*/
virtual int checkRows(doublereal & valueSmall) const;
//! Check to see if we have any zero columns in the jacobian
/*!
* This utility routine checks to see if any columns are zero.
* The smallest column is returned along with the largest coefficient in that column
*
* @param valueSmall OUTPUT value of the largest coefficient in the smallest column
*
* @return index of the column that is most nearly zero
*/
virtual int checkColumns(doublereal & valueSmall) const;
protected:
//! the factor flag
int m_factored;
public:
//! Work vector for QR algorithm
vector_fp tau;
@ -141,11 +276,10 @@ namespace Cantera {
//! Integer work vector for QR algorithms
std::vector<int> iwork_;
protected:
//! 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
int useQR_;
};
@ -153,5 +287,3 @@ namespace Cantera {
#endif

43
Cantera/src/numerics/ctlapack.h
View file

@ -29,6 +29,7 @@
#define _DGETRS_ dgetrs
#define _DGETRI_ dgetri
#define _DGELSS_ dgelss
#define _DGBCON_ dgbcon
#define _DGBSV_ dgbsv
#define _DGBTRF_ dgbtrf
#define _DGBTRS_ dgbtrs
@ -51,6 +52,7 @@
#define _DGETRS_ dgetrs_
#define _DGETRI_ dgetri_
#define _DGELSS_ dgelss_
#define _DGBCON_ dgbcon_
#define _DGBSV_ dgbsv_
#define _DGBTRF_ dgbtrf_
#define _DGBTRS_ dgbtrs_
@ -222,6 +224,16 @@ extern "C" {
#endif
#ifdef LAPACK_FTN_STRING_LEN_AT_END
int _DGBCON_(const char *norm, const integer* n, integer *kl, integer *ku, doublereal* ab, const integer* ldab,
const integer *ipiv, const doublereal *anorm, const doublereal *rcond,
doublereal* work, const integer* iwork, integer *info, ftnlen nosize);
#else
int _DGBCON_(const char *norm, ftnlen nosize, const integer* n, integer *kl, integer *ku, doublereal* ab, const integer* ldab,
const integer *ipiv, const doublereal *anorm, 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);
@ -535,6 +547,37 @@ namespace Cantera {
#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_dgbcon(const char norm, int n, int kl, int ku, doublereal* a, int ldab, int *ipiv, doublereal anorm,
doublereal* work, int *iwork, int &info) {
char cnorm = '1';
if (norm) {
cnorm = norm;
}
integer f_n = n;
integer f_kl = kl;
integer f_ku = ku;
integer f_ldab = ldab;
integer f_info = info;
doublereal rcond;
#ifdef NO_FTN_STRING_LEN_AT_END
_DGBCON_(&cnorm, &f_n , &f_kl, &f_ku, a, &f_ldab, ipiv, &anorm, &rcond, work, iwork, &f_info);
#else
ftnlen trsize = 1;
#ifdef LAPACK_FTN_STRING_LEN_AT_END
_DGBCON_(&cnorm, &f_n, &f_kl, &f_ku, a, &f_ldab, ipiv, &anorm, &rcond, work, iwork, &f_info, trsize);
#else
_DGBCON_(&cnorm, trsize, &f_n, &f_kl, &f_ku, a, &f_ldab, ipiv, &anorm, &rcond, work, iwork, &f_info);
#endif
#endif
info = f_info;
return rcond;

2
Cantera/src/oneD/MultiNewton.cpp
View file

@ -133,7 +133,7 @@ namespace Cantera {
}
#endif
iok = jac.solve(sz, step, step);
iok = jac.solve(step, step);
// if iok is non-zero, then solve failed
if (iok > 0) {

10
Cantera/src/thermo/Makefile.in
View file

@ -79,11 +79,13 @@ ELECTRO_H = MolalityVPSSTP.h VPStandardStateTP.h \
endif
ifeq ($(do_issp),1)
ISSP_OBJ = IdealSolidSolnPhase.o StoichSubstanceSSTP.o SingleSpeciesTP.o MineralEQ3.o \
GibbsExcessVPSSTP.o PseudoBinaryVPSSTP.o MargulesVPSSTP.o \
IonsFromNeutralVPSSTP.o PDSS_IonsFromNeutral.o FixedChemPotSSTP.o
GibbsExcessVPSSTP.o MolarityIonicVPSSTP.o MargulesVPSSTP.o \
IonsFromNeutralVPSSTP.o PDSS_IonsFromNeutral.o FixedChemPotSSTP.o \
MixedSolventElectrolyte.o
ISSP_H = IdealSolidSolnPhase.h StoichSubstanceSSTP.h SingleSpeciesTP.h MineralEQ3.h \
GibbsExcessVPSSTP.h PseudoBinaryVPSSTP.h MargulesVPSSTP.h \
IonsFromNeutralVPSSTP.h PDSS_IonsFromNeutral.h FixedChemPotSSTP.h
GibbsExcessVPSSTP.h MolarityIonicVPSSTP.h MargulesVPSSTP.h \
IonsFromNeutralVPSSTP.h PDSS_IonsFromNeutral.h FixedChemPotSSTP.h \
MixedSolventElectrolyte.h
endif
CATHERMO_OBJ = $(THERMO_OBJ) $(ELECTRO_OBJ) $(ISSP_OBJ)

2
Cantera/src/thermo/MargulesVPSSTP.h
View file

@ -23,7 +23,7 @@
#ifndef CT_MARGULESVPSSTP_H
#define CT_MARGULESVPSSTP_H
#include "PseudoBinaryVPSSTP.h"
#include "GibbsExcessVPSSTP.h"
namespace Cantera {

16
Cantera/src/thermo/MixedSolventElectrolyte.cpp
View file

@ -31,7 +31,7 @@ namespace Cantera {
*
*/
MixedSolventElectrolyte::MixedSolventElectrolyte() :
GibbsExcessVPSSTP(),
MolarityIonicVPSSTP(),
numBinaryInteractions_(0),
formMargules_(0),
formTempModel_(0)
@ -48,7 +48,7 @@ namespace Cantera {
*/
MixedSolventElectrolyte::MixedSolventElectrolyte(std::string inputFile, std::string id) :
GibbsExcessVPSSTP(),
MolarityIonicVPSSTP(),
numBinaryInteractions_(0),
formMargules_(0),
formTempModel_(0)
@ -57,7 +57,7 @@ namespace Cantera {
}
MixedSolventElectrolyte::MixedSolventElectrolyte(XML_Node& phaseRoot, std::string id) :
GibbsExcessVPSSTP(),
MolarityIonicVPSSTP(),
numBinaryInteractions_(0),
formMargules_(0),
formTempModel_(0)
@ -73,7 +73,7 @@ namespace Cantera {
* has a working copy constructor
*/
MixedSolventElectrolyte::MixedSolventElectrolyte(const MixedSolventElectrolyte &b) :
GibbsExcessVPSSTP()
MolarityIonicVPSSTP()
{
MixedSolventElectrolyte::operator=(b);
}
@ -90,7 +90,7 @@ namespace Cantera {
return *this;
}
GibbsExcessVPSSTP::operator=(b);
MolarityIonicVPSSTP::operator=(b);
numBinaryInteractions_ = b.numBinaryInteractions_ ;
m_HE_b_ij = b.m_HE_b_ij;
@ -141,7 +141,7 @@ namespace Cantera {
*
*/
MixedSolventElectrolyte::MixedSolventElectrolyte(int testProb) :
GibbsExcessVPSSTP(),
MolarityIonicVPSSTP(),
numBinaryInteractions_(0),
formMargules_(0),
formTempModel_(0)
@ -659,7 +659,7 @@ namespace Cantera {
*/
void MixedSolventElectrolyte::initThermo() {
initLengths();
GibbsExcessVPSSTP::initThermo();
MolarityIonicVPSSTP::initThermo();
}
@ -741,7 +741,7 @@ namespace Cantera {
/*
* Go down the chain
*/
GibbsExcessVPSSTP::initThermoXML(phaseNode, id);
MolarityIonicVPSSTP::initThermoXML(phaseNode, id);
}

4
Cantera/src/thermo/MixedSolventElectrolyte.h
View file

@ -23,7 +23,7 @@
#ifndef CT_MIXEDSOLVENTELECTROLYTEVPSSTP_H
#define CT_MIXEDSOLVENTELECTROLYTEVPSSTP_H
#include "PseudoBinaryVPSSTP.h"
#include "MolarityIonicVPSSTP.h"
namespace Cantera {
@ -311,7 +311,7 @@ namespace Cantera {
*
*/
class MixedSolventElectrolyte : public GibbsExcessVPSSTP {
class MixedSolventElectrolyte : public MolarityIonicVPSSTP {
public:

179
Cantera/src/thermo/PseudoBinaryVPSSTP.cpp → Cantera/src/thermo/MolarityIonicVPSSTP.cpp
View file

@ -1,9 +1,9 @@
/**
* @file PseudoBinaryVPSSTP.cpp
* @file MolarityIonicVPSSTP.cpp
* Definitions for intermediate ThermoPhase object for phases which
* employ excess gibbs free energy formulations
* (see \ref thermoprops
* and class \link Cantera::PseudoBinaryVPSSTP PseudoBinaryVPSSTP\endlink).
* and class \link Cantera::MolarityIonicVPSSTP MolarityIonicVPSSTP\endlink).
*
* Header file for a derived class of ThermoPhase that handles
* variable pressure standard state methods for calculating
@ -17,24 +17,24 @@
* U.S. Government retains certain rights in this software.
*/
/*
* $Date$
* $Revision$
* $Date: 2009-11-09 16:36:49 -0700 (Mon, 09 Nov 2009) $
* $Revision: 255 $
*/
#include "PseudoBinaryVPSSTP.h"
#include "MolarityIonicVPSSTP.h"
#include <cmath>
using namespace std;
namespace Cantera {
//====================================================================================================================
/*
* Default constructor.
*
*/
PseudoBinaryVPSSTP::PseudoBinaryVPSSTP() :
MolarityIonicVPSSTP::MolarityIonicVPSSTP() :
GibbsExcessVPSSTP(),
PBType_(PBTYPE_PASSTHROUGH),
numPBSpecies_(m_kk),
@ -42,19 +42,17 @@ namespace Cantera {
numCationSpecies_(0),
numAnionSpecies_(0),
numPassThroughSpecies_(0),
neutralPBindexStart(0),
cationPhase_(0),
anionPhase_(0)
neutralPBindexStart(0)
{
}
//====================================================================================================================
/*
* Copy Constructor:
*
* Note this stuff will not work until the underlying phase
* has a working copy constructor
*/
PseudoBinaryVPSSTP::PseudoBinaryVPSSTP(const PseudoBinaryVPSSTP &b) :
MolarityIonicVPSSTP::MolarityIonicVPSSTP(const MolarityIonicVPSSTP &b) :
GibbsExcessVPSSTP(),
PBType_(PBTYPE_PASSTHROUGH),
numPBSpecies_(m_kk),
@ -62,21 +60,19 @@ namespace Cantera {
numCationSpecies_(0),
numAnionSpecies_(0),
numPassThroughSpecies_(0),
neutralPBindexStart(0),
cationPhase_(0),
anionPhase_(0)
neutralPBindexStart(0)
{
*this = operator=(b);
}
//====================================================================================================================
/*
* operator=()
*
* Note this stuff will not work until the underlying phase
* has a working assignment operator
*/
PseudoBinaryVPSSTP& PseudoBinaryVPSSTP::
operator=(const PseudoBinaryVPSSTP &b) {
MolarityIonicVPSSTP& MolarityIonicVPSSTP::
operator=(const MolarityIonicVPSSTP &b) {
if (&b != this) {
GibbsExcessVPSSTP::operator=(b);
}
@ -92,21 +88,19 @@ namespace Cantera {
passThroughList_ = b.passThroughList_;
numPassThroughSpecies_ = b.numPassThroughSpecies_;
neutralPBindexStart = b.neutralPBindexStart;
cationPhase_ = b.cationPhase_;
anionPhase_ = b.anionPhase_;
moleFractionsTmp_ = b.moleFractionsTmp_;
return *this;
}
//====================================================================================================================
/**
*
* ~PseudoBinaryVPSSTP(): (virtual)
* ~MolarityIonicVPSSTP(): (virtual)
*
* Destructor: does nothing:
*
*/
PseudoBinaryVPSSTP::~PseudoBinaryVPSSTP() {
MolarityIonicVPSSTP::~MolarityIonicVPSSTP() {
}
/*
@ -114,25 +108,25 @@ namespace Cantera {
* a pointer to ThermoPhase.
*/
ThermoPhase*
PseudoBinaryVPSSTP::duplMyselfAsThermoPhase() const {
PseudoBinaryVPSSTP* mtp = new PseudoBinaryVPSSTP(*this);
MolarityIonicVPSSTP::duplMyselfAsThermoPhase() const {
MolarityIonicVPSSTP* mtp = new MolarityIonicVPSSTP(*this);
return (ThermoPhase *) mtp;
}
/*
* -------------- Utilities -------------------------------
*/
//====================================================================================================================
// Equation of state type flag.
/*
* The ThermoPhase base class returns
* zero. Subclasses should define this to return a unique
* non-zero value. Known constants defined for this purpose are
* listed in mix_defs.h. The PseudoBinaryVPSSTP class also returns
* listed in mix_defs.h. The MolarityIonicVPSSTP class also returns
* zero, as it is a non-complete class.
*/
int PseudoBinaryVPSSTP::eosType() const {
int MolarityIonicVPSSTP::eosType() const {
return 0;
}
@ -147,32 +141,35 @@ namespace Cantera {
* - Activities, Standard States, Activity Concentrations -----------
*/
doublereal PseudoBinaryVPSSTP::standardConcentration(int k) const {
//====================================================================================================================
doublereal MolarityIonicVPSSTP::standardConcentration(int k) const {
err("standardConcentration");
return -1.0;
}
doublereal PseudoBinaryVPSSTP::logStandardConc(int k) const {
//====================================================================================================================
doublereal MolarityIonicVPSSTP::logStandardConc(int k) const {
err("logStandardConc");
return -1.0;
}
//====================================================================================================================
void PseudoBinaryVPSSTP::getElectrochemPotentials(doublereal* mu) const {
void MolarityIonicVPSSTP::getElectrochemPotentials(doublereal* mu) const {
getChemPotentials(mu);
double ve = Faraday * electricPotential();
for (int k = 0; k < m_kk; k++) {
mu[k] += ve*charge(k);
}
}
void PseudoBinaryVPSSTP::calcPseudoBinaryMoleFractions() const {
//====================================================================================================================
void MolarityIonicVPSSTP::calcPseudoBinaryMoleFractions() const {
int k;
int kCat;
int kMax;
doublereal sumCat;
doublereal sumAnion;
doublereal chP, chM;
doublereal sum = 0.0;
doublereal sumMax;
switch (PBType_) {
case PBTYPE_PASSTHROUGH:
for (k = 0; k < m_kk; k++) {
@ -185,26 +182,43 @@ namespace Cantera {
for (k = 0; k < m_kk; k++) {
moleFractionsTmp_[k] = moleFractions_[k];
}
kMax = -1;
sumMax = 0.0;
for (k = 0; k < (int) cationList_.size(); k++) {
sumCat += moleFractions_[cationList_[k]];
kCat = cationList_[k];
chP = m_speciesCharge[kCat];
if (moleFractions_[kCat] > sumMax) {
kMax = k;
sumMax = moleFractions_[kCat];
}
sumCat += chP * moleFractions_[kCat];
}
k = anionList_[0];
chM = m_speciesCharge[k];
sumAnion = moleFractions_[k] * chM;
sum = sumCat - sumAnion;
if (fabs(sum) > 1.0E-16) {
moleFractionsTmp_[cationList_[kMax]] -= sum / m_speciesCharge[kMax];
sum = 0.0;
for (k = 0; k < numCationSpecies_; k++) {
sum += moleFractionsTmp_[k];
}
for (k = 0; k < numCationSpecies_; k++) {
moleFractionsTmp_[k]/= sum;
}
}
sumAnion = moleFractions_[anionList_[k]];
PBMoleFractions_[0] = sumCat -sumAnion;
moleFractionsTmp_[indexSpecialSpecies_] -= PBMoleFractions_[0];
for (k = 0; k < numCationSpecies_; k++) {
PBMoleFractions_[1+k] = moleFractionsTmp_[cationList_[k]];
PBMoleFractions_[k] = moleFractionsTmp_[cationList_[k]];
}
for (k = 0; k < numPassThroughSpecies_; k++) {
PBMoleFractions_[neutralPBindexStart + k] =
moleFractions_[cationList_[k]];
PBMoleFractions_[neutralPBindexStart + k] = moleFractions_[passThroughList_[k]];
}
sum = fmaxx(0.0, PBMoleFractions_[0]);
for (k = 1; k < numPBSpecies_; k++) {
sum += PBMoleFractions_[k];
}
for (k = 0; k < numPBSpecies_; k++) {
PBMoleFractions_[k] /= sum;
@ -226,19 +240,17 @@ namespace Cantera {
}
}
//====================================================================================================================
/*
* ------------ Partial Molar Properties of the Solution ------------
*/
doublereal PseudoBinaryVPSSTP::err(std::string msg) const {
throw CanteraError("PseudoBinaryVPSSTP","Base class method "
//====================================================================================================================
doublereal MolarityIonicVPSSTP::err(std::string msg) const {
throw CanteraError("MolarityIonicVPSSTP","Base class method "
+msg+" called. Equation of state type: "+int2str(eosType()));
return 0;
}
//====================================================================================================================
/*
* @internal Initialize. This method is provided to allow
* subclasses to perform any initialization required after all
@ -252,19 +264,49 @@ namespace Cantera {
*
* @see importCTML.cpp
*/
void PseudoBinaryVPSSTP::initThermo() {
initLengths();
void MolarityIonicVPSSTP::initThermo() {
GibbsExcessVPSSTP::initThermo();
initLengths();
/*
* Go find the list of cations and anions
*/
double ch;
numCationSpecies_ = 0.0;
cationList_.clear();
anionList_.clear();
passThroughList_.clear();
for (int k = 0; k < m_kk; k++) {
ch = m_speciesCharge[k];
if (ch > 0.0) {
cationList_.push_back(k);
numCationSpecies_++;
} else if (ch < 0.0) {
anionList_.push_back(k);
numAnionSpecies_++;
} else {
passThroughList_.push_back(k);
numPassThroughSpecies_++;
}
}
numPBSpecies_ = numCationSpecies_ + numAnionSpecies_ - 1;
neutralPBindexStart = numPBSpecies_;
PBType_ = PBTYPE_MULTICATIONANION;
if (numAnionSpecies_ == 1) {
PBType_ = PBTYPE_SINGLEANION;
} else if (numCationSpecies_ == 1) {
PBType_ = PBTYPE_SINGLECATION;
}
if (numAnionSpecies_ == 0 && numCationSpecies_ == 0) {
PBType_ = PBTYPE_PASSTHROUGH;
}
}
// Initialize lengths of local variables after all species have
// been identified.
void PseudoBinaryVPSSTP::initLengths() {
//====================================================================================================================
// Initialize lengths of local variables after all species have been identified.
void MolarityIonicVPSSTP::initLengths() {
m_kk = nSpecies();
moleFractions_.resize(m_kk);
moleFractionsTmp_.resize(m_kk);
}
//====================================================================================================================
/*
* initThermoXML() (virtual from ThermoPhase)
* Import and initialize a ThermoPhase object
@ -280,18 +322,16 @@ namespace Cantera {
* to see if phaseNode is pointing to the phase
* with the correct id.
*/
void PseudoBinaryVPSSTP::initThermoXML(XML_Node& phaseNode, std::string id) {
void MolarityIonicVPSSTP::initThermoXML(XML_Node& phaseNode, std::string id) {
GibbsExcessVPSSTP::initThermoXML(phaseNode, id);
}
/**
//====================================================================================================================
/*
* Format a summary of the mixture state for output.
*/
std::string PseudoBinaryVPSSTP::report(bool show_thermo) const {
std::string MolarityIonicVPSSTP::report(bool show_thermo) const {
char p[800];
string s = "";
try {
@ -364,7 +404,6 @@ namespace Cantera {
}
return s;
}
//====================================================================================================================
}

66
Cantera/src/thermo/PseudoBinaryVPSSTP.h → Cantera/src/thermo/MolarityIonicVPSSTP.h
View file

@ -1,15 +1,15 @@
/**
* @file PseudoBinaryVPSSTP.h
* @file MolarityIonicVPSSTP.h
* Header for intermediate ThermoPhase object for phases which
* employ gibbs excess free energy based formulations
* (see \ref thermoprops
* and class \link Cantera::gibbsExcessVPSSTP gibbsExcessVPSSTP\endlink).
* and class \link Cantera::MolarityIonicVPSSTP MolarityIonicVPSSTP\endlink).
*
* Header file for a derived class of ThermoPhase that handles
* variable pressure standard state methods for calculating
* thermodynamic properties that are further based upon activities
* based on the molality scale. These include most of the methods for
* calculating liquid electrolyte thermodynamics.
* based on the molarity scale. In this class, we expect that there are
* ions, but they are treated on the molarity scale.
*/
/*
* Copywrite (2006) Sandia Corporation. Under the terms of
@ -17,11 +17,11 @@
* U.S. Government retains certain rights in this software.
*/
/*
* $Id$
* $Id: MolarityIonicVPSSTP.h 255 2009-11-09 23:36:49Z hkmoffa $
*/
#ifndef CT_PSEUDOBINARYVPSSTP_H
#define CT_PSEUDOBINARYVPSSTP_H
#ifndef CT_MOLARITYIONICVPSSTP_H
#define CT_MOLARITYIONICVPSSTP_H
#include "GibbsExcessVPSSTP.h"
@ -32,22 +32,14 @@ namespace Cantera {
*/
/*!
* PseudoBinaryVPSSTP is a derived class of ThermoPhase
* MolarityIonicVPSSTP is a derived class of ThermoPhase
* GibbsExcessVPSSTP that handles
* variable pressure standard state methods for calculating
* thermodynamic properties that are further based on
* expressing the Excess Gibbs free energy as a function of
* the mole fractions (or pseudo mole fractions) of consitituents.
* This category is the workhorse for describing molten salts,
* solid-phase mixtures of semiconductors, and mixtures of miscible
* and semi-miscible compounds.
*
* It includes
* . regular solutions
* . Margueles expansions
* . NTRL equation
* . Wilson's equation
* . UNIQUAC equation of state.
* the mole fractions (or pseudo mole fractions) of the consitituents.
* This category is the workhorse for describing ionic systems which
* are not on the molality scale.
*
* This class adds additional functions onto the %ThermoPhase interface
* that handles the calculation of the excess Gibbs free energy. The %ThermoPhase
@ -60,15 +52,14 @@ namespace Cantera {
* symmetrical formulations.
*
* This layer will massage the mole fraction vector to implement
* cation and anion based mole numbers in an optional manner
*
* The way that it collects the cation and anion based mole numbers
* is via holding two extra ThermoPhase objects. These
* can include standard states for salts.
*
* cation and anion based mole numbers in an optional manner, such that
* it is expected that there exists a charge balance at all times.
* One of the ions must be a "special ion" in the sense that its' thermodynamic
* functions are set to zero, and the thermo functions of all other
* ions are based on a valuation relative to the special ion.
*
*/
class PseudoBinaryVPSSTP : public GibbsExcessVPSSTP {
class MolarityIonicVPSSTP : public GibbsExcessVPSSTP {
public:
@ -81,7 +72,7 @@ namespace Cantera {
* density conservation and therefore element conservation
* is the more important principle to follow.
*/
PseudoBinaryVPSSTP();
MolarityIonicVPSSTP();
//! Copy constructor
/*!
@ -90,17 +81,17 @@ namespace Cantera {
*
* @param b class to be copied
*/
PseudoBinaryVPSSTP(const PseudoBinaryVPSSTP&b);
MolarityIonicVPSSTP(const MolarityIonicVPSSTP&b);
/// Assignment operator
/*!
*
* @param b class to be copied.
*/
PseudoBinaryVPSSTP& operator=(const PseudoBinaryVPSSTP&b);
MolarityIonicVPSSTP& operator=(const MolarityIonicVPSSTP&b);
/// Destructor.
virtual ~PseudoBinaryVPSSTP();
virtual ~MolarityIonicVPSSTP();
//! Duplication routine for objects which inherit from ThermoPhase.
/*!
@ -344,6 +335,13 @@ namespace Cantera {
protected:
// Pseudobinary type
/*!
* PBTYPE_PASSTHROUGH All species are passthrough species
* PBTYPE_SINGLEANION there is only one anion in the mixture
* PBTYPE_SINGLECATION there is only one cation in the mixture
* PBTYPE_MULTICATIONANION Complex mixture
*/
int PBType_;
//! Number of pseudo binary species
@ -354,20 +352,20 @@ namespace Cantera {
mutable std::vector<doublereal> PBMoleFractions_;
//! Vector of cation indecises in the mixture
std::vector<int> cationList_;
//! Number of cations in the mixture
int numCationSpecies_;
std::vector<int>anionList_;
std::vector<int> anionList_;
int numAnionSpecies_;
std::vector<int> passThroughList_;
int numPassThroughSpecies_;
int neutralPBindexStart;
ThermoPhase *cationPhase_;
ThermoPhase *anionPhase_;
mutable std::vector<doublereal> moleFractionsTmp_;
private:

2
Cantera/src/thermo/PureFluidPhase.cpp
View file

@ -354,7 +354,7 @@ namespace Cantera {
void PureFluidPhase::getEnthalpy_RT_ref(doublereal *hrt) const {
double psave = pressure();
double t = temperature();
double pref = m_spthermo->refPressure();
//double pref = m_spthermo->refPressure();
double plow = 1.0E-8;
Set(tpx::TP, t, plow);
getEnthalpy_RT(hrt);

10
Cantera/src/thermo/ThermoFactory.cpp
View file

@ -78,6 +78,8 @@
#include "HMWSoln.h"
#include "DebyeHuckel.h"
#include "IdealMolalSoln.h"
#include "MolarityIonicVPSSTP.h"
#include "MixedSolventElectrolyte.h"
#endif
#include "IdealSolnGasVPSS.h"
@ -97,7 +99,7 @@ namespace Cantera {
/*!
* @deprecated This entire structure could be replaced with a std::map
*/
static int ntypes = 20;
static int ntypes = 22;
//! Define the string name of the %ThermoPhase types that are handled by this factory routine
static string _types[] = {"IdealGas", "Incompressible",
@ -106,7 +108,8 @@ namespace Cantera {
"HMW", "IdealSolidSolution", "DebyeHuckel",
"IdealMolalSolution", "IdealGasVPSS",
"MineralEQ3", "MetalSHEelectrons", "Margules", "PhaseCombo_Interaction",
"IonsFromNeutralMolecule", "FixedChemPot"
"IonsFromNeutralMolecule", "FixedChemPot", "MolarityIonicVPSSTP",
"MixedSolventElectrolyte"
};
//! Define the integer id of the %ThermoPhase types that are handled by this factory routine
@ -116,7 +119,8 @@ namespace Cantera {
cHMW, cIdealSolidSolnPhase, cDebyeHuckel,
cIdealMolalSoln, cVPSS_IdealGas,
cMineralEQ3, cMetalSHEelectrons,
cMargulesVPSSTP, cPhaseCombo_Interaction, cIonsFromNeutral, cFixedChemPot
cMargulesVPSSTP, cPhaseCombo_Interaction, cIonsFromNeutral, cFixedChemPot,
cMolarityIonicVPSSTP, cMixedSolventElectrolyte
};
/*

3
Cantera/src/thermo/mix_defs.h
View file

@ -79,6 +79,9 @@ namespace Cantera {
const int cMargulesVPSSTP = 301;
const int cMolarityIonicVPSSTP = 401;
const int cMixedSolventElectrolyte = 402;
const int cPhaseCombo_Interaction = 305;
const int cIonsFromNeutral = 2000;

1548
configure vendored
View file

File diff suppressed because it is too large Load diff

8
docs/Cantera.cfg.in
View file

@ -798,7 +798,13 @@ FILE_PATTERNS = Kinetics.h Kinetics.cpp \
RootFind.h \
RootFind.cpp \
NonlinearSolver.h \
NonlinearSolver.cpp
NonlinearSolver.cpp \
BandMatrix.h \
BandMatrix.cpp \
GeneralMatrix.h \
GeneralMatrix.cpp \
SquareMatrix.h \
SquareMatrix.cpp
# The RECURSIVE tag can be used to turn specify whether or not subdirectories
# should be searched for input files as well. Possible values are YES and NO.

3
docs/install_examples/linux.64_sierra_gcc444_python264_numpy
View file

@ -85,6 +85,9 @@ export CC
F77='/sierra/Sntools/extras/compilers/gcc-4.4.4/bin/gfortran'
export F77
FFLAGS="-g -fno-second-underscore"
export FFLAGS
CFLAGS="-g -Wall"
export CFLAGS

3
ext/f2c_lapack/Makefile.in
View file

@ -111,6 +111,9 @@ dtrtrs.o \
dgecon.o \
dgeequ.o \
dgerfs.o \
dgbcon.o \
dgbequ.o \
dlatbs.o \
ieeeck.o \
ilaenv.o

283
ext/f2c_lapack/dgbcon.c Normal file
View file

@ -0,0 +1,283 @@
/* dgbcon.f -- translated by f2c (version 20031025).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
/* Table of constant values */
static integer c__1 = 1;
/* Subroutine */ int dgbcon_(char *norm, integer *n, integer *kl, integer *ku,
doublereal *ab, integer *ldab, integer *ipiv, doublereal *anorm,
doublereal *rcond, doublereal *work, integer *iwork, integer *info,
ftnlen norm_len)
{
/* System generated locals */
integer ab_dim1, ab_offset, i__1, i__2, i__3;
doublereal d__1;
/* Local variables */
static integer j;
static doublereal t;
static integer kd, lm, jp, ix, kase;
extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
integer *);
static integer kase1;
static doublereal scale;
extern logical lsame_(char *, char *, ftnlen, ftnlen);
extern /* Subroutine */ int drscl_(integer *, doublereal *, doublereal *,
integer *);
static logical lnoti;
extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *,
integer *, doublereal *, integer *);
extern doublereal dlamch_(char *, ftnlen);
extern /* Subroutine */ int dlacon_(integer *, doublereal *, doublereal *,
integer *, doublereal *, integer *);
extern integer idamax_(integer *, doublereal *, integer *);
extern /* Subroutine */ int dlatbs_(char *, char *, char *, char *,
integer *, integer *, doublereal *, integer *, doublereal *,
doublereal *, doublereal *, integer *, ftnlen, ftnlen, ftnlen,
ftnlen), xerbla_(char *, integer *, ftnlen);
static doublereal ainvnm;
static logical onenrm;
static char normin[1];
static doublereal smlnum;
/* -- LAPACK routine (version 3.0) -- */
/* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
/* Courant Institute, Argonne National Lab, and Rice University */
/* September 30, 1994 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* DGBCON estimates the reciprocal of the condition number of a real */
/* general band matrix A, in either the 1-norm or the infinity-norm, */
/* using the LU factorization computed by DGBTRF. */
/* An estimate is obtained for norm(inv(A)), and the reciprocal of the */
/* condition number is computed as */
/* RCOND = 1 / ( norm(A) * norm(inv(A)) ). */
/* Arguments */
/* ========= */
/* NORM (input) CHARACTER*1 */
/* Specifies whether the 1-norm condition number or the */
/* infinity-norm condition number is required: */
/* = '1' or 'O': 1-norm; */
/* = 'I': Infinity-norm. */
/* N (input) INTEGER */
/* The order of the matrix A. N >= 0. */
/* KL (input) INTEGER */
/* The number of subdiagonals within the band of A. KL >= 0. */
/* KU (input) INTEGER */
/* The number of superdiagonals within the band of A. KU >= 0. */
/* AB (input) DOUBLE PRECISION array, dimension (LDAB,N) */
/* Details of the LU factorization of the band matrix A, as */
/* computed by DGBTRF. U is stored as an upper triangular band */
/* matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and */
/* the multipliers used during the factorization are stored in */
/* rows KL+KU+2 to 2*KL+KU+1. */
/* LDAB (input) INTEGER */
/* The leading dimension of the array AB. LDAB >= 2*KL+KU+1. */
/* IPIV (input) INTEGER array, dimension (N) */
/* The pivot indices; for 1 <= i <= N, row i of the matrix was */
/* interchanged with row IPIV(i). */
/* ANORM (input) DOUBLE PRECISION */
/* If NORM = '1' or 'O', the 1-norm of the original matrix A. */
/* If NORM = 'I', the infinity-norm of the original matrix A. */
/* RCOND (output) DOUBLE PRECISION */
/* The reciprocal of the condition number of the matrix A, */
/* computed as RCOND = 1/(norm(A) * norm(inv(A))). */
/* WORK (workspace) DOUBLE PRECISION array, dimension (3*N) */
/* IWORK (workspace) INTEGER array, dimension (N) */
/* INFO (output) INTEGER */
/* = 0: successful exit */
/* < 0: if INFO = -i, the i-th argument had an illegal value */
/* ===================================================================== */
/* .. Parameters .. */
/* .. */
/* .. Local Scalars .. */
/* .. */
/* .. External Functions .. */
/* .. */
/* .. External Subroutines .. */
/* .. */
/* .. Intrinsic Functions .. */
/* .. */
/* .. Executable Statements .. */
/* Test the input parameters. */
/* Parameter adjustments */
ab_dim1 = *ldab;
ab_offset = 1 + ab_dim1;
ab -= ab_offset;
--ipiv;
--work;
--iwork;
/* Function Body */
*info = 0;
onenrm = *(unsigned char *)norm == '1' || lsame_(norm, "O", (ftnlen)1, (
ftnlen)1);
if (! onenrm && ! lsame_(norm, "I", (ftnlen)1, (ftnlen)1)) {
*info = -1;
} else if (*n < 0) {
*info = -2;
} else if (*kl < 0) {
*info = -3;
} else if (*ku < 0) {
*info = -4;
} else if (*ldab < (*kl << 1) + *ku + 1) {
*info = -6;
} else if (*anorm < 0.) {
*info = -8;
}
if (*info != 0) {
i__1 = -(*info);
xerbla_("DGBCON", &i__1, (ftnlen)6);
return 0;
}
/* Quick return if possible */
*rcond = 0.;
if (*n == 0) {
*rcond = 1.;
return 0;
} else if (*anorm == 0.) {
return 0;
}
smlnum = dlamch_("Safe minimum", (ftnlen)12);
/* Estimate the norm of inv(A). */
ainvnm = 0.;
*(unsigned char *)normin = 'N';
if (onenrm) {
kase1 = 1;
} else {
kase1 = 2;
}
kd = *kl + *ku + 1;
lnoti = *kl > 0;
kase = 0;
L10:
dlacon_(n, &work[*n + 1], &work[1], &iwork[1], &ainvnm, &kase);
if (kase != 0) {
if (kase == kase1) {
/* Multiply by inv(L). */
if (lnoti) {
i__1 = *n - 1;
for (j = 1; j <= i__1; ++j) {
/* Computing MIN */
i__2 = *kl, i__3 = *n - j;
lm = min(i__2,i__3);
jp = ipiv[j];
t = work[jp];
if (jp != j) {
work[jp] = work[j];
work[j] = t;
}
d__1 = -t;
daxpy_(&lm, &d__1, &ab[kd + 1 + j * ab_dim1], &c__1, &
work[j + 1], &c__1);
/* L20: */
}
}
/* Multiply by inv(U). */
i__1 = *kl + *ku;
dlatbs_("Upper", "No transpose", "Non-unit", normin, n, &i__1, &
ab[ab_offset], ldab, &work[1], &scale, &work[(*n << 1) +
1], info, (ftnlen)5, (ftnlen)12, (ftnlen)8, (ftnlen)1);
} else {
/* Multiply by inv(U'). */
i__1 = *kl + *ku;
dlatbs_("Upper", "Transpose", "Non-unit", normin, n, &i__1, &ab[
ab_offset], ldab, &work[1], &scale, &work[(*n << 1) + 1],
info, (ftnlen)5, (ftnlen)9, (ftnlen)8, (ftnlen)1);
/* Multiply by inv(L'). */
if (lnoti) {
for (j = *n - 1; j >= 1; --j) {
/* Computing MIN */
i__1 = *kl, i__2 = *n - j;
lm = min(i__1,i__2);
work[j] -= ddot_(&lm, &ab[kd + 1 + j * ab_dim1], &c__1, &
work[j + 1], &c__1);
jp = ipiv[j];
if (jp != j) {
t = work[jp];
work[jp] = work[j];
work[j] = t;
}
/* L30: */
}
}
}
/* Divide X by 1/SCALE if doing so will not cause overflow. */
*(unsigned char *)normin = 'Y';
if (scale != 1.) {
ix = idamax_(n, &work[1], &c__1);
if (scale < (d__1 = work[ix], abs(d__1)) * smlnum || scale == 0.)
{
goto L40;
}
drscl_(n, &scale, &work[1], &c__1);
}
goto L10;
}
/* Compute the estimate of the reciprocal condition number. */
if (ainvnm != 0.) {
*rcond = 1. / ainvnm / *anorm;
}
L40:
return 0;
/* End of DGBCON */
} /* dgbcon_ */

321
ext/f2c_lapack/dgbequ.c Normal file
View file

@ -0,0 +1,321 @@
/* dgbequ.f -- translated by f2c (version 20031025).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
/* Subroutine */ int dgbequ_(integer *m, integer *n, integer *kl, integer *ku,
doublereal *ab, integer *ldab, doublereal *r__, doublereal *c__,
doublereal *rowcnd, doublereal *colcnd, doublereal *amax, integer *
info)
{
/* System generated locals */
integer ab_dim1, ab_offset, i__1, i__2, i__3, i__4;
doublereal d__1, d__2, d__3;
/* Local variables */
static integer i__, j, kd;
static doublereal rcmin, rcmax;
extern doublereal dlamch_(char *, ftnlen);
extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
static doublereal bignum, smlnum;
/* -- LAPACK routine (version 3.0) -- */
/* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
/* Courant Institute, Argonne National Lab, and Rice University */
/* March 31, 1993 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* DGBEQU computes row and column scalings intended to equilibrate an */
/* M-by-N band matrix A and reduce its condition number. R returns the */
/* row scale factors and C the column scale factors, chosen to try to */
/* make the largest element in each row and column of the matrix B with */
/* elements B(i,j)=R(i)*A(i,j)*C(j) have absolute value 1. */
/* R(i) and C(j) are restricted to be between SMLNUM = smallest safe */
/* number and BIGNUM = largest safe number. Use of these scaling */
/* factors is not guaranteed to reduce the condition number of A but */
/* works well in practice. */
/* Arguments */
/* ========= */
/* M (input) INTEGER */
/* The number of rows of the matrix A. M >= 0. */
/* N (input) INTEGER */
/* The number of columns of the matrix A. N >= 0. */
/* KL (input) INTEGER */
/* The number of subdiagonals within the band of A. KL >= 0. */
/* KU (input) INTEGER */
/* The number of superdiagonals within the band of A. KU >= 0. */
/* AB (input) DOUBLE PRECISION array, dimension (LDAB,N) */
/* The band matrix A, stored in rows 1 to KL+KU+1. The j-th */
/* column of A is stored in the j-th column of the array AB as */
/* follows: */
/* AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl). */
/* LDAB (input) INTEGER */
/* The leading dimension of the array AB. LDAB >= KL+KU+1. */
/* R (output) DOUBLE PRECISION array, dimension (M) */
/* If INFO = 0, or INFO > M, R contains the row scale factors */
/* for A. */
/* C (output) DOUBLE PRECISION array, dimension (N) */
/* If INFO = 0, C contains the column scale factors for A. */
/* ROWCND (output) DOUBLE PRECISION */
/* If INFO = 0 or INFO > M, ROWCND contains the ratio of the */
/* smallest R(i) to the largest R(i). If ROWCND >= 0.1 and */
/* AMAX is neither too large nor too small, it is not worth */
/* scaling by R. */
/* COLCND (output) DOUBLE PRECISION */
/* If INFO = 0, COLCND contains the ratio of the smallest */
/* C(i) to the largest C(i). If COLCND >= 0.1, it is not */
/* worth scaling by C. */
/* AMAX (output) DOUBLE PRECISION */
/* Absolute value of largest matrix element. If AMAX is very */
/* close to overflow or very close to underflow, the matrix */
/* should be scaled. */
/* INFO (output) INTEGER */
/* = 0: successful exit */
/* < 0: if INFO = -i, the i-th argument had an illegal value */
/* > 0: if INFO = i, and i is */
/* <= M: the i-th row of A is exactly zero */
/* > M: the (i-M)-th column of A is exactly zero */
/* ===================================================================== */
/* .. Parameters .. */
/* .. */
/* .. Local Scalars .. */
/* .. */
/* .. External Functions .. */
/* .. */
/* .. External Subroutines .. */
/* .. */
/* .. Intrinsic Functions .. */
/* .. */
/* .. Executable Statements .. */
/* Test the input parameters */
/* Parameter adjustments */
ab_dim1 = *ldab;
ab_offset = 1 + ab_dim1;
ab -= ab_offset;
--r__;
--c__;
/* Function Body */
*info = 0;
if (*m < 0) {
*info = -1;
} else if (*n < 0) {
*info = -2;
} else if (*kl < 0) {
*info = -3;
} else if (*ku < 0) {
*info = -4;
} else if (*ldab < *kl + *ku + 1) {
*info = -6;
}
if (*info != 0) {
i__1 = -(*info);
xerbla_("DGBEQU", &i__1, (ftnlen)6);
return 0;
}
/* Quick return if possible */
if (*m == 0 || *n == 0) {
*rowcnd = 1.;
*colcnd = 1.;
*amax = 0.;
return 0;
}
/* Get machine constants. */
smlnum = dlamch_("S", (ftnlen)1);
bignum = 1. / smlnum;
/* Compute row scale factors. */
i__1 = *m;
for (i__ = 1; i__ <= i__1; ++i__) {
r__[i__] = 0.;
/* L10: */
}
/* Find the maximum element in each row. */
kd = *ku + 1;
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
/* Computing MAX */
i__2 = j - *ku;
/* Computing MIN */
i__4 = j + *kl;
i__3 = min(i__4,*m);
for (i__ = max(i__2,1); i__ <= i__3; ++i__) {
/* Computing MAX */
d__2 = r__[i__], d__3 = (d__1 = ab[kd + i__ - j + j * ab_dim1],
abs(d__1));
r__[i__] = max(d__2,d__3);
/* L20: */
}
/* L30: */
}
/* Find the maximum and minimum scale factors. */
rcmin = bignum;
rcmax = 0.;
i__1 = *m;
for (i__ = 1; i__ <= i__1; ++i__) {
/* Computing MAX */
d__1 = rcmax, d__2 = r__[i__];
rcmax = max(d__1,d__2);
/* Computing MIN */
d__1 = rcmin, d__2 = r__[i__];
rcmin = min(d__1,d__2);
/* L40: */
}
*amax = rcmax;
if (rcmin == 0.) {
/* Find the first zero scale factor and return an error code. */
i__1 = *m;
for (i__ = 1; i__ <= i__1; ++i__) {
if (r__[i__] == 0.) {
*info = i__;
return 0;
}
/* L50: */
}
} else {
/* Invert the scale factors. */
i__1 = *m;
for (i__ = 1; i__ <= i__1; ++i__) {
/* Computing MIN */
/* Computing MAX */
d__2 = r__[i__];
d__1 = max(d__2,smlnum);
r__[i__] = 1. / min(d__1,bignum);
/* L60: */
}
/* Compute ROWCND = min(R(I)) / max(R(I)) */
*rowcnd = max(rcmin,smlnum) / min(rcmax,bignum);
}
/* Compute column scale factors */
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
c__[j] = 0.;
/* L70: */
}
/* Find the maximum element in each column, */
/* assuming the row scaling computed above. */
kd = *ku + 1;
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
/* Computing MAX */
i__3 = j - *ku;
/* Computing MIN */
i__4 = j + *kl;
i__2 = min(i__4,*m);
for (i__ = max(i__3,1); i__ <= i__2; ++i__) {
/* Computing MAX */
d__2 = c__[j], d__3 = (d__1 = ab[kd + i__ - j + j * ab_dim1], abs(
d__1)) * r__[i__];
c__[j] = max(d__2,d__3);
/* L80: */
}
/* L90: */
}
/* Find the maximum and minimum scale factors. */
rcmin = bignum;
rcmax = 0.;
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
/* Computing MIN */
d__1 = rcmin, d__2 = c__[j];
rcmin = min(d__1,d__2);
/* Computing MAX */
d__1 = rcmax, d__2 = c__[j];
rcmax = max(d__1,d__2);
/* L100: */
}
if (rcmin == 0.) {
/* Find the first zero scale factor and return an error code. */
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
if (c__[j] == 0.) {
*info = *m + j;
return 0;
}
/* L110: */
}
} else {
/* Invert the scale factors. */
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
/* Computing MIN */
/* Computing MAX */
d__2 = c__[j];
d__1 = max(d__2,smlnum);
c__[j] = 1. / min(d__1,bignum);
/* L120: */
}
/* Compute COLCND = min(C(J)) / max(C(J)) */
*colcnd = max(rcmin,smlnum) / min(rcmax,bignum);
}
return 0;
/* End of DGBEQU */
} /* dgbequ_ */

855
ext/f2c_lapack/dlatbs.c Normal file
View file

@ -0,0 +1,855 @@
/* dlatbs.f -- translated by f2c (version 20031025).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
/* Table of constant values */
static integer c__1 = 1;
static doublereal c_b36 = .5;
/* Subroutine */ int dlatbs_(char *uplo, char *trans, char *diag, char *
normin, integer *n, integer *kd, doublereal *ab, integer *ldab,
doublereal *x, doublereal *scale, doublereal *cnorm, integer *info,
ftnlen uplo_len, ftnlen trans_len, ftnlen diag_len, ftnlen normin_len)
{
/* System generated locals */
integer ab_dim1, ab_offset, i__1, i__2, i__3, i__4;
doublereal d__1, d__2, d__3;
/* Local variables */
static integer i__, j;
static doublereal xj, rec, tjj;
static integer jinc, jlen;
extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
integer *);
static doublereal xbnd;
static integer imax;
static doublereal tmax, tjjs, xmax, grow, sumj;
extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
integer *);
static integer maind;
extern logical lsame_(char *, char *, ftnlen, ftnlen);
static doublereal tscal, uscal;
extern doublereal dasum_(integer *, doublereal *, integer *);
static integer jlast;
extern /* Subroutine */ int dtbsv_(char *, char *, char *, integer *,
integer *, doublereal *, integer *, doublereal *, integer *,
ftnlen, ftnlen, ftnlen), daxpy_(integer *, doublereal *,
doublereal *, integer *, doublereal *, integer *);
static logical upper;
extern doublereal dlamch_(char *, ftnlen);
extern integer idamax_(integer *, doublereal *, integer *);
extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
static doublereal bignum;
static logical notran;
static integer jfirst;
static doublereal smlnum;
static logical nounit;
/* -- LAPACK auxiliary routine (version 3.0) -- */
/* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
/* Courant Institute, Argonne National Lab, and Rice University */
/* June 30, 1992 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* DLATBS solves one of the triangular systems */
/* A *x = s*b or A'*x = s*b */
/* with scaling to prevent overflow, where A is an upper or lower */
/* triangular band matrix. Here A' denotes the transpose of A, x and b */
/* are n-element vectors, and s is a scaling factor, usually less than */
/* or equal to 1, chosen so that the components of x will be less than */
/* the overflow threshold. If the unscaled problem will not cause */
/* overflow, the Level 2 BLAS routine DTBSV is called. If the matrix A */
/* is singular (A(j,j) = 0 for some j), then s is set to 0 and a */
/* non-trivial solution to A*x = 0 is returned. */
/* Arguments */
/* ========= */
/* UPLO (input) CHARACTER*1 */
/* Specifies whether the matrix A is upper or lower triangular. */
/* = 'U': Upper triangular */
/* = 'L': Lower triangular */
/* TRANS (input) CHARACTER*1 */
/* Specifies the operation applied to A. */
/* = 'N': Solve A * x = s*b (No transpose) */
/* = 'T': Solve A'* x = s*b (Transpose) */
/* = 'C': Solve A'* x = s*b (Conjugate transpose = Transpose) */
/* DIAG (input) CHARACTER*1 */
/* Specifies whether or not the matrix A is unit triangular. */
/* = 'N': Non-unit triangular */
/* = 'U': Unit triangular */
/* NORMIN (input) CHARACTER*1 */
/* Specifies whether CNORM has been set or not. */
/* = 'Y': CNORM contains the column norms on entry */
/* = 'N': CNORM is not set on entry. On exit, the norms will */
/* be computed and stored in CNORM. */
/* N (input) INTEGER */
/* The order of the matrix A. N >= 0. */
/* KD (input) INTEGER */
/* The number of subdiagonals or superdiagonals in the */
/* triangular matrix A. KD >= 0. */
/* AB (input) DOUBLE PRECISION array, dimension (LDAB,N) */
/* The upper or lower triangular band matrix A, stored in the */
/* first KD+1 rows of the array. The j-th column of A is stored */
/* in the j-th column of the array AB as follows: */
/* if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; */
/* if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). */
/* LDAB (input) INTEGER */
/* The leading dimension of the array AB. LDAB >= KD+1. */
/* X (input/output) DOUBLE PRECISION array, dimension (N) */
/* On entry, the right hand side b of the triangular system. */
/* On exit, X is overwritten by the solution vector x. */
/* SCALE (output) DOUBLE PRECISION */
/* The scaling factor s for the triangular system */
/* A * x = s*b or A'* x = s*b. */
/* If SCALE = 0, the matrix A is singular or badly scaled, and */
/* the vector x is an exact or approximate solution to A*x = 0. */
/* CNORM (input or output) DOUBLE PRECISION array, dimension (N) */
/* If NORMIN = 'Y', CNORM is an input argument and CNORM(j) */
/* contains the norm of the off-diagonal part of the j-th column */
/* of A. If TRANS = 'N', CNORM(j) must be greater than or equal */
/* to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j) */
/* must be greater than or equal to the 1-norm. */
/* If NORMIN = 'N', CNORM is an output argument and CNORM(j) */
/* returns the 1-norm of the offdiagonal part of the j-th column */
/* of A. */
/* INFO (output) INTEGER */
/* = 0: successful exit */
/* < 0: if INFO = -k, the k-th argument had an illegal value */
/* Further Details */
/* ======= ======= */
/* A rough bound on x is computed; if that is less than overflow, DTBSV */
/* is called, otherwise, specific code is used which checks for possible */
/* overflow or divide-by-zero at every operation. */
/* A columnwise scheme is used for solving A*x = b. The basic algorithm */
/* if A is lower triangular is */
/* x[1:n] := b[1:n] */
/* for j = 1, ..., n */
/* x(j) := x(j) / A(j,j) */
/* x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j] */
/* end */
/* Define bounds on the components of x after j iterations of the loop: */
/* M(j) = bound on x[1:j] */
/* G(j) = bound on x[j+1:n] */
/* Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}. */
/* Then for iteration j+1 we have */
/* M(j+1) <= G(j) / | A(j+1,j+1) | */
/* G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] | */
/* <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | ) */
/* where CNORM(j+1) is greater than or equal to the infinity-norm of */
/* column j+1 of A, not counting the diagonal. Hence */
/* G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | ) */
/* 1<=i<=j */
/* and */
/* |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| ) */
/* 1<=i< j */
/* Since |x(j)| <= M(j), we use the Level 2 BLAS routine DTBSV if the */
/* reciprocal of the largest M(j), j=1,..,n, is larger than */
/* max(underflow, 1/overflow). */
/* The bound on x(j) is also used to determine when a step in the */
/* columnwise method can be performed without fear of overflow. If */
/* the computed bound is greater than a large constant, x is scaled to */
/* prevent overflow, but if the bound overflows, x is set to 0, x(j) to */
/* 1, and scale to 0, and a non-trivial solution to A*x = 0 is found. */
/* Similarly, a row-wise scheme is used to solve A'*x = b. The basic */
/* algorithm for A upper triangular is */
/* for j = 1, ..., n */
/* x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j) */
/* end */
/* We simultaneously compute two bounds */
/* G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j */
/* M(j) = bound on x(i), 1<=i<=j */
/* The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we */
/* add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1. */
/* Then the bound on x(j) is */
/* M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) | */
/* <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| ) */
/* 1<=i<=j */
/* and we can safely call DTBSV if 1/M(n) and 1/G(n) are both greater */
/* than max(underflow, 1/overflow). */
/* ===================================================================== */
/* .. Parameters .. */
/* .. */
/* .. Local Scalars .. */
/* .. */
/* .. External Functions .. */
/* .. */
/* .. External Subroutines .. */
/* .. */
/* .. Intrinsic Functions .. */
/* .. */
/* .. Executable Statements .. */
/* Parameter adjustments */
ab_dim1 = *ldab;
ab_offset = 1 + ab_dim1;
ab -= ab_offset;
--x;
--cnorm;
/* Function Body */
*info = 0;
upper = lsame_(uplo, "U", (ftnlen)1, (ftnlen)1);
notran = lsame_(trans, "N", (ftnlen)1, (ftnlen)1);
nounit = lsame_(diag, "N", (ftnlen)1, (ftnlen)1);
/* Test the input parameters. */
if (! upper && ! lsame_(uplo, "L", (ftnlen)1, (ftnlen)1)) {
*info = -1;
} else if (! notran && ! lsame_(trans, "T", (ftnlen)1, (ftnlen)1) && !
lsame_(trans, "C", (ftnlen)1, (ftnlen)1)) {
*info = -2;
} else if (! nounit && ! lsame_(diag, "U", (ftnlen)1, (ftnlen)1)) {
*info = -3;
} else if (! lsame_(normin, "Y", (ftnlen)1, (ftnlen)1) && ! lsame_(normin,
"N", (ftnlen)1, (ftnlen)1)) {
*info = -4;
} else if (*n < 0) {
*info = -5;
} else if (*kd < 0) {
*info = -6;
} else if (*ldab < *kd + 1) {
*info = -8;
}
if (*info != 0) {
i__1 = -(*info);
xerbla_("DLATBS", &i__1, (ftnlen)6);
return 0;
}
/* Quick return if possible */
if (*n == 0) {
return 0;
}
/* Determine machine dependent parameters to control overflow. */
smlnum = dlamch_("Safe minimum", (ftnlen)12) / dlamch_("Precision", (
ftnlen)9);
bignum = 1. / smlnum;
*scale = 1.;
if (lsame_(normin, "N", (ftnlen)1, (ftnlen)1)) {
/* Compute the 1-norm of each column, not including the diagonal. */
if (upper) {
/* A is upper triangular. */
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
/* Computing MIN */
i__2 = *kd, i__3 = j - 1;
jlen = min(i__2,i__3);
cnorm[j] = dasum_(&jlen, &ab[*kd + 1 - jlen + j * ab_dim1], &
c__1);
/* L10: */
}
} else {
/* A is lower triangular. */
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
/* Computing MIN */
i__2 = *kd, i__3 = *n - j;
jlen = min(i__2,i__3);
if (jlen > 0) {
cnorm[j] = dasum_(&jlen, &ab[j * ab_dim1 + 2], &c__1);
} else {
cnorm[j] = 0.;
}
/* L20: */
}
}
}
/* Scale the column norms by TSCAL if the maximum element in CNORM is */
/* greater than BIGNUM. */
imax = idamax_(n, &cnorm[1], &c__1);
tmax = cnorm[imax];
if (tmax <= bignum) {
tscal = 1.;
} else {
tscal = 1. / (smlnum * tmax);
dscal_(n, &tscal, &cnorm[1], &c__1);
}
/* Compute a bound on the computed solution vector to see if the */
/* Level 2 BLAS routine DTBSV can be used. */
j = idamax_(n, &x[1], &c__1);
xmax = (d__1 = x[j], abs(d__1));
xbnd = xmax;
if (notran) {
/* Compute the growth in A * x = b. */
if (upper) {
jfirst = *n;
jlast = 1;
jinc = -1;
maind = *kd + 1;
} else {
jfirst = 1;
jlast = *n;
jinc = 1;
maind = 1;
}
if (tscal != 1.) {
grow = 0.;
goto L50;
}
if (nounit) {
/* A is non-unit triangular. */
/* Compute GROW = 1/G(j) and XBND = 1/M(j). */
/* Initially, G(0) = max{x(i), i=1,...,n}. */
grow = 1. / max(xbnd,smlnum);
xbnd = grow;
i__1 = jlast;
i__2 = jinc;
for (j = jfirst; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
/* Exit the loop if the growth factor is too small. */
if (grow <= smlnum) {
goto L50;
}
/* M(j) = G(j-1) / abs(A(j,j)) */
tjj = (d__1 = ab[maind + j * ab_dim1], abs(d__1));
/* Computing MIN */
d__1 = xbnd, d__2 = min(1.,tjj) * grow;
xbnd = min(d__1,d__2);
if (tjj + cnorm[j] >= smlnum) {
/* G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) ) */
grow *= tjj / (tjj + cnorm[j]);
} else {
/* G(j) could overflow, set GROW to 0. */
grow = 0.;
}
/* L30: */
}
grow = xbnd;
} else {
/* A is unit triangular. */
/* Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}. */
/* Computing MIN */
d__1 = 1., d__2 = 1. / max(xbnd,smlnum);
grow = min(d__1,d__2);
i__2 = jlast;
i__1 = jinc;
for (j = jfirst; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
/* Exit the loop if the growth factor is too small. */
if (grow <= smlnum) {
goto L50;
}
/* G(j) = G(j-1)*( 1 + CNORM(j) ) */
grow *= 1. / (cnorm[j] + 1.);
/* L40: */
}
}
L50:
;
} else {
/* Compute the growth in A' * x = b. */
if (upper) {
jfirst = 1;
jlast = *n;
jinc = 1;
maind = *kd + 1;
} else {
jfirst = *n;
jlast = 1;
jinc = -1;
maind = 1;
}
if (tscal != 1.) {
grow = 0.;
goto L80;
}
if (nounit) {
/* A is non-unit triangular. */
/* Compute GROW = 1/G(j) and XBND = 1/M(j). */
/* Initially, M(0) = max{x(i), i=1,...,n}. */
grow = 1. / max(xbnd,smlnum);
xbnd = grow;
i__1 = jlast;
i__2 = jinc;
for (j = jfirst; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
/* Exit the loop if the growth factor is too small. */
if (grow <= smlnum) {
goto L80;
}
/* G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) ) */
xj = cnorm[j] + 1.;
/* Computing MIN */
d__1 = grow, d__2 = xbnd / xj;
grow = min(d__1,d__2);
/* M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j)) */
tjj = (d__1 = ab[maind + j * ab_dim1], abs(d__1));
if (xj > tjj) {
xbnd *= tjj / xj;
}
/* L60: */
}
grow = min(grow,xbnd);
} else {
/* A is unit triangular. */
/* Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}. */
/* Computing MIN */
d__1 = 1., d__2 = 1. / max(xbnd,smlnum);
grow = min(d__1,d__2);
i__2 = jlast;
i__1 = jinc;
for (j = jfirst; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
/* Exit the loop if the growth factor is too small. */
if (grow <= smlnum) {
goto L80;
}
/* G(j) = ( 1 + CNORM(j) )*G(j-1) */
xj = cnorm[j] + 1.;
grow /= xj;
/* L70: */
}
}
L80:
;
}
if (grow * tscal > smlnum) {
/* Use the Level 2 BLAS solve if the reciprocal of the bound on */
/* elements of X is not too small. */
dtbsv_(uplo, trans, diag, n, kd, &ab[ab_offset], ldab, &x[1], &c__1, (
ftnlen)1, (ftnlen)1, (ftnlen)1);
} else {
/* Use a Level 1 BLAS solve, scaling intermediate results. */
if (xmax > bignum) {
/* Scale X so that its components are less than or equal to */
/* BIGNUM in absolute value. */
*scale = bignum / xmax;
dscal_(n, scale, &x[1], &c__1);
xmax = bignum;
}
if (notran) {
/* Solve A * x = b */
i__1 = jlast;
i__2 = jinc;
for (j = jfirst; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
/* Compute x(j) = b(j) / A(j,j), scaling x if necessary. */
xj = (d__1 = x[j], abs(d__1));
if (nounit) {
tjjs = ab[maind + j * ab_dim1] * tscal;
} else {
tjjs = tscal;
if (tscal == 1.) {
goto L100;
}
}
tjj = abs(tjjs);
if (tjj > smlnum) {
/* abs(A(j,j)) > SMLNUM: */
if (tjj < 1.) {
if (xj > tjj * bignum) {
/* Scale x by 1/b(j). */
rec = 1. / xj;
dscal_(n, &rec, &x[1], &c__1);
*scale *= rec;
xmax *= rec;
}
}
x[j] /= tjjs;
xj = (d__1 = x[j], abs(d__1));
} else if (tjj > 0.) {
/* 0 < abs(A(j,j)) <= SMLNUM: */
if (xj > tjj * bignum) {
/* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM */
/* to avoid overflow when dividing by A(j,j). */
rec = tjj * bignum / xj;
if (cnorm[j] > 1.) {
/* Scale by 1/CNORM(j) to avoid overflow when */
/* multiplying x(j) times column j. */
rec /= cnorm[j];
}
dscal_(n, &rec, &x[1], &c__1);
*scale *= rec;
xmax *= rec;
}
x[j] /= tjjs;
xj = (d__1 = x[j], abs(d__1));
} else {
/* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and */
/* scale = 0, and compute a solution to A*x = 0. */
i__3 = *n;
for (i__ = 1; i__ <= i__3; ++i__) {
x[i__] = 0.;
/* L90: */
}
x[j] = 1.;
xj = 1.;
*scale = 0.;
xmax = 0.;
}
L100:
/* Scale x if necessary to avoid overflow when adding a */
/* multiple of column j of A. */
if (xj > 1.) {
rec = 1. / xj;
if (cnorm[j] > (bignum - xmax) * rec) {
/* Scale x by 1/(2*abs(x(j))). */
rec *= .5;
dscal_(n, &rec, &x[1], &c__1);
*scale *= rec;
}
} else if (xj * cnorm[j] > bignum - xmax) {
/* Scale x by 1/2. */
dscal_(n, &c_b36, &x[1], &c__1);
*scale *= .5;
}
if (upper) {
if (j > 1) {
/* Compute the update */
/* x(max(1,j-kd):j-1) := x(max(1,j-kd):j-1) - */
/* x(j)* A(max(1,j-kd):j-1,j) */
/* Computing MIN */
i__3 = *kd, i__4 = j - 1;
jlen = min(i__3,i__4);
d__1 = -x[j] * tscal;
daxpy_(&jlen, &d__1, &ab[*kd + 1 - jlen + j * ab_dim1]
, &c__1, &x[j - jlen], &c__1);
i__3 = j - 1;
i__ = idamax_(&i__3, &x[1], &c__1);
xmax = (d__1 = x[i__], abs(d__1));
}
} else if (j < *n) {
/* Compute the update */
/* x(j+1:min(j+kd,n)) := x(j+1:min(j+kd,n)) - */
/* x(j) * A(j+1:min(j+kd,n),j) */
/* Computing MIN */
i__3 = *kd, i__4 = *n - j;
jlen = min(i__3,i__4);
if (jlen > 0) {
d__1 = -x[j] * tscal;
daxpy_(&jlen, &d__1, &ab[j * ab_dim1 + 2], &c__1, &x[
j + 1], &c__1);
}
i__3 = *n - j;
i__ = j + idamax_(&i__3, &x[j + 1], &c__1);
xmax = (d__1 = x[i__], abs(d__1));
}
/* L110: */
}
} else {
/* Solve A' * x = b */
i__2 = jlast;
i__1 = jinc;
for (j = jfirst; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
/* Compute x(j) = b(j) - sum A(k,j)*x(k). */
/* k<>j */
xj = (d__1 = x[j], abs(d__1));
uscal = tscal;
rec = 1. / max(xmax,1.);
if (cnorm[j] > (bignum - xj) * rec) {
/* If x(j) could overflow, scale x by 1/(2*XMAX). */
rec *= .5;
if (nounit) {
tjjs = ab[maind + j * ab_dim1] * tscal;
} else {
tjjs = tscal;
}
tjj = abs(tjjs);
if (tjj > 1.) {
/* Divide by A(j,j) when scaling x if A(j,j) > 1. */
/* Computing MIN */
d__1 = 1., d__2 = rec * tjj;
rec = min(d__1,d__2);
uscal /= tjjs;
}
if (rec < 1.) {
dscal_(n, &rec, &x[1], &c__1);
*scale *= rec;
xmax *= rec;
}
}
sumj = 0.;
if (uscal == 1.) {
/* If the scaling needed for A in the dot product is 1, */
/* call DDOT to perform the dot product. */
if (upper) {
/* Computing MIN */
i__3 = *kd, i__4 = j - 1;
jlen = min(i__3,i__4);
sumj = ddot_(&jlen, &ab[*kd + 1 - jlen + j * ab_dim1],
&c__1, &x[j - jlen], &c__1);
} else {
/* Computing MIN */
i__3 = *kd, i__4 = *n - j;
jlen = min(i__3,i__4);
if (jlen > 0) {
sumj = ddot_(&jlen, &ab[j * ab_dim1 + 2], &c__1, &
x[j + 1], &c__1);
}
}
} else {
/* Otherwise, use in-line code for the dot product. */
if (upper) {
/* Computing MIN */
i__3 = *kd, i__4 = j - 1;
jlen = min(i__3,i__4);
i__3 = jlen;
for (i__ = 1; i__ <= i__3; ++i__) {
sumj += ab[*kd + i__ - jlen + j * ab_dim1] *
uscal * x[j - jlen - 1 + i__];
/* L120: */
}
} else {
/* Computing MIN */
i__3 = *kd, i__4 = *n - j;
jlen = min(i__3,i__4);
i__3 = jlen;
for (i__ = 1; i__ <= i__3; ++i__) {
sumj += ab[i__ + 1 + j * ab_dim1] * uscal * x[j +
i__];
/* L130: */
}
}
}
if (uscal == tscal) {
/* Compute x(j) := ( x(j) - sumj ) / A(j,j) if 1/A(j,j) */
/* was not used to scale the dotproduct. */
x[j] -= sumj;
xj = (d__1 = x[j], abs(d__1));
if (nounit) {
/* Compute x(j) = x(j) / A(j,j), scaling if necessary. */
tjjs = ab[maind + j * ab_dim1] * tscal;
} else {
tjjs = tscal;
if (tscal == 1.) {
goto L150;
}
}
tjj = abs(tjjs);
if (tjj > smlnum) {
/* abs(A(j,j)) > SMLNUM: */
if (tjj < 1.) {
if (xj > tjj * bignum) {
/* Scale X by 1/abs(x(j)). */
rec = 1. / xj;
dscal_(n, &rec, &x[1], &c__1);
*scale *= rec;
xmax *= rec;
}
}
x[j] /= tjjs;
} else if (tjj > 0.) {
/* 0 < abs(A(j,j)) <= SMLNUM: */
if (xj > tjj * bignum) {
/* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM. */
rec = tjj * bignum / xj;
dscal_(n, &rec, &x[1], &c__1);
*scale *= rec;
xmax *= rec;
}
x[j] /= tjjs;
} else {
/* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and */
/* scale = 0, and compute a solution to A'*x = 0. */
i__3 = *n;
for (i__ = 1; i__ <= i__3; ++i__) {
x[i__] = 0.;
/* L140: */
}
x[j] = 1.;
*scale = 0.;
xmax = 0.;
}
L150:
;
} else {
/* Compute x(j) := x(j) / A(j,j) - sumj if the dot */
/* product has already been divided by 1/A(j,j). */
x[j] = x[j] / tjjs - sumj;
}
/* Computing MAX */
d__2 = xmax, d__3 = (d__1 = x[j], abs(d__1));
xmax = max(d__2,d__3);
/* L160: */
}
}
*scale /= tscal;
}
/* Scale the column norms by 1/TSCAL for return. */
if (tscal != 1.) {
d__1 = 1. / tscal;
dscal_(n, &d__1, &cnorm[1], &c__1);
}
return 0;
/* End of DLATBS */
} /* dlatbs_ */

2
ext/lapack/Makefile.in
View file

@ -18,6 +18,7 @@ F_FLAGS = @FFLAGS@ $(PIC_FLAG)
OBJS = \
dbdsqr.o \
dgbcon.o \
dgbtrf.o \
dgbtf2.o \
dgbtrs.o \
@ -59,6 +60,7 @@ dlasrt.o \
dlassq.o \
dlasv2.o \
dlaswp.o \
dlatbs.o \
dorg2r.o \
dorgbr.o \
dorgl2.o \

222
ext/lapack/dgbcon.f Normal file
View file

@ -0,0 +1,222 @@
SUBROUTINE DGBCON( NORM, N, KL, KU, AB, LDAB, IPIV, ANORM, RCOND,
$ WORK, IWORK, INFO )
*
* -- LAPACK routine (version 3.0) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* September 30, 1994
*
* .. Scalar Arguments ..
CHARACTER NORM
INTEGER INFO, KL, KU, LDAB, N
DOUBLE PRECISION ANORM, RCOND
* ..
* .. Array Arguments ..
INTEGER IPIV( * ), IWORK( * )
DOUBLE PRECISION AB( LDAB, * ), WORK( * )
* ..
*
* Purpose
* =======
*
* DGBCON estimates the reciprocal of the condition number of a real
* general band matrix A, in either the 1-norm or the infinity-norm,
* using the LU factorization computed by DGBTRF.
*
* An estimate is obtained for norm(inv(A)), and the reciprocal of the
* condition number is computed as
* RCOND = 1 / ( norm(A) * norm(inv(A)) ).
*
* Arguments
* =========
*
* NORM (input) CHARACTER*1
* Specifies whether the 1-norm condition number or the
* infinity-norm condition number is required:
* = '1' or 'O': 1-norm;
* = 'I': Infinity-norm.
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* KL (input) INTEGER
* The number of subdiagonals within the band of A. KL >= 0.
*
* KU (input) INTEGER
* The number of superdiagonals within the band of A. KU >= 0.
*
* AB (input) DOUBLE PRECISION array, dimension (LDAB,N)
* Details of the LU factorization of the band matrix A, as
* computed by DGBTRF. U is stored as an upper triangular band
* matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and
* the multipliers used during the factorization are stored in
* rows KL+KU+2 to 2*KL+KU+1.
*
* LDAB (input) INTEGER
* The leading dimension of the array AB. LDAB >= 2*KL+KU+1.
*
* IPIV (input) INTEGER array, dimension (N)
* The pivot indices; for 1 <= i <= N, row i of the matrix was
* interchanged with row IPIV(i).
*
* ANORM (input) DOUBLE PRECISION
* If NORM = '1' or 'O', the 1-norm of the original matrix A.
* If NORM = 'I', the infinity-norm of the original matrix A.
*
* RCOND (output) DOUBLE PRECISION
* The reciprocal of the condition number of the matrix A,
* computed as RCOND = 1/(norm(A) * norm(inv(A))).
*
* WORK (workspace) DOUBLE PRECISION array, dimension (3*N)
*
* IWORK (workspace) INTEGER array, dimension (N)
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
*
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ONE, ZERO
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
* ..
* .. Local Scalars ..
LOGICAL LNOTI, ONENRM
CHARACTER NORMIN
INTEGER IX, J, JP, KASE, KASE1, KD, LM
DOUBLE PRECISION AINVNM, SCALE, SMLNUM, T
* ..
* .. External Functions ..
LOGICAL LSAME
INTEGER IDAMAX
DOUBLE PRECISION DDOT, DLAMCH
EXTERNAL LSAME, IDAMAX, DDOT, DLAMCH
* ..
* .. External Subroutines ..
EXTERNAL DAXPY, DLACON, DLATBS, DRSCL, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MIN
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' )
IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN
INFO = -1
ELSE IF( N.LT.0 ) THEN
INFO = -2
ELSE IF( KL.LT.0 ) THEN
INFO = -3
ELSE IF( KU.LT.0 ) THEN
INFO = -4
ELSE IF( LDAB.LT.2*KL+KU+1 ) THEN
INFO = -6
ELSE IF( ANORM.LT.ZERO ) THEN
INFO = -8
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DGBCON', -INFO )
RETURN
END IF
*
* Quick return if possible
*
RCOND = ZERO
IF( N.EQ.0 ) THEN
RCOND = ONE
RETURN
ELSE IF( ANORM.EQ.ZERO ) THEN
RETURN
END IF
*
SMLNUM = DLAMCH( 'Safe minimum' )
*
* Estimate the norm of inv(A).
*
AINVNM = ZERO
NORMIN = 'N'
IF( ONENRM ) THEN
KASE1 = 1
ELSE
KASE1 = 2
END IF
KD = KL + KU + 1
LNOTI = KL.GT.0
KASE = 0
10 CONTINUE
CALL DLACON( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE )
IF( KASE.NE.0 ) THEN
IF( KASE.EQ.KASE1 ) THEN
*
* Multiply by inv(L).
*
IF( LNOTI ) THEN
DO 20 J = 1, N - 1
LM = MIN( KL, N-J )
JP = IPIV( J )
T = WORK( JP )
IF( JP.NE.J ) THEN
WORK( JP ) = WORK( J )
WORK( J ) = T
END IF
CALL DAXPY( LM, -T, AB( KD+1, J ), 1, WORK( J+1 ), 1 )
20 CONTINUE
END IF
*
* Multiply by inv(U).
*
CALL DLATBS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
$ KL+KU, AB, LDAB, WORK, SCALE, WORK( 2*N+1 ),
$ INFO )
ELSE
*
* Multiply by inv(U').
*
CALL DLATBS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N,
$ KL+KU, AB, LDAB, WORK, SCALE, WORK( 2*N+1 ),
$ INFO )
*
* Multiply by inv(L').
*
IF( LNOTI ) THEN
DO 30 J = N - 1, 1, -1
LM = MIN( KL, N-J )
WORK( J ) = WORK( J ) - DDOT( LM, AB( KD+1, J ), 1,
$ WORK( J+1 ), 1 )
JP = IPIV( J )
IF( JP.NE.J ) THEN
T = WORK( JP )
WORK( JP ) = WORK( J )
WORK( J ) = T
END IF
30 CONTINUE
END IF
END IF
*
* Divide X by 1/SCALE if doing so will not cause overflow.
*
NORMIN = 'Y'
IF( SCALE.NE.ONE ) THEN
IX = IDAMAX( N, WORK, 1 )
IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
$ GO TO 40
CALL DRSCL( N, SCALE, WORK, 1 )
END IF
GO TO 10
END IF
*
* Compute the estimate of the reciprocal condition number.
*
IF( AINVNM.NE.ZERO )
$ RCOND = ( ONE / AINVNM ) / ANORM
*
40 CONTINUE
RETURN
*
* End of DGBCON
*
END

240
ext/lapack/dgbequ.f Normal file
View file

@ -0,0 +1,240 @@
SUBROUTINE DGBEQU( M, N, KL, KU, AB, LDAB, R, C, ROWCND, COLCND,
$ AMAX, INFO )
*
* -- LAPACK routine (version 3.0) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* March 31, 1993
*
* .. Scalar Arguments ..
INTEGER INFO, KL, KU, LDAB, M, N
DOUBLE PRECISION AMAX, COLCND, ROWCND
* ..
* .. Array Arguments ..
DOUBLE PRECISION AB( LDAB, * ), C( * ), R( * )
* ..
*
* Purpose
* =======
*
* DGBEQU computes row and column scalings intended to equilibrate an
* M-by-N band matrix A and reduce its condition number. R returns the
* row scale factors and C the column scale factors, chosen to try to
* make the largest element in each row and column of the matrix B with
* elements B(i,j)=R(i)*A(i,j)*C(j) have absolute value 1.
*
* R(i) and C(j) are restricted to be between SMLNUM = smallest safe
* number and BIGNUM = largest safe number. Use of these scaling
* factors is not guaranteed to reduce the condition number of A but
* works well in practice.
*
* Arguments
* =========
*
* M (input) INTEGER
* The number of rows of the matrix A. M >= 0.
*
* N (input) INTEGER
* The number of columns of the matrix A. N >= 0.
*
* KL (input) INTEGER
* The number of subdiagonals within the band of A. KL >= 0.
*
* KU (input) INTEGER
* The number of superdiagonals within the band of A. KU >= 0.
*
* AB (input) DOUBLE PRECISION array, dimension (LDAB,N)
* The band matrix A, stored in rows 1 to KL+KU+1. The j-th
* column of A is stored in the j-th column of the array AB as
* follows:
* AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl).
*
* LDAB (input) INTEGER
* The leading dimension of the array AB. LDAB >= KL+KU+1.
*
* R (output) DOUBLE PRECISION array, dimension (M)
* If INFO = 0, or INFO > M, R contains the row scale factors
* for A.
*
* C (output) DOUBLE PRECISION array, dimension (N)
* If INFO = 0, C contains the column scale factors for A.
*
* ROWCND (output) DOUBLE PRECISION
* If INFO = 0 or INFO > M, ROWCND contains the ratio of the
* smallest R(i) to the largest R(i). If ROWCND >= 0.1 and
* AMAX is neither too large nor too small, it is not worth
* scaling by R.
*
* COLCND (output) DOUBLE PRECISION
* If INFO = 0, COLCND contains the ratio of the smallest
* C(i) to the largest C(i). If COLCND >= 0.1, it is not
* worth scaling by C.
*
* AMAX (output) DOUBLE PRECISION
* Absolute value of largest matrix element. If AMAX is very
* close to overflow or very close to underflow, the matrix
* should be scaled.
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
* > 0: if INFO = i, and i is
* <= M: the i-th row of A is exactly zero
* > M: the (i-M)-th column of A is exactly zero
*
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ONE, ZERO
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
* ..
* .. Local Scalars ..
INTEGER I, J, KD
DOUBLE PRECISION BIGNUM, RCMAX, RCMIN, SMLNUM
* ..
* .. External Functions ..
DOUBLE PRECISION DLAMCH
EXTERNAL DLAMCH
* ..
* .. External Subroutines ..
EXTERNAL XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN
* ..
* .. Executable Statements ..
*
* Test the input parameters
*
INFO = 0
IF( M.LT.0 ) THEN
INFO = -1
ELSE IF( N.LT.0 ) THEN
INFO = -2
ELSE IF( KL.LT.0 ) THEN
INFO = -3
ELSE IF( KU.LT.0 ) THEN
INFO = -4
ELSE IF( LDAB.LT.KL+KU+1 ) THEN
INFO = -6
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DGBEQU', -INFO )
RETURN
END IF
*
* Quick return if possible
*
IF( M.EQ.0 .OR. N.EQ.0 ) THEN
ROWCND = ONE
COLCND = ONE
AMAX = ZERO
RETURN
END IF
*
* Get machine constants.
*
SMLNUM = DLAMCH( 'S' )
BIGNUM = ONE / SMLNUM
*
* Compute row scale factors.
*
DO 10 I = 1, M
R( I ) = ZERO
10 CONTINUE
*
* Find the maximum element in each row.
*
KD = KU + 1
DO 30 J = 1, N
DO 20 I = MAX( J-KU, 1 ), MIN( J+KL, M )
R( I ) = MAX( R( I ), ABS( AB( KD+I-J, J ) ) )
20 CONTINUE
30 CONTINUE
*
* Find the maximum and minimum scale factors.
*
RCMIN = BIGNUM
RCMAX = ZERO
DO 40 I = 1, M
RCMAX = MAX( RCMAX, R( I ) )
RCMIN = MIN( RCMIN, R( I ) )
40 CONTINUE
AMAX = RCMAX
*
IF( RCMIN.EQ.ZERO ) THEN
*
* Find the first zero scale factor and return an error code.
*
DO 50 I = 1, M
IF( R( I ).EQ.ZERO ) THEN
INFO = I
RETURN
END IF
50 CONTINUE
ELSE
*
* Invert the scale factors.
*
DO 60 I = 1, M
R( I ) = ONE / MIN( MAX( R( I ), SMLNUM ), BIGNUM )
60 CONTINUE
*
* Compute ROWCND = min(R(I)) / max(R(I))
*
ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
END IF
*
* Compute column scale factors
*
DO 70 J = 1, N
C( J ) = ZERO
70 CONTINUE
*
* Find the maximum element in each column,
* assuming the row scaling computed above.
*
KD = KU + 1
DO 90 J = 1, N
DO 80 I = MAX( J-KU, 1 ), MIN( J+KL, M )
C( J ) = MAX( C( J ), ABS( AB( KD+I-J, J ) )*R( I ) )
80 CONTINUE
90 CONTINUE
*
* Find the maximum and minimum scale factors.
*
RCMIN = BIGNUM
RCMAX = ZERO
DO 100 J = 1, N
RCMIN = MIN( RCMIN, C( J ) )
RCMAX = MAX( RCMAX, C( J ) )
100 CONTINUE
*
IF( RCMIN.EQ.ZERO ) THEN
*
* Find the first zero scale factor and return an error code.
*
DO 110 J = 1, N
IF( C( J ).EQ.ZERO ) THEN
INFO = M + J
RETURN
END IF
110 CONTINUE
ELSE
*
* Invert the scale factors.
*
DO 120 J = 1, N
C( J ) = ONE / MIN( MAX( C( J ), SMLNUM ), BIGNUM )
120 CONTINUE
*
* Compute COLCND = min(C(J)) / max(C(J))
*
COLCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
END IF
*
RETURN
*
* End of DGBEQU
*
END

724
ext/lapack/dlatbs.f Normal file
View file

@ -0,0 +1,724 @@
SUBROUTINE DLATBS( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X,
$ SCALE, CNORM, INFO )
*
* -- LAPACK auxiliary routine (version 3.0) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* June 30, 1992
*
* .. Scalar Arguments ..
CHARACTER DIAG, NORMIN, TRANS, UPLO
INTEGER INFO, KD, LDAB, N
DOUBLE PRECISION SCALE
* ..
* .. Array Arguments ..
DOUBLE PRECISION AB( LDAB, * ), CNORM( * ), X( * )
* ..
*
* Purpose
* =======
*
* DLATBS solves one of the triangular systems
*
* A *x = s*b or A'*x = s*b
*
* with scaling to prevent overflow, where A is an upper or lower
* triangular band matrix. Here A' denotes the transpose of A, x and b
* are n-element vectors, and s is a scaling factor, usually less than
* or equal to 1, chosen so that the components of x will be less than
* the overflow threshold. If the unscaled problem will not cause
* overflow, the Level 2 BLAS routine DTBSV is called. If the matrix A
* is singular (A(j,j) = 0 for some j), then s is set to 0 and a
* non-trivial solution to A*x = 0 is returned.
*
* Arguments
* =========
*
* UPLO (input) CHARACTER*1
* Specifies whether the matrix A is upper or lower triangular.
* = 'U': Upper triangular
* = 'L': Lower triangular
*
* TRANS (input) CHARACTER*1
* Specifies the operation applied to A.
* = 'N': Solve A * x = s*b (No transpose)
* = 'T': Solve A'* x = s*b (Transpose)
* = 'C': Solve A'* x = s*b (Conjugate transpose = Transpose)
*
* DIAG (input) CHARACTER*1
* Specifies whether or not the matrix A is unit triangular.
* = 'N': Non-unit triangular
* = 'U': Unit triangular
*
* NORMIN (input) CHARACTER*1
* Specifies whether CNORM has been set or not.
* = 'Y': CNORM contains the column norms on entry
* = 'N': CNORM is not set on entry. On exit, the norms will
* be computed and stored in CNORM.
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* KD (input) INTEGER
* The number of subdiagonals or superdiagonals in the
* triangular matrix A. KD >= 0.
*
* AB (input) DOUBLE PRECISION array, dimension (LDAB,N)
* The upper or lower triangular band matrix A, stored in the
* first KD+1 rows of the array. The j-th column of A is stored
* in the j-th column of the array AB as follows:
* if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
* if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
*
* LDAB (input) INTEGER
* The leading dimension of the array AB. LDAB >= KD+1.
*
* X (input/output) DOUBLE PRECISION array, dimension (N)
* On entry, the right hand side b of the triangular system.
* On exit, X is overwritten by the solution vector x.
*
* SCALE (output) DOUBLE PRECISION
* The scaling factor s for the triangular system
* A * x = s*b or A'* x = s*b.
* If SCALE = 0, the matrix A is singular or badly scaled, and
* the vector x is an exact or approximate solution to A*x = 0.
*
* CNORM (input or output) DOUBLE PRECISION array, dimension (N)
*
* If NORMIN = 'Y', CNORM is an input argument and CNORM(j)
* contains the norm of the off-diagonal part of the j-th column
* of A. If TRANS = 'N', CNORM(j) must be greater than or equal
* to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)
* must be greater than or equal to the 1-norm.
*
* If NORMIN = 'N', CNORM is an output argument and CNORM(j)
* returns the 1-norm of the offdiagonal part of the j-th column
* of A.
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -k, the k-th argument had an illegal value
*
* Further Details
* ======= =======
*
* A rough bound on x is computed; if that is less than overflow, DTBSV
* is called, otherwise, specific code is used which checks for possible
* overflow or divide-by-zero at every operation.
*
* A columnwise scheme is used for solving A*x = b. The basic algorithm
* if A is lower triangular is
*
* x[1:n] := b[1:n]
* for j = 1, ..., n
* x(j) := x(j) / A(j,j)
* x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
* end
*
* Define bounds on the components of x after j iterations of the loop:
* M(j) = bound on x[1:j]
* G(j) = bound on x[j+1:n]
* Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
*
* Then for iteration j+1 we have
* M(j+1) <= G(j) / | A(j+1,j+1) |
* G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
* <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
*
* where CNORM(j+1) is greater than or equal to the infinity-norm of
* column j+1 of A, not counting the diagonal. Hence
*
* G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
* 1<=i<=j
* and
*
* |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
* 1<=i< j
*
* Since |x(j)| <= M(j), we use the Level 2 BLAS routine DTBSV if the
* reciprocal of the largest M(j), j=1,..,n, is larger than
* max(underflow, 1/overflow).
*
* The bound on x(j) is also used to determine when a step in the
* columnwise method can be performed without fear of overflow. If
* the computed bound is greater than a large constant, x is scaled to
* prevent overflow, but if the bound overflows, x is set to 0, x(j) to
* 1, and scale to 0, and a non-trivial solution to A*x = 0 is found.
*
* Similarly, a row-wise scheme is used to solve A'*x = b. The basic
* algorithm for A upper triangular is
*
* for j = 1, ..., n
* x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j)
* end
*
* We simultaneously compute two bounds
* G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j
* M(j) = bound on x(i), 1<=i<=j
*
* The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we
* add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.
* Then the bound on x(j) is
*
* M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
*
* <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
* 1<=i<=j
*
* and we can safely call DTBSV if 1/M(n) and 1/G(n) are both greater
* than max(underflow, 1/overflow).
*
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO, HALF, ONE
PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0 )
* ..
* .. Local Scalars ..
LOGICAL NOTRAN, NOUNIT, UPPER
INTEGER I, IMAX, J, JFIRST, JINC, JLAST, JLEN, MAIND
DOUBLE PRECISION BIGNUM, GROW, REC, SMLNUM, SUMJ, TJJ, TJJS,
$ TMAX, TSCAL, USCAL, XBND, XJ, XMAX
* ..
* .. External Functions ..
LOGICAL LSAME
INTEGER IDAMAX
DOUBLE PRECISION DASUM, DDOT, DLAMCH
EXTERNAL LSAME, IDAMAX, DASUM, DDOT, DLAMCH
* ..
* .. External Subroutines ..
EXTERNAL DAXPY, DSCAL, DTBSV, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN
* ..
* .. Executable Statements ..
*
INFO = 0
UPPER = LSAME( UPLO, 'U' )
NOTRAN = LSAME( TRANS, 'N' )
NOUNIT = LSAME( DIAG, 'N' )
*
* Test the input parameters.
*
IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
INFO = -1
ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
$ LSAME( TRANS, 'C' ) ) THEN
INFO = -2
ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
INFO = -3
ELSE IF( .NOT.LSAME( NORMIN, 'Y' ) .AND. .NOT.
$ LSAME( NORMIN, 'N' ) ) THEN
INFO = -4
ELSE IF( N.LT.0 ) THEN
INFO = -5
ELSE IF( KD.LT.0 ) THEN
INFO = -6
ELSE IF( LDAB.LT.KD+1 ) THEN
INFO = -8
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DLATBS', -INFO )
RETURN
END IF
*
* Quick return if possible
*
IF( N.EQ.0 )
$ RETURN
*
* Determine machine dependent parameters to control overflow.
*
SMLNUM = DLAMCH( 'Safe minimum' ) / DLAMCH( 'Precision' )
BIGNUM = ONE / SMLNUM
SCALE = ONE
*
IF( LSAME( NORMIN, 'N' ) ) THEN
*
* Compute the 1-norm of each column, not including the diagonal.
*
IF( UPPER ) THEN
*
* A is upper triangular.
*
DO 10 J = 1, N
JLEN = MIN( KD, J-1 )
CNORM( J ) = DASUM( JLEN, AB( KD+1-JLEN, J ), 1 )
10 CONTINUE
ELSE
*
* A is lower triangular.
*
DO 20 J = 1, N
JLEN = MIN( KD, N-J )
IF( JLEN.GT.0 ) THEN
CNORM( J ) = DASUM( JLEN, AB( 2, J ), 1 )
ELSE
CNORM( J ) = ZERO
END IF
20 CONTINUE
END IF
END IF
*
* Scale the column norms by TSCAL if the maximum element in CNORM is
* greater than BIGNUM.
*
IMAX = IDAMAX( N, CNORM, 1 )
TMAX = CNORM( IMAX )
IF( TMAX.LE.BIGNUM ) THEN
TSCAL = ONE
ELSE
TSCAL = ONE / ( SMLNUM*TMAX )
CALL DSCAL( N, TSCAL, CNORM, 1 )
END IF
*
* Compute a bound on the computed solution vector to see if the
* Level 2 BLAS routine DTBSV can be used.
*
J = IDAMAX( N, X, 1 )
XMAX = ABS( X( J ) )
XBND = XMAX
IF( NOTRAN ) THEN
*
* Compute the growth in A * x = b.
*
IF( UPPER ) THEN
JFIRST = N
JLAST = 1
JINC = -1
MAIND = KD + 1
ELSE
JFIRST = 1
JLAST = N
JINC = 1
MAIND = 1
END IF
*
IF( TSCAL.NE.ONE ) THEN
GROW = ZERO
GO TO 50
END IF
*
IF( NOUNIT ) THEN
*
* A is non-unit triangular.
*
* Compute GROW = 1/G(j) and XBND = 1/M(j).
* Initially, G(0) = max{x(i), i=1,...,n}.
*
GROW = ONE / MAX( XBND, SMLNUM )
XBND = GROW
DO 30 J = JFIRST, JLAST, JINC
*
* Exit the loop if the growth factor is too small.
*
IF( GROW.LE.SMLNUM )
$ GO TO 50
*
* M(j) = G(j-1) / abs(A(j,j))
*
TJJ = ABS( AB( MAIND, J ) )
XBND = MIN( XBND, MIN( ONE, TJJ )*GROW )
IF( TJJ+CNORM( J ).GE.SMLNUM ) THEN
*
* G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
*
GROW = GROW*( TJJ / ( TJJ+CNORM( J ) ) )
ELSE
*
* G(j) could overflow, set GROW to 0.
*
GROW = ZERO
END IF
30 CONTINUE
GROW = XBND
ELSE
*
* A is unit triangular.
*
* Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
*
GROW = MIN( ONE, ONE / MAX( XBND, SMLNUM ) )
DO 40 J = JFIRST, JLAST, JINC
*
* Exit the loop if the growth factor is too small.
*
IF( GROW.LE.SMLNUM )
$ GO TO 50
*
* G(j) = G(j-1)*( 1 + CNORM(j) )
*
GROW = GROW*( ONE / ( ONE+CNORM( J ) ) )
40 CONTINUE
END IF
50 CONTINUE
*
ELSE
*
* Compute the growth in A' * x = b.
*
IF( UPPER ) THEN
JFIRST = 1
JLAST = N
JINC = 1
MAIND = KD + 1
ELSE
JFIRST = N
JLAST = 1
JINC = -1
MAIND = 1
END IF
*
IF( TSCAL.NE.ONE ) THEN
GROW = ZERO
GO TO 80
END IF
*
IF( NOUNIT ) THEN
*
* A is non-unit triangular.
*
* Compute GROW = 1/G(j) and XBND = 1/M(j).
* Initially, M(0) = max{x(i), i=1,...,n}.
*
GROW = ONE / MAX( XBND, SMLNUM )
XBND = GROW
DO 60 J = JFIRST, JLAST, JINC
*
* Exit the loop if the growth factor is too small.
*
IF( GROW.LE.SMLNUM )
$ GO TO 80
*
* G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
*
XJ = ONE + CNORM( J )
GROW = MIN( GROW, XBND / XJ )
*
* M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
*
TJJ = ABS( AB( MAIND, J ) )
IF( XJ.GT.TJJ )
$ XBND = XBND*( TJJ / XJ )
60 CONTINUE
GROW = MIN( GROW, XBND )
ELSE
*
* A is unit triangular.
*
* Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
*
GROW = MIN( ONE, ONE / MAX( XBND, SMLNUM ) )
DO 70 J = JFIRST, JLAST, JINC
*
* Exit the loop if the growth factor is too small.
*
IF( GROW.LE.SMLNUM )
$ GO TO 80
*
* G(j) = ( 1 + CNORM(j) )*G(j-1)
*
XJ = ONE + CNORM( J )
GROW = GROW / XJ
70 CONTINUE
END IF
80 CONTINUE
END IF
*
IF( ( GROW*TSCAL ).GT.SMLNUM ) THEN
*
* Use the Level 2 BLAS solve if the reciprocal of the bound on
* elements of X is not too small.
*
CALL DTBSV( UPLO, TRANS, DIAG, N, KD, AB, LDAB, X, 1 )
ELSE
*
* Use a Level 1 BLAS solve, scaling intermediate results.
*
IF( XMAX.GT.BIGNUM ) THEN
*
* Scale X so that its components are less than or equal to
* BIGNUM in absolute value.
*
SCALE = BIGNUM / XMAX
CALL DSCAL( N, SCALE, X, 1 )
XMAX = BIGNUM
END IF
*
IF( NOTRAN ) THEN
*
* Solve A * x = b
*
DO 110 J = JFIRST, JLAST, JINC
*
* Compute x(j) = b(j) / A(j,j), scaling x if necessary.
*
XJ = ABS( X( J ) )
IF( NOUNIT ) THEN
TJJS = AB( MAIND, J )*TSCAL
ELSE
TJJS = TSCAL
IF( TSCAL.EQ.ONE )
$ GO TO 100
END IF
TJJ = ABS( TJJS )
IF( TJJ.GT.SMLNUM ) THEN
*
* abs(A(j,j)) > SMLNUM:
*
IF( TJJ.LT.ONE ) THEN
IF( XJ.GT.TJJ*BIGNUM ) THEN
*
* Scale x by 1/b(j).
*
REC = ONE / XJ
CALL DSCAL( N, REC, X, 1 )
SCALE = SCALE*REC
XMAX = XMAX*REC
END IF
END IF
X( J ) = X( J ) / TJJS
XJ = ABS( X( J ) )
ELSE IF( TJJ.GT.ZERO ) THEN
*
* 0 < abs(A(j,j)) <= SMLNUM:
*
IF( XJ.GT.TJJ*BIGNUM ) THEN
*
* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM
* to avoid overflow when dividing by A(j,j).
*
REC = ( TJJ*BIGNUM ) / XJ
IF( CNORM( J ).GT.ONE ) THEN
*
* Scale by 1/CNORM(j) to avoid overflow when
* multiplying x(j) times column j.
*
REC = REC / CNORM( J )
END IF
CALL DSCAL( N, REC, X, 1 )
SCALE = SCALE*REC
XMAX = XMAX*REC
END IF
X( J ) = X( J ) / TJJS
XJ = ABS( X( J ) )
ELSE
*
* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
* scale = 0, and compute a solution to A*x = 0.
*
DO 90 I = 1, N
X( I ) = ZERO
90 CONTINUE
X( J ) = ONE
XJ = ONE
SCALE = ZERO
XMAX = ZERO
END IF
100 CONTINUE
*
* Scale x if necessary to avoid overflow when adding a
* multiple of column j of A.
*
IF( XJ.GT.ONE ) THEN
REC = ONE / XJ
IF( CNORM( J ).GT.( BIGNUM-XMAX )*REC ) THEN
*
* Scale x by 1/(2*abs(x(j))).
*
REC = REC*HALF
CALL DSCAL( N, REC, X, 1 )
SCALE = SCALE*REC
END IF
ELSE IF( XJ*CNORM( J ).GT.( BIGNUM-XMAX ) ) THEN
*
* Scale x by 1/2.
*
CALL DSCAL( N, HALF, X, 1 )
SCALE = SCALE*HALF
END IF
*
IF( UPPER ) THEN
IF( J.GT.1 ) THEN
*
* Compute the update
* x(max(1,j-kd):j-1) := x(max(1,j-kd):j-1) -
* x(j)* A(max(1,j-kd):j-1,j)
*
JLEN = MIN( KD, J-1 )
CALL DAXPY( JLEN, -X( J )*TSCAL,
$ AB( KD+1-JLEN, J ), 1, X( J-JLEN ), 1 )
I = IDAMAX( J-1, X, 1 )
XMAX = ABS( X( I ) )
END IF
ELSE IF( J.LT.N ) THEN
*
* Compute the update
* x(j+1:min(j+kd,n)) := x(j+1:min(j+kd,n)) -
* x(j) * A(j+1:min(j+kd,n),j)
*
JLEN = MIN( KD, N-J )
IF( JLEN.GT.0 )
$ CALL DAXPY( JLEN, -X( J )*TSCAL, AB( 2, J ), 1,
$ X( J+1 ), 1 )
I = J + IDAMAX( N-J, X( J+1 ), 1 )
XMAX = ABS( X( I ) )
END IF
110 CONTINUE
*
ELSE
*
* Solve A' * x = b
*
DO 160 J = JFIRST, JLAST, JINC
*
* Compute x(j) = b(j) - sum A(k,j)*x(k).
* k<>j
*
XJ = ABS( X( J ) )
USCAL = TSCAL
REC = ONE / MAX( XMAX, ONE )
IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN
*
* If x(j) could overflow, scale x by 1/(2*XMAX).
*
REC = REC*HALF
IF( NOUNIT ) THEN
TJJS = AB( MAIND, J )*TSCAL
ELSE
TJJS = TSCAL
END IF
TJJ = ABS( TJJS )
IF( TJJ.GT.ONE ) THEN
*
* Divide by A(j,j) when scaling x if A(j,j) > 1.
*
REC = MIN( ONE, REC*TJJ )
USCAL = USCAL / TJJS
END IF
IF( REC.LT.ONE ) THEN
CALL DSCAL( N, REC, X, 1 )
SCALE = SCALE*REC
XMAX = XMAX*REC
END IF
END IF
*
SUMJ = ZERO
IF( USCAL.EQ.ONE ) THEN
*
* If the scaling needed for A in the dot product is 1,
* call DDOT to perform the dot product.
*
IF( UPPER ) THEN
JLEN = MIN( KD, J-1 )
SUMJ = DDOT( JLEN, AB( KD+1-JLEN, J ), 1,
$ X( J-JLEN ), 1 )
ELSE
JLEN = MIN( KD, N-J )
IF( JLEN.GT.0 )
$ SUMJ = DDOT( JLEN, AB( 2, J ), 1, X( J+1 ), 1 )
END IF
ELSE
*
* Otherwise, use in-line code for the dot product.
*
IF( UPPER ) THEN
JLEN = MIN( KD, J-1 )
DO 120 I = 1, JLEN
SUMJ = SUMJ + ( AB( KD+I-JLEN, J )*USCAL )*
$ X( J-JLEN-1+I )
120 CONTINUE
ELSE
JLEN = MIN( KD, N-J )
DO 130 I = 1, JLEN
SUMJ = SUMJ + ( AB( I+1, J )*USCAL )*X( J+I )
130 CONTINUE
END IF
END IF
*
IF( USCAL.EQ.TSCAL ) THEN
*
* Compute x(j) := ( x(j) - sumj ) / A(j,j) if 1/A(j,j)
* was not used to scale the dotproduct.
*
X( J ) = X( J ) - SUMJ
XJ = ABS( X( J ) )
IF( NOUNIT ) THEN
*
* Compute x(j) = x(j) / A(j,j), scaling if necessary.
*
TJJS = AB( MAIND, J )*TSCAL
ELSE
TJJS = TSCAL
IF( TSCAL.EQ.ONE )
$ GO TO 150
END IF
TJJ = ABS( TJJS )
IF( TJJ.GT.SMLNUM ) THEN
*
* abs(A(j,j)) > SMLNUM:
*
IF( TJJ.LT.ONE ) THEN
IF( XJ.GT.TJJ*BIGNUM ) THEN
*
* Scale X by 1/abs(x(j)).
*
REC = ONE / XJ
CALL DSCAL( N, REC, X, 1 )
SCALE = SCALE*REC
XMAX = XMAX*REC
END IF
END IF
X( J ) = X( J ) / TJJS
ELSE IF( TJJ.GT.ZERO ) THEN
*
* 0 < abs(A(j,j)) <= SMLNUM:
*
IF( XJ.GT.TJJ*BIGNUM ) THEN
*
* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
*
REC = ( TJJ*BIGNUM ) / XJ
CALL DSCAL( N, REC, X, 1 )
SCALE = SCALE*REC
XMAX = XMAX*REC
END IF
X( J ) = X( J ) / TJJS
ELSE
*
* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
* scale = 0, and compute a solution to A'*x = 0.
*
DO 140 I = 1, N
X( I ) = ZERO
140 CONTINUE
X( J ) = ONE
SCALE = ZERO
XMAX = ZERO
END IF
150 CONTINUE
ELSE
*
* Compute x(j) := x(j) / A(j,j) - sumj if the dot
* product has already been divided by 1/A(j,j).
*
X( J ) = X( J ) / TJJS - SUMJ
END IF
XMAX = MAX( XMAX, ABS( X( J ) ) )
160 CONTINUE
END IF
SCALE = SCALE / TSCAL
END IF
*
* Scale the column norms by 1/TSCAL for return.
*
IF( TSCAL.NE.ONE ) THEN
CALL DSCAL( N, ONE / TSCAL, CNORM, 1 )
END IF
*
RETURN
*
* End of DLATBS
*
END

2
test_problems/cathermo/HMW_test_1/output_blessed.txt
View file

@ -1,7 +1,7 @@
Index Name MoleF MolalityCropped Charge
0 H2O(L) 8.1992706e-01 5.5508435e+01 0.0
1 Cl- 9.0036447e-02 6.0953986e+00 -1.0
2 H+ 3.1947185e-11 2.1628000e-09 1.0
2 H+ 3.1947189e-11 2.1628000e-09 1.0
3 Na+ 9.0036468e-02 6.0954000e+00 1.0
4 OH- 2.0645728e-08 1.3977000e-06 -1.0

274
test_problems/cathermo/HMW_test_1/output_noD_blessed.txt
View file

@ -1,7 +1,7 @@
Index Name MoleF MolalityCropped Charge
0 H2O(L) 8.1992706e-01 5.5508435e+01 0.0
1 Cl- 9.0036447e-02 6.0953986e+00 -1.0
2 H+ 3.1947185e-11 2.1628000e-09 1.0
2 H+ 3.1947189e-11 2.1628000e-09 1.0
3 Na+ 9.0036468e-02 6.0954000e+00 1.0
4 OH- 2.0645728e-08 1.3977000e-06 -1.0
@ -28,6 +28,278 @@ Index Name MoleF MolalityCropped Charge
OH- Na+ Cl- -0.00600
a1 = 3.04284e-10
a2 = 3.04284e-10
Debugging information from hmw_act
Step 1:
ionic strenth = 6.0997000e+00
total molar charge = 1.2199400e+01
Is = 6.0997
ij = 1, elambda = 0.0454012, elambda1 = -0.00306854
ij = 2, elambda = 0.200776, elambda1 = -0.014532
ij = 3, elambda = 0.47109, elambda1 = -0.0351127
ij = 4, elambda = 0.857674, elambda1 = -0.0650149
ij = 4, elambda = 0.857674, elambda1 = -0.0650149
ij = 6, elambda = 1.98206, elambda1 = -0.153152
ij = 8, elambda = 3.57685, elambda1 = -0.279391
ij = 9, elambda = 4.55112, elambda1 = -0.356872
ij = 12, elambda = 8.18289, elambda1 = -0.646977
ij = 16, elambda = 14.6822, elambda1 = -1.16875
Step 2:
z1= 1 z2= 1 E-theta(I) = 0.000000, E-thetaprime(I) = 0.000000
z1= 1 z2= 2 E-theta(I) = -0.059044, E-thetaprime(I) = 0.004790
z1= 1 z2= 3 E-theta(I) = -0.355533, E-thetaprime(I) = 0.028969
z1= 1 z2= 4 E-theta(I) = -1.068400, E-thetaprime(I) = 0.087216
z1= 2 z2= 1 E-theta(I) = -0.059044, E-thetaprime(I) = 0.004790
z1= 2 z2= 2 E-theta(I) = 0.000000, E-thetaprime(I) = 0.000000
z1= 2 z2= 3 E-theta(I) = -0.178237, E-thetaprime(I) = 0.014566
z1= 2 z2= 4 E-theta(I) = -0.951372, E-thetaprime(I) = 0.077813
z1= 3 z2= 1 E-theta(I) = -0.355533, E-thetaprime(I) = 0.028969
z1= 3 z2= 2 E-theta(I) = -0.178237, E-thetaprime(I) = 0.014566
z1= 3 z2= 3 E-theta(I) = 0.000000, E-thetaprime(I) = 0.000000
z1= 3 z2= 4 E-theta(I) = -0.357010, E-thetaprime(I) = 0.029220
z1= 4 z2= 1 E-theta(I) = -1.068400, E-thetaprime(I) = 0.087216
z1= 4 z2= 2 E-theta(I) = -0.951372, E-thetaprime(I) = 0.077813
z1= 4 z2= 3 E-theta(I) = -0.357010, E-thetaprime(I) = 0.029220
z1= 4 z2= 4 E-theta(I) = 0.000000, E-thetaprime(I) = 0.000000
Step 3:
Species Species g(x) hfunc(x)
Cl- H+ 0.07849 -0.07133
Cl- Na+ 0.07849 -0.07133
Cl- OH- 0.00000 0.00000
H+ Na+ 0.00000 0.00000
H+ OH- 0.07849 -0.07133
Na+ OH- 0.07849 -0.07133
Step 4:
Species Species BMX BprimeMX BphiMX
1 0.200614: 0.1775 0.2945 0 0.0784862
Cl- H+ 0.2006142 -0.0034438 0.1796081
2 0.0974087: 0.0765 0.2664 0 0.0784862
Cl- Na+ 0.0974087 -0.0031152 0.0784069
Cl- OH- 0.0000000 0.0000000 0.0000000
H+ Na+ 0.0000000 0.0000000 0.0000000
5 0: 0 0 0 0.0784862
H+ OH- 0.0000000 0.0000000 0.0000000
6 0.106257: 0.0864 0.253 0 0.0784862
Na+ OH- 0.1062570 -0.0029585 0.0882110
Step 5:
Species Species CMX
Cl- H+ 0.0004000
Cl- Na+ 0.0006350
Cl- OH- 0.0000000
H+ Na+ 0.0000000
H+ OH- 0.0000000
Na+ OH- 0.0022000
Step 6:
Species Species Phi_ij Phiprime_ij Phi^phi_ij
Cl- H+ 0.000000 0.000000 0.000000
Cl- Na+ 0.000000 0.000000 0.000000
Cl- OH- -0.050000 0.000000 -0.050000
H+ Na+ 0.036000 0.000000 0.036000
H+ OH- 0.000000 0.000000 0.000000
Na+ OH- 0.000000 0.000000 0.000000
Step 7:
initial value of F = -1.143942
F = -1.143942
F = -1.259847
F = -1.259847
F = -1.259847
F = -1.259847
F = -1.259847
Step 8: Summing in All Contributions to Activity Coefficients
Contributions to ln(ActCoeff_Cl-):
Unary term: z*z*F = -1.25985
Tern CMX term on Cl-,H+: abs(z_i) m_j m_k CMX = 0.00000
Tern CMX term on Cl-,Na+: abs(z_i) m_j m_k CMX = 0.02363
Bin term with H+: 2 m_j BMX = 0.00000
m_j Z CMX = 0.00000
Psi term on H+,Na+: m_j m_k psi_ijk = -0.00000
Bin term with Na+: 2 m_j BMX = 1.18833
m_j Z CMX = 0.04725
Phi term with OH-: 2 m_j Phi_aa = -0.00000
Psi term on OH-,Na+: m_j m_k psi_ijk = -0.00000
Tern CMX term on OH-,Na+: abs(z_i) m_j m_k CMX = 0.00000
Net Cl- lngamma[i] = -0.00064 gamma[i]= 0.999359
Contributions to ln(ActCoeff_H+):
Unary term: z*z*F = -1.25985
Bin term with Cl-: 2 m_j BMX = 2.44737
m_j Z CMX = 0.02977
Tern CMX term on H+,Cl-: abs(z_i) m_j m_k CMX = 0.00000
Phi term with Na+: 2 m_j Phi_cc = 0.43918
Psi term on Na+,Cl-: m_j m_k psi_ijk = -0.14883
Tern CMX term on Na+,Cl-: abs(z_i) m_j m_k CMX = 0.02363
Tern CMX term on Na+,OH-: abs(z_i) m_j m_k CMX = 0.00000
Bin term with OH-: 2 m_j BMX = 0.00000
m_j Z CMX = 0.00000
Net H+ lngamma[i] = 1.53127 gamma[i]= 4.624042
Contributions to ln(ActCoeff_Na+):
Unary term: z*z*F = -1.25985
Bin term with Cl-: 2 m_j BMX = 1.18833
m_j Z CMX = 0.04725
Psi term on Cl-,OH-: m_j m_k psi_ijk = -0.00000
Phi term with H+: 2 m_j Phi_cc = 0.00000
Psi term on H+,Cl-: m_j m_k psi_ijk = -0.00000
Tern CMX term on H+,Cl-: abs(z_i) m_j m_k CMX = 0.00000
Tern CMX term on Na+,Cl-: abs(z_i) m_j m_k CMX = 0.02363
Tern CMX term on Na+,OH-: abs(z_i) m_j m_k CMX = 0.00000
Bin term with OH-: 2 m_j BMX = 0.00000
m_j Z CMX = 0.00000
Net Na+ lngamma[i] = -0.00064 gamma[i]= 0.999359
Contributions to ln(ActCoeff_OH-):
Unary term: z*z*F = -1.25985
Phi term with Cl-: 2 m_j Phi_aa = -0.60997
Tern CMX term on Cl-,H+: abs(z_i) m_j m_k CMX = 0.00000
Psi term on Cl-,Na+: m_j m_k psi_ijk = -0.22324
Tern CMX term on Cl-,Na+: abs(z_i) m_j m_k CMX = 0.02363
Bin term with H+: 2 m_j BMX = 0.00000
m_j Z CMX = 0.00000
Bin term with Na+: 2 m_j BMX = 1.29627
m_j Z CMX = 0.16371
Tern CMX term on OH-,Na+: abs(z_i) m_j m_k CMX = 0.00000
Net OH- lngamma[i] = -0.60945 gamma[i]= 0.543650
Step 9:
term1= -1.489777 sum1= 3.205458 sum2= 0.000000 sum3= -0.000001 sum4= 0.000000 sum5= 0.000000
sum_m_phi_minus_1= 3.431360 osmotic_coef= 1.281273
Step 10:
Weight of Solvent = 18.01528
molalitySumUncropped = 12.1994
ln_a_water= -0.281593 a_water= 0.754581
Debugging information from hmw_act
Step 1:
ionic strenth = 6.0997000e+00
total molar charge = 1.2199400e+01
Is = 6.0997
ij = 1, elambda = 0.0454012, elambda1 = -0.00306854
ij = 2, elambda = 0.200776, elambda1 = -0.014532
ij = 3, elambda = 0.47109, elambda1 = -0.0351127
ij = 4, elambda = 0.857674, elambda1 = -0.0650149
ij = 4, elambda = 0.857674, elambda1 = -0.0650149
ij = 6, elambda = 1.98206, elambda1 = -0.153152
ij = 8, elambda = 3.57685, elambda1 = -0.279391
ij = 9, elambda = 4.55112, elambda1 = -0.356872
ij = 12, elambda = 8.18289, elambda1 = -0.646977
ij = 16, elambda = 14.6822, elambda1 = -1.16875
Step 2:
z1= 1 z2= 1 E-theta(I) = 0.000000, E-thetaprime(I) = 0.000000
z1= 1 z2= 2 E-theta(I) = -0.059044, E-thetaprime(I) = 0.004790
z1= 1 z2= 3 E-theta(I) = -0.355533, E-thetaprime(I) = 0.028969
z1= 1 z2= 4 E-theta(I) = -1.068400, E-thetaprime(I) = 0.087216
z1= 2 z2= 1 E-theta(I) = -0.059044, E-thetaprime(I) = 0.004790
z1= 2 z2= 2 E-theta(I) = 0.000000, E-thetaprime(I) = 0.000000
z1= 2 z2= 3 E-theta(I) = -0.178237, E-thetaprime(I) = 0.014566
z1= 2 z2= 4 E-theta(I) = -0.951372, E-thetaprime(I) = 0.077813
z1= 3 z2= 1 E-theta(I) = -0.355533, E-thetaprime(I) = 0.028969
z1= 3 z2= 2 E-theta(I) = -0.178237, E-thetaprime(I) = 0.014566
z1= 3 z2= 3 E-theta(I) = 0.000000, E-thetaprime(I) = 0.000000
z1= 3 z2= 4 E-theta(I) = -0.357010, E-thetaprime(I) = 0.029220
z1= 4 z2= 1 E-theta(I) = -1.068400, E-thetaprime(I) = 0.087216
z1= 4 z2= 2 E-theta(I) = -0.951372, E-thetaprime(I) = 0.077813
z1= 4 z2= 3 E-theta(I) = -0.357010, E-thetaprime(I) = 0.029220
z1= 4 z2= 4 E-theta(I) = 0.000000, E-thetaprime(I) = 0.000000
Step 3:
Species Species g(x) hfunc(x)
Cl- H+ 0.07849 -0.07133
Cl- Na+ 0.07849 -0.07133
Cl- OH- 0.00000 0.00000
H+ Na+ 0.00000 0.00000
H+ OH- 0.07849 -0.07133
Na+ OH- 0.07849 -0.07133
Step 4:
Species Species BMX BprimeMX BphiMX
1 0.200614: 0.1775 0.2945 0 0.0784862
Cl- H+ 0.2006142 -0.0034438 0.1796081
2 0.0974087: 0.0765 0.2664 0 0.0784862
Cl- Na+ 0.0974087 -0.0031152 0.0784069
Cl- OH- 0.0000000 0.0000000 0.0000000
H+ Na+ 0.0000000 0.0000000 0.0000000
5 0: 0 0 0 0.0784862
H+ OH- 0.0000000 0.0000000 0.0000000
6 0.106257: 0.0864 0.253 0 0.0784862
Na+ OH- 0.1062570 -0.0029585 0.0882110
Step 5:
Species Species CMX
Cl- H+ 0.0004000
Cl- Na+ 0.0006350
Cl- OH- 0.0000000
H+ Na+ 0.0000000
H+ OH- 0.0000000
Na+ OH- 0.0022000
Step 6:
Species Species Phi_ij Phiprime_ij Phi^phi_ij
Cl- H+ 0.000000 0.000000 0.000000
Cl- Na+ 0.000000 0.000000 0.000000
Cl- OH- -0.050000 0.000000 -0.050000
H+ Na+ 0.036000 0.000000 0.036000
H+ OH- 0.000000 0.000000 0.000000
Na+ OH- 0.000000 0.000000 0.000000
Step 7:
initial value of F = -1.143942
F = -1.143942
F = -1.259847
F = -1.259847
F = -1.259847
F = -1.259847
F = -1.259847
Step 8: Summing in All Contributions to Activity Coefficients
Contributions to ln(ActCoeff_Cl-):
Unary term: z*z*F = -1.25985
Tern CMX term on Cl-,H+: abs(z_i) m_j m_k CMX = 0.00000
Tern CMX term on Cl-,Na+: abs(z_i) m_j m_k CMX = 0.02363
Bin term with H+: 2 m_j BMX = 0.00000
m_j Z CMX = 0.00000
Psi term on H+,Na+: m_j m_k psi_ijk = -0.00000
Bin term with Na+: 2 m_j BMX = 1.18833
m_j Z CMX = 0.04725
Phi term with OH-: 2 m_j Phi_aa = -0.00000
Psi term on OH-,Na+: m_j m_k psi_ijk = -0.00000
Tern CMX term on OH-,Na+: abs(z_i) m_j m_k CMX = 0.00000
Net Cl- lngamma[i] = -0.00064 gamma[i]= 0.999359
Contributions to ln(ActCoeff_H+):
Unary term: z*z*F = -1.25985
Bin term with Cl-: 2 m_j BMX = 2.44737
m_j Z CMX = 0.02977
Tern CMX term on H+,Cl-: abs(z_i) m_j m_k CMX = 0.00000
Phi term with Na+: 2 m_j Phi_cc = 0.43918
Psi term on Na+,Cl-: m_j m_k psi_ijk = -0.14883
Tern CMX term on Na+,Cl-: abs(z_i) m_j m_k CMX = 0.02363
Tern CMX term on Na+,OH-: abs(z_i) m_j m_k CMX = 0.00000
Bin term with OH-: 2 m_j BMX = 0.00000
m_j Z CMX = 0.00000
Net H+ lngamma[i] = 1.53127 gamma[i]= 4.624042
Contributions to ln(ActCoeff_Na+):
Unary term: z*z*F = -1.25985
Bin term with Cl-: 2 m_j BMX = 1.18833
m_j Z CMX = 0.04725
Psi term on Cl-,OH-: m_j m_k psi_ijk = -0.00000
Phi term with H+: 2 m_j Phi_cc = 0.00000
Psi term on H+,Cl-: m_j m_k psi_ijk = -0.00000
Tern CMX term on H+,Cl-: abs(z_i) m_j m_k CMX = 0.00000
Tern CMX term on Na+,Cl-: abs(z_i) m_j m_k CMX = 0.02363
Tern CMX term on Na+,OH-: abs(z_i) m_j m_k CMX = 0.00000
Bin term with OH-: 2 m_j BMX = 0.00000
m_j Z CMX = 0.00000
Net Na+ lngamma[i] = -0.00064 gamma[i]= 0.999359
Contributions to ln(ActCoeff_OH-):
Unary term: z*z*F = -1.25985
Phi term with Cl-: 2 m_j Phi_aa = -0.60997
Tern CMX term on Cl-,H+: abs(z_i) m_j m_k CMX = 0.00000
Psi term on Cl-,Na+: m_j m_k psi_ijk = -0.22324
Tern CMX term on Cl-,Na+: abs(z_i) m_j m_k CMX = 0.02363
Bin term with H+: 2 m_j BMX = 0.00000
m_j Z CMX = 0.00000
Bin term with Na+: 2 m_j BMX = 1.29627
m_j Z CMX = 0.16371
Tern CMX term on OH-,Na+: abs(z_i) m_j m_k CMX = 0.00000
Net OH- lngamma[i] = -0.60945 gamma[i]= 0.543650
Step 9:
term1= -1.489777 sum1= 3.205458 sum2= 0.000000 sum3= -0.000001 sum4= 0.000000 sum5= 0.000000
sum_m_phi_minus_1= 3.431360 osmotic_coef= 1.281273
Step 10:
Weight of Solvent = 18.01528
molalitySumUncropped = 12.1994
ln_a_water= -0.281593 a_water= 0.754581
Name Activity ActCoeffMolal MoleFract Molality
H2O(L) 0.754581 0.92042 0.819823 55.5084
Cl- 6.09579 0.999359 0.0900885 6.0997