Skip to content

Commit

Permalink
Simplify matrix stuff and fix a weird bug where rhs = mat * solution
Browse files Browse the repository at this point in the history
  • Loading branch information
@alkino
alkino committed Dec 12, 2023
1 parent 72f7558 commit 8eab595
Show file tree
Hide file tree
Showing 2 changed files with 51 additions and 53 deletions.
99 changes: 48 additions & 51 deletions src/nrniv/kschan.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -2111,42 +2103,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");
}
Expand All @@ -2156,39 +2140,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()) {

Check warning on line 2155 in src/nrniv/kschan.cpp

View workflow job for this annotation

GitHub Actions / ubuntu-22.04 - cmake (-DNRN_ENABLE_CORENEURON=ON -DNRN_ENABLE_INTERVIEWS=OFF -DNMODL_SANITIZERS=undefinedundefined) 

enumeration value 'Success' not handled in switch [-Wswitch]
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];
}
}
Expand Down
5 changes: 3 additions & 2 deletions src/nrniv/kschan.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
#include "ivocvect.h"
#include "nrnunits.h"

#include "ocmatrix.h"
#include <Eigen/Sparse>

// extern double dt;
extern double celsius;
Expand Down Expand Up @@ -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
Expand Down

0 comments on commit 8eab595

Please sign in to comment.