From 881fcbaf6ffd68643265fde00b67a65fd362d27e Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Sun, 27 Mar 2016 22:52:49 -0400 Subject: [PATCH] Improve use of local variables in src/numerics --- src/numerics/BandMatrix.cpp | 27 ++++++----------------- src/numerics/CVodesIntegrator.cpp | 29 ++++++------------------- src/numerics/DenseMatrix.cpp | 11 ++++------ src/numerics/Func1.cpp | 1 - src/numerics/IDA_Solver.cpp | 36 +++++++++++-------------------- 5 files changed, 29 insertions(+), 75 deletions(-) diff --git a/src/numerics/BandMatrix.cpp b/src/numerics/BandMatrix.cpp index e91d68352..99379a53d 100644 --- a/src/numerics/BandMatrix.cpp +++ b/src/numerics/BandMatrix.cpp @@ -285,10 +285,8 @@ vector_fp::const_iterator BandMatrix::end() const ostream& operator<<(ostream& s, const BandMatrix& m) { - size_t nr = m.nRows(); - size_t nc = m.nColumns(); - for (size_t i = 0; i < nr; i++) { - for (size_t j = 0; j < nc; j++) { + for (size_t i = 0; i < m.nRows(); i++) { + for (size_t j = 0; j < m.nColumns(); j++) { s << m(i,j) << ", "; } s << endl; @@ -298,30 +296,19 @@ ostream& operator<<(ostream& s, const BandMatrix& m) doublereal BandMatrix::rcond(doublereal a1norm) { - int printLevel = 0; - int useReturnErrorCode = 0; - if (iwork_.size() < m_n) { - iwork_.resize(m_n); - } - if (work_.size() < 3 * m_n) { - work_.resize(3 * m_n); - } - doublereal rcond = 0.0; + iwork_.resize(m_n); + work_.resize(3 * m_n); + if (m_factored != 1) { throw CanteraError("BandMatrix::rcond()", "matrix isn't factored correctly"); } size_t ldab = (2 *m_kl + m_ku + 1); int rinfo = 0; - rcond = ct_dgbcon('1', m_n, m_kl, m_ku, ludata.data(), ldab, m_ipiv.data(), a1norm, work_.data(), + double rcond = ct_dgbcon('1', m_n, m_kl, m_ku, ludata.data(), ldab, m_ipiv.data(), a1norm, work_.data(), iwork_.data(), rinfo); if (rinfo != 0) { - if (printLevel) { - writelogf("BandMatrix::rcond(): DGBCON returned INFO = %d\n", rinfo); - } - if (! useReturnErrorCode) { - throw CanteraError("BandMatrix::rcond()", "DGBCON returned INFO = {}", rinfo); - } + throw CanteraError("BandMatrix::rcond()", "DGBCON returned INFO = {}", rinfo); } return rcond; } diff --git a/src/numerics/CVodesIntegrator.cpp b/src/numerics/CVodesIntegrator.cpp index 95429305d..671f3f648 100644 --- a/src/numerics/CVodesIntegrator.cpp +++ b/src/numerics/CVodesIntegrator.cpp @@ -58,15 +58,9 @@ extern "C" { void* f_data) { try { - double* ydata = NV_DATA_S(y); - double* ydotdata = NV_DATA_S(ydot); FuncData* d = (FuncData*)f_data; FuncEval* f = d->m_func; - if (d->m_pars.size() == 0) { - f->eval(t, ydata, ydotdata, NULL); - } else { - f->eval(t, ydata, ydotdata, d->m_pars.data()); - } + f->eval(t, NV_DATA_S(y), NV_DATA_S(ydot), d->m_pars.data()); } catch (CanteraError& err) { std::cerr << err.what() << std::endl; return 1; // possibly recoverable error @@ -232,18 +226,12 @@ void CVodesIntegrator::setIterator(IterType t) void CVodesIntegrator::sensInit(double t0, FuncEval& func) { m_np = func.nparams(); - size_t nv = func.neq(); m_sens_ok = false; - doublereal* data; - N_Vector y; - y = N_VNew_Serial(static_cast(nv)); + N_Vector y = N_VNew_Serial(static_cast(func.neq())); m_yS = N_VCloneVectorArray_Serial(static_cast(m_np), y); for (size_t n = 0; n < m_np; n++) { - data = NV_DATA_S(m_yS[n]); - for (size_t j = 0; j < nv; j++) { - data[j] =0.0; - } + N_VConst(0.0, m_yS[n]); } int flag = CVodeSensInit(m_cvode_mem, static_cast(m_np), @@ -253,9 +241,7 @@ void CVodesIntegrator::sensInit(double t0, FuncEval& func) throw CanteraError("CVodesIntegrator::sensInit", "Error in CVodeSensMalloc"); } vector_fp atol(m_np, m_abstolsens); - double rtol = m_reltolsens; - flag = CVodeSensSStolerances(m_cvode_mem, rtol, atol.data()); - + flag = CVodeSensSStolerances(m_cvode_mem, m_reltolsens, atol.data()); } void CVodesIntegrator::initialize(double t0, FuncEval& func) @@ -268,9 +254,7 @@ void CVodesIntegrator::initialize(double t0, FuncEval& func) N_VDestroy_Serial(m_y); // free solution vector if already allocated } m_y = N_VNew_Serial(static_cast(m_neq)); // allocate solution vector - for (size_t i = 0; i < m_neq; i++) { - NV_Ith_S(m_y, i) = 0.0; - } + N_VConst(0.0, m_y); // check abs tolerance array size if (m_itol == CV_SV && m_nabs < m_neq) { throw CanteraError("CVodesIntegrator::initialize", @@ -347,8 +331,7 @@ void CVodesIntegrator::reinitialize(double t0, FuncEval& func) m_time = t0; func.getState(NV_DATA_S(m_y)); - int result; - result = CVodeReInit(m_cvode_mem, m_t0, m_y); + int result = CVodeReInit(m_cvode_mem, m_t0, m_y); if (result != CV_SUCCESS) { throw CanteraError("CVodesIntegrator::reinitialize", "CVodeReInit failed. result = {}", result); diff --git a/src/numerics/DenseMatrix.cpp b/src/numerics/DenseMatrix.cpp index db86b41bf..6a200e930 100644 --- a/src/numerics/DenseMatrix.cpp +++ b/src/numerics/DenseMatrix.cpp @@ -108,12 +108,9 @@ void DenseMatrix::mult(const DenseMatrix& B, DenseMatrix& prod) const void DenseMatrix::leftMult(const double* const b, double* const prod) const { - size_t nc = nColumns(); - size_t nr = nRows(); - double sum = 0.0; - for (size_t n = 0; n < nc; n++) { - sum = 0.0; - for (size_t i = 0; i < nr; i++) { + for (size_t n = 0; n < nColumns(); n++) { + double sum = 0.0; + for (size_t i = 0; i < nRows(); i++) { sum += value(i,n)*b[i]; } prod[n] = sum; @@ -127,13 +124,13 @@ vector_int& DenseMatrix::ipiv() int solve(DenseMatrix& A, double* b, size_t nrhs, size_t ldb) { - int info = 0; if (A.nColumns() != A.nRows()) { if (A.m_printLevel) { writelogf("solve(DenseMatrix& A, double* b): Can only solve a square matrix\n"); } throw CanteraError("solve(DenseMatrix& A, double* b)", "Can only solve a square matrix"); } + int info = 0; ct_dgetrf(A.nRows(), A.nColumns(), A.ptrColumn(0), A.nRows(), &A.ipiv()[0], info); if (info > 0) { diff --git a/src/numerics/Func1.cpp b/src/numerics/Func1.cpp index d18d53263..0d5ef34f2 100644 --- a/src/numerics/Func1.cpp +++ b/src/numerics/Func1.cpp @@ -216,7 +216,6 @@ string Func1::write(const std::string& arg) const string Pow1::write(const std::string& arg) const { - string c = ""; if (m_c == 0.5) { return "\\sqrt{" + arg + "}"; } diff --git a/src/numerics/IDA_Solver.cpp b/src/numerics/IDA_Solver.cpp index d26c6735d..cd23d3216 100644 --- a/src/numerics/IDA_Solver.cpp +++ b/src/numerics/IDA_Solver.cpp @@ -63,21 +63,19 @@ extern "C" { */ static int ida_resid(realtype t, N_Vector y, N_Vector ydot, N_Vector r, void* f_data) { - double* ydata = NV_DATA_S(y); - double* ydotdata = NV_DATA_S(ydot); - double* rdata = NV_DATA_S(r); Cantera::ResidData* d = (Cantera::ResidData*) f_data; Cantera::ResidJacEval* f = d->m_func; Cantera::IDA_Solver* s = d->m_solver; double delta_t = s->getCurrentStepFromIDA(); // TODO evaluate evalType. Assumed to be Base_ResidEval - int retn = 0; - int flag = f->evalResidNJ(t, delta_t, ydata, ydotdata, rdata); + int flag = f->evalResidNJ(t, delta_t, NV_DATA_S(y), NV_DATA_S(ydot), + NV_DATA_S(r)); if (flag < 0) { // This signals to IDA that a nonrecoverable error has occurred. - retn = flag; + return flag; + } else { + return 0; } - return retn; } //! Function called by by IDA to evaluate the Jacobian, given y and ydot. @@ -98,15 +96,12 @@ extern "C" { 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) { - doublereal* ydata = NV_DATA_S(y); - doublereal* ydotdata = NV_DATA_S(ydot); - doublereal* rdata = NV_DATA_S(r); Cantera::ResidData* d = (Cantera::ResidData*) f_data; Cantera::ResidJacEval* f = d->m_func; - doublereal* const* colPts = Jac->cols; Cantera::IDA_Solver* s = d->m_solver; double delta_t = s->getCurrentStepFromIDA(); - f->evalJacobianDP(t, delta_t, c_j, ydata, ydotdata, colPts, rdata); + f->evalJacobianDP(t, delta_t, c_j, NV_DATA_S(y), NV_DATA_S(ydot), + Jac->cols, NV_DATA_S(r)); return 0; } } @@ -465,7 +460,6 @@ void IDA_Solver::init(doublereal t0) void IDA_Solver::correctInitial_Y_given_Yp(doublereal* y, doublereal* yp, doublereal tout) { - int icopt = IDA_Y_INIT; doublereal tout1 = tout; if (tout == 0.0) { double h0 = 1.0E-5; @@ -475,7 +469,7 @@ void IDA_Solver::correctInitial_Y_given_Yp(doublereal* y, doublereal* yp, double tout1 = m_t0 + h0; } - int flag = IDACalcIC(m_ida_mem, icopt, tout1); + int flag = IDACalcIC(m_ida_mem, IDA_Y_INIT, tout1); if (flag != IDA_SUCCESS) { throw CanteraError("IDA_Solver::correctInitial_Y_given_Yp", "IDACalcIC failed: error = {}", flag); @@ -486,12 +480,9 @@ void IDA_Solver::correctInitial_Y_given_Yp(doublereal* y, doublereal* yp, double throw CanteraError("IDA_Solver::correctInitial_Y_given_Yp", "IDAGetSolution failed: error = {}", flag); } - doublereal* yy = NV_DATA_S(m_y); - doublereal* yyp = NV_DATA_S(m_ydot); - for (int i = 0; i < m_neq; i++) { - y[i] = yy[i]; - yp[i] = yyp[i]; + y[i] = NV_Ith_S(m_y, i); + yp[i] = NV_Ith_S(m_ydot, i); } } @@ -518,12 +509,9 @@ void IDA_Solver::correctInitial_YaYp_given_Yd(doublereal* y, doublereal* yp, dou throw CanteraError("IDA_Solver::correctInitial_YaYp_given_Yd", "IDAGetSolution failed: error = {}", flag); } - doublereal* yy = NV_DATA_S(m_y); - doublereal* yyp = NV_DATA_S(m_ydot); - for (int i = 0; i < m_neq; i++) { - y[i] = yy[i]; - yp[i] = yyp[i]; + y[i] = NV_Ith_S(m_y, i); + yp[i] = NV_Ith_S(m_ydot, i); } }