From 6cba0bd7d28f2a7bc2e6d0ab546478d32e20a7c7 Mon Sep 17 00:00:00 2001
From: Ioannis Magkanaris <iomagkanaris@gmail.com>
Date: Thu, 12 Sep 2019 16:20:06 +0200
Subject: [PATCH 01/14] Initial commit

---
 coreneuron/nrniv/main1.cpp         | 11 +++++
 coreneuron/nrniv/nrnoptarg.cpp     |  1 +
 coreneuron/nrnoc/capac.cpp         |  7 +++
 coreneuron/nrnoc/fadvance_core.cpp | 29 +++++++++++++
 coreneuron/nrnoc/finitialize.cpp   |  7 +++
 coreneuron/nrnoc/multicore.cpp     | 69 ++++++++++++++++++++++++++++++
 coreneuron/nrnoc/multicore.h       |  8 ++++
 coreneuron/nrnoc/treeset_core.cpp  | 28 +++++++++++-
 8 files changed, 159 insertions(+), 1 deletion(-)

diff --git a/coreneuron/nrniv/main1.cpp b/coreneuron/nrniv/main1.cpp
index a7ad6c0f2..cd37002aa 100644
--- a/coreneuron/nrniv/main1.cpp
+++ b/coreneuron/nrniv/main1.cpp
@@ -416,6 +416,13 @@ static void trajectory_return() {
     }
 }
 
+bool nrn_use_fast_imem;
+extern void nrn_fast_imem_alloc();
+static void use_fast_imem() {
+    nrn_use_fast_imem = true;
+    nrn_fast_imem_alloc();
+}
+
 }  // namespace coreneuron
 
 /// The following high-level functions are marked as "extern C"
@@ -501,6 +508,10 @@ extern "C" int run_solve_core(int argc, char** argv) {
             (*nrn2core_part2_clean_)();
         }
 
+        if (nrnopt_get_flag("--fast_imem")) {
+            use_fast_imem();
+        }
+
         // TODO : if some ranks are empty then restore will go in deadlock
         // phase (as some ranks won't have restored anything and hence return
         // false in checkpoint_initialize
