Skip to content

Commit

Permalink
To complex
Browse files Browse the repository at this point in the history
  • Loading branch information
Nicolas Cornu committed Dec 1, 2023
1 parent f71b138 commit 817edb2
Showing 1 changed file with 35 additions and 95 deletions.
130 changes: 35 additions & 95 deletions src/nrniv/nonlinz.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,16 +6,15 @@
#include "nrniv_mf.h"
#include "nrnoc2iv.h"
#include "nrnmpi.h"
#include "cspmatrix.h"
#include "membfunc.h"

#include <Eigen/Sparse>
#include <complex>

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_setup(int mode);
Expand All @@ -36,7 +35,7 @@ class NonLinImpRep {
int n_v_, n_ext_, n_lin_, n_ode_, neq_v_, neq_;
std::vector<neuron::container::data_handle<double>> pv_, pvdot_;
int* v_index_;
double* rv_;
std::vector<std::complex<double>> rv_;
double* jv_;
std::vector<std::complex<double>*> diag_;
double* deltavec_; // just like cvode.atol*cvode.atolscale for ode's
Expand Down Expand Up @@ -66,9 +65,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_->rv_[vloc]);
}
double NonLinImp::input_amp(int curloc) {
if (nrnmpi_numprocs > 1 && nrnthread_v_transfer_) {
Expand All @@ -80,9 +77,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_->rv_[curloc]);
}
double NonLinImp::transfer_phase(int curloc, int vloc) {
if (nrnmpi_numprocs > 1 && nrnthread_v_transfer_ && curloc != rep_->iloc_) {
Expand All @@ -92,9 +87,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_->rv_[vloc]);
}
double NonLinImp::input_phase(int curloc) {
if (nrnmpi_numprocs > 1 && nrnthread_v_transfer_) {
Expand All @@ -106,9 +99,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_->rv_[curloc]);
}
double NonLinImp::ratio_amp(int clmploc, int vloc) {
if (nrnmpi_numprocs > 1 && nrnthread_v_transfer_) {
Expand All @@ -120,15 +111,7 @@ 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_->rv_[vloc] * std::conj(rep_->rv_[clmploc]) / std::norm(rep_->rv_[clmploc]));
}
void NonLinImp::compute(double omega, double deltafac, int maxiter) {
v_setup_vectors();
Expand Down Expand Up @@ -175,29 +158,19 @@ 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_->rv_[i] = std::complex<double>(0, 0);
}
if (curloc >= 0) {
rep_->rv_[curloc] = 1.e2 / NODEAREA(_nt->_v_node[curloc]);
}
if (nrnthread_v_transfer_) {
rval = rep_->gapsolve();
} else {
std::vector<std::complex<double>> 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::SparseLU<decltype(rep_->m_)>(rep_->m_);
lu.analyzePattern(rep_->m_);
lu.factorize(rep_->m_);
auto rjv_ = Eigen::Map<Eigen::Vector<std::complex<double>, 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::Vector<std::complex<double>, Eigen::Dynamic>>(rep_->rv_.data(), rep_->rv_.size());
rv = lu.solve(rv);
}
}
return rval;
Expand Down Expand Up @@ -248,10 +221,7 @@ NonLinImpRep::NonLinImpRep() {
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;
rv_.resize(neq_);
diag_.resize(neq_);
deltavec_ = new double[neq_];

Expand All @@ -273,8 +243,6 @@ NonLinImpRep::NonLinImpRep() {

NonLinImpRep::~NonLinImpRep() {
delete[] v_index_;
delete[](rv_ - 1);
delete[](jv_ - 1);
delete[] deltavec_;
}

Expand Down Expand Up @@ -381,27 +349,25 @@ 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
NODERHS(nd) = 0;
// compute rhs
current(i, ml, in);
// save rhs
x2[in] = NODERHS(nd);
rv_[in].imag() = NODERHS(nd);

Check failure on line 359 in src/nrniv/nonlinz.cpp

View workflow job for this annotation

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

expression is not assignable
// each state incremented separately and restored
for (iis = 0; iis < cnt; ++iis) {
is = ieq + in * cnt + iis;
// save s
x1[is] = *pv_[is];
rv_[is].real() = *pv_[is];

Check failure on line 364 in src/nrniv/nonlinz.cpp

View workflow job for this annotation

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

expression is not assignable
// 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] = rv_[is].real(); // restore s
double g = (NODERHS(nd) - rv_[in].imag()) / deltavec_[is];
if (g != 0.) {
m_.coeffRef(v_index_[nd->v_node_index], is + 1) = -g;
}
Expand All @@ -428,22 +394,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);
rv_[in].real() = NODEV(nd);

Check failure on line 403 in src/nrniv/nonlinz.cpp

View workflow job for this annotation

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

expression is not assignable
}
// 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 (rv_[in].real() == v) {
nd->v() = v + delta_;
}
}
Expand All @@ -453,18 +417,18 @@ 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];
rv_[is].imag() = *pvdot_[is];

