From 44c53a6f2712ca5a0b9c15eaa5b28f3f296ce30b Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Mon, 24 Mar 2014 04:02:18 +0000 Subject: [PATCH] [Kinetics] Abstract use of LAPACK in solveSP and solveProb --- include/cantera/kinetics/solveSP.h | 7 ++----- include/cantera/numerics/solveProb.h | 10 ++-------- src/kinetics/solveSP.cpp | 12 +++--------- src/numerics/solveProb.cpp | 10 ++-------- 4 files changed, 9 insertions(+), 30 deletions(-) diff --git a/include/cantera/kinetics/solveSP.h b/include/cantera/kinetics/solveSP.h index b7704c7b1..e850b8707 100644 --- a/include/cantera/kinetics/solveSP.h +++ b/include/cantera/kinetics/solveSP.h @@ -14,7 +14,7 @@ #include "ImplicitSurfChem.h" #include "cantera/kinetics/InterfaceKinetics.h" -#include "cantera/base/Array.h" +#include "cantera/numerics/SquareMatrix.h" //! @defgroup solvesp_methods Surface Problem Solver Methods //! @{ @@ -533,9 +533,6 @@ private: //! Vector of mole fractions. length m_maxTotSpecies vector_fp m_XMolKinSpecies; - //! pivots. length MAX(1, m_neq) - vector_int m_ipiv; - //! Vector of pointers to the top of the columns of the Jacobian /*! * The "dim" by "dim" computed Jacobian matrix for the @@ -545,7 +542,7 @@ private: //! Jacobian. m_neq by m_neq computed Jacobian matrix for the local //! Newton's method. - Array2D m_Jac; + SquareMatrix m_Jac; public: int m_ioflag; diff --git a/include/cantera/numerics/solveProb.h b/include/cantera/numerics/solveProb.h index cf3cd6607..0fef8f9dd 100644 --- a/include/cantera/numerics/solveProb.h +++ b/include/cantera/numerics/solveProb.h @@ -17,7 +17,7 @@ * @defgroup solverGroup Solvers for Equation Systems */ -#include "cantera/base/Array.h" +#include "cantera/numerics/SquareMatrix.h" #include "ResidEval.h" //! Solution Methods @@ -415,12 +415,6 @@ private: */ vector_fp m_resid; - //! pivots - /*! - * length MAX(1, m_neq) - */ - vector_int m_ipiv; - //! Vector of pointers to the top of the columns of the jacobians /*! * The "dim" by "dim" computed Jacobian matrix for the @@ -432,7 +426,7 @@ private: /*! * m_neq by m_neq computed Jacobian matrix for the local Newton's method. */ - Array2D m_Jac; + SquareMatrix m_Jac; //! Top bounds for the solution vector /*! diff --git a/src/kinetics/solveSP.cpp b/src/kinetics/solveSP.cpp index 435474536..8afa8eeae 100644 --- a/src/kinetics/solveSP.cpp +++ b/src/kinetics/solveSP.cpp @@ -10,7 +10,6 @@ #include "cantera/kinetics/solveSP.h" #include "cantera/base/clockWC.h" -#include "cantera/numerics/ctlapack.h" #include @@ -125,8 +124,6 @@ solveSP::solveSP(ImplicitSurfChem* surfChemPtr, int bulkFunc) : m_wtResid.resize(dim1, 0.0); m_wtSpecies.resize(dim1, 0.0); m_resid.resize(dim1, 0.0); - m_ipiv.resize(dim1, 0); - m_Jac.resize(dim1, dim1, 0.0); m_JacCol.resize(dim1, 0); for (size_t k = 0; k < dim1; k++) { @@ -152,7 +149,6 @@ int solveSP::solveSurfProb(int ifunc, doublereal time_scale, doublereal TKelvin, doublereal label_factor = 1.0; int iter=0; // iteration number on numlinear solver int iter_max=1000; // maximum number of nonlinear iterations - int nrhs=1; doublereal deltaT = 1.0E-10; // Delta time step doublereal damp=1.0, tmp; // Weighted L2 norm of the residual. Currently, this is only @@ -291,13 +287,11 @@ int solveSP::solveSurfProb(int ifunc, doublereal time_scale, doublereal TKelvin, DATA_PTR(m_resid), m_neq); /* - * Solve Linear system (with LAPACK). The solution is in resid[] + * Solve Linear system. The solution is in resid[] */ - ct_dgetrf(m_neq, m_neq, m_JacCol[0], m_neq, DATA_PTR(m_ipiv), info); + info = m_Jac.factor(); if (info==0) { - ct_dgetrs(ctlapack::NoTranspose, m_neq, nrhs, m_JacCol[0], - m_neq, DATA_PTR(m_ipiv), DATA_PTR(m_resid), m_neq, - info); + m_Jac.solve(&m_resid[0]); } /* * Force convergence if residual is small to avoid diff --git a/src/numerics/solveProb.cpp b/src/numerics/solveProb.cpp index 7e5238761..bebb1cf1e 100644 --- a/src/numerics/solveProb.cpp +++ b/src/numerics/solveProb.cpp @@ -11,7 +11,6 @@ #include "cantera/numerics/solveProb.h" #include "cantera/base/clockWC.h" -#include "cantera/numerics/ctlapack.h" #include "cantera/base/stringUtils.h" using namespace std; @@ -48,7 +47,6 @@ solveProb::solveProb(ResidEval* resid) : m_wtResid.resize(dim1, 0.0); m_wtSpecies.resize(dim1, 0.0); m_resid.resize(dim1, 0.0); - m_ipiv.resize(dim1, 0); m_topBounds.resize(dim1, 1.0); m_botBounds.resize(dim1, 0.0); @@ -84,7 +82,6 @@ int solveProb::solve(int ifunc, doublereal time_scale, doublereal label_factor = 1.0; int iter=0; // iteration number on numlinear solver int iter_max=1000; // maximum number of nonlinear iterations - int nrhs=1; doublereal deltaT = 1.0E-10; // Delta time step doublereal damp=1.0, tmp; // Weighted L2 norm of the residual. Currently, this is only @@ -238,12 +235,9 @@ int solveProb::solve(int ifunc, doublereal time_scale, /* * Solve Linear system (with LAPACK). The solution is in resid[] */ - - ct_dgetrf(m_neq, m_neq, m_JacCol[0], m_neq, DATA_PTR(m_ipiv), info); + info = m_Jac.factor(); if (info==0) { - ct_dgetrs(ctlapack::NoTranspose, m_neq, nrhs, m_JacCol[0], - m_neq, DATA_PTR(m_ipiv), DATA_PTR(m_resid), m_neq, - info); + m_Jac.solve(&m_resid[0]); } /* * Force convergence if residual is small to avoid