Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Review complexity of sympy generated (scary?) code for KINETIC block #930

Closed
pramodk opened this issue Sep 9, 2022 · 5 comments
Closed

Review complexity of sympy generated (scary?) code for KINETIC block #930

pramodk opened this issue Sep 9, 2022 · 5 comments
Assignees
@cattabiani
Labels
codegen Code generation backend performance Related to performance improvement solver Solver and numerical methods

Comments

@pramodk
Copy link
Contributor

pramodk commented Sep 9, 2022

While fixing the issues in #888, I came across glia__dbbs_mod_collection__cdp5__CR.mod file (use attached version!). By using following nmodl command:

./bin/nmodl glia__dbbs_mod_collection__cdp5__CR.mod host --c passes --inline sympy --analytic

The generated code from sympy for KINETIC block looks like:

                functor(NrnThread* nt, glia__dbbs_mod_collection__cdp5__CR_Instance* inst, int id, int pnodecount, double v, Datum* indexes, double* data, ThreadDatum* thread) : nt{nt}, inst{inst}, id{id}, pnodecount{pnodecount}, v{v}, indexes{indexes}, data{data}, thread{thread} {}
                void operator()(const Eigen::Matrix<double, 21, 1>& nmodl_eigen_xm, Eigen::Matrix<double, 21, 1>& nmodl_eigen_fm, Eigen::Matrix<double, 21, 21>& nmodl_eigen_jm) const {
                    const double* nmodl_eigen_x = nmodl_eigen_xm.data();
                    double* nmodl_eigen_j = nmodl_eigen_jm.data();
                    double* nmodl_eigen_f = nmodl_eigen_fm.data();
                    nmodl_eigen_f[static_cast<int>(0)] = 0.5 * (2.0 * FARADAY * pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id] * ( -nmodl_eigen_x[static_cast<int>(0)] + old_ca) + FARADAY * nt->_dt * ( -2.0 * nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(11)] * inst->global->c1 * inst->dsqvol[id] - 2.0 * nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(13)] * inst->dsqvol[id] * inst->global->nT1 - 2.0 * nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(14)] * inst->dsqvol[id] * inst->global->nR1 - 4.0 * nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(15)] * inst->dsqvol[id] * inst->global->nR1 - 2.0 * nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(16)] * inst->dsqvol[id] * inst->global->nR1 - 2.0 * nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(16)] * inst->dsqvol[id] * inst->global->nT1 - 2.0 * nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(17)] * inst->dsqvol[id] * inst->global->nT1 - 2.0 * nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(18)] * inst->dsqvol[id] * inst->global->nR1 - 20000000000.0 * nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(19)] * inst->global->kpmp1 * inst->parea[id] - 4.0 * nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(1)] * inst->dsqvol[id] * inst->global->nT1 - 2.0 * nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(1)] * inst->dsqvol[id] * inst->global->nV1 - 2.0 * nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(2)] * inst->dsqvol[id] * inst->global->nR1 - 2.0 * nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(2)] * inst->dsqvol[id] * inst->global->nT1 - 2.0 * nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(5)] * inst->dsqvol[id] * inst->global->rf1 - 2.0 * nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(7)] * inst->dsqvol[id] * inst->rf3[id] - 2.0 * nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(9)] * inst->global->b1 * inst->dsqvol[id] + 2.0 * nmodl_eigen_x[static_cast<int>(10)] * inst->global->b2 * inst->dsqvol[id] + 2.0 * nmodl_eigen_x[static_cast<int>(12)] * inst->global->c2 * inst->dsqvol[id] + 2.0 * nmodl_eigen_x[static_cast<int>(13)] * inst->dsqvol[id] * inst->global->nR2 + 2.0 * nmodl_eigen_x[static_cast<int>(14)] * inst->dsqvol[id] * inst->global->nR2 + 2.0 * nmodl_eigen_x[static_cast<int>(14)] * inst->dsqvol[id] * inst->global->nT2 + 4.0 * nmodl_eigen_x[static_cast<int>(15)] * inst->dsqvol[id] * inst->global->nT2 + 2.0 * nmodl_eigen_x[static_cast<int>(16)] * inst->dsqvol[id] * inst->global->nT2 + 2.0 * nmodl_eigen_x[static_cast<int>(17)] * inst->dsqvol[id] * inst->global->nR2 + 2.0 * nmodl_eigen_x[static_cast<int>(18)] * inst->dsqvol[id] * inst->global->nR2 + 2.0 * nmodl_eigen_x[static_cast<int>(18)] * inst->dsqvol[id] * inst->global->nT2 + 20000000000.0 * nmodl_eigen_x[static_cast<int>(20)] * inst->global->kpmp2 * inst->parea[id] + 2.0 * nmodl_eigen_x[static_cast<int>(2)] * inst->dsqvol[id] * inst->global->nT2 + 4.0 * nmodl_eigen_x[static_cast<int>(3)] * inst->dsqvol[id] * inst->global->nR2 + 2.0 * nmodl_eigen_x[static_cast<int>(4)] * inst->dsqvol[id] * inst->global->nV2 + 2.0 * nmodl_eigen_x[static_cast<int>(6)] * inst->dsqvol[id] * inst->global->rf2 + 2.0 * nmodl_eigen_x[static_cast<int>(8)] * inst->dsqvol[id] * inst->rf4[id]) - PI * inst->diam[indexes[4*pnodecount + id]] * nt->_dt * inst->ica[id]) / (FARADAY * pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(0)] =  -(nmodl_eigen_x[static_cast<int>(11)] * inst->global->c1 * inst->dsqvol[id] * nt->_dt + nmodl_eigen_x[static_cast<int>(13)] * inst->dsqvol[id] * nt->_dt * inst->global->nT1 + nmodl_eigen_x[static_cast<int>(14)] * inst->dsqvol[id] * nt->_dt * inst->global->nR1 + 2.0 * nmodl_eigen_x[static_cast<int>(15)] * inst->dsqvol[id] * nt->_dt * inst->global->nR1 + nmodl_eigen_x[static_cast<int>(16)] * inst->dsqvol[id] * nt->_dt * inst->global->nR1 + nmodl_eigen_x[static_cast<int>(16)] * inst->dsqvol[id] * nt->_dt * inst->global->nT1 + nmodl_eigen_x[static_cast<int>(17)] * inst->dsqvol[id] * nt->_dt * inst->global->nT1 + nmodl_eigen_x[static_cast<int>(18)] * inst->dsqvol[id] * nt->_dt * inst->global->nR1 + 10000000000.0 * nmodl_eigen_x[static_cast<int>(19)] * nt->_dt * inst->global->kpmp1 * inst->parea[id] + 2.0 * nmodl_eigen_x[static_cast<int>(1)] * inst->dsqvol[id] * nt->_dt * inst->global->nT1 + nmodl_eigen_x[static_cast<int>(1)] * inst->dsqvol[id] * nt->_dt * inst->global->nV1 + nmodl_eigen_x[static_cast<int>(2)] * inst->dsqvol[id] * nt->_dt * inst->global->nR1 + nmodl_eigen_x[static_cast<int>(2)] * inst->dsqvol[id] * nt->_dt * inst->global->nT1 + nmodl_eigen_x[static_cast<int>(5)] * inst->dsqvol[id] * nt->_dt * inst->global->rf1 + nmodl_eigen_x[static_cast<int>(7)] * inst->dsqvol[id] * nt->_dt * inst->rf3[id] + nmodl_eigen_x[static_cast<int>(9)] * inst->global->b1 * inst->dsqvol[id] * nt->_dt + pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(21)] =  -nmodl_eigen_x[static_cast<int>(0)] * inst->dsqvol[id] * nt->_dt * (2.0 * inst->global->nT1 + inst->global->nV1) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(42)] = inst->dsqvol[id] * nt->_dt * ( -nmodl_eigen_x[static_cast<int>(0)] * inst->global->nR1 - nmodl_eigen_x[static_cast<int>(0)] * inst->global->nT1 + inst->global->nT2) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(63)] = 2.0 * inst->dsqvol[id] * nt->_dt * inst->global->nR2 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(84)] = inst->dsqvol[id] * nt->_dt * inst->global->nV2 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(105)] =  -nmodl_eigen_x[static_cast<int>(0)] * inst->dsqvol[id] * nt->_dt * inst->global->rf1 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(126)] = inst->dsqvol[id] * nt->_dt * inst->global->rf2 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(147)] =  -nmodl_eigen_x[static_cast<int>(0)] * inst->dsqvol[id] * nt->_dt * inst->rf3[id] / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(168)] = inst->dsqvol[id] * nt->_dt * inst->rf4[id] / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(189)] =  -nmodl_eigen_x[static_cast<int>(0)] * inst->global->b1 * inst->dsqvol[id] * nt->_dt / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(210)] = inst->global->b2 * inst->dsqvol[id] * nt->_dt / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(231)] =  -nmodl_eigen_x[static_cast<int>(0)] * inst->global->c1 * inst->dsqvol[id] * nt->_dt / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(252)] = inst->global->c2 * inst->dsqvol[id] * nt->_dt / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(273)] = inst->dsqvol[id] * nt->_dt * ( -nmodl_eigen_x[static_cast<int>(0)] * inst->global->nT1 + inst->global->nR2) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(294)] = inst->dsqvol[id] * nt->_dt * ( -nmodl_eigen_x[static_cast<int>(0)] * inst->global->nR1 + inst->global->nR2 + inst->global->nT2) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(315)] = 2.0 * inst->dsqvol[id] * nt->_dt * ( -nmodl_eigen_x[static_cast<int>(0)] * inst->global->nR1 + inst->global->nT2) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(336)] = inst->dsqvol[id] * nt->_dt * ( -nmodl_eigen_x[static_cast<int>(0)] * inst->global->nR1 - nmodl_eigen_x[static_cast<int>(0)] * inst->global->nT1 + inst->global->nT2) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(357)] = inst->dsqvol[id] * nt->_dt * ( -nmodl_eigen_x[static_cast<int>(0)] * inst->global->nT1 + inst->global->nR2) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(378)] = inst->dsqvol[id] * nt->_dt * ( -nmodl_eigen_x[static_cast<int>(0)] * inst->global->nR1 + inst->global->nR2 + inst->global->nT2) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(399)] =  -10000000000.0 * nmodl_eigen_x[static_cast<int>(0)] * nt->_dt * inst->global->kpmp1 * inst->parea[id] / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(420)] = 10000000000.0 * nt->_dt * inst->global->kpmp2 * inst->parea[id] / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_f[static_cast<int>(1)] = ( -2.0 * nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(1)] * inst->dsqvol[id] * nt->_dt * inst->global->nT1 - nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(1)] * inst->dsqvol[id] * nt->_dt * inst->global->nV1 + nmodl_eigen_x[static_cast<int>(16)] * inst->dsqvol[id] * nt->_dt * inst->global->nT2 + nmodl_eigen_x[static_cast<int>(2)] * inst->dsqvol[id] * nt->_dt * inst->global->nT2 + nmodl_eigen_x[static_cast<int>(4)] * inst->dsqvol[id] * nt->_dt * inst->global->nV2 + pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id] * ( -nmodl_eigen_x[static_cast<int>(1)] + old_CR)) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(1)] =  -nmodl_eigen_x[static_cast<int>(1)] * inst->dsqvol[id] * nt->_dt * (2.0 * inst->global->nT1 + inst->global->nV1) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(22)] =  -(2.0 * nmodl_eigen_x[static_cast<int>(0)] * inst->dsqvol[id] * nt->_dt * inst->global->nT1 + nmodl_eigen_x[static_cast<int>(0)] * inst->dsqvol[id] * nt->_dt * inst->global->nV1 + pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(43)] = inst->dsqvol[id] * nt->_dt * inst->global->nT2 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(64)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(85)] = inst->dsqvol[id] * nt->_dt * inst->global->nV2 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(106)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(127)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(148)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(169)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(190)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(211)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(232)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(253)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(274)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(295)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(316)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(337)] = inst->dsqvol[id] * nt->_dt * inst->global->nT2 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(358)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(379)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(400)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(421)] = 0.0;
                    nmodl_eigen_f[static_cast<int>(2)] = (nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(1)] * inst->dsqvol[id] * nt->_dt * inst->global->nT1 - nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(2)] * inst->dsqvol[id] * nt->_dt * inst->global->nR1 - nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(2)] * inst->dsqvol[id] * nt->_dt * inst->global->nT1 + nmodl_eigen_x[static_cast<int>(13)] * inst->dsqvol[id] * nt->_dt * inst->global->nR2 + nmodl_eigen_x[static_cast<int>(15)] * inst->dsqvol[id] * nt->_dt * inst->global->nT2 - nmodl_eigen_x[static_cast<int>(2)] * inst->dsqvol[id] * nt->_dt * inst->global->nT2 + pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id] * ( -nmodl_eigen_x[static_cast<int>(2)] + old_CR_1C_0N)) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(2)] = inst->dsqvol[id] * nt->_dt * (nmodl_eigen_x[static_cast<int>(1)] * inst->global->nT1 - nmodl_eigen_x[static_cast<int>(2)] * inst->global->nR1 - nmodl_eigen_x[static_cast<int>(2)] * inst->global->nT1) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(23)] = nmodl_eigen_x[static_cast<int>(0)] * inst->dsqvol[id] * nt->_dt * inst->global->nT1 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(44)] =  -(nmodl_eigen_x[static_cast<int>(0)] * inst->dsqvol[id] * nt->_dt * inst->global->nR1 + nmodl_eigen_x[static_cast<int>(0)] * inst->dsqvol[id] * nt->_dt * inst->global->nT1 + pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id] + inst->dsqvol[id] * nt->_dt * inst->global->nT2) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(65)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(86)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(107)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(128)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(149)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(170)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(191)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(212)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(233)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(254)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(275)] = inst->dsqvol[id] * nt->_dt * inst->global->nR2 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(296)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(317)] = inst->dsqvol[id] * nt->_dt * inst->global->nT2 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(338)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(359)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(380)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(401)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(422)] = 0.0;
                    nmodl_eigen_f[static_cast<int>(3)] = (nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(14)] * inst->dsqvol[id] * nt->_dt * inst->global->nR1 + nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(18)] * inst->dsqvol[id] * nt->_dt * inst->global->nR1 - 2.0 * nmodl_eigen_x[static_cast<int>(3)] * inst->dsqvol[id] * nt->_dt * inst->global->nR2 + pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id] * ( -nmodl_eigen_x[static_cast<int>(3)] + old_CR_2C_2N)) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(3)] = inst->dsqvol[id] * nt->_dt * inst->global->nR1 * (nmodl_eigen_x[static_cast<int>(14)] + nmodl_eigen_x[static_cast<int>(18)]) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(24)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(45)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(66)] =  -1.0 - 2.0 * inst->dsqvol[id] * nt->_dt * inst->global->nR2 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(87)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(108)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(129)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(150)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(171)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(192)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(213)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(234)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(255)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(276)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(297)] = nmodl_eigen_x[static_cast<int>(0)] * inst->dsqvol[id] * nt->_dt * inst->global->nR1 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(318)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(339)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(360)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(381)] = nmodl_eigen_x[static_cast<int>(0)] * inst->dsqvol[id] * nt->_dt * inst->global->nR1 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(402)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(423)] = 0.0;
                    nmodl_eigen_f[static_cast<int>(4)] = nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(1)] * inst->dsqvol[id] * nt->_dt * inst->global->nV1 - nmodl_eigen_x[static_cast<int>(4)] * inst->dsqvol[id] * nt->_dt * inst->global->nV2 - nmodl_eigen_x[static_cast<int>(4)] + old_CR_1V;
                    nmodl_eigen_j[static_cast<int>(4)] = nmodl_eigen_x[static_cast<int>(1)] * inst->dsqvol[id] * nt->_dt * inst->global->nV1;
                    nmodl_eigen_j[static_cast<int>(25)] = nmodl_eigen_x[static_cast<int>(0)] * inst->dsqvol[id] * nt->_dt * inst->global->nV1;
                    nmodl_eigen_j[static_cast<int>(46)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(67)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(88)] =  -inst->dsqvol[id] * nt->_dt * inst->global->nV2 - 1.0;
                    nmodl_eigen_j[static_cast<int>(109)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(130)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(151)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(172)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(193)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(214)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(235)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(256)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(277)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(298)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(319)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(340)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(361)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(382)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(403)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(424)] = 0.0;
                    nmodl_eigen_f[static_cast<int>(5)] = ( -nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(5)] * inst->dsqvol[id] * nt->_dt * inst->global->rf1 + nmodl_eigen_x[static_cast<int>(6)] * inst->dsqvol[id] * nt->_dt * inst->global->rf2 + pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id] * ( -nmodl_eigen_x[static_cast<int>(5)] + old_Buff1)) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(5)] =  -nmodl_eigen_x[static_cast<int>(5)] * inst->dsqvol[id] * nt->_dt * inst->global->rf1 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(26)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(47)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(68)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(89)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(110)] =  -nmodl_eigen_x[static_cast<int>(0)] * inst->dsqvol[id] * nt->_dt * inst->global->rf1 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]) - 1.0;
                    nmodl_eigen_j[static_cast<int>(131)] = inst->dsqvol[id] * nt->_dt * inst->global->rf2 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(152)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(173)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(194)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(215)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(236)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(257)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(278)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(299)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(320)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(341)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(362)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(383)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(404)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(425)] = 0.0;
                    nmodl_eigen_f[static_cast<int>(6)] = (nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(5)] * inst->dsqvol[id] * nt->_dt * inst->global->rf1 - nmodl_eigen_x[static_cast<int>(6)] * inst->dsqvol[id] * nt->_dt * inst->global->rf2 + pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id] * ( -nmodl_eigen_x[static_cast<int>(6)] + old_Buff1_ca)) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(6)] = nmodl_eigen_x[static_cast<int>(5)] * inst->dsqvol[id] * nt->_dt * inst->global->rf1 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(27)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(48)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(69)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(90)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(111)] = nmodl_eigen_x[static_cast<int>(0)] * inst->dsqvol[id] * nt->_dt * inst->global->rf1 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(132)] =  -1.0 - inst->dsqvol[id] * nt->_dt * inst->global->rf2 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(153)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(174)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(195)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(216)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(237)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(258)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(279)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(300)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(321)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(342)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(363)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(384)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(405)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(426)] = 0.0;
                    nmodl_eigen_f[static_cast<int>(7)] = ( -nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(7)] * inst->dsqvol[id] * nt->_dt * inst->rf3[id] + nmodl_eigen_x[static_cast<int>(8)] * inst->dsqvol[id] * nt->_dt * inst->rf4[id] + pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id] * ( -nmodl_eigen_x[static_cast<int>(7)] + old_Buff2)) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(7)] =  -nmodl_eigen_x[static_cast<int>(7)] * inst->dsqvol[id] * nt->_dt * inst->rf3[id] / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(28)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(49)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(70)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(91)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(112)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(133)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(154)] =  -nmodl_eigen_x[static_cast<int>(0)] * inst->dsqvol[id] * nt->_dt * inst->rf3[id] / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]) - 1.0;
                    nmodl_eigen_j[static_cast<int>(175)] = inst->dsqvol[id] * nt->_dt * inst->rf4[id] / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(196)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(217)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(238)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(259)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(280)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(301)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(322)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(343)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(364)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(385)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(406)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(427)] = 0.0;
                    nmodl_eigen_f[static_cast<int>(8)] = (nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(7)] * inst->dsqvol[id] * nt->_dt * inst->rf3[id] - nmodl_eigen_x[static_cast<int>(8)] * inst->dsqvol[id] * nt->_dt * inst->rf4[id] + pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id] * ( -nmodl_eigen_x[static_cast<int>(8)] + old_Buff2_ca)) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(8)] = nmodl_eigen_x[static_cast<int>(7)] * inst->dsqvol[id] * nt->_dt * inst->rf3[id] / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(29)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(50)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(71)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(92)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(113)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(134)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(155)] = nmodl_eigen_x[static_cast<int>(0)] * inst->dsqvol[id] * nt->_dt * inst->rf3[id] / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(176)] =  -1.0 - inst->dsqvol[id] * nt->_dt * inst->rf4[id] / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(197)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(218)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(239)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(260)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(281)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(302)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(323)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(344)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(365)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(386)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(407)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(428)] = 0.0;
                    nmodl_eigen_f[static_cast<int>(9)] = ( -nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(9)] * inst->global->b1 * inst->dsqvol[id] * nt->_dt + nmodl_eigen_x[static_cast<int>(10)] * inst->global->b2 * inst->dsqvol[id] * nt->_dt + pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id] * ( -nmodl_eigen_x[static_cast<int>(9)] + old_BTC)) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(9)] =  -nmodl_eigen_x[static_cast<int>(9)] * inst->global->b1 * inst->dsqvol[id] * nt->_dt / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(30)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(51)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(72)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(93)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(114)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(135)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(156)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(177)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(198)] =  -nmodl_eigen_x[static_cast<int>(0)] * inst->global->b1 * inst->dsqvol[id] * nt->_dt / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]) - 1.0;
                    nmodl_eigen_j[static_cast<int>(219)] = inst->global->b2 * inst->dsqvol[id] * nt->_dt / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(240)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(261)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(282)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(303)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(324)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(345)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(366)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(387)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(408)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(429)] = 0.0;
                    nmodl_eigen_f[static_cast<int>(10)] = (nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(9)] * inst->global->b1 * inst->dsqvol[id] * nt->_dt - nmodl_eigen_x[static_cast<int>(10)] * inst->global->b2 * inst->dsqvol[id] * nt->_dt + pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id] * ( -nmodl_eigen_x[static_cast<int>(10)] + old_BTC_ca)) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(10)] = nmodl_eigen_x[static_cast<int>(9)] * inst->global->b1 * inst->dsqvol[id] * nt->_dt / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(31)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(52)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(73)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(94)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(115)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(136)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(157)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(178)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(199)] = nmodl_eigen_x[static_cast<int>(0)] * inst->global->b1 * inst->dsqvol[id] * nt->_dt / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(220)] =  -inst->global->b2 * inst->dsqvol[id] * nt->_dt / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]) - 1.0;
                    nmodl_eigen_j[static_cast<int>(241)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(262)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(283)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(304)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(325)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(346)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(367)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(388)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(409)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(430)] = 0.0;
                    nmodl_eigen_f[static_cast<int>(11)] = ( -nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(11)] * inst->global->c1 * inst->dsqvol[id] * nt->_dt + nmodl_eigen_x[static_cast<int>(12)] * inst->global->c2 * inst->dsqvol[id] * nt->_dt + pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id] * ( -nmodl_eigen_x[static_cast<int>(11)] + old_DMNPE)) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(11)] =  -nmodl_eigen_x[static_cast<int>(11)] * inst->global->c1 * inst->dsqvol[id] * nt->_dt / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(32)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(53)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(74)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(95)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(116)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(137)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(158)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(179)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(200)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(221)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(242)] =  -nmodl_eigen_x[static_cast<int>(0)] * inst->global->c1 * inst->dsqvol[id] * nt->_dt / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]) - 1.0;
                    nmodl_eigen_j[static_cast<int>(263)] = inst->global->c2 * inst->dsqvol[id] * nt->_dt / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(284)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(305)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(326)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(347)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(368)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(389)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(410)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(431)] = 0.0;
                    nmodl_eigen_f[static_cast<int>(12)] = (nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(11)] * inst->global->c1 * inst->dsqvol[id] * nt->_dt - nmodl_eigen_x[static_cast<int>(12)] * inst->global->c2 * inst->dsqvol[id] * nt->_dt + pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id] * ( -nmodl_eigen_x[static_cast<int>(12)] + old_DMNPE_ca)) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(12)] = nmodl_eigen_x[static_cast<int>(11)] * inst->global->c1 * inst->dsqvol[id] * nt->_dt / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(33)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(54)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(75)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(96)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(117)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(138)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(159)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(180)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(201)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(222)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(243)] = nmodl_eigen_x[static_cast<int>(0)] * inst->global->c1 * inst->dsqvol[id] * nt->_dt / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(264)] =  -inst->global->c2 * inst->dsqvol[id] * nt->_dt / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]) - 1.0;
                    nmodl_eigen_j[static_cast<int>(285)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(306)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(327)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(348)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(369)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(390)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(411)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(432)] = 0.0;
                    nmodl_eigen_f[static_cast<int>(13)] = ( -nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(13)] * inst->dsqvol[id] * nt->_dt * inst->global->nT1 + nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(2)] * inst->dsqvol[id] * nt->_dt * inst->global->nR1 - nmodl_eigen_x[static_cast<int>(13)] * inst->dsqvol[id] * nt->_dt * inst->global->nR2 + nmodl_eigen_x[static_cast<int>(14)] * inst->dsqvol[id] * nt->_dt * inst->global->nT2 + pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id] * ( -nmodl_eigen_x[static_cast<int>(13)] + old_CR_2C_0N)) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(13)] = inst->dsqvol[id] * nt->_dt * ( -nmodl_eigen_x[static_cast<int>(13)] * inst->global->nT1 + nmodl_eigen_x[static_cast<int>(2)] * inst->global->nR1) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(34)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(55)] = nmodl_eigen_x[static_cast<int>(0)] * inst->dsqvol[id] * nt->_dt * inst->global->nR1 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(76)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(97)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(118)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(139)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(160)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(181)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(202)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(223)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(244)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(265)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(286)] =  -(nmodl_eigen_x[static_cast<int>(0)] * inst->dsqvol[id] * nt->_dt * inst->global->nT1 + pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id] + inst->dsqvol[id] * nt->_dt * inst->global->nR2) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(307)] = inst->dsqvol[id] * nt->_dt * inst->global->nT2 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(328)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(349)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(370)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(391)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(412)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(433)] = 0.0;
                    nmodl_eigen_f[static_cast<int>(14)] = (nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(13)] * inst->dsqvol[id] * nt->_dt * inst->global->nT1 - nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(14)] * inst->dsqvol[id] * nt->_dt * inst->global->nR1 + nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(15)] * inst->dsqvol[id] * nt->_dt * inst->global->nR1 - nmodl_eigen_x[static_cast<int>(14)] * inst->dsqvol[id] * nt->_dt * inst->global->nR2 - nmodl_eigen_x[static_cast<int>(14)] * inst->dsqvol[id] * nt->_dt * inst->global->nT2 + nmodl_eigen_x[static_cast<int>(3)] * inst->dsqvol[id] * nt->_dt * inst->global->nR2 + pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id] * ( -nmodl_eigen_x[static_cast<int>(14)] + old_CR_2C_1N)) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(14)] = inst->dsqvol[id] * nt->_dt * (nmodl_eigen_x[static_cast<int>(13)] * inst->global->nT1 - nmodl_eigen_x[static_cast<int>(14)] * inst->global->nR1 + nmodl_eigen_x[static_cast<int>(15)] * inst->global->nR1) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(35)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(56)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(77)] = inst->dsqvol[id] * nt->_dt * inst->global->nR2 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(98)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(119)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(140)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(161)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(182)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(203)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(224)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(245)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(266)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(287)] = nmodl_eigen_x[static_cast<int>(0)] * inst->dsqvol[id] * nt->_dt * inst->global->nT1 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(308)] =  -(nmodl_eigen_x[static_cast<int>(0)] * inst->dsqvol[id] * nt->_dt * inst->global->nR1 + pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id] + inst->dsqvol[id] * nt->_dt * inst->global->nR2 + inst->dsqvol[id] * nt->_dt * inst->global->nT2) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(329)] = nmodl_eigen_x[static_cast<int>(0)] * inst->dsqvol[id] * nt->_dt * inst->global->nR1 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(350)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(371)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(392)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(413)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(434)] = 0.0;
                    nmodl_eigen_f[static_cast<int>(15)] = ( -2.0 * nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(15)] * inst->dsqvol[id] * nt->_dt * inst->global->nR1 + nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(16)] * inst->dsqvol[id] * nt->_dt * inst->global->nT1 + nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(2)] * inst->dsqvol[id] * nt->_dt * inst->global->nT1 + nmodl_eigen_x[static_cast<int>(14)] * inst->dsqvol[id] * nt->_dt * inst->global->nR2 - 2.0 * nmodl_eigen_x[static_cast<int>(15)] * inst->dsqvol[id] * nt->_dt * inst->global->nT2 + nmodl_eigen_x[static_cast<int>(18)] * inst->dsqvol[id] * nt->_dt * inst->global->nR2 + pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id] * ( -nmodl_eigen_x[static_cast<int>(15)] + old_CR_1C_1N)) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(15)] = inst->dsqvol[id] * nt->_dt * ( -2.0 * nmodl_eigen_x[static_cast<int>(15)] * inst->global->nR1 + nmodl_eigen_x[static_cast<int>(16)] * inst->global->nT1 + nmodl_eigen_x[static_cast<int>(2)] * inst->global->nT1) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(36)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(57)] = nmodl_eigen_x[static_cast<int>(0)] * inst->dsqvol[id] * nt->_dt * inst->global->nT1 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(78)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(99)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(120)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(141)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(162)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(183)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(204)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(225)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(246)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(267)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(288)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(309)] = inst->dsqvol[id] * nt->_dt * inst->global->nR2 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(330)] =  -(2.0 * nmodl_eigen_x[static_cast<int>(0)] * inst->dsqvol[id] * nt->_dt * inst->global->nR1 + pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id] + 2.0 * inst->dsqvol[id] * nt->_dt * inst->global->nT2) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(351)] = nmodl_eigen_x[static_cast<int>(0)] * inst->dsqvol[id] * nt->_dt * inst->global->nT1 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(372)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(393)] = inst->dsqvol[id] * nt->_dt * inst->global->nR2 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(414)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(435)] = 0.0;
                    nmodl_eigen_f[static_cast<int>(16)] = ( -nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(16)] * inst->dsqvol[id] * nt->_dt * inst->global->nR1 - nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(16)] * inst->dsqvol[id] * nt->_dt * inst->global->nT1 + nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(1)] * inst->dsqvol[id] * nt->_dt * inst->global->nT1 + nmodl_eigen_x[static_cast<int>(15)] * inst->dsqvol[id] * nt->_dt * inst->global->nT2 - nmodl_eigen_x[static_cast<int>(16)] * inst->dsqvol[id] * nt->_dt * inst->global->nT2 + nmodl_eigen_x[static_cast<int>(17)] * inst->dsqvol[id] * nt->_dt * inst->global->nR2 + pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id] * ( -nmodl_eigen_x[static_cast<int>(16)] + old_CR_0C_1N)) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(16)] = inst->dsqvol[id] * nt->_dt * ( -nmodl_eigen_x[static_cast<int>(16)] * inst->global->nR1 - nmodl_eigen_x[static_cast<int>(16)] * inst->global->nT1 + nmodl_eigen_x[static_cast<int>(1)] * inst->global->nT1) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(37)] = nmodl_eigen_x[static_cast<int>(0)] * inst->dsqvol[id] * nt->_dt * inst->global->nT1 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(58)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(79)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(100)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(121)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(142)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(163)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(184)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(205)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(226)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(247)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(268)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(289)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(310)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(331)] = inst->dsqvol[id] * nt->_dt * inst->global->nT2 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(352)] =  -(nmodl_eigen_x[static_cast<int>(0)] * inst->dsqvol[id] * nt->_dt * inst->global->nR1 + nmodl_eigen_x[static_cast<int>(0)] * inst->dsqvol[id] * nt->_dt * inst->global->nT1 + pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id] + inst->dsqvol[id] * nt->_dt * inst->global->nT2) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(373)] = inst->dsqvol[id] * nt->_dt * inst->global->nR2 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(394)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(415)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(436)] = 0.0;
                    nmodl_eigen_f[static_cast<int>(17)] = (nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(16)] * inst->dsqvol[id] * nt->_dt * inst->global->nR1 - nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(17)] * inst->dsqvol[id] * nt->_dt * inst->global->nT1 - nmodl_eigen_x[static_cast<int>(17)] * inst->dsqvol[id] * nt->_dt * inst->global->nR2 + nmodl_eigen_x[static_cast<int>(18)] * inst->dsqvol[id] * nt->_dt * inst->global->nT2 + pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id] * ( -nmodl_eigen_x[static_cast<int>(17)] + old_CR_0C_2N)) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(17)] = inst->dsqvol[id] * nt->_dt * (nmodl_eigen_x[static_cast<int>(16)] * inst->global->nR1 - nmodl_eigen_x[static_cast<int>(17)] * inst->global->nT1) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(38)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(59)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(80)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(101)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(122)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(143)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(164)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(185)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(206)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(227)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(248)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(269)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(290)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(311)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(332)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(353)] = nmodl_eigen_x[static_cast<int>(0)] * inst->dsqvol[id] * nt->_dt * inst->global->nR1 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(374)] =  -(nmodl_eigen_x[static_cast<int>(0)] * inst->dsqvol[id] * nt->_dt * inst->global->nT1 + pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id] + inst->dsqvol[id] * nt->_dt * inst->global->nR2) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(395)] = inst->dsqvol[id] * nt->_dt * inst->global->nT2 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(416)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(437)] = 0.0;
                    nmodl_eigen_f[static_cast<int>(18)] = (nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(15)] * inst->dsqvol[id] * nt->_dt * inst->global->nR1 + nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(17)] * inst->dsqvol[id] * nt->_dt * inst->global->nT1 - nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(18)] * inst->dsqvol[id] * nt->_dt * inst->global->nR1 - nmodl_eigen_x[static_cast<int>(18)] * inst->dsqvol[id] * nt->_dt * inst->global->nR2 - nmodl_eigen_x[static_cast<int>(18)] * inst->dsqvol[id] * nt->_dt * inst->global->nT2 + nmodl_eigen_x[static_cast<int>(3)] * inst->dsqvol[id] * nt->_dt * inst->global->nR2 + pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id] * ( -nmodl_eigen_x[static_cast<int>(18)] + old_CR_1C_2N)) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(18)] = inst->dsqvol[id] * nt->_dt * (nmodl_eigen_x[static_cast<int>(15)] * inst->global->nR1 + nmodl_eigen_x[static_cast<int>(17)] * inst->global->nT1 - nmodl_eigen_x[static_cast<int>(18)] * inst->global->nR1) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(39)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(60)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(81)] = inst->dsqvol[id] * nt->_dt * inst->global->nR2 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(102)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(123)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(144)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(165)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(186)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(207)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(228)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(249)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(270)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(291)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(312)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(333)] = nmodl_eigen_x[static_cast<int>(0)] * inst->dsqvol[id] * nt->_dt * inst->global->nR1 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(354)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(375)] = nmodl_eigen_x[static_cast<int>(0)] * inst->dsqvol[id] * nt->_dt * inst->global->nT1 / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(396)] =  -(nmodl_eigen_x[static_cast<int>(0)] * inst->dsqvol[id] * nt->_dt * inst->global->nR1 + pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id] + inst->dsqvol[id] * nt->_dt * inst->global->nR2 + inst->dsqvol[id] * nt->_dt * inst->global->nT2) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(417)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(438)] = 0.0;
                    nmodl_eigen_f[static_cast<int>(19)] =  -1.0 * nmodl_eigen_x[static_cast<int>(0)] * nmodl_eigen_x[static_cast<int>(19)] * nt->_dt * inst->global->kpmp1 - 1.0 * nmodl_eigen_x[static_cast<int>(19)] + 1.0 * nmodl_eigen_x[static_cast<int>(20)] * nt->_dt * inst->global->kpmp2 + 1.0 * nmodl_eigen_x[static_cast<int>(20)] * nt->_dt * inst->global->kpmp3 + 1.0 * old_pump;
                    nmodl_eigen_j[static_cast<int>(19)] =  -1.0 * nmodl_eigen_x[static_cast<int>(19)] * nt->_dt * inst->global->kpmp1;
                    nmodl_eigen_j[static_cast<int>(40)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(61)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(82)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(103)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(124)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(145)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(166)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(187)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(208)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(229)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(250)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(271)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(292)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(313)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(334)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(355)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(376)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(397)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(418)] =  -1.0 * nmodl_eigen_x[static_cast<int>(0)] * nt->_dt * inst->global->kpmp1 - 1.0;
                    nmodl_eigen_j[static_cast<int>(439)] = 1.0 * nt->_dt * (inst->global->kpmp2 + inst->global->kpmp3);
                    nmodl_eigen_f[static_cast<int>(20)] = 1.0 * inst->TotalPump[id] - 1.0 * nmodl_eigen_x[static_cast<int>(19)] - 1.0 * nmodl_eigen_x[static_cast<int>(20)];
                    nmodl_eigen_j[static_cast<int>(20)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(41)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(62)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(83)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(104)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(125)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(146)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(167)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(188)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(209)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(230)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(251)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(272)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(293)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(314)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(335)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(356)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(377)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(398)] = 0.0;
                    nmodl_eigen_j[static_cast<int>(419)] =  -1.0;
                    nmodl_eigen_j[static_cast<int>(440)] =  -1.0;
                }