Check failure on line 420 in src/nrniv/nonlinz.cpp

View workflow job for this annotation

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

expression is not assignable
*pvdot_[is] = 0;
}
nd->v() = x1[in];
nd->v() = rv_[in].real();
}
// compute the rhs(v)
ode(i, ml);
// fill the ds/dv elements
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 = (rv_[is].imag() - *pvdot_[is]) / delta_;
if (ds != 0.) {
m_.coeffRef(is + 1, v_index_[nd->v_node_index]) = -ds;
}
Expand All @@ -491,21 +455,19 @@ 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];
rv_[is].real() = *pv_[is];

Check failure on line 462 in src/nrniv/nonlinz.cpp

View workflow job for this annotation

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

expression is not assignable
}
}
// compute rhs. this is the rhs(s)
ode(i, ml);
// save rhs
for (in = 0; in < ml->nodecount; ++in) {
for (is = ieq + in * cnt, iis = 0; iis < cnt; ++iis, ++is) {
x2[is] = *pvdot_[is];
rv_[is].imag() = *pvdot_[is];

Check failure on line 470 in src/nrniv/nonlinz.cpp

View workflow job for this annotation

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

expression is not assignable
}
}
// iterate over the states
Expand All @@ -525,11 +487,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] - rv_[is].imag()) / deltavec_[is];
if (ds != 0.) {
m_.coeffRef(is, ks) = -ds;
}
*pv_[ks] = x1[ks];
*pv_[ks] = rv_[ks].real();
}
}
// perhaps not necessary but ensures the last computation is with
Expand Down Expand Up @@ -580,22 +542,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<std::complex<double>> rx(neq_);
std::vector<std::complex<double>> rx1(neq_);
std::vector<std::complex<double>> rb(rv_.begin(), rv_.end());

// iterate till change in x is small
double tol = 1e-9;
Expand All @@ -614,19 +563,17 @@ int NonLinImpRep::gapsolve() {
auto lu = Eigen::SparseLU<decltype(m_)>(m_);
lu.analyzePattern(m_);
lu.factorize(m_);
auto rjv_ = Eigen::Map<Eigen::Vector<std::complex<double>, 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::Vector<std::complex<double>, Eigen::Dynamic>>(rb.data(), rb.size());
auto rx1_ = Eigen::Map<Eigen::Vector<std::complex<double>, 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;
}
Expand All @@ -640,20 +587,13 @@ int NonLinImpRep::gapsolve() {
}
#endif
if (success) {
for (int i = 0; i < neq_; ++i) {
rv_[i] = rx1[i];
jv_[i] = jx1[i];
}
rv_ = 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 = rv_;
pargap_jacobi_rhs(rb, rx);

Check failure on line 597 in src/nrniv/nonlinz.cpp

View workflow job for this annotation

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

no matching function for call to 'pargap_jacobi_rhs'
pargap_jacobi_rhs(jb, jx);

Check failure on line 598 in src/nrniv/nonlinz.cpp

View workflow job for this annotation

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

use of undeclared identifier 'jb'

Check failure on line 598 in src/nrniv/nonlinz.cpp

View workflow job for this annotation

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

use of undeclared identifier 'jx'
}
Expand Down

0 comments on commit 817edb2

Please sign in to comment.