diff --git a/Cantera/src/numerics/BandMatrix.cpp b/Cantera/src/numerics/BandMatrix.cpp index b643fe1b2..0571df81c 100644 --- a/Cantera/src/numerics/BandMatrix.cpp +++ b/Cantera/src/numerics/BandMatrix.cpp @@ -22,160 +22,243 @@ using namespace std; namespace Cantera { + //==================================================================================================================== + BandMatrix::BandMatrix() : + m_factored(false), + m_n(0), + m_kl(0), + m_ku(0), + m_zero(0.0) + { + data.clear(); + ludata.clear(); + } + //==================================================================================================================== + BandMatrix::BandMatrix(int n, int kl, int ku, doublereal v) : + m_factored(false), + m_n(n), + m_kl(kl), + m_ku(ku), + m_zero(0.0) + { + data.resize(n*(2*kl + ku + 1)); + ludata.resize(n*(2*kl + ku + 1)); + fill(data.begin(), data.end(), v); + fill(ludata.begin(), ludata.end(), 0.0); + m_ipiv.resize(m_n); + } + //==================================================================================================================== + BandMatrix::BandMatrix(const BandMatrix& y) : + m_factored(false), + m_n(0), + m_kl(0), + m_ku(0), + m_zero(0.0) + { + m_n = y.m_n; + m_kl = y.m_kl; + m_ku = y.m_ku; + data = y.data; + ludata = y.ludata; + m_factored = y.m_factored; + m_ipiv = y.m_ipiv; + } + //==================================================================================================================== + BandMatrix::~BandMatrix() { - /// Default constructor. - BandMatrix::BandMatrix() : m_factored(false), m_n(0), - m_kl(0), m_ku(0), m_zero(0.0) { - data.clear(); ludata.clear(); + } + //==================================================================================================================== + BandMatrix& BandMatrix::operator=(const BandMatrix & y) { + if (&y == this) return *this; + m_n = y.m_n; + m_kl = y.m_kl; + m_ku = y.m_ku; + m_ipiv = y.m_ipiv; + data = y.data; + ludata = y.ludata; + m_factored = y.m_factored; + return *this; + } + //==================================================================================================================== + void BandMatrix::resize(int n, int kl, int ku, doublereal v) { + m_n = n; + m_kl = kl; + m_ku = ku; + data.resize(n*(2*kl + ku + 1)); + ludata.resize(n*(2*kl + ku + 1)); + m_ipiv.resize(m_n); + fill(data.begin(), data.end(), v); + m_factored = false; + } + //==================================================================================================================== + void BandMatrix::bfill(doublereal v) { + std::fill(data.begin(), data.end(), v); + m_factored = false; + } + //==================================================================================================================== + doublereal& BandMatrix::operator()(int i, int j) { + return value(i,j); + } + //==================================================================================================================== + doublereal BandMatrix::operator()(int i, int j) const { + return value(i,j); + } + //==================================================================================================================== + doublereal& BandMatrix::value(int i, int j) { + m_factored = false; + if (i < j - m_ku || i > j + m_kl) { + return m_zero; } - - - /** - * Constructor. Create an n by n banded matrix. - * @param n number of rows and columns - * @param kl number of subdiagonals - * @param ku number of superdiagonals - * @param v initial value (default = 0.0) - */ - BandMatrix::BandMatrix(int n, int kl, int ku, doublereal v) - : m_factored(false), m_n(n), m_kl(kl), m_ku(ku) { - data.resize(n*(2*kl + ku + 1)); - ludata.resize(n*(2*kl + ku + 1)); - fill(data.begin(), data.end(), v); - fill(ludata.begin(), ludata.end(), 0.0); - m_ipiv.resize(m_n); + return data[index(i,j)]; + } + //==================================================================================================================== + doublereal BandMatrix::value( int i, int j) const { + if (i < j - m_ku || i > j + m_kl) return 0.0; + return data[index(i,j)]; + } + //==================================================================================================================== + int BandMatrix::index(int i, int j) const { + int rw = m_kl + m_ku + i - j; + return (2*m_kl + m_ku + 1)*j + rw; + } + //==================================================================================================================== + doublereal BandMatrix::_value(int i, int j) const { + return data[index(i,j)]; + } + //==================================================================================================================== + // Number of rows + int BandMatrix::nRows() const { + return m_n; + } + //==================================================================================================================== + // Number of columns + int BandMatrix::nColumns() const { + return m_n; + } + //==================================================================================================================== + // Number of subdiagonals + int BandMatrix::nSubDiagonals() const { + return m_kl; + } + //==================================================================================================================== + // Number of superdiagonals + int BandMatrix::nSuperDiagonals() const { + return m_ku; + } + //==================================================================================================================== + int BandMatrix::ldim() const { + return 2*m_kl + m_ku + 1; + } + //==================================================================================================================== + vector_int & BandMatrix::ipiv() { + return m_ipiv; + } + //==================================================================================================================== + /* + * Multiply A*b and write result to \c prod. + */ + void BandMatrix::mult(const doublereal * const b, doublereal * const prod) const { + int nr = nRows(); + doublereal sum = 0.0; + for (int m = 0; m < nr; m++) { + sum = 0.0; + for (int j = m - m_kl; j <= m + m_ku; j++) { + if (j >= 0 && j < m_n) + sum += _value(m,j) * b[j]; + } + prod[m] = sum; } - - /// copy constructor - BandMatrix::BandMatrix(const BandMatrix& y) { - m_n = y.m_n; - m_kl = y.m_kl; - m_ku = y.m_ku; - data = y.data; - ludata = y.ludata; - m_factored = y.m_factored; - m_ipiv = y.m_ipiv; + } + //==================================================================================================================== + /* + * Multiply b*A and write result to \c prod. + */ + void BandMatrix::leftMult(const doublereal * const b, doublereal * const prod) const { + int nc = nColumns(); + doublereal sum = 0.0; + for (int n = 0; n < nc; n++) { + sum = 0.0; + for (int i = n - m_ku; i <= n + m_kl; i++) { + if (i >= 0 && i < m_n) + sum += _value(i,n) * b[i]; + } + prod[n] = sum; } + } + //==================================================================================================================== + /* + * Perform an LU decomposition. LAPACK routine DGBTRF is used. + * The factorization is saved in ludata. + */ + int BandMatrix::factor() { + int info=0; + copy(data.begin(), data.end(), ludata.begin()); + ct_dgbtrf(nRows(), nColumns(), nSubDiagonals(), nSuperDiagonals(), + DATA_PTR(ludata), ldim(), DATA_PTR(ipiv()), info); - BandMatrix& BandMatrix::operator=(const BandMatrix& y) { - if (&y == this) return *this; - m_n = y.m_n; - m_kl = y.m_kl; - m_ku = y.m_ku; - m_ipiv = y.m_ipiv; - data = y.data; - ludata = y.ludata; - m_factored = y.m_factored; - return *this; + // if info = 0, LU decomp succeeded. + if (info == 0) { + m_factored = true; + } else { + m_factored = false; + ofstream fout("bandmatrix.csv"); + fout << *this << endl; + fout.close(); } + 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(int n, doublereal* b) { + int info = 0; + if (!m_factored) info = factor(); + if (info == 0) + ct_dgbtrs(ctlapack::NoTranspose, nColumns(), nSubDiagonals(), + nSuperDiagonals(), 1, DATA_PTR(ludata), ldim(), + DATA_PTR(ipiv()), b, nColumns(), info); - void BandMatrix::resize(int n, int kl, int ku, doublereal v) { - m_n = n; - m_kl = kl; - m_ku = ku; - data.resize(n*(2*kl + ku + 1)); - ludata.resize(n*(2*kl + ku + 1)); - m_ipiv.resize(m_n); - fill(data.begin(), data.end(), v); - fill(data.begin(), data.end(), 0.0); - m_factored = false; + // error handling + if (info != 0) { + ofstream fout("bandmatrix.csv"); + fout << *this << endl; + fout.close(); } - - - /** - * Multiply A*b and write result to \c prod. - */ - void BandMatrix::mult(const double* b, double* prod) const { - int nr = rows(); - int m, j; - double sum = 0.0; - for (m = 0; m < nr; m++) { - sum = 0.0; - for (j = m - m_kl; j <= m + m_ku; j++) { - if (j >= 0 && j < m_n) - sum += _value(m,j)*b[j]; - } - prod[m] = sum; - } - } - - - /** - * Multiply b*A and write result to \c prod. - */ - void BandMatrix::leftMult(const double* b, double* prod) const { - int nc = columns(); - int n, i; - double sum = 0.0; - for (n = 0; n < nc; n++) { - sum = 0.0; - for (i = n - m_ku; i <= n + m_kl; i++) { - if (i >= 0 && i < m_n) - sum += _value(i,n)*b[i]; - } - prod[n] = sum; - } - } - - - /** - * Perform an LU decomposition. LAPACK routine DGBTRF is used. - * The factorization is saved in ludata. - */ - int BandMatrix::factor() { - int info=0; - copy(data.begin(), data.end(), ludata.begin()); - ct_dgbtrf(rows(), columns(), nSubDiagonals(), nSuperDiagonals(), - DATA_PTR(ludata), ldim(), DATA_PTR(ipiv()), info); - - // if info = 0, LU decomp succeeded. - if (info == 0) { - m_factored = true; - } - else { - m_factored = false; - ofstream fout("bandmatrix.csv"); - fout << *this << endl; - fout.close(); - } - return info; - } - - - - int BandMatrix::solve(int n, const doublereal* b, doublereal* x) { - copy(b, b+n, x); - return solve(n, x); - } - - int BandMatrix::solve(int n, doublereal* b) { - int info = 0; - if (!m_factored) info = factor(); - if (info == 0) - ct_dgbtrs(ctlapack::NoTranspose, columns(), nSubDiagonals(), - nSuperDiagonals(), 1, DATA_PTR(ludata), ldim(), - DATA_PTR(ipiv()), b, columns(), info); - - // error handling - if (info != 0) { - ofstream fout("bandmatrix.csv"); - fout << *this << endl; - fout.close(); - } - return info; - } - - ostream& operator<<(ostream& s, const BandMatrix& m) { - int nr = m.rows(); - int nc = m.columns(); - int i,j; - for (i = 0; i < nr; i++) { - for (j = 0; j < nc; j++) { - s << m(i,j) << ", "; - } - s << endl; - } - return s; + return info; + } + //==================================================================================================================== + vector_fp::iterator BandMatrix::begin() { + m_factored = false; + return data.begin(); + } + //==================================================================================================================== + vector_fp::iterator BandMatrix::end() { + m_factored = false; + return data.end(); + } + //==================================================================================================================== + vector_fp::const_iterator BandMatrix::begin() const { + return data.begin(); + } + //==================================================================================================================== + vector_fp::const_iterator BandMatrix::end() const { + return data.end(); + } + //==================================================================================================================== + ostream& operator<<(ostream& s, const BandMatrix& m) { + int nr = m.nRows(); + int nc = m.nColumns(); + for (int i = 0; i < nr; i++) { + for (int j = 0; j < nc; j++) { + s << m(i,j) << ", "; + } + s << endl; } + return s; + } + //==================================================================================================================== } diff --git a/Cantera/src/numerics/BandMatrix.h b/Cantera/src/numerics/BandMatrix.h index d3f7fdd3e..148029c95 100644 --- a/Cantera/src/numerics/BandMatrix.h +++ b/Cantera/src/numerics/BandMatrix.h @@ -23,130 +23,258 @@ namespace Cantera { - /** - * A class for banded matrices. + /** + * A class for banded matrices. This class has matrix inversion processes. + * The class is based upon the LAPACK banded storage matrix format. + */ + class BandMatrix { + + public: + + //! Base Constructor + /*! + * * Create an \c 0 by \c 0 matrix, and initialize all elements to \c 0. */ - class BandMatrix { + BandMatrix(); - public: + //! Creates a banded matrix and sets all elements to zero + /*! + * Create an \c n by \c n banded matrix, and initialize all elements to \c v. + * + * @param n size of the square matrix + * @param kl band size on the lower portion of the matrix + * @param ku band size on the upper portion of the matrix + * @param v intial value of all matrix components. + */ + BandMatrix(int n, int kl, int ku, doublereal v = 0.0); - BandMatrix(); - BandMatrix(int n, int kl, int ku, doublereal v = 0.0); + //! Copy constructor + /*! + * @param y Matrix to be copied + */ + BandMatrix(const BandMatrix& y); - /// copy constructor - BandMatrix(const BandMatrix& y); + //! Destructor. Does nothing. + virtual ~BandMatrix(); - /// Destructor. Does nothing. - virtual ~BandMatrix(){} + //! assignment operator + /*! + * @param y reference to the matrix to be copied + */ + BandMatrix& operator=(const BandMatrix& y); - /// assignment. - BandMatrix& operator=(const BandMatrix& y); + //! Resize the matrix problem + /*! + * All data is lost + * + * @param n size of the square matrix + * @param kl band size on the lower portion of the matrix + * @param ku band size on the upper portion of the matrix + * @param v intial value of all matrix components. + */ + void resize(int n, int kl, int ku, doublereal v = 0.0); - void resize(int n, int kl, int ku, doublereal v = 0.0); + //! Fill or zero the matrix + /*! + * @param v Fill value, defaults to zero. + */ + void bfill(doublereal v = 0.0); - void bfill(doublereal v) { - std::fill(data.begin(), data.end(), v); - m_factored = false; - } - - doublereal& operator()( int i, int j) { - return value(i,j); - } - - doublereal operator() ( int i, int j) const { - return value(i,j); - } - - /// Return a reference to element (i,j). Since this method may - /// alter the element value, it may need to be refactored, so - /// the flag m_factored is set to false. - doublereal& value( int i, int j) { - m_factored = false; - if (i < j - m_ku || i > j + m_kl) { - m_zero = 0.0; - return m_zero; - } - return data[index(i,j)]; - } - - /// Return the value of element (i,j). This method does not - /// alter the array. - doublereal value( int i, int j) const { - if (i < j - m_ku || i > j + m_kl) return 0.0; - return data[index(i,j)]; - } - - /// Return the location in the internal 1D array corresponding to - /// the (i,j) element in the banded array. - int index(int i, int j) const { - int rw = m_kl + m_ku + i - j; - return (2*m_kl + m_ku + 1)*j + rw; - } - - /// Return the value of the (i,j) element for (i,j) within the - /// bandwidth. For efficiency, this method does not check that - /// (i,j) are within the bandwidth; it is up to the calling - /// program to insure that this is true. - doublereal _value(int i, int j) const { - return data[index(i,j)]; - } - - /// Number of rows - int nRows() const { return m_n; } - /// @deprecated Redundant. - int rows() const { return m_n; } - - /// Number of columns - int nColumns() const { return m_n; } - /// @deprecated Redundant. - int columns() const { return m_n; } - - /// Number of subdiagonals - int nSubDiagonals() const { return m_kl; } - - /// Number of superdiagonals - int nSuperDiagonals() const { return m_ku; } - - int ldim() const { return 2*m_kl + m_ku + 1; } - vector_int& ipiv() { return m_ipiv; } - - /// Multiply A*b and write result to prod. - void mult(const double* b, double* prod) const; - - /// Multiply b*A and write result to prod. - void leftMult(const double* b, double* prod) const; - - int factor(); - - //void solve(const vector_fp& b, vector_fp& x); - - int solve(int n, const doublereal* b, doublereal* x); - int solve(int n, doublereal* b); - - vector_fp::iterator begin() { - m_factored = false; - return data.begin(); - } - vector_fp::iterator end() { - m_factored = false; - return data.end(); - } - vector_fp::const_iterator begin() const { return data.begin(); } - vector_fp::const_iterator end() const { return data.end(); } - - protected: - vector_fp data; - vector_fp ludata; - bool m_factored; + //! Index into the (i,j) element + /*! + * @param i row + * @param j column + * + * Returns a changeable reference to the matrix entry + */ + doublereal& operator()(int i, int j); - int m_n, m_kl, m_ku; - doublereal m_zero; - vector_int m_ipiv; + //! Constant Index into the (i,j) element + /*! + * @param i row + * @param j column + * + * Returns an unchangeable reference to the matrix entry + */ + doublereal operator() (int i, int j) const; - }; + //! Return a changeable reference to element (i,j). + /*! + * Since this method may alter the element value, it may need to be refactored, so + * the flag m_factored is set to false. + * + * @param i row + * @param j column + * + * @return Returns a reference to the value of the matrix entry + */ + doublereal& value( int i, int j); - std::ostream& operator<<(std::ostream& s, const BandMatrix& m); + + //! Return the value of element (i,j). + /*! + * This method does not alter the array. + * @param i row + * @param j column + * + * @return Returns the value of the matrix entry + */ + doublereal value( int i, int j) const; + + //! Returns the location in the internal 1D array corresponding to the (i,j) element in the banded array + /*! + * @param i row + * @param j column + * + * @return Returns the index of the matrix entry + */ + int index(int i, int j) const; + + //! Return the value of the (i,j) element for (i,j) within the bandwidth. + /*! + * For efficiency, this method does not check that (i,j) are within the bandwidth; it is up to the calling + * program to insure that this is true. + * + * @param i row + * @param j column + * + * @return Returns the value of the matrix entry + */ + doublereal _value(int i, int j) const; + + //! Returns the number of rows + int nRows() const; + + //! Number of columns + int nColumns() const; + + //! Number of subdiagonals + int nSubDiagonals() const; + + //! Number of superdiagonals + int nSuperDiagonals() const; + + //! Return the number of rows of storage needed for the band storage + int ldim() const; + + //! Return a reference to the pivot vector + /*! + * @return return a reference to the pivot vector + */ + vector_int& ipiv(); + + //! 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 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; + + //! Perform an LU decomposition, the LAPACK routine DGBTRF is used. + /*! + * + * The factorization is saved in ludata. + * + * @return Return a success flag. + * 0 indicates a success + * ~0 Some error occurred, see the LAPACK documentation + */ + int factor(); + + + //! 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 + * + * @return Return a success flag + * 0 indicates a success + * ~0 Some error occurred, see the LAPACK documentation + */ + int solve(int n, 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 + * + * @return Return a success flag + * 0 indicates a success + * ~0 Some error occurred, see the LAPACK documentation + */ + int solve(int n, 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(); + + //! Returns an iterator for the end of the band storage data + /*! + * Iterator points to the end of the data, and it is changeable. + */ + vector_fp::iterator end(); + + //! Returns a const iterator for the start of the band storage data + /*! + * Iterator points to the beginning of the data, and it is not changeable. + */ + vector_fp::const_iterator begin() const; + + //! Returns a const iterator for the end of the band storage data + /*! + * Iterator points to the end of the data, and it is not changeable. + */ + vector_fp::const_iterator end() const; + + protected: + + //! Matrix data + vector_fp data; + + //! Factorized data + vector_fp ludata; + + //! Boolean indicating whether the matrix is factored + bool m_factored; + + //! Number of rows and columns of the matrix + int m_n; + + //! Number of subdiagonals of the matrix + int m_kl; + + //! Number of super diagonals of the matrix + int m_ku; + + //! value of zero + doublereal m_zero; + + //! Pivot vector + vector_int m_ipiv; + + }; + + //! Utility routine to print out the matrix + /*! + * @param s ostream to print the matrix out to + * @param m Matrix to be printed + * + * @return Returns a reference to the ostream + */ + std::ostream& operator<<(std::ostream& s, const BandMatrix& m); } diff --git a/Cantera/src/numerics/NonlinearSolver.cpp b/Cantera/src/numerics/NonlinearSolver.cpp index 91fa49d4b..962389c64 100644 --- a/Cantera/src/numerics/NonlinearSolver.cpp +++ b/Cantera/src/numerics/NonlinearSolver.cpp @@ -716,8 +716,8 @@ 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, double* y_comm, double* ydot_comm, doublereal time_curr, - int num_newt_its) + void NonlinearSolver::scaleMatrix(SquareMatrix& jac, doublereal * const y_comm, doublereal * const ydot_comm, + doublereal time_curr, int num_newt_its) { int irow, jcol; /* @@ -2319,10 +2319,10 @@ namespace Cantera { * NSOLN_RETN_FAIL_DAMPSTEP * Unsuccessful step. We can not find a damping factor that is suitable. */ - int NonlinearSolver::dampStep(const doublereal time_curr, const doublereal * y_n_curr, - const doublereal * ydot_n_curr, doublereal * const step_1, - doublereal * const y_n_1, double* const ydot_n_1, doublereal * const step_2, - double& stepNorm_2, SquareMatrix& jac, bool writetitle, int& num_backtracks) + 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) { int j, m; int info = 0; @@ -2729,7 +2729,7 @@ namespace Cantera { */ int NonlinearSolver::decideStep(const doublereal time_curr, int leg, doublereal alpha, const doublereal * const y_n_curr, const doublereal * const ydot_n_curr, const std::vector & step_1, - const doublereal * const y_n_1, const double* const ydot_n_1, doublereal trustDeltaOld) + const doublereal * const y_n_1, const doublereal * const ydot_n_1, doublereal trustDeltaOld) { int retn = 2; bool goodStep = false; @@ -2906,8 +2906,8 @@ namespace Cantera { * @return A positive value indicates a successful convergence * -1 Failed convergence */ - int NonlinearSolver::solve_nonlinear_problem(int SolnType, double* y_comm, double* ydot_comm, doublereal CJ, - doublereal time_curr, SquareMatrix& jac, + int NonlinearSolver::solve_nonlinear_problem(int SolnType, doublereal * const y_comm, doublereal * const ydot_comm, + doublereal CJ, doublereal time_curr, SquareMatrix& jac, int &num_newt_its, int &num_linear_solves, int &num_backtracks, int loglevelInput) { diff --git a/Cantera/src/numerics/NonlinearSolver.h b/Cantera/src/numerics/NonlinearSolver.h index a2861dddc..582c3fa6c 100644 --- a/Cantera/src/numerics/NonlinearSolver.h +++ b/Cantera/src/numerics/NonlinearSolver.h @@ -311,6 +311,7 @@ namespace Cantera { */ void setDefaultDeltaBoundsMagnitudes(); + //! Adjust the step minimums void adjustUpStepMinimums(); //! Set the delta Bounds magnitudes by hand @@ -506,10 +507,10 @@ namespace Cantera { * * @return returns an integer indicating what happened. */ - int dampStep(const doublereal time_curr, const double* y_n_curr, - const doublereal *ydot_n_curr, double * const step_1, - double* const y_n_1, double* const ydot_n_1, double* step_2, - double& stepNorm_2, SquareMatrix& jac, bool writetitle, + 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, int& num_backtracks); //! Find the solution to F(X) = 0 by damped Newton iteration. @@ -539,7 +540,7 @@ namespace Cantera { * @return A positive value indicates a successful convergence * -1 Failed convergence */ - int solve_nonlinear_problem(int SolnType, double* y_comm,double* ydot_comm, doublereal CJ, + 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, int &num_linear_solves, int &num_backtracks, int loglevelInput); @@ -586,8 +587,10 @@ namespace Cantera { * @param y_comm Current value of the solution vector * @param ydot_comm Current value of the time derivative of the solution vector * @param time_curr current value of the time + * @param num_newt_its Current value of the number of newt its */ - void scaleMatrix(SquareMatrix& jac, double* y_comm, double* ydot_comm, doublereal time_curr, int num_newt_its); + void scaleMatrix(SquareMatrix& jac, doublereal * const y_comm, doublereal * const ydot_comm, + doublereal time_curr, int num_newt_its); //! Print solution norm contribution /*! @@ -802,11 +805,11 @@ namespace Cantera { * @param time_curr INPUT Current value of the time * @param leg INPUT Leg of the dogleg that we are on * @param alpha INPUT Distance down that leg that we are on - * @param y0 INPUT Current value of the solution vector - * @param ydot0 INPUT Current value of the derivative of the solution vector - * @param step0 INPUT Trial step - * @param y1 OUTPUT Solution values at the conditions which are evalulated for success - * @param ydot1 OUTPUT Time derivates of solution at the conditions which are evalulated for success + * @param y_n_curr INPUT Current value of the solution vector + * @param ydot_n_curr INPUT Current value of the derivative of the solution vector + * @param step_1 INPUT Trial step + * @param y_n_1 OUTPUT Solution values at the conditions which are evalulated for success + * @param ydot_n_1 OUTPUT Time derivates of solution at the conditions which are evalulated for success * @param trustDeltaOld INPUT Value of the trust length at the old conditions * * @@ -819,10 +822,10 @@ namespace Cantera { * -2 Current value of the solution vector caused a residual error in its evaluation. * Step is a failure, and the step size must be reduced in order to proceed further. */ - int decideStep(const doublereal time_curr, int leg, doublereal alpha, const doublereal * const y0, - const doublereal * const ydot0, - const std::vector & step0, - const doublereal * const y1, const doublereal * const ydot1, doublereal trustDeltaOld); + int decideStep(const doublereal time_curr, int leg, doublereal alpha, const doublereal * const y_n_curr, + const doublereal * const ydot_n_curr, + const std::vector & step_1, + const doublereal * const y_n_1, const doublereal * const ydot_n_1, doublereal trustDeltaOld); //! Calculated the expected residual along the double dogleg curve. /*! @@ -923,6 +926,7 @@ namespace Cantera { //! Value of the solution time derivative at the new point that is to be considered std::vector m_ydot_n_1; + //! Value of the step to be taken in the solution std::vector m_step_1; //! Vector of column scaling factors