From 9181bf2f145dcf1d85d29485494495315c00dfa8 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 30 Nov 2023 12:01:28 +0100 Subject: [PATCH 01/16] nrnlinz.cpp: use eigen instead of sparse13 --- src/nrniv/nonlinz.cpp | 87 +++++++++++++++++++++---------------------- 1 file changed, 43 insertions(+), 44 deletions(-) diff --git a/src/nrniv/nonlinz.cpp b/src/nrniv/nonlinz.cpp index 65bb00dd5f..040646567a 100644 --- a/src/nrniv/nonlinz.cpp +++ b/src/nrniv/nonlinz.cpp @@ -9,6 +9,8 @@ #include "cspmatrix.h" #include "membfunc.h" +#include + extern void v_setup_vectors(); extern int nrndae_extra_eqn_count(); extern Symlist* hoc_built_in_symlist; @@ -29,14 +31,14 @@ class NonLinImpRep { void dsds(); int gapsolve(); - char* m_; + Eigen::SparseMatrix, Eigen::RowMajor> m_{}; int scnt_; // structure_change int n_v_, n_ext_, n_lin_, n_ode_, neq_v_, neq_; std::vector> pv_, pvdot_; int* v_index_; double* rv_; double* jv_; - double** diag_; + std::vector*> diag_; double* deltavec_; // just like cvode.atol*cvode.atolscale for ode's double delta_; // slightly more efficient and easier for v. void current(int, Memb_list*, int); @@ -152,7 +154,7 @@ void NonLinImp::compute(double omega, double deltafac, int maxiter) { rep_->omega_ = 1000. * omega; rep_->delta(deltafac); // fill matrix - cmplx_spClear(rep_->m_); + rep_->m_.setZero(); rep_->didv(); rep_->dsds(); #if 1 // when 0 equivalent to standard method @@ -160,20 +162,6 @@ void NonLinImp::compute(double omega, double deltafac, int maxiter) { rep_->dsdv(); #endif - // cmplx_spPrint(rep_->m_, 0, 1, 1); - // for (int i=0; i < rep_->neq_; ++i) { - // printf("i=%d %g %g\n", i, rep_->diag_[i][0], rep_->diag_[i][1]); - // } - int e = cmplx_spFactor(rep_->m_); - switch (e) { - case spZERO_DIAG: - hoc_execerror("cmplx_spFactor error:", "Zero Diagonal"); - case spNO_MEMORY: - hoc_execerror("cmplx_spFactor error:", "No Memory"); - case spSINGULAR: - hoc_execerror("cmplx_spFactor error:", "Singular"); - } - rep_->iloc_ = -2; } @@ -196,8 +184,20 @@ int NonLinImp::solve(int curloc) { if (nrnthread_v_transfer_) { rval = rep_->gapsolve(); } else { - assert(rep_->m_); - cmplx_spSolve(rep_->m_, rep_->rv_ - 1, rep_->rv_ - 1, rep_->jv_ - 1, rep_->jv_ - 1); + std::vector> rjv{}; + rjv.reserve(rep_->neq_); + for (size_t i = 0; i < rep_->neq_; ++i) { + rjv.emplace_back(rep_->rv_[i], rep_->jv_[i]); + } + auto lu = Eigen::SparseLUm_)>(rep_->m_); + lu.analyzePattern(rep_->m_); + lu.factorize(rep_->m_); + auto rjv_ = Eigen::Map, Eigen::Dynamic>>(rjv.data(), rjv.size()); + rjv_ = lu.solve(rjv_); + for (size_t i = 0; i < rep_->neq_; ++i) { + rep_->rv_[i] = rjv[i].real(); + rep_->jv_[i] = rjv[i].imag(); + } } } return rval; @@ -211,7 +211,6 @@ NonLinImpRep::NonLinImpRep() { int i, j, ieq, cnt; NrnThread* _nt = nrn_threads; maxiter_ = 500; - m_ = NULL; vsymtol_ = NULL; Symbol* vsym = hoc_table_lookup("v", hoc_built_in_symlist); @@ -245,8 +244,7 @@ NonLinImpRep::NonLinImpRep() { if (neq_ == 0) { return; } - m_ = cmplx_spCreate(neq_, 1, &err); - assert(err == spOKAY); + m_ = Eigen::SparseMatrix, Eigen::RowMajor>(neq_, neq_); pv_.resize(neq_); pvdot_.resize(neq_); v_index_ = new int[n_v_]; @@ -254,7 +252,7 @@ NonLinImpRep::NonLinImpRep() { rv_ += 1; jv_ = new double[neq_ + 1]; jv_ += 1; - diag_ = new double*[neq_]; + diag_.resize(neq_); deltavec_ = new double[neq_]; for (i = 0; i < n_v_; ++i) { @@ -265,23 +263,18 @@ NonLinImpRep::NonLinImpRep() { v_index_[i] = i + 1; } for (i = 0; i < n_v_; ++i) { - diag_[i] = cmplx_spGetElement(m_, v_index_[i], v_index_[i]); + diag_[i] = &m_.coeffRef(v_index_[i], v_index_[i]); } for (i = neq_v_; i < neq_; ++i) { - diag_[i] = cmplx_spGetElement(m_, i + 1, i + 1); + diag_[i] = &m_.coeffRef(i, i); } scnt_ = structure_change_cnt; } NonLinImpRep::~NonLinImpRep() { - if (!m_) { - return; - } - cmplx_spDestroy(m_); delete[] v_index_; delete[](rv_ - 1); delete[](jv_ - 1); - delete[] diag_; delete[] deltavec_; } @@ -317,10 +310,8 @@ void NonLinImpRep::didv() { for (i = _nt->ncell; i < n_v_; ++i) { nd = _nt->_v_node[i]; ip = _nt->_v_parent[i]->v_node_index; - double* a = cmplx_spGetElement(m_, v_index_[ip], v_index_[i]); - double* b = cmplx_spGetElement(m_, v_index_[i], v_index_[ip]); - *a += NODEA(nd); - *b += NODEB(nd); + m_.coeffRef(v_index_[ip], v_index_[i]) += NODEA(nd); + m_.coeffRef(v_index_[i], v_index_[ip]) += NODEB(nd); *diag_[i] -= NODEB(nd); *diag_[ip] -= NODEA(nd); } @@ -412,9 +403,7 @@ void NonLinImpRep::dids() { *pv_[is] = x1[is]; // restore s double g = (NODERHS(nd) - x2[in]) / deltavec_[is]; if (g != 0.) { - double* elm = - cmplx_spGetElement(m_, v_index_[nd->v_node_index], is + 1); - elm[0] = -g; + m_.coeffRef(v_index_[nd->v_node_index], is + 1) = -g; } } // don't know if this is necessary but make sure last @@ -477,9 +466,7 @@ void NonLinImpRep::dsdv() { for (is = ieq + in * cnt, iis = 0; iis < cnt; ++iis, ++is) { double ds = (x2[is] - *pvdot_[is]) / delta_; if (ds != 0.) { - double* elm = - cmplx_spGetElement(m_, is + 1, v_index_[nd->v_node_index]); - elm[0] = -ds; + m_.coeffRef(is + 1, v_index_[nd->v_node_index]) = -ds; } } } @@ -494,7 +481,7 @@ void NonLinImpRep::dsds() { NrnThread* nt = nrn_threads; // jw term for (i = neq_v_; i < neq_; ++i) { - diag_[i][1] += omega_; + *diag_[i] += std::complex(0, omega_); } ieq = neq_v_; for (NrnThreadMembList* tml = nt->tml; tml; tml = tml->next) { @@ -540,8 +527,7 @@ void NonLinImpRep::dsds() { for (is = ieq + in * cnt, iis = 0; iis < cnt; ++iis, ++is) { double ds = (*pvdot_[is] - x2[is]) / deltavec_[is]; if (ds != 0.) { - double* elm = cmplx_spGetElement(m_, is + 1, ks + 1); - elm[0] = -ds; + m_.coeffRef(is, ks) = -ds; } *pv_[ks] = x1[ks]; } @@ -620,7 +606,20 @@ int NonLinImpRep::gapsolve() { for (iter = 1; iter <= maxiter_; ++iter) { if (neq_) { - cmplx_spSolve(m_, rb - 1, rx1 - 1, jb - 1, jx1 - 1); + std::vector> rjv{}; + rjv.reserve(neq_); + for (size_t i = 0; i < neq_; ++i) { + rjv.emplace_back(rb[i], rx1[i]); + } + auto lu = Eigen::SparseLU(m_); + lu.analyzePattern(m_); + lu.factorize(m_); + auto rjv_ = Eigen::Map, Eigen::Dynamic>>(rjv.data(), rjv.size()); + rjv_ = lu.solve(rjv_); + for (size_t i = 0; i < neq_; ++i) { + jb[i] = rjv[i].real(); + jx1[i] = rjv[i].imag(); + } } // if any change in x > tol, then do another iteration. From 3dec29609fb7f88f7724535662235cc94a377b32 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Fri, 1 Dec 2023 09:05:08 +0100 Subject: [PATCH 02/16] To complex --- src/nrniv/nonlinz.cpp | 194 +++++++++++------------------------------ src/nrniv/nonlinz.h | 7 +- src/nrniv/partrans.cpp | 86 ++++++++++-------- 3 files changed, 105 insertions(+), 182 deletions(-) diff --git a/src/nrniv/nonlinz.cpp b/src/nrniv/nonlinz.cpp index 040646567a..94af378ef2 100644 --- a/src/nrniv/nonlinz.cpp +++ b/src/nrniv/nonlinz.cpp @@ -6,24 +6,22 @@ #include "nrniv_mf.h" #include "nrnoc2iv.h" #include "nrnmpi.h" -#include "cspmatrix.h" #include "membfunc.h" #include +#include extern void v_setup_vectors(); extern int nrndae_extra_eqn_count(); extern Symlist* hoc_built_in_symlist; extern void (*nrnthread_v_transfer_)(NrnThread*); -extern spREAL* spGetElement(char*, int, int); -extern void pargap_jacobi_rhs(double*, double*); +extern void pargap_jacobi_rhs(std::vector>&, const std::vector>&); extern void pargap_jacobi_setup(int mode); -class NonLinImpRep { +class NonLinImpRep final { public: NonLinImpRep(); - virtual ~NonLinImpRep(); void delta(double); void didv(); void dids(); @@ -32,14 +30,14 @@ class NonLinImpRep { int gapsolve(); Eigen::SparseMatrix, Eigen::RowMajor> m_{}; + std::unique_ptr> lu_{}; int scnt_; // structure_change int n_v_, n_ext_, n_lin_, n_ode_, neq_v_, neq_; std::vector> pv_, pvdot_; - int* v_index_; - double* rv_; - double* jv_; + std::vector v_index_; + std::vector> v_; std::vector*> diag_; - double* deltavec_; // just like cvode.atol*cvode.atolscale for ode's + std::vector deltavec_; // just like cvode.atol*cvode.atolscale for ode's double delta_; // slightly more efficient and easier for v. void current(int, Memb_list*, int); void ode(int, Memb_list*); @@ -50,14 +48,10 @@ class NonLinImpRep { int maxiter_; }; -NonLinImp::NonLinImp() { - rep_ = NULL; -} NonLinImp::~NonLinImp() { - if (rep_) { - delete rep_; - } + delete rep_; } + double NonLinImp::transfer_amp(int curloc, int vloc) { if (nrnmpi_numprocs > 1 && nrnthread_v_transfer_ && curloc != rep_->iloc_) { hoc_execerror( @@ -66,9 +60,7 @@ double NonLinImp::transfer_amp(int curloc, int vloc) { if (curloc != rep_->iloc_) { solve(curloc); } - double x = rep_->rv_[vloc]; - double y = rep_->jv_[vloc]; - return sqrt(x * x + y * y); + return std::abs(rep_->v_[vloc]); } double NonLinImp::input_amp(int curloc) { if (nrnmpi_numprocs > 1 && nrnthread_v_transfer_) { @@ -80,9 +72,7 @@ double NonLinImp::input_amp(int curloc) { if (curloc < 0) { return 0.0; } - double x = rep_->rv_[curloc]; - double y = rep_->jv_[curloc]; - return sqrt(x * x + y * y); + return std::abs(rep_->v_[curloc]); } double NonLinImp::transfer_phase(int curloc, int vloc) { if (nrnmpi_numprocs > 1 && nrnthread_v_transfer_ && curloc != rep_->iloc_) { @@ -92,9 +82,7 @@ double NonLinImp::transfer_phase(int curloc, int vloc) { if (curloc != rep_->iloc_) { solve(curloc); } - double x = rep_->rv_[vloc]; - double y = rep_->jv_[vloc]; - return atan2(y, x); + return std::arg(rep_->v_[vloc]); } double NonLinImp::input_phase(int curloc) { if (nrnmpi_numprocs > 1 && nrnthread_v_transfer_) { @@ -106,9 +94,7 @@ double NonLinImp::input_phase(int curloc) { if (curloc < 0) { return 0.0; } - double x = rep_->rv_[curloc]; - double y = rep_->jv_[curloc]; - return atan2(y, x); + return std::arg(rep_->v_[curloc]); } double NonLinImp::ratio_amp(int clmploc, int vloc) { if (nrnmpi_numprocs > 1 && nrnthread_v_transfer_) { @@ -120,22 +106,13 @@ double NonLinImp::ratio_amp(int clmploc, int vloc) { if (clmploc != rep_->iloc_) { solve(clmploc); } - double ax, bx, cx, ay, by, cy, bb; - ax = rep_->rv_[vloc]; - ay = rep_->jv_[vloc]; - bx = rep_->rv_[clmploc]; - by = rep_->jv_[clmploc]; - bb = bx * bx + by * by; - cx = (ax * bx + ay * by) / bb; - cy = (ay * bx - ax * by) / bb; - return sqrt(cx * cx + cy * cy); + return std::abs(rep_->v_[vloc] * std::conj(rep_->v_[clmploc]) / std::norm(rep_->v_[clmploc])); } void NonLinImp::compute(double omega, double deltafac, int maxiter) { v_setup_vectors(); nrn_rhs(nrn_ensure_model_data_are_sorted(), nrn_threads[0]); if (rep_ && rep_->scnt_ != structure_change_cnt) { - delete rep_; - rep_ = NULL; + rep_ = nullptr; } if (!rep_) { rep_ = new NonLinImpRep(); @@ -162,6 +139,10 @@ void NonLinImp::compute(double omega, double deltafac, int maxiter) { rep_->dsdv(); #endif + rep_->lu_ = std::make_uniquem_)>>(rep_->m_); + rep_->lu_->analyzePattern(rep_->m_); + rep_->lu_->factorize(rep_->m_); + rep_->iloc_ = -2; } @@ -175,29 +156,16 @@ int NonLinImp::solve(int curloc) { int i; rep_->iloc_ = curloc; for (i = 0; i < rep_->neq_; ++i) { - rep_->rv_[i] = 0; - rep_->jv_[i] = 0; + rep_->v_[i] = std::complex(0, 0); } if (curloc >= 0) { - rep_->rv_[curloc] = 1.e2 / NODEAREA(_nt->_v_node[curloc]); + rep_->v_[curloc] = 1.e2 / NODEAREA(_nt->_v_node[curloc]); } if (nrnthread_v_transfer_) { rval = rep_->gapsolve(); } else { - std::vector> rjv{}; - rjv.reserve(rep_->neq_); - for (size_t i = 0; i < rep_->neq_; ++i) { - rjv.emplace_back(rep_->rv_[i], rep_->jv_[i]); - } - auto lu = Eigen::SparseLUm_)>(rep_->m_); - lu.analyzePattern(rep_->m_); - lu.factorize(rep_->m_); - auto rjv_ = Eigen::Map, Eigen::Dynamic>>(rjv.data(), rjv.size()); - rjv_ = lu.solve(rjv_); - for (size_t i = 0; i < rep_->neq_; ++i) { - rep_->rv_[i] = rjv[i].real(); - rep_->jv_[i] = rjv[i].imag(); - } + auto rv = Eigen::Map, Eigen::Dynamic>>(rep_->v_.data(), rep_->v_.size()); + rv = rep_->lu_->solve(rv); } } return rval; @@ -247,13 +215,10 @@ NonLinImpRep::NonLinImpRep() { m_ = Eigen::SparseMatrix, Eigen::RowMajor>(neq_, neq_); pv_.resize(neq_); pvdot_.resize(neq_); - v_index_ = new int[n_v_]; - rv_ = new double[neq_ + 1]; - rv_ += 1; - jv_ = new double[neq_ + 1]; - jv_ += 1; + v_index_.resize(n_v_); + v_.resize(neq_); diag_.resize(neq_); - deltavec_ = new double[neq_]; + deltavec_.resize(neq_); for (i = 0; i < n_v_; ++i) { // utilize nd->eqn_index in case of use_sparse13 later @@ -271,13 +236,6 @@ NonLinImpRep::NonLinImpRep() { scnt_ = structure_change_cnt; } -NonLinImpRep::~NonLinImpRep() { - delete[] v_index_; - delete[](rv_ - 1); - delete[](jv_ - 1); - delete[] deltavec_; -} - void NonLinImpRep::delta(double deltafac) { // also defines pv_,pvdot_ map for ode int i, nc, cnt, ieq; NrnThread* nt = nrn_threads; @@ -293,7 +251,7 @@ void NonLinImpRep::delta(double deltafac) { // also defines pv_,pvdot_ map for nrn_ode_map_t ode_map = memb_func[i].ode_map; for (auto j = 0; j < nc; ++j) { ode_map( - ml->prop[j], ieq, pv_.data() + ieq, pvdot_.data() + ieq, deltavec_ + ieq, i); + ml->prop[j], ieq, pv_.data() + ieq, pvdot_.data() + ieq, deltavec_.data() + ieq, i); ieq += cnt; } } @@ -381,8 +339,6 @@ void NonLinImpRep::dids() { nrn_ode_count_t s = memb_func[i].ode_count; int cnt = (*s)(i); if (memb_func[i].current) { - double* x1 = rv_; // use as temporary storage - double* x2 = jv_; for (in = 0; in < ml->nodecount; ++in) { Node* nd = ml->nodelist[in]; // zero rhs @@ -390,18 +346,18 @@ void NonLinImpRep::dids() { // compute rhs current(i, ml, in); // save rhs - x2[in] = NODERHS(nd); + v_[in].imag(NODERHS(nd)); // each state incremented separately and restored for (iis = 0; iis < cnt; ++iis) { is = ieq + in * cnt + iis; // save s - x1[is] = *pv_[is]; + v_[is].real(*pv_[is]); // increment s and zero rhs *pv_[is] += deltavec_[is]; NODERHS(nd) = 0; current(i, ml, in); - *pv_[is] = x1[is]; // restore s - double g = (NODERHS(nd) - x2[in]) / deltavec_[is]; + *pv_[is] = v_[is].real(); // restore s + double g = (NODERHS(nd) - v_[in].imag()) / deltavec_[is]; if (g != 0.) { m_.coeffRef(v_index_[nd->v_node_index], is + 1) = -g; } @@ -428,22 +384,20 @@ void NonLinImpRep::dsdv() { nrn_ode_count_t s = memb_func[i].ode_count; int cnt = (*s)(i); if (memb_func[i].current) { - double* x1 = rv_; // use as temporary storage - double* x2 = jv_; // zero rhs, save v for (in = 0; in < ml->nodecount; ++in) { Node* nd = ml->nodelist[in]; for (is = ieq + in * cnt, iis = 0; iis < cnt; ++iis, ++is) { *pvdot_[is] = 0.; } - x1[in] = NODEV(nd); + v_[in].real(NODEV(nd)); } // increment v only once in case there are multiple // point processes at the same location for (in = 0; in < ml->nodecount; ++in) { Node* nd = ml->nodelist[in]; auto const v = nd->v(); - if (x1[in] == v) { + if (v_[in].real() == v) { nd->v() = v + delta_; } } @@ -453,10 +407,10 @@ void NonLinImpRep::dsdv() { for (in = 0; in < ml->nodecount; ++in) { Node* nd = ml->nodelist[in]; for (is = ieq + in * cnt, iis = 0; iis < cnt; ++iis, ++is) { - x2[is] = *pvdot_[is]; + v_[is].imag(*pvdot_[is]); *pvdot_[is] = 0; } - nd->v() = x1[in]; + nd->v() = v_[in].real(); } // compute the rhs(v) ode(i, ml); @@ -464,7 +418,7 @@ void NonLinImpRep::dsdv() { for (in = 0; in < ml->nodecount; ++in) { Node* nd = ml->nodelist[in]; for (is = ieq + in * cnt, iis = 0; iis < cnt; ++iis, ++is) { - double ds = (x2[is] - *pvdot_[is]) / delta_; + double ds = (v_[is].imag() - *pvdot_[is]) / delta_; if (ds != 0.) { m_.coeffRef(is + 1, v_index_[nd->v_node_index]) = -ds; } @@ -491,13 +445,11 @@ void NonLinImpRep::dsds() { int nc = ml->nodecount; nrn_ode_count_t s = memb_func[i].ode_count; int cnt = (*s)(i); - double* x1 = rv_; // use as temporary storage - double* x2 = jv_; // zero rhs, save s for (in = 0; in < ml->nodecount; ++in) { for (is = ieq + in * cnt, iis = 0; iis < cnt; ++iis, ++is) { *pvdot_[is] = 0.; - x1[is] = *pv_[is]; + v_[is].real(*pv_[is]); } } // compute rhs. this is the rhs(s) @@ -505,7 +457,7 @@ void NonLinImpRep::dsds() { // save rhs for (in = 0; in < ml->nodecount; ++in) { for (is = ieq + in * cnt, iis = 0; iis < cnt; ++iis, ++is) { - x2[is] = *pvdot_[is]; + v_[is].imag(*pvdot_[is]); } } // iterate over the states @@ -525,11 +477,11 @@ void NonLinImpRep::dsds() { Node* nd = ml->nodelist[in]; ks = ieq + in * cnt + kks; for (is = ieq + in * cnt, iis = 0; iis < cnt; ++iis, ++is) { - double ds = (*pvdot_[is] - x2[is]) / deltavec_[is]; + double ds = (*pvdot_[is] - v_[is].imag()) / deltavec_[is]; if (ds != 0.) { m_.coeffRef(is, ks) = -ds; } - *pv_[ks] = x1[ks]; + *pv_[ks] = v_[ks].real(); } } // perhaps not necessary but ensures the last computation is with @@ -560,8 +512,8 @@ void NonLinImpRep::ode(int im, Memb_list* ml) { // assume there is in fact an o } int NonLinImpRep::gapsolve() { - // On entry, rv_ and jv_ contain the complex b for A*x = b. - // On return rv_ and jv_ contain complex solution, x. + // On entry, v_ contains the complex b for A*x = b. + // On return v_ contains complex solution, x. // m_ is the factored matrix for the trees without gap junctions // Jacobi method (easy for parallel) // A = D + R @@ -580,22 +532,9 @@ int NonLinImpRep::gapsolve() { pargap_jacobi_setup(0); - double *rx, *jx, *rx1, *jx1, *rb, *jb; - if (neq_) { - rx = new double[neq_]; - jx = new double[neq_]; - rx1 = new double[neq_]; - jx1 = new double[neq_]; - rb = new double[neq_]; - jb = new double[neq_]; - } - - // initialize for first iteration - for (int i = 0; i < neq_; ++i) { - rx[i] = jx[i] = 0.0; - rb[i] = rv_[i]; - jb[i] = jv_[i]; - } + std::vector> rx(neq_); + std::vector> rx1(neq_); + std::vector> rb(v_.begin(), v_.end()); // iterate till change in x is small double tol = 1e-9; @@ -606,27 +545,17 @@ int NonLinImpRep::gapsolve() { for (iter = 1; iter <= maxiter_; ++iter) { if (neq_) { - std::vector> rjv{}; - rjv.reserve(neq_); - for (size_t i = 0; i < neq_; ++i) { - rjv.emplace_back(rb[i], rx1[i]); - } - auto lu = Eigen::SparseLU(m_); - lu.analyzePattern(m_); - lu.factorize(m_); - auto rjv_ = Eigen::Map, Eigen::Dynamic>>(rjv.data(), rjv.size()); - rjv_ = lu.solve(rjv_); - for (size_t i = 0; i < neq_; ++i) { - jb[i] = rjv[i].real(); - jx1[i] = rjv[i].imag(); - } + auto rb_ = Eigen::Map, Eigen::Dynamic>>(rb.data(), rb.size()); + auto rx1_ = Eigen::Map, Eigen::Dynamic>>(rx1.data(), rx1.size()); + rx1_ = lu_->solve(rb_); } // if any change in x > tol, then do another iteration. success = 1; delta = 0.0; for (int i = 0; i < neq_; ++i) { - double err = fabs(rx1[i] - rx[i]) + fabs(jx1[i] - jx[i]); + auto diff = rx1[i] - rx[i]; + double err = std::abs(diff.real()) + std::abs(diff.imag()); if (err > tol) { success = 0; } @@ -640,35 +569,18 @@ int NonLinImpRep::gapsolve() { } #endif if (success) { - for (int i = 0; i < neq_; ++i) { - rv_[i] = rx1[i]; - jv_[i] = jx1[i]; - } + v_ = rx1; break; } // setup for next iteration - for (int i = 0; i < neq_; ++i) { - rx[i] = rx1[i]; - jx[i] = jx1[i]; - rb[i] = rv_[i]; - jb[i] = jv_[i]; - } + rx = rx1; + rb = v_; pargap_jacobi_rhs(rb, rx); - pargap_jacobi_rhs(jb, jx); } pargap_jacobi_setup(1); // tear down - if (neq_) { - delete[] rx; - delete[] jx; - delete[] rx1; - delete[] jx1; - delete[] rb; - delete[] jb; - } - if (!success) { char buf[256]; Sprintf(buf, diff --git a/src/nrniv/nonlinz.h b/src/nrniv/nonlinz.h index 115a191b27..ca7a568f37 100644 --- a/src/nrniv/nonlinz.h +++ b/src/nrniv/nonlinz.h @@ -3,10 +3,9 @@ class NonLinImpRep; -class NonLinImp { +class NonLinImp final { public: - NonLinImp(); - virtual ~NonLinImp(); + ~NonLinImp(); void compute(double omega, double deltafac, int maxiter); double transfer_amp(int curloc, int vloc); // v_node[arg] is the node double transfer_phase(int curloc, int vloc); @@ -16,7 +15,7 @@ class NonLinImp { int solve(int curloc); private: - NonLinImpRep* rep_; + NonLinImpRep* rep_{}; }; #endif diff --git a/src/nrniv/partrans.cpp b/src/nrniv/partrans.cpp index 1f34f9bbf2..d902446d9b 100644 --- a/src/nrniv/partrans.cpp +++ b/src/nrniv/partrans.cpp @@ -13,6 +13,7 @@ #include #include +#include #include // Replaces NrnHash for MapSgid2Int #include #include @@ -873,47 +874,58 @@ void pargap_jacobi_setup(int mode) { } } -void pargap_jacobi_rhs(double* b, double* x) { - // helper for complex impedance with parallel gap junctions - // b = b - R*x R are the off diagonal gap elements of the jacobian. - // we presume 1 thread. First nrn_thread[0].end equations are in node order. - if (!nrnthread_v_transfer_) { - return; - } +void pargap_jacobi_rhs(std::vector>& b, const std::vector>& x) { + // First loop for real, second for imag + for (int real_imag = 0; real_imag < 2; ++real_imag) { + // helper for complex impedance with parallel gap junctions + // b = b - R*x R are the off diagonal gap elements of the jacobian. + // we presume 1 thread. First nrn_thread[0].end equations are in node order. + if (!nrnthread_v_transfer_) { + return; + } - NrnThread* _nt = nrn_threads; + NrnThread* _nt = nrn_threads; - // transfer gap node voltages to gap vpre - for (size_t i = 0; i < visources_.size(); ++i) { - Node* nd = visources_[i]; - nd->v() = x[nd->v_node_index]; - } - mpi_transfer(); - thread_transfer(_nt); + // transfer gap node voltages to gap vpre + for (size_t i = 0; i < visources_.size(); ++i) { + Node* nd = visources_[i]; + if (real_imag == 0) { + nd->v() = x[nd->v_node_index].real(); + } else { + nd->v() = x[nd->v_node_index].imag(); + } + } + mpi_transfer(); + thread_transfer(_nt); - // set gap node voltages to 0 so we can use nrn_cur to set rhs - for (size_t i = 0; i < visources_.size(); ++i) { - Node* nd = visources_[i]; - nd->v() = 0.0; - } - auto const sorted_token = nrn_ensure_model_data_are_sorted(); - auto* const vec_rhs = _nt->node_rhs_storage(); - // Initialize rhs to 0. - for (int i = 0; i < _nt->end; ++i) { - vec_rhs[i] = 0.0; - } - for (int k = 0; k < imped_current_type_count_; ++k) { - int type = imped_current_type_[k]; - Memb_list* ml = imped_current_ml_[k]; - memb_func[type].current(sorted_token, _nt, ml, type); - } + // set gap node voltages to 0 so we can use nrn_cur to set rhs + for (size_t i = 0; i < visources_.size(); ++i) { + Node* nd = visources_[i]; + nd->v() = 0.0; + } + auto const sorted_token = nrn_ensure_model_data_are_sorted(); + auto* const vec_rhs = _nt->node_rhs_storage(); + // Initialize rhs to 0. + for (int i = 0; i < _nt->end; ++i) { + vec_rhs[i] = 0.0; + } + for (int k = 0; k < imped_current_type_count_; ++k) { + int type = imped_current_type_[k]; + Memb_list* ml = imped_current_ml_[k]; + memb_func[type].current(sorted_token, _nt, ml, type); + } - // possibly many gap junctions in same node (and possible even different - // types) but rhs is the accumulation of all those instances at each node - // so ... The only thing that can go wrong is if there are intances of - // gap junctions that are not being used (not in the target list). - for (int i = 0; i < _nt->end; ++i) { - b[i] += vec_rhs[i]; + // possibly many gap junctions in same node (and possible even different + // types) but rhs is the accumulation of all those instances at each node + // so ... The only thing that can go wrong is if there are intances of + // gap junctions that are not being used (not in the target list). + for (int i = 0; i < _nt->end; ++i) { + if (real_imag == 0) { + b[i] += vec_rhs[i]; + } else { + b[i] += std::complex(0, vec_rhs[i]); + } + } } } From ad3155087242dfc5d4b61005e360094c17209d7f Mon Sep 17 00:00:00 2001 From: Luc Grosheintz Date: Fri, 1 Dec 2023 12:09:20 +0100 Subject: [PATCH 03/16] formatting. --- src/nrniv/nonlinz.cpp | 23 ++++++++++++++++------- src/nrniv/partrans.cpp | 3 ++- 2 files changed, 18 insertions(+), 8 deletions(-) diff --git a/src/nrniv/nonlinz.cpp b/src/nrniv/nonlinz.cpp index 94af378ef2..23bc09e566 100644 --- a/src/nrniv/nonlinz.cpp +++ b/src/nrniv/nonlinz.cpp @@ -16,7 +16,8 @@ extern int nrndae_extra_eqn_count(); extern Symlist* hoc_built_in_symlist; extern void (*nrnthread_v_transfer_)(NrnThread*); -extern void pargap_jacobi_rhs(std::vector>&, const std::vector>&); +extern void pargap_jacobi_rhs(std::vector>&, + const std::vector>&); extern void pargap_jacobi_setup(int mode); class NonLinImpRep final { @@ -38,7 +39,7 @@ class NonLinImpRep final { std::vector> v_; std::vector*> diag_; std::vector deltavec_; // just like cvode.atol*cvode.atolscale for ode's - double delta_; // slightly more efficient and easier for v. + double delta_; // slightly more efficient and easier for v. void current(int, Memb_list*, int); void ode(int, Memb_list*); @@ -164,7 +165,9 @@ int NonLinImp::solve(int curloc) { if (nrnthread_v_transfer_) { rval = rep_->gapsolve(); } else { - auto rv = Eigen::Map, Eigen::Dynamic>>(rep_->v_.data(), rep_->v_.size()); + auto rv = + Eigen::Map, Eigen::Dynamic>>(rep_->v_.data(), + rep_->v_.size()); rv = rep_->lu_->solve(rv); } } @@ -250,8 +253,12 @@ void NonLinImpRep::delta(double deltafac) { // also defines pv_,pvdot_ map for if (nrn_ode_count_t s = memb_func[i].ode_count; s && (cnt = s(i)) > 0) { nrn_ode_map_t ode_map = memb_func[i].ode_map; for (auto j = 0; j < nc; ++j) { - ode_map( - ml->prop[j], ieq, pv_.data() + ieq, pvdot_.data() + ieq, deltavec_.data() + ieq, i); + ode_map(ml->prop[j], + ieq, + pv_.data() + ieq, + pvdot_.data() + ieq, + deltavec_.data() + ieq, + i); ieq += cnt; } } @@ -545,8 +552,10 @@ int NonLinImpRep::gapsolve() { for (iter = 1; iter <= maxiter_; ++iter) { if (neq_) { - auto rb_ = Eigen::Map, Eigen::Dynamic>>(rb.data(), rb.size()); - auto rx1_ = Eigen::Map, Eigen::Dynamic>>(rx1.data(), rx1.size()); + auto rb_ = Eigen::Map, Eigen::Dynamic>>(rb.data(), + rb.size()); + auto rx1_ = Eigen::Map, Eigen::Dynamic>>(rx1.data(), + rx1.size()); rx1_ = lu_->solve(rb_); } diff --git a/src/nrniv/partrans.cpp b/src/nrniv/partrans.cpp index d902446d9b..3f5ea46182 100644 --- a/src/nrniv/partrans.cpp +++ b/src/nrniv/partrans.cpp @@ -874,7 +874,8 @@ void pargap_jacobi_setup(int mode) { } } -void pargap_jacobi_rhs(std::vector>& b, const std::vector>& x) { +void pargap_jacobi_rhs(std::vector>& b, + const std::vector>& x) { // First loop for real, second for imag for (int real_imag = 0; real_imag < 2; ++real_imag) { // helper for complex impedance with parallel gap junctions From 6e249c23c55ad4f0ae935be6d3addfd620a8e9fc Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Fri, 1 Dec 2023 15:49:20 +0100 Subject: [PATCH 04/16] Fix --- src/nrniv/nonlinz.cpp | 59 +++++++++++++++++-------------------------- 1 file changed, 23 insertions(+), 36 deletions(-) diff --git a/src/nrniv/nonlinz.cpp b/src/nrniv/nonlinz.cpp index 23bc09e566..3c1c2c5a70 100644 --- a/src/nrniv/nonlinz.cpp +++ b/src/nrniv/nonlinz.cpp @@ -35,9 +35,7 @@ class NonLinImpRep final { int scnt_; // structure_change int n_v_, n_ext_, n_lin_, n_ode_, neq_v_, neq_; std::vector> pv_, pvdot_; - std::vector v_index_; std::vector> v_; - std::vector*> diag_; std::vector deltavec_; // just like cvode.atol*cvode.atolscale for ode's double delta_; // slightly more efficient and easier for v. void current(int, Memb_list*, int); @@ -45,8 +43,8 @@ class NonLinImpRep final { double omega_; int iloc_; // current injection site of last solve - float* vsymtol_; - int maxiter_; + float* vsymtol_{}; + int maxiter_{500}; }; NonLinImp::~NonLinImp() { @@ -113,6 +111,7 @@ void NonLinImp::compute(double omega, double deltafac, int maxiter) { v_setup_vectors(); nrn_rhs(nrn_ensure_model_data_are_sorted(), nrn_threads[0]); if (rep_ && rep_->scnt_ != structure_change_cnt) { + delete rep_; rep_ = nullptr; } if (!rep_) { @@ -131,8 +130,10 @@ void NonLinImp::compute(double omega, double deltafac, int maxiter) { rep_->omega_ = 1000. * omega; rep_->delta(deltafac); - // fill matrix + rep_->m_.setZero(); + + // fill matrix rep_->didv(); rep_->dsds(); #if 1 // when 0 equivalent to standard method @@ -154,21 +155,18 @@ int NonLinImp::solve(int curloc) { hoc_execerror("Must call Impedance.compute first", 0); } if (rep_->iloc_ != curloc) { - int i; rep_->iloc_ = curloc; - for (i = 0; i < rep_->neq_; ++i) { - rep_->v_[i] = std::complex(0, 0); - } + rep_->v_ = std::vector>(rep_->neq_); if (curloc >= 0) { - rep_->v_[curloc] = 1.e2 / NODEAREA(_nt->_v_node[curloc]); + rep_->v_[curloc].real(1.e2 / NODEAREA(_nt->_v_node[curloc])); } if (nrnthread_v_transfer_) { rval = rep_->gapsolve(); } else { - auto rv = + auto v = Eigen::Map, Eigen::Dynamic>>(rep_->v_.data(), rep_->v_.size()); - rv = rep_->lu_->solve(rv); + v = rep_->lu_->solve(v); } } return rval; @@ -178,12 +176,9 @@ int NonLinImp::solve(int curloc) { // mapping is already done there. NonLinImpRep::NonLinImpRep() { - int err; - int i, j, ieq, cnt; + int i, cnt; NrnThread* _nt = nrn_threads; - maxiter_ = 500; - vsymtol_ = NULL; Symbol* vsym = hoc_table_lookup("v", hoc_built_in_symlist); if (vsym->extra) { vsymtol_ = &vsym->extra->tolerance; @@ -218,9 +213,7 @@ NonLinImpRep::NonLinImpRep() { m_ = Eigen::SparseMatrix, Eigen::RowMajor>(neq_, neq_); pv_.resize(neq_); pvdot_.resize(neq_); - v_index_.resize(n_v_); v_.resize(neq_); - diag_.resize(neq_); deltavec_.resize(neq_); for (i = 0; i < n_v_; ++i) { @@ -228,13 +221,6 @@ NonLinImpRep::NonLinImpRep() { Node* nd = _nt->_v_node[i]; pv_[i] = nd->v_handle(); pvdot_[i] = nd->rhs_handle(); - v_index_[i] = i + 1; - } - for (i = 0; i < n_v_; ++i) { - diag_[i] = &m_.coeffRef(v_index_[i], v_index_[i]); - } - for (i = neq_v_; i < neq_; ++i) { - diag_[i] = &m_.coeffRef(i, i); } scnt_ = structure_change_cnt; } @@ -275,17 +261,17 @@ void NonLinImpRep::didv() { for (i = _nt->ncell; i < n_v_; ++i) { nd = _nt->_v_node[i]; ip = _nt->_v_parent[i]->v_node_index; - m_.coeffRef(v_index_[ip], v_index_[i]) += NODEA(nd); - m_.coeffRef(v_index_[i], v_index_[ip]) += NODEB(nd); - *diag_[i] -= NODEB(nd); - *diag_[ip] -= NODEA(nd); + m_.coeffRef(ip, i) += NODEA(nd); + m_.coeffRef(i, ip) += NODEB(nd); + m_.coeffRef(i, i) -= NODEB(nd); + m_.coeffRef(ip, ip) -= NODEA(nd); } // jwC term Memb_list* mlc = _nt->tml->ml; int n = mlc->nodecount; for (i = 0; i < n; ++i) { j = mlc->nodelist[i]->v_node_index; - diag_[v_index_[j] - 1][1] += .001 * mlc->data(i, 0) * omega_; + m_.coeffRef(j, j) += std::complex(0, .001 * mlc->data(i, 0) * omega_); } // di/dv terms // because there may be several point processes of the same type @@ -326,7 +312,8 @@ void NonLinImpRep::didv() { current(i, ml, j); // conductance // add to matrix - *diag_[v_index_[nd->v_node_index] - 1] -= (x2 - NODERHS(nd)) / delta_; + m_.coeffRef(nd->v_node_index, + nd->v_node_index) -= std::complex((x2 - NODERHS(nd)) / delta_, 0); } } } @@ -366,7 +353,7 @@ void NonLinImpRep::dids() { *pv_[is] = v_[is].real(); // restore s double g = (NODERHS(nd) - v_[in].imag()) / deltavec_[is]; if (g != 0.) { - m_.coeffRef(v_index_[nd->v_node_index], is + 1) = -g; + m_.coeffRef(nd->v_node_index, is) = -g; } } // don't know if this is necessary but make sure last @@ -427,7 +414,7 @@ void NonLinImpRep::dsdv() { for (is = ieq + in * cnt, iis = 0; iis < cnt; ++iis, ++is) { double ds = (v_[is].imag() - *pvdot_[is]) / delta_; if (ds != 0.) { - m_.coeffRef(is + 1, v_index_[nd->v_node_index]) = -ds; + m_.coeffRef(is, nd->v_node_index) = -ds; } } } @@ -442,7 +429,7 @@ void NonLinImpRep::dsds() { NrnThread* nt = nrn_threads; // jw term for (i = neq_v_; i < neq_; ++i) { - *diag_[i] += std::complex(0, omega_); + m_.coeffRef(i, i) += std::complex(0, omega_); } ieq = neq_v_; for (NrnThreadMembList* tml = nt->tml; tml; tml = tml->next) { @@ -541,11 +528,11 @@ int NonLinImpRep::gapsolve() { std::vector> rx(neq_); std::vector> rx1(neq_); - std::vector> rb(v_.begin(), v_.end()); + std::vector> rb(v_); // iterate till change in x is small double tol = 1e-9; - double delta; + double delta{}; int success = 0; int iter; From 9c5e21041a6aaad0d25cee0bcdb4d0add0174208 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Mon, 4 Dec 2023 16:29:06 +0100 Subject: [PATCH 05/16] Simplifications and errors --- src/nrniv/nonlinz.cpp | 28 ++++++++++++++++++++++------ 1 file changed, 22 insertions(+), 6 deletions(-) diff --git a/src/nrniv/nonlinz.cpp b/src/nrniv/nonlinz.cpp index 3c1c2c5a70..d9f53a8f4e 100644 --- a/src/nrniv/nonlinz.cpp +++ b/src/nrniv/nonlinz.cpp @@ -11,6 +11,8 @@ #include #include +using namespace std::complex_literals; + extern void v_setup_vectors(); extern int nrndae_extra_eqn_count(); extern Symlist* hoc_built_in_symlist; @@ -141,9 +143,23 @@ void NonLinImp::compute(double omega, double deltafac, int maxiter) { rep_->dsdv(); #endif + // Now that the matrix is filled we can compressed it (mandatory for SparseLU) + rep_->m_.makeCompressed(); + rep_->lu_ = std::make_uniquem_)>>(rep_->m_); - rep_->lu_->analyzePattern(rep_->m_); - rep_->lu_->factorize(rep_->m_); + rep_->lu_->compute(rep_->m_); + auto info = rep_->lu_->info(); + switch (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 invalid", nullptr); + break; + }; rep_->iloc_ = -2; } @@ -158,7 +174,7 @@ int NonLinImp::solve(int curloc) { rep_->iloc_ = curloc; rep_->v_ = std::vector>(rep_->neq_); if (curloc >= 0) { - rep_->v_[curloc].real(1.e2 / NODEAREA(_nt->_v_node[curloc])); + rep_->v_[curloc] = 1.e2 / NODEAREA(_nt->_v_node[curloc]); } if (nrnthread_v_transfer_) { rval = rep_->gapsolve(); @@ -271,7 +287,7 @@ void NonLinImpRep::didv() { int n = mlc->nodecount; for (i = 0; i < n; ++i) { j = mlc->nodelist[i]->v_node_index; - m_.coeffRef(j, j) += std::complex(0, .001 * mlc->data(i, 0) * omega_); + m_.coeffRef(j, j) += .001 * mlc->data(i, 0) * omega_ * 1i; } // di/dv terms // because there may be several point processes of the same type @@ -313,7 +329,7 @@ void NonLinImpRep::didv() { // conductance // add to matrix m_.coeffRef(nd->v_node_index, - nd->v_node_index) -= std::complex((x2 - NODERHS(nd)) / delta_, 0); + nd->v_node_index) -= (x2 - NODERHS(nd)) / delta_; } } } @@ -429,7 +445,7 @@ void NonLinImpRep::dsds() { NrnThread* nt = nrn_threads; // jw term for (i = neq_v_; i < neq_; ++i) { - m_.coeffRef(i, i) += std::complex(0, omega_); + m_.coeffRef(i, i) += omega_ * 1i; } ieq = neq_v_; for (NrnThreadMembList* tml = nt->tml; tml; tml = tml->next) { From c27defa412ad6563531f5b038e2495169e6b662f Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Mon, 4 Dec 2023 16:41:44 +0100 Subject: [PATCH 06/16] Format --- src/nrniv/nonlinz.cpp | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/src/nrniv/nonlinz.cpp b/src/nrniv/nonlinz.cpp index d9f53a8f4e..41c492c1f4 100644 --- a/src/nrniv/nonlinz.cpp +++ b/src/nrniv/nonlinz.cpp @@ -150,15 +150,17 @@ void NonLinImp::compute(double omega, double deltafac, int maxiter) { rep_->lu_->compute(rep_->m_); auto info = rep_->lu_->info(); switch (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 invalid", 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 invalid", nullptr); + break; }; rep_->iloc_ = -2; @@ -328,8 +330,7 @@ void NonLinImpRep::didv() { current(i, ml, j); // conductance // add to matrix - m_.coeffRef(nd->v_node_index, - nd->v_node_index) -= (x2 - NODERHS(nd)) / delta_; + m_.coeffRef(nd->v_node_index, nd->v_node_index) -= (x2 - NODERHS(nd)) / delta_; } } } From 76b7e84fe28bf79f340d423886197d6c5f91ac70 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Wed, 6 Dec 2023 16:28:32 +0100 Subject: [PATCH 07/16] decltype is broken with nvhpc --- src/nrniv/nonlinz.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/nrniv/nonlinz.cpp b/src/nrniv/nonlinz.cpp index 41c492c1f4..edc689d70a 100644 --- a/src/nrniv/nonlinz.cpp +++ b/src/nrniv/nonlinz.cpp @@ -32,8 +32,8 @@ class NonLinImpRep final { void dsds(); int gapsolve(); - Eigen::SparseMatrix, Eigen::RowMajor> m_{}; - std::unique_ptr> lu_{}; + Eigen::SparseMatrix> m_{}; + std::unique_ptr>>> lu_{}; int scnt_; // structure_change int n_v_, n_ext_, n_lin_, n_ode_, neq_v_, neq_; std::vector> pv_, pvdot_; @@ -146,8 +146,8 @@ void NonLinImp::compute(double omega, double deltafac, int maxiter) { // Now that the matrix is filled we can compressed it (mandatory for SparseLU) rep_->m_.makeCompressed(); - rep_->lu_ = std::make_uniquem_)>>(rep_->m_); - rep_->lu_->compute(rep_->m_); + rep_->lu_ = std::make_unique>>>( + rep_->m_); auto info = rep_->lu_->info(); switch (info) { case Eigen::NumericalIssue: From f53ba212bdc597366ec19dc654569cca711bbf83 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Wed, 6 Dec 2023 18:10:03 +0100 Subject: [PATCH 08/16] Simplify --- src/nrniv/nonlinz.cpp | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/src/nrniv/nonlinz.cpp b/src/nrniv/nonlinz.cpp index edc689d70a..c9fec203ec 100644 --- a/src/nrniv/nonlinz.cpp +++ b/src/nrniv/nonlinz.cpp @@ -33,7 +33,7 @@ class NonLinImpRep final { int gapsolve(); Eigen::SparseMatrix> m_{}; - std::unique_ptr>>> lu_{}; + Eigen::SparseLU>> lu_{}; int scnt_; // structure_change int n_v_, n_ext_, n_lin_, n_ode_, neq_v_, neq_; std::vector> pv_, pvdot_; @@ -146,10 +146,8 @@ void NonLinImp::compute(double omega, double deltafac, int maxiter) { // Now that the matrix is filled we can compressed it (mandatory for SparseLU) rep_->m_.makeCompressed(); - rep_->lu_ = std::make_unique>>>( - rep_->m_); - auto info = rep_->lu_->info(); - switch (info) { + rep_->lu_.compute(rep_->m_); + switch (rep_->lu_.info()) { case Eigen::NumericalIssue: hoc_execerror( "NumericalIssue: The matrix is not valid following what expect Eigen SparseLu", @@ -184,7 +182,7 @@ int NonLinImp::solve(int curloc) { auto v = Eigen::Map, Eigen::Dynamic>>(rep_->v_.data(), rep_->v_.size()); - v = rep_->lu_->solve(v); + v = rep_->lu_.solve(v); } } return rval; @@ -560,7 +558,7 @@ int NonLinImpRep::gapsolve() { rb.size()); auto rx1_ = Eigen::Map, Eigen::Dynamic>>(rx1.data(), rx1.size()); - rx1_ = lu_->solve(rb_); + rx1_ = lu_.solve(rb_); } // if any change in x > tol, then do another iteration. From 624c7dde1b5a61131f8c592b8f2a3d2ad256d097 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 7 Dec 2023 09:51:49 +0100 Subject: [PATCH 09/16] Remove all RowMajor --- src/nrniv/nonlinz.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/nrniv/nonlinz.cpp b/src/nrniv/nonlinz.cpp index c9fec203ec..2e84e8fd0d 100644 --- a/src/nrniv/nonlinz.cpp +++ b/src/nrniv/nonlinz.cpp @@ -226,7 +226,7 @@ NonLinImpRep::NonLinImpRep() { if (neq_ == 0) { return; } - m_ = Eigen::SparseMatrix, Eigen::RowMajor>(neq_, neq_); + m_ = Eigen::SparseMatrix>(neq_, neq_); pv_.resize(neq_); pvdot_.resize(neq_); v_.resize(neq_); From fb2a698a5721ed581d7925f10edfad8473f49ac5 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Tue, 12 Dec 2023 17:45:14 +0100 Subject: [PATCH 10/16] Avoid a useless copy --- src/nrniv/nonlinz.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/nrniv/nonlinz.cpp b/src/nrniv/nonlinz.cpp index 2e84e8fd0d..9c645a5661 100644 --- a/src/nrniv/nonlinz.cpp +++ b/src/nrniv/nonlinz.cpp @@ -226,7 +226,7 @@ NonLinImpRep::NonLinImpRep() { if (neq_ == 0) { return; } - m_ = Eigen::SparseMatrix>(neq_, neq_); + m_.resize(neq_, neq_); pv_.resize(neq_); pvdot_.resize(neq_); v_.resize(neq_); From e06535b06e1da968696fbf88768a7d505fccd26a Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Wed, 13 Dec 2023 09:59:11 +0100 Subject: [PATCH 11/16] No more complex sparse13 after that --- cmake/NeuronFileLists.cmake | 2 -- src/sparse13/CMakeLists.txt | 8 +----- src/sparse13/cspalloc.cpp | 2 -- src/sparse13/cspbuild.cpp | 2 -- src/sparse13/cspfactor.cpp | 2 -- src/sparse13/cspmatrix.h | 2 -- src/sparse13/cspoutput.cpp | 2 -- src/sparse13/cspredef.h | 49 ------------------------------------- src/sparse13/cspsolve.cpp | 2 -- src/sparse13/csputils.cpp | 2 -- 10 files changed, 1 insertion(+), 72 deletions(-) delete mode 100644 src/sparse13/cspalloc.cpp delete mode 100644 src/sparse13/cspbuild.cpp delete mode 100644 src/sparse13/cspfactor.cpp delete mode 100644 src/sparse13/cspmatrix.h delete mode 100644 src/sparse13/cspoutput.cpp delete mode 100644 src/sparse13/cspredef.h delete mode 100644 src/sparse13/cspsolve.cpp delete mode 100644 src/sparse13/csputils.cpp diff --git a/cmake/NeuronFileLists.cmake b/cmake/NeuronFileLists.cmake index 8cd9a06553..a3e88a85f9 100644 --- a/cmake/NeuronFileLists.cmake +++ b/cmake/NeuronFileLists.cmake @@ -68,8 +68,6 @@ set(HEADER_FILES_TO_INSTALL scopmath/sparse_thread.hpp scopmath/ssimplic.hpp scopmath/ssimplic_thread.hpp - sparse13/cspmatrix.h - sparse13/cspredef.h sparse13/spconfig.h sparse13/spmatrix.h) diff --git a/src/sparse13/CMakeLists.txt b/src/sparse13/CMakeLists.txt index 69f9e1a955..d7a58f6702 100644 --- a/src/sparse13/CMakeLists.txt +++ b/src/sparse13/CMakeLists.txt @@ -5,11 +5,5 @@ add_library( spfactor.cpp spoutput.cpp spsolve.cpp - sputils.cpp - cspalloc.cpp - cspbuild.cpp - cspfactor.cpp - cspoutput.cpp - cspsolve.cpp - csputils.cpp) + sputils.cpp) set_property(TARGET sparse13 PROPERTY POSITION_INDEPENDENT_CODE ON) diff --git a/src/sparse13/cspalloc.cpp b/src/sparse13/cspalloc.cpp deleted file mode 100644 index 29b7ee44c8..0000000000 --- a/src/sparse13/cspalloc.cpp +++ /dev/null @@ -1,2 +0,0 @@ -#define cmplx_spPrefix -#include "spalloc.cpp" diff --git a/src/sparse13/cspbuild.cpp b/src/sparse13/cspbuild.cpp deleted file mode 100644 index b002c1e5c2..0000000000 --- a/src/sparse13/cspbuild.cpp +++ /dev/null @@ -1,2 +0,0 @@ -#define cmplx_spPrefix -#include "spbuild.cpp" diff --git a/src/sparse13/cspfactor.cpp b/src/sparse13/cspfactor.cpp deleted file mode 100644 index 2a9d91aaaa..0000000000 --- a/src/sparse13/cspfactor.cpp +++ /dev/null @@ -1,2 +0,0 @@ -#define cmplx_spPrefix -#include "spfactor.cpp" diff --git a/src/sparse13/cspmatrix.h b/src/sparse13/cspmatrix.h deleted file mode 100644 index 7200d07b55..0000000000 --- a/src/sparse13/cspmatrix.h +++ /dev/null @@ -1,2 +0,0 @@ -#define cmplx_spPrefix -#include "spmatrix.h" diff --git a/src/sparse13/cspoutput.cpp b/src/sparse13/cspoutput.cpp deleted file mode 100644 index 59325c189a..0000000000 --- a/src/sparse13/cspoutput.cpp +++ /dev/null @@ -1,2 +0,0 @@ -#define cmplx_spPrefix -#include "spoutput.cpp" diff --git a/src/sparse13/cspredef.h b/src/sparse13/cspredef.h deleted file mode 100644 index 2df7687583..0000000000 --- a/src/sparse13/cspredef.h +++ /dev/null @@ -1,49 +0,0 @@ -/* mostly generated from -cat temp | /usr/bin/tr -cs "[:alpha:]" "[\n*]" | sort | uniq | grep '^sp' > sp_redef.h -where temp is the last part of spmatrix.h -*/ -#define spClear cmplx_spClear -#define spCondition cmplx_spCondition -#define spCreate cmplx_spCreate -#define spDeleteRowAndCol cmplx_spDeleteRowAndCol -#define spDestroy cmplx_spDestroy -#define spDeterminant cmplx_spDeterminant -#define spElementCount cmplx_spElementCount -#define spError cmplx_spError -#define spFactor cmplx_spFactor -#define spFileMatrix cmplx_spFileMatrix -#define spFileStats cmplx_spFileStats -#define spFileVector cmplx_spFileVector -#define spFillinCount cmplx_spFillinCount -#define spGetAdmittance cmplx_spGetAdmittance -#define spGetElement cmplx_spGetElement -#define spGetInitInfo cmplx_spGetInitInfo -#define spGetOnes cmplx_spGetOnes -#define spGetQuad cmplx_spGetQuad -#define spGetSize cmplx_spGetSize -#define spInitialize cmplx_spInitialize -#define spInstallInitInfo cmplx_spInstallInitInfo -#define spLargestElement cmplx_spLargestElement -#define spMNA_Preorder cmplx_spMNA_Preorder -#define spMultTransposed cmplx_spMultTransposed -#define spMultiply cmplx_spMultiply -#define spNorm cmplx_spNorm -#define spOrderAndFactor cmplx_spOrderAndFactor -#define spPartition cmplx_spPartition -#define spPrint cmplx_spPrint -#define spPseudoCondition cmplx_spPseudoCondition -#define spRoundoff cmplx_spRoundoff -#define spScale cmplx_spScale -#define spSetComplex cmplx_spSetComplex -#define spSetReal cmplx_spSetReal -#define spSolve cmplx_spSolve -#define spSolveTransposed cmplx_spSolveTransposed -#define spStripFills cmplx_spStripFills -#define spWhereSingular cmplx_spWhereSingular -#define spcGetFillin cmplx_spcGetFillin -#define spcGetElement cmplx_spcGetElement -#define spcLinkRows cmplx_spcLinkRows -#define spcFindElementInCol cmplx_spcFindElementInCol -#define spcCreateElement cmplx_spcCreateElement -#define spcRowExchange cmplx_spcRowExchange -#define spcColExchange cmplx_spcColExchange diff --git a/src/sparse13/cspsolve.cpp b/src/sparse13/cspsolve.cpp deleted file mode 100644 index f1aeba7d67..0000000000 --- a/src/sparse13/cspsolve.cpp +++ /dev/null @@ -1,2 +0,0 @@ -#define cmplx_spPrefix -#include "spsolve.cpp" diff --git a/src/sparse13/csputils.cpp b/src/sparse13/csputils.cpp deleted file mode 100644 index 0f0363b8ae..0000000000 --- a/src/sparse13/csputils.cpp +++ /dev/null @@ -1,2 +0,0 @@ -#define cmplx_spPrefix -#include "sputils.cpp" From b895c2fdd51a72aaefe701bddb67ad443acbbe43 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Wed, 13 Dec 2023 10:01:00 +0100 Subject: [PATCH 12/16] Add Success for silence warnings --- src/nrniv/nonlinz.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/nrniv/nonlinz.cpp b/src/nrniv/nonlinz.cpp index 9c645a5661..10989d097a 100644 --- a/src/nrniv/nonlinz.cpp +++ b/src/nrniv/nonlinz.cpp @@ -148,6 +148,9 @@ void NonLinImp::compute(double omega, double deltafac, int maxiter) { rep_->lu_.compute(rep_->m_); switch (rep_->lu_.info()) { + case Eigen::Success: + // Everything fine + break; case Eigen::NumericalIssue: hoc_execerror( "NumericalIssue: The matrix is not valid following what expect Eigen SparseLu", @@ -159,7 +162,7 @@ void NonLinImp::compute(double omega, double deltafac, int maxiter) { case Eigen::InvalidInput: hoc_execerror("InvalidInput: the inputs are invalid", nullptr); break; - }; + } rep_->iloc_ = -2; } From a8f33b3fd1ff5c802153821d9ca55416df0b6f3f Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Wed, 13 Dec 2023 10:39:27 +0100 Subject: [PATCH 13/16] Format + sonar --- src/nrniv/nonlinz.cpp | 7 +++---- src/sparse13/CMakeLists.txt | 10 ++-------- 2 files changed, 5 insertions(+), 12 deletions(-) diff --git a/src/nrniv/nonlinz.cpp b/src/nrniv/nonlinz.cpp index 10989d097a..747cdbcb82 100644 --- a/src/nrniv/nonlinz.cpp +++ b/src/nrniv/nonlinz.cpp @@ -195,7 +195,6 @@ int NonLinImp::solve(int curloc) { // mapping is already done there. NonLinImpRep::NonLinImpRep() { - int i, cnt; NrnThread* _nt = nrn_threads; Symbol* vsym = hoc_table_lookup("v", hoc_built_in_symlist); @@ -217,10 +216,10 @@ NonLinImpRep::NonLinImpRep() { n_ode_ = 0; for (NrnThreadMembList* tml = _nt->tml; tml; tml = tml->next) { Memb_list* ml = tml->ml; - i = tml->index; + int i = tml->index; nrn_ode_count_t s = memb_func[i].ode_count; if (s) { - cnt = (*s)(i); + int cnt = (*s)(i); n_ode_ += cnt * ml->nodecount; } } @@ -235,7 +234,7 @@ NonLinImpRep::NonLinImpRep() { v_.resize(neq_); deltavec_.resize(neq_); - for (i = 0; i < n_v_; ++i) { + for (int i = 0; i < n_v_; ++i) { // utilize nd->eqn_index in case of use_sparse13 later Node* nd = _nt->_v_node[i]; pv_[i] = nd->v_handle(); diff --git a/src/sparse13/CMakeLists.txt b/src/sparse13/CMakeLists.txt index d7a58f6702..b2adc31dc2 100644 --- a/src/sparse13/CMakeLists.txt +++ b/src/sparse13/CMakeLists.txt @@ -1,9 +1,3 @@ -add_library( - sparse13 STATIC - spalloc.cpp - spbuild.cpp - spfactor.cpp - spoutput.cpp - spsolve.cpp - sputils.cpp) +add_library(sparse13 STATIC spalloc.cpp spbuild.cpp spfactor.cpp spoutput.cpp spsolve.cpp + sputils.cpp) set_property(TARGET sparse13 PROPERTY POSITION_INDEPENDENT_CODE ON) From 27af30aca901e96ce99533733d0f522face3f2b5 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Mon, 18 Dec 2023 16:10:27 +0100 Subject: [PATCH 14/16] Follow review --- src/nrniv/nonlinz.cpp | 18 +++++++++++++----- src/nrniv/nonlinz.h | 3 +++ 2 files changed, 16 insertions(+), 5 deletions(-) diff --git a/src/nrniv/nonlinz.cpp b/src/nrniv/nonlinz.cpp index 747cdbcb82..17c3f0c2fb 100644 --- a/src/nrniv/nonlinz.cpp +++ b/src/nrniv/nonlinz.cpp @@ -32,7 +32,9 @@ class NonLinImpRep final { void dsds(); int gapsolve(); + // Matrix containing the non linear system to solve. Eigen::SparseMatrix> m_{}; + // The solver of the matrix using the LU decomposition method. Eigen::SparseLU>> lu_{}; int scnt_; // structure_change int n_v_, n_ext_, n_lin_, n_ode_, neq_v_, neq_; @@ -143,7 +145,7 @@ void NonLinImp::compute(double omega, double deltafac, int maxiter) { rep_->dsdv(); #endif - // Now that the matrix is filled we can compressed it (mandatory for SparseLU) + // Now that the matrix is filled we can compress it (mandatory for SparseLU) rep_->m_.makeCompressed(); rep_->lu_.compute(rep_->m_); @@ -153,14 +155,20 @@ void NonLinImp::compute(double omega, double deltafac, int maxiter) { break; case Eigen::NumericalIssue: hoc_execerror( - "NumericalIssue: The matrix is not valid following what expect Eigen SparseLu", - nullptr); + "Eigen Sparse LU factorization failed with Eigen::NumericalIssue, please check the " + "input matrix:", + rep_->lu_.lastErrorMessage().c_str()); break; case Eigen::NoConvergence: - hoc_execerror("NoConvergence: The matrix did not converge", nullptr); + hoc_execerror( + "Eigen Sparse LU factorization reports Eigen::NonConvergence after calling compute():", + rep_->lu_.lastErrorMessage().c_str()); break; case Eigen::InvalidInput: - hoc_execerror("InvalidInput: the inputs are invalid", nullptr); + hoc_execerror( + "Eigen Sparse LU factorization failed with Eigen::InvalidInput, the input matrix seems " + "invalid:", + rep_->lu_.lastErrorMessage().c_str()); break; } diff --git a/src/nrniv/nonlinz.h b/src/nrniv/nonlinz.h index ca7a568f37..96e5125aa2 100644 --- a/src/nrniv/nonlinz.h +++ b/src/nrniv/nonlinz.h @@ -3,9 +3,12 @@ class NonLinImpRep; +// A solve for non linear equation of complex numbers. +// Matrix should be squared. class NonLinImp final { public: ~NonLinImp(); + // Prepare the matrix before solving it. void compute(double omega, double deltafac, int maxiter); double transfer_amp(int curloc, int vloc); // v_node[arg] is the node double transfer_phase(int curloc, int vloc); From 532b8c64d1ec1a55219f8c2ed96cba52faacd9f5 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Tue, 19 Dec 2023 11:26:01 +0100 Subject: [PATCH 15/16] Remove final --- src/nrniv/nonlinz.cpp | 6 +++++- src/nrniv/nonlinz.h | 8 +++++--- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/src/nrniv/nonlinz.cpp b/src/nrniv/nonlinz.cpp index 17c3f0c2fb..17f6ad4271 100644 --- a/src/nrniv/nonlinz.cpp +++ b/src/nrniv/nonlinz.cpp @@ -22,14 +22,17 @@ extern void pargap_jacobi_rhs(std::vector>&, const std::vector>&); extern void pargap_jacobi_setup(int mode); -class NonLinImpRep final { +class NonLinImpRep { public: NonLinImpRep(); void delta(double); + + // Functions to fill the matrix void didv(); void dids(); void dsdv(); void dsds(); + int gapsolve(); // Matrix containing the non linear system to solve. @@ -148,6 +151,7 @@ void NonLinImp::compute(double omega, double deltafac, int maxiter) { // Now that the matrix is filled we can compress it (mandatory for SparseLU) rep_->m_.makeCompressed(); + // Factorize the matrix so this is ready to solve rep_->lu_.compute(rep_->m_); switch (rep_->lu_.info()) { case Eigen::Success: diff --git a/src/nrniv/nonlinz.h b/src/nrniv/nonlinz.h index 96e5125aa2..b477bf783a 100644 --- a/src/nrniv/nonlinz.h +++ b/src/nrniv/nonlinz.h @@ -3,18 +3,20 @@ class NonLinImpRep; -// A solve for non linear equation of complex numbers. +// A solver for non linear equation of complex numbers. // Matrix should be squared. -class NonLinImp final { +class NonLinImp { public: ~NonLinImp(); // Prepare the matrix before solving it. void compute(double omega, double deltafac, int maxiter); - double transfer_amp(int curloc, int vloc); // v_node[arg] is the node + + double transfer_amp(int curloc, int vloc); double transfer_phase(int curloc, int vloc); double input_amp(int curloc); double input_phase(int curloc); double ratio_amp(int clmploc, int vloc); + int solve(int curloc); private: From 73123806bd09752215f58647e686dc16d5adf9bc Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Tue, 19 Dec 2023 14:40:42 +0100 Subject: [PATCH 16/16] More doc --- src/nrniv/nonlinz.cpp | 40 +++++++++++++++++++++------------------- 1 file changed, 21 insertions(+), 19 deletions(-) diff --git a/src/nrniv/nonlinz.cpp b/src/nrniv/nonlinz.cpp index 17f6ad4271..f0a24eac8d 100644 --- a/src/nrniv/nonlinz.cpp +++ b/src/nrniv/nonlinz.cpp @@ -534,6 +534,8 @@ void NonLinImpRep::ode(int im, Memb_list* ml) { // assume there is in fact an o memb_func[im].ode_spec(nrn_ensure_model_data_are_sorted(), nrn_threads, ml, im); } +// This function compute a solution of a converging system by iteration. +// The value returned is the number of iterations to reach a precision of "tol" (1e-9). int NonLinImpRep::gapsolve() { // On entry, v_ contains the complex b for A*x = b. // On return v_ contains complex solution, x. @@ -553,11 +555,11 @@ int NonLinImpRep::gapsolve() { } #endif - pargap_jacobi_setup(0); + pargap_jacobi_setup(0); // 0 means 'setup' - std::vector> rx(neq_); - std::vector> rx1(neq_); - std::vector> rb(v_); + std::vector> x_old(neq_); + std::vector> x(neq_); + std::vector> b(v_); // iterate till change in x is small double tol = 1e-9; @@ -568,25 +570,25 @@ int NonLinImpRep::gapsolve() { for (iter = 1; iter <= maxiter_; ++iter) { if (neq_) { - auto rb_ = Eigen::Map, Eigen::Dynamic>>(rb.data(), - rb.size()); - auto rx1_ = Eigen::Map, Eigen::Dynamic>>(rx1.data(), - rx1.size()); - rx1_ = lu_.solve(rb_); + auto b_ = Eigen::Map, Eigen::Dynamic>>(b.data(), + b.size()); + auto x_ = Eigen::Map, Eigen::Dynamic>>(x.data(), + x.size()); + x_ = lu_.solve(b_); } // if any change in x > tol, then do another iteration. success = 1; delta = 0.0; + // Do the substraction of the previous result (x_old) and current result (x). + // If all differences are < tol stop the loop, otherwise continue to iterate for (int i = 0; i < neq_; ++i) { - auto diff = rx1[i] - rx[i]; + auto diff = x[i] - x_old[i]; double err = std::abs(diff.real()) + std::abs(diff.imag()); if (err > tol) { success = 0; } - if (delta < err) { - delta = err; - } + delta = std::max(err, delta); } #if NRNMPI if (nrnmpi_numprocs > 1) { @@ -594,17 +596,17 @@ int NonLinImpRep::gapsolve() { } #endif if (success) { - v_ = rx1; + v_ = x; break; } // setup for next iteration - rx = rx1; - rb = v_; - pargap_jacobi_rhs(rb, rx); + x_old = x; + b = v_; + pargap_jacobi_rhs(b, x_old); } - pargap_jacobi_setup(1); // tear down + pargap_jacobi_setup(1); // 1 means 'tear down' if (!success) { char buf[256]; @@ -614,7 +616,7 @@ int NonLinImpRep::gapsolve() { maxiter_, delta, tol); - execerror(buf, 0); + hoc_execerror(buf, nullptr); } return iter; }