[glia__dbbs_mod_collection__cdp5__CR.mod.zip](https://github.com/BlueBrain/nmodl/files/9536176/glia__dbbs_mod_collection__cdp5__CR.mod.zip)


                void finalize() {
                }
            };
@pramodk pramodk added codegen Code generation backend solver Solver and numerical methods performance Related to performance improvement labels Sep 9, 2022
@pramodk pramodk mentioned this issue Sep 9, 2022
Open
6 tasks
@cattabiani
Copy link
Contributor

I do not see the problem. We have a matrix represented as a vector and we are filling it with the values. It is scary because it is ~20x20 (afaik). With 20 different variables the problem is deemed to be complex.

@ohm314
Copy link
Contributor

ohm314 commented Sep 14, 2022

I was also gonna say, this looks fine to me? One potential improvement could be to first initialize all with 0 (when allocating) and then only set the non-zero values.

@1uc
Copy link
Collaborator

1uc commented Aug 25, 2023

More easily visible link to the mod file:
https://github.com/BlueBrain/nmodl/files/9536176/glia__dbbs_mod_collection__cdp5__CR.mod.zip

@1uc
Copy link
Collaborator

1uc commented Aug 29, 2023 •
edited
Loading

The generated code will inhibit the compiler from optimizing. For example the assembly of the following lines doesn't extract any common subexpressions:

                    nmodl_eigen_j[static_cast<int>(336)] = inst->dsqvol[id] * nt->_dt * ( -nmodl_eigen_x[static_cast<int>(0)] * inst->global->nR1 - nmodl_eigen_x[static_cast<int>(0)] * inst->global->nT1 + inst->global->nT2) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(357)] = inst->dsqvol[id] * nt->_dt * ( -nmodl_eigen_x[static_cast<int>(0)] * inst->global->nT1 + inst->global->nR2) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);
                    nmodl_eigen_j[static_cast<int>(378)] = inst->dsqvol[id] * nt->_dt * ( -nmodl_eigen_x[static_cast<int>(0)] * inst->global->nR1 + inst->global->nR2 + inst->global->nT2) / (pow(inst->diam[indexes[4*pnodecount + id]], 2.0) * inst->vrat[id]);

This is due to the indirection in the terms, e.g. inst->dsqvol[id], nt->_dt, inst->vrat[id], etc. Due to the possibility of aliasing, the compiler wont optimize. (Confirmed by manually eliminating the possibility of aliasing; and comparing the generated assembly.) By additionally eliminating the division term; adding a memset and removing the assignment of zero; and sorting the lines leads to substantial reduction in runtime (nearly 2x); by further adding -ffast-math when compiling the mod file gives another boost (just over 2x). (Speedup is of course only for this one function.)

Despite this, the runtime of this MOD file is dominated by the Newton solver which calls this function; the overall improvement to runtime is very small for this MOD file.

@1uc
Copy link
Collaborator

1uc commented Aug 22, 2024

The details of how KINETIC blocks work have changed enough to turn this issue stale.

@1uc 1uc closed this as completed Aug 22, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees

@cattabiani cattabiani

Labels
codegen Code generation backend performance Related to performance improvement solver Solver and numerical methods
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
4 participants
@ohm314 @pramodk @cattabiani @1uc