From 75d9ef93c5a751ac104f1c17b9b4f36a749c363e Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 11 Nov 2015 21:04:12 -0500 Subject: [PATCH] Clean up Doxygen docs and comments for numerics classes --- include/cantera/numerics/BandMatrix.h | 165 ++++++------ include/cantera/numerics/DAE_Solver.h | 29 +-- include/cantera/numerics/DenseMatrix.h | 172 ++++++------ include/cantera/numerics/Func1.h | 18 +- include/cantera/numerics/GeneralMatrix.h | 88 +++---- include/cantera/numerics/IDA_Solver.h | 39 +-- include/cantera/numerics/Integrator.h | 13 +- include/cantera/numerics/ResidEval.h | 33 ++- include/cantera/numerics/ResidJacEval.h | 101 ++++---- include/cantera/numerics/RootFind.h | 316 +++++++++++------------ include/cantera/numerics/SquareMatrix.h | 26 +- include/cantera/numerics/funcs.h | 12 +- include/cantera/numerics/polyfit.h | 96 ++++--- src/numerics/BandMatrix.cpp | 6 +- src/numerics/CVodesIntegrator.cpp | 26 +- src/numerics/DAE_solvers.cpp | 4 +- src/numerics/DenseMatrix.cpp | 3 +- src/numerics/GeneralMatrix.cpp | 4 +- src/numerics/IDA_Solver.cpp | 51 ++-- src/numerics/ResidJacEval.cpp | 4 +- src/numerics/RootFind.cpp | 155 ++++------- src/numerics/SquareMatrix.cpp | 23 +- src/numerics/funcs.cpp | 6 +- 23 files changed, 638 insertions(+), 752 deletions(-) diff --git a/include/cantera/numerics/BandMatrix.h b/include/cantera/numerics/BandMatrix.h index 0e827628e..9cf72c2f0 100644 --- a/include/cantera/numerics/BandMatrix.h +++ b/include/cantera/numerics/BandMatrix.h @@ -15,33 +15,33 @@ namespace Cantera { -//! A class for banded matrices, involving 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. + * 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. + * 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, within - * the oneD directory, the class is used as is, without reference to the GeneralMatrix - * base type. + * This class is a derived class of the base class GeneralMatrix. However, + * within the oneD directory, the class is used as is, without reference to the + * GeneralMatrix base type. */ class BandMatrix : public GeneralMatrix { public: //! Base Constructor /*! - * * Create an \c 0 by \c 0 matrix, and initialize all elements to \c 0. + * Create an \c 0 by \c 0 matrix, and initialize all elements to \c 0. */ BandMatrix(); //! 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. + * 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 @@ -50,23 +50,14 @@ public: */ BandMatrix(size_t n, size_t kl, size_t ku, doublereal v = 0.0); - //! Copy constructor - /*! - * @param y Matrix to be copied - */ BandMatrix(const BandMatrix& y); - - //! assignment operator - /*! - * @param y reference to the matrix to be copied - */ BandMatrix& operator=(const BandMatrix& y); //! Resize the matrix problem /*! - * All data is lost + * All data is lost * - * @param n size of the square matrix + * @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 initial value of all matrix components. @@ -84,40 +75,41 @@ public: //! 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. + * 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 + * @param i row + * @param j column + * @returns a reference to the value of the matrix entry */ doublereal& value(size_t i, size_t j); //! 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 + * This method does not alter the array. + * @param i row + * @param j column + * @returns the value of the matrix entry */ doublereal value(size_t i, size_t j) const; - //! Returns the location in the internal 1D array corresponding to the (i,j) element in the banded array + //! 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 + * @param i row + * @param j column + * @returns the index of the matrix entry */ size_t index(size_t i, size_t 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. + * 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 + * @param i row + * @param j column + * @returns the value of the matrix entry */ doublereal _value(size_t i, size_t j) const; @@ -125,10 +117,12 @@ public: //! Return the size and structure of the matrix /*! - * @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. + * @param iStruct OUTPUT Pointer to a vector of ints that describe the + * structure of the matrix. + * + * istruct[0] = kl + * istruct[1] = ku + * @returns the number of rows and columns in the matrix. */ virtual size_t nRowsAndStruct(size_t* const iStruct = 0) const; @@ -155,55 +149,52 @@ public: /*! * The factorization is saved in ludata. * - * @return Return a success flag. - * 0 indicates a success - * ~0 Some error occurred, see the LAPACK documentation + * @return Return a success flag. 0 indicates a success; ~0 indicates some + * error occurred, see the LAPACK documentation */ int factor(); //! Solve the matrix problem Ax = b /*! - * @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 + * @param b INPUT RHS of the problem + * @param x OUTPUT solution to the problem + * @return a success flag. 0 indicates a success; ~0 indicates some error + * occurred, see the LAPACK documentation */ int solve(const doublereal* const b, doublereal* const x); //! Solve the matrix problem Ax = b /*! - * @param b INPUT RHS of the problem - * OUTPUT solution to the problem - * @param nrhs Number of right hand sides to solve - * @param ldb Leading dimension of `b`. Default is nColumns() - * @return Return a success flag - * 0 indicates a success - * ~0 Some error occurred, see the LAPACK documentation + * @param b INPUT RHS of the problem + * OUTPUT solution to the problem + * @param nrhs Number of right hand sides to solve + * @param ldb Leading dimension of `b`. Default is nColumns() + * @returns a success flag. 0 indicates a success; ~0 indicates some error + * occurred, see the LAPACK documentation */ int solve(doublereal* b, size_t nrhs=1, size_t ldb=0); //! Returns an iterator for the start of the band storage data /*! - * Iterator points to the beginning of the data, and it is changeable. + * Iterator points to the beginning of the data, and it is changeable. */ virtual 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. + * 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. + * 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. + * Iterator points to the end of the data, and it is not changeable. */ vector_fp::const_iterator end() const; @@ -211,10 +202,10 @@ public: //! Returns an estimate of the inverse of the condition number for the matrix /*! - * The matrix must have been previously factored using the LU algorithm + * 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 + * @returns the inverse of the condition number */ virtual doublereal rcond(doublereal a1norm); @@ -227,42 +218,43 @@ public: virtual GeneralMatrix* duplMyselfAsGeneralMatrix() const; - //! Return a pointer to the top of column j, column values are assumed to be contiguous in memory + //! Return a pointer to the top of column j /*! - * The LAPACK bandstructure has column values which are contiguous in memory: + * Column values are assumed to be contiguous in memory (LAPACK band matrix + * structure) * - * 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) + * 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: * - * This routine returns the position of AB(1,j) (fortran-1 indexing) in the above format + * AB(KL + KU + 1 + i - j,j) = A(i,j) for max(1, j - KU) <= i <= min(m, j + KL) * - * So to address the (i,j) position, you use the following indexing: + * 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 + * @returns a pointer to the top of the column */ virtual doublereal* ptrColumn(size_t 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. + * 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. + * @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 + * 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 @@ -271,8 +263,8 @@ public: //! 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 + * 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 @@ -315,8 +307,7 @@ protected: /*! * @param s ostream to print the matrix out to * @param m Matrix to be printed - * - * @return Returns a reference to the ostream + * @returns a reference to the ostream */ std::ostream& operator<<(std::ostream& s, const BandMatrix& m); diff --git a/include/cantera/numerics/DAE_Solver.h b/include/cantera/numerics/DAE_Solver.h index adef4fbea..53f647074 100644 --- a/include/cantera/numerics/DAE_Solver.h +++ b/include/cantera/numerics/DAE_Solver.h @@ -141,24 +141,22 @@ public: //! Calculate consistent value of the starting solution given the starting //! solution derivatives /** - * This method may be called if the initial conditions do not - * satisfy the residual equation F = 0. Given the derivatives - * of all variables, this method computes the initial y - * values. + * This method may be called if the initial conditions do not satisfy the + * residual equation F = 0. Given the derivatives of all variables, this + * method computes the initial y values. */ virtual void correctInitial_Y_given_Yp(doublereal* y, doublereal* yp, doublereal tout) { warn("correctInitial_Y_given_Yp"); } - //! Calculate consistent value of the algebraic constraints and - //! derivatives at the start of the problem + //! Calculate consistent value of the algebraic constraints and derivatives + //! at the start of the problem /** - * This method may be called if the initial conditions do not - * satisfy the residual equation F = 0. Given the initial - * values of all differential variables, it computes the - * initial values of all algebraic variables and the initial - * derivatives of all differential variables. + * This method may be called if the initial conditions do not satisfy the + * residual equation F = 0. Given the initial values of all differential + * variables, it computes the initial values of all algebraic variables and + * the initial derivatives of all differential variables. * @param y Calculated value of the solution vector after the procedure ends * @param yp Calculated value of the solution derivative after the procedure * @param tout The first value of t at which a soluton will be @@ -252,12 +250,9 @@ private: //! Factor method for choosing a DAE solver /*! - * - * @param itype String identifying the type - * (IDA is the only option) - * @param f Residual function to be solved by the DAE algorithm - * - * @return Returns a point to the instantiated DAE_Solver object + * @param itype String identifying the type (IDA is the only option) + * @param f Residual function to be solved by the DAE algorithm + * @returns a point to the instantiated DAE_Solver object */ DAE_Solver* newDAE_Solver(const std::string& itype, ResidJacEval& f); diff --git a/include/cantera/numerics/DenseMatrix.h b/include/cantera/numerics/DenseMatrix.h index 0f2bc667b..22f1ee7e9 100644 --- a/include/cantera/numerics/DenseMatrix.h +++ b/include/cantera/numerics/DenseMatrix.h @@ -19,50 +19,52 @@ namespace Cantera /** * @defgroup numerics Numerical Utilities within Cantera * - * Cantera contains some capabilities for solving nonlinear equations and - * integrating both ODE and DAE equation systems in time. This section describes these - * capabilities. + * Cantera contains some capabilities for solving nonlinear equations and + * integrating both ODE and DAE equation systems in time. This section describes + * these capabilities. * */ -//! Exception thrown when an LAPACK error is encountered associated with inverting or solving a matrix +//! Exception thrown when an LAPACK error is encountered associated with +//! inverting or solving a matrix /*! - * A named error condition is used so that the calling code may differentiate this type of error - * from other error conditions. + * A named error condition is used so that the calling code may differentiate + * this type of error from other error conditions. */ class CELapackError : public CanteraError { public: //! Constructor passes through to main Cantera error handler /*! - * @param routine Name of calling routine - * @param msg Informative message + * @param routine Name of calling routine + * @param msg Informative message */ CELapackError(const std::string& routine, const std::string& msg) : CanteraError(routine + " LAPACK ERROR", msg) { } }; -//! A class for full (non-sparse) matrices with Fortran-compatible -//! data storage, which adds matrix operations to class Array2D. +//! A class for full (non-sparse) matrices with Fortran-compatible data storage, +//! which adds matrix operations to class Array2D. /*! - * The dense matrix class adds matrix operations onto the Array2D class. - * These matrix operations are carried out by the appropriate BLAS and LAPACK routines + * The dense matrix class adds matrix operations onto the Array2D class. These + * matrix operations are carried out by the appropriate BLAS and LAPACK routines * - * Error handling from BLAS and LAPACK are handled via the following formulation. - * Depending on a variable, a singular matrix or other terminal error condition from - * LAPACK is handled by either throwing an exception of type, CELapackError, or by - * returning the error code condition to the calling routine. + * Error handling from BLAS and LAPACK are handled via the following + * formulation. Depending on a variable, a singular matrix or other terminal + * error condition from LAPACK is handled by either throwing an exception of + * type, CELapackError, or by returning the error code condition to the calling + * routine. * - * The int variable, m_useReturnErrorCode, determines which method is used. - * The default value of zero means that an exception is thrown. A value of 1 - * means that a return code is used. + * The int variable, m_useReturnErrorCode, determines which method is used. The + * default value of zero means that an exception is thrown. A value of 1 means + * that a return code is used. * - * Reporting of these LAPACK error conditions is handled by the class variable - * m_printLevel. The default is for no reporting. If m_printLevel is nonzero, - * the error condition is reported to Cantera's log file. + * Reporting of these LAPACK error conditions is handled by the class variable + * m_printLevel. The default is for no reporting. If m_printLevel is nonzero, + * the error condition is reported to Cantera's log file. * - * @ingroup numerics + * @ingroup numerics */ class DenseMatrix : public Array2D { @@ -72,33 +74,24 @@ public: //! Constructor. /*! - * Create an \c n by \c m matrix, and initialize all elements to \c v. + * Create an \c n by \c m matrix, and initialize all elements to \c v. * - * @param n New number of rows - * @param m New number of columns - * @param v Default fill value. defaults to zero. + * @param n New number of rows + * @param m New number of columns + * @param v Default fill value. defaults to zero. */ DenseMatrix(size_t n, size_t m, doublereal v = 0.0); - //! Copy constructor - /*! - * @param y Object to be copied - */ DenseMatrix(const DenseMatrix& y); - - //! Assignment operator - /*! - * @param y Object to be copied - */ DenseMatrix& operator=(const DenseMatrix& y); //! Resize the matrix /*! - * Resize the matrix to n rows by m cols. + * Resize the matrix to n rows by m cols. * - * @param n New number of rows - * @param m New number of columns - * @param v Default fill value. defaults to zero. + * @param n New number of rows + * @param m New number of columns + * @param v Default fill value. defaults to zero. */ void resize(size_t n, size_t m, doublereal v = 0.0); @@ -106,11 +99,10 @@ public: //! Return a const vector of const pointers to the columns /*! - * Note, the Jacobian can not be altered by this routine, and - * therefore the member function is const. + * Note, the Jacobian can not be altered by this routine, and therefore the + * member function is const. * - * @return returns a vector of pointers to the top of the columns - * of the matrices. + * @returns a vector of pointers to the top of the columns of the matrices. */ const doublereal* const* const_colPts() const; @@ -118,29 +110,29 @@ public: //! Multiply A*B and write result to \c prod. /*! - * @param b input DenseMatrix B of size NxN - * @param prod output output DenseMatrix prod size NxN + * @param[in] b DenseMatrix B of size NxN + * @param[out] prod DenseMatrix prod size NxN */ virtual void mult(const DenseMatrix& b, DenseMatrix& prod) const; //! Left-multiply the matrix by transpose(b), and write the result to prod. /*! - * @param b left multiply by this vector. The length must be equal to n - * the number of rows in the matrix. - * @param prod Resulting vector. This is of length m, the number of columns - * in the matrix + * @param b left multiply by this vector. The length must be equal to n + * the number of rows in the matrix. + * @param prod Resulting vector. This is of length m, the number of columns + * in the matrix */ virtual void leftMult(const double* const b, double* const prod) const; //! Return a changeable value of the pivot vector /*! - * @return Returns a reference to the pivot vector as a vector_int + * @returns a reference to the pivot vector as a vector_int */ vector_int& ipiv(); //! Return a changeable value of the pivot vector /*! - * @return Returns a reference to the pivot vector as a vector_int + * @returns a reference to the pivot vector as a vector_int */ const vector_int& ipiv() const { return m_ipiv; @@ -156,22 +148,25 @@ protected: public: //! Error Handling Flag /*! - * The default is to set this to 0. In this case, if a factorization is requested and can't be achieved, - * a CESingularMatrix exception is triggered. No return code is used, because an exception is thrown. - * If this is set to 1, then an exception is not thrown. Routines return with an error code, that is up - * to the calling routine to handle correctly. Negative return codes always throw an exception. + * The default is to set this to 0. In this case, if a factorization is + * requested and can't be achieved, a CESingularMatrix exception is + * triggered. No return code is used, because an exception is thrown. If + * this is set to 1, then an exception is not thrown. Routines return with + * an error code, that is up to the calling routine to handle correctly. + * Negative return codes always throw an exception. */ int m_useReturnErrorCode; //! Print Level /*! - * Printing is done to the log file using the routine writelogf(). + * Printing is done to the log file using the routine writelogf(). * - * Level of printing that is carried out. Only error conditions are printed out, if this value is nonzero. + * Level of printing that is carried out. Only error conditions are printed + * out, if this value is nonzero. */ int m_printLevel; - // Listing of friend functions which are defined below + // Listing of friend functions which are defined below friend int solve(DenseMatrix& A, double* b, size_t nrhs, size_t ldb); friend int solve(DenseMatrix& A, DenseMatrix& b); @@ -181,60 +176,61 @@ public: //! Solve Ax = b. Array b is overwritten on exit with x. /*! - * The solve class uses the LAPACK routine dgetrf to invert the m xy n matrix. + * The solve function uses the LAPACK routine dgetrf to invert the m xy n matrix. + * + * The factorization has the form * - * The factorization has the form * A = P * L * U - * where P is a permutation matrix, L is lower triangular with unit - * diagonal elements (lower trapezoidal if m > n), and U is upper - * triangular (upper trapezoidal if m < n). * - * The system is then solved using the LAPACK routine dgetrs + * where P is a permutation matrix, L is lower triangular with unit diagonal + * elements (lower trapezoidal if m > n), and U is upper triangular (upper + * trapezoidal if m < n). * - * @param A Dense matrix to be factored - * @param b RHS(s) to be solved. - * @param nrhs Number of right hand sides to solve - * @param ldb Leading dimension of b, if nrhs > 1 + * The system is then solved using the LAPACK routine dgetrs + * + * @param A Dense matrix to be factored + * @param b RHS(s) to be solved. + * @param nrhs Number of right hand sides to solve + * @param ldb Leading dimension of b, if nrhs > 1 */ int solve(DenseMatrix& A, double* b, size_t nrhs=1, size_t ldb=0); -//! Solve Ax = b for multiple right-hand-side vectors. +//! Solve Ax = b for multiple right-hand-side vectors. /*! - * @param A Dense matrix to be factored - * @param b Dense matrix of RHS's. Each column is a RHS + * @param A Dense matrix to be factored + * @param b Dense matrix of RHS's. Each column is a RHS */ int solve(DenseMatrix& A, DenseMatrix& b); //! Multiply \c A*b and return the result in \c prod. Uses BLAS routine DGEMV. /*! - * \f[ - * prod_i = sum^N_{j = 1}{A_{ij} b_j} - * \f] + * \f[ + * prod_i = sum^N_{j = 1}{A_{ij} b_j} + * \f] * - * @param A input Dense Matrix A with M rows and N columns - * @param b input vector b with length N - * @param prod output output vector prod length = M + * @param[in] A Dense Matrix A with M rows and N columns + * @param[in] b vector b with length N + * @param[out] prod vector prod length = M */ void multiply(const DenseMatrix& A, const double* const b, double* const prod); //! Multiply \c A*b and add it to the result in \c prod. Uses BLAS routine DGEMV. /*! - * \f[ - * prod_i += sum^N_{j = 1}{A_{ij} b_j} - * \f] + * \f[ + * prod_i += sum^N_{j = 1}{A_{ij} b_j} + * \f] * - * @param A input Dense Matrix A with M rows and N columns - * @param b input vector b with length N - * @param prod output output vector prod length = M + * @param[in] A Dense Matrix A with M rows and N columns + * @param[in] b vector b with length N + * @param[out] prod vector prod length = M */ void increment(const DenseMatrix& A, const double* const b, double* const prod); //! invert A. A is overwritten with A^-1. /*! - * @param A Invert the matrix A and store it back in place - * - * @param nn Size of A. This defaults to -1, which means that the number - * of rows is used as the default size of n + * @param A Invert the matrix A and store it back in place + * @param nn Size of A. This defaults to -1, which means that the number of + * rows is used as the default size of n */ int invert(DenseMatrix& A, size_t nn=npos); diff --git a/include/cantera/numerics/Func1.h b/include/cantera/numerics/Func1.h index c5505223c..41c087220 100644 --- a/include/cantera/numerics/Func1.h +++ b/include/cantera/numerics/Func1.h @@ -50,8 +50,8 @@ public: //! Duplicate the current function. /*! - * This duplicates the current function, returning a - * reference to the new malloced function. + * This duplicates the current function, returning a reference to the new + * malloced function. */ virtual Func1& duplicate() const; @@ -65,17 +65,16 @@ public: //! Creates a derivative to the current function /*! - * This will malloc a derivative function and - * return a reference to the function. + * This will malloc a derivative function and return a reference to the + * function. */ virtual Func1& derivative() const; //! Routine to determine if two functions are the same. /*! - * Two functions are the same if they are the same function. - * This means that the ID and stored constant is the same. - * This means that the m_f1 and m_f2 are identical if they - * are non-null. + * Two functions are the same if they are the same function. This means + * that the ID and stored constant is the same. This means that the m_f1 + * and m_f2 are identical if they are non-null. */ bool isIdentical(Func1& other) const; @@ -987,8 +986,7 @@ protected: }; /** - * Periodic function. Takes any function and makes it - * periodic with period T. + * Periodic function. Takes any function and makes it periodic with period T. */ class Periodic1 : public Func1 { diff --git a/include/cantera/numerics/GeneralMatrix.h b/include/cantera/numerics/GeneralMatrix.h index 4871c67cb..ff0689230 100644 --- a/include/cantera/numerics/GeneralMatrix.h +++ b/include/cantera/numerics/GeneralMatrix.h @@ -32,20 +32,15 @@ public: */ GeneralMatrix(int matType); - //! Copy Constructor GeneralMatrix(const GeneralMatrix& right); - - //! Assignment operator GeneralMatrix& operator=(const GeneralMatrix& right); - - //! Destructor. Does nothing. virtual ~GeneralMatrix() {} //! Duplicator member function /*! - * This function will duplicate the matrix given a generic GeneralMatrix + * This function will duplicate the matrix given a generic GeneralMatrix * - * @return Returns a pointer to the malloced object + * @returns a pointer to the malloced object */ virtual GeneralMatrix* duplMyselfAsGeneralMatrix() const = 0; @@ -54,21 +49,21 @@ public: //! Multiply A*b and write result to prod. /*! - * @param b Vector to do the rh multiplication - * @param prod OUTPUT vector to receive the result + * @param b Vector to do the rh multiplication + * @param prod OUTPUT vector to receive the result */ virtual void mult(const doublereal* b, doublereal* prod) const = 0; //! Multiply b*A and write result to prod. /*! - * @param b Vector to do the lh multiplication - * @param prod OUTPUT vector to receive the result + * @param b Vector to do the lh multiplication + * @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. + /*! + * We flip m_factored boolean to indicate that the matrix is now A-1. */ virtual int factor() = 0; @@ -76,7 +71,7 @@ public: /*! * we set m_factored to 2 to indicate the matrix is now QR factored * - * @return Returns the info variable from LAPACK + * @returns the info variable from LAPACK */ virtual int factorQR() { throw NotImplementedError("GeneralMatrix::factorQR"); @@ -84,9 +79,9 @@ public: //! Returns an estimate of the inverse of the condition number for the matrix /*! - * The matrix must have been previously factored using the QR algorithm + * The matrix must have been previously factored using the QR algorithm * - * @return returns the inverse of the condition number + * @returns the inverse of the condition number */ virtual doublereal rcondQR() { throw NotImplementedError("GeneralMatrix::rcondQR"); @@ -94,10 +89,10 @@ public: //! Returns an estimate of the inverse of the condition number for the matrix /*! - * The matrix must have been previously factored using the LU algorithm + * 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 + * @returns the inverse of the condition number */ virtual doublereal rcond(doublereal a1norm) = 0; @@ -122,8 +117,9 @@ public: //! Return the size and structure of the matrix /*! - * @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. + * @param iStruct OUTPUT Pointer to a vector of ints that describe the + * structure of the matrix. + * @returns the number of rows and columns in the matrix. */ virtual size_t nRowsAndStruct(size_t* const iStruct = 0) const = 0; @@ -134,10 +130,10 @@ public: //! Solves the Ax = b system returning x in the b spot. /*! - * @param b Vector for the RHS of the equation system - * @param nrhs Number of right-hand sides to solve, default 1 - * @param ldb Leading dimension of the right-hand side array. - * Defaults to nRows() + * @param b Vector for the RHS of the equation system + * @param nrhs Number of right-hand sides to solve, default 1 + * @param ldb Leading dimension of the right-hand side array. Defaults to + * nRows() */ virtual int solve(doublereal* b, size_t nrhs=1, size_t ldb=0) = 0; @@ -146,57 +142,55 @@ public: return (m_factored != 0); } - //! Return a pointer to the top of column j, columns are assumed to be contiguous in memory + //! 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 + * @param j Value of the column + * @returns a pointer to the top of the column */ virtual doublereal* ptrColumn(size_t j) = 0; //! Index into the (i,j) element /*! - * @param i row - * @param j column - * - * Returns a changeable reference to the matrix entry + * @param i row + * @param j column + * @returns a changeable reference to the matrix entry */ virtual doublereal& operator()(size_t i, size_t j) = 0; //! Constant Index into the (i,j) element /*! - * @param i row - * @param j column - * - * Returns an unchangeable reference to the matrix entry + * @param i row + * @param j column + * @returns an unchangeable reference to the matrix entry */ virtual doublereal operator()(size_t i, size_t j) const = 0; //! Return an iterator pointing to the first element /*! - * We might drop this later + * 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 + * 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. + * 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. + * @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 + * 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 @@ -205,8 +199,8 @@ public: //! 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 + * 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 @@ -221,8 +215,8 @@ public: int matrixType_; protected: - //! Indicates whether the matrix is factored. 0 for unfactored; Non-zero values - //! indicate a particular factorization (LU=1, QR=2). + //! Indicates whether the matrix is factored. 0 for unfactored; Non-zero + //! values indicate a particular factorization (LU=1, QR=2). int m_factored; }; diff --git a/include/cantera/numerics/IDA_Solver.h b/include/cantera/numerics/IDA_Solver.h index 62a26a02d..36f4f2ec1 100644 --- a/include/cantera/numerics/IDA_Solver.h +++ b/include/cantera/numerics/IDA_Solver.h @@ -42,9 +42,10 @@ class IDA_Solver : public DAE_Solver public: //! Constructor. /*! - * Default settings: dense Jacobian, no user-supplied Jacobian function, Newton iteration. + * Default settings: dense Jacobian, no user-supplied Jacobian function, + * Newton iteration. * - * @param f Function that will supply the time dependent residual to be solved + * @param f Function that will supply the time dependent residual to be solved */ IDA_Solver(ResidJacEval& f); @@ -62,8 +63,8 @@ public: //! Set up the problem to use a band solver /*! - * @param m_upper upper band width of the matrix - * @param m_lower lower band width of the matrix + * @param m_upper upper band width of the matrix + * @param m_lower lower band width of the matrix */ virtual void setBandedLinearSolver(int m_upper, int m_lower); @@ -71,34 +72,34 @@ public: //! Set the maximum number of time steps /*! - * @param n input of maximum number of time steps + * @param n input of maximum number of time steps */ virtual void setMaxNumSteps(int n); //! Sset the initial step size /*! - * @param h0 initial step size value + * @param h0 initial step size value */ virtual void setInitialStepSize(doublereal h0); //! Set the stop time /*! - * @param tstop the independent variable value past which the solution is not to proceed. + * @param tstop the independent variable value past which the solution is + * not to proceed. */ virtual void setStopTime(doublereal tstop); //! Get the current step size from IDA via a call /*! - * @return Returns the current step size. + * @returns the current step size. */ virtual double getCurrentStepFromIDA(); - //! Set the form of the Jacobian /*! - * @param formJac Form of the Jacobian - * 0 numerical Jacobian - * 1 analytical Jacobian given by the evalJacobianDP() function + * @param formJac Form of the Jacobian + * 0 numerical Jacobian + * 1 analytical Jacobian given by the evalJacobianDP() function */ virtual void setJacobianType(int formJac); @@ -113,7 +114,7 @@ public: //! Set the maximum number of nonlinear solver convergence failures /*! - * @param n Value of nonlin failures. If value is exceeded, the calculation terminates. + * @param n Value of nonlin failures. If value is exceeded, the calculation terminates. */ virtual void setMaxNonlinConvFailures(int n); @@ -132,12 +133,12 @@ public: //! Step the system to a final value of the time /*! - * @param tout Final value of the time - * @return Returns the IDASolve() return flag + * @param tout Final value of the time + * @returns the IDASolve() return flag * - * The return values for IDASolve are described below. - * (The numerical return values are defined above in this file.) - * All unsuccessful returns give a negative return value. + * The return values for IDASolve are described below. (The numerical return + * values are defined above in this file.) All unsuccessful returns give a + * negative return value. * * IDA_SUCCESS * IDASolve succeeded and no roots were found. @@ -226,7 +227,7 @@ protected: //! Initial value of the time doublereal m_t0; - //! Current value of the solution vector + //! Current value of the solution vector N_Vector m_y; //! Current value of the derivative of the solution vector diff --git a/include/cantera/numerics/Integrator.h b/include/cantera/numerics/Integrator.h index f90adf66e..6b2eb150a 100644 --- a/include/cantera/numerics/Integrator.h +++ b/include/cantera/numerics/Integrator.h @@ -44,7 +44,6 @@ enum IterType { Functional_Iter }; - //! Abstract base class for ODE system integrators. /*! * @ingroup odeGroup @@ -60,7 +59,7 @@ public: virtual ~Integrator() { } - //! Set error tolerances. + //! Set error tolerances. /*! * @param reltol scalar relative tolerance * @param n Number of equations @@ -71,7 +70,7 @@ public: warn("setTolerances"); } - //! Set error tolerances. + //! Set error tolerances. /*! * @param reltol scalar relative tolerance * @param abstol scalar absolute tolerance @@ -80,7 +79,7 @@ public: warn("setTolerances"); } - //! Set the sensitivity error tolerances + //! Set the sensitivity error tolerances /*! * @param reltol scalar relative tolerance * @param abstol scalar absolute tolerance @@ -88,7 +87,7 @@ public: virtual void setSensitivityTolerances(doublereal reltol, doublereal abstol) { } - //! Set the problem type. + //! Set the problem type. /*! * @param probtype Type of the problem */ @@ -97,8 +96,8 @@ public: } /** - * Initialize the integrator for a new problem. Call after - * all options have been set. + * Initialize the integrator for a new problem. Call after all options have + * been set. * @param t0 initial time * @param func RHS evaluator object for system of equations. */ diff --git a/include/cantera/numerics/ResidEval.h b/include/cantera/numerics/ResidEval.h index 08285558a..9f9ad7018 100644 --- a/include/cantera/numerics/ResidEval.h +++ b/include/cantera/numerics/ResidEval.h @@ -1,6 +1,4 @@ -/** - * @file ResidEval.h - */ +//! @file ResidEval.h // Copyright 2006 California Institute of Technology @@ -60,9 +58,8 @@ public: } /** - * Specify that solution component k is purely algebraic - - * that is, the derivative of this component does not appear - * in the residual function. + * Specify that solution component k is purely algebraic - that is, the + * derivative of this component does not appear in the residual function. */ virtual void setAlgebraic(const int k) { if ((int) m_alg.size() < (k+1)) { @@ -108,11 +105,11 @@ public: /*! * Values for both the solution and the value of ydot may be provided. * - * @param[in] t0 Time - * @param[out] y Solution vector - * @param[out] ydot Rate of change of solution vector. + * @param[in] t0 Time + * @param[out] y Solution vector + * @param[out] ydot Rate of change of solution vector. * - * @return Returns a flag to indicate that operation is successful. + * @returns a flag to indicate that operation is successful. * 1 Means a successful operation * -0 or neg value Means an unsuccessful operation */ @@ -126,10 +123,10 @@ public: //! Return the number of equations in the equation system virtual int nEquations() const = 0; - //! Write out to a file or to standard output the current solution + //! Write out to a file or to standard output the current solution /*! - * ievent is a description of the event that caused this - * function to be called. + * ievent is a description of the event that caused this function to be + * called. */ virtual void writeSolution(int ievent, const double time, const double deltaT, @@ -151,10 +148,10 @@ public: } } - //! Return the number of parameters in the calculation + //! Return the number of parameters in the calculation /*! - * This is the number of parameters in the sensitivity calculation. We have - * set this to zero and have included it for later expansion + * This is the number of parameters in the sensitivity calculation. We have + * set this to zero and have included it for later expansion */ int nparams() const { return 0; @@ -163,8 +160,8 @@ public: protected: //! Mapping vector that stores whether a degree of freedom is a DAE or not /*! - * The first index is the equation number. The second index is 1 if it is a DAE, - * and zero if it is not. + * The first index is the equation number. The second index is 1 if it is a + * DAE, and zero if it is not. */ vector_int m_alg; std::map m_constrain; diff --git a/include/cantera/numerics/ResidJacEval.h b/include/cantera/numerics/ResidJacEval.h index f4624baf8..4bbf1450b 100644 --- a/include/cantera/numerics/ResidJacEval.h +++ b/include/cantera/numerics/ResidJacEval.h @@ -35,42 +35,39 @@ enum ResidEval_Type_Enum { Base_ShowSolution, //! Base residual calculation containing any lagged components /*! - * We use this to calculate residuals when doing line searches along - * directions determined by Jacobians that are missing contributions - * from lagged entries. + * We use this to calculate residuals when doing line searches along + * irections determined by Jacobians that are missing contributions from + * lagged entries. */ Base_LaggedSolutionComponents }; //! Wrappers for the function evaluators for Nonlinear solvers and Time steppers /*! - * A class for full (non-sparse dense matrices with Fortran-compatible data storage. - * The class adds support for identifying what types of calls are made to the residual - * evaluator by adding the ResidEval_Type_Enum class. + * A class for full (non-sparse dense matrices with Fortran-compatible data + * storage. The class adds support for identifying what types of calls are made + * to the residual evaluator by adding the ResidEval_Type_Enum class. */ class ResidJacEval : public ResidEval { public: //!Default constructor /*! - * @param atol Initial value of the global tolerance (defaults to 1.0E-13) + * @param atol Initial value of the global tolerance (defaults to 1.0E-13) */ ResidJacEval(doublereal atol = 1.0e-13); - //!Copy Constructor ResidJacEval(const ResidJacEval& right); - - //! Assignment operator ResidJacEval& operator=(const ResidJacEval& right); //! Duplication routine for objects derived from residJacEval /*! - * This virtual routine can be used to duplicate objects which inherit - * from ResidJacEval even if the application only has a pointer to - * ResidJacEval to work with. + * This virtual routine can be used to duplicate objects which inherit from + * ResidJacEval even if the application only has a pointer to ResidJacEval + * to work with. * - * These routines are basically wrappers around the derived copy - * constructor. + * These routines are basically wrappers around the derived copy + * constructor. */ virtual ResidJacEval* duplMyselfAsResidJacEval() const; @@ -90,7 +87,7 @@ public: * differenced or that the residual doesn't take this issue into account) * @param delta_x Value of the delta used in the numerical differencing * - * @return Returns a flag to indicate that operation is successful. + * @returns a flag to indicate that operation is successful. * 1 Means a successful operation * -0 or neg value Means an unsuccessful operation */ @@ -110,8 +107,8 @@ public: //! Filter the solution predictions /*! - * Codes might provide a predicted step change. This routine filters the predicted - * solution vector eliminating illegal directions. + * Codes might provide a predicted step change. This routine filters the + * predicted solution vector eliminating illegal directions. * * @param t Time (input) * @param ybase Solution vector (input, output) @@ -123,12 +120,12 @@ public: //! Filter the solution predictions /*! - * Codes might provide a predicted solution vector. This routine filters the predicted - * solution vector. + * Codes might provide a predicted solution vector. This routine filters the + * predicted solution vector. * * @param t Time (input) * @param y Solution vector (input, output) - * @return Return the norm of the amount of filtering + * @returns the norm of the amount of filtering */ virtual doublereal filterSolnPrediction(const doublereal t, doublereal* const y); @@ -140,15 +137,16 @@ public: //! Evaluate the time tracking equations, if any /*! - * Evaluate time integrated quantities that are calculated at the - * end of every successful time step. This call is made once at the end of every successful - * time step that advances the time. It's also made once at the start of the time stepping. + * Evaluate time integrated quantities that are calculated at the end of + * every successful time step. This call is made once at the end of every + * successful time step that advances the time. It's also made once at the + * start of the time stepping. * * @param t Time (input) * @param delta_t The current value of the time step (input) * @param y Solution vector (input, do not modify) * @param ydot Rate of change of solution vector. (input, do not modify) - * @return Returns a flag to indicate that operation is successful. + * @returns a flag to indicate that operation is successful. * 1 Means a successful operation * -0 or neg value Means an unsuccessful operation */ @@ -157,11 +155,12 @@ public: //! Evaluate any stopping criteria other than a final time limit /*! - * If we are to stop the time integration for any reason other than reaching a final time limit, tout, - * provide a test here. This call is made at the end of every successful time step iteration + * If we are to stop the time integration for any reason other than reaching + * a final time limit, tout, provide a test here. This call is made at the + * end of every successful time step iteration * - * @return If true, the the time stepping is stopped. If false, then time stepping is stopped if t >= tout - * Defaults to false. + * @return If true, the the time stepping is stopped. If false, then time + * stepping is stopped if t >= tout Defaults to false. * @param t Time (input) * @param delta_t The current value of the time step (input) * @param y Solution vector (input, do not modify) @@ -174,17 +173,17 @@ public: //! Return a vector of delta y's for calculation of the numerical Jacobian /*! - * There is a default algorithm provided. + * There is a default algorithm provided. * - * delta_y[i] = atol[i] + 1.0E-6 ysoln[i] - * delta_y[i] = atol[i] + MAX(1.0E-6 ysoln[i] * 0.01 * solnWeights[i]) + * delta_y[i] = atol[i] + 1.0E-6 ysoln[i] + * delta_y[i] = atol[i] + MAX(1.0E-6 ysoln[i] * 0.01 * solnWeights[i]) * * @param t Time (input) * @param y Solution vector (input, do not modify) * @param ydot Rate of change of solution vector. (input, do not modify) * @param delta_y Value of the delta to be used in calculating the numerical Jacobian * @param solnWeights Value of the solution weights that are used in determining convergence (default = 0) - * @return Returns a flag to indicate that operation is successful. + * @returns a flag to indicate that operation is successful. * 1 Means a successful operation * -0 or neg value Means an unsuccessful operation */ @@ -195,9 +194,10 @@ public: doublereal* const delta_y, const doublereal* const solnWeights = 0); - //! Returns a vector of column scale factors that can be used to column scale Jacobians. + //! Returns a vector of column scale factors that can be used to column + //! scale Jacobians. /*! - * Default to yScales[] = 1.0 + * Default to yScales[] = 1.0 * * @param t Time (input) * @param y Solution vector (input, do not modify) @@ -207,7 +207,8 @@ public: virtual void calcSolnScales(const doublereal t, const doublereal* const y, const doublereal* const y_old, doublereal* const yScales); - //! This function may be used to create output at various points in the execution of an application. + //! This function may be used to create output at various points in the + //! execution of an application. /*! * @param ifunc identity of the call * 0 Initial call @@ -224,7 +225,8 @@ public: const doublereal* const y, const doublereal* const ydot); - //! This function may be used to create output at various points in the execution of an application. + //! This function may be used to create output at various points in the + //! execution of an application. /*! * This routine calls user_out2(). * @@ -239,27 +241,28 @@ public: //! Multiply the matrix by another matrix that leads to better conditioning /*! - * Provide a left sided matrix that will multiply the current Jacobian, after scaling - * and lead to a better conditioned system. - * This routine is called just before the matrix is factored. + * Provide a left sided matrix that will multiply the current Jacobian, + * after scaling and lead to a better conditioned system. This routine is + * called just before the matrix is factored. * - * Original Problem: + * Original Problem: * J delta_x = - Resid * - * New problem: + * New problem: * M (J delta_x) = - M Resid * - * @param matrix Pointer to the current Jacobian (if zero, it's already been factored) - * @param nrows offsets for the matrix - * @param rhs residual vector. This also needs to be LHS multiplied by M - * @return Returns a flag to indicate that operation is successful. + * @param matrix Pointer to the current Jacobian (if zero, it's already been factored) + * @param nrows offsets for the matrix + * @param rhs residual vector. This also needs to be LHS multiplied by M + * @returns a flag to indicate that operation is successful. * 1 Means a successful operation * -0 or neg value Means an unsuccessful operation */ 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 * @@ -270,7 +273,7 @@ public: * @param ydot Rate of change of solution vector. (input, do not modify) * @param J Reference to the SquareMatrix object to be calculated (output) * @param resid Value of the residual that is computed (output) - * @return Returns a flag to indicate that operation is successful. + * @returns a flag to indicate that operation is successful. * 1 Means a successful operation * -0 or neg value Means an unsuccessful operation */ @@ -290,7 +293,7 @@ public: * @param jacobianColPts Pointer to the vector of pts to columns of the SquareMatrix * object to be calculated (output) * @param resid Value of the residual that is computed (output) - * @return Returns a flag to indicate that operation is successful. + * @returns a flag to indicate that operation is successful. * 1 Means a successful operation * -0 or neg value Means an unsuccessful operation */ diff --git a/include/cantera/numerics/RootFind.h b/include/cantera/numerics/RootFind.h index 9da61e89b..fa1af6ee0 100644 --- a/include/cantera/numerics/RootFind.h +++ b/include/cantera/numerics/RootFind.h @@ -27,105 +27,109 @@ namespace Cantera //! This means that the root solver was a success #define ROOTFIND_SUCCESS 0 -//! This return value means that the root finder resolved a solution in the x coordinate -//! However, convergence in F was not achieved. +//! This return value means that the root finder resolved a solution in the x +//! coordinate, however, convergence in F was not achieved. /*! - * A common situation for this to happen is that f(x) is discontinuous about f(x) = f_0, - * where we seek the x where the function is equal to f_0. f(x) spans the - * f_0 while not being equal to f_0 anywhere. + * A common situation for this to happen is that f(x) is discontinuous about + * f(x) = f_0, where we seek the x where the function is equal to f_0. f(x) + * spans the f_0 while not being equal to f_0 anywhere. */ #define ROOTFIND_SUCCESS_XCONVERGENCEONLY 1 -//! This means that the root solver failed to achieve convergence +//! This means that the root solver failed to achieve convergence #define ROOTFIND_FAILEDCONVERGENCE -1 -//! This means that the input to the root solver was defective +//! This means that the input to the root solver was defective #define ROOTFIND_BADINPUT -2 -//! This means that the rootfinder believes the solution is lower than xmin +//! This means that the rootfinder believes the solution is lower than xmin #define ROOTFIND_SOLNLOWERTHANXMIN -3 -//! This means that the rootfinder believes the solution is higher than xmax +//! This means that the rootfinder believes the solution is higher than xmax #define ROOTFIND_SOLNHIGHERTHANXMAX -4 //@} //! Root finder for 1D problems /*! - * The root finder solves a single nonlinear equation described below. + * The root finder solves a single nonlinear equation described below. * - * \f[ - * f(x) = f_0 - * \f] + * \f[ + * f(x) = f_0 + * \f] * - * \f$ f(x) \f$ is assumed to be single valued as a function of x.\f$ f(x) \f$ is not assumed to be continuous nor is - * its derivative assumed to be well formed. + * \f$ f(x) \f$ is assumed to be single valued as a function of x.\f$ f(x) \f$ + * is not assumed to be continuous nor is its derivative assumed to be well + * formed. * - * Root finders are significantly different in the sense that do not have to rely - * solely on Newton's method to find the answer to the problem. Instead they use a method to bound - * the solution between high and low values and then use a method to refine that bound. The eventual - * solution to the problem is presented as x_best and as a bound, delta_X, on the solution - * component. Because of this, they are far more stable for functions and Jacobians that have discontinuities - * or noise associated with them. + * Root finders are significantly different in the sense that do not have to + * rely solely on Newton's method to find the answer to the problem. Instead + * they use a method to bound the solution between high and low values and then + * use a method to refine that bound. The eventual solution to the problem is + * presented as x_best and as a bound, delta_X, on the solution component. + * Because of this, they are far more stable for functions and Jacobians that + * have discontinuities or noise associated with them. * - * The algorithm is a convolution of a local Secant method with an approach of finding a straddle in x. - * The Jacobian is never required. + * The algorithm is a convolution of a local Secant method with an approach of + * finding a straddle in x. The Jacobian is never required. * - * There is a general breakdown of the algorithm into stages. The first stage seeks to find a straddle of the - * function. The second stage seeks to reduce the bounds in x and f in order to satisfy the specification of the - * stopping criteria. In the last stage the algorithm seeks to find the base value of x that satisfies the - * original equation given what it current knows about the function. + * There is a general breakdown of the algorithm into stages. The first stage + * seeks to find a straddle of the function. The second stage seeks to reduce + * the bounds in x and f in order to satisfy the specification of the stopping + * criteria. In the last stage the algorithm seeks to find the base value of x + * that satisfies the original equation given what it current knows about the + * function. * - * Globalization strategy + * Globalization strategy * - * Specifying the General Changes in x + * Specifying the General Changes in x * - * Supplying Hints with General Function Behavior Flags + * Supplying Hints with General Function Behavior Flags * - * Stopping Criteria + * Stopping Criteria * - * Specification of the Stopping Criteria + * Specification of the Stopping Criteria * - * Additional constraints + * Additional constraints * - * Bounds Criteria For the Routine + * Bounds Criteria For the Routine * - * Example + * Example * - * @code - * // Define a residual. The definition of a residual involves a lot more work than is shown here. - * ResidEval * ec; - * // Instantiate the root finder with the residual to be solved, ec. - * RootFind rf(&ec); - * // Set the relative and absolute tolerancess for f and x. - * rf.setTol(1.0E-5, 1.0E-10, 1.0E-5, 1.0E-11); - * // Give a hint about the function's dependence on x. This is needed, for example, if the function has - * // flat regions. - * rf.setFuncIsGenerallyIncreasing(true); - * rf.setDeltaX(0.01); - * // Supply an initial guess for the solution - * double xbest = phiM; - * double oldP = printLvl_; - * // Set the print level for the solver. Zero produces no output. Two produces a summary table of each iteration. - * rf.setPrintLvl(2); - * // Define a minimum and maximum for the independent variable. - * double phimin = 1.3; - * double phimax = 2.2; - * // Define a maximum iteration number - * int itmax = 100; - * // Define the f_0 value, and on return will contain the actual value of f(x) obtained - * double currentObtained; - * // Call the solver - * status = rf.solve(phimin, phimax, 100, currentObtained, &xbest); - * if (status == 0) { - * if (printLvl_ > 1) { - * printf("Electrode::integrateConstantCurrent(): Volts (%g amps) = %g\n", currentObtained, xbest); - * } - * } else { - * if (printLvl_) { - * printf("Electrode::integrateConstantCurrent(): bad status = %d Volts (%g amps) = %g\n", - * status, currentObtained, xbest); - * } - * } - * @endcode + * @code + * // Define a residual. The definition of a residual involves a lot more work than is shown here. + * ResidEval * ec; + * // Instantiate the root finder with the residual to be solved, ec. + * RootFind rf(&ec); + * // Set the relative and absolute tolerancess for f and x. + * rf.setTol(1.0E-5, 1.0E-10, 1.0E-5, 1.0E-11); + * // Give a hint about the function's dependence on x. This is needed, for example, if the function has + * // flat regions. + * rf.setFuncIsGenerallyIncreasing(true); + * rf.setDeltaX(0.01); + * // Supply an initial guess for the solution + * double xbest = phiM; + * double oldP = printLvl_; + * // Set the print level for the solver. Zero produces no output. Two produces a summary table of each iteration. + * rf.setPrintLvl(2); + * // Define a minimum and maximum for the independent variable. + * double phimin = 1.3; + * double phimax = 2.2; + * // Define a maximum iteration number + * int itmax = 100; + * // Define the f_0 value, and on return will contain the actual value of f(x) obtained + * double currentObtained; + * // Call the solver + * status = rf.solve(phimin, phimax, 100, currentObtained, &xbest); + * if (status == 0) { + * if (printLvl_ > 1) { + * printf("Electrode::integrateConstantCurrent(): Volts (%g amps) = %g\n", currentObtained, xbest); + * } + * } else { + * if (printLvl_) { + * printf("Electrode::integrateConstantCurrent(): bad status = %d Volts (%g amps) = %g\n", + * status, currentObtained, xbest); + * } + * } + * @endcode * - * @todo Noise - * @todo General Search to be done when all else fails + * @todo Noise + * @todo General Search to be done when all else fails */ class RootFind { @@ -136,19 +140,15 @@ public: */ RootFind(ResidEval* resid); - //! Copy constructor RootFind(const RootFind& r); - ~RootFind() {} - - //! Assignment operator RootFind& operator=(const RootFind& right); private: //! Calculate a deltaX from an input value of x /*! - * This routine ensure that the deltaX will be greater or equal to DeltaXNorm_ - * or 1.0E-14 x + * This routine ensure that the deltaX will be greater or equal to + * DeltaXNorm_ or 1.0E-14 x * * @param x1 input value of x */ @@ -156,8 +156,8 @@ private: //! Calculate a deltaX from an input value of x /*! - * This routine ensure that the deltaX will be greater or equal to DeltaXNorm_ - * or 1.0E-14 x or deltaXConverged_. + * This routine ensure that the deltaX will be greater or equal to + * DeltaXNorm_ or 1.0E-14 x or deltaXConverged_. * * @param x1 input value of x */ @@ -165,49 +165,50 @@ private: //! Calculate a controlled, nonzero delta between two numbers /*! - * The delta is designed to be greater than or equal to delXMeaningful(x) defined above - * with the same sign as the original delta. Therefore if you subtract it from either - * of the two original numbers, you get a different number. + * The delta is designed to be greater than or equal to delXMeaningful(x) + * defined above with the same sign as the original delta. Therefore if you + * subtract it from either of the two original numbers, you get a different + * number. * - * @param x2 first number - * @param x1 second number + * @param x2 first number + * @param x1 second number */ doublereal deltaXControlled(doublereal x2, doublereal x1) const; //! Function to decide whether two real numbers are the same or not /*! - * A comparison is made between the two numbers to decide whether they - * are close to one another. This is defined as being within factor * delXMeaningful() of each other. + * A comparison is made between the two numbers to decide whether they are + * close to one another. This is defined as being within factor * + * delXMeaningful() of each other. * - * The basic premise here is that if the two numbers are too close, the noise - * will prevent an accurate calculation of the function and its slope. + * The basic premise here is that if the two numbers are too close, the + * noise will prevent an accurate calculation of the function and its slope. * * @param x1 First number * @param x2 second number * @param factor Multiplicative factor to multiple deltaX with - * @return Returns a boolean indicating whether the two numbers are the same or not. + * @returns a boolean indicating whether the two numbers are the same or not. */ bool theSame(doublereal x2, doublereal x1, doublereal factor = 1.0) const; public: - //! Using a line search method, find the root of a 1D function + //! Using a line search method, find the root of a 1D function /*! - * This routine solves the following equation. + * This routine solves the following equation. * - * \f[ - * R(x) = f(x) - f_o = 0 - * \f] + * \f[ + * R(x) = f(x) - f_o = 0 + * \f] * - * @param xmin Minimum value of x to be used. - * @param xmax Maximum value of x to be used - * @param itmax maximum number of iterations. Usually, it can be less than 50. - * @param funcTargetValue - * Value of \f$ f_o \f$ in the equation. - * On return, it contains the value of the function actually obtained. - * @param xbest Returns the x that satisfies the function - * On input, xbest should contain the best estimate of the solution. - * An attempt to find the solution near xbest is made. - * @return: + * @param xmin Minimum value of x to be used. + * @param xmax Maximum value of x to be used + * @param itmax maximum number of iterations. Usually, it can be less than 50. + * @param funcTargetValue Value of \f$ f_o \f$ in the equation. On return, + * it contains the value of the function actually obtained. + * @param xbest Returns the x that satisfies the function On input, xbest + * should contain the best estimate of the solution. An + * attempt to find the solution near xbest is made. + * @return: * 0 = ROOTFIND_SUCCESS Found function * -1 = ROOTFIND_FAILEDCONVERGENCE Failed to find the answer * -2 = ROOTFIND_BADINPUT Bad input was detected @@ -222,37 +223,36 @@ public: * R(x) = f(x) - f_o = 0 * \f] * - * @param x Value of the independent variable + * @param x Value of the independent variable * - * @return The routine returns the value of \f$ R(x) \f$ + * @return The routine returns the value of \f$ R(x) \f$ */ doublereal func(doublereal x); //! Set the tolerance parameters for the rootfinder /*! - * These tolerance parameters are used on the function value and the independent value - * to determine convergence + * These tolerance parameters are used on the function value and the + * independent value to determine convergence * * @param rtolf Relative tolerance. The default is 10^-5 * @param atolf absolute tolerance. The default is 10^-11 - * @param rtolx Relative tolerance. The default is 10^-5 - * Default parameter is 0.0, in which case rtolx is set equal to rtolf - * @param atolx absolute tolerance. The default is 10^-11 - * Default parameter is 0.0, in which case atolx is set equal to atolf + * @param rtolx Relative tolerance. The default is 10^-5. Default parameter + * is 0.0, in which case rtolx is set equal to rtolf + * @param atolx absolute tolerance. The default is 10^-11. Default + * parameter is 0.0, in which case atolx is set equal to + * atolf */ void setTol(doublereal rtolf, doublereal atolf, doublereal rtolx = 0.0, doublereal atolx = 0.0); //! Set the print level from the rootfinder /*! - * 0 -> absolutely nothing is printed for a single time step. - * 1 -> One line summary per solve_nonlinear call - * 2 -> short description, points of interest: Table of nonlinear solve - one line per iteration - * 3 -> Table is included -> More printing per nonlinear iteration (default) that occurs during the table - * 4 -> Summaries of the nonlinear solve iteration as they are occurring -> table no longer printed - * 5 -> Algorithm information on the nonlinear iterates are printed out - * 6 -> Additional info on the nonlinear iterates are printed out - * 7 -> Additional info on the linear solve is printed out. - * 8 -> Info on a per iterate of the linear solve is printed out. + * - 0: No printing of any kind + * - 1: Single print line indicating success or failure of the routine. + * - 2: Summary table printed at the end of the routine, with a convergence + * history + * - 3: Printouts during the iteration are added. Summary table is printed + * out at the end. if writeLogAllowed_ is turned on, a file is written + * out with the convergence history. * * @param printLvl integer value */ @@ -260,12 +260,12 @@ public: //! Set the function behavior flag /*! - * If this is true, the function is generally an increasing function of x. - * In particular, if the algorithm is seeking a higher value of f, it will look - * in the positive x direction. + * If this is true, the function is generally an increasing function of x. + * In particular, if the algorithm is seeking a higher value of f, it will + * look in the positive x direction. * - * This type of function is needed because this algorithm must deal with regions of f(x) where - * f is not changing with x. + * This type of function is needed because this algorithm must deal with + * regions of f(x) where f is not changing with x. * * @param value boolean value */ @@ -273,12 +273,12 @@ public: //! Set the function behavior flag /*! - * If this is true, the function is generally a decreasing function of x. - * In particular, if the algorithm is seeking a higher value of f, it will look - * in the negative x direction. + * If this is true, the function is generally a decreasing function of x. In + * particular, if the algorithm is seeking a higher value of f, it will look + * in the negative x direction. * - * This type of function is needed because this algorithm must deal with regions of f(x) where - * f is not changing with x. + * This type of function is needed because this algorithm must deal with + * regions of f(x) where f is not changing with x. * * @param value boolean value */ @@ -286,16 +286,12 @@ public: //! Set the minimum value of deltaX /*! - * This sets the value of deltaXNorm_ - * - * @param deltaXNorm + * @param deltaXNorm */ void setDeltaX(doublereal deltaXNorm); //! Set the maximum value of deltaX /*! - * This sets the value of deltaXMax_ - * * @param deltaX */ void setDeltaXMax(doublereal deltaX); @@ -304,36 +300,30 @@ public: void printTable(); public: - //! Pointer to the residual function evaluator + //! Pointer to the residual function evaluator ResidEval* m_residFunc; - //! Target value for the function. We seek the value of f that is equal to this value + //! Target value for the function. We seek the value of f that is equal to + //! this value doublereal m_funcTargetValue; - //! Absolute tolerance for the value of f + //! Absolute tolerance for the value of f doublereal m_atolf; - //! Absolute tolerance for the value of x + //! Absolute tolerance for the value of x doublereal m_atolx; - //! Relative tolerance for the value of f and x + //! Relative tolerance for the value of f and x doublereal m_rtolf; - //! Relative tolerance for the value of x + //! Relative tolerance for the value of x doublereal m_rtolx; - //! Maximum number of step sizes + //! Maximum number of step sizes doublereal m_maxstep; protected: - //! Print level - /*! - * 0 No printing of any kind - * 1 Single print line indicating success or failure of the routine. - * 2 Summary table printed at the end of the routine, with a convergence history - * 3 Printouts during the iteration are added. Summary table is printed out at the end. - * if writeLogAllowed_ is turned on, a file is written out with the convergence history. - */ + //! Print level. @see setPrintLvl int printLvl; public: @@ -341,31 +331,35 @@ public: bool writeLogAllowed_; protected: - //! Delta X norm. This is the nominal value of deltaX that will be used by the program + //! Delta X norm. This is the nominal value of deltaX that will be used by + //! the program doublereal DeltaXnorm_; - //! Boolean indicating whether DeltaXnorm_ has been specified by the user or not + //! Boolean indicating whether DeltaXnorm_ has been specified by the user or + //! not int specifiedDeltaXnorm_; - //! Delta X Max. This is the maximum value of deltaX that will be used by the program + //! Delta X Max. /*! - * Sometimes a large change in x causes problems. + * This is the maximum value of deltaX that will be used by the program. + * Sometimes a large change in x causes problems. */ doublereal DeltaXMax_; - //! Boolean indicating whether DeltaXMax_ has been specified by the user or not + //! Boolean indicating whether DeltaXMax_ has been specified by the user or + //! not int specifiedDeltaXMax_; //! Boolean indicating whether the function is an increasing with x bool FuncIsGenerallyIncreasing_; - //! Boolean indicating whether the function is decreasing with x + //! Boolean indicating whether the function is decreasing with x bool FuncIsGenerallyDecreasing_; //! Value of delta X that is needed for convergence /*! - * X will be considered as converged if we are within deltaXConverged_ of the solution - * The default is zero. + * X will be considered as converged if we are within deltaXConverged_ of + * the solution The default is zero. */ doublereal deltaXConverged_; diff --git a/include/cantera/numerics/SquareMatrix.h b/include/cantera/numerics/SquareMatrix.h index 1e9b1d42f..18b0e4c90 100644 --- a/include/cantera/numerics/SquareMatrix.h +++ b/include/cantera/numerics/SquareMatrix.h @@ -1,7 +1,4 @@ -/** - * @file SquareMatrix.h - * Dense, Square (not sparse) matrices. - */ +//! @file SquareMatrix.h Dense, Square (not sparse) matrices. /* * Copyright 2004 Sandia Corporation. Under the terms of Contract @@ -20,8 +17,8 @@ namespace Cantera { /** - * A class for full (non-sparse) matrices with Fortran-compatible - * data storage. Adds matrix inversion operations to this class from DenseMatrix. + * 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, public GeneralMatrix { @@ -38,10 +35,7 @@ public: */ SquareMatrix(size_t n, doublereal v = 0.0); - //! Copy Constructor SquareMatrix(const SquareMatrix& right); - - //! Assignment operator SquareMatrix& operator=(const SquareMatrix& right); int solve(doublereal* b, size_t nrhs=1, size_t ldb=0); @@ -63,7 +57,8 @@ public: virtual doublereal oneNorm() const; - //! Solves the linear problem Ax=b using the QR algorithm returning x in the b spot + //! Solves the linear problem Ax=b using the QR algorithm returning x in the + //! b spot /*! * @param b RHS to be solved. */ @@ -99,10 +94,10 @@ public: /*! * This is inherited from GeneralMatrix * - * @param iStruct OUTPUT Pointer to a vector of ints that describe the structure of the matrix. - * not used + * @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. + * @returns the number of rows and columns in the matrix. */ size_t nRowsAndStruct(size_t* const iStruct = 0) const; @@ -125,10 +120,11 @@ public: //! Integer work vector for QR algorithms vector_int iwork_; protected: - //! 1-norm of the matrix. This is determined immediately before every factorization + //! 1-norm of the matrix. This is determined immediately before every + //! factorization doublereal a1norm_; - //! Use the QR algorithm to factor and invert the matrix + //! Use the QR algorithm to factor and invert the matrix int useQR_; }; } diff --git a/include/cantera/numerics/funcs.h b/include/cantera/numerics/funcs.h index 9b4d69ede..460956315 100644 --- a/include/cantera/numerics/funcs.h +++ b/include/cantera/numerics/funcs.h @@ -16,17 +16,15 @@ namespace Cantera //! Linearly interpolate a function defined on a discrete grid. /*! - * Vector xpts contains a monotonic sequence of grid points, and - * vector fpts contains function values defined at these points. - * The value returned is the linear interpolate at point x. - * If x is outside the range of xpts, the value of fpts at the - * nearest end is returned. + * Vector xpts contains a monotonic sequence of grid points, and vector fpts + * contains function values defined at these points. The value returned is the + * linear interpolate at point x. If x is outside the range of xpts, the value + * of fpts at the nearest end is returned. * * @param x value of the x coordinate * @param xpts value of the grid points * @param fpts value of the interpolant at the grid points - * @return Returned value is the value of of the interpolated - * function at x. + * @returns the value of of the interpolated function at x. */ doublereal linearInterp(doublereal x, const vector_fp& xpts, const vector_fp& fpts); diff --git a/include/cantera/numerics/polyfit.h b/include/cantera/numerics/polyfit.h index 6d56682dc..c7cf472b2 100644 --- a/include/cantera/numerics/polyfit.h +++ b/include/cantera/numerics/polyfit.h @@ -14,62 +14,54 @@ namespace Cantera { -//! Fits a polynomial function to a set of data points +//! Fits a polynomial function to a set of data points /*! - * Given a collection of points X(I) and a set of values Y(I) which - * correspond to some function or measurement at each of the X(I), - * subroutine DPOLFT computes the weighted least-squares polynomial - * fits of all degrees up to some degree either specified by the user - * or determined by the routine. The fits thus obtained are in - * orthogonal polynomial form. Subroutine DP1VLU may then be - * called to evaluate the fitted polynomials and any of their - * derivatives at any point. The subroutine DPCOEF may be used to - * express the polynomial fits as powers of (X-C) for any specified - * point C. + * Given a collection of points X(I) and a set of values Y(I) which correspond + * to some function or measurement at each of the X(I), subroutine DPOLFT + * computes the weighted least-squares polynomial fits of all degrees up to some + * degree either specified by the user or determined by the routine. The fits + * thus obtained are in orthogonal polynomial form. Subroutine DP1VLU may then + * be called to evaluate the fitted polynomials and any of their derivatives at + * any point. The subroutine DPCOEF may be used to express the polynomial fits + * as powers of (X-C) for any specified point C. * - * @param n The number of data points. - * @param x A set of grid points on which the data is specified. - * The array of values of the independent variable. These - * values may appear in any order and need not all be - * distinct. There are n of them. - * @param y array of corresponding function values. There are n of them - * @param w array of positive values to be used as weights. If - * W[0] is negative, DPOLFT will set all the weights - * to 1.0, which means unweighted least squares error - * will be minimized. To minimize relative error, the - * user should set the weights to: W(I) = 1.0/Y(I)**2, - * I = 1,...,N . - * @param maxdeg maximum degree to be allowed for polynomial fit. - * MAXDEG may be any non-negative integer less than N. - * Note -- MAXDEG cannot be equal to N-1 when a - * statistical test is to be used for degree selection, - * i.e., when input value of EPS is negative. - * @param ndeg output degree of the fit computed. - * @param eps Specifies the criterion to be used in determining - * the degree of fit to be computed. - * (1) If EPS is input negative, DPOLFT chooses the - * degree based on a statistical F test of - * significance. One of three possible - * significance levels will be used: .01, .05 or - * .10. If EPS=-1.0 , the routine will - * automatically select one of these levels based - * on the number of data points and the maximum - * degree to be considered. If EPS is input as - * -.01, -.05, or -.10, a significance level of - * .01, .05, or .10, respectively, will be used. - * (2) If EPS is set to 0., DPOLFT computes the - * polynomials of degrees 0 through MAXDEG . - * (3) If EPS is input positive, EPS is the RMS - * error tolerance which must be satisfied by the - * fitted polynomial. DPOLFT will increase the - * degree of fit until this criterion is met or - * until the maximum degree is reached. - * @param r Output vector containing the first ndeg+1 Taylor coefficients + * @param n The number of data points. + * @param x A set of grid points on which the data is specified. The array of + * values of the independent variable. These values may appear in + * any order and need not all be distinct. There are n of them. + * @param y array of corresponding function values. There are n of them + * @param w array of positive values to be used as weights. If W[0] is + * negative, DPOLFT will set all the weights to 1.0, which means + * unweighted least squares error will be minimized. To minimize + * relative error, the user should set the weights to: W(I) = + * 1.0/Y(I)**2, I = 1,...,N . + * @param maxdeg maximum degree to be allowed for polynomial fit. MAXDEG may be + * any non-negative integer less than N. Note -- MAXDEG cannot be + * equal to N-1 when a statistical test is to be used for degree + * selection, i.e., when input value of EPS is negative. + * @param ndeg output degree of the fit computed. + * @param eps Specifies the criterion to be used in determining the degree of + * fit to be computed. + * 1. If EPS is input negative, DPOLFT chooses the degree based on a + * statistical F test of significance. One of three possible + * significance levels will be used: .01, .05 or .10. If + * EPS=-1.0 , the routine will automatically select one of these + * levels based on the number of data points and the maximum + * degree to be considered. If EPS is input as -.01, -.05, or + * -.10, a significance level of .01, .05, or .10, respectively, + * will be used. + * 2. If EPS is set to 0., DPOLFT computes the polynomials of degrees + * 0 through MAXDEG. + * 3. If EPS is input positive, EPS is the RMS error tolerance which + * must be satisfied by the fitted polynomial. DPOLFT will + * increase the degree of fit until this criterion is met or until + * the maximum degree is reached. * - * P(X) = r[0] + r[1]*(X-C) + ... + r[ndeg] * (X-C)**ndeg + * @param r Output vector containing the first ndeg+1 Taylor coefficients + * + * P(X) = r[0] + r[1]*(X-C) + ... + r[ndeg] * (X-C)**ndeg * ( here C = 0.0) - * @return Returned value is the value of the rms of the interpolated - * function at x. + * @returns value of the rms of the interpolated function at x. */ doublereal polyfit(int n, doublereal* x, doublereal* y, doublereal* w, int maxdeg, int& ndeg, doublereal eps, doublereal* r); diff --git a/src/numerics/BandMatrix.cpp b/src/numerics/BandMatrix.cpp index f3bba7ec7..134c60442 100644 --- a/src/numerics/BandMatrix.cpp +++ b/src/numerics/BandMatrix.cpp @@ -1,8 +1,4 @@ -/** - * @file BandMatrix.cpp - * - * Banded matrices. - */ +//! @file BandMatrix.cpp Banded matrices. // Copyright 2001 California Institute of Technology diff --git a/src/numerics/CVodesIntegrator.cpp b/src/numerics/CVodesIntegrator.cpp index 8ae2526fa..056dd47d1 100644 --- a/src/numerics/CVodesIntegrator.cpp +++ b/src/numerics/CVodesIntegrator.cpp @@ -1,6 +1,4 @@ -/** - * @file CVodesIntegrator.cpp - */ +//! @file CVodesIntegrator.cpp // Copyright 2001 California Institute of Technology #include "cantera/numerics/CVodesIntegrator.h" @@ -51,13 +49,12 @@ public: extern "C" { /** - * Function called by cvodes to evaluate ydot given y. The CVODE - * integrator allows passing in a void* pointer to access - * external data. This pointer is cast to a pointer to a instance - * of class FuncEval. The equations to be integrated should be - * specified by deriving a class from FuncEval that evaluates the - * desired equations. - * @ingroup odeGroup + * Function called by cvodes to evaluate ydot given y. The CVODE integrator + * allows passing in a void* pointer to access external data. This pointer + * is cast to a pointer to a instance of class FuncEval. The equations to be + * integrated should be specified by deriving a class from FuncEval that + * evaluates the desired equations. + * @ingroup odeGroup */ static int cvodes_rhs(realtype t, N_Vector y, N_Vector ydot, void* f_data) @@ -287,12 +284,9 @@ void CVodesIntegrator::initialize(double t0, FuncEval& func) CVodeFree(&m_cvode_mem); } - /* - * Specify the method and the iteration type: - * Cantera Defaults: - * CV_BDF - Use BDF methods - * CV_NEWTON - use Newton's method - */ + //! Specify the method and the iteration type. Cantera Defaults: + //! CV_BDF - Use BDF methods + //! CV_NEWTON - use Newton's method m_cvode_mem = CVodeCreate(m_method, m_iter); if (!m_cvode_mem) { throw CVodesErr("CVodeCreate failed."); diff --git a/src/numerics/DAE_solvers.cpp b/src/numerics/DAE_solvers.cpp index fa2981d78..6ca6edf5b 100644 --- a/src/numerics/DAE_solvers.cpp +++ b/src/numerics/DAE_solvers.cpp @@ -1,6 +1,4 @@ -/** - * @file DAE_solvers.cpp Factory routine for picking the DAE solver package - */ +//! @file DAE_solvers.cpp Factory routine for picking the DAE solver package /* * Copyright 2004 Sandia Corporation. Under the terms of Contract diff --git a/src/numerics/DenseMatrix.cpp b/src/numerics/DenseMatrix.cpp index 5fb5e868e..430aa72b7 100644 --- a/src/numerics/DenseMatrix.cpp +++ b/src/numerics/DenseMatrix.cpp @@ -100,7 +100,8 @@ void DenseMatrix::mult(const DenseMatrix& B, DenseMatrix& prod) const const doublereal* const* bcols = B.const_colPts(); doublereal* const* prodcols = prod.colPts(); for (size_t col=0; col < m_ncols; ++col) { - // Loop over ncols multiplying A*column of B and storing in corresponding prod column + // Loop over ncols multiplying A*column of B and storing in + // corresponding prod column mult(bcols[col], prodcols[col]); } } diff --git a/src/numerics/GeneralMatrix.cpp b/src/numerics/GeneralMatrix.cpp index 370ef9295..937b4d3aa 100644 --- a/src/numerics/GeneralMatrix.cpp +++ b/src/numerics/GeneralMatrix.cpp @@ -1,6 +1,4 @@ -/** - * @file GeneralMatrix.cpp - */ +//! @file GeneralMatrix.cpp /* * Copyright 2004 Sandia Corporation. Under the terms of Contract diff --git a/src/numerics/IDA_Solver.cpp b/src/numerics/IDA_Solver.cpp index 3b732defa..2af60b189 100644 --- a/src/numerics/IDA_Solver.cpp +++ b/src/numerics/IDA_Solver.cpp @@ -1,6 +1,4 @@ -/** - * @file IDA_Solver.cpp - */ +//! @file IDA_Solver.cpp // Copyright 2006 California Institute of Technology @@ -27,10 +25,8 @@ typedef long int sd_size_t; namespace Cantera { -/** - * A simple class to hold an array of parameter values and a pointer to - * an instance of a subclass of ResidEval. - */ +//! A simple class to hold an array of parameter values and a pointer to an +//! instance of a subclass of ResidEval. class ResidData { public: @@ -48,20 +44,23 @@ public: } extern "C" { - //! Function called by IDA to evaluate the residual, given y and ydot. + //! Function called by IDA to evaluate the residual, given y and ydot. /*! - * IDA allows passing in a void* pointer to access external data. Instead of requiring the user to provide a - * residual function directly to IDA (which would require using - * the sundials data types N_Vector, etc.), we define this function as the single function that IDA always calls. The - * real evaluation of the residual is done by an instance of a subclass of ResidEval, passed in to this - * function as a pointer in the parameters. + * IDA allows passing in a void* pointer to access external data. Instead of + * requiring the user to provide a residual function directly to IDA (which + * would require using the sundials data types N_Vector, etc.), we define + * this function as the single function that IDA always calls. The real + * evaluation of the residual is done by an instance of a subclass of + * ResidEval, passed in to this function as a pointer in the parameters. * - * FROM IDA WRITEUP -> What the IDA solver expects as a return flag from its residual routines ------ - * A IDAResFn res should return a value of 0 if successful, a positive - * value if a recoverable error occured (e.g. yy has an illegal value), - * or a negative value if a nonrecoverable error occured. In the latter - * case, the program halts. If a recoverable error occured, the integrator - * will attempt to correct and retry. + * FROM IDA WRITEUP -> What the IDA solver expects as a return flag from its + * residual routines: + * + * A IDAResFn res should return a value of 0 if successful, a positive value + * if a recoverable error occured (e.g. yy has an illegal value), or a + * negative value if a nonrecoverable error occured. In the latter case, the + * program halts. If a recoverable error occured, the integrator will + * attempt to correct and retry. */ static int ida_resid(realtype t, N_Vector y, N_Vector ydot, N_Vector r, void* f_data) { @@ -90,11 +89,12 @@ extern "C" { * N_Vector tmp1, N_Vector tmp2, N_Vector tmp3); * * A IDADlsDenseJacFn should return - * 0 if successful, - * a positive int if a recoverable error occurred, or - * a negative int if a nonrecoverable error occurred. - * In the case of a recoverable error return, the integrator will - * attempt to recover by reducing the stepsize (which changes cj). + * - 0 if successful, + * - a positive int if a recoverable error occurred, or + * - a negative int if a nonrecoverable error occurred. + * + * In the case of a recoverable error return, the integrator will attempt to + * recover by reducing the stepsize (which changes cj). */ static int ida_jacobian(sd_size_t nrows, realtype t, realtype c_j, N_Vector y, N_Vector ydot, N_Vector r, DlsMat Jac, void* f_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) @@ -371,10 +371,7 @@ void IDA_Solver::init(doublereal t0) } } - //----------------------------------- // set the linear solver type - //----------------------------------- - if (m_type == 1 || m_type == 0) { long int N = m_neq; flag = IDADense(m_ida_mem, N); diff --git a/src/numerics/ResidJacEval.cpp b/src/numerics/ResidJacEval.cpp index c933b3157..5a5ce2ea9 100644 --- a/src/numerics/ResidJacEval.cpp +++ b/src/numerics/ResidJacEval.cpp @@ -1,6 +1,4 @@ -/** - * @file ResidJacEval.cpp - */ +//! @file ResidJacEval.cpp /* * Copyright 2004 Sandia Corporation. Under the terms of Contract diff --git a/src/numerics/RootFind.cpp b/src/numerics/RootFind.cpp index 237cd799b..1e5fe2fcf 100644 --- a/src/numerics/RootFind.cpp +++ b/src/numerics/RootFind.cpp @@ -1,6 +1,4 @@ -/** - * @file: RootFind.cpp root finder for 1D problems - */ +//! @file: RootFind.cpp root finder for 1D problems /* * Copyright 2004 Sandia Corporation. Under the terms of Contract @@ -25,7 +23,7 @@ using namespace std; namespace Cantera { -//! Print out a form for the current function evaluation +//! Print out a form for the current function evaluation /*! * @param fp Pointer to the FILE object * @param xval Current value of x @@ -170,11 +168,8 @@ bool RootFind::theSame(doublereal x2, doublereal x1, doublereal factor) const int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& funcTargetValue, doublereal* xbest) { - /* - * We store the function target and then actually calculate a modified functional - * - * func = eval(x1) - m_funcTargetValue = 0 - */ + // We store the function target and then actually calculate a modified + // functional, func = eval(x1) - m_funcTargetValue = 0 m_funcTargetValue = funcTargetValue; static int callNum = 0; @@ -196,7 +191,7 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun doublereal fPosF = 1.0E300; doublereal xNegF = 0.0; doublereal fNegF = -1.0E300; - doublereal fnorm; /* A valid norm for the making the function value dimensionless */ + doublereal fnorm; // A valid norm for the making the function value dimensionless doublereal xDelMin; doublereal sgn; doublereal fnoise = 0.0; @@ -220,10 +215,8 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun return ROOTFIND_BADINPUT; } - /* - * If the maximum step size has not been specified, set it here to 1/5 of the - * domain range of x. - */ + // If the maximum step size has not been specified, set it here to 1/5 of + // the domain range of x. if (!specifiedDeltaXMax_) { DeltaXMax_ = 0.2 *(xmax - xmin); } @@ -240,9 +233,7 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun } } - /* - * Calculate an initial value of deltaXConverged_ - */ + // Calculate an initial value of deltaXConverged_ deltaXConverged_ = m_rtolx * (*xbest) + m_atolx; if (DeltaXnorm_ < deltaXConverged_) { writelogf("%s DeltaXnorm_, %g, is too small compared to tols, increasing to %g\n", @@ -250,10 +241,8 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun DeltaXnorm_ = deltaXConverged_; } - /* - * Find the first function value f1 = func(x1), by using the value entered into xbest. - * Process it - */ + // Find the first function value f1 = func(x1), by using the value entered + // into xbest. Process it x1 = *xbest; if (x1 < xmin || x1 > xmax) { x1 = (xmin + xmax) / 2.0; @@ -296,10 +285,9 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun rfHistory_.push_back(rfT); rfT.clear(); - /* - * Now, this is actually a tricky part of the algorithm - Find the x value for - * the second point. It's tricky because we don't have a valid idea of the scale of x yet - */ + // Now, this is actually a tricky part of the algorithm - Find the x value + // for the second point. It's tricky because we don't have a valid idea of + // the scale of x yet rfT.reasoning = "Second Point: "; if (x1 == 0.0) { x2 = x1 + 0.01 * DeltaXnorm_; @@ -313,9 +301,7 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun rfT.reasoning += " - But adjusted to be within bounds"; } - /* - * Find the second function value f2 = func(x2), Process it - */ + // Find the second function value f2 = func(x2), Process it deltaX2 = x2 - x1; its++; f2 = func(x2); @@ -324,10 +310,9 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun fprintf(fp, "%-5d %-5d %-15.5E %-15.5E", -1, 0, x2, f2); } - /* - * Calculate the norm of the function, this is the nominal value of f. We try - * to reduce the nominal value of f by rtolf, this is the main convergence requirement. - */ + // Calculate the norm of the function, this is the nominal value of f. We + // try to reduce the nominal value of f by rtolf, this is the main + // convergence requirement. if (m_funcTargetValue != 0.0) { fnorm = m_atolf + fabs(m_funcTargetValue); } else { @@ -357,9 +342,7 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun rfT.foundPos = foundPosF; rfT.foundNeg = foundNegF; - /* - * See if we have already achieved a straddle - */ + // See if we have already achieved a straddle foundStraddle = foundPosF && foundNegF; if (foundStraddle) { if (xPosF > xNegF) { @@ -371,14 +354,11 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun bool useNextStrat = false; bool slopePointingToHigher = true; - // --------------------------------------------------------------------------------------------- - // MAIN LOOP - // --------------------------------------------------------------------------------------------- + + // MAIN LOOP while (!converged && its < itmax) { - /* - * Find an estimate of the next point, xnew, to try based on - * a linear approximation from the last two points. - */ + // Find an estimate of the next point, xnew, to try based on a linear + // approximation from the last two points. if (DEBUG_MODE_ENABLED && fabs(x2 - x1) < 1.0E-14) { writelogf(" RootFind: we are here x2 = %g x1 = %g\n", x2, x1); } @@ -411,9 +391,8 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun fprintf(fp, " | xlin = %-11.5E", xnew); } deltaXnew = xnew - x2; - /* - * If the suggested step size is too big, throw out step - */ + + // If the suggested step size is too big, throw out step if (!foundStraddle) { if (fabs(xnew - x2) > DeltaXMax_) { useNextStrat = true; @@ -424,9 +403,9 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun xnew = x2 + deltaXnew; } } - /* - * If the slope can't be trusted using a different strategy for picking the next point - */ + + // If the slope can't be trusted using a different strategy for picking + // the next point if (useNextStrat) { rfT.reasoning += "Using DeltaXnorm, " + fp2str(DeltaXnorm_) + " and FuncIsGenerallyIncreasing hints. "; if (f2 < 0.0) { @@ -472,10 +451,9 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun } } - /* - * Here, if we have a straddle, we purposefully overshoot the smaller side by 5%. Yes it does lead to - * more iterations. However, we're interested in bounding x, and not just doing Newton's method. - */ + // Here, if we have a straddle, we purposefully overshoot the smaller + // side by 5%. Yes it does lead to more iterations. However, we're + // interested in bounding x, and not just doing Newton's method. if (foundStraddle) { double delta = fabs(x2 - x1); if (fabs(xnew - x1) < .01 * delta) { @@ -490,16 +468,13 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun } } } - /* - * OK, we have an estimate xnew. - * Put heuristic bounds on the step jump - */ + + // OK, we have an estimate xnew. Put heuristic bounds on the step jump if ((xnew > x1 && xnew < x2) || (xnew < x1 && xnew > x2)) { - /* - * If we are doing a jump in between the two previous points, make sure - * the new trial is no closer that 10% of the distances between x2-x1 to - * any of the original points. This is an important part of finding a good bound. - */ + // If we are doing a jump in between the two previous points, make + // sure the new trial is no closer that 10% of the distances between + // x2-x1 to any of the original points. This is an important part of + // finding a good bound. xDelMin = fabs(x2 - x1) / 10.; if (fabs(xnew - x1) < xDelMin) { xnew = x1 + sign(xnew-x1) * xDelMin; @@ -514,11 +489,9 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun } } } else { - /* - * If we are venturing into new ground, only allow the step jump - * to increase by 50% at each iteration, unless the step jump is less than - * the user has said that it is ok to take - */ + // If we are venturing into new ground, only allow the step jump to + // increase by 50% at each iteration, unless the step jump is less + // than the user has said that it is ok to take doublereal xDelMax = 1.5 * fabs(x2 - x1); if (specifiedDeltaXnorm_ && 0.5 * DeltaXnorm_ > xDelMax) { xDelMax = 0.5 *DeltaXnorm_; @@ -529,11 +502,11 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun fprintf(fp, " | xlimitsize = %-11.5E", xnew); } } - /* - * If we are doing a jump outside the two previous points, make sure - * the new trial is no closer that 10% of the distances between x2-x1 to - * any of the original points. This is an important part of finding a good bound. - */ + + // If we are doing a jump outside the two previous points, make sure + // the new trial is no closer that 10% of the distances between + // x2-x1 to any of the original points. This is an important part of + // finding a good bound. xDelMin = 0.1 * fabs(x2 - x1); if (fabs(xnew - x2) < xDelMin) { xnew = x2 + sign(xnew - x2) * xDelMin; @@ -548,9 +521,8 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun } } } - /* - * HKM -> Not sure this section is needed - */ + + // HKM -> Not sure this section is needed if (foundStraddle) { double xorig = xnew; if (posStraddle) { @@ -591,9 +563,7 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun } } - /* - * Enforce a minimum stepsize if we haven't found a straddle. - */ + // Enforce a minimum stepsize if we haven't found a straddle. deltaXnew = xnew - x2; if (fabs(deltaXnew) < 1.2 * delXMeaningful(xnew) && !foundStraddle) { sgn = 1.0; @@ -606,9 +576,7 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun xnew = x2 + deltaXnew; } - /* - * Guard against going above xmax or below xmin - */ + // Guard against going above xmax or below xmin if (xnew > xmax) { topBump++; if (topBump < 3) { @@ -732,11 +700,9 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun x2 = xnew; f2 = fnew; - /* - * As we go on to new data points, we make sure that - * we have the best straddle of the solution with the choice of F1 and F2 when - * we do have a straddle to work with. - */ + // As we go on to new data points, we make sure that we have the best + // straddle of the solution with the choice of F1 and F2 when we do have + // a straddle to work with. if (foundStraddle) { bool foundBetterPos = false; bool foundBetterNeg = false; @@ -831,9 +797,7 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun rfT.delX = std::max(rfT.delX, xmax - x2); } } - /* - * Section To Determine CONVERGENCE criteria - */ + // Section To Determine CONVERGENCE criteria doFinalFuncCall = 0; if ((fabs(fnew / fnorm) < m_rtolf) && foundStraddle) { if (fabs(deltaX2) < deltaXConverged_ && fabs(deltaXnew) < deltaXConverged_) { @@ -850,9 +814,7 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun } } - /* - * Check for excess convergence in the x coordinate - */ + // Check for excess convergence in the x coordinate if (!converged && foundStraddle) { doublereal denom = fabs(x1 - x2); if (denom < 1.0E-200) { @@ -867,9 +829,7 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun } } } else { - /* - * We are here when F is not converged, but we may want to end anyway - */ + // We are here when F is not converged, but we may want to end anyway if (!converged && foundStraddle) { doublereal denom = fabs(x1 - x2); if (denom < 1.0E-200) { @@ -877,10 +837,9 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun converged = true; rfT.reasoning += "FNotConverged but X1X2Identical"; } - /* - * The premise here is that if x1 and x2 get close to one another, - * then the accuracy of the calculation gets destroyed. - */ + + // The premise here is that if x1 and x2 get close to one + // another, then the accuracy of the calculation gets destroyed. if (theSame(x2, x1, 1.0E-5)) { converged = true; retn = ROOTFIND_SUCCESS_XCONVERGENCEONLY; diff --git a/src/numerics/SquareMatrix.cpp b/src/numerics/SquareMatrix.cpp index 36cc329ff..58c00e537 100644 --- a/src/numerics/SquareMatrix.cpp +++ b/src/numerics/SquareMatrix.cpp @@ -1,6 +1,4 @@ -/** - * @file SquareMatrix.cpp - */ +//! @file SquareMatrix.cpp /* * Copyright 2004 Sandia Corporation. Under the terms of Contract @@ -59,9 +57,8 @@ int SquareMatrix::solve(doublereal* b, size_t nrhs, size_t ldb) return solveQR(b); } int info=0; - /* - * Check to see whether the matrix has been factored. - */ + + // Check to see whether the matrix has been factored. if (!m_factored) { int retn = factor(); if (retn) { @@ -71,9 +68,8 @@ int SquareMatrix::solve(doublereal* b, size_t nrhs, size_t ldb) if (ldb == 0) { ldb = nColumns(); } - /* - * Solve the factored system - */ + + // Solve the factored system ct_dgetrs(ctlapack::NoTranspose, static_cast(nRows()), nrhs, &*begin(), static_cast(nRows()), ipiv().data(), b, ldb, info); @@ -168,9 +164,8 @@ int SquareMatrix::factorQR() int SquareMatrix::solveQR(doublereal* b) { int info=0; - /* - * Check to see whether the matrix has been factored. - */ + + // Check to see whether the matrix has been factored. if (!m_factored) { int retn = factorQR(); if (retn) { @@ -184,9 +179,7 @@ int SquareMatrix::solveQR(doublereal* b) lwork = 8 * m_nrows; } - /* - * Solve the factored system - */ + // Solve the factored system ct_dormqr(ctlapack::Left, ctlapack::Transpose, m_nrows, 1, m_nrows, &*begin(), m_nrows, tau.data(), b, m_nrows, work.data(), lwork, info); if (info != 0) { diff --git a/src/numerics/funcs.cpp b/src/numerics/funcs.cpp index a228a84dd..6a9cb2a7e 100644 --- a/src/numerics/funcs.cpp +++ b/src/numerics/funcs.cpp @@ -1,7 +1,5 @@ -/** - * @file funcs.cpp file containing miscellaneous - * numerical functions. - */ +//! @file funcs.cpp file containing miscellaneous numerical functions. + /* * Copyright 2001-2003 California Institute of Technology * See file License.txt for licensing information