research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Improve use of local variables in src/numerics

Browse source
This commit is contained in:
Ray Speth 2016-03-27 22:52:49 -04:00
parent 0995725a13
commit 881fcbaf6f
5 changed files with 29 additions and 75 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

27
src/numerics/BandMatrix.cpp
View file

@ -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;
}

29
src/numerics/CVodesIntegrator.cpp
View file

@ -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<sd_size_t>(nv));
N_Vector y = N_VNew_Serial(static_cast<sd_size_t>(func.neq()));
m_yS = N_VCloneVectorArray_Serial(static_cast<sd_size_t>(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<sd_size_t>(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<sd_size_t>(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);

11
src/numerics/DenseMatrix.cpp
View file

@ -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) {

1
src/numerics/Func1.cpp
View file

@ -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 + "}";
}

36
src/numerics/IDA_Solver.cpp
View file

@ -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);
}
}