diff --git a/Cantera/fortran/f77demos/Makefile.in b/Cantera/fortran/f77demos/Makefile.in index d99410041..8ba8599e5 100644 --- a/Cantera/fortran/f77demos/Makefile.in +++ b/Cantera/fortran/f77demos/Makefile.in @@ -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@ diff --git a/Cantera/src/numerics/BEulerInt.cpp b/Cantera/src/numerics/BEulerInt.cpp index 78753267d..a2abc55fb 100644 --- a/Cantera/src/numerics/BEulerInt.cpp +++ b/Cantera/src/numerics/BEulerInt.cpp @@ -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 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, diff --git a/Cantera/src/numerics/BEulerInt.h b/Cantera/src/numerics/BEulerInt.h index 7481558fc..3a02a1368 100644 --- a/Cantera/src/numerics/BEulerInt.h +++ b/Cantera/src/numerics/BEulerInt.h @@ -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. diff --git a/Cantera/src/numerics/BandMatrix.cpp b/Cantera/src/numerics/BandMatrix.cpp index 0571df81c..44a045dfb 100644 --- a/Cantera/src/numerics/BandMatrix.cpp +++ b/Cantera/src/numerics/BandMatrix.cpp @@ -17,6 +17,7 @@ #include "ctexceptions.h" #include "stringUtils.h" #include "global.h" +#include 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(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(&y); + (void) memcpy(DATA_PTR(data), yyPtr->ptrColumn(0), n); + } + //==================================================================================================================== + /* + * clear the factored flag + */ + void BandMatrix::clearFactorFlag() { + m_factored = 0; + } + //==================================================================================================================== + //==================================================================================================================== } diff --git a/Cantera/src/numerics/BandMatrix.h b/Cantera/src/numerics/BandMatrix.h index 148029c95..6a756172e 100644 --- a/Cantera/src/numerics/BandMatrix.h +++ b/Cantera/src/numerics/BandMatrix.h @@ -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 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 diff --git a/Cantera/src/numerics/DenseMatrix.cpp b/Cantera/src/numerics/DenseMatrix.cpp index c56f242fb..134acfe72 100644 --- a/Cantera/src/numerics/DenseMatrix.cpp +++ b/Cantera/src/numerics/DenseMatrix.cpp @@ -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) diff --git a/Cantera/src/numerics/DenseMatrix.h b/Cantera/src/numerics/DenseMatrix.h index b6d37af1c..c43018d4a 100644 --- a/Cantera/src/numerics/DenseMatrix.h +++ b/Cantera/src/numerics/DenseMatrix.h @@ -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 /*! diff --git a/Cantera/src/numerics/GeneralMatrix.cpp b/Cantera/src/numerics/GeneralMatrix.cpp new file mode 100644 index 000000000..5a76c5e6b --- /dev/null +++ b/Cantera/src/numerics/GeneralMatrix.cpp @@ -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() + { + } + //==================================================================================================================== +} diff --git a/Cantera/src/numerics/GeneralMatrix.h b/Cantera/src/numerics/GeneralMatrix.h new file mode 100644 index 000000000..8781ff524 --- /dev/null +++ b/Cantera/src/numerics/GeneralMatrix.h @@ -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 diff --git a/Cantera/src/numerics/Makefile.in b/Cantera/src/numerics/Makefile.in index 8787bb236..e83028c68 100644 --- a/Cantera/src/numerics/Makefile.in +++ b/Cantera/src/numerics/Makefile.in @@ -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 diff --git a/Cantera/src/numerics/NonlinearSolver.cpp b/Cantera/src/numerics/NonlinearSolver.cpp index 962389c64..81c92ec00 100644 --- a/Cantera/src/numerics/NonlinearSolver.cpp +++ b/Cantera/src/numerics/NonlinearSolver.cpp @@ -19,6 +19,7 @@ #include #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 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 & 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; } diff --git a/Cantera/src/numerics/NonlinearSolver.h b/Cantera/src/numerics/NonlinearSolver.h index 582c3fa6c..6972934ac 100644 --- a/Cantera/src/numerics/NonlinearSolver.h +++ b/Cantera/src/numerics/NonlinearSolver.h @@ -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 & 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 diff --git a/Cantera/src/numerics/ResidJacEval.cpp b/Cantera/src/numerics/ResidJacEval.cpp index 62ca74b79..a4fc0ac22 100644 --- a/Cantera/src/numerics/ResidJacEval.cpp +++ b/Cantera/src/numerics/ResidJacEval.cpp @@ -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:: diff --git a/Cantera/src/numerics/ResidJacEval.h b/Cantera/src/numerics/ResidJacEval.h index 9514a1635..ccbaf4b98 100644 --- a/Cantera/src/numerics/ResidJacEval.h +++ b/Cantera/src/numerics/ResidJacEval.h @@ -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. /*! diff --git a/Cantera/src/numerics/SquareMatrix.cpp b/Cantera/src/numerics/SquareMatrix.cpp index 1724f8c93..29e16bef0 100644 --- a/Cantera/src/numerics/SquareMatrix.cpp +++ b/Cantera/src/numerics/SquareMatrix.cpp @@ -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(& 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(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; + } + //===================================================================================================================== + } diff --git a/Cantera/src/numerics/SquareMatrix.h b/Cantera/src/numerics/SquareMatrix.h index b95d3b35b..0eb03dbe8 100644 --- a/Cantera/src/numerics/SquareMatrix.h +++ b/Cantera/src/numerics/SquareMatrix.h @@ -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 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 - - diff --git a/Cantera/src/numerics/ctlapack.h b/Cantera/src/numerics/ctlapack.h index ea7083c1f..b1fae29b0 100644 --- a/Cantera/src/numerics/ctlapack.h +++ b/Cantera/src/numerics/ctlapack.h @@ -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; diff --git a/Cantera/src/oneD/MultiNewton.cpp b/Cantera/src/oneD/MultiNewton.cpp index e85f4da54..1f79e125f 100644 --- a/Cantera/src/oneD/MultiNewton.cpp +++ b/Cantera/src/oneD/MultiNewton.cpp @@ -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) { diff --git a/Cantera/src/thermo/Makefile.in b/Cantera/src/thermo/Makefile.in index 9a933701f..4cdc7580e 100644 --- a/Cantera/src/thermo/Makefile.in +++ b/Cantera/src/thermo/Makefile.in @@ -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) diff --git a/Cantera/src/thermo/MargulesVPSSTP.h b/Cantera/src/thermo/MargulesVPSSTP.h index 93adb082c..ebeb9d5e2 100644 --- a/Cantera/src/thermo/MargulesVPSSTP.h +++ b/Cantera/src/thermo/MargulesVPSSTP.h @@ -23,7 +23,7 @@ #ifndef CT_MARGULESVPSSTP_H #define CT_MARGULESVPSSTP_H -#include "PseudoBinaryVPSSTP.h" +#include "GibbsExcessVPSSTP.h" namespace Cantera { diff --git a/Cantera/src/thermo/MixedSolventElectrolyte.cpp b/Cantera/src/thermo/MixedSolventElectrolyte.cpp index 64d4b9bd1..fc6edbfa4 100644 --- a/Cantera/src/thermo/MixedSolventElectrolyte.cpp +++ b/Cantera/src/thermo/MixedSolventElectrolyte.cpp @@ -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); } diff --git a/Cantera/src/thermo/MixedSolventElectrolyte.h b/Cantera/src/thermo/MixedSolventElectrolyte.h index e91caaca2..7ffad1f05 100644 --- a/Cantera/src/thermo/MixedSolventElectrolyte.h +++ b/Cantera/src/thermo/MixedSolventElectrolyte.h @@ -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: diff --git a/Cantera/src/thermo/PseudoBinaryVPSSTP.cpp b/Cantera/src/thermo/MolarityIonicVPSSTP.cpp similarity index 60% rename from Cantera/src/thermo/PseudoBinaryVPSSTP.cpp rename to Cantera/src/thermo/MolarityIonicVPSSTP.cpp index f909098a6..e657aea3d 100644 --- a/Cantera/src/thermo/PseudoBinaryVPSSTP.cpp +++ b/Cantera/src/thermo/MolarityIonicVPSSTP.cpp @@ -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 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; } - - + //==================================================================================================================== } diff --git a/Cantera/src/thermo/PseudoBinaryVPSSTP.h b/Cantera/src/thermo/MolarityIonicVPSSTP.h similarity index 84% rename from Cantera/src/thermo/PseudoBinaryVPSSTP.h rename to Cantera/src/thermo/MolarityIonicVPSSTP.h index 3cc9a9d5b..dd6c339e8 100644 --- a/Cantera/src/thermo/PseudoBinaryVPSSTP.h +++ b/Cantera/src/thermo/MolarityIonicVPSSTP.h @@ -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 PBMoleFractions_; + //! Vector of cation indecises in the mixture std::vector cationList_; + + //! Number of cations in the mixture int numCationSpecies_; - std::vectoranionList_; + std::vector anionList_; int numAnionSpecies_; std::vector passThroughList_; int numPassThroughSpecies_; int neutralPBindexStart; - ThermoPhase *cationPhase_; - ThermoPhase *anionPhase_; - mutable std::vector moleFractionsTmp_; private: diff --git a/Cantera/src/thermo/PureFluidPhase.cpp b/Cantera/src/thermo/PureFluidPhase.cpp index 7f1266f9e..0c30fe03e 100644 --- a/Cantera/src/thermo/PureFluidPhase.cpp +++ b/Cantera/src/thermo/PureFluidPhase.cpp @@ -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); diff --git a/Cantera/src/thermo/ThermoFactory.cpp b/Cantera/src/thermo/ThermoFactory.cpp index 9c3a13f52..ff9480d97 100644 --- a/Cantera/src/thermo/ThermoFactory.cpp +++ b/Cantera/src/thermo/ThermoFactory.cpp @@ -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 }; /* diff --git a/Cantera/src/thermo/mix_defs.h b/Cantera/src/thermo/mix_defs.h index 3a4dec238..1aef185ea 100644 --- a/Cantera/src/thermo/mix_defs.h +++ b/Cantera/src/thermo/mix_defs.h @@ -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; diff --git a/configure b/configure index ed3d01168..eafb4a6a6 100755 --- a/configure +++ b/configure @@ -308,7 +308,7 @@ ac_includes_default="\ # include #endif" -ac_subst_vars='SHELL PATH_SEPARATOR PACKAGE_NAME PACKAGE_TARNAME PACKAGE_VERSION PACKAGE_STRING PACKAGE_BUGREPORT exec_prefix prefix program_transform_name bindir sbindir libexecdir datadir sysconfdir sharedstatedir localstatedir libdir includedir oldincludedir infodir mandir build_alias host_alias target_alias DEFS ECHO_C ECHO_N ECHO_T LIBS BITCOMPILE BITHARDWARE BITCHANGE ldemulationarg CVF_LIBDIR USE_CLIB_DLL local_inst local_python_inst python_prefix python_win_prefix ctversion homedir ct_libdir ct_bindir ct_incdir ct_incroot ct_datadir ct_demodir ct_templdir ct_tutdir ct_docdir ct_dir ct_mandir build build_cpu build_vendor build_os host host_cpu host_vendor host_os target target_cpu target_vendor target_os username ctroot buildinc buildlib buildbin MAKE GRAPHVIZDIR ARCHIVE DO_RANLIB RANLIB CXX_DEPENDS USERDIR INCL_USER_CODE CXX CXXFLAGS LDFLAGS CPPFLAGS ac_ct_CXX EXEEXT OBJEXT use_sundials CVODE_LIBS IDA_LIBS sundials_include sundials_lib_dir sundials_lib sundials_lib_dep CANTERA_DEBUG_MODE COMPILE_PURE_FLUIDS phase_object_files phase_header_files COMPILE_IDEAL_SOLUTIONS COMPILE_ELECTROLYTES NEED_CATHERMO COMPILE_KINETICS COMPILE_HETEROKIN COMPILE_RXNPATH WITH_REACTORS KERNEL KERNEL_OBJ BUILD_CK LIB_DIR COMPILE_VCSNONIDEAL COMPILE_H298MODIFY_CAPABILITY COMPILE_INTERMEDIATE_ZEROED_KINETICS BOOST_INCLUDE BOOST_LIB PURIFY build_lapack build_blas BLAS_LAPACK_LIBS BLAS_LAPACK_LINK BLAS_LAPACK_DIR build_with_f2c build_f2c_lib F2C_SYSTEMLIB BOOST_LIB_DIR LOCAL_LIB_DIRS LOCAL_LIBS LOCAL_LIBS_DEP INSTALL_LIBS_DEP RAW_LIBS_DEP CANTERA_CORE_LIBS CANTERA_CORE_LIBS_DEP CT_SHARED_LIB PYTHON_CMD BUILD_PYTHON NUMPY_INC_DIR NUMPY_HOME NUMARRAY_INC_DIR NUMARRAY_HOME CANTERA_PYTHON_HOME CVSTAG MATLAB_CMD BUILD_MATLAB BUILD_CLIB export_name PIC INSTALL_PROGRAM INSTALL_SCRIPT INSTALL_DATA CC CFLAGS ac_ct_CC CXXCPP EGREP SOEXT SHARED CXX_INCLUDES LCXX_FLAGS LCXX_END_LIBS HAVE_STRIPSYMBOLS F77 FFLAGS ac_ct_F77 FLIBS F90 BUILD_F90 F90FLAGS F90BUILDFLAGS F90LIBS LCXX_FLIBS precompile_headers OS_IS_DARWIN OS_IS_WIN OS_IS_CYGWIN SHARED_CTLIB mex_ext F77_EXT CXX_EXT OBJ_EXT EXE_EXT math_libs SO LDSHARED EXTRA_LINK TSCOMPARE_abs INSTALL_abs INSTALL_VERBOSE LIBOBJS LTLIBOBJS' +ac_subst_vars='SHELL PATH_SEPARATOR PACKAGE_NAME PACKAGE_TARNAME PACKAGE_VERSION PACKAGE_STRING PACKAGE_BUGREPORT exec_prefix prefix program_transform_name bindir sbindir libexecdir datadir sysconfdir sharedstatedir localstatedir libdir includedir oldincludedir infodir mandir build_alias host_alias target_alias DEFS ECHO_C ECHO_N ECHO_T LIBS F77 FFLAGS LDFLAGS ac_ct_F77 EXEEXT OBJEXT FLIBS BITCOMPILE BITHARDWARE BITCHANGE ldemulationarg CVF_LIBDIR USE_CLIB_DLL local_inst local_python_inst python_prefix python_win_prefix ctversion homedir ct_libdir ct_bindir ct_incdir ct_incroot ct_datadir ct_demodir ct_templdir ct_tutdir ct_docdir ct_dir ct_mandir build build_cpu build_vendor build_os host host_cpu host_vendor host_os target target_cpu target_vendor target_os username ctroot buildinc buildlib buildbin MAKE GRAPHVIZDIR ARCHIVE DO_RANLIB RANLIB CXX_DEPENDS USERDIR INCL_USER_CODE CXX CXXFLAGS CPPFLAGS ac_ct_CXX use_sundials CVODE_LIBS IDA_LIBS sundials_include sundials_lib_dir sundials_lib sundials_lib_dep CANTERA_DEBUG_MODE COMPILE_PURE_FLUIDS phase_object_files phase_header_files COMPILE_IDEAL_SOLUTIONS COMPILE_ELECTROLYTES NEED_CATHERMO COMPILE_KINETICS COMPILE_HETEROKIN COMPILE_RXNPATH WITH_REACTORS KERNEL KERNEL_OBJ BUILD_CK LIB_DIR COMPILE_VCSNONIDEAL COMPILE_H298MODIFY_CAPABILITY COMPILE_INTERMEDIATE_ZEROED_KINETICS BOOST_INCLUDE BOOST_LIB PURIFY build_lapack build_blas BLAS_LAPACK_LIBS BLAS_LAPACK_LINK BLAS_LAPACK_DIR build_with_f2c build_f2c_lib F2C_SYSTEMLIB BOOST_LIB_DIR LOCAL_LIB_DIRS LOCAL_LIBS LOCAL_LIBS_DEP INSTALL_LIBS_DEP RAW_LIBS_DEP CANTERA_CORE_LIBS CANTERA_CORE_LIBS_DEP CT_SHARED_LIB PYTHON_CMD BUILD_PYTHON NUMPY_INC_DIR NUMPY_HOME NUMARRAY_INC_DIR NUMARRAY_HOME CANTERA_PYTHON_HOME CVSTAG MATLAB_CMD BUILD_MATLAB BUILD_CLIB export_name PIC INSTALL_PROGRAM INSTALL_SCRIPT INSTALL_DATA CC CFLAGS ac_ct_CC CXXCPP EGREP SOEXT SHARED CXX_INCLUDES LCXX_FLAGS LCXX_END_LIBS HAVE_STRIPSYMBOLS FC FCFLAGS ac_ct_FC FCLIBS F90 BUILD_F90 F90FLAGS F90BUILDFLAGS F90LIBS LCXX_FLIBS precompile_headers OS_IS_DARWIN OS_IS_WIN OS_IS_CYGWIN SHARED_CTLIB mex_ext F77_EXT CXX_EXT OBJ_EXT EXE_EXT math_libs SO LDSHARED EXTRA_LINK TSCOMPARE_abs INSTALL_abs INSTALL_VERBOSE LIBOBJS LTLIBOBJS' ac_subst_files='' # Initialize some variables set by options. @@ -749,6 +749,18 @@ ac_env_target_alias_set=${target_alias+set} ac_env_target_alias_value=$target_alias ac_cv_env_target_alias_set=${target_alias+set} ac_cv_env_target_alias_value=$target_alias +ac_env_F77_set=${F77+set} +ac_env_F77_value=$F77 +ac_cv_env_F77_set=${F77+set} +ac_cv_env_F77_value=$F77 +ac_env_FFLAGS_set=${FFLAGS+set} +ac_env_FFLAGS_value=$FFLAGS +ac_cv_env_FFLAGS_set=${FFLAGS+set} +ac_cv_env_FFLAGS_value=$FFLAGS +ac_env_LDFLAGS_set=${LDFLAGS+set} +ac_env_LDFLAGS_value=$LDFLAGS +ac_cv_env_LDFLAGS_set=${LDFLAGS+set} +ac_cv_env_LDFLAGS_value=$LDFLAGS ac_env_CXX_set=${CXX+set} ac_env_CXX_value=$CXX ac_cv_env_CXX_set=${CXX+set} @@ -757,10 +769,6 @@ ac_env_CXXFLAGS_set=${CXXFLAGS+set} ac_env_CXXFLAGS_value=$CXXFLAGS ac_cv_env_CXXFLAGS_set=${CXXFLAGS+set} ac_cv_env_CXXFLAGS_value=$CXXFLAGS -ac_env_LDFLAGS_set=${LDFLAGS+set} -ac_env_LDFLAGS_value=$LDFLAGS -ac_cv_env_LDFLAGS_set=${LDFLAGS+set} -ac_cv_env_LDFLAGS_value=$LDFLAGS ac_env_CPPFLAGS_set=${CPPFLAGS+set} ac_env_CPPFLAGS_value=$CPPFLAGS ac_cv_env_CPPFLAGS_set=${CPPFLAGS+set} @@ -777,14 +785,14 @@ ac_env_CXXCPP_set=${CXXCPP+set} ac_env_CXXCPP_value=$CXXCPP ac_cv_env_CXXCPP_set=${CXXCPP+set} ac_cv_env_CXXCPP_value=$CXXCPP -ac_env_F77_set=${F77+set} -ac_env_F77_value=$F77 -ac_cv_env_F77_set=${F77+set} -ac_cv_env_F77_value=$F77 -ac_env_FFLAGS_set=${FFLAGS+set} -ac_env_FFLAGS_value=$FFLAGS -ac_cv_env_FFLAGS_set=${FFLAGS+set} -ac_cv_env_FFLAGS_value=$FFLAGS +ac_env_FC_set=${FC+set} +ac_env_FC_value=$FC +ac_cv_env_FC_set=${FC+set} +ac_cv_env_FC_value=$FC +ac_env_FCFLAGS_set=${FCFLAGS+set} +ac_env_FCFLAGS_value=$FCFLAGS +ac_cv_env_FCFLAGS_set=${FCFLAGS+set} +ac_cv_env_FCFLAGS_value=$FCFLAGS # # Report the --help message. @@ -860,17 +868,19 @@ if test -n "$ac_init_help"; then cat <<\_ACEOF Some influential environment variables: - CXX C++ compiler command - CXXFLAGS C++ compiler flags + F77 Fortran 77 compiler command + FFLAGS Fortran 77 compiler flags LDFLAGS linker flags, e.g. -L if you have libraries in a nonstandard directory + CXX C++ compiler command + CXXFLAGS C++ compiler flags CPPFLAGS C/C++ preprocessor flags, e.g. -I if you have headers in a nonstandard directory CC C compiler command CFLAGS C compiler flags CXXCPP C++ preprocessor - F77 Fortran 77 compiler command - FFLAGS Fortran 77 compiler flags + FC Fortran compiler command + FCFLAGS Fortran compiler flags Use these variables to override the choices made by `configure' or to help it to find libraries and programs with nonstandard names/locations. @@ -1303,6 +1313,726 @@ ac_compiler_gnu=$ac_cv_c_compiler_gnu + + +ac_ext=f +ac_compile='$F77 -c $FFLAGS conftest.$ac_ext >&5' +ac_link='$F77 -o conftest$ac_exeext $FFLAGS $LDFLAGS conftest.$ac_ext $LIBS >&5' +ac_compiler_gnu=$ac_cv_f77_compiler_gnu +if test -n "$ac_tool_prefix"; then + for ac_prog in g77 f77 xlf frt pgf77 fort77 fl32 af77 f90 xlf90 pgf90 epcf90 f95 fort xlf95 ifc efc pgf95 lf95 gfortran + do + # Extract the first word of "$ac_tool_prefix$ac_prog", so it can be a program name with args. +set dummy $ac_tool_prefix$ac_prog; ac_word=$2 +echo "$as_me:$LINENO: checking for $ac_word" >&5 +echo $ECHO_N "checking for $ac_word... $ECHO_C" >&6 +if test "${ac_cv_prog_F77+set}" = set; then + echo $ECHO_N "(cached) $ECHO_C" >&6 +else + if test -n "$F77"; then + ac_cv_prog_F77="$F77" # Let the user override the test. +else +as_save_IFS=$IFS; IFS=$PATH_SEPARATOR +for as_dir in $PATH +do + IFS=$as_save_IFS + test -z "$as_dir" && as_dir=. + for ac_exec_ext in '' $ac_executable_extensions; do + if $as_executable_p "$as_dir/$ac_word$ac_exec_ext"; then + ac_cv_prog_F77="$ac_tool_prefix$ac_prog" + echo "$as_me:$LINENO: found $as_dir/$ac_word$ac_exec_ext" >&5 + break 2 + fi +done +done + +fi +fi +F77=$ac_cv_prog_F77 +if test -n "$F77"; then + echo "$as_me:$LINENO: result: $F77" >&5 +echo "${ECHO_T}$F77" >&6 +else + echo "$as_me:$LINENO: result: no" >&5 +echo "${ECHO_T}no" >&6 +fi + + test -n "$F77" && break + done +fi +if test -z "$F77"; then + ac_ct_F77=$F77 + for ac_prog in g77 f77 xlf frt pgf77 fort77 fl32 af77 f90 xlf90 pgf90 epcf90 f95 fort xlf95 ifc efc pgf95 lf95 gfortran +do + # Extract the first word of "$ac_prog", so it can be a program name with args. +set dummy $ac_prog; ac_word=$2 +echo "$as_me:$LINENO: checking for $ac_word" >&5 +echo $ECHO_N "checking for $ac_word... $ECHO_C" >&6 +if test "${ac_cv_prog_ac_ct_F77+set}" = set; then + echo $ECHO_N "(cached) $ECHO_C" >&6 +else + if test -n "$ac_ct_F77"; then + ac_cv_prog_ac_ct_F77="$ac_ct_F77" # Let the user override the test. +else +as_save_IFS=$IFS; IFS=$PATH_SEPARATOR +for as_dir in $PATH +do + IFS=$as_save_IFS + test -z "$as_dir" && as_dir=. + for ac_exec_ext in '' $ac_executable_extensions; do + if $as_executable_p "$as_dir/$ac_word$ac_exec_ext"; then + ac_cv_prog_ac_ct_F77="$ac_prog" + echo "$as_me:$LINENO: found $as_dir/$ac_word$ac_exec_ext" >&5 + break 2 + fi +done +done + +fi +fi +ac_ct_F77=$ac_cv_prog_ac_ct_F77 +if test -n "$ac_ct_F77"; then + echo "$as_me:$LINENO: result: $ac_ct_F77" >&5 +echo "${ECHO_T}$ac_ct_F77" >&6 +else + echo "$as_me:$LINENO: result: no" >&5 +echo "${ECHO_T}no" >&6 +fi + + test -n "$ac_ct_F77" && break +done + + F77=$ac_ct_F77 +fi + + +# Provide some information about the compiler. +echo "$as_me:1410:" \ + "checking for Fortran 77 compiler version" >&5 +ac_compiler=`set X $ac_compile; echo $2` +{ (eval echo "$as_me:$LINENO: \"$ac_compiler --version &5\"") >&5 + (eval $ac_compiler --version &5) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } +{ (eval echo "$as_me:$LINENO: \"$ac_compiler -v &5\"") >&5 + (eval $ac_compiler -v &5) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } +{ (eval echo "$as_me:$LINENO: \"$ac_compiler -V &5\"") >&5 + (eval $ac_compiler -V &5) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } +rm -f a.out + +cat >conftest.$ac_ext <<_ACEOF + program main + + end +_ACEOF +ac_clean_files_save=$ac_clean_files +ac_clean_files="$ac_clean_files a.out a.exe b.out" +# Try to create an executable without -o first, disregard a.out. +# It will help us diagnose broken compilers, and finding out an intuition +# of exeext. +echo "$as_me:$LINENO: checking for Fortran 77 compiler default output file name" >&5 +echo $ECHO_N "checking for Fortran 77 compiler default output file name... $ECHO_C" >&6 +ac_link_default=`echo "$ac_link" | sed 's/ -o *conftest[^ ]*//'` +if { (eval echo "$as_me:$LINENO: \"$ac_link_default\"") >&5 + (eval $ac_link_default) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); }; then + # Find the output, starting from the most likely. This scheme is +# not robust to junk in `.', hence go to wildcards (a.*) only as a last +# resort. + +# Be careful to initialize this variable, since it used to be cached. +# Otherwise an old cache value of `no' led to `EXEEXT = no' in a Makefile. +ac_cv_exeext= +# b.out is created by i960 compilers. +for ac_file in a_out.exe a.exe conftest.exe a.out conftest a.* conftest.* b.out +do + test -f "$ac_file" || continue + case $ac_file in + *.$ac_ext | *.xcoff | *.tds | *.d | *.pdb | *.xSYM | *.bb | *.bbg | *.o | *.obj ) + ;; + conftest.$ac_ext ) + # This is the source file. + ;; + [ab].out ) + # We found the default executable, but exeext='' is most + # certainly right. + break;; + *.* ) + ac_cv_exeext=`expr "$ac_file" : '[^.]*\(\..*\)'` + # FIXME: I believe we export ac_cv_exeext for Libtool, + # but it would be cool to find out if it's true. Does anybody + # maintain Libtool? --akim. + export ac_cv_exeext + break;; + * ) + break;; + esac +done +else + echo "$as_me: failed program was:" >&5 +sed 's/^/| /' conftest.$ac_ext >&5 + +{ { echo "$as_me:$LINENO: error: Fortran 77 compiler cannot create executables +See \`config.log' for more details." >&5 +echo "$as_me: error: Fortran 77 compiler cannot create executables +See \`config.log' for more details." >&2;} + { (exit 77); exit 77; }; } +fi + +ac_exeext=$ac_cv_exeext +echo "$as_me:$LINENO: result: $ac_file" >&5 +echo "${ECHO_T}$ac_file" >&6 + +# Check the compiler produces executables we can run. If not, either +# the compiler is broken, or we cross compile. +echo "$as_me:$LINENO: checking whether the Fortran 77 compiler works" >&5 +echo $ECHO_N "checking whether the Fortran 77 compiler works... $ECHO_C" >&6 +# FIXME: These cross compiler hacks should be removed for Autoconf 3.0 +# If not cross compiling, check that we can run a simple program. +if test "$cross_compiling" != yes; then + if { ac_try='./$ac_file' + { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5 + (eval $ac_try) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); }; }; then + cross_compiling=no + else + if test "$cross_compiling" = maybe; then + cross_compiling=yes + else + { { echo "$as_me:$LINENO: error: cannot run Fortran 77 compiled programs. +If you meant to cross compile, use \`--host'. +See \`config.log' for more details." >&5 +echo "$as_me: error: cannot run Fortran 77 compiled programs. +If you meant to cross compile, use \`--host'. +See \`config.log' for more details." >&2;} + { (exit 1); exit 1; }; } + fi + fi +fi +echo "$as_me:$LINENO: result: yes" >&5 +echo "${ECHO_T}yes" >&6 + +rm -f a.out a.exe conftest$ac_cv_exeext b.out +ac_clean_files=$ac_clean_files_save +# Check the compiler produces executables we can run. If not, either +# the compiler is broken, or we cross compile. +echo "$as_me:$LINENO: checking whether we are cross compiling" >&5 +echo $ECHO_N "checking whether we are cross compiling... $ECHO_C" >&6 +echo "$as_me:$LINENO: result: $cross_compiling" >&5 +echo "${ECHO_T}$cross_compiling" >&6 + +echo "$as_me:$LINENO: checking for suffix of executables" >&5 +echo $ECHO_N "checking for suffix of executables... $ECHO_C" >&6 +if { (eval echo "$as_me:$LINENO: \"$ac_link\"") >&5 + (eval $ac_link) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); }; then + # If both `conftest.exe' and `conftest' are `present' (well, observable) +# catch `conftest.exe'. For instance with Cygwin, `ls conftest' will +# work properly (i.e., refer to `conftest.exe'), while it won't with +# `rm'. +for ac_file in conftest.exe conftest conftest.*; do + test -f "$ac_file" || continue + case $ac_file in + *.$ac_ext | *.xcoff | *.tds | *.d | *.pdb | *.xSYM | *.bb | *.bbg | *.o | *.obj ) ;; + *.* ) ac_cv_exeext=`expr "$ac_file" : '[^.]*\(\..*\)'` + export ac_cv_exeext + break;; + * ) break;; + esac +done +else + { { echo "$as_me:$LINENO: error: cannot compute suffix of executables: cannot compile and link +See \`config.log' for more details." >&5 +echo "$as_me: error: cannot compute suffix of executables: cannot compile and link +See \`config.log' for more details." >&2;} + { (exit 1); exit 1; }; } +fi + +rm -f conftest$ac_cv_exeext +echo "$as_me:$LINENO: result: $ac_cv_exeext" >&5 +echo "${ECHO_T}$ac_cv_exeext" >&6 + +rm -f conftest.$ac_ext +EXEEXT=$ac_cv_exeext +ac_exeext=$EXEEXT +echo "$as_me:$LINENO: checking for suffix of object files" >&5 +echo $ECHO_N "checking for suffix of object files... $ECHO_C" >&6 +if test "${ac_cv_objext+set}" = set; then + echo $ECHO_N "(cached) $ECHO_C" >&6 +else + cat >conftest.$ac_ext <<_ACEOF + program main + + end +_ACEOF +rm -f conftest.o conftest.obj +if { (eval echo "$as_me:$LINENO: \"$ac_compile\"") >&5 + (eval $ac_compile) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); }; then + for ac_file in `(ls conftest.o conftest.obj; ls conftest.*) 2>/dev/null`; do + case $ac_file in + *.$ac_ext | *.xcoff | *.tds | *.d | *.pdb | *.xSYM | *.bb | *.bbg ) ;; + *) ac_cv_objext=`expr "$ac_file" : '.*\.\(.*\)'` + break;; + esac +done +else + echo "$as_me: failed program was:" >&5 +sed 's/^/| /' conftest.$ac_ext >&5 + +{ { echo "$as_me:$LINENO: error: cannot compute suffix of object files: cannot compile +See \`config.log' for more details." >&5 +echo "$as_me: error: cannot compute suffix of object files: cannot compile +See \`config.log' for more details." >&2;} + { (exit 1); exit 1; }; } +fi + +rm -f conftest.$ac_cv_objext conftest.$ac_ext +fi +echo "$as_me:$LINENO: result: $ac_cv_objext" >&5 +echo "${ECHO_T}$ac_cv_objext" >&6 +OBJEXT=$ac_cv_objext +ac_objext=$OBJEXT +# If we don't use `.F' as extension, the preprocessor is not run on the +# input file. (Note that this only needs to work for GNU compilers.) +ac_save_ext=$ac_ext +ac_ext=F +echo "$as_me:$LINENO: checking whether we are using the GNU Fortran 77 compiler" >&5 +echo $ECHO_N "checking whether we are using the GNU Fortran 77 compiler... $ECHO_C" >&6 +if test "${ac_cv_f77_compiler_gnu+set}" = set; then + echo $ECHO_N "(cached) $ECHO_C" >&6 +else + cat >conftest.$ac_ext <<_ACEOF + program main +#ifndef __GNUC__ + choke me +#endif + + end +_ACEOF +rm -f conftest.$ac_objext +if { (eval echo "$as_me:$LINENO: \"$ac_compile\"") >&5 + (eval $ac_compile) 2>conftest.er1 + ac_status=$? + grep -v '^ *+' conftest.er1 >conftest.err + rm -f conftest.er1 + cat conftest.err >&5 + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } && + { ac_try='test -z "$ac_f77_werror_flag" + || test ! -s conftest.err' + { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5 + (eval $ac_try) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); }; } && + { ac_try='test -s conftest.$ac_objext' + { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5 + (eval $ac_try) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); }; }; then + ac_compiler_gnu=yes +else + echo "$as_me: failed program was:" >&5 +sed 's/^/| /' conftest.$ac_ext >&5 + +ac_compiler_gnu=no +fi +rm -f conftest.err conftest.$ac_objext conftest.$ac_ext +ac_cv_f77_compiler_gnu=$ac_compiler_gnu + +fi +echo "$as_me:$LINENO: result: $ac_cv_f77_compiler_gnu" >&5 +echo "${ECHO_T}$ac_cv_f77_compiler_gnu" >&6 +ac_ext=$ac_save_ext +ac_test_FFLAGS=${FFLAGS+set} +ac_save_FFLAGS=$FFLAGS +FFLAGS= +echo "$as_me:$LINENO: checking whether $F77 accepts -g" >&5 +echo $ECHO_N "checking whether $F77 accepts -g... $ECHO_C" >&6 +if test "${ac_cv_prog_f77_g+set}" = set; then + echo $ECHO_N "(cached) $ECHO_C" >&6 +else + FFLAGS=-g +cat >conftest.$ac_ext <<_ACEOF + program main + + end +_ACEOF +rm -f conftest.$ac_objext +if { (eval echo "$as_me:$LINENO: \"$ac_compile\"") >&5 + (eval $ac_compile) 2>conftest.er1 + ac_status=$? + grep -v '^ *+' conftest.er1 >conftest.err + rm -f conftest.er1 + cat conftest.err >&5 + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } && + { ac_try='test -z "$ac_f77_werror_flag" + || test ! -s conftest.err' + { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5 + (eval $ac_try) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); }; } && + { ac_try='test -s conftest.$ac_objext' + { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5 + (eval $ac_try) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); }; }; then + ac_cv_prog_f77_g=yes +else + echo "$as_me: failed program was:" >&5 +sed 's/^/| /' conftest.$ac_ext >&5 + +ac_cv_prog_f77_g=no +fi +rm -f conftest.err conftest.$ac_objext conftest.$ac_ext + +fi +echo "$as_me:$LINENO: result: $ac_cv_prog_f77_g" >&5 +echo "${ECHO_T}$ac_cv_prog_f77_g" >&6 +if test "$ac_test_FFLAGS" = set; then + FFLAGS=$ac_save_FFLAGS +elif test $ac_cv_prog_f77_g = yes; then + if test "x$ac_cv_f77_compiler_gnu" = xyes; then + FFLAGS="-g -O2" + else + FFLAGS="-g" + fi +else + if test "x$ac_cv_f77_compiler_gnu" = xyes; then + FFLAGS="-O2" + else + FFLAGS= + fi +fi + +G77=`test $ac_compiler_gnu = yes && echo yes` +ac_ext=c +ac_cpp='$CPP $CPPFLAGS' +ac_compile='$CC -c $CFLAGS $CPPFLAGS conftest.$ac_ext >&5' +ac_link='$CC -o conftest$ac_exeext $CFLAGS $CPPFLAGS $LDFLAGS conftest.$ac_ext $LIBS >&5' +ac_compiler_gnu=$ac_cv_c_compiler_gnu + + +ac_ext=f +ac_compile='$F77 -c $FFLAGS conftest.$ac_ext >&5' +ac_link='$F77 -o conftest$ac_exeext $FFLAGS $LDFLAGS conftest.$ac_ext $LIBS >&5' +ac_compiler_gnu=$ac_cv_f77_compiler_gnu +echo "$as_me:$LINENO: checking how to get verbose linking output from $F77" >&5 +echo $ECHO_N "checking how to get verbose linking output from $F77... $ECHO_C" >&6 +if test "${ac_cv_prog_f77_v+set}" = set; then + echo $ECHO_N "(cached) $ECHO_C" >&6 +else + cat >conftest.$ac_ext <<_ACEOF + program main + + end +_ACEOF +rm -f conftest.$ac_objext +if { (eval echo "$as_me:$LINENO: \"$ac_compile\"") >&5 + (eval $ac_compile) 2>conftest.er1 + ac_status=$? + grep -v '^ *+' conftest.er1 >conftest.err + rm -f conftest.er1 + cat conftest.err >&5 + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } && + { ac_try='test -z "$ac_f77_werror_flag" + || test ! -s conftest.err' + { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5 + (eval $ac_try) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); }; } && + { ac_try='test -s conftest.$ac_objext' + { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5 + (eval $ac_try) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); }; }; then + ac_cv_prog_f77_v= +# Try some options frequently used verbose output +for ac_verb in -v -verbose --verbose -V -\#\#\#; do + cat >conftest.$ac_ext <<_ACEOF + program main + + end +_ACEOF + +# Compile and link our simple test program by passing a flag (argument +# 1 to this macro) to the Fortran compiler in order to get +# "verbose" output that we can then parse for the Fortran linker +# flags. +ac_save_FFLAGS=$FFLAGS +FFLAGS="$FFLAGS $ac_verb" +(eval echo $as_me:1787: \"$ac_link\") >&5 +ac_f77_v_output=`eval $ac_link 5>&1 2>&1 | grep -v 'Driving:'` +echo "$ac_f77_v_output" >&5 +FFLAGS=$ac_save_FFLAGS + +rm -f conftest* + +# On HP/UX there is a line like: "LPATH is: /foo:/bar:/baz" where +# /foo, /bar, and /baz are search directories for the Fortran linker. +# Here, we change these into -L/foo -L/bar -L/baz (and put it first): +ac_f77_v_output="`echo $ac_f77_v_output | + grep 'LPATH is:' | + sed 's,.*LPATH is\(: *[^ ]*\).*,\1,;s,: */, -L/,g'` $ac_f77_v_output" + +case $ac_f77_v_output in + # If we are using xlf then replace all the commas with spaces. + *xlfentry*) + ac_f77_v_output=`echo $ac_f77_v_output | sed 's/,/ /g'` ;; + + # With Intel ifc, ignore the quoted -mGLOB_options_string stuff (quoted + # $LIBS confuse us, and the libraries appear later in the output anyway). + *mGLOB_options_string*) + ac_f77_v_output=`echo $ac_f77_v_output | sed 's/\"-mGLOB[^\"]*\"/ /g'` ;; + + # If we are using Cray Fortran then delete quotes. + # Use "\"" instead of '"' for font-lock-mode. + # FIXME: a more general fix for quoted arguments with spaces? + *cft90*) + ac_f77_v_output=`echo $ac_f77_v_output | sed "s/\"//g"` ;; +esac + + + # look for -l* and *.a constructs in the output + for ac_arg in $ac_f77_v_output; do + case $ac_arg in + [\\/]*.a | ?:[\\/]*.a | -[lLRu]*) + ac_cv_prog_f77_v=$ac_verb + break 2 ;; + esac + done +done +if test -z "$ac_cv_prog_f77_v"; then + { echo "$as_me:$LINENO: WARNING: cannot determine how to obtain linking information from $F77" >&5 +echo "$as_me: WARNING: cannot determine how to obtain linking information from $F77" >&2;} +fi +else + echo "$as_me: failed program was:" >&5 +sed 's/^/| /' conftest.$ac_ext >&5 + +{ echo "$as_me:$LINENO: WARNING: compilation failed" >&5 +echo "$as_me: WARNING: compilation failed" >&2;} +fi +rm -f conftest.err conftest.$ac_objext conftest.$ac_ext + +fi +echo "$as_me:$LINENO: result: $ac_cv_prog_f77_v" >&5 +echo "${ECHO_T}$ac_cv_prog_f77_v" >&6 +echo "$as_me:$LINENO: checking for Fortran libraries of $F77" >&5 +echo $ECHO_N "checking for Fortran libraries of $F77... $ECHO_C" >&6 +if test "${ac_cv_f77_libs+set}" = set; then + echo $ECHO_N "(cached) $ECHO_C" >&6 +else + if test "x$FLIBS" != "x"; then + ac_cv_f77_libs="$FLIBS" # Let the user override the test. +else + +cat >conftest.$ac_ext <<_ACEOF + program main + + end +_ACEOF + +# Compile and link our simple test program by passing a flag (argument +# 1 to this macro) to the Fortran compiler in order to get +# "verbose" output that we can then parse for the Fortran linker +# flags. +ac_save_FFLAGS=$FFLAGS +FFLAGS="$FFLAGS $ac_cv_prog_f77_v" +(eval echo $as_me:1865: \"$ac_link\") >&5 +ac_f77_v_output=`eval $ac_link 5>&1 2>&1 | grep -v 'Driving:'` +echo "$ac_f77_v_output" >&5 +FFLAGS=$ac_save_FFLAGS + +rm -f conftest* + +# On HP/UX there is a line like: "LPATH is: /foo:/bar:/baz" where +# /foo, /bar, and /baz are search directories for the Fortran linker. +# Here, we change these into -L/foo -L/bar -L/baz (and put it first): +ac_f77_v_output="`echo $ac_f77_v_output | + grep 'LPATH is:' | + sed 's,.*LPATH is\(: *[^ ]*\).*,\1,;s,: */, -L/,g'` $ac_f77_v_output" + +case $ac_f77_v_output in + # If we are using xlf then replace all the commas with spaces. + *xlfentry*) + ac_f77_v_output=`echo $ac_f77_v_output | sed 's/,/ /g'` ;; + + # With Intel ifc, ignore the quoted -mGLOB_options_string stuff (quoted + # $LIBS confuse us, and the libraries appear later in the output anyway). + *mGLOB_options_string*) + ac_f77_v_output=`echo $ac_f77_v_output | sed 's/\"-mGLOB[^\"]*\"/ /g'` ;; + + # If we are using Cray Fortran then delete quotes. + # Use "\"" instead of '"' for font-lock-mode. + # FIXME: a more general fix for quoted arguments with spaces? + *cft90*) + ac_f77_v_output=`echo $ac_f77_v_output | sed "s/\"//g"` ;; +esac + + + +ac_cv_f77_libs= + +# Save positional arguments (if any) +ac_save_positional="$@" + +set X $ac_f77_v_output +while test $# != 1; do + shift + ac_arg=$1 + case $ac_arg in + [\\/]*.a | ?:[\\/]*.a) + ac_exists=false + for ac_i in $ac_cv_f77_libs; do + if test x"$ac_arg" = x"$ac_i"; then + ac_exists=true + break + fi + done + + if test x"$ac_exists" = xtrue; then + : +else + ac_cv_f77_libs="$ac_cv_f77_libs $ac_arg" +fi + + ;; + -bI:*) + ac_exists=false + for ac_i in $ac_cv_f77_libs; do + if test x"$ac_arg" = x"$ac_i"; then + ac_exists=true + break + fi + done + + if test x"$ac_exists" = xtrue; then + : +else + if test "$ac_compiler_gnu" = yes; then + for ac_link_opt in $ac_arg; do + ac_cv_f77_libs="$ac_cv_f77_libs -Xlinker $ac_link_opt" + done +else + ac_cv_f77_libs="$ac_cv_f77_libs $ac_arg" +fi +fi + + ;; + # Ignore these flags. + -lang* | -lcrt[01].o | -lcrtbegin.o | -lc | -lgcc | -libmil | -LANG:=*) + ;; + -lkernel32) + test x"$CYGWIN" != xyes && ac_cv_f77_libs="$ac_cv_f77_libs $ac_arg" + ;; + -[LRuY]) + # These flags, when seen by themselves, take an argument. + # We remove the space between option and argument and re-iterate + # unless we find an empty arg or a new option (starting with -) + case $2 in + "" | -*);; + *) + ac_arg="$ac_arg$2" + shift; shift + set X $ac_arg "$@" + ;; + esac + ;; + -YP,*) + for ac_j in `echo $ac_arg | sed -e 's/-YP,/-L/;s/:/ -L/g'`; do + ac_exists=false + for ac_i in $ac_cv_f77_libs; do + if test x"$ac_j" = x"$ac_i"; then + ac_exists=true + break + fi + done + + if test x"$ac_exists" = xtrue; then + : +else + ac_arg="$ac_arg $ac_j" + ac_cv_f77_libs="$ac_cv_f77_libs $ac_j" +fi + + done + ;; + -[lLR]*) + ac_exists=false + for ac_i in $ac_cv_f77_libs; do + if test x"$ac_arg" = x"$ac_i"; then + ac_exists=true + break + fi + done + + if test x"$ac_exists" = xtrue; then + : +else + ac_cv_f77_libs="$ac_cv_f77_libs $ac_arg" +fi + + ;; + # Ignore everything else. + esac +done +# restore positional arguments +set X $ac_save_positional; shift + +# We only consider "LD_RUN_PATH" on Solaris systems. If this is seen, +# then we insist that the "run path" must be an absolute path (i.e. it +# must begin with a "/"). +case `(uname -sr) 2>/dev/null` in + "SunOS 5"*) + ac_ld_run_path=`echo $ac_f77_v_output | + sed -n 's,^.*LD_RUN_PATH *= *\(/[^ ]*\).*$,-R\1,p'` + test "x$ac_ld_run_path" != x && + if test "$ac_compiler_gnu" = yes; then + for ac_link_opt in $ac_ld_run_path; do + ac_cv_f77_libs="$ac_cv_f77_libs -Xlinker $ac_link_opt" + done +else + ac_cv_f77_libs="$ac_cv_f77_libs $ac_ld_run_path" +fi + ;; +esac +fi # test "x$[]_AC_LANG_PREFIX[]LIBS" = "x" + +fi +echo "$as_me:$LINENO: result: $ac_cv_f77_libs" >&5 +echo "${ECHO_T}$ac_cv_f77_libs" >&6 +FLIBS="$ac_cv_f77_libs" + + +ac_ext=c +ac_cpp='$CPP $CPPFLAGS' +ac_compile='$CC -c $CFLAGS $CPPFLAGS conftest.$ac_ext >&5' +ac_link='$CC -o conftest$ac_exeext $CFLAGS $CPPFLAGS $LDFLAGS conftest.$ac_ext $LIBS >&5' +ac_compiler_gnu=$ac_cv_c_compiler_gnu ac_config_headers="$ac_config_headers config.h" @@ -1888,207 +2618,6 @@ ac_compiler=`set X $ac_compile; echo $2` echo "$as_me:$LINENO: \$? = $ac_status" >&5 (exit $ac_status); } -cat >conftest.$ac_ext <<_ACEOF -/* confdefs.h. */ -_ACEOF -cat confdefs.h >>conftest.$ac_ext -cat >>conftest.$ac_ext <<_ACEOF -/* end confdefs.h. */ - -int -main () -{ - - ; - return 0; -} -_ACEOF -ac_clean_files_save=$ac_clean_files -ac_clean_files="$ac_clean_files a.out a.exe b.out" -# Try to create an executable without -o first, disregard a.out. -# It will help us diagnose broken compilers, and finding out an intuition -# of exeext. -echo "$as_me:$LINENO: checking for C++ compiler default output file name" >&5 -echo $ECHO_N "checking for C++ compiler default output file name... $ECHO_C" >&6 -ac_link_default=`echo "$ac_link" | sed 's/ -o *conftest[^ ]*//'` -if { (eval echo "$as_me:$LINENO: \"$ac_link_default\"") >&5 - (eval $ac_link_default) 2>&5 - ac_status=$? - echo "$as_me:$LINENO: \$? = $ac_status" >&5 - (exit $ac_status); }; then - # Find the output, starting from the most likely. This scheme is -# not robust to junk in `.', hence go to wildcards (a.*) only as a last -# resort. - -# Be careful to initialize this variable, since it used to be cached. -# Otherwise an old cache value of `no' led to `EXEEXT = no' in a Makefile. -ac_cv_exeext= -# b.out is created by i960 compilers. -for ac_file in a_out.exe a.exe conftest.exe a.out conftest a.* conftest.* b.out -do - test -f "$ac_file" || continue - case $ac_file in - *.$ac_ext | *.xcoff | *.tds | *.d | *.pdb | *.xSYM | *.bb | *.bbg | *.o | *.obj ) - ;; - conftest.$ac_ext ) - # This is the source file. - ;; - [ab].out ) - # We found the default executable, but exeext='' is most - # certainly right. - break;; - *.* ) - ac_cv_exeext=`expr "$ac_file" : '[^.]*\(\..*\)'` - # FIXME: I believe we export ac_cv_exeext for Libtool, - # but it would be cool to find out if it's true. Does anybody - # maintain Libtool? --akim. - export ac_cv_exeext - break;; - * ) - break;; - esac -done -else - echo "$as_me: failed program was:" >&5 -sed 's/^/| /' conftest.$ac_ext >&5 - -{ { echo "$as_me:$LINENO: error: C++ compiler cannot create executables -See \`config.log' for more details." >&5 -echo "$as_me: error: C++ compiler cannot create executables -See \`config.log' for more details." >&2;} - { (exit 77); exit 77; }; } -fi - -ac_exeext=$ac_cv_exeext -echo "$as_me:$LINENO: result: $ac_file" >&5 -echo "${ECHO_T}$ac_file" >&6 - -# Check the compiler produces executables we can run. If not, either -# the compiler is broken, or we cross compile. -echo "$as_me:$LINENO: checking whether the C++ compiler works" >&5 -echo $ECHO_N "checking whether the C++ compiler works... $ECHO_C" >&6 -# FIXME: These cross compiler hacks should be removed for Autoconf 3.0 -# If not cross compiling, check that we can run a simple program. -if test "$cross_compiling" != yes; then - if { ac_try='./$ac_file' - { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5 - (eval $ac_try) 2>&5 - ac_status=$? - echo "$as_me:$LINENO: \$? = $ac_status" >&5 - (exit $ac_status); }; }; then - cross_compiling=no - else - if test "$cross_compiling" = maybe; then - cross_compiling=yes - else - { { echo "$as_me:$LINENO: error: cannot run C++ compiled programs. -If you meant to cross compile, use \`--host'. -See \`config.log' for more details." >&5 -echo "$as_me: error: cannot run C++ compiled programs. -If you meant to cross compile, use \`--host'. -See \`config.log' for more details." >&2;} - { (exit 1); exit 1; }; } - fi - fi -fi -echo "$as_me:$LINENO: result: yes" >&5 -echo "${ECHO_T}yes" >&6 - -rm -f a.out a.exe conftest$ac_cv_exeext b.out -ac_clean_files=$ac_clean_files_save -# Check the compiler produces executables we can run. If not, either -# the compiler is broken, or we cross compile. -echo "$as_me:$LINENO: checking whether we are cross compiling" >&5 -echo $ECHO_N "checking whether we are cross compiling... $ECHO_C" >&6 -echo "$as_me:$LINENO: result: $cross_compiling" >&5 -echo "${ECHO_T}$cross_compiling" >&6 - -echo "$as_me:$LINENO: checking for suffix of executables" >&5 -echo $ECHO_N "checking for suffix of executables... $ECHO_C" >&6 -if { (eval echo "$as_me:$LINENO: \"$ac_link\"") >&5 - (eval $ac_link) 2>&5 - ac_status=$? - echo "$as_me:$LINENO: \$? = $ac_status" >&5 - (exit $ac_status); }; then - # If both `conftest.exe' and `conftest' are `present' (well, observable) -# catch `conftest.exe'. For instance with Cygwin, `ls conftest' will -# work properly (i.e., refer to `conftest.exe'), while it won't with -# `rm'. -for ac_file in conftest.exe conftest conftest.*; do - test -f "$ac_file" || continue - case $ac_file in - *.$ac_ext | *.xcoff | *.tds | *.d | *.pdb | *.xSYM | *.bb | *.bbg | *.o | *.obj ) ;; - *.* ) ac_cv_exeext=`expr "$ac_file" : '[^.]*\(\..*\)'` - export ac_cv_exeext - break;; - * ) break;; - esac -done -else - { { echo "$as_me:$LINENO: error: cannot compute suffix of executables: cannot compile and link -See \`config.log' for more details." >&5 -echo "$as_me: error: cannot compute suffix of executables: cannot compile and link -See \`config.log' for more details." >&2;} - { (exit 1); exit 1; }; } -fi - -rm -f conftest$ac_cv_exeext -echo "$as_me:$LINENO: result: $ac_cv_exeext" >&5 -echo "${ECHO_T}$ac_cv_exeext" >&6 - -rm -f conftest.$ac_ext -EXEEXT=$ac_cv_exeext -ac_exeext=$EXEEXT -echo "$as_me:$LINENO: checking for suffix of object files" >&5 -echo $ECHO_N "checking for suffix of object files... $ECHO_C" >&6 -if test "${ac_cv_objext+set}" = set; then - echo $ECHO_N "(cached) $ECHO_C" >&6 -else - cat >conftest.$ac_ext <<_ACEOF -/* confdefs.h. */ -_ACEOF -cat confdefs.h >>conftest.$ac_ext -cat >>conftest.$ac_ext <<_ACEOF -/* end confdefs.h. */ - -int -main () -{ - - ; - return 0; -} -_ACEOF -rm -f conftest.o conftest.obj -if { (eval echo "$as_me:$LINENO: \"$ac_compile\"") >&5 - (eval $ac_compile) 2>&5 - ac_status=$? - echo "$as_me:$LINENO: \$? = $ac_status" >&5 - (exit $ac_status); }; then - for ac_file in `(ls conftest.o conftest.obj; ls conftest.*) 2>/dev/null`; do - case $ac_file in - *.$ac_ext | *.xcoff | *.tds | *.d | *.pdb | *.xSYM | *.bb | *.bbg ) ;; - *) ac_cv_objext=`expr "$ac_file" : '.*\.\(.*\)'` - break;; - esac -done -else - echo "$as_me: failed program was:" >&5 -sed 's/^/| /' conftest.$ac_ext >&5 - -{ { echo "$as_me:$LINENO: error: cannot compute suffix of object files: cannot compile -See \`config.log' for more details." >&5 -echo "$as_me: error: cannot compute suffix of object files: cannot compile -See \`config.log' for more details." >&2;} - { (exit 1); exit 1; }; } -fi - -rm -f conftest.$ac_cv_objext conftest.$ac_ext -fi -echo "$as_me:$LINENO: result: $ac_cv_objext" >&5 -echo "${ECHO_T}$ac_cv_objext" >&6 -OBJEXT=$ac_cv_objext -ac_objext=$OBJEXT echo "$as_me:$LINENO: checking whether we are using the GNU C++ compiler" >&5 echo $ECHO_N "checking whether we are using the GNU C++ compiler... $ECHO_C" >&6 if test "${ac_cv_cxx_compiler_gnu+set}" = set; then @@ -2991,6 +3520,7 @@ fi case $ac_sys_system in Linux) NEED_F2C=1 ;; esac +NEED_F2C= # # Create a variable build_f2c_lib that determines whether @@ -3002,7 +3532,7 @@ if test -n "$NEED_F2C" ; then build_f2c_lib=1 else case $ac_sys_system in - Linux) F2C_SYSTEMLIB="-lg2c" + Linux) F2C_SYSTEMLIB="" esac fi @@ -3142,7 +3672,7 @@ then RAW_LIBS_DEP=$RAW_LIBS_DEP' 'libctf2c.a else case $ac_sys_system in - Linux) LOCAL_LIBS=$LOCAL_LIBS' '-lg2c;; + Linux) LOCAL_LIBS=$LOCAL_LIBS;; esac fi # Darwin*) LOCAL_LIBS=$LOCAL_LIBS' '-lg2c;; @@ -9975,7 +10505,11 @@ echo 'checking for a strip symbol command ... ' $HAVE_STRIPSYMBOLS # Fortran #--------------------------------------------------------------------------- -#if test x"$build_with_f2c" = "x0"; then +# +# This macro sets the substitution variable, @F77@ and @G77@ +# +echo " this is $F77" +echo "this is a test" ac_ext=f ac_compile='$F77 -c $FFLAGS conftest.$ac_ext >&5' ac_link='$F77 -o conftest$ac_exeext $FFLAGS $LDFLAGS conftest.$ac_ext $LIBS >&5' @@ -10068,7 +10602,7 @@ fi # Provide some information about the compiler. -echo "$as_me:10071:" \ +echo "$as_me:10605:" \ "checking for Fortran 77 compiler version" >&5 ac_compiler=`set X $ac_compile; echo $2` { (eval echo "$as_me:$LINENO: \"$ac_compiler --version &5\"") >&5 @@ -10222,8 +10756,7 @@ else FFLAGS=$FFLAGS' -fno-second-underscore' fi fi - - +#F77LDRCLIBS= ac_ext=f ac_compile='$F77 -c $FFLAGS conftest.$ac_ext >&5' ac_link='$F77 -o conftest$ac_exeext $FFLAGS $LDFLAGS conftest.$ac_ext $LIBS >&5' @@ -10275,7 +10808,7 @@ _ACEOF # flags. ac_save_FFLAGS=$FFLAGS FFLAGS="$FFLAGS $ac_verb" -(eval echo $as_me:10278: \"$ac_link\") >&5 +(eval echo $as_me:10811: \"$ac_link\") >&5 ac_f77_v_output=`eval $ac_link 5>&1 2>&1 | grep -v 'Driving:'` echo "$ac_f77_v_output" >&5 FFLAGS=$ac_save_FFLAGS @@ -10353,7 +10886,7 @@ _ACEOF # flags. ac_save_FFLAGS=$FFLAGS FFLAGS="$FFLAGS $ac_cv_prog_f77_v" -(eval echo $as_me:10356: \"$ac_link\") >&5 +(eval echo $as_me:10889: \"$ac_link\") >&5 ac_f77_v_output=`eval $ac_link 5>&1 2>&1 | grep -v 'Driving:'` echo "$ac_f77_v_output" >&5 FFLAGS=$ac_save_FFLAGS @@ -10526,6 +11059,550 @@ ac_link='$CXX -o conftest$ac_exeext $CXXFLAGS $CPPFLAGS $LDFLAGS conftest.$ac_ex ac_compiler_gnu=$ac_cv_cxx_compiler_gnu +echo " Macro returned with FLIBS defined as " $FLIBS + + +# +ac_ext=${FC_SRCEXT-f} +ac_compile='$FC -c $FCFLAGS $FCFLAGS_SRCEXT conftest.$ac_ext >&5' +ac_link='$FC -o conftest$ac_exeext $FCFLAGS $LDFLAGS $FCFLAGS_SRCEXT conftest.$ac_ext $LIBS >&5' +ac_compiler_gnu=$ac_cv_fc_compiler_gnu +if test -n "$ac_tool_prefix"; then + for ac_prog in f95 fort xlf95 ifc efc pgf95 lf95 gfortran f90 xlf90 pgf90 epcf90 g77 f77 xlf frt pgf77 fort77 fl32 af77 + do + # Extract the first word of "$ac_tool_prefix$ac_prog", so it can be a program name with args. +set dummy $ac_tool_prefix$ac_prog; ac_word=$2 +echo "$as_me:$LINENO: checking for $ac_word" >&5 +echo $ECHO_N "checking for $ac_word... $ECHO_C" >&6 +if test "${ac_cv_prog_FC+set}" = set; then + echo $ECHO_N "(cached) $ECHO_C" >&6 +else + if test -n "$FC"; then + ac_cv_prog_FC="$FC" # Let the user override the test. +else +as_save_IFS=$IFS; IFS=$PATH_SEPARATOR +for as_dir in $PATH +do + IFS=$as_save_IFS + test -z "$as_dir" && as_dir=. + for ac_exec_ext in '' $ac_executable_extensions; do + if $as_executable_p "$as_dir/$ac_word$ac_exec_ext"; then + ac_cv_prog_FC="$ac_tool_prefix$ac_prog" + echo "$as_me:$LINENO: found $as_dir/$ac_word$ac_exec_ext" >&5 + break 2 + fi +done +done + +fi +fi +FC=$ac_cv_prog_FC +if test -n "$FC"; then + echo "$as_me:$LINENO: result: $FC" >&5 +echo "${ECHO_T}$FC" >&6 +else + echo "$as_me:$LINENO: result: no" >&5 +echo "${ECHO_T}no" >&6 +fi + + test -n "$FC" && break + done +fi +if test -z "$FC"; then + ac_ct_FC=$FC + for ac_prog in f95 fort xlf95 ifc efc pgf95 lf95 gfortran f90 xlf90 pgf90 epcf90 g77 f77 xlf frt pgf77 fort77 fl32 af77 +do + # Extract the first word of "$ac_prog", so it can be a program name with args. +set dummy $ac_prog; ac_word=$2 +echo "$as_me:$LINENO: checking for $ac_word" >&5 +echo $ECHO_N "checking for $ac_word... $ECHO_C" >&6 +if test "${ac_cv_prog_ac_ct_FC+set}" = set; then + echo $ECHO_N "(cached) $ECHO_C" >&6 +else + if test -n "$ac_ct_FC"; then + ac_cv_prog_ac_ct_FC="$ac_ct_FC" # Let the user override the test. +else +as_save_IFS=$IFS; IFS=$PATH_SEPARATOR +for as_dir in $PATH +do + IFS=$as_save_IFS + test -z "$as_dir" && as_dir=. + for ac_exec_ext in '' $ac_executable_extensions; do + if $as_executable_p "$as_dir/$ac_word$ac_exec_ext"; then + ac_cv_prog_ac_ct_FC="$ac_prog" + echo "$as_me:$LINENO: found $as_dir/$ac_word$ac_exec_ext" >&5 + break 2 + fi +done +done + +fi +fi +ac_ct_FC=$ac_cv_prog_ac_ct_FC +if test -n "$ac_ct_FC"; then + echo "$as_me:$LINENO: result: $ac_ct_FC" >&5 +echo "${ECHO_T}$ac_ct_FC" >&6 +else + echo "$as_me:$LINENO: result: no" >&5 +echo "${ECHO_T}no" >&6 +fi + + test -n "$ac_ct_FC" && break +done + + FC=$ac_ct_FC +fi + + +# Provide some information about the compiler. +echo "$as_me:11158:" \ + "checking for Fortran compiler version" >&5 +ac_compiler=`set X $ac_compile; echo $2` +{ (eval echo "$as_me:$LINENO: \"$ac_compiler --version &5\"") >&5 + (eval $ac_compiler --version &5) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } +{ (eval echo "$as_me:$LINENO: \"$ac_compiler -v &5\"") >&5 + (eval $ac_compiler -v &5) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } +{ (eval echo "$as_me:$LINENO: \"$ac_compiler -V &5\"") >&5 + (eval $ac_compiler -V &5) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } +rm -f a.out + +# If we don't use `.F' as extension, the preprocessor is not run on the +# input file. (Note that this only needs to work for GNU compilers.) +ac_save_ext=$ac_ext +ac_ext=F +echo "$as_me:$LINENO: checking whether we are using the GNU Fortran compiler" >&5 +echo $ECHO_N "checking whether we are using the GNU Fortran compiler... $ECHO_C" >&6 +if test "${ac_cv_fc_compiler_gnu+set}" = set; then + echo $ECHO_N "(cached) $ECHO_C" >&6 +else + cat >conftest.$ac_ext <<_ACEOF + program main +#ifndef __GNUC__ + choke me +#endif + + end +_ACEOF +rm -f conftest.$ac_objext +if { (eval echo "$as_me:$LINENO: \"$ac_compile\"") >&5 + (eval $ac_compile) 2>conftest.er1 + ac_status=$? + grep -v '^ *+' conftest.er1 >conftest.err + rm -f conftest.er1 + cat conftest.err >&5 + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } && + { ac_try='test -z "$ac_fc_werror_flag" + || test ! -s conftest.err' + { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5 + (eval $ac_try) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); }; } && + { ac_try='test -s conftest.$ac_objext' + { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5 + (eval $ac_try) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); }; }; then + ac_compiler_gnu=yes +else + echo "$as_me: failed program was:" >&5 +sed 's/^/| /' conftest.$ac_ext >&5 + +ac_compiler_gnu=no +fi +rm -f conftest.err conftest.$ac_objext conftest.$ac_ext +ac_cv_fc_compiler_gnu=$ac_compiler_gnu + +fi +echo "$as_me:$LINENO: result: $ac_cv_fc_compiler_gnu" >&5 +echo "${ECHO_T}$ac_cv_fc_compiler_gnu" >&6 +ac_ext=$ac_save_ext +ac_test_FFLAGS=${FCFLAGS+set} +ac_save_FFLAGS=$FCFLAGS +FCFLAGS= +echo "$as_me:$LINENO: checking whether $FC accepts -g" >&5 +echo $ECHO_N "checking whether $FC accepts -g... $ECHO_C" >&6 +if test "${ac_cv_prog_fc_g+set}" = set; then + echo $ECHO_N "(cached) $ECHO_C" >&6 +else + FCFLAGS=-g +cat >conftest.$ac_ext <<_ACEOF + program main + + end +_ACEOF +rm -f conftest.$ac_objext +if { (eval echo "$as_me:$LINENO: \"$ac_compile\"") >&5 + (eval $ac_compile) 2>conftest.er1 + ac_status=$? + grep -v '^ *+' conftest.er1 >conftest.err + rm -f conftest.er1 + cat conftest.err >&5 + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } && + { ac_try='test -z "$ac_fc_werror_flag" + || test ! -s conftest.err' + { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5 + (eval $ac_try) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); }; } && + { ac_try='test -s conftest.$ac_objext' + { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5 + (eval $ac_try) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); }; }; then + ac_cv_prog_fc_g=yes +else + echo "$as_me: failed program was:" >&5 +sed 's/^/| /' conftest.$ac_ext >&5 + +ac_cv_prog_fc_g=no +fi +rm -f conftest.err conftest.$ac_objext conftest.$ac_ext + +fi +echo "$as_me:$LINENO: result: $ac_cv_prog_fc_g" >&5 +echo "${ECHO_T}$ac_cv_prog_fc_g" >&6 +if test "$ac_test_FFLAGS" = set; then + FCFLAGS=$ac_save_FFLAGS +elif test $ac_cv_prog_fc_g = yes; then + if test "x$ac_cv_fc_compiler_gnu" = xyes; then + FCFLAGS="-g -O2" + else + FCFLAGS="-g" + fi +else + if test "x$ac_cv_fc_compiler_gnu" = xyes; then + FCFLAGS="-O2" + else + FCFLAGS= + fi +fi + +ac_ext=cc +ac_cpp='$CXXCPP $CPPFLAGS' +ac_compile='$CXX -c $CXXFLAGS $CPPFLAGS conftest.$ac_ext >&5' +ac_link='$CXX -o conftest$ac_exeext $CXXFLAGS $CPPFLAGS $LDFLAGS conftest.$ac_ext $LIBS >&5' +ac_compiler_gnu=$ac_cv_cxx_compiler_gnu + + +ac_ext=${FC_SRCEXT-f} +ac_compile='$FC -c $FCFLAGS $FCFLAGS_SRCEXT conftest.$ac_ext >&5' +ac_link='$FC -o conftest$ac_exeext $FCFLAGS $LDFLAGS $FCFLAGS_SRCEXT conftest.$ac_ext $LIBS >&5' +ac_compiler_gnu=$ac_cv_fc_compiler_gnu +echo "$as_me:$LINENO: checking how to get verbose linking output from $FC" >&5 +echo $ECHO_N "checking how to get verbose linking output from $FC... $ECHO_C" >&6 +if test "${ac_cv_prog_fc_v+set}" = set; then + echo $ECHO_N "(cached) $ECHO_C" >&6 +else + cat >conftest.$ac_ext <<_ACEOF + program main + + end +_ACEOF +rm -f conftest.$ac_objext +if { (eval echo "$as_me:$LINENO: \"$ac_compile\"") >&5 + (eval $ac_compile) 2>conftest.er1 + ac_status=$? + grep -v '^ *+' conftest.er1 >conftest.err + rm -f conftest.er1 + cat conftest.err >&5 + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } && + { ac_try='test -z "$ac_fc_werror_flag" + || test ! -s conftest.err' + { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5 + (eval $ac_try) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); }; } && + { ac_try='test -s conftest.$ac_objext' + { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5 + (eval $ac_try) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); }; }; then + ac_cv_prog_fc_v= +# Try some options frequently used verbose output +for ac_verb in -v -verbose --verbose -V -\#\#\#; do + cat >conftest.$ac_ext <<_ACEOF + program main + + end +_ACEOF + +# Compile and link our simple test program by passing a flag (argument +# 1 to this macro) to the Fortran compiler in order to get +# "verbose" output that we can then parse for the Fortran linker +# flags. +ac_save_FFLAGS=$FCFLAGS +FCFLAGS="$FCFLAGS $ac_verb" +(eval echo $as_me:11353: \"$ac_link\") >&5 +ac_fc_v_output=`eval $ac_link 5>&1 2>&1 | grep -v 'Driving:'` +echo "$ac_fc_v_output" >&5 +FCFLAGS=$ac_save_FFLAGS + +rm -f conftest* + +# On HP/UX there is a line like: "LPATH is: /foo:/bar:/baz" where +# /foo, /bar, and /baz are search directories for the Fortran linker. +# Here, we change these into -L/foo -L/bar -L/baz (and put it first): +ac_fc_v_output="`echo $ac_fc_v_output | + grep 'LPATH is:' | + sed 's,.*LPATH is\(: *[^ ]*\).*,\1,;s,: */, -L/,g'` $ac_fc_v_output" + +case $ac_fc_v_output in + # If we are using xlf then replace all the commas with spaces. + *xlfentry*) + ac_fc_v_output=`echo $ac_fc_v_output | sed 's/,/ /g'` ;; + + # With Intel ifc, ignore the quoted -mGLOB_options_string stuff (quoted + # $LIBS confuse us, and the libraries appear later in the output anyway). + *mGLOB_options_string*) + ac_fc_v_output=`echo $ac_fc_v_output | sed 's/\"-mGLOB[^\"]*\"/ /g'` ;; + + # If we are using Cray Fortran then delete quotes. + # Use "\"" instead of '"' for font-lock-mode. + # FIXME: a more general fix for quoted arguments with spaces? + *cft90*) + ac_fc_v_output=`echo $ac_fc_v_output | sed "s/\"//g"` ;; +esac + + + # look for -l* and *.a constructs in the output + for ac_arg in $ac_fc_v_output; do + case $ac_arg in + [\\/]*.a | ?:[\\/]*.a | -[lLRu]*) + ac_cv_prog_fc_v=$ac_verb + break 2 ;; + esac + done +done +if test -z "$ac_cv_prog_fc_v"; then + { echo "$as_me:$LINENO: WARNING: cannot determine how to obtain linking information from $FC" >&5 +echo "$as_me: WARNING: cannot determine how to obtain linking information from $FC" >&2;} +fi +else + echo "$as_me: failed program was:" >&5 +sed 's/^/| /' conftest.$ac_ext >&5 + +{ echo "$as_me:$LINENO: WARNING: compilation failed" >&5 +echo "$as_me: WARNING: compilation failed" >&2;} +fi +rm -f conftest.err conftest.$ac_objext conftest.$ac_ext + +fi +echo "$as_me:$LINENO: result: $ac_cv_prog_fc_v" >&5 +echo "${ECHO_T}$ac_cv_prog_fc_v" >&6 +echo "$as_me:$LINENO: checking for Fortran libraries of $FC" >&5 +echo $ECHO_N "checking for Fortran libraries of $FC... $ECHO_C" >&6 +if test "${ac_cv_fc_libs+set}" = set; then + echo $ECHO_N "(cached) $ECHO_C" >&6 +else + if test "x$FCLIBS" != "x"; then + ac_cv_fc_libs="$FCLIBS" # Let the user override the test. +else + +cat >conftest.$ac_ext <<_ACEOF + program main + + end +_ACEOF + +# Compile and link our simple test program by passing a flag (argument +# 1 to this macro) to the Fortran compiler in order to get +# "verbose" output that we can then parse for the Fortran linker +# flags. +ac_save_FFLAGS=$FCFLAGS +FCFLAGS="$FCFLAGS $ac_cv_prog_fc_v" +(eval echo $as_me:11431: \"$ac_link\") >&5 +ac_fc_v_output=`eval $ac_link 5>&1 2>&1 | grep -v 'Driving:'` +echo "$ac_fc_v_output" >&5 +FCFLAGS=$ac_save_FFLAGS + +rm -f conftest* + +# On HP/UX there is a line like: "LPATH is: /foo:/bar:/baz" where +# /foo, /bar, and /baz are search directories for the Fortran linker. +# Here, we change these into -L/foo -L/bar -L/baz (and put it first): +ac_fc_v_output="`echo $ac_fc_v_output | + grep 'LPATH is:' | + sed 's,.*LPATH is\(: *[^ ]*\).*,\1,;s,: */, -L/,g'` $ac_fc_v_output" + +case $ac_fc_v_output in + # If we are using xlf then replace all the commas with spaces. + *xlfentry*) + ac_fc_v_output=`echo $ac_fc_v_output | sed 's/,/ /g'` ;; + + # With Intel ifc, ignore the quoted -mGLOB_options_string stuff (quoted + # $LIBS confuse us, and the libraries appear later in the output anyway). + *mGLOB_options_string*) + ac_fc_v_output=`echo $ac_fc_v_output | sed 's/\"-mGLOB[^\"]*\"/ /g'` ;; + + # If we are using Cray Fortran then delete quotes. + # Use "\"" instead of '"' for font-lock-mode. + # FIXME: a more general fix for quoted arguments with spaces? + *cft90*) + ac_fc_v_output=`echo $ac_fc_v_output | sed "s/\"//g"` ;; +esac + + + +ac_cv_fc_libs= + +# Save positional arguments (if any) +ac_save_positional="$@" + +set X $ac_fc_v_output +while test $# != 1; do + shift + ac_arg=$1 + case $ac_arg in + [\\/]*.a | ?:[\\/]*.a) + ac_exists=false + for ac_i in $ac_cv_fc_libs; do + if test x"$ac_arg" = x"$ac_i"; then + ac_exists=true + break + fi + done + + if test x"$ac_exists" = xtrue; then + : +else + ac_cv_fc_libs="$ac_cv_fc_libs $ac_arg" +fi + + ;; + -bI:*) + ac_exists=false + for ac_i in $ac_cv_fc_libs; do + if test x"$ac_arg" = x"$ac_i"; then + ac_exists=true + break + fi + done + + if test x"$ac_exists" = xtrue; then + : +else + if test "$ac_compiler_gnu" = yes; then + for ac_link_opt in $ac_arg; do + ac_cv_fc_libs="$ac_cv_fc_libs -Xlinker $ac_link_opt" + done +else + ac_cv_fc_libs="$ac_cv_fc_libs $ac_arg" +fi +fi + + ;; + # Ignore these flags. + -lang* | -lcrt[01].o | -lcrtbegin.o | -lc | -lgcc | -libmil | -LANG:=*) + ;; + -lkernel32) + test x"$CYGWIN" != xyes && ac_cv_fc_libs="$ac_cv_fc_libs $ac_arg" + ;; + -[LRuY]) + # These flags, when seen by themselves, take an argument. + # We remove the space between option and argument and re-iterate + # unless we find an empty arg or a new option (starting with -) + case $2 in + "" | -*);; + *) + ac_arg="$ac_arg$2" + shift; shift + set X $ac_arg "$@" + ;; + esac + ;; + -YP,*) + for ac_j in `echo $ac_arg | sed -e 's/-YP,/-L/;s/:/ -L/g'`; do + ac_exists=false + for ac_i in $ac_cv_fc_libs; do + if test x"$ac_j" = x"$ac_i"; then + ac_exists=true + break + fi + done + + if test x"$ac_exists" = xtrue; then + : +else + ac_arg="$ac_arg $ac_j" + ac_cv_fc_libs="$ac_cv_fc_libs $ac_j" +fi + + done + ;; + -[lLR]*) + ac_exists=false + for ac_i in $ac_cv_fc_libs; do + if test x"$ac_arg" = x"$ac_i"; then + ac_exists=true + break + fi + done + + if test x"$ac_exists" = xtrue; then + : +else + ac_cv_fc_libs="$ac_cv_fc_libs $ac_arg" +fi + + ;; + # Ignore everything else. + esac +done +# restore positional arguments +set X $ac_save_positional; shift + +# We only consider "LD_RUN_PATH" on Solaris systems. If this is seen, +# then we insist that the "run path" must be an absolute path (i.e. it +# must begin with a "/"). +case `(uname -sr) 2>/dev/null` in + "SunOS 5"*) + ac_ld_run_path=`echo $ac_fc_v_output | + sed -n 's,^.*LD_RUN_PATH *= *\(/[^ ]*\).*$,-R\1,p'` + test "x$ac_ld_run_path" != x && + if test "$ac_compiler_gnu" = yes; then + for ac_link_opt in $ac_ld_run_path; do + ac_cv_fc_libs="$ac_cv_fc_libs -Xlinker $ac_link_opt" + done +else + ac_cv_fc_libs="$ac_cv_fc_libs $ac_ld_run_path" +fi + ;; +esac +fi # test "x$[]_AC_LANG_PREFIX[]LIBS" = "x" + +fi +echo "$as_me:$LINENO: result: $ac_cv_fc_libs" >&5 +echo "${ECHO_T}$ac_cv_fc_libs" >&6 +FCLIBS="$ac_cv_fc_libs" + + +ac_ext=cc +ac_cpp='$CXXCPP $CPPFLAGS' +ac_compile='$CXX -c $CXXFLAGS $CPPFLAGS conftest.$ac_ext >&5' +ac_link='$CXX -o conftest$ac_exeext $CXXFLAGS $CPPFLAGS $LDFLAGS conftest.$ac_ext $LIBS >&5' +ac_compiler_gnu=$ac_cv_cxx_compiler_gnu + + +echo 'FLIBS = ' $FLIBS + override_f77_libs=0; #case $ac_sys_system in @@ -10704,7 +11781,6 @@ F90BUILDFLAGS=${f90buildopts}' '${savef90flags} - # filename extensions for Fortran 77 if test -z "$F77_EXT"; then F77_EXT=f; fi @@ -11577,6 +12653,13 @@ s,@ECHO_C@,$ECHO_C,;t t s,@ECHO_N@,$ECHO_N,;t t s,@ECHO_T@,$ECHO_T,;t t s,@LIBS@,$LIBS,;t t +s,@F77@,$F77,;t t +s,@FFLAGS@,$FFLAGS,;t t +s,@LDFLAGS@,$LDFLAGS,;t t +s,@ac_ct_F77@,$ac_ct_F77,;t t +s,@EXEEXT@,$EXEEXT,;t t +s,@OBJEXT@,$OBJEXT,;t t +s,@FLIBS@,$FLIBS,;t t s,@BITCOMPILE@,$BITCOMPILE,;t t s,@BITHARDWARE@,$BITHARDWARE,;t t s,@BITCHANGE@,$BITCHANGE,;t t @@ -11627,11 +12710,8 @@ s,@USERDIR@,$USERDIR,;t t s,@INCL_USER_CODE@,$INCL_USER_CODE,;t t s,@CXX@,$CXX,;t t s,@CXXFLAGS@,$CXXFLAGS,;t t -s,@LDFLAGS@,$LDFLAGS,;t t s,@CPPFLAGS@,$CPPFLAGS,;t t s,@ac_ct_CXX@,$ac_ct_CXX,;t t -s,@EXEEXT@,$EXEEXT,;t t -s,@OBJEXT@,$OBJEXT,;t t s,@use_sundials@,$use_sundials,;t t s,@CVODE_LIBS@,$CVODE_LIBS,;t t s,@IDA_LIBS@,$IDA_LIBS,;t t @@ -11704,10 +12784,10 @@ s,@CXX_INCLUDES@,$CXX_INCLUDES,;t t s,@LCXX_FLAGS@,$LCXX_FLAGS,;t t s,@LCXX_END_LIBS@,$LCXX_END_LIBS,;t t s,@HAVE_STRIPSYMBOLS@,$HAVE_STRIPSYMBOLS,;t t -s,@F77@,$F77,;t t -s,@FFLAGS@,$FFLAGS,;t t -s,@ac_ct_F77@,$ac_ct_F77,;t t -s,@FLIBS@,$FLIBS,;t t +s,@FC@,$FC,;t t +s,@FCFLAGS@,$FCFLAGS,;t t +s,@ac_ct_FC@,$ac_ct_FC,;t t +s,@FCLIBS@,$FCLIBS,;t t s,@F90@,$F90,;t t s,@BUILD_F90@,$BUILD_F90,;t t s,@F90FLAGS@,$F90FLAGS,;t t diff --git a/docs/Cantera.cfg.in b/docs/Cantera.cfg.in index 789deaf44..f3f3a4dab 100755 --- a/docs/Cantera.cfg.in +++ b/docs/Cantera.cfg.in @@ -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. diff --git a/docs/install_examples/linux.64_sierra_gcc444_python264_numpy b/docs/install_examples/linux.64_sierra_gcc444_python264_numpy index 32760d503..7ab813602 100755 --- a/docs/install_examples/linux.64_sierra_gcc444_python264_numpy +++ b/docs/install_examples/linux.64_sierra_gcc444_python264_numpy @@ -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 diff --git a/ext/f2c_lapack/Makefile.in b/ext/f2c_lapack/Makefile.in index 16c2f1a14..4b0e77110 100755 --- a/ext/f2c_lapack/Makefile.in +++ b/ext/f2c_lapack/Makefile.in @@ -111,6 +111,9 @@ dtrtrs.o \ dgecon.o \ dgeequ.o \ dgerfs.o \ +dgbcon.o \ +dgbequ.o \ +dlatbs.o \ ieeeck.o \ ilaenv.o diff --git a/ext/f2c_lapack/dgbcon.c b/ext/f2c_lapack/dgbcon.c new file mode 100644 index 000000000..72083f759 --- /dev/null +++ b/ext/f2c_lapack/dgbcon.c @@ -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_ */ + diff --git a/ext/f2c_lapack/dgbequ.c b/ext/f2c_lapack/dgbequ.c new file mode 100644 index 000000000..1c26aadfd --- /dev/null +++ b/ext/f2c_lapack/dgbequ.c @@ -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_ */ + diff --git a/ext/f2c_lapack/dlatbs.c b/ext/f2c_lapack/dlatbs.c new file mode 100644 index 000000000..eaea649f4 --- /dev/null +++ b/ext/f2c_lapack/dlatbs.c @@ -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_ */ + diff --git a/ext/lapack/Makefile.in b/ext/lapack/Makefile.in index e0f242113..b71f30b4f 100755 --- a/ext/lapack/Makefile.in +++ b/ext/lapack/Makefile.in @@ -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 \ diff --git a/ext/lapack/dgbcon.f b/ext/lapack/dgbcon.f new file mode 100644 index 000000000..ba613b735 --- /dev/null +++ b/ext/lapack/dgbcon.f @@ -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 diff --git a/ext/lapack/dgbequ.f b/ext/lapack/dgbequ.f new file mode 100644 index 000000000..309798c98 --- /dev/null +++ b/ext/lapack/dgbequ.f @@ -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 diff --git a/ext/lapack/dlatbs.f b/ext/lapack/dlatbs.f new file mode 100644 index 000000000..fe7dd11cb --- /dev/null +++ b/ext/lapack/dlatbs.f @@ -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 diff --git a/test_problems/cathermo/HMW_test_1/output_blessed.txt b/test_problems/cathermo/HMW_test_1/output_blessed.txt index 3d63f73b0..296945da9 100644 --- a/test_problems/cathermo/HMW_test_1/output_blessed.txt +++ b/test_problems/cathermo/HMW_test_1/output_blessed.txt @@ -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 diff --git a/test_problems/cathermo/HMW_test_1/output_noD_blessed.txt b/test_problems/cathermo/HMW_test_1/output_noD_blessed.txt index 966d76de0..296945da9 100644 --- a/test_problems/cathermo/HMW_test_1/output_noD_blessed.txt +++ b/test_problems/cathermo/HMW_test_1/output_noD_blessed.txt @@ -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