Simplify matrix stuff and fix a weird bug where rhs = mat * solution

src/nrniv/kschan.cpp

@@ -9,7 +9,7 @@
 #include "ocnotify.h"
 #include "parse.hpp"
 #include "nrniv_mf.h"
-#include "ocmatrix.h"
+#include <Eigen/Sparse>
 
 #define NSingleIndex 0
 #if defined(__MWERKS__) && !defined(_MSC_VER)
@@ -879,9 +879,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 +1272,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 +2104,34 @@ void KSChan::setupmat() {
     if (!nksstate_) {
         return;
     }
-    mat_ = std::make_unique<OcFullMatrix>(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<double>(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,39 +2141,52 @@ 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<double> 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<Eigen::Vector<double, Eigen::Dynamic>>(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];
     }
 }
 
 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<double> s(nksstate_);
+    std::vector<double> 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<Eigen::Vector<double, Eigen::Dynamic>>(s.data(), s.size());
+    auto ds_ = Eigen::Map<Eigen::Vector<double, Eigen::Dynamic>>(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];
     }
 }

src/nrniv/kschan.h

@@ -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<OcFullMatrix> mat_;
+    Eigen::SparseMatrix<double> mat_{};
+    Eigen::SparseLU<Eigen::SparseMatrix<double>> lu_{};
     double** elms_;
     double** diag_;
     int dsize_;       // size of prop->dparam