diff --git a/coreneuron/nrniv/nrnoptarg.cpp b/coreneuron/nrniv/nrnoptarg.cpp
index 0099e476b..93f9ecb8e 100644
--- a/coreneuron/nrniv/nrnoptarg.cpp
+++ b/coreneuron/nrniv/nrnoptarg.cpp
@@ -98,6 +98,7 @@ static param_flag param_flag_args[] = {
     {"--multisend", "Use Multisend spike exchange instead of Allgather."},
     {"--binqueue", "Use bin queue."},
     {"--skip-mpi-finalize", "Do not call mpi finalize."},
+    {"--fast_imem", "Use fast_i_membrane."},
     {NULL, NULL}};
 
 static param_str param_str_args[] = {
diff --git a/coreneuron/nrnoc/capac.cpp b/coreneuron/nrnoc/capac.cpp
index e7b07e69f..7fafee179 100644
--- a/coreneuron/nrnoc/capac.cpp
+++ b/coreneuron/nrnoc/capac.cpp
@@ -206,6 +206,13 @@ void nrn_div_capacity(NrnThread* _nt, Memb_list* ml, int type) {
         VEC_RHS(ni[_iml]) /= 1.e-3 * cm;
         // fprintf(stderr, "== nrn_div_cap: RHS[%d]=%.12f\n", ni[_iml], VEC_RHS(ni[_iml])) ;
     }
+    /*if (_nt->_nrn_fast_imem) {
+        Node **vnode = ml->nodelist;
+        double* p = _nt->_nrn_fast_imem->_nrn_sav_rhs;
+        for (int i=0; i <_cntml_actual; ++i) {
+            p[vnode[i]->v_node_index] += i_cap;
+        }
+    }*/
 }
 
 void nrn_mul_capacity(NrnThread* _nt, Memb_list* ml, int type) {
diff --git a/coreneuron/nrnoc/fadvance_core.cpp b/coreneuron/nrnoc/fadvance_core.cpp
index 25b703af1..a3f28d61d 100644
--- a/coreneuron/nrnoc/fadvance_core.cpp
+++ b/coreneuron/nrnoc/fadvance_core.cpp
@@ -39,6 +39,8 @@ THE POSSIBILITY OF SUCH DAMAGE.
 
 namespace coreneuron {
 
+extern bool nrn_use_fast_imem;
+
 static void* nrn_fixed_step_thread(NrnThread*);
 static void* nrn_fixed_step_group_thread(NrnThread*);
 
@@ -153,6 +155,32 @@ static void* nrn_fixed_step_group_thread(NrnThread* nth) {
     return (void*)0;
 }
 
+void nrn_calc_fast_imem(NrnThread* _nt) {
+    int i;
+    int i1 = 0;
+    int i3 = _nt->end;
+
+    double* vec_rhs = &(VEC_RHS(0));
+    double* vec_area = &(VEC_AREA(0));
+
+    double* pd = _nt->_nrn_fast_imem->_nrn_sav_d;
+    double* prhs = _nt->_nrn_fast_imem->_nrn_sav_rhs;
+    FILE *fp_rhs, *fp_d;
+    fp_rhs = fopen("rhs.CORENEURON", "a");
+    fp_d = fopen("d.CORENEURON", "a");
+    fprintf(fp_rhs, "%.8e time\n", _nt->_t);
+    fprintf(fp_rhs, "%.8e time\n", _nt->_t);
+    for (i = i1; i < i3 ; ++i) {
+        prhs[i] = (pd[i]*vec_rhs[i] + prhs[i])*vec_area[i]*0.01;
+        fprintf(fp_rhs, "%.8e, ", prhs[i]);
+    }
+    fprintf(fp_rhs, "\n");
+    for (i = i1; i < i3 ; ++i) {
+        fprintf(fp_d, "%.8e, ", pd[i]);
+    }
+    fprintf(fp_d, "\n");
+}
+
 void update(NrnThread* _nt) {
     int i, i1, i2;
     i1 = 0;
@@ -190,6 +218,7 @@ void update(NrnThread* _nt) {
         assert(_nt->tml->index == CAP);
         nrn_cur_capacitance(_nt, _nt->tml->ml, _nt->tml->index);
     }
+    if (nrn_use_fast_imem) { nrn_calc_fast_imem(_nt); }
 }
 
 void nonvint(NrnThread* _nt) {
diff --git a/coreneuron/nrnoc/finitialize.cpp b/coreneuron/nrnoc/finitialize.cpp
index 2ef266d03..0e2fdfe04 100644
--- a/coreneuron/nrnoc/finitialize.cpp
+++ b/coreneuron/nrnoc/finitialize.cpp
@@ -32,6 +32,10 @@ THE POSSIBILITY OF SUCH DAMAGE.
 #include "coreneuron/nrniv/profiler_interface.h"
 
 namespace coreneuron {
+
+extern bool nrn_use_fast_imem;
+extern void nrn_calc_fast_imem(NrnThread* _nt);
+
 void nrn_finitialize(int setv, double v) {
     int i;
     NrnThread* _nt;
@@ -98,6 +102,9 @@ void nrn_finitialize(int setv, double v) {
     }
     for (i = 0; i < nrn_nthread; ++i) {
         setup_tree_matrix_minimal(nrn_threads + i);
+        if (nrn_use_fast_imem) {
+            nrn_calc_fast_imem(nrn_threads + i);
+        }
     }
     for (i = 0; i < nrn_nthread; ++i) {
         nrn_ba(nrn_threads + i, BEFORE_STEP);
diff --git a/coreneuron/nrnoc/multicore.cpp b/coreneuron/nrnoc/multicore.cpp
index a4ebe6b16..b98513998 100644
--- a/coreneuron/nrnoc/multicore.cpp
+++ b/coreneuron/nrnoc/multicore.cpp
@@ -137,6 +137,8 @@ void nrn_threads_create(int n) {
                 nt->_watch_types = NULL;
                 nt->mapping = NULL;
                 nt->trajec_requests = NULL;
+
+                nt->_nrn_fast_imem = 0;
             }
         }
         v_structure_change = 1;
@@ -145,6 +147,73 @@ void nrn_threads_create(int n) {
     /*printf("nrn_threads_create %d %d\n", nrn_nthread, nrn_thread_parallel_);*/
 }
 
+/*
+Avoid invalidating pointers to i_membrane_ unless the number of compartments
+in a thread has changed.
+*/
+extern bool nrn_use_fast_imem;
+static int fast_imem_nthread_ = 0;
+static int* fast_imem_size_ = NULL;
+static _nrn_Fast_Imem* fast_imem_;
+
+static void fast_imem_free() {
+    int i;
+    for (i = 0; i < nrn_nthread; ++i) {
+        nrn_threads[i]._nrn_fast_imem = NULL;
+    }
+    for (i = 0; i < fast_imem_nthread_; ++i) {
+        if (fast_imem_size_[i] > 0) {
+            free(fast_imem_[i]._nrn_sav_rhs);
+            free(fast_imem_[i]._nrn_sav_d);
+        }
+    }
+    if (fast_imem_nthread_) {
+        free(fast_imem_size_);
+        free(fast_imem_);
+        fast_imem_nthread_ = 0;
+        fast_imem_size_ = NULL;
+        fast_imem_ = NULL;
+    }
+}
+
+static void fast_imem_alloc() {
+    int i;
+    if (fast_imem_nthread_ != nrn_nthread) {
+        fast_imem_free();
+        fast_imem_nthread_ = nrn_nthread;
+        fast_imem_size_ = (int*)ecalloc(nrn_nthread, sizeof(int));
+        fast_imem_ = (_nrn_Fast_Imem*)ecalloc(nrn_nthread, sizeof(_nrn_Fast_Imem));
+    }
+    for (i=0; i < nrn_nthread; ++i) {
+        NrnThread* nt = nrn_threads + i;
+        int n = nt->end;
+        _nrn_Fast_Imem* fi = fast_imem_ + i;
+        if (n != fast_imem_size_[i]) {
+            if (fast_imem_size_[i] > 0) {
+                free(fi->_nrn_sav_rhs);
+                free(fi->_nrn_sav_d);
+            }
+            if (n > 0) {
+                fi->_nrn_sav_rhs = (double*)emalloc_align(n * sizeof(double));
+                fi->_nrn_sav_d = (double*)emalloc_align(n * sizeof(double));
+            }
+            fast_imem_size_[i] = n;
+        }
+    }
+}
+
+void nrn_fast_imem_alloc() {
+    if (nrn_use_fast_imem) {
+        int i;
+        fast_imem_alloc();
+        for (i=0; i < nrn_nthread; ++i) {
+            nrn_threads[i]._nrn_fast_imem = fast_imem_ + i;
+        }
+    }else{
+        fast_imem_free();
+    }
+}
+
 void nrn_threads_free() {
     if (nrn_nthread) {
         free_memory((void*)nrn_threads);
diff --git a/coreneuron/nrnoc/multicore.h b/coreneuron/nrnoc/multicore.h
index 7b9112080..7272fb5bb 100644
--- a/coreneuron/nrnoc/multicore.h
+++ b/coreneuron/nrnoc/multicore.h
@@ -59,6 +59,11 @@ struct NrnThreadBAList {
     NrnThreadBAList* next;
 };
 
+typedef struct _nrn_Fast_Imem {
+    double* _nrn_sav_rhs;
+    double* _nrn_sav_d;
+} _nrn_Fast_Imem;
+
 struct TrajectoryRequests {
     void** vpr;       /* PlayRecord Objects known by NEURON */
     double** scatter; /* if bsize == 0, each time step */
@@ -116,6 +121,9 @@ struct NrnThread {
                              compartment */
     double* _shadow_d;    /* Not pointer into _data. Avoid race for multiple POINT_PROCESS in same
                              compartment */
+
+    _nrn_Fast_Imem* _nrn_fast_imem;
+
     int* _v_parent_index;
     int* _permute;
     char* _sp13mat;              /* handle to general sparse matrix */
diff --git a/coreneuron/nrnoc/treeset_core.cpp b/coreneuron/nrnoc/treeset_core.cpp
index 52ebc4eb4..fe27a0304 100644
--- a/coreneuron/nrnoc/treeset_core.cpp
+++ b/coreneuron/nrnoc/treeset_core.cpp
@@ -67,6 +67,13 @@ static void nrn_rhs(NrnThread* _nt) {
         vec_d[i] = 0.;
     }
 
+    if (_nt->_nrn_fast_imem) {
+        for (i = i1; i < i3; ++i) {
+            _nt->_nrn_fast_imem->_nrn_sav_rhs[i] = 0.;
+            _nt->_nrn_fast_imem->_nrn_sav_d[i] = 0.;
+        }
+    }
+
     nrn_ba(_nt, BEFORE_BREAKPOINT);
     /* note that CAP has no current */
     for (tml = _nt->tml; tml; tml = tml->next)
@@ -82,7 +89,15 @@ static void nrn_rhs(NrnThread* _nt) {
             }
 #endif
         }
-
+    if (_nt->_nrn_fast_imem) {
+        /* _nrn_save_rhs has only the contribution of electrode current
+           so here we transform so it only has membrane current contribution
+        */
+        double* p = _nt->_nrn_fast_imem->_nrn_sav_rhs;
+        for (i = i1; i < i3; ++i) {
+            p[i] -= vec_rhs[i];
+        }
+    }
 /* now the internal axial currents.
 The extracellular mechanism contribution is already done.
         rhs += ai_j*(vi_j - vi)
@@ -151,6 +166,17 @@ static void nrn_lhs(NrnThread* _nt) {
     double* vec_b = &(VEC_B(0));
     int* parent_index = _nt->_v_parent_index;
 
+
+    if (_nt->_nrn_fast_imem) {
+        /* _nrn_save_d has only the contribution of electrode current
+           so here we transform so it only has membrane current contribution
+        */
+        double* p = _nt->_nrn_fast_imem->_nrn_sav_d;
+        for (i = i1; i < i3; ++i) {
+            p[i] += vec_d[i];
+        }
+    }
+
 /* now add the axial currents */
 // clang-format off
     #pragma acc parallel loop present(          \

From b38ddf5461ed2f6dbfac086ffe90d7321b54627e Mon Sep 17 00:00:00 2001
From: Ioannis Magkanaris <iomagkanaris@gmail.com>
Date: Mon, 16 Sep 2019 16:20:40 +0200
Subject: [PATCH 02/14] Code refactoring

- Turned the struct _nrn_Fast_Imem of pointers to double arrays in
  NrnThread
- Moved fast_imem calculation related functions to seperate files
---
 coreneuron/nrniv/main1.cpp         |   8 +-
 coreneuron/nrnoc/capac.cpp         |   7 --
 coreneuron/nrnoc/fadvance_core.cpp |  25 +----
 coreneuron/nrnoc/fast_imem.cpp     | 147 +++++++++++++++++++++++++++++
 coreneuron/nrnoc/fast_imem.h       |  43 +++++++++
 coreneuron/nrnoc/finitialize.cpp   |   2 +-
 coreneuron/nrnoc/multicore.cpp     |  70 +-------------
 coreneuron/nrnoc/multicore.h       |   9 +-
 coreneuron/nrnoc/treeset_core.cpp  |  14 +--
 9 files changed, 206 insertions(+), 119 deletions(-)
 create mode 100644 coreneuron/nrnoc/fast_imem.cpp
 create mode 100644 coreneuron/nrnoc/fast_imem.h

diff --git a/coreneuron/nrniv/main1.cpp b/coreneuron/nrniv/main1.cpp
index cd37002aa..a764f4812 100644
--- a/coreneuron/nrniv/main1.cpp
+++ b/coreneuron/nrniv/main1.cpp
@@ -37,6 +37,7 @@ THE POSSIBILITY OF SUCH DAMAGE.
 #include "coreneuron/engine.h"
 #include "coreneuron/utils/randoms/nrnran123.h"
 #include "coreneuron/nrnconf.h"
+#include "coreneuron/nrnoc/fast_imem.h"
 #include "coreneuron/nrnoc/multicore.h"
 #include "coreneuron/nrnoc/nrnoc_decl.h"
 #include "coreneuron/nrnmpi/nrnmpi.h"
@@ -416,12 +417,6 @@ static void trajectory_return() {
     }
 }
 
-bool nrn_use_fast_imem;
-extern void nrn_fast_imem_alloc();
-static void use_fast_imem() {
-    nrn_use_fast_imem = true;
-    nrn_fast_imem_alloc();
-}
 
 }  // namespace coreneuron
 
@@ -509,6 +504,7 @@ extern "C" int run_solve_core(int argc, char** argv) {
         }
 
         if (nrnopt_get_flag("--fast_imem")) {
+            std::cout << "USING FAST I MEM" << std::endl;
             use_fast_imem();
         }
 
diff --git a/coreneuron/nrnoc/capac.cpp b/coreneuron/nrnoc/capac.cpp
index 7fafee179..e7b07e69f 100644
--- a/coreneuron/nrnoc/capac.cpp
+++ b/coreneuron/nrnoc/capac.cpp
@@ -206,13 +206,6 @@ void nrn_div_capacity(NrnThread* _nt, Memb_list* ml, int type) {
         VEC_RHS(ni[_iml]) /= 1.e-3 * cm;
         // fprintf(stderr, "== nrn_div_cap: RHS[%d]=%.12f\n", ni[_iml], VEC_RHS(ni[_iml])) ;
     }
-    /*if (_nt->_nrn_fast_imem) {
-        Node **vnode = ml->nodelist;
-        double* p = _nt->_nrn_fast_imem->_nrn_sav_rhs;
-        for (int i=0; i <_cntml_actual; ++i) {
-            p[vnode[i]->v_node_index] += i_cap;
-        }
-    }*/
 }
 
 void nrn_mul_capacity(NrnThread* _nt, Memb_list* ml, int type) {
diff --git a/coreneuron/nrnoc/fadvance_core.cpp b/coreneuron/nrnoc/fadvance_core.cpp
index a3f28d61d..b50417779 100644
--- a/coreneuron/nrnoc/fadvance_core.cpp
+++ b/coreneuron/nrnoc/fadvance_core.cpp
@@ -30,6 +30,7 @@ THE POSSIBILITY OF SUCH DAMAGE.
 #include "coreneuron/nrnconf.h"
 #include "coreneuron/nrnoc/multicore.h"
 #include "coreneuron/nrnmpi/nrnmpi.h"
+#include "coreneuron/nrnoc/fast_imem.h"
 #include "coreneuron/nrnoc/nrnoc_decl.h"
 #include "coreneuron/nrniv/nrn_acc_manager.h"
 #include "coreneuron/utils/reports/nrnreport.h"
@@ -155,31 +156,7 @@ static void* nrn_fixed_step_group_thread(NrnThread* nth) {
     return (void*)0;
 }
 
-void nrn_calc_fast_imem(NrnThread* _nt) {
-    int i;
-    int i1 = 0;
-    int i3 = _nt->end;
 
-    double* vec_rhs = &(VEC_RHS(0));
-    double* vec_area = &(VEC_AREA(0));
-
-    double* pd = _nt->_nrn_fast_imem->_nrn_sav_d;
-    double* prhs = _nt->_nrn_fast_imem->_nrn_sav_rhs;
-    FILE *fp_rhs, *fp_d;
-    fp_rhs = fopen("rhs.CORENEURON", "a");
-    fp_d = fopen("d.CORENEURON", "a");
-    fprintf(fp_rhs, "%.8e time\n", _nt->_t);
-    fprintf(fp_rhs, "%.8e time\n", _nt->_t);
-    for (i = i1; i < i3 ; ++i) {
-        prhs[i] = (pd[i]*vec_rhs[i] + prhs[i])*vec_area[i]*0.01;
-        fprintf(fp_rhs, "%.8e, ", prhs[i]);
-    }
-    fprintf(fp_rhs, "\n");
-    for (i = i1; i < i3 ; ++i) {
-        fprintf(fp_d, "%.8e, ", pd[i]);
-    }
-    fprintf(fp_d, "\n");
-}
 
 void update(NrnThread* _nt) {
     int i, i1, i2;
diff --git a/coreneuron/nrnoc/fast_imem.cpp b/coreneuron/nrnoc/fast_imem.cpp
new file mode 100644
index 000000000..8f7c09040
--- /dev/null
+++ b/coreneuron/nrnoc/fast_imem.cpp
@@ -0,0 +1,147 @@
+/*
+Copyright (c) 2019, Blue Brain Project
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification,
+are permitted provided that the following conditions are met:
+1. Redistributions of source code must retain the above copyright notice,
+   this list of conditions and the following disclaimer.
+2. Redistributions in binary form must reproduce the above copyright notice,
+   this list of conditions and the following disclaimer in the documentation
+   and/or other materials provided with the distribution.
+3. Neither the name of the copyright holder nor the names of its contributors
+   may be used to endorse or promote products derived from this software
+   without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
+LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
+THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+#include "coreneuron/nrnoc/fast_imem.h"
+
+namespace coreneuron {
+
+extern int nrn_nthread;
+extern NrnThread *nrn_threads;
+bool nrn_use_fast_imem;
+static int fast_imem_nthread_ = 0;
+static int *fast_imem_size_ = nullptr;
+static double **_nrn_Fast_Imem_rhs = nullptr;
+static double **_nrn_Fast_Imem_d = nullptr;
+
+static void fast_imem_free() {
+    int i;
+    for (i = 0; i < nrn_nthread; ++i) {
+        nrn_threads[i]._nrn_fast_imem_sav_rhs = nullptr;
+        nrn_threads[i]._nrn_fast_imem_sav_d = nullptr;
+    }
+
+    for (i = 0; i < fast_imem_nthread_; ++i) {
+        if (fast_imem_size_[i] > 0) {
+            free(_nrn_Fast_Imem_rhs[i]);
+            free(_nrn_Fast_Imem_d[i]);
+        }
+    }
+    if (fast_imem_nthread_) {
+        free(fast_imem_size_);
+        free(_nrn_Fast_Imem_rhs);
+        free(_nrn_Fast_Imem_d);
+        fast_imem_nthread_ = 0;
+        fast_imem_size_ = nullptr;
+        _nrn_Fast_Imem_rhs = nullptr;
+        _nrn_Fast_Imem_rhs = nullptr;
+    }
+}
+
+/*
+Avoid invalidating pointers to i_membrane_ unless the number of compartments
+in a thread has changed.
+*/
+
+static void fast_imem_alloc() {
+    int i;
+    if (fast_imem_nthread_ != nrn_nthread) {
+        fast_imem_free();
+        fast_imem_nthread_ = nrn_nthread;
+        fast_imem_size_ = (int *) ecalloc(nrn_nthread, sizeof(int));
+        _nrn_Fast_Imem_rhs = (double **) ecalloc(nrn_nthread, sizeof(double));
+        _nrn_Fast_Imem_d = (double **) ecalloc(nrn_nthread, sizeof(double *));
+    }
+    double **fi_rhs = &_nrn_Fast_Imem_rhs[i];
+    double **fi_d = &_nrn_Fast_Imem_d[i];
+    for (i = 0; i < nrn_nthread; ++i) {
+        NrnThread *nt = nrn_threads + i;
+        int n = nt->end;
+        if (n != fast_imem_size_[i]) {
+            if (fast_imem_size_[i] > 0) {
+                free(*fi_rhs);
+                free(*fi_d);
+            }
+            if (n > 0) {
+                *fi_rhs = (double *) emalloc_align(n * sizeof(double));
+                *fi_d = (double *) emalloc_align(n * sizeof(double));
+            }
+            fast_imem_size_[i] = n;
+        }
+    }
+}
+
+void nrn_fast_imem_alloc() {
+    if (nrn_use_fast_imem) {
+        int i;
+        fast_imem_alloc();
+        for (i = 0; i < nrn_nthread; ++i) {
+            nrn_threads[i]._nrn_fast_imem_sav_rhs = _nrn_Fast_Imem_rhs[i];
+            nrn_threads[i]._nrn_fast_imem_sav_d = _nrn_Fast_Imem_d[i];
+        }
+    } else {
+        fast_imem_free();
+    }
+}
+
+void use_fast_imem() {
+    nrn_use_fast_imem = true;
+    nrn_fast_imem_alloc();
+}
+
+void nrn_calc_fast_imem(NrnThread* _nt) {
+    int i;
+    int i1 = 0;
+    int i3 = _nt->end;
+
+    double* vec_rhs = &(VEC_RHS(0));
+    double* vec_area = &(VEC_AREA(0));
+
+    double* pd = _nt->_nrn_fast_imem_sav_d;
+    double* prhs = _nt->_nrn_fast_imem_sav_rhs;
+    FILE *fp_rhs, *fp_d;
+    char rhs_filename[20], d_filename[20];
+    sprintf(rhs_filename,"rhs.CORENEURON.%d",nrnmpi_myid);
+    sprintf(d_filename,"d.CORENEURON.%d",nrnmpi_myid);
+    fp_rhs = fopen(rhs_filename, "a");
+    fp_d = fopen(d_filename, "a");
+    fprintf(fp_rhs, "\n%.8e time\n", _nt->_t);
+    fprintf(fp_d, "\n%.8e time\n", _nt->_t);
+    for (i = i1; i < i3 ; ++i) {
+        prhs[i] = (pd[i]*vec_rhs[i] + prhs[i])*vec_area[i]*0.01;
+        fprintf(fp_rhs, "%.8e, ", prhs[i]);
+    }
+    fprintf(fp_rhs, "\n");
+    for (i = i1; i < i3 ; ++i) {
+        fprintf(fp_d, "%.8e, ", pd[i]);
+    }
+    fprintf(fp_d, "\n");
+    fclose(fp_rhs);
+    fclose(fp_d);
+}
+
+}
\ No newline at end of file
diff --git a/coreneuron/nrnoc/fast_imem.h b/coreneuron/nrnoc/fast_imem.h
new file mode 100644
index 000000000..31c3b8e17
--- /dev/null
+++ b/coreneuron/nrnoc/fast_imem.h
@@ -0,0 +1,43 @@
+//
+// Created by magkanar on 9/12/19.
+//
+
+#ifndef fast_imem_h
+#define fast_imem_h
+
+#include "coreneuron/nrnconf.h"
+#include "coreneuron/nrniv/memory.h"
+#include "coreneuron/nrnmpi/nrnmpi.h"
+#include "coreneuron/nrnoc/multicore.h"
+
+namespace coreneuron {
+
+/* Enables fast membrane curent culculation and allocates required
+ * memory.
+ * Found in src/nrncvode/cvodeobj.cpp in NEURON.
+ */
+void use_fast_imem();
+
+/* Free memory allocated for the fast current membrane calculation.
+ * Found in src/nrnoc/multicore.c in NEURON.
+ */
+static void fast_imem_free();
+
+/* Allocate memory for the rhs and d arrays needed for the fast
+ * current membrane calculation.
+ * Found in src/nrnoc/multicore.c in NEURON.
+ */
+static void fast_imem_alloc();
+
+/* fast_imem_alloc() wrapper.
+ * Found in src/nrnoc/multicore.c in NEURON.
+ */
+void nrn_fast_imem_alloc();
+
+/* Calculate the new values of rhs array at every timestep.
+ * Found in src/nrnoc/fadvance.c in NEURON.
+ */
+void nrn_calc_fast_imem(NrnThread* _nt);
+
+}  // namespace coreneuron
+#endif //fast_imem_h
diff --git a/coreneuron/nrnoc/finitialize.cpp b/coreneuron/nrnoc/finitialize.cpp
index 0e2fdfe04..448ff66bf 100644
--- a/coreneuron/nrnoc/finitialize.cpp
+++ b/coreneuron/nrnoc/finitialize.cpp
@@ -27,6 +27,7 @@ THE POSSIBILITY OF SUCH DAMAGE.
 */
 
 #include "coreneuron/nrnconf.h"
+#include "coreneuron/nrnoc/fast_imem.h"
 #include "coreneuron/nrnoc/multicore.h"
 #include "coreneuron/nrnoc/nrnoc_decl.h"
 #include "coreneuron/nrniv/profiler_interface.h"
@@ -34,7 +35,6 @@ THE POSSIBILITY OF SUCH DAMAGE.
 namespace coreneuron {
 
 extern bool nrn_use_fast_imem;
-extern void nrn_calc_fast_imem(NrnThread* _nt);
 
 void nrn_finitialize(int setv, double v) {
     int i;
diff --git a/coreneuron/nrnoc/multicore.cpp b/coreneuron/nrnoc/multicore.cpp
index b98513998..8a473ff3e 100644
--- a/coreneuron/nrnoc/multicore.cpp
+++ b/coreneuron/nrnoc/multicore.cpp
@@ -138,7 +138,8 @@ void nrn_threads_create(int n) {
                 nt->mapping = NULL;
                 nt->trajec_requests = NULL;
 
-                nt->_nrn_fast_imem = 0;
+                nt->_nrn_fast_imem_sav_rhs = nullptr;
+                nt->_nrn_fast_imem_sav_d = nullptr;
             }
         }
         v_structure_change = 1;
@@ -147,73 +148,6 @@ void nrn_threads_create(int n) {
     /*printf("nrn_threads_create %d %d\n", nrn_nthread, nrn_thread_parallel_);*/
 }
 
-/*
-Avoid invalidating pointers to i_membrane_ unless the number of compartments
-in a thread has changed.
-*/
-extern bool nrn_use_fast_imem;
-static int fast_imem_nthread_ = 0;
-static int* fast_imem_size_ = NULL;
-static _nrn_Fast_Imem* fast_imem_;
-
-static void fast_imem_free() {
-    int i;
-    for (i = 0; i < nrn_nthread; ++i) {
-        nrn_threads[i]._nrn_fast_imem = NULL;
-    }
-    for (i = 0; i < fast_imem_nthread_; ++i) {
-        if (fast_imem_size_[i] > 0) {
-            free(fast_imem_[i]._nrn_sav_rhs);
-            free(fast_imem_[i]._nrn_sav_d);
-        }
-    }
-    if (fast_imem_nthread_) {
-        free(fast_imem_size_);
-        free(fast_imem_);
-        fast_imem_nthread_ = 0;
-        fast_imem_size_ = NULL;
-        fast_imem_ = NULL;
-    }
-}
-
-static void fast_imem_alloc() {
-    int i;
-    if (fast_imem_nthread_ != nrn_nthread) {
-        fast_imem_free();
-        fast_imem_nthread_ = nrn_nthread;
-        fast_imem_size_ = (int*)ecalloc(nrn_nthread, sizeof(int));
-        fast_imem_ = (_nrn_Fast_Imem*)ecalloc(nrn_nthread, sizeof(_nrn_Fast_Imem));
-    }
-    for (i=0; i < nrn_nthread; ++i) {
-        NrnThread* nt = nrn_threads + i;
-        int n = nt->end;
-        _nrn_Fast_Imem* fi = fast_imem_ + i;
-        if (n != fast_imem_size_[i]) {
-            if (fast_imem_size_[i] > 0) {
-                free(fi->_nrn_sav_rhs);
-                free(fi->_nrn_sav_d);
-            }
-            if (n > 0) {
-                fi->_nrn_sav_rhs = (double*)emalloc_align(n * sizeof(double));
-                fi->_nrn_sav_d = (double*)emalloc_align(n * sizeof(double));
-            }
-            fast_imem_size_[i] = n;
-        }
-    }
-}
-
-void nrn_fast_imem_alloc() {
-    if (nrn_use_fast_imem) {
-        int i;
-        fast_imem_alloc();
-        for (i=0; i < nrn_nthread; ++i) {
-            nrn_threads[i]._nrn_fast_imem = fast_imem_ + i;
-        }
-    }else{
-        fast_imem_free();
-    }
-}
-
 void nrn_threads_free() {
     if (nrn_nthread) {
         free_memory((void*)nrn_threads);
diff --git a/coreneuron/nrnoc/multicore.h b/coreneuron/nrnoc/multicore.h
index 7272fb5bb..e25e6d5b1 100644
--- a/coreneuron/nrnoc/multicore.h
+++ b/coreneuron/nrnoc/multicore.h
@@ -59,11 +59,6 @@ struct NrnThreadBAList {
     NrnThreadBAList* next;
 };
 
-typedef struct _nrn_Fast_Imem {
-    double* _nrn_sav_rhs;
-    double* _nrn_sav_d;
-} _nrn_Fast_Imem;
-
 struct TrajectoryRequests {
     void** vpr;       /* PlayRecord Objects known by NEURON */
     double** scatter; /* if bsize == 0, each time step */
@@ -122,7 +117,9 @@ struct NrnThread {
     double* _shadow_d;    /* Not pointer into _data. Avoid race for multiple POINT_PROCESS in same
                              compartment */
 
-    _nrn_Fast_Imem* _nrn_fast_imem;
+    /* Fast membrane current calculation arrays */
+    double* _nrn_fast_imem_sav_rhs; /* Save equation right hand side array */
+    double* _nrn_fast_imem_sav_d;   /* Save diagonal array */
 
     int* _v_parent_index;
     int* _permute;
diff --git a/coreneuron/nrnoc/treeset_core.cpp b/coreneuron/nrnoc/treeset_core.cpp
index fe27a0304..fa0c3bb62 100644
--- a/coreneuron/nrnoc/treeset_core.cpp
+++ b/coreneuron/nrnoc/treeset_core.cpp
@@ -67,10 +67,10 @@ static void nrn_rhs(NrnThread* _nt) {
         vec_d[i] = 0.;
     }
 
-    if (_nt->_nrn_fast_imem) {
+    if (_nt->_nrn_fast_imem_sav_rhs && _nt->_nrn_fast_imem_sav_d) {
         for (i = i1; i < i3; ++i) {
-            _nt->_nrn_fast_imem->_nrn_sav_rhs[i] = 0.;
-            _nt->_nrn_fast_imem->_nrn_sav_d[i] = 0.;
+            _nt->_nrn_fast_imem_sav_rhs[i] = 0.;
+            _nt->_nrn_fast_imem_sav_d[i] = 0.;
         }
     }
 
@@ -89,11 +89,11 @@ static void nrn_rhs(NrnThread* _nt) {
             }
 #endif
         }
-    if (_nt->_nrn_fast_imem) {
+    if (_nt->_nrn_fast_imem_sav_rhs) {
         /* _nrn_save_rhs has only the contribution of electrode current
            so here we transform so it only has membrane current contribution
         */
-        double* p = _nt->_nrn_fast_imem->_nrn_sav_rhs;
+        double* p = _nt->_nrn_fast_imem_sav_rhs;
         for (i = i1; i < i3; ++i) {
             p[i] -= vec_rhs[i];
         }
@@ -167,11 +167,11 @@ static void nrn_lhs(NrnThread* _nt) {
     int* parent_index = _nt->_v_parent_index;
 
 
-    if (_nt->_nrn_fast_imem) {
+    if (_nt->_nrn_fast_imem_sav_d) {
         /* _nrn_save_d has only the contribution of electrode current
            so here we transform so it only has membrane current contribution
         */
-        double* p = _nt->_nrn_fast_imem->_nrn_sav_d;
+        double* p = _nt->_nrn_fast_imem_sav_d;
         for (i = i1; i < i3; ++i) {
             p[i] += vec_d[i];
         }

From 768e9beb55c81e2950ed74fc203a0eaae09410ea Mon Sep 17 00:00:00 2001
From: Ioannis Magkanaris <iomagkanaris@gmail.com>
Date: Tue, 17 Sep 2019 15:43:17 +0200
Subject: [PATCH 03/14] Reverted back to NrnFastImem struct instead of two
 double arrays

---
 coreneuron/nrniv/main1.cpp                    |  5 --
 coreneuron/nrnoc/fast_imem.cpp                | 51 +++++++++----------
 coreneuron/nrnoc/fast_imem.h                  | 33 +++++++++---
 coreneuron/nrnoc/multicore.cpp                |  3 +-
 coreneuron/nrnoc/multicore.h                  | 10 ++--
 coreneuron/nrnoc/treeset_core.cpp             | 14 ++---
 .../reports/report_configuration_parser.cpp   |  5 ++
 7 files changed, 70 insertions(+), 51 deletions(-)

diff --git a/coreneuron/nrniv/main1.cpp b/coreneuron/nrniv/main1.cpp
index a764f4812..f579463a4 100644
--- a/coreneuron/nrniv/main1.cpp
+++ b/coreneuron/nrniv/main1.cpp
@@ -503,11 +503,6 @@ extern "C" int run_solve_core(int argc, char** argv) {
             (*nrn2core_part2_clean_)();
         }
 
-        if (nrnopt_get_flag("--fast_imem")) {
-            std::cout << "USING FAST I MEM" << std::endl;
-            use_fast_imem();
-        }
-
         // TODO : if some ranks are empty then restore will go in deadlock
         // phase (as some ranks won't have restored anything and hence return
         // false in checkpoint_initialize
diff --git a/coreneuron/nrnoc/fast_imem.cpp b/coreneuron/nrnoc/fast_imem.cpp
index 8f7c09040..019eb9d2d 100644
--- a/coreneuron/nrnoc/fast_imem.cpp
+++ b/coreneuron/nrnoc/fast_imem.cpp
@@ -26,7 +26,10 @@ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
 THE POSSIBILITY OF SUCH DAMAGE.
 */
 
+#include "coreneuron/nrnconf.h"
 #include "coreneuron/nrnoc/fast_imem.h"
+#include "coreneuron/nrniv/memory.h"
+#include "coreneuron/nrnmpi/nrnmpi.h"
 
 namespace coreneuron {
 
@@ -35,30 +38,25 @@ extern NrnThread *nrn_threads;
 bool nrn_use_fast_imem;
 static int fast_imem_nthread_ = 0;
 static int *fast_imem_size_ = nullptr;
-static double **_nrn_Fast_Imem_rhs = nullptr;
-static double **_nrn_Fast_Imem_d = nullptr;
+static NrnFastImem* fast_imem_;
 
 static void fast_imem_free() {
     int i;
     for (i = 0; i < nrn_nthread; ++i) {
-        nrn_threads[i]._nrn_fast_imem_sav_rhs = nullptr;
-        nrn_threads[i]._nrn_fast_imem_sav_d = nullptr;
+        nrn_threads[i].nrn_fast_imem = NULL;
     }
-
     for (i = 0; i < fast_imem_nthread_; ++i) {
         if (fast_imem_size_[i] > 0) {
-            free(_nrn_Fast_Imem_rhs[i]);
-            free(_nrn_Fast_Imem_d[i]);
+            free(fast_imem_[i].nrn_sav_rhs);
+            free(fast_imem_[i].nrn_sav_d);
         }
     }
     if (fast_imem_nthread_) {
         free(fast_imem_size_);
-        free(_nrn_Fast_Imem_rhs);
-        free(_nrn_Fast_Imem_d);
+        free(fast_imem_);
         fast_imem_nthread_ = 0;
         fast_imem_size_ = nullptr;
-        _nrn_Fast_Imem_rhs = nullptr;
-        _nrn_Fast_Imem_rhs = nullptr;
+        fast_imem_ = nullptr;
     }
 }
 
@@ -72,23 +70,21 @@ static void fast_imem_alloc() {
     if (fast_imem_nthread_ != nrn_nthread) {
         fast_imem_free();
         fast_imem_nthread_ = nrn_nthread;
-        fast_imem_size_ = (int *) ecalloc(nrn_nthread, sizeof(int));
-        _nrn_Fast_Imem_rhs = (double **) ecalloc(nrn_nthread, sizeof(double));
-        _nrn_Fast_Imem_d = (double **) ecalloc(nrn_nthread, sizeof(double *));
+        fast_imem_size_ = (int*)ecalloc(nrn_nthread, sizeof(int));
+        fast_imem_ = (NrnFastImem*)ecalloc(nrn_nthread, sizeof(NrnFastImem));
     }
-    double **fi_rhs = &_nrn_Fast_Imem_rhs[i];
-    double **fi_d = &_nrn_Fast_Imem_d[i];
-    for (i = 0; i < nrn_nthread; ++i) {
-        NrnThread *nt = nrn_threads + i;
+    for (i=0; i < nrn_nthread; ++i) {
+        NrnThread* nt = nrn_threads + i;
         int n = nt->end;
+        NrnFastImem* fi = fast_imem_ + i;
         if (n != fast_imem_size_[i]) {
             if (fast_imem_size_[i] > 0) {
-                free(*fi_rhs);
-                free(*fi_d);
+                free(fi->nrn_sav_rhs);
+                free(fi->nrn_sav_d);
             }
             if (n > 0) {
-                *fi_rhs = (double *) emalloc_align(n * sizeof(double));
-                *fi_d = (double *) emalloc_align(n * sizeof(double));
+                fi->nrn_sav_rhs = (double*)emalloc_align(n * sizeof(double));
+                fi->nrn_sav_d = (double*)emalloc_align(n * sizeof(double));
             }
             fast_imem_size_[i] = n;
         }
@@ -99,11 +95,10 @@ void nrn_fast_imem_alloc() {
     if (nrn_use_fast_imem) {
         int i;
         fast_imem_alloc();
-        for (i = 0; i < nrn_nthread; ++i) {
-            nrn_threads[i]._nrn_fast_imem_sav_rhs = _nrn_Fast_Imem_rhs[i];
-            nrn_threads[i]._nrn_fast_imem_sav_d = _nrn_Fast_Imem_d[i];
+        for (i=0; i < nrn_nthread; ++i) {
+            nrn_threads[i].nrn_fast_imem = fast_imem_ + i;
         }
-    } else {
+    }else{
         fast_imem_free();
     }
 }
@@ -121,8 +116,8 @@ void nrn_calc_fast_imem(NrnThread* _nt) {
     double* vec_rhs = &(VEC_RHS(0));
     double* vec_area = &(VEC_AREA(0));
 
-    double* pd = _nt->_nrn_fast_imem_sav_d;
-    double* prhs = _nt->_nrn_fast_imem_sav_rhs;
+    double* pd = _nt->nrn_fast_imem->nrn_sav_d;
+    double* prhs = _nt->nrn_fast_imem->nrn_sav_rhs;
     FILE *fp_rhs, *fp_d;
     char rhs_filename[20], d_filename[20];
     sprintf(rhs_filename,"rhs.CORENEURON.%d",nrnmpi_myid);
diff --git a/coreneuron/nrnoc/fast_imem.h b/coreneuron/nrnoc/fast_imem.h
index 31c3b8e17..4e54c1512 100644
--- a/coreneuron/nrnoc/fast_imem.h
+++ b/coreneuron/nrnoc/fast_imem.h
@@ -1,13 +1,34 @@
-//
-// Created by magkanar on 9/12/19.
-//
+/*
+Copyright (c) 2019, Blue Brain Project
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification,
+are permitted provided that the following conditions are met:
+1. Redistributions of source code must retain the above copyright notice,
+   this list of conditions and the following disclaimer.
+2. Redistributions in binary form must reproduce the above copyright notice,
+   this list of conditions and the following disclaimer in the documentation
+   and/or other materials provided with the distribution.
+3. Neither the name of the copyright holder nor the names of its contributors
+   may be used to endorse or promote products derived from this software
+   without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
+LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
+THE POSSIBILITY OF SUCH DAMAGE.
+*/
 
 #ifndef fast_imem_h
 #define fast_imem_h
 
-#include "coreneuron/nrnconf.h"
-#include "coreneuron/nrniv/memory.h"
-#include "coreneuron/nrnmpi/nrnmpi.h"
 #include "coreneuron/nrnoc/multicore.h"
 
 namespace coreneuron {
diff --git a/coreneuron/nrnoc/multicore.cpp b/coreneuron/nrnoc/multicore.cpp
index 8a473ff3e..e3bc03e03 100644
--- a/coreneuron/nrnoc/multicore.cpp
+++ b/coreneuron/nrnoc/multicore.cpp
@@ -138,8 +138,7 @@ void nrn_threads_create(int n) {
                 nt->mapping = NULL;
                 nt->trajec_requests = NULL;
 
-                nt->_nrn_fast_imem_sav_rhs = nullptr;
-                nt->_nrn_fast_imem_sav_d = nullptr;
+                nt->nrn_fast_imem = nullptr;
             }
         }
         v_structure_change = 1;
diff --git a/coreneuron/nrnoc/multicore.h b/coreneuron/nrnoc/multicore.h
index e25e6d5b1..9aaec8bf0 100644
--- a/coreneuron/nrnoc/multicore.h
+++ b/coreneuron/nrnoc/multicore.h
@@ -59,6 +59,11 @@ struct NrnThreadBAList {
     NrnThreadBAList* next;
 };
 
+typedef struct NrnFastImem {
+    double* nrn_sav_rhs;
+    double* nrn_sav_d;
+} NrnFastImem;
+
 struct TrajectoryRequests {
     void** vpr;       /* PlayRecord Objects known by NEURON */
     double** scatter; /* if bsize == 0, each time step */
@@ -117,9 +122,8 @@ struct NrnThread {
     double* _shadow_d;    /* Not pointer into _data. Avoid race for multiple POINT_PROCESS in same
                              compartment */
 
-    /* Fast membrane current calculation arrays */
-    double* _nrn_fast_imem_sav_rhs; /* Save equation right hand side array */
-    double* _nrn_fast_imem_sav_d;   /* Save diagonal array */
+    /* Fast membrane current calculation struct */
+    NrnFastImem* nrn_fast_imem;
 
     int* _v_parent_index;
     int* _permute;
diff --git a/coreneuron/nrnoc/treeset_core.cpp b/coreneuron/nrnoc/treeset_core.cpp
index fa0c3bb62..bfdb8cd0a 100644
--- a/coreneuron/nrnoc/treeset_core.cpp
+++ b/coreneuron/nrnoc/treeset_core.cpp
@@ -67,10 +67,10 @@ static void nrn_rhs(NrnThread* _nt) {
         vec_d[i] = 0.;
     }
 
-    if (_nt->_nrn_fast_imem_sav_rhs && _nt->_nrn_fast_imem_sav_d) {
+    if (_nt->nrn_fast_imem) {
         for (i = i1; i < i3; ++i) {
-            _nt->_nrn_fast_imem_sav_rhs[i] = 0.;
-            _nt->_nrn_fast_imem_sav_d[i] = 0.;
+            _nt->nrn_fast_imem->nrn_sav_rhs[i] = 0.;
+            _nt->nrn_fast_imem->nrn_sav_d[i] = 0.;
         }
     }
 
@@ -89,11 +89,11 @@ static void nrn_rhs(NrnThread* _nt) {
             }
 #endif
         }
-    if (_nt->_nrn_fast_imem_sav_rhs) {
+    if (_nt->nrn_fast_imem) {
         /* _nrn_save_rhs has only the contribution of electrode current
            so here we transform so it only has membrane current contribution
         */
-        double* p = _nt->_nrn_fast_imem_sav_rhs;
+        double* p = _nt->nrn_fast_imem->nrn_sav_rhs;
         for (i = i1; i < i3; ++i) {
             p[i] -= vec_rhs[i];
         }
@@ -167,11 +167,11 @@ static void nrn_lhs(NrnThread* _nt) {
     int* parent_index = _nt->_v_parent_index;
 
 
-    if (_nt->_nrn_fast_imem_sav_d) {
+    if (_nt->nrn_fast_imem) {
         /* _nrn_save_d has only the contribution of electrode current
            so here we transform so it only has membrane current contribution
         */
-        double* p = _nt->_nrn_fast_imem_sav_d;
+        double* p = _nt->nrn_fast_imem->nrn_sav_d;
         for (i = i1; i < i3; ++i) {
             p[i] += vec_d[i];
         }
diff --git a/coreneuron/utils/reports/report_configuration_parser.cpp b/coreneuron/utils/reports/report_configuration_parser.cpp
index 901469e52..ad0fb61e4 100644
--- a/coreneuron/utils/reports/report_configuration_parser.cpp
+++ b/coreneuron/utils/reports/report_configuration_parser.cpp
@@ -28,6 +28,7 @@
 
 #include "coreneuron/nrniv/nrn_assert.h"
 #include "coreneuron/utils/reports/nrnreport.h"
+#include "coreneuron/nrnoc/fast_imem.h"
 #include "coreneuron/nrnoc/mech_mapping.hpp"
 #include <stdio.h>
 #include <stdlib.h>
@@ -101,6 +102,10 @@ std::vector<ReportConfiguration> create_report_configurations(const char* conf_f
             abort();
         }
 
+        if (strcmp(report_on, "i_membrane_")) {
+            use_fast_imem();
+        }
+
         if (report.type == SynapseReport)
             parse_filter_string(report_on, report);
 

From 3d904d0caf91f71af88592437108e56bf23d0f39 Mon Sep 17 00:00:00 2001
From: Ioannis Magkanaris <ioannis.magkanaris@epfl.ch>
Date: Wed, 18 Sep 2019 12:28:17 +0200
Subject: [PATCH 04/14] Small code changes

- Related to enabling fast_imem calculation from report parser
- Removed some debugging code
---
 coreneuron/nrniv/nrn_setup.cpp |  4 ++++
 coreneuron/nrnoc/fast_imem.cpp | 19 ++-----------------
 2 files changed, 6 insertions(+), 17 deletions(-)

diff --git a/coreneuron/nrniv/nrn_setup.cpp b/coreneuron/nrniv/nrn_setup.cpp
index dbb3605d3..031742e0c 100644
--- a/coreneuron/nrniv/nrn_setup.cpp
+++ b/coreneuron/nrniv/nrn_setup.cpp
@@ -34,6 +34,7 @@ THE POSSIBILITY OF SUCH DAMAGE.
 #include "coreneuron/nrnconf.h"
 #include "coreneuron/nrnoc/multicore.h"
 #include "coreneuron/nrniv/nrniv_decl.h"
+#include "coreneuron/nrnoc/fast_imem.h"
 #include "coreneuron/nrnoc/nrnoc_decl.h"
 #include "coreneuron/nrniv/vrecitem.h"
 #include "coreneuron/nrniv/multisend.h"
@@ -795,6 +796,9 @@ void nrn_setup(const char* filesdat,
     mk_cell_indices();
 #endif
 
+    /// Allocate memory for fast_imem calculation
+    nrn_fast_imem_alloc();
+
     /// Generally, tables depend on a few parameters. And if those parameters change,
     /// then the table needs to be recomputed. This is obviously important in NEURON
     /// since the user can change those parameters at any time. However, there is no
diff --git a/coreneuron/nrnoc/fast_imem.cpp b/coreneuron/nrnoc/fast_imem.cpp
index 019eb9d2d..d545cff34 100644
--- a/coreneuron/nrnoc/fast_imem.cpp
+++ b/coreneuron/nrnoc/fast_imem.cpp
@@ -118,25 +118,10 @@ void nrn_calc_fast_imem(NrnThread* _nt) {
 
     double* pd = _nt->nrn_fast_imem->nrn_sav_d;
     double* prhs = _nt->nrn_fast_imem->nrn_sav_rhs;
-    FILE *fp_rhs, *fp_d;
-    char rhs_filename[20], d_filename[20];
-    sprintf(rhs_filename,"rhs.CORENEURON.%d",nrnmpi_myid);
-    sprintf(d_filename,"d.CORENEURON.%d",nrnmpi_myid);
-    fp_rhs = fopen(rhs_filename, "a");
-    fp_d = fopen(d_filename, "a");
-    fprintf(fp_rhs, "\n%.8e time\n", _nt->_t);
-    fprintf(fp_d, "\n%.8e time\n", _nt->_t);
     for (i = i1; i < i3 ; ++i) {
         prhs[i] = (pd[i]*vec_rhs[i] + prhs[i])*vec_area[i]*0.01;
-        fprintf(fp_rhs, "%.8e, ", prhs[i]);
     }
-    fprintf(fp_rhs, "\n");
-    for (i = i1; i < i3 ; ++i) {
-        fprintf(fp_d, "%.8e, ", pd[i]);
-    }
-    fprintf(fp_d, "\n");
-    fclose(fp_rhs);
-    fclose(fp_d);
 }
 
-}
\ No newline at end of file
+}
+

From abf194f758405c2a9e02870f7f7f3a969266a1fe Mon Sep 17 00:00:00 2001
From: Ioannis Magkanaris <ioannis.magkanaris@epfl.ch>
Date: Wed, 18 Sep 2019 13:27:22 +0200
Subject: [PATCH 05/14] Little code polishing

---
 coreneuron/nrniv/main1.cpp                               | 2 --
 coreneuron/nrniv/nrnoptarg.cpp                           | 1 -
 coreneuron/nrnoc/fadvance_core.cpp                       | 6 +++---
 coreneuron/nrnoc/treeset_core.cpp                        | 3 ++-
 coreneuron/utils/reports/report_configuration_parser.cpp | 5 -----
 5 files changed, 5 insertions(+), 12 deletions(-)

diff --git a/coreneuron/nrniv/main1.cpp b/coreneuron/nrniv/main1.cpp
index f579463a4..a7ad6c0f2 100644
--- a/coreneuron/nrniv/main1.cpp
+++ b/coreneuron/nrniv/main1.cpp
@@ -37,7 +37,6 @@ THE POSSIBILITY OF SUCH DAMAGE.
 #include "coreneuron/engine.h"
 #include "coreneuron/utils/randoms/nrnran123.h"
 #include "coreneuron/nrnconf.h"
-#include "coreneuron/nrnoc/fast_imem.h"
 #include "coreneuron/nrnoc/multicore.h"
 #include "coreneuron/nrnoc/nrnoc_decl.h"
 #include "coreneuron/nrnmpi/nrnmpi.h"
@@ -417,7 +416,6 @@ static void trajectory_return() {
     }
 }
 
-
 }  // namespace coreneuron
 
 /// The following high-level functions are marked as "extern C"
diff --git a/coreneuron/nrniv/nrnoptarg.cpp b/coreneuron/nrniv/nrnoptarg.cpp
index 93f9ecb8e..0099e476b 100644
--- a/coreneuron/nrniv/nrnoptarg.cpp
+++ b/coreneuron/nrniv/nrnoptarg.cpp
@@ -98,7 +98,6 @@ static param_flag param_flag_args[] = {
     {"--multisend", "Use Multisend spike exchange instead of Allgather."},
     {"--binqueue", "Use bin queue."},
     {"--skip-mpi-finalize", "Do not call mpi finalize."},
-    {"--fast_imem", "Use fast_i_membrane."},
     {NULL, NULL}};
 
 static param_str param_str_args[] = {
diff --git a/coreneuron/nrnoc/fadvance_core.cpp b/coreneuron/nrnoc/fadvance_core.cpp
index b50417779..7579d0c1d 100644
--- a/coreneuron/nrnoc/fadvance_core.cpp
+++ b/coreneuron/nrnoc/fadvance_core.cpp
@@ -156,8 +156,6 @@ static void* nrn_fixed_step_group_thread(NrnThread* nth) {
     return (void*)0;
 }
 
-
-
 void update(NrnThread* _nt) {
     int i, i1, i2;
     i1 = 0;
@@ -195,7 +193,9 @@ void update(NrnThread* _nt) {
         assert(_nt->tml->index == CAP);
         nrn_cur_capacitance(_nt, _nt->tml->ml, _nt->tml->index);
     }
-    if (nrn_use_fast_imem) { nrn_calc_fast_imem(_nt); }
+    if (nrn_use_fast_imem) { 
+        nrn_calc_fast_imem(_nt);
+    }
 }
 
 void nonvint(NrnThread* _nt) {
diff --git a/coreneuron/nrnoc/treeset_core.cpp b/coreneuron/nrnoc/treeset_core.cpp
index bfdb8cd0a..135979854 100644
--- a/coreneuron/nrnoc/treeset_core.cpp
+++ b/coreneuron/nrnoc/treeset_core.cpp
@@ -89,6 +89,7 @@ static void nrn_rhs(NrnThread* _nt) {
             }
 #endif
         }
+    
     if (_nt->nrn_fast_imem) {
         /* _nrn_save_rhs has only the contribution of electrode current
            so here we transform so it only has membrane current contribution
@@ -98,6 +99,7 @@ static void nrn_rhs(NrnThread* _nt) {
             p[i] -= vec_rhs[i];
         }
     }
+
 /* now the internal axial currents.
 The extracellular mechanism contribution is already done.
         rhs += ai_j*(vi_j - vi)
@@ -166,7 +168,6 @@ static void nrn_lhs(NrnThread* _nt) {
     double* vec_b = &(VEC_B(0));
     int* parent_index = _nt->_v_parent_index;
 
-
     if (_nt->nrn_fast_imem) {
         /* _nrn_save_d has only the contribution of electrode current
            so here we transform so it only has membrane current contribution
diff --git a/coreneuron/utils/reports/report_configuration_parser.cpp b/coreneuron/utils/reports/report_configuration_parser.cpp
index ad0fb61e4..901469e52 100644
--- a/coreneuron/utils/reports/report_configuration_parser.cpp
+++ b/coreneuron/utils/reports/report_configuration_parser.cpp
@@ -28,7 +28,6 @@
 
 #include "coreneuron/nrniv/nrn_assert.h"
 #include "coreneuron/utils/reports/nrnreport.h"
-#include "coreneuron/nrnoc/fast_imem.h"
 #include "coreneuron/nrnoc/mech_mapping.hpp"
 #include <stdio.h>
 #include <stdlib.h>
@@ -102,10 +101,6 @@ std::vector<ReportConfiguration> create_report_configurations(const char* conf_f
             abort();
         }
 
-        if (strcmp(report_on, "i_membrane_")) {
-            use_fast_imem();
-        }
-
         if (report.type == SynapseReport)
             parse_filter_string(report_on, report);
 

From 891dd1a9301327d64fc94d93e56683c34248ef7e Mon Sep 17 00:00:00 2001
From: Ioannis Magkanaris <iomagkanaris@gmail.com>
Date: Fri, 20 Sep 2019 14:04:24 +0200
Subject: [PATCH 06/14] Added i_membrane to neuron direct memory transfer and
 refactored fast_imem memory allocation

---
 coreneuron/nrniv/main1.cpp         |  6 ++-
 coreneuron/nrniv/nrn_setup.cpp     | 18 +++++++-
 coreneuron/nrniv/nrniv_decl.h      |  4 ++
 coreneuron/nrnoc/fadvance_core.cpp |  2 -
 coreneuron/nrnoc/fast_imem.cpp     | 71 ++++++------------------------
 coreneuron/nrnoc/fast_imem.h       |  9 ++--
 coreneuron/nrnoc/finitialize.cpp   |  2 -
 7 files changed, 42 insertions(+), 70 deletions(-)

diff --git a/coreneuron/nrniv/main1.cpp b/coreneuron/nrniv/main1.cpp
index a7ad6c0f2..fa8d83690 100644
--- a/coreneuron/nrniv/main1.cpp
+++ b/coreneuron/nrniv/main1.cpp
@@ -37,6 +37,7 @@ THE POSSIBILITY OF SUCH DAMAGE.
 #include "coreneuron/engine.h"
 #include "coreneuron/utils/randoms/nrnran123.h"
 #include "coreneuron/nrnconf.h"
+#include "coreneuron/nrnoc/fast_imem.h"
 #include "coreneuron/nrnoc/multicore.h"
 #include "coreneuron/nrnoc/nrnoc_decl.h"
 #include "coreneuron/nrnmpi/nrnmpi.h"
@@ -126,10 +127,13 @@ char* prepare_args(int& argc, char**& argv, int use_mpi, const char* arg) {
     return first;
 }
 
-int corenrn_embedded_run(int nthread, int have_gaps, int use_mpi, const char* arg) {
+int corenrn_embedded_run(int nthread, int have_gaps, int use_mpi, int use_fast_imem, const char* arg) {
     corenrn_embedded = 1;
     corenrn_embedded_nthread = nthread;
     coreneuron::nrn_have_gaps = have_gaps;
+    if (use_fast_imem) {
+        coreneuron::nrn_use_fast_imem = true;
+    }
 
     set_openmp_threads(nthread);
     int argc = 0;
diff --git a/coreneuron/nrniv/nrn_setup.cpp b/coreneuron/nrniv/nrn_setup.cpp
index 031742e0c..be9c835b5 100644
--- a/coreneuron/nrniv/nrn_setup.cpp
+++ b/coreneuron/nrniv/nrn_setup.cpp
@@ -51,6 +51,8 @@ THE POSSIBILITY OF SUCH DAMAGE.
 #include "coreneuron/utils/reports/nrnsection_mapping.h"
 
 // callbacks into nrn/src/nrniv/nrnbbcore_write.cpp
+#include "coreneuron/nrnoc/fast_imem.h"
+#include "coreneuron/nrniv/nrniv_decl.h"
 #include "coreneuron/nrniv/nrn2core_direct.h"
 
 int corenrn_embedded;
@@ -917,9 +919,9 @@ int nrn_i_layout(int icnt, int cnt, int isz, int sz, int layout) {
 }
 
 // take into account alignment, layout, permutation
-// only voltage or mechanism data index allowed. (mtype 0 means time)
+// only voltage, i_membrane or mechanism data index allowed. (mtype 0 means time)
 double* stdindex2ptr(int mtype, int index, NrnThread& nt) {
-    if (mtype == -5) {  // voltage
+    if (mtype == voltage) {  // voltage
         int v0 = nt._actual_v - nt._data;
         int ix = index;  // relative to _actual_v
         nrn_assert((ix >= 0) && (ix < nt.end));
@@ -927,6 +929,14 @@ double* stdindex2ptr(int mtype, int index, NrnThread& nt) {
             node_permute(&ix, 1, nt._permute);
         }
         return nt._data + (v0 + ix);                // relative to nt._data
+    } else if (mtype == i_membrane) {  // membrane current from fast_imem calculation
+            int i_mem = nt.nrn_fast_imem->nrn_sav_rhs - nt._data;
+            int ix = index;  // relative to _actual_v
+            nrn_assert((ix >= 0) && (ix < nt.end));
+            if (nt._permute) {
+                node_permute(&ix, 1, nt._permute);
+            }
+            return nt._data + (i_mem + ix);         // relative to nt._data
     } else if (mtype > 0 && mtype < n_memb_func) {  //
         Memb_list* ml = nt._ml_list[mtype];
         nrn_assert(ml);
@@ -1154,6 +1164,10 @@ void nrn_cleanup(bool clean_ion_global_map) {
         }
 
         free_memory(nt->_ml_list);
+
+        if (nt->nrn_fast_imem) {
+            fast_imem_free();
+        }
     }
 
 #if NRN_MULTISEND
diff --git a/coreneuron/nrniv/nrniv_decl.h b/coreneuron/nrniv/nrniv_decl.h
index dcb4c9a53..5be1c1b9b 100644
--- a/coreneuron/nrniv/nrniv_decl.h
+++ b/coreneuron/nrniv/nrniv_decl.h
@@ -35,6 +35,10 @@ THE POSSIBILITY OF SUCH DAMAGE.
 #include "coreneuron/utils/endianness.h"
 #include "coreneuron/nrniv/nrnoptarg.h"
 namespace coreneuron {
+
+/// Mechanism type to be used from stdindex2ptr and nrn_dblpntr2nrncore (in Neuron)
+enum mech_type {voltage = -1, i_membrane = -2};
+
 extern int cvode_active_;
 /// Vector of maps for negative presyns
 extern std::vector<std::map<int, PreSyn*> > neg_gid2out;
diff --git a/coreneuron/nrnoc/fadvance_core.cpp b/coreneuron/nrnoc/fadvance_core.cpp
index 7579d0c1d..179354cd1 100644
--- a/coreneuron/nrnoc/fadvance_core.cpp
+++ b/coreneuron/nrnoc/fadvance_core.cpp
@@ -40,8 +40,6 @@ THE POSSIBILITY OF SUCH DAMAGE.
 
 namespace coreneuron {
 
-extern bool nrn_use_fast_imem;
-
 static void* nrn_fixed_step_thread(NrnThread*);
 static void* nrn_fixed_step_group_thread(NrnThread*);
 
diff --git a/coreneuron/nrnoc/fast_imem.cpp b/coreneuron/nrnoc/fast_imem.cpp
index d545cff34..68c36603d 100644
--- a/coreneuron/nrnoc/fast_imem.cpp
+++ b/coreneuron/nrnoc/fast_imem.cpp
@@ -40,74 +40,29 @@ static int fast_imem_nthread_ = 0;
 static int *fast_imem_size_ = nullptr;
 static NrnFastImem* fast_imem_;
 
-static void fast_imem_free() {
-    int i;
-    for (i = 0; i < nrn_nthread; ++i) {
-        nrn_threads[i].nrn_fast_imem = NULL;
-    }
-    for (i = 0; i < fast_imem_nthread_; ++i) {
-        if (fast_imem_size_[i] > 0) {
-            free(fast_imem_[i].nrn_sav_rhs);
-            free(fast_imem_[i].nrn_sav_d);
-        }
-    }
-    if (fast_imem_nthread_) {
-        free(fast_imem_size_);
-        free(fast_imem_);
-        fast_imem_nthread_ = 0;
-        fast_imem_size_ = nullptr;
-        fast_imem_ = nullptr;
-    }
-}
-
-/*
-Avoid invalidating pointers to i_membrane_ unless the number of compartments
-in a thread has changed.
-*/
-
-static void fast_imem_alloc() {
-    int i;
-    if (fast_imem_nthread_ != nrn_nthread) {
-        fast_imem_free();
-        fast_imem_nthread_ = nrn_nthread;
-        fast_imem_size_ = (int*)ecalloc(nrn_nthread, sizeof(int));
-        fast_imem_ = (NrnFastImem*)ecalloc(nrn_nthread, sizeof(NrnFastImem));
-    }
-    for (i=0; i < nrn_nthread; ++i) {
-        NrnThread* nt = nrn_threads + i;
-        int n = nt->end;
-        NrnFastImem* fi = fast_imem_ + i;
-        if (n != fast_imem_size_[i]) {
-            if (fast_imem_size_[i] > 0) {
-                free(fi->nrn_sav_rhs);
-                free(fi->nrn_sav_d);
-            }
-            if (n > 0) {
-                fi->nrn_sav_rhs = (double*)emalloc_align(n * sizeof(double));
-                fi->nrn_sav_d = (double*)emalloc_align(n * sizeof(double));
-            }
-            fast_imem_size_[i] = n;
+void fast_imem_free() {
+    for (int i = 0; i < nrn_nthread; ++i) {
+        if (nrn_threads[i].nrn_fast_imem) {
+            free(nrn_threads[i].nrn_fast_imem->nrn_sav_rhs);
+            free(nrn_threads[i].nrn_fast_imem->nrn_sav_d);
+            free(nrn_threads[i].nrn_fast_imem);
+            nrn_threads[i].nrn_fast_imem = nullptr;
         }
     }
 }
 
 void nrn_fast_imem_alloc() {
     if (nrn_use_fast_imem) {
-        int i;
-        fast_imem_alloc();
-        for (i=0; i < nrn_nthread; ++i) {
-            nrn_threads[i].nrn_fast_imem = fast_imem_ + i;
-        }
-    }else{
         fast_imem_free();
+        for (int i=0; i < nrn_nthread; ++i) {
+            int n = nrn_threads[i].end;
+            nrn_threads[i].nrn_fast_imem = (NrnFastImem*)ecalloc(1, sizeof(NrnFastImem));
+            nrn_threads[i].nrn_fast_imem->nrn_sav_rhs = (double*)emalloc_align(n * sizeof(double));
+            nrn_threads[i].nrn_fast_imem->nrn_sav_d = (double*)emalloc_align(n * sizeof(double));
+        }
     }
 }
 
-void use_fast_imem() {
-    nrn_use_fast_imem = true;
-    nrn_fast_imem_alloc();
-}
-
 void nrn_calc_fast_imem(NrnThread* _nt) {
     int i;
     int i1 = 0;
diff --git a/coreneuron/nrnoc/fast_imem.h b/coreneuron/nrnoc/fast_imem.h
index 4e54c1512..449021ecb 100644
--- a/coreneuron/nrnoc/fast_imem.h
+++ b/coreneuron/nrnoc/fast_imem.h
@@ -33,16 +33,15 @@ THE POSSIBILITY OF SUCH DAMAGE.
 
 namespace coreneuron {
 
-/* Enables fast membrane curent culculation and allocates required
- * memory.
- * Found in src/nrncvode/cvodeobj.cpp in NEURON.
+/* Bool global variable to define if the fast_imem
+ * calculations should be enabled.
  */
-void use_fast_imem();
+extern bool nrn_use_fast_imem;
 
 /* Free memory allocated for the fast current membrane calculation.
  * Found in src/nrnoc/multicore.c in NEURON.
  */
-static void fast_imem_free();
+void fast_imem_free();
 
 /* Allocate memory for the rhs and d arrays needed for the fast
  * current membrane calculation.
diff --git a/coreneuron/nrnoc/finitialize.cpp b/coreneuron/nrnoc/finitialize.cpp
index 448ff66bf..9f2deb8bf 100644
--- a/coreneuron/nrnoc/finitialize.cpp
+++ b/coreneuron/nrnoc/finitialize.cpp
@@ -34,8 +34,6 @@ THE POSSIBILITY OF SUCH DAMAGE.
 
 namespace coreneuron {
 
-extern bool nrn_use_fast_imem;
-
 void nrn_finitialize(int setv, double v) {
     int i;
     NrnThread* _nt;

From c07b04bf261ed221cebc0eb1d0ecdc1cd57c6709 Mon Sep 17 00:00:00 2001
From: Ioannis Magkanaris <iomagkanaris@gmail.com>
Date: Fri, 20 Sep 2019 15:43:14 +0200
Subject: [PATCH 07/14] Renamed i_membrane enum to i_membrane_

---
 coreneuron/nrniv/nrn_setup.cpp | 4 ++--
 coreneuron/nrniv/nrniv_decl.h  | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/coreneuron/nrniv/nrn_setup.cpp b/coreneuron/nrniv/nrn_setup.cpp
index be9c835b5..2d749fe98 100644
--- a/coreneuron/nrniv/nrn_setup.cpp
+++ b/coreneuron/nrniv/nrn_setup.cpp
@@ -919,7 +919,7 @@ int nrn_i_layout(int icnt, int cnt, int isz, int sz, int layout) {
 }
 
 // take into account alignment, layout, permutation
-// only voltage, i_membrane or mechanism data index allowed. (mtype 0 means time)
+// only voltage, i_membrane_ or mechanism data index allowed. (mtype 0 means time)
 double* stdindex2ptr(int mtype, int index, NrnThread& nt) {
     if (mtype == voltage) {  // voltage
         int v0 = nt._actual_v - nt._data;
@@ -929,7 +929,7 @@ double* stdindex2ptr(int mtype, int index, NrnThread& nt) {
             node_permute(&ix, 1, nt._permute);
         }
         return nt._data + (v0 + ix);                // relative to nt._data
-    } else if (mtype == i_membrane) {  // membrane current from fast_imem calculation
+    } else if (mtype == i_membrane_) {  // membrane current from fast_imem calculation
             int i_mem = nt.nrn_fast_imem->nrn_sav_rhs - nt._data;
             int ix = index;  // relative to _actual_v
             nrn_assert((ix >= 0) && (ix < nt.end));
diff --git a/coreneuron/nrniv/nrniv_decl.h b/coreneuron/nrniv/nrniv_decl.h
index 5be1c1b9b..e2d491218 100644
--- a/coreneuron/nrniv/nrniv_decl.h
+++ b/coreneuron/nrniv/nrniv_decl.h
@@ -37,7 +37,7 @@ THE POSSIBILITY OF SUCH DAMAGE.
 namespace coreneuron {
 
 /// Mechanism type to be used from stdindex2ptr and nrn_dblpntr2nrncore (in Neuron)
-enum mech_type {voltage = -1, i_membrane = -2};
+enum mech_type {voltage = -1, i_membrane_ = -2};
 
 extern int cvode_active_;
 /// Vector of maps for negative presyns

From 9abff72e249b2e872a372fea91e6c44dc0b6a20e Mon Sep 17 00:00:00 2001
From: Ioannis Magkanaris <iomagkanaris@gmail.com>
Date: Mon, 23 Sep 2019 15:29:21 +0200
Subject: [PATCH 08/14] Enable fast_imem calculation through report parsing and
 generate reports for the fast_imem

---
 coreneuron/utils/reports/nrnreport.cpp        | 20 +++++++++++--------
 coreneuron/utils/reports/nrnreport.h          |  2 +-
 .../reports/report_configuration_parser.cpp   |  6 +++++-
 external/mod2c                                |  2 +-
 4 files changed, 19 insertions(+), 11 deletions(-)

diff --git a/coreneuron/utils/reports/nrnreport.cpp b/coreneuron/utils/reports/nrnreport.cpp
index 93568eb7c..a85396fb2 100644
--- a/coreneuron/utils/reports/nrnreport.cpp
+++ b/coreneuron/utils/reports/nrnreport.cpp
@@ -99,7 +99,7 @@ class ReportEvent : public DiscreteEvent {
     }
 };
 
-VarsToReport get_soma_vars_to_report(NrnThread& nt, std::set<int>& target) {
+VarsToReport get_soma_vars_to_report(NrnThread& nt, std::set<int>& target, double* report_variable) {
     VarsToReport vars_to_report;
     NrnThreadMappingInfo* mapinfo = (NrnThreadMappingInfo*)nt.mapping;
 
@@ -123,15 +123,15 @@ VarsToReport get_soma_vars_to_report(NrnThread& nt, std::set<int>& target) {
             /** 1st key is section-id and 1st value is segment of soma */
             int section_id = s->secmap.begin()->first;
             int idx = s->secmap.begin()->second.front();
-            double* v = nt._actual_v + idx;
-            to_report.push_back(VarWithMapping(section_id, v));
+            double* variable = report_variable + idx;
+            to_report.push_back(VarWithMapping(section_id, variable));
             vars_to_report[gid] = to_report;
         }
     }
     return vars_to_report;
 }
 
-VarsToReport get_compartment_vars_to_report(NrnThread& nt, std::set<int>& target) {
+VarsToReport get_compartment_vars_to_report(NrnThread& nt, std::set<int>& target, double* report_variable) {
     VarsToReport vars_to_report;
     NrnThreadMappingInfo* mapinfo = (NrnThreadMappingInfo*)nt.mapping;
     if (!mapinfo) {
@@ -160,8 +160,8 @@ VarsToReport get_compartment_vars_to_report(NrnThread& nt, std::set<int>& target
                     for (size_t k = 0; k < vec.size(); k++) {
                         int idx = vec[k];
                         /** corresponding voltage in coreneuron voltage array */
-                        double* v = nt._actual_v + idx;
-                        to_report.push_back(VarWithMapping(compartment_id, v));
+                        double* variable = report_variable + idx;
+                        to_report.push_back(VarWithMapping(compartment_id, variable));
                     }
                 }
             }
@@ -413,11 +413,15 @@ void register_report(double dt, double tstop, double delay, ReportConfiguration&
         VarsToReport vars_to_report;
         switch (report.type) {
             case SomaReport:
-                vars_to_report = get_soma_vars_to_report(nt, report.target);
+                vars_to_report = get_soma_vars_to_report(nt, report.target, nt._actual_v);
                 register_soma_report(nt, report, vars_to_report);
                 break;
             case CompartmentReport:
-                vars_to_report = get_compartment_vars_to_report(nt, report.target);
+                vars_to_report = get_compartment_vars_to_report(nt, report.target, nt._actual_v);
+                register_compartment_report(nt, report, vars_to_report);
+                break;
+            case IMembraneReport:
+                vars_to_report = get_compartment_vars_to_report(nt, report.target, nt.nrn_fast_imem->nrn_sav_rhs);
                 register_compartment_report(nt, report, vars_to_report);
                 break;
             default:
diff --git a/coreneuron/utils/reports/nrnreport.h b/coreneuron/utils/reports/nrnreport.h
index c0c523779..3cef03bb7 100644
--- a/coreneuron/utils/reports/nrnreport.h
+++ b/coreneuron/utils/reports/nrnreport.h
@@ -51,7 +51,7 @@ namespace coreneuron {
 /// name of the variable in mod file used for setting synapse id
 #define SYNAPSE_ID_MOD_NAME "synapseID"
 
-enum ReportType { SomaReport, CompartmentReport, SynapseReport };
+enum ReportType { SomaReport, CompartmentReport, SynapseReport, IMembraneReport };
 
 struct ReportConfiguration {
     char name[MAX_REPORT_NAME_LEN];         // name of the report
diff --git a/coreneuron/utils/reports/report_configuration_parser.cpp b/coreneuron/utils/reports/report_configuration_parser.cpp
index 901469e52..984019839 100644
--- a/coreneuron/utils/reports/report_configuration_parser.cpp
+++ b/coreneuron/utils/reports/report_configuration_parser.cpp
@@ -28,6 +28,7 @@
 
 #include "coreneuron/nrniv/nrn_assert.h"
 #include "coreneuron/utils/reports/nrnreport.h"
+#include "coreneuron/nrnoc/fast_imem.h"
 #include "coreneuron/nrnoc/mech_mapping.hpp"
 #include <stdio.h>
 #include <stdlib.h>
@@ -90,7 +91,10 @@ std::vector<ReportConfiguration> create_report_configurations(const char* conf_f
         }
         sprintf(report.output_path, "%s/%s", output_dir, report.name);
         if (strcmp(report.type_str, "compartment") == 0) {
-            if (is_soma)
+            if (strcmp(report_on, "i_membrane") == 0) {
+                nrn_use_fast_imem = true;
+                report.type = IMembraneReport;
+            } else if (is_soma)
                 report.type = SomaReport;
             else
                 report.type = CompartmentReport;
diff --git a/external/mod2c b/external/mod2c
index 788af0ef1..20a743613 160000
--- a/external/mod2c
+++ b/external/mod2c
@@ -1 +1 @@
-Subproject commit 788af0ef16821bb51e9b80b5140f5d90a780070a
+Subproject commit 20a74361318147ba032c826d92996281883e103a

From e96273555bd3164448322120713d0b7577579501 Mon Sep 17 00:00:00 2001
From: Ioannis Magkanaris <iomagkanaris@gmail.com>
Date: Wed, 25 Sep 2019 18:25:14 +0200
Subject: [PATCH 09/14] Added fast_imem calculation to direct memory test

---
 tests/jenkins/neuron_direct.py | 12 ++++++++++--
 1 file changed, 10 insertions(+), 2 deletions(-)

diff --git a/tests/jenkins/neuron_direct.py b/tests/jenkins/neuron_direct.py
index a47f072b8..e43b3a964 100644
--- a/tests/jenkins/neuron_direct.py
+++ b/tests/jenkins/neuron_direct.py
@@ -10,17 +10,22 @@
 ic.dur = 0.1
 ic.amp = 0.3
 
+h.cvode.use_fast_imem(1)
+h.cvode.cache_efficient(1)
+
 v = h.Vector()
 v.record(h.soma(.5)._ref_v, sec = h.soma)
+i_mem = h.Vector()
+i_mem.record(h.soma(.5)._ref_i_membrane_, sec = h.soma)
 tv = h.Vector()
 tv.record(h._ref_t, sec=h.soma)
 h.run()
 vstd = v.cl()
 tvstd = tv.cl()
+i_memstd = i_mem.cl()
 
 #h.CoreNeuronRun[0].run()
 pc = h.ParallelContext()
-h.cvode.cache_efficient(1)
 h.stdinit()
 pc.nrncore_run("-e %g"%h.tstop, 0)
 
@@ -30,6 +35,9 @@
 if v.cl().sub(vstd).abs().max() >= 1e-10:
     print("Voltage difference greater than or equal to 1e-10")
     sys.exit(-1)
+if i_mem.cl().sub(i_memstd).abs().max() >= 1e-10:
+    print("i_mem difference greater than or equal to 1e-10")
+    sys.exit(-1)
 
-print("Voltage times are same and difference is less than 1e-10")
+print("Voltage times and i_membrane_ are same and difference is less than 1e-10")
 h.quit()

From 2304e9a30767ddaf70d4a0e6837090e511d30137 Mon Sep 17 00:00:00 2001
From: Ioannis Magkanaris <iomagkanaris@gmail.com>
Date: Thu, 26 Sep 2019 13:50:02 +0200
Subject: [PATCH 10/14] Added small explanation about the mech_type enum

---
 coreneuron/nrniv/nrniv_decl.h          |  2 ++
 coreneuron/nrniv/tqueue.h              |  2 +-
 coreneuron/nrnoc/fast_imem.cpp         | 16 ++++++++++++++++
 coreneuron/utils/reports/nrnreport.cpp |  1 +
 4 files changed, 20 insertions(+), 1 deletion(-)

diff --git a/coreneuron/nrniv/nrniv_decl.h b/coreneuron/nrniv/nrniv_decl.h
index e2d491218..578d0045d 100644
--- a/coreneuron/nrniv/nrniv_decl.h
+++ b/coreneuron/nrniv/nrniv_decl.h
@@ -37,6 +37,8 @@ THE POSSIBILITY OF SUCH DAMAGE.
 namespace coreneuron {
 
 /// Mechanism type to be used from stdindex2ptr and nrn_dblpntr2nrncore (in Neuron)
+/// Values of the mechanism types should be negative numbers to avoid any conflict with
+/// mechanism types of Memb_list(>0) or time(0) passed from Neuron
 enum mech_type {voltage = -1, i_membrane_ = -2};
 
 extern int cvode_active_;
diff --git a/coreneuron/nrniv/tqueue.h b/coreneuron/nrniv/tqueue.h
index 0271d564e..baef104bb 100644
--- a/coreneuron/nrniv/tqueue.h
+++ b/coreneuron/nrniv/tqueue.h
@@ -60,7 +60,7 @@ namespace coreneuron {
 #define STRCMP(a, b) (a - b)
 
 class TQItem;
-#define SPBLK TQItem
+#define SPBLK SPC
 #define leftlink left_
 #define rightlink right_
 #define uplink parent_
diff --git a/coreneuron/nrnoc/fast_imem.cpp b/coreneuron/nrnoc/fast_imem.cpp
index 68c36603d..984e8830e 100644
--- a/coreneuron/nrnoc/fast_imem.cpp
+++ b/coreneuron/nrnoc/fast_imem.cpp
@@ -73,9 +73,25 @@ void nrn_calc_fast_imem(NrnThread* _nt) {
 
     double* pd = _nt->nrn_fast_imem->nrn_sav_d;
     double* prhs = _nt->nrn_fast_imem->nrn_sav_rhs;
+    FILE *fp_rhs, *fp_d;
+    char rhs_filename[20], d_filename[20];
+    sprintf(rhs_filename,"rhs.CORENEURON.%d",nrnmpi_myid);
+    sprintf(d_filename,"d.CORENEURON.%d",nrnmpi_myid);
+    fp_rhs = fopen(rhs_filename, "a");
+    fp_d = fopen(d_filename, "a");
+    fprintf(fp_rhs, "\n%.8e time\n", _nt->_t);
+    fprintf(fp_d, "\n%.8e time\n", _nt->_t);
     for (i = i1; i < i3 ; ++i) {
         prhs[i] = (pd[i]*vec_rhs[i] + prhs[i])*vec_area[i]*0.01;
+        fprintf(fp_rhs, "%.8e, ", prhs[i]);
     }
+    fprintf(fp_rhs, "\n");
+    for (i = i1; i < i3 ; ++i) {
+        fprintf(fp_d, "%.8e, ", pd[i]);
+    }
+    fprintf(fp_d, "\n");
+    fclose(fp_rhs);
+    fclose(fp_d);
 }
 
 }
diff --git a/coreneuron/utils/reports/nrnreport.cpp b/coreneuron/utils/reports/nrnreport.cpp
index a85396fb2..9b6e6ff9f 100644
--- a/coreneuron/utils/reports/nrnreport.cpp
+++ b/coreneuron/utils/reports/nrnreport.cpp
@@ -88,6 +88,7 @@ class ReportEvent : public DiscreteEvent {
 #pragma omp critical
         {
             // each thread needs to know its own step
+            printf("Real time: %lf\n",t);
             records_nrec(step, gids_to_report.size(), &gids_to_report[0], report_path);
             send(t + dt, nc, nt);
             step++;

From f0d2c3e36027af44576c34deb534e17cc91148cc Mon Sep 17 00:00:00 2001
From: Ioannis Magkanaris <iomagkanaris@gmail.com>
Date: Thu, 26 Sep 2019 13:52:42 +0200
Subject: [PATCH 11/14] Revert "Added small explanation about the mech_type
 enum"

This reverts commit 2304e9a30767ddaf70d4a0e6837090e511d30137 and adds
only the proper files
---
 coreneuron/nrniv/tqueue.h              |  2 +-
 coreneuron/nrnoc/fast_imem.cpp         | 16 ----------------
 coreneuron/utils/reports/nrnreport.cpp |  1 -
 3 files changed, 1 insertion(+), 18 deletions(-)

diff --git a/coreneuron/nrniv/tqueue.h b/coreneuron/nrniv/tqueue.h
index baef104bb..0271d564e 100644
--- a/coreneuron/nrniv/tqueue.h
+++ b/coreneuron/nrniv/tqueue.h
@@ -60,7 +60,7 @@ namespace coreneuron {
 #define STRCMP(a, b) (a - b)
 
 class TQItem;
-#define SPBLK SPC
+#define SPBLK TQItem
 #define leftlink left_
 #define rightlink right_
 #define uplink parent_
diff --git a/coreneuron/nrnoc/fast_imem.cpp b/coreneuron/nrnoc/fast_imem.cpp
index 984e8830e..68c36603d 100644
--- a/coreneuron/nrnoc/fast_imem.cpp
+++ b/coreneuron/nrnoc/fast_imem.cpp
@@ -73,25 +73,9 @@ void nrn_calc_fast_imem(NrnThread* _nt) {
 
     double* pd = _nt->nrn_fast_imem->nrn_sav_d;
     double* prhs = _nt->nrn_fast_imem->nrn_sav_rhs;
-    FILE *fp_rhs, *fp_d;
-    char rhs_filename[20], d_filename[20];
-    sprintf(rhs_filename,"rhs.CORENEURON.%d",nrnmpi_myid);
-    sprintf(d_filename,"d.CORENEURON.%d",nrnmpi_myid);
-    fp_rhs = fopen(rhs_filename, "a");
-    fp_d = fopen(d_filename, "a");
-    fprintf(fp_rhs, "\n%.8e time\n", _nt->_t);
-    fprintf(fp_d, "\n%.8e time\n", _nt->_t);
     for (i = i1; i < i3 ; ++i) {
         prhs[i] = (pd[i]*vec_rhs[i] + prhs[i])*vec_area[i]*0.01;
-        fprintf(fp_rhs, "%.8e, ", prhs[i]);
     }
-    fprintf(fp_rhs, "\n");
-    for (i = i1; i < i3 ; ++i) {
-        fprintf(fp_d, "%.8e, ", pd[i]);
-    }
-    fprintf(fp_d, "\n");
-    fclose(fp_rhs);
-    fclose(fp_d);
 }
 
 }
diff --git a/coreneuron/utils/reports/nrnreport.cpp b/coreneuron/utils/reports/nrnreport.cpp
index 9b6e6ff9f..a85396fb2 100644
--- a/coreneuron/utils/reports/nrnreport.cpp
+++ b/coreneuron/utils/reports/nrnreport.cpp
@@ -88,7 +88,6 @@ class ReportEvent : public DiscreteEvent {
 #pragma omp critical
         {
             // each thread needs to know its own step
-            printf("Real time: %lf\n",t);
             records_nrec(step, gids_to_report.size(), &gids_to_report[0], report_path);
             send(t + dt, nc, nt);
             step++;

From 36bcbb5501e0ea0462c7c0407d2293526fd459d5 Mon Sep 17 00:00:00 2001
From: Ioannis Magkanaris <iomagkanaris@gmail.com>
Date: Thu, 26 Sep 2019 17:24:22 +0200
Subject: [PATCH 12/14] Addressesd review comments

- Added documentation in stdindex2ptr function
- Changed fast_imem (de)allocation functions with "CoreNeuron" for loop
- Small identation fix
---
 coreneuron/nrniv/nrn_setup.cpp |  7 +++++--
 coreneuron/nrnoc/fast_imem.cpp | 22 +++++++++++-----------
 2 files changed, 16 insertions(+), 13 deletions(-)

diff --git a/coreneuron/nrniv/nrn_setup.cpp b/coreneuron/nrniv/nrn_setup.cpp
index 2d749fe98..0eb3c4a9b 100644
--- a/coreneuron/nrniv/nrn_setup.cpp
+++ b/coreneuron/nrniv/nrn_setup.cpp
@@ -918,6 +918,9 @@ int nrn_i_layout(int icnt, int cnt, int isz, int sz, int layout) {
     return 0;
 }
 
+// This function is related to nrn_dblpntr2nrncore in Neuron to determine which values should
+// be transferred from CoreNeuron. Types correspond to the value to be transferred based on
+// mech_type enum or non-artificial cell mechanisms.
 // take into account alignment, layout, permutation
 // only voltage, i_membrane_ or mechanism data index allowed. (mtype 0 means time)
 double* stdindex2ptr(int mtype, int index, NrnThread& nt) {
@@ -929,9 +932,9 @@ double* stdindex2ptr(int mtype, int index, NrnThread& nt) {
             node_permute(&ix, 1, nt._permute);
         }
         return nt._data + (v0 + ix);                // relative to nt._data
-    } else if (mtype == i_membrane_) {  // membrane current from fast_imem calculation
+    } else if (mtype == i_membrane_) {              // membrane current from fast_imem calculation
             int i_mem = nt.nrn_fast_imem->nrn_sav_rhs - nt._data;
-            int ix = index;  // relative to _actual_v
+            int ix = index;  // relative to nrn_fast_imem->nrn_sav_rhs
             nrn_assert((ix >= 0) && (ix < nt.end));
             if (nt._permute) {
                 node_permute(&ix, 1, nt._permute);
diff --git a/coreneuron/nrnoc/fast_imem.cpp b/coreneuron/nrnoc/fast_imem.cpp
index 68c36603d..141eee889 100644
--- a/coreneuron/nrnoc/fast_imem.cpp
+++ b/coreneuron/nrnoc/fast_imem.cpp
@@ -41,12 +41,12 @@ static int *fast_imem_size_ = nullptr;
 static NrnFastImem* fast_imem_;
 
 void fast_imem_free() {
-    for (int i = 0; i < nrn_nthread; ++i) {
-        if (nrn_threads[i].nrn_fast_imem) {
-            free(nrn_threads[i].nrn_fast_imem->nrn_sav_rhs);
-            free(nrn_threads[i].nrn_fast_imem->nrn_sav_d);
-            free(nrn_threads[i].nrn_fast_imem);
-            nrn_threads[i].nrn_fast_imem = nullptr;
+    for (NrnThread* nt = nrn_threads; nt < nrn_threads + nrn_nthread; ++nt) {
+        if (nt->nrn_fast_imem) {
+            free(nt->nrn_fast_imem->nrn_sav_rhs);
+            free(nt->nrn_fast_imem->nrn_sav_d);
+            free(nt->nrn_fast_imem);
+            nt->nrn_fast_imem = nullptr;
         }
     }
 }
@@ -54,11 +54,11 @@ void fast_imem_free() {
 void nrn_fast_imem_alloc() {
     if (nrn_use_fast_imem) {
         fast_imem_free();
-        for (int i=0; i < nrn_nthread; ++i) {
-            int n = nrn_threads[i].end;
-            nrn_threads[i].nrn_fast_imem = (NrnFastImem*)ecalloc(1, sizeof(NrnFastImem));
-            nrn_threads[i].nrn_fast_imem->nrn_sav_rhs = (double*)emalloc_align(n * sizeof(double));
-            nrn_threads[i].nrn_fast_imem->nrn_sav_d = (double*)emalloc_align(n * sizeof(double));
+        for (NrnThread* nt = nrn_threads; nt < nrn_threads + nrn_nthread; ++nt) {
+            int n = nt->end;
+            nt->nrn_fast_imem = (NrnFastImem*)ecalloc(1, sizeof(NrnFastImem));
+            nt->nrn_fast_imem->nrn_sav_rhs = (double*)emalloc_align(n * sizeof(double));
+            nt->nrn_fast_imem->nrn_sav_d = (double*)emalloc_align(n * sizeof(double));
         }
     }
 }

From b3c7b03ea4287096685a7cef7339acb9cf188ea9 Mon Sep 17 00:00:00 2001
From: Ioannis Magkanaris <iomagkanaris@gmail.com>
Date: Mon, 30 Sep 2019 10:59:27 +0200
Subject: [PATCH 13/14] Addressed review comments

- Removed unused variables
- Improvements on used variables
- Removed change in mod2c submodule
---
 coreneuron/nrnoc/fast_imem.cpp | 20 ++++++++------------
 coreneuron/nrnoc/multicore.h   |  4 ++--
 external/mod2c                 |  2 +-
 3 files changed, 11 insertions(+), 15 deletions(-)

diff --git a/coreneuron/nrnoc/fast_imem.cpp b/coreneuron/nrnoc/fast_imem.cpp
index 141eee889..99614cb42 100644
--- a/coreneuron/nrnoc/fast_imem.cpp
+++ b/coreneuron/nrnoc/fast_imem.cpp
@@ -36,9 +36,6 @@ namespace coreneuron {
 extern int nrn_nthread;
 extern NrnThread *nrn_threads;
 bool nrn_use_fast_imem;
-static int fast_imem_nthread_ = 0;
-static int *fast_imem_size_ = nullptr;
-static NrnFastImem* fast_imem_;
 
 void fast_imem_free() {
     for (NrnThread* nt = nrn_threads; nt < nrn_threads + nrn_nthread; ++nt) {
@@ -63,18 +60,17 @@ void nrn_fast_imem_alloc() {
     }
 }
 
-void nrn_calc_fast_imem(NrnThread* _nt) {
-    int i;
+void nrn_calc_fast_imem(NrnThread* nt) {
     int i1 = 0;
-    int i3 = _nt->end;
+    int i3 = nt->end;
 
-    double* vec_rhs = &(VEC_RHS(0));
-    double* vec_area = &(VEC_AREA(0));
+    double* vec_rhs = nt->_actual_rhs;
+    double* vec_area = nt->_actual_area;
 
-    double* pd = _nt->nrn_fast_imem->nrn_sav_d;
-    double* prhs = _nt->nrn_fast_imem->nrn_sav_rhs;
-    for (i = i1; i < i3 ; ++i) {
-        prhs[i] = (pd[i]*vec_rhs[i] + prhs[i])*vec_area[i]*0.01;
+    double* fast_imem_d = nt->nrn_fast_imem->nrn_sav_d;
+    double* fast_imem_rhs = nt->nrn_fast_imem->nrn_sav_rhs;
+    for (int i = i1; i < i3 ; ++i) {
+        fast_imem_rhs[i] = (fast_imem_d[i]*vec_rhs[i] + fast_imem_rhs[i])*vec_area[i]*0.01;
     }
 }
 
diff --git a/coreneuron/nrnoc/multicore.h b/coreneuron/nrnoc/multicore.h
index 9aaec8bf0..9d879fbdf 100644
--- a/coreneuron/nrnoc/multicore.h
+++ b/coreneuron/nrnoc/multicore.h
@@ -59,10 +59,10 @@ struct NrnThreadBAList {
     NrnThreadBAList* next;
 };
 
-typedef struct NrnFastImem {
+struct NrnFastImem {
     double* nrn_sav_rhs;
     double* nrn_sav_d;
-} NrnFastImem;
+};
 
 struct TrajectoryRequests {
     void** vpr;       /* PlayRecord Objects known by NEURON */
diff --git a/external/mod2c b/external/mod2c
index 20a743613..788af0ef1 160000
--- a/external/mod2c
+++ b/external/mod2c
@@ -1 +1 @@
-Subproject commit 20a74361318147ba032c826d92996281883e103a
+Subproject commit 788af0ef16821bb51e9b80b5140f5d90a780070a

From 6b280f10b03c4d7a5647bd196d83c2dc0c25c00f Mon Sep 17 00:00:00 2001
From: Ioannis Magkanaris <iomagkanaris@gmail.com>
Date: Mon, 30 Sep 2019 13:16:57 +0200
Subject: [PATCH 14/14] Updated mod2c submodule commit

---
 external/mod2c | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/external/mod2c b/external/mod2c
index 788af0ef1..0a0bee5e7 160000
--- a/external/mod2c
+++ b/external/mod2c
@@ -1 +1 @@
-Subproject commit 788af0ef16821bb51e9b80b5140f5d90a780070a
+Subproject commit 0a0bee5e79fbb533e8bb766c7e0efeb3b6daafea