From ca688b0e268a7f67fc263a33a25db4390b5fb52d Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Tue, 5 Sep 2023 17:46:50 +0200 Subject: [PATCH 1/8] One less sparse13 --- src/nrniv/kschan.cpp | 69 ++++++++++++++------------------------------ src/nrniv/kschan.h | 4 +-- 2 files changed, 24 insertions(+), 49 deletions(-) diff --git a/src/nrniv/kschan.cpp b/src/nrniv/kschan.cpp index 655d4b376a..639678e758 100644 --- a/src/nrniv/kschan.cpp +++ b/src/nrniv/kschan.cpp @@ -10,6 +10,7 @@ #include "ocnotify.h" #include "parse.hpp" #include "nrniv_mf.h" +#include "ocmatrix.h" #define NSingleIndex 0 #if defined(__MWERKS__) && !defined(_MSC_VER) @@ -22,7 +23,6 @@ static KSChanList* channels; extern char* hoc_symbol_units(Symbol*, const char*); extern void nrn_mk_table_check(); -extern spREAL* spGetElement(char*, int, int); static Symbol* ksstate_sym; static Symbol* ksgate_sym; @@ -880,7 +880,7 @@ KSChan::KSChan(Object* obj, bool is_p) { Sprintf(buf, "Chan%d", obj_->index); name_ = buf; ion_ = "NonSpecific"; - mat_ = NULL; + mat_ = nullptr; elms_ = NULL; diag_ = NULL; gmax_deflt_ = 0.; @@ -1277,10 +1277,9 @@ void KSChan::free1() { ligands_ = NULL; } if (mat_) { - spDestroy(mat_); + mat_ = nullptr; delete[] elms_; delete[] diag_; - mat_ = NULL; } ngate_ = 0; nstate_ = 0; @@ -2109,44 +2108,30 @@ void KSChan::freesym(Symbol* s, Symbol* top) { } void KSChan::setupmat() { - int i, j, err; // printf("KSChan::setupmat nksstate=%d\n", nksstate_); - if (mat_) { - spDestroy(mat_); - delete[] elms_; - delete[] diag_; - mat_ = NULL; - } if (!nksstate_) { return; } - mat_ = spCreate(nksstate_, 0, &err); - if (err != spOKAY) { - hoc_execerror("Couldn't create sparse matrix", 0); - } - spFactor(mat_); // will fail but creates an internal vector needed by - // mulmat which might be called prior to initialization - // when switching to cvode active. + mat_ = std::make_unique(nksstate_, nksstate_); elms_ = new double*[4 * (ntrans_ - ivkstrans_)]; diag_ = new double*[nksstate_]; - for (i = ivkstrans_, j = 0; i < ntrans_; ++i) { - int s, t; - s = trans_[i].src_ - nhhstate_ + 1; - t = trans_[i].target_ - nhhstate_ + 1; - elms_[j++] = spGetElement(mat_, s, s); - elms_[j++] = spGetElement(mat_, s, t); - elms_[j++] = spGetElement(mat_, t, t); - elms_[j++] = spGetElement(mat_, t, s); + for (int i = ivkstrans_, j = 0; i < ntrans_; ++i) { + int s = trans_[i].src_ - nhhstate_; + int t = trans_[i].target_ - nhhstate_; + elms_[j++] = mat_->mep(s, s); + elms_[j++] = mat_->mep(s, t); + elms_[j++] = mat_->mep(t, t); + elms_[j++] = mat_->mep(t, s); } - for (i = 0; i < nksstate_; ++i) { - diag_[i] = spGetElement(mat_, i + 1, i + 1); + for (int i = 0; i < nksstate_; ++i) { + diag_[i] = mat_->mep(i, i); } } void KSChan::fillmat(double v, Datum* pd) { int i, j; double a, b; - spClear(mat_); + mat_->zero(); for (i = ivkstrans_, j = 0; i < iligtrans_; ++i) { trans_[i].ab(v, a, b); // printf("trans %d v=%g a=%g b=%g\n", i, v, a, b); @@ -2164,7 +2149,6 @@ void KSChan::fillmat(double v, Datum* pd) { *elms_[j++] += a; } // printf("after fill\n"); - // spPrint(mat_, 0, 1, 0); } void KSChan::mat_dt(double dt, Memb_list* ml, std::size_t instance, std::size_t offset) { @@ -2181,26 +2165,15 @@ void KSChan::mat_dt(double dt, Memb_list* ml, std::size_t instance, std::size_t void KSChan::solvemat(Memb_list* ml, std::size_t instance, std::size_t offset) { // spSolve seems to require that the parameters are contiguous, which // they're not anymore in the real NEURON data structure - std::vector s(nksstate_ + 1); // +1 so the pointer arithmetic to account for 1-based - // indexing is valid + std::vector s(nksstate_); for (auto j = 0; j < nksstate_; ++j) { - s[j + 1] = ml->data(instance, offset + j); - } - auto const e = spFactor(mat_); - if (e != spOKAY) { - switch (e) { - case spZERO_DIAG: - hoc_execerror("spFactor error:", "Zero Diagonal"); - case spNO_MEMORY: - hoc_execerror("spFactor error:", "No Memory"); - case spSINGULAR: - hoc_execerror("spFactor error:", "Singular"); - } + s[j] = ml->data(instance, offset + j); } - spSolve(mat_, s.data(), s.data()); + IvocVect in{}; in.vec() = s; + mat_->solv(&in, &in, false); // Propgate the solution back to the mechanism data for (auto j = 0; j < nksstate_; ++j) { - ml->data(instance, offset + j) = s[j + 1]; + ml->data(instance, offset + j) = s[j]; } } @@ -2216,7 +2189,9 @@ void KSChan::mulmat(Memb_list* ml, s[j + 1] = ml->data(instance, offset_s + j); ds[j + 1] = ml->data(instance, offset_ds + j); } - spMultiply(mat_, ds.data(), s.data()); + IvocVect in{}; in.vec() = ds; + IvocVect out{}; out.vec() = s; + mat_->mulv(&in, &out); // Propagate the results for (auto j = 0; j < nksstate_; ++j) { ml->data(instance, offset_s + j) = s[j + 1]; diff --git a/src/nrniv/kschan.h b/src/nrniv/kschan.h index de49789be5..5cc4ba997e 100644 --- a/src/nrniv/kschan.h +++ b/src/nrniv/kschan.h @@ -6,7 +6,7 @@ #include "ivocvect.h" #include "nrnunits.h" -#include "spmatrix.h" +#include "ocmatrix.h" // extern double dt; extern double celsius; @@ -468,7 +468,7 @@ class KSChan { int cvode_ieq_; Symbol* mechsym_; // the top level symbol (insert sym or new sym) Symbol* rlsym_; // symbol with the range list (= mechsym_ when density) - char* mat_; + std::unique_ptr mat_; double** elms_; double** diag_; int dsize_; // size of prop->dparam From 8a05c9e74981882277c92a48d91f5b7f5a2dcbfa Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Tue, 28 Nov 2023 08:14:20 +0100 Subject: [PATCH 2/8] Debug --- src/nrniv/kschan.cpp | 40 +++++++++++++++++++++------------------- src/nrniv/kschan.h | 2 +- 2 files changed, 22 insertions(+), 20 deletions(-) diff --git a/src/nrniv/kschan.cpp b/src/nrniv/kschan.cpp index 639678e758..b339c9d77e 100644 --- a/src/nrniv/kschan.cpp +++ b/src/nrniv/kschan.cpp @@ -2112,7 +2112,7 @@ void KSChan::setupmat() { if (!nksstate_) { return; } - mat_ = std::make_unique(nksstate_, nksstate_); + mat_ = std::make_unique(nksstate_, nksstate_); elms_ = new double*[4 * (ntrans_ - ivkstrans_)]; diag_ = new double*[nksstate_]; for (int i = ivkstrans_, j = 0; i < ntrans_; ++i) { @@ -2129,10 +2129,10 @@ void KSChan::setupmat() { } void KSChan::fillmat(double v, Datum* pd) { - int i, j; - double a, b; mat_->zero(); - for (i = ivkstrans_, j = 0; i < iligtrans_; ++i) { + int j = 0; + for (int i = ivkstrans_; i < iligtrans_; ++i) { + double a, b; trans_[i].ab(v, a, b); // printf("trans %d v=%g a=%g b=%g\n", i, v, a, b); *elms_[j++] -= a; @@ -2140,9 +2140,9 @@ void KSChan::fillmat(double v, Datum* pd) { *elms_[j++] -= b; *elms_[j++] += a; } - for (i = iligtrans_; i < ntrans_; ++i) { - a = trans_[i].alpha(pd); - b = trans_[i].beta(); + for (int i = iligtrans_; i < ntrans_; ++i) { + double a = trans_[i].alpha(pd); + double b = trans_[i].beta(); *elms_[j++] -= a; *elms_[j++] += b; *elms_[j++] -= b; @@ -2154,7 +2154,6 @@ void KSChan::fillmat(double v, Datum* pd) { void KSChan::mat_dt(double dt, Memb_list* ml, std::size_t instance, std::size_t offset) { // y' = m*y this part add the dt for the form ynew/dt - yold/dt =m*ynew // the matrix ends up as (m-1/dt)ynew = -1/dt*yold - int i; double dt1 = -1. / dt; for (int i = 0; i < nksstate_; ++i) { *(diag_[i]) += dt1; @@ -2165,15 +2164,18 @@ void KSChan::mat_dt(double dt, Memb_list* ml, std::size_t instance, std::size_t void KSChan::solvemat(Memb_list* ml, std::size_t instance, std::size_t offset) { // spSolve seems to require that the parameters are contiguous, which // they're not anymore in the real NEURON data structure - std::vector s(nksstate_); + std::vector s(nksstate_ + 1); for (auto j = 0; j < nksstate_; ++j) { - s[j] = ml->data(instance, offset + j); + s[j + 1] = ml->data(instance, offset + j); } - IvocVect in{}; in.vec() = s; - mat_->solv(&in, &in, false); + IvocVect inout{}; + inout.vec() = std::vector(s.begin(), s.end() - 1); + mat_->solv(&inout, &inout, false); + std::copy(inout.vec().begin(), inout.vec().end(), s.begin()); + // Propgate the solution back to the mechanism data for (auto j = 0; j < nksstate_; ++j) { - ml->data(instance, offset + j) = s[j]; + ml->data(instance, offset + j) = s[j + 1]; } } @@ -2181,17 +2183,17 @@ void KSChan::mulmat(Memb_list* ml, std::size_t instance, std::size_t offset_s, std::size_t offset_ds) { - std::vector s, ds; - s.resize(nksstate_ + 1); // +1 so the pointer arithmetic to account for 1-based indexing is - // valid - ds.resize(nksstate_ + 1); + std::vector s(nksstate_+1), ds(nksstate_+1); for (auto j = 0; j < nksstate_; ++j) { s[j + 1] = ml->data(instance, offset_s + j); ds[j + 1] = ml->data(instance, offset_ds + j); } - IvocVect in{}; in.vec() = ds; - IvocVect out{}; out.vec() = s; + IvocVect in{}, out{}; + in.vec() = std::vector(ds.begin(), ds.end() - 1); + out.vec() = std::vector(s.begin(), s.end() - 1); mat_->mulv(&in, &out); + std::copy(in.vec().begin(), in.vec().end(), ds.begin()); + std::copy(out.vec().begin(), out.vec().end(), s.begin()); // Propagate the results for (auto j = 0; j < nksstate_; ++j) { ml->data(instance, offset_s + j) = s[j + 1]; diff --git a/src/nrniv/kschan.h b/src/nrniv/kschan.h index 5cc4ba997e..8f274ce39f 100644 --- a/src/nrniv/kschan.h +++ b/src/nrniv/kschan.h @@ -468,7 +468,7 @@ class KSChan { int cvode_ieq_; Symbol* mechsym_; // the top level symbol (insert sym or new sym) Symbol* rlsym_; // symbol with the range list (= mechsym_ when density) - std::unique_ptr mat_; + std::unique_ptr mat_; double** elms_; double** diag_; int dsize_; // size of prop->dparam From 72f7558728bfd01a4f0c9f0eede609bbf85ca103 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Tue, 28 Nov 2023 17:55:16 +0100 Subject: [PATCH 3/8] Fix --- src/ivoc/ocmatrix.cpp | 2 +- src/nrniv/kschan.cpp | 55 +++++++++++++++++++------------------------ 2 files changed, 25 insertions(+), 32 deletions(-) diff --git a/src/ivoc/ocmatrix.cpp b/src/ivoc/ocmatrix.cpp index 2823c58b50..a702470395 100644 --- a/src/ivoc/ocmatrix.cpp +++ b/src/ivoc/ocmatrix.cpp @@ -223,7 +223,7 @@ void OcFullMatrix::solv(Vect* in, Vect* out, bool use_lu) { } auto v1 = Vect2VEC(in); auto v2 = Vect2VEC(out); - v2 = lu_->solve(v1); + v2 = lu_->solve(v1).eval(); } double OcFullMatrix::det(int* e) { diff --git a/src/nrniv/kschan.cpp b/src/nrniv/kschan.cpp index a57f574308..c28bb01f31 100644 --- a/src/nrniv/kschan.cpp +++ b/src/nrniv/kschan.cpp @@ -2115,12 +2115,13 @@ void KSChan::setupmat() { elms_ = new double*[4 * (ntrans_ - ivkstrans_)]; diag_ = new double*[nksstate_]; for (int i = ivkstrans_, j = 0; i < ntrans_; ++i) { - int s = trans_[i].src_ - nhhstate_; - int t = trans_[i].target_ - nhhstate_; - elms_[j++] = mat_->mep(s, s); - elms_[j++] = mat_->mep(s, t); - elms_[j++] = mat_->mep(t, t); - elms_[j++] = mat_->mep(t, s); + int s, t; + s = trans_[i].src_ - nhhstate_; + t = trans_[i].target_ - nhhstate_; + elms_[j++] = mat_->mep(s, s); + elms_[j++] = mat_->mep(s, t); + elms_[j++] = mat_->mep(t, t); + elms_[j++] = mat_->mep(t, s); } for (int i = 0; i < nksstate_; ++i) { diag_[i] = mat_->mep(i, i); @@ -2128,10 +2129,10 @@ void KSChan::setupmat() { } void KSChan::fillmat(double v, Datum* pd) { + int i, j; + double a, b; mat_->zero(); - int j = 0; - for (int i = ivkstrans_; i < iligtrans_; ++i) { - double a, b; + for (i = ivkstrans_, j = 0; i < iligtrans_; ++i) { trans_[i].ab(v, a, b); // printf("trans %d v=%g a=%g b=%g\n", i, v, a, b); *elms_[j++] -= a; @@ -2139,9 +2140,9 @@ void KSChan::fillmat(double v, Datum* pd) { *elms_[j++] -= b; *elms_[j++] += a; } - for (int i = iligtrans_; i < ntrans_; ++i) { - double a = trans_[i].alpha(pd); - double b = trans_[i].beta(); + for (i = iligtrans_; i < ntrans_; ++i) { + a = trans_[i].alpha(pd); + b = trans_[i].beta(); *elms_[j++] -= a; *elms_[j++] += b; *elms_[j++] -= b; @@ -2163,18 +2164,15 @@ void KSChan::mat_dt(double dt, Memb_list* ml, std::size_t instance, std::size_t void KSChan::solvemat(Memb_list* ml, std::size_t instance, std::size_t offset) { // spSolve seems to require that the parameters are contiguous, which // they're not anymore in the real NEURON data structure - std::vector s(nksstate_ + 1); + IvocVect s(nksstate_), out(nksstate_); for (auto j = 0; j < nksstate_; ++j) { - s[j + 1] = ml->data(instance, offset + j); + s[j] = ml->data(instance, offset + j); } - IvocVect inout{}; - inout.vec() = std::vector(s.begin(), s.end() - 1); - mat_->solv(&inout, &inout, false); - std::copy(inout.vec().begin(), inout.vec().end(), s.begin()); + mat_->solv(&s, &out, false); // Propgate the solution back to the mechanism data for (auto j = 0; j < nksstate_; ++j) { - ml->data(instance, offset + j) = s[j + 1]; + ml->data(instance, offset + j) = out[j]; } } @@ -2182,21 +2180,16 @@ void KSChan::mulmat(Memb_list* ml, std::size_t instance, std::size_t offset_s, std::size_t offset_ds) { - std::vector s(nksstate_+1), ds(nksstate_+1); + IvocVect s(nksstate_), ds(nksstate_); for (auto j = 0; j < nksstate_; ++j) { - s[j + 1] = ml->data(instance, offset_s + j); - ds[j + 1] = ml->data(instance, offset_ds + j); - } - IvocVect in{}, out{}; - in.vec() = std::vector(ds.begin(), ds.end() - 1); - out.vec() = std::vector(s.begin(), s.end() - 1); - mat_->mulv(&in, &out); - std::copy(in.vec().begin(), in.vec().end(), ds.begin()); - std::copy(out.vec().begin(), out.vec().end(), s.begin()); + s[j] = ml->data(instance, offset_s + j); + ds[j] = ml->data(instance, offset_ds + j); + } + mat_->mulv(&ds, &s); // Propagate the results for (auto j = 0; j < nksstate_; ++j) { - ml->data(instance, offset_s + j) = s[j + 1]; - ml->data(instance, offset_ds + j) = ds[j + 1]; + ml->data(instance, offset_s + j) = s[j]; + ml->data(instance, offset_ds + j) = ds[j]; } } From 8eab5953af2358d8a7464e7c7146b4cfbbc7b1aa Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Tue, 12 Dec 2023 17:10:20 +0100 Subject: [PATCH 4/8] Simplify matrix stuff and fix a weird bug where rhs = mat * solution --- src/nrniv/kschan.cpp | 99 +++++++++++++++++++++----------------------- src/nrniv/kschan.h | 5 ++- 2 files changed, 51 insertions(+), 53 deletions(-) diff --git a/src/nrniv/kschan.cpp b/src/nrniv/kschan.cpp index c28bb01f31..f787bdaa8e 100644 --- a/src/nrniv/kschan.cpp +++ b/src/nrniv/kschan.cpp @@ -9,7 +9,6 @@ #include "ocnotify.h" #include "parse.hpp" #include "nrniv_mf.h" -#include "ocmatrix.h" #define NSingleIndex 0 #if defined(__MWERKS__) && !defined(_MSC_VER) @@ -879,9 +878,6 @@ KSChan::KSChan(Object* obj, bool is_p) { Sprintf(buf, "Chan%d", obj_->index); name_ = buf; ion_ = "NonSpecific"; - mat_ = nullptr; - elms_ = NULL; - diag_ = NULL; gmax_deflt_ = 0.; erev_deflt_ = 0.; soffset_ = 4; // gmax, e, g, i before the first state in p array @@ -1275,11 +1271,7 @@ void KSChan::free1() { delete[] ligands_; ligands_ = NULL; } - if (mat_) { - mat_ = nullptr; - delete[] elms_; - delete[] diag_; - } + mat_.setZero(); ngate_ = 0; nstate_ = 0; ntrans_ = 0; @@ -2111,42 +2103,34 @@ void KSChan::setupmat() { if (!nksstate_) { return; } - mat_ = std::make_unique(nksstate_, nksstate_); - elms_ = new double*[4 * (ntrans_ - ivkstrans_)]; - diag_ = new double*[nksstate_]; - for (int i = ivkstrans_, j = 0; i < ntrans_; ++i) { - int s, t; - s = trans_[i].src_ - nhhstate_; - t = trans_[i].target_ - nhhstate_; - elms_[j++] = mat_->mep(s, s); - elms_[j++] = mat_->mep(s, t); - elms_[j++] = mat_->mep(t, t); - elms_[j++] = mat_->mep(t, s); - } - for (int i = 0; i < nksstate_; ++i) { - diag_[i] = mat_->mep(i, i); - } + mat_ = Eigen::SparseMatrix(nksstate_, nksstate_); } void KSChan::fillmat(double v, Datum* pd) { - int i, j; - double a, b; - mat_->zero(); - for (i = ivkstrans_, j = 0; i < iligtrans_; ++i) { - trans_[i].ab(v, a, b); + mat_.setZero(); + int j = 0; + for (int i = ivkstrans_; i < iligtrans_; ++i, ++j) { + auto& trans = trans_[i]; + double a, b; + trans.ab(v, a, b); // printf("trans %d v=%g a=%g b=%g\n", i, v, a, b); - *elms_[j++] -= a; - *elms_[j++] += b; - *elms_[j++] -= b; - *elms_[j++] += a; - } - for (i = iligtrans_; i < ntrans_; ++i) { - a = trans_[i].alpha(pd); - b = trans_[i].beta(); - *elms_[j++] -= a; - *elms_[j++] += b; - *elms_[j++] -= b; - *elms_[j++] += a; + auto s = trans.src_ - nhhstate_; + auto t = trans.target_ - nhhstate_; + mat_.coeffRef(s, s) -= a; + mat_.coeffRef(s, t) += b; + mat_.coeffRef(t, t) -= b; + mat_.coeffRef(t, s) += a; + } + for (int i = iligtrans_; i < ntrans_; ++i, ++j) { + auto& trans = trans_[i]; + auto s = trans.src_ - nhhstate_; + auto t = trans.target_ - nhhstate_; + double a = trans.alpha(pd); + double b = trans.beta(); + mat_.coeffRef(s, s) -= a; + mat_.coeffRef(s, t) += b; + mat_.coeffRef(t, t) -= b; + mat_.coeffRef(t, s) += a; } // printf("after fill\n"); } @@ -2156,23 +2140,35 @@ void KSChan::mat_dt(double dt, Memb_list* ml, std::size_t instance, std::size_t // the matrix ends up as (m-1/dt)ynew = -1/dt*yold double dt1 = -1. / dt; for (int i = 0; i < nksstate_; ++i) { - *(diag_[i]) += dt1; + mat_.coeffRef(i, i) += dt1; ml->data(instance, offset + i) *= dt1; } } void KSChan::solvemat(Memb_list* ml, std::size_t instance, std::size_t offset) { - // spSolve seems to require that the parameters are contiguous, which - // they're not anymore in the real NEURON data structure - IvocVect s(nksstate_), out(nksstate_); + std::vector s(nksstate_); for (auto j = 0; j < nksstate_; ++j) { s[j] = ml->data(instance, offset + j); } - mat_->solv(&s, &out, false); + mat_.makeCompressed(); + lu_.compute(mat_); + switch (lu_.info()) { + case Eigen::NumericalIssue: + hoc_execerror("NumericalIssue: The matrix is not valid following what expect Eigen SparseLU", nullptr); + break; + case Eigen::NoConvergence: + hoc_execerror("NoConvergence: The matrix did not converge", nullptr); + break; + case Eigen::InvalidInput: + hoc_execerror("InvalidInput: the inputs are invliad", nullptr); + break; + } + auto s_ = Eigen::Map>(s.data(), s.size()); + s_ = lu_.solve(s_); // Propgate the solution back to the mechanism data for (auto j = 0; j < nksstate_; ++j) { - ml->data(instance, offset + j) = out[j]; + ml->data(instance, offset + j) = s[j]; } } @@ -2180,15 +2176,16 @@ void KSChan::mulmat(Memb_list* ml, std::size_t instance, std::size_t offset_s, std::size_t offset_ds) { - IvocVect s(nksstate_), ds(nksstate_); + std::vector s(nksstate_); + std::vector ds(nksstate_); for (auto j = 0; j < nksstate_; ++j) { s[j] = ml->data(instance, offset_s + j); - ds[j] = ml->data(instance, offset_ds + j); } - mat_->mulv(&ds, &s); + auto s_ = Eigen::Map>(s.data(), s.size()); + auto ds_ = Eigen::Map>(ds.data(), ds.size()); + ds_ = mat_ * s_; // Propagate the results for (auto j = 0; j < nksstate_; ++j) { - ml->data(instance, offset_s + j) = s[j]; ml->data(instance, offset_ds + j) = ds[j]; } } diff --git a/src/nrniv/kschan.h b/src/nrniv/kschan.h index 8f274ce39f..58bad8e2f4 100644 --- a/src/nrniv/kschan.h +++ b/src/nrniv/kschan.h @@ -6,7 +6,7 @@ #include "ivocvect.h" #include "nrnunits.h" -#include "ocmatrix.h" +#include // extern double dt; extern double celsius; @@ -468,7 +468,8 @@ class KSChan { int cvode_ieq_; Symbol* mechsym_; // the top level symbol (insert sym or new sym) Symbol* rlsym_; // symbol with the range list (= mechsym_ when density) - std::unique_ptr mat_; + Eigen::SparseMatrix mat_{}; + Eigen::SparseLU> lu_{}; double** elms_; double** diag_; int dsize_; // size of prop->dparam From 5815e28e0a44c6bdbf10f9af46e2ad98768a99f8 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Tue, 12 Dec 2023 17:27:09 +0100 Subject: [PATCH 5/8] Resize instead of creating a new one --- src/nrniv/kschan.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/nrniv/kschan.cpp b/src/nrniv/kschan.cpp index f787bdaa8e..f2bb3fdf84 100644 --- a/src/nrniv/kschan.cpp +++ b/src/nrniv/kschan.cpp @@ -2103,7 +2103,7 @@ void KSChan::setupmat() { if (!nksstate_) { return; } - mat_ = Eigen::SparseMatrix(nksstate_, nksstate_); + mat_.resize(nksstate_, nksstate_); } void KSChan::fillmat(double v, Datum* pd) { From 15ec90ec751f88cf764ec9d4fbe779c142ecc0aa Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Tue, 12 Dec 2023 17:28:15 +0100 Subject: [PATCH 6/8] Format --- src/nrniv/kschan.cpp | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/src/nrniv/kschan.cpp b/src/nrniv/kschan.cpp index f2bb3fdf84..f9b0efd2d5 100644 --- a/src/nrniv/kschan.cpp +++ b/src/nrniv/kschan.cpp @@ -2153,15 +2153,17 @@ void KSChan::solvemat(Memb_list* ml, std::size_t instance, std::size_t offset) { mat_.makeCompressed(); lu_.compute(mat_); switch (lu_.info()) { - case Eigen::NumericalIssue: - hoc_execerror("NumericalIssue: The matrix is not valid following what expect Eigen SparseLU", nullptr); - break; - case Eigen::NoConvergence: - hoc_execerror("NoConvergence: The matrix did not converge", nullptr); - break; - case Eigen::InvalidInput: - hoc_execerror("InvalidInput: the inputs are invliad", nullptr); - break; + case Eigen::NumericalIssue: + hoc_execerror( + "NumericalIssue: The matrix is not valid following what expect Eigen SparseLU", + nullptr); + break; + case Eigen::NoConvergence: + hoc_execerror("NoConvergence: The matrix did not converge", nullptr); + break; + case Eigen::InvalidInput: + hoc_execerror("InvalidInput: the inputs are invliad", nullptr); + break; } auto s_ = Eigen::Map>(s.data(), s.size()); s_ = lu_.solve(s_); From 6226178b35dce094ae7418909eff711160cf3a54 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Wed, 13 Dec 2023 11:40:53 +0100 Subject: [PATCH 7/8] last clean? --- src/nrniv/kschan.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/nrniv/kschan.cpp b/src/nrniv/kschan.cpp index f9b0efd2d5..aab3c3c2e4 100644 --- a/src/nrniv/kschan.cpp +++ b/src/nrniv/kschan.cpp @@ -2146,13 +2146,16 @@ void KSChan::mat_dt(double dt, Memb_list* ml, std::size_t instance, std::size_t } void KSChan::solvemat(Memb_list* ml, std::size_t instance, std::size_t offset) { - std::vector s(nksstate_); + Eigen::VectorXd s(nksstate_); for (auto j = 0; j < nksstate_; ++j) { s[j] = ml->data(instance, offset + j); } mat_.makeCompressed(); lu_.compute(mat_); switch (lu_.info()) { + case Eigen::Success: + // Everything fine, at least no warning + break; case Eigen::NumericalIssue: hoc_execerror( "NumericalIssue: The matrix is not valid following what expect Eigen SparseLU", @@ -2165,8 +2168,7 @@ void KSChan::solvemat(Memb_list* ml, std::size_t instance, std::size_t offset) { hoc_execerror("InvalidInput: the inputs are invliad", nullptr); break; } - auto s_ = Eigen::Map>(s.data(), s.size()); - s_ = lu_.solve(s_); + s = lu_.solve(s); // Propgate the solution back to the mechanism data for (auto j = 0; j < nksstate_; ++j) { @@ -2178,14 +2180,12 @@ void KSChan::mulmat(Memb_list* ml, std::size_t instance, std::size_t offset_s, std::size_t offset_ds) { - std::vector s(nksstate_); - std::vector ds(nksstate_); + Eigen::VectorXd s(nksstate_); + Eigen::VectorXd ds(nksstate_); for (auto j = 0; j < nksstate_; ++j) { s[j] = ml->data(instance, offset_s + j); } - auto s_ = Eigen::Map>(s.data(), s.size()); - auto ds_ = Eigen::Map>(ds.data(), ds.size()); - ds_ = mat_ * s_; + ds = mat_ * s; // Propagate the results for (auto j = 0; j < nksstate_; ++j) { ml->data(instance, offset_ds + j) = ds[j]; From 211a7cbc818f7b3ea77d344dd9908c3e75f36cd9 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Mon, 8 Jan 2024 15:21:41 +0100 Subject: [PATCH 8/8] Change the tolerance to see if everything is magically fixed --- test/hoctests/tests/test_kschan.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/hoctests/tests/test_kschan.py b/test/hoctests/tests/test_kschan.py index 3774b37e56..00c7ef1184 100644 --- a/test/hoctests/tests/test_kschan.py +++ b/test/hoctests/tests/test_kschan.py @@ -328,7 +328,7 @@ def test_2(): h.cvode_active(1) # At least executes KSChan::mulmat hrun( - "kchan without single cvode=True", t_tol=2e-7, v_tol=1e-11, v_tol_per_time=5e-7 + "kchan without single cvode=True", t_tol=4e-7, v_tol=1e-11, v_tol_per_time=5e-7 ) h.cvode_active(0)