text: #include <../../nrnconf.h>
  
  #include <nrnmpi.h>
  #include <stdlib.h>
  #include <errno.h>
  #include "neuron.h"
  #include "section.h"
  #include "nrniv_mf.h"
  #include "multisplit.h"
  #define nrnoc_fadvance_c
  #include "utils/profile/profiler_interface.h"
  #include "nonvintblock.h"
  #include "nrncvode.h"
  #include "spmatrix.h"
  #include <vector>
  
  /*
   after an fadvance from t-dt to t, v is defined at t
   states that depend on v are defined at t+dt/2
   states, such as concentration that depend on current are defined at t
   and currents are defined at t-dt/2.
   */
  /*
   fcurrent is used to set up all assigned variables  without changing any
   states.  It assumes all states have their values at time t which is
   only first order correct. It determines the nonvint assignments first
   so that ena, etc. are correct for the current determinations. We
   demand that nonvint assignments can be done without changing states when
   dt is temporarily set to 0.
  
   It turned out to be a bad idea to do this. Many methods (KINETIC, PROCEDURE)
   for solving are not even called when dt =0. Also it plays havoc with tables
   that are recomputed several times when dt is changed then changed back.
   Therefore we are no longer trying to initialize assigned variables within
   SOLVE'd blocks. Instead we are making initialization more convenient,
   at least within the normal situations, by introducing a finitialize()
   with an optional argument vinit. finitialize will call the INITIALIZE block
   for all mechanisms in all segments. (Default initialization sets all states
   to state0.) The user can set things up specially in a models INITIAL
   block. INITIAL blocks can make use of v. With care they can make use of
   ionic concentrations just like breakpoint and solve blocks. When the argument
   is present, v for all segments are set to that value.
   */
  
  #if !defined(NRNMPI) || NRNMPI == 0
  extern double nrnmpi_wtime();
  #endif
  extern double* nrn_mech_wtime_;
  extern double t, dt;
  extern double chkarg(int, double low, double high);
  extern void nrn_fixed_step();
  extern void nrn_fixed_step_group(int);
  static void* nrn_fixed_step_thread(NrnThread*);
  static void* nrn_fixed_step_group_thread(NrnThread* nth);
  static void* nrn_fixed_step_lastpart(NrnThread*);
  extern void* setup_tree_matrix(NrnThread*);
  extern void nrn_solve(NrnThread*);
  extern void nonvint(NrnThread* nt);
  extern void nrncvode_set_t(double t);
  
  static void* nrn_ms_treeset_through_triang(NrnThread*);
  static void* nrn_ms_reduce_solve(NrnThread*);
  static void* nrn_ms_bksub(NrnThread*);
  static void* nrn_ms_bksub_through_triang(NrnThread*);
  extern void* nrn_multisplit_triang(NrnThread*);
  extern void* nrn_multisplit_reduce_solve(NrnThread*);
  extern void* nrn_multisplit_bksub(NrnThread*);
  extern void (*nrn_multisplit_setup_)();
  void (*nrn_allthread_handle)();
  
  extern int tree_changed;
  extern int diam_changed;
  extern int state_discon_allowed_;
  extern double hoc_epsilon;
  
  static void update(NrnThread*);
  
  #if 0
  /*
    1 to save space, but must worry about other uses of the memb_list
    such as in netcvode.cpp in fornetcon_prepare
  */
  extern short* nrn_is_artificial_;
  extern cTemplate** nrn_pnt_template_;
  #endif
  
  #define NRNCTIME 1
  #define NONVINT_ODE_COUNT 5
  
  #if NRNCTIME
  #define CTBEGIN double wt = nrnmpi_wtime();
  #define CTADD nth->_ctime += nrnmpi_wtime() - wt;
  #else
  #define CTBEGIN /**/
  #define CTADD /**/
  #endif
  
  #define ELIMINATE_T_ROUNDOFF 0
  #if ELIMINATE_T_ROUNDOFF
  /* in order to simplify and as much as possible avoid the handling
  of round-off error due to repeated t += dt, we use
  t = nrn_tbase_ + nrn_ndt_ * nrn_dt_;
  to advance time. The problems that must be overcome are when the
  user changes dt in the middle of a run or switches from cvode or
  abruptly and arbitrarily sets a new value of t and continues from
  there (hence nrn_tbase_)
  */
  double nrn_ndt_, nrn_tbase_, nrn_dt_;
  void nrn_chk_ndt() {
  	if (dt != nrn_dt_ || t != nrn_tbase_ + nrn_ndt_ * nrn_dt_) {
  if (nrnmpi_myid == 0) Printf("nrn_chk_ndt t=%g dt=%g old nrn_tbase_=%g nrn_ndt_=%g nrn_dt_=%g\n",
  t, dt, nrn_tbase_, nrn_ndt_, nrn_dt_);
  		nrn_dt_ = dt;
  		nrn_tbase_ = t;
  		nrn_ndt_ = 0.;
  	}
  }
  #endif /* ELIMINATE_T_ROUNDOFF */
  
  /*
  There are (too) many variants of nrn_fixed_step depending on
  nrnmpi_numprocs 1 or > 1, nrn_nthread 1 or > 1,
  nrnmpi_v_transfer nil or callable, nrn_multisplit_setup nil or callable,
  and whether one step with fadvance
  or possibly many with ParallelContext.psolve before synchronizing with
  NetParEvent. The combination of simultaneous nrnmpi_numprocs > 1 and
  nrn_nthread > 1 with parallel transfer
  requires some refactoring of the use of the old
  (*nrnmpi_v_transfer_)() from within nonvint(nt) which handled either mpi
  or threads but cannot handle both simultaneously. (Unless the first
  thread that arrives in nrnmpi_v_transfer is the one that accomplishes
  the mpi transfer). Instead, we replace with
  nrnthread_v_transfer(nt) to handle the per thread copying of source
  data into the target space owned by the thread. Between update (called by
  nrn_fixed_step_thread and nrn_ms_bksub), and nonvint (called by
  nrn_fixed_step_lastpart, we need to have thread 0 do the interprocessor
  transfer of source voltages to a either a staging area for later
  copying to the target by a thread (or directly to the target if only one
  thread). This will happen with (*nrnmpi_v_transfer_)(). (Look Ma, no threads).
  This deals properly with the necessary mpi synchronization. And leaves
  thread handling where it was before. Also, writing to the staging area
  is only done by thread 0. Fixed step and global variable step
  logic is limited to the case where an nrnmpi_v_transfer requires
  existence of nrnthread_v_transfer (even if one thread).
  */
  #if 1 || PARANEURON
  void (*nrnmpi_v_transfer_)(); /* called by thread 0 */
  void (*nrnthread_v_transfer_)(NrnThread* nt);
  /* if at least one gap junction has a source voltage with extracellular inserted */
  void (*nrnthread_vi_compute_)(NrnThread* nt);
  #endif
  
  #if VECTORIZE
  extern "C" int v_structure_change;
  #endif
  
  #if CVODE
  int cvode_active_;
  #endif
  
  int stoprun;
  int nrn_use_fast_imem;
  
  #define PROFILE 0
  #include "profile.h"
  
  void fadvance(void) {
      nrn::Instrumentor::phase p_fadvance("fadvance");
      tstopunset;
  #if CVODE
  	if (cvode_active_) {
  		cvode_fadvance(-1.);
  		tstopunset;
  		hoc_retpushx(1.);
  		return;
  	}
  #endif
  	if (tree_changed) {
  		setup_topology();
  	}
  	if (v_structure_change) {
  		v_setup_vectors();
  	}
  	if (diam_changed) {
  		recalc_diam();
  	}
  	nrn_fixed_step();
  	tstopunset;
  	hoc_retpushx(1.);
  }
  
  /*
    batch_save() initializes list of variables
    batch_save(&varname, ...) adds variable names to list for saving
    batch_save("varname", ...) adds variable names to the list and name
    				will appear in header.
    batch_run(tstop, tstep, "file") saves variables in file every tstep
  */
  
  static void batch_out(), batch_open(), batch_close();
  
  static FILE* batch_file;
  static int batch_size;
  static int batch_n;
  static double** batch_var;
  
  static void  batch_open(char* name, double tstop, double tstep, char* comment){
      if (batch_file) {
          batch_close();
      }
      if (!name) {
          return;
      }
      batch_file = fopen(name, "w");
      if (!batch_file) {
          hoc_execerror("Couldn't open batch file", name);
      }
      fprintf(batch_file, "%s\nbatch_run from t = %g to %g in steps of %g with dt = %g\n", comment, t, tstop, tstep, dt);
      #if 0
      fprintf(batch_file, "variable names --\n");
          if (!batch_var) {
              batch_var = hoc_newlist();
          }
          count = 0;
          ITERATE(q, batch_varname) {
              fprintf(batch_file, "%s\n", STR(q));
              ++count;
          }
          if (count != batch_n) {
              if (batch_var) {
                  free(batch_var);
              }
              batch_n = count;
              batch_var = (double**)ecalloc(batch_n, sizeof(double*));
          }
          count = 0;
          ITERATE(q, batch_varname) {
              batch_var[count++] = hoc_val_pointer(STR(q));
          }
      #endif
  }
  
  void batch_run(void) /* avoid interpreter overhead */
  {
  	double tstop, tstep, tnext;
  	char* filename;
  	char* comment;
  
  	tstopunset;
  	tstop = chkarg(1,0.,1e20);
  	tstep = chkarg(2, 0., 1e20);
  	if (ifarg(3)) {
  		filename = gargstr(3);
  	}else{
  		filename = 0;
  	}
  	if (ifarg(4)) {
  		comment = gargstr(4);
  	}else{
  		comment = "";
  	}
  	
  	if (tree_changed) {
  		setup_topology();
  	}
  #if VECTORIZE
  	if (v_structure_change) {
  		v_setup_vectors();
  	}
  #endif
  	batch_open(filename, tstop, tstep, comment);
  	batch_out();
  	if (cvode_active_) {
  		while (t < tstop) {
  			cvode_fadvance(t+tstep);
  			batch_out();
  		}
  	}else{
  		tstep -= dt/4.;
  		tstop -= dt/4.;
  		tnext = t + tstep;
  		while (t < tstop) {
  			nrn_fixed_step();
  			if (t > tnext) {
  				batch_out();
  				tnext = t + tstep;
  			}
  			if (stoprun) { tstopunset; break; }
  		}
  	}
  	batch_close();
  	hoc_retpushx(1.);
  }
  
  static void dt2thread(double adt) {
      if (adt != nrn_threads[0]._dt) {
  	int i;
  	for (i=0; i < nrn_nthread; ++i) {
  		NrnThread* nt = nrn_threads + i;
  		nt->_t = t;
  		nt->_dt = dt;
  		if (secondorder) {
  			nt->cj = 2.0/dt;
  		}else{
  			nt->cj = 1.0/dt;
  		}
  	}
      }
  }
  
  static int _upd;
  static void* daspk_init_step_thread(NrnThread* nt) {
  	setup_tree_matrix(nt);
  	nrn_solve(nt);
  	if (_upd) {
  		update(nt);
  	}
  	return nullptr;
  }
  
  void nrn_daspk_init_step(double tt, double dteps, int upd){
  	int i;
  	double dtsav = nrn_threads->_dt;
  	int so = secondorder;
  	dt = dteps;
  	t = tt;
  	secondorder = 0;
  	dt2thread(dteps);
  	nrn_thread_table_check();
  	_upd = upd;
  	nrn_multithread_job(daspk_init_step_thread);
  	dt = dtsav;
  	secondorder = so;
  	dt2thread(dtsav);
  	nrn_thread_table_check();
  }
  
  void nrn_fixed_step() {
      nrn::Instrumentor::phase p_timestep("timestep");
      int i;
  #if ELIMINATE_T_ROUNDOFF
  	nrn_chk_ndt();
  #endif
  	if (t != nrn_threads->_t) {
  		dt2thread(-1.);
  	}else{
  		dt2thread(dt);
  	}
  	nrn_thread_table_check();
  	if (nrn_multisplit_setup_) {
  		nrn_multithread_job(nrn_ms_treeset_through_triang);
  		// remove to avoid possible deadlock where some ranks do a
  		// v transfer and others do a spike exchange.
  		// i.e. must complete the full multisplit time step.
  		//if (!nrn_allthread_handle) {
  			nrn_multithread_job(nrn_ms_reduce_solve);
  			nrn_multithread_job(nrn_ms_bksub);
  			/* see comment below */
  			if (nrnthread_v_transfer_) {
  				if (nrnmpi_v_transfer_) {
                      nrn::Instrumentor::phase p_gap("gap-v-transfer");
                      (*nrnmpi_v_transfer_)();
  				}
  				nrn_multithread_job(nrn_fixed_step_lastpart);
  			}
          //}
  	}else{
  		nrn_multithread_job(nrn_fixed_step_thread);
  /* if there is no nrnthread_v_transfer then there cannot be
     a nrnmpi_v_transfer and lastpart
     will be done in above call.
  */
  		if (nrnthread_v_transfer_) {
  			if (nrnmpi_v_transfer_) {
                  nrn::Instrumentor::phase p_gap("gap-v-transfer");
                  (*nrnmpi_v_transfer_)();
  			}
  			nrn_multithread_job(nrn_fixed_step_lastpart);
  		}
  	}
  	t = nrn_threads[0]._t;
  	if (nrn_allthread_handle) { (*nrn_allthread_handle)(); }
  }
  
  /* better cache efficiency since a thread can do an entire minimum delay
  integration interval before joining
  */
  static int step_group_n;
  static int step_group_begin;
  static int step_group_end;
  
  void nrn_fixed_step_group(int n) {
  	int i;
  #if ELIMINATE_T_ROUNDOFF
  	nrn_chk_ndt();
  #endif
  	dt2thread(dt);
  	nrn_thread_table_check();
  	if (nrn_multisplit_setup_) {
  		int b = 0;
  		nrn_multithread_job(nrn_ms_treeset_through_triang);
  		step_group_n = 0; /* abort at bksub flag */
  		for (i=1; i < n; ++i) {
  			nrn_multithread_job(nrn_ms_reduce_solve);
  			nrn_multithread_job(nrn_ms_bksub_through_triang);
  			if (step_group_n) {
  				step_group_n = 0;
  				if (nrn_allthread_handle) {
  					(*nrn_allthread_handle)();
  				}
  				/* aborted step at bksub, so if not stopped
  				   must do the triang*/
  				b = 1;
  				if (!stoprun) {
  					nrn_multithread_job(nrn_ms_treeset_through_triang);
  				}
  			}
  			if (stoprun) { break; }
  			b = 0;
  		}
  		if (!b) {
  			nrn_multithread_job(nrn_ms_reduce_solve);
  			nrn_multithread_job(nrn_ms_bksub);
  		}
  		if (nrn_allthread_handle) { (*nrn_allthread_handle)(); }
  	}else{
  		step_group_n = n;
  		step_group_begin = 0;
  		step_group_end = 0;
  		while(step_group_end < step_group_n) {
  /*printf("step_group_end=%d step_group_n=%d\n", step_group_end, step_group_n);*/
  			nrn_multithread_job(nrn_fixed_step_group_thread);
  			if (nrn_allthread_handle) { (*nrn_allthread_handle)(); }
  			if (stoprun) { break; }
  			step_group_begin = step_group_end;
  		}
  	}
  	t = nrn_threads[0]._t;
  }
  
  void* nrn_fixed_step_group_thread(NrnThread* nth) {
  	int i;
  	nth->_stop_stepping = 0;
  	for (i = step_group_begin; i < step_group_n; ++i) {
          nrn::Instrumentor::phase p_timestep("timestep");
          nrn_fixed_step_thread(nth);
  		if (nth->_stop_stepping) {
  			if (nth->id == 0) { step_group_end = i + 1; }
  			nth->_stop_stepping = 0;
  			return nullptr;
  		}
  	}
  	if (nth->id == 0) { step_group_end = step_group_n; }
  	return nullptr;
  }
  
  void* nrn_fixed_step_thread(NrnThread* nth) {
  	double wt;
      {
          nrn::Instrumentor::phase p("deliver-events");
          deliver_net_events(nth);
      }
  
      wt = nrnmpi_wtime();
  	nrn_random_play();
  #if ELIMINATE_T_ROUNDOFF
  	nth->nrn_ndt_ += .5;
  	nth->_t = nrn_tbase_ + nth->nrn_ndt_ * nrn_dt_;
  #else
  	nth->_t += .5 * nth->_dt;
  #endif
  	fixed_play_continuous(nth);
  	setup_tree_matrix(nth);
      {
          nrn::Instrumentor::phase p("matrix-solver");
          nrn_solve(nth);
      }
      {
          nrn::Instrumentor::phase p("second-order-cur");
          second_order_cur(nth);
      }
      {
          nrn::Instrumentor::phase p("update");
          update(nth);
      }
      CTADD
  /*
    To simplify the logic,
    if there is no nrnthread_v_transfer then there cannot be an nrnmpi_v_transfer.
  */
  	if (!nrnthread_v_transfer_) {
  		nrn_fixed_step_lastpart(nth);
  	}
      return nullptr;
  }
  
  extern void nrn_extra_scatter_gather(int direction, int tid);
  extern void nrn_ba(NrnThread*, int);
  
  void* nrn_fixed_step_lastpart(NrnThread* nth) {
  	CTBEGIN
  #if ELIMINATE_T_ROUNDOFF
  	nth->nrn_ndt_ += .5;
  	nth->_t = nrn_tbase_ + nth->nrn_ndt_ * nrn_dt_;
  #else
  	nth->_t += .5 * nth->_dt;
  #endif
  	fixed_play_continuous(nth);
  	nrn_extra_scatter_gather(0, nth->id);
  	nonvint(nth);
  	nrn_ba(nth, AFTER_SOLVE);
  	fixed_record_continuous(nth);
      CTADD
      {
          nrn::Instrumentor::phase p("deliver-events");
          nrn_deliver_events(nth); /* up to but not past texit */
      }
      return nullptr;
  }
  
  /* nrn_fixed_step_thread is split into three pieces */
  
  void* nrn_ms_treeset_through_triang(NrnThread* nth) {
  	double wt;
  	deliver_net_events(nth);
  	wt = nrnmpi_wtime();
  	nrn_random_play();
  #if ELIMINATE_T_ROUNDOFF
  	nth->nrn_ndt_ += .5;
  	nth->_t = nrn_tbase_ + nth->nrn_ndt_ * nrn_dt_;
  #else
  	nth->_t += .5 * nth->_dt;
  #endif
  	fixed_play_continuous(nth);
  	setup_tree_matrix(nth);
  	nrn_multisplit_triang(nth);
  	CTADD
  	return nullptr;
  }
  void* nrn_ms_reduce_solve(NrnThread* nth) {
  	nrn_multisplit_reduce_solve(nth);
  	return nullptr;
  }
  void* nrn_ms_bksub(NrnThread* nth) {
  	CTBEGIN
  	nrn_multisplit_bksub(nth);
  	second_order_cur(nth);
  	update(nth);
  	CTADD
  /* see above comment in nrn_fixed_step_thread */
  	if (!nrnthread_v_transfer_) {
  		nrn_fixed_step_lastpart(nth);
  	}
  	return nullptr;
  }
  void* nrn_ms_bksub_through_triang(NrnThread* nth) {
  	nrn_ms_bksub(nth);
  	if (nth->_stop_stepping) {
  		nth->_stop_stepping = 0;
  		if (nth == nrn_threads) {
  			step_group_n = 1;
  		}
  		return nullptr;
  	}
  	nrn_ms_treeset_through_triang(nth);
  	return nullptr;
  }
  
  
  static void update(NrnThread* _nt)
  {
  	int i, i1, i2;
  	i1 = 0;
  	i2 = _nt->end;
  #if CACHEVEC
      if (use_cachevec) {
  	/* do not need to worry about linmod or extracellular*/
  	if (secondorder) {
  		for (i=i1; i < i2; ++i) {
  			VEC_V(i) += 2.*VEC_RHS(i);
  		}
  	}else{
  		for (i=i1; i < i2; ++i) {
  			VEC_V(i) += VEC_RHS(i);
  		}
  	}
      }else
  #endif
      {	/* use original non-vectorized update */
   	if (secondorder) {
  #if _CRAY
  #pragma _CRI ivdep
  #endif
  		for (i=i1; i < i2; ++i) {
  			NODEV(_nt->_v_node[i]) += 2.*NODERHS(_nt->_v_node[i]);
  		}
  	}else{
  #if _CRAY
  #pragma _CRI ivdep
  #endif
  		for (i=i1; i < i2; ++i) {
  			NODEV(_nt->_v_node[i]) += NODERHS(_nt->_v_node[i]);
  		}
  		if (use_sparse13) {
  			nrndae_update();
  		}
  	}
      } /* end of non-vectorized update */
  
  #if EXTRACELLULAR
  	nrn_update_2d(_nt);
  #endif
  	if (nrnthread_vi_compute_) { (*nrnthread_vi_compute_)(_nt); }
  #if I_MEMBRANE
  	if (_nt->tml) {
  		assert(_nt->tml->index == CAP);
  		nrn_capacity_current(_nt, _nt->tml->ml);
  	}
  #endif
  	if (nrn_use_fast_imem) { nrn_calc_fast_imem(_nt); }
  }
  
  void nrn_calc_fast_imem(NrnThread* _nt) {
  	int i;
  	int i1 = 0;
  	int i3 = _nt->end;
  	double* pd = _nt->_nrn_fast_imem->_nrn_sav_d;
  	double* prhs = _nt->_nrn_fast_imem->_nrn_sav_rhs;
      if (use_cachevec) {
  	for (i = i1; i < i3 ; ++i) {
  		prhs[i] = (pd[i]*VEC_RHS(i) + prhs[i])*VEC_AREA(i)*0.01;
  	}
      }else{
  	for (i = i1; i < i3 ; ++i) {
  		Node* nd = _nt->_v_node[i];
  		prhs[i] = (pd[i]*NODERHS(nd) + prhs[i])*NODEAREA(nd)*0.01;
  	}
      }
  }
  
  void nrn_calc_fast_imem_fixedstep_init(NrnThread* _nt) {
    // nrn_rhs() is called near end of nrn_finitialize() via setup_tree_matrix()
    // and nrn_rhs() sets up _nrn_sav_rhs as -total ionic_current (without
    // ELECTRODE_CURRENT contributions) and RHS as axial + ionic + stim currents
    // So sum, scaled by area, is i_membrane_ in nA.
    // (Note: capacitance does not appear on rhs because delta_v is the
    // variable in the current balance equations set up by setup_tree_matrix.)
    // Warning: Have not thought deeply about extracellular or LinearMechanism.
    //          But there is a good chance things are ok. But needs testing.
    // I don't believe this is used by Cvode or IDA.
  	int i;
  	int i1 = 0;
  	int i3 = _nt->end;
  	double* prhs = _nt->_nrn_fast_imem->_nrn_sav_rhs;
      if (use_cachevec) {
  	for (i = i1; i < i3 ; ++i) {
  		prhs[i] = (VEC_RHS(i) + prhs[i])*VEC_AREA(i)*0.01;
  	}
      }else{
  	for (i = i1; i < i3 ; ++i) {
  		Node* nd = _nt->_v_node[i];
  		prhs[i] = (NODERHS(nd) + prhs[i])*NODEAREA(nd)*0.01;
  	}
      }
  }
  
  void fcurrent(void)
  {
  	int i;
  	if (tree_changed) {
  		setup_topology();
  	}
  	if (v_structure_change) {
  		v_setup_vectors();
  	}
  	if (diam_changed) {
  		recalc_diam();
  	}
  
  	dt2thread(-1.);
  	nrn_thread_table_check();
  	state_discon_allowed_ = 0;
  	nrn_multithread_job(setup_tree_matrix);
  	state_discon_allowed_ = 1;
  	hoc_retpushx(1.);
  }
  
  void nrn_print_matrix(NrnThread* _nt) {
  	extern int section_count;
  	extern Section** secorder;
  	int isec, inode;
  	Section* sec;
  	Node* nd;
      if (use_sparse13) {
  	if(ifarg(1) && chkarg(1, 0., 1.) == 0.) {
  		spPrint(_nt->_sp13mat, 1, 0, 1);
  	}else{
  		int i, n = spGetSize(_nt->_sp13mat, 0);
  		spPrint(_nt->_sp13mat, 1, 1, 1);
  		for (i=1; i <= n; ++i) {
  			Printf("%d %g\n", i, _nt->_actual_rhs[i]);
  		}
  	}
      }else if (_nt) {
  	for (inode = 0; inode <  _nt->end; ++inode) {
  		nd = _nt->_v_node[inode];
  Printf("%d %g %g %g %g\n", inode, ClassicalNODEB(nd), ClassicalNODEA(nd), NODED(nd), NODERHS(nd));
  	}		
      }else{
  	for (isec = 0; isec < section_count; ++isec) {
  		sec = secorder[isec];
  		for (inode = 0; inode < sec->nnode; ++inode) {
  			nd = sec->pnode[inode];
  Printf("%d %d %g %g %g %g\n", isec, inode, ClassicalNODEB(nd), ClassicalNODEA(nd), NODED(nd), NODERHS(nd));
  		}
  	}
      }
  }
  
  void fmatrix(void) {
  	if (ifarg(1)) {
  		double x;
  		Section* sec;
  		int id;
  		Node* nd;
  		NrnThread* _nt;
  		nrn_seg_or_x_arg(1, &sec, &x);
  		id = (int)chkarg(2, 1., 4.);
  		nd = node_exact(sec, x);
  		_nt = nd->_nt;
  		switch (id) {
  		case 1: hoc_retpushx(NODEA(nd)); break;
  		case 2: hoc_retpushx(NODED(nd)); break;
  		case 3: hoc_retpushx(NODEB(nd)); break;
  		case 4: hoc_retpushx(NODERHS(nd)); break;
  		}
  		return;
  	}
  	nrn_print_matrix(nrn_threads);
  	hoc_retpushx(1.);
  	return;
  }
  
  void nonvint(NrnThread* _nt)
  {
  #if VECTORIZE
  	int i=0;
  	double w;
  	int measure = 0;
  	NrnThreadMembList* tml;
  #if 1 || PARANEURON
  	/* nrnmpi_v_transfer if needed was done earlier */
      if (nrnthread_v_transfer_) {
          nrn::Instrumentor::phase p_gap("gap-v-transfer");
          (*nrnthread_v_transfer_)(_nt);
      }
  #endif
      nrn::Instrumentor::phase_begin("state-update");
      if (_nt->id == 0 && nrn_mech_wtime_) { measure = 1; }
  	errno = 0;
  	for (tml = _nt->tml; tml; tml = tml->next) if (memb_func[tml->index].state) {
              std::string mechname("state-");
              mechname += memb_func[tml->index].sym->name;
              nrn::Instrumentor::phase_begin(mechname.c_str());
              Pvmi s = memb_func[tml->index].state;
              nrn::Instrumentor::phase_end(mechname.c_str());
              if (measure) {
                  w = nrnmpi_wtime(); }
  		(*s)(_nt, tml->ml, tml->index);
  		if (measure) { nrn_mech_wtime_[tml->index] += nrnmpi_wtime() - w; }
  		if (errno) {
  			if (nrn_errno_check(i)) {
  hoc_warning("errno set during calculation of states", (char*)0);
  			}
  		}
    	  }
  	long_difus_solve(0, _nt); /* if any longitudinal diffusion */
  	nrn_nonvint_block_fixed_step_solve(_nt->id);
      nrn::Instrumentor::phase_end("state-update");
  #endif
  }
  
  #if VECTORIZE
  int nrn_errno_check(int i)
  {
  	int ierr;
  	ierr = hoc_errno_check();
  	if (ierr) {
  		fprintf(stderr, 
  "%d errno=%d at t=%g during call to mechanism %s\n", nrnmpi_myid, ierr, t, memb_func[i].sym->name);
  	}
  	return ierr;
  }
  #else
  int nrn_errno_check(Prop* p, int inode, Section* sec)
  {
  	int ierr;
  	char* secname();
  	ierr = hoc_errno_check();
  	if (ierr) {
  		fprintf(stderr, 
  "%d errno set at t=%g during call to mechanism %s at node %d in section %s\n",
  nrnmpi_myid, t, memb_func[p->type].sym->name, inode, secname(sec));
  	}
  	return ierr;
  }
  #endif
  
  void frecord_init(void) { /* useful when changing states after an finitialize() */
  	int i;
  	dt2thread(-1);
  	nrn_record_init();
  	if (!cvode_active_) {
  		for (i=0; i < nrn_nthread; ++i) {
  			fixed_record_continuous(nrn_threads + i);
  		}
  	}
  	hoc_retpushx(1.);
  }
  
  void verify_structure(void) {
  	if (tree_changed) {
  		setup_topology();
  	}
  	if (v_structure_change) {
  		v_setup_vectors();
  	}
  	if (diam_changed) {
  		recalc_diam();
  	}
  	nrn_solver_prepare(); /* cvode ready to be used */
  }
  	
  void nrn_finitialize(int setv, double v) {
      int iord, i;
      NrnThread *_nt;
      extern int _ninits;
      extern short *nrn_is_artificial_;
      ++_ninits;
  
      nrn::Instrumentor::phase_begin("finitialize");
      nrn_fihexec(3); /* model structure changes can be made */
      verify_structure();
  #if ELIMINATE_T_ROUNDOFF
      nrn_ndt_ = 0.; nrn_dt_ = dt; nrn_tbase_ = 0.;
  #else
      t = 0.;
      dt2thread(-1.);
  #endif
      if (cvode_active_) {
          nrncvode_set_t(t);
      }
      nrn_thread_table_check();
      clear_event_queue();
      nrn_spike_exchange_init();
      nrn_random_play();
  #if VECTORIZE
      nrn_play_init(); /* Vector.play */
      for (i = 0; i < nrn_nthread; ++i) {
          nrn_deliver_events(nrn_threads + i); /* The play events at t=0 */
      }
      if (setv) {
  #if _CRAY
  #pragma _CRI ivdep
  #endif
          FOR_THREADS(_nt)
          for (i = 0; i < _nt->end; ++i) {
              NODEV(_nt->_v_node[i]) = v;
          }
      }
  #if 1 || PARANEURON
      if (nrnthread_vi_compute_) FOR_THREADS(_nt)
      {
          (*nrnthread_vi_compute_)(_nt);
      }
      {
          nrn::Instrumentor::phase p_gap("gap-v-transfer");
          if (nrnmpi_v_transfer_) {
              (nrnmpi_v_transfer_)();
          }
          if (nrnthread_v_transfer_)
              FOR_THREADS(_nt) {
                  (*nrnthread_v_transfer_)(_nt);
              }
      }
  #endif
      nrn_fihexec(0); /* after v is set but before INITIAL blocks are called*/
      for (i = 0; i < nrn_nthread; ++i) {
          nrn_ba(nrn_threads + i, BEFORE_INITIAL);
      }
  #if _CRAY
      cray_node_init();
  #endif
      /* the INITIAL blocks are ordered so that mechanisms that write
         concentrations are after ions and before mechanisms that read
         concentrations.
      */
      /* the memblist list in NrnThread is already so ordered */
  #if MULTICORE
      for (i = 0; i < nrn_nthread; ++i) {
          NrnThread *nt = nrn_threads + i;
          nrn_nonvint_block_init(nt->id);
          NrnThreadMembList *tml;
          for (tml = nt->tml; tml; tml = tml->next) {
              Pvmi s = memb_func[tml->index].initialize;
              if (s) {
                  (*s)(nt, tml->ml, tml->index);
              }
          }
      }
  #endif
      for (iord = 0; iord < n_memb_func; ++iord) {
          i = memb_order_[iord];
          /* first clause due to MULTICORE */
          if (nrn_is_artificial_[i])
              if (memb_func[i].initialize) {
                  Pvmi s = memb_func[i].initialize;
  #if 0
                  if (nrn_is_artificial_[i]) {
                  /*
                       I hope the space saving of the memb_list arrays is worth
                       doing this specifically. And if art cells are needed
                       in the memb_list anywhere else we will have do do something
                       similar to this. This gives up vectorization but this is only
                       initialization and all the other use of artcell should be
                       event driven.
                  */
                  Prop* p;
                  hoc_Item* q;
                  hoc_List* list = nrn_pnt_template_[i]->olist;
                  ITERATE(q, list) {
                      Object* obj = OBJ(q);
                      Prop* p = ((Point_process*)obj->u.this_pointer)->prop;
                      (*s)((Node*)0, p->param, p->dparam);
                  }
                  }else if (memb_list[i].nodecount){
  #else
                  if (memb_list[i].nodecount) {
  #endif
                      (*s)(nrn_threads, memb_list + i, i);
                  }
                  if (errno) {
                      if (nrn_errno_check(i)) {
                          hoc_warning("errno set during call to INITIAL block", (char *) 0);
                      }
                  }
              }
      }
  #endif
      if (use_sparse13) {
          nrndae_init();
      }
  
      init_net_events();
      for (i = 0; i < nrn_nthread; ++i) {
          nrn_ba(nrn_threads + i, AFTER_INITIAL);
      }
      nrn_fihexec(1); /* after INITIAL blocks, before fcurrent*/
  
      for (i = 0; i < nrn_nthread; ++i) {
          nrn_deliver_events(nrn_threads + i); /* The INITIAL sent events at t=0 */
      }
      if (cvode_active_) {
          cvode_finitialize(t);
          nrn_record_init();
      } else {
          state_discon_allowed_ = 0;
          for (i = 0; i < nrn_nthread; ++i) {
              setup_tree_matrix(nrn_threads + i);
              if (nrn_use_fast_imem) {
                  nrn_calc_fast_imem_fixedstep_init(nrn_threads + i);
              }
          }
          state_discon_allowed_ = 1;
          nrn_record_init();
          for (i = 0; i < nrn_nthread; ++i) {
              fixed_record_continuous(nrn_threads + i);
          }
      }
      for (i = 0; i < nrn_nthread; ++i) {
          nrn_deliver_events(nrn_threads + i); /* The record events at t=0 */
      }
  #if NRNMPI
      nrn_spike_exchange(nrn_threads);
  #endif
      if (nrn_allthread_handle) { (*nrn_allthread_handle)(); }
  
      nrn_fihexec(2); /* just before return */
      nrn::Instrumentor::phase_end("finitialize");
  }
  
  void finitialize(void) {
      int setv;
      double v = 0.0;
      setv = 0;
      if (ifarg(1)) {
          v = *getarg(1);
          setv = 1;
      }
      tstopunset;
      nrn_finitialize(setv, v);
      tstopunset;
      hoc_retpushx(1.);
  }
  
  
  static void batch_close() {
      if (batch_file) {
          fclose(batch_file);
          batch_file = 0;
      }
  }
  
  static void batch_out() {
      if (batch_file) {
          int i;
          for (i = 0; i < batch_n; ++i) {
              fprintf(batch_file, " %g", *batch_var[i]);
          }
          fprintf(batch_file, "\n");
      }
  }
  
  void batch_save(void) {
      int i;
      if (!ifarg(1)) {
          batch_n = 0;
      } else {
          for (i = 1; ifarg(i); ++i) {
              if (batch_size <= batch_n) {
                  batch_size += 20;
                  batch_var = (double **) erealloc(batch_var, batch_size * sizeof(double *));
              }
              batch_var[batch_n] = hoc_pgetarg(i);
              ++batch_n;
          }
      }
      hoc_retpushx(1.);
  }
  
  void nrn_ba(NrnThread *nt, int bat) {
      NrnThreadBAList *tbl;
      int i;
      for (tbl = nt->tbl[bat]; tbl; tbl = tbl->next) {
          nrn_bamech_t f = tbl->bam->f;
          int type = tbl->bam->type;
          Memb_list *ml = tbl->ml;
          for (i = 0; i < ml->nodecount; ++i) {
              (*f)(ml->nodelist[i], ml->data[i], ml->pdata[i], ml->_thread, nt);
          }
      }
  }
  
  typedef int (*NonVintBlockItem)(int method, int size, double *pd1, double *pd2, int tid);
  /* list to store the nrn_nonvint_block functions */
  static std::vector<NonVintBlockItem> nonvint_block_list;
  
  int nrn_nonvint_block_exe(int method, int size, double *pd1, double *pd2, int tid) {
      /* execute all functions in nonvint_block_list and return the sum of the
       * return values
       */
      int rval, sum = 0;
  
      for (auto& func: nonvint_block_list) {
          rval = (*func)(method, size, pd1, pd2, tid);
          if (rval == -1) {
              hoc_execerror("nrn_nonvint_block error", 0);
          } else {
              sum += rval;
          }
          if (method == NONVINT_ODE_COUNT) {
              size += rval;
          }
      }
  
      return sum;
  }
  
  extern "C" int set_nonvint_block(NonVintBlockItem func) {
  
      /* store new_nrn_nonvint_block functions in a list */
      nonvint_block_list.push_back(func);
  
      /* could this be set directly in nrn_nonvint_block_helper? */
      nrn_nonvint_block = &nrn_nonvint_block_exe;
  
      return 0;
  }
  
  extern "C" int unset_nonvint_block(NonVintBlockItem func) {
      // remove new_nrn_nonvint_block functions from the list
      // in c++-20 one could say std::erase(nonvint_block_list, func);
      auto n = nonvint_block_list.size();
      for (size_t i = 0; i < n; ++i) {
          if (func == nonvint_block_list[i]) {
              nonvint_block_list.erase(nonvint_block_list.begin() + i);
              break;
          }
      }
      if (nonvint_block_list.empty()) {
          nrn_nonvint_block = NULL;
      }
      return 0;
  }
  
  int nrn_nonvint_block_helper(int method, int size, double *pd1, double *pd2, int tid) {
      assert(nrn_nonvint_block);
      int rval = (*nrn_nonvint_block)(method, size, pd1, pd2, tid);
      if (rval == -1) {
          hoc_execerror("nrn_nonvint_block error", 0);
      }
      return rval;
  }
  
  /*
     Derived from scopmath/euler.cpp. Here because scopmath does not know
     about NrnThread
  */
  #include "nrniv_mf.h"
  
  #undef SUCCESS
  #define SUCCESS 0
  #define der_(arg)  p[der[arg]]
  #define var_(arg)  p[var[arg]]
  
  /* ARGSUSED */
  extern "C" int euler_thread(int neqn, int *var, int *der, double *p,
                   int (*func)(double *, Datum *, Datum *, NrnThread *),
                   Datum *ppvar, Datum *thread, NrnThread *nt) {
      int i;
      double dt = nt->_dt;
  
      /* Calculate the derivatives */
  
      (*func)(p, ppvar, thread, nt);
  
      /* Update dependent variables --- note defines in euler above*/
  
      for (i = 0; i < neqn; i++)
          var_(i) += dt * (der_(i));
  
      return (SUCCESS);
  }
  
  