libpspm
Public Member Functions | Static Public Member Functions | Public Attributes | Private Attributes | List of all members
LSODA Class Reference

#include <lsoda.h>

Collaboration diagram for LSODA:
[legend]

Public Member Functions

 LSODA ()
 
 ~LSODA ()
 
size_t idamax1 (const std::vector< double > &dx, const size_t n, const size_t offset)
 
void dscal1 (const double da, std::vector< double > &dx, const size_t n, const size_t offset)
 
double ddot1 (const std::vector< double > &a, const std::vector< double > &b, const size_t n, const size_t offsetA, const size_t offsetB)
 
void daxpy1 (const double da, const std::vector< double > &dx, std::vector< double > &dy, const size_t n, const size_t offsetX, const size_t offsetY)
 
void dgesl (const std::vector< std::vector< double >> &a, const size_t n, std::vector< int > &ipvt, std::vector< double > &b, const size_t job)
 
void dgefa (std::vector< std::vector< double >> &a, const size_t n, std::vector< int > &ipvt, size_t *const info)
 
template<class Functor >
void prja (const size_t neq, std::vector< double > &y, Functor &derivs, void *_data)
 
template<class Functor , class AfterStep >
void lsoda (Functor &derivs, AfterStep &after_step, const size_t neq, std::vector< double > &y, double *t, double tout, int itask, int *istate, int iopt, int jt, std::array< int, 7 > &iworks, std::array< double, 4 > &rworks, void *_data)
 
template<class Functor >
void correction (const size_t neq, std::vector< double > &y, Functor &derivs, size_t *corflag, double pnorm, double *del, double *delp, double *told, size_t *ncf, double *rh, size_t *m, void *_data)
 
template<class Functor >
void stoda (const size_t neq, std::vector< double > &y, Functor &derivs, void *_data)
 
template<class Functor , class AfterStep >
void lsoda_update (Functor &derivs, AfterStep &after_step, const size_t neq, std::vector< double > &y, double *t, const double tout, void *const _data, double rtol=1e-6, double atol=1e-6)
 
template<class AfterStep >
void successreturn (AfterStep &after_step, std::vector< double > &y, double *t, int itask, int ihit, double tcrit, int *istate)
 
void terminate (int *istate)
 
void terminate2 (std::vector< double > &y, double *t)
 
void _freevectors (void)
 
void ewset (const std::vector< double > &ycur)
 
void resetcoeff (void)
 
void solsy (std::vector< double > &y)
 
void endstoda (void)
 
void orderswitch (double *rhup, double dsm, double *pdh, double *rh, size_t *orderflag)
 
void intdy (double t, int k, std::vector< double > &dky, int *iflag)
 
void corfailure (double *told, double *rh, size_t *ncf, size_t *corflag)
 
void methodswitch (double dsm, double pnorm, double *pdh, double *rh)
 
void cfode (int meth_)
 
void scaleh (double *rh, double *pdh)
 
double fnorm (int n, const std::vector< std::vector< double >> &a, const std::vector< double > &w)
 
double vmnorm (const size_t n, const std::vector< double > &v, const std::vector< double > &w)
 
void resize_system (int nyh, int lenyh)
 
void set_istate (int n)
 
int get_istate () const
 
int get_fncalls ()
 

Static Public Member Functions

static bool abs_compare (double a, double b)
 

Public Attributes

void * param = nullptr
 

Private Attributes

double ETA = 2.2204460492503131e-16
 
size_t ml
 
size_t mu
 
size_t imxer
 
double sqrteta
 
std::array< size_t, 3 > mord
 
std::array< double, 13 > sm1
 
std::array< double, 14 > el
 
std::array< double, 13 > cm1
 
std::array< double, 6 > cm2
 
std::array< std::array< double, 14 >, 13 > elco
 
std::array< std::array< double, 4 >, 13 > tesco
 
size_t illin
 
size_t init
 
size_t ierpj
 
size_t iersl
 
size_t jcur
 
size_t l
 
size_t miter
 
size_t maxord
 
size_t maxcor
 
size_t msbp
 
size_t mxncf
 
int kflag
 
int jstart
 
size_t ixpr = 0
 
size_t jtyp
 
size_t mused
 
size_t mxordn
 
size_t mxords = 12
 
size_t meth_
 
size_t n
 
size_t nq
 
size_t nst
 
size_t nfe
 
size_t nje
 
size_t nqu
 
size_t mxstep
 
size_t mxhnil
 
size_t nslast
 
size_t nhnil
 
size_t ntrep
 
size_t nyh
 
double ccmax
 
double el0
 
double h_ = .0
 
double hmin
 
double hmxi
 
double hu
 
double rc
 
double tn_ = 0.0
 
double tsw
 
double pdnorm
 
double conit
 
double crate
 
double hold
 
double rmax
 
size_t ialth
 
size_t ipup
 
size_t lmax
 
size_t nslp
 
double pdest
 
double pdlast
 
double ratio
 
int icount
 
int irflag
 
std::vector< double > ewt
 
std::vector< double > savf
 
std::vector< double > acor
 
std::vector< std::vector< double > > yh_
 
std::vector< std::vector< double > > wm_
 
std::vector< int > ipvt
 
int itol_ = 2
 
std::vector< double > rtol_
 
std::vector< double > atol_
 
int iState = 1
 
std::vector< double > yout
 

Detailed Description

Type definition of LSODA ode system. See the file test_LSODA.cpp for an example.

time, double y, array of double. dydt, array of double data, void*

void

Definition at line 38 of file lsoda.h.

Constructor & Destructor Documentation

◆ LSODA()

LSODA::LSODA ( )

Definition at line 34 of file lsoda.cpp.

35 {
36  // Initialize arrays.
37  mord = {{12, 5}};
38  sm1 = {{0., 0.5, 0.575, 0.55, 0.45, 0.35, 0.25, 0.2, 0.15, 0.1, 0.075, 0.05, 0.025}};
39  el = {{0}};
40  cm1 = {{0}};
41  cm2 = {{0}};
42 }
LSODA::sm1
std::array< double, 13 > sm1
Definition: lsoda.h:121
LSODA::mord
std::array< size_t, 3 > mord
Definition: lsoda.h:120
LSODA::cm1
std::array< double, 13 > cm1
Definition: lsoda.h:124
LSODA::cm2
std::array< double, 6 > cm2
Definition: lsoda.h:125
LSODA::el
std::array< double, 14 > el
Definition: lsoda.h:123

◆ ~LSODA()

LSODA::~LSODA ( )

Definition at line 44 of file lsoda.cpp.

45 {
46 }

Member Function Documentation

◆ _freevectors()

void LSODA::_freevectors ( void  )

Definition at line 831 of file lsoda.cpp.

832 {
833  // Does nothing. USE c++ memory mechanism here.
834 }

◆ abs_compare()

bool LSODA::abs_compare ( double  a,
double  b 
)
static

Definition at line 48 of file lsoda.cpp.

49 {
50  return (std::abs(a) < std::abs(b));
51 }

◆ cfode()

void LSODA::cfode ( int  meth_)

Definition at line 331 of file lsoda.cpp.

332 {
333  int i, nq, nqm1, nqp1;
334  double agamq, fnq, fnqm1, pc[13], pint, ragq, rqfac, rq1fac, tsign, xpin;
335  /*
336  cfode is called by the integrator routine to set coefficients
337  needed there. The coefficients for the current method, as
338  given by the value of meth_, are set for all orders and saved.
339  The maximum order assumed here is 12 if meth_ = 1 and 5 if meth_ = 2.
340  ( A smaller value of the maximum order is also allowed. )
341  cfode is called once at the beginning of the problem, and
342  is not called again unless and until meth_ is changed.
343 
344  The elco array contains the basic method coefficients.
345  The coefficients el[i], 1 < i < nq+1, for the method of
346  order nq are stored in elco[nq][i]. They are given by a generating
347  polynomial, i.e.,
348 
349  l(x) = el[1] + el[2]*x + ... + el[nq+1]*x^nq.
350 
351  For the implicit Adams method, l(x) is given by
352 
353  dl/dx = (x+1)*(x+2)*...*(x+nq-1)/factorial(nq-1), l(-1) = 0.
354 
355  For the bdf methods, l(x) is given by
356 
357  l(x) = (x+1)*(x+2)*...*(x+nq)/k,
358 
359  where k = factorial(nq)*(1+1/2+...+1/nq).
360 
361  The tesco array contains test constants used for the
362  local error test and the selection of step size and/or order.
363  At order nq, tesco[nq][k] is used for the selection of step
364  size at order nq-1 if k = 1, at order nq if k = 2, and at order
365  nq+1 if k = 3.
366  */
367  if(meth_ == 1) {
368  elco[1][1] = 1.;
369  elco[1][2] = 1.;
370  tesco[1][1] = 0.;
371  tesco[1][2] = 2.;
372  tesco[2][1] = 1.;
373  tesco[12][3] = 0.;
374  pc[1] = 1.;
375  rqfac = 1.;
376  for(nq = 2; nq <= 12; nq++) {
377  /*
378  The pc array will contain the coefficients of the polynomial
379 
380  p(x) = (x+1)*(x+2)*...*(x+nq-1).
381 
382  Initially, p(x) = 1.
383  */
384  rq1fac = rqfac;
385  rqfac = rqfac / (double)nq;
386  nqm1 = nq - 1;
387  fnqm1 = (double)nqm1;
388  nqp1 = nq + 1;
389  /*
390  Form coefficients of p(x)*(x+nq-1).
391  */
392  pc[nq] = 0.;
393  for(i = nq; i >= 2; i--)
394  pc[i] = pc[i - 1] + fnqm1 * pc[i];
395  pc[1] = fnqm1 * pc[1];
396  /*
397  Compute integral, -1 to 0, of p(x) and x*p(x).
398  */
399  pint = pc[1];
400  xpin = pc[1] / 2.;
401  tsign = 1.;
402  for(i = 2; i <= nq; i++) {
403  tsign = -tsign;
404  pint += tsign * pc[i] / (double)i;
405  xpin += tsign * pc[i] / (double)(i + 1);
406  }
407  /*
408  Store coefficients in elco and tesco.
409  */
410  elco[nq][1] = pint * rq1fac;
411  elco[nq][2] = 1.;
412  for(i = 2; i <= nq; i++)
413  elco[nq][i + 1] = rq1fac * pc[i] / (double)i;
414  agamq = rqfac * xpin;
415  ragq = 1. / agamq;
416  tesco[nq][2] = ragq;
417  if(nq < 12)
418  tesco[nqp1][1] = ragq * rqfac / (double)nqp1;
419  tesco[nqm1][3] = ragq;
420  } /* end for */
421  return;
422  } /* end if ( meth_ == 1 ) */
423 
424  /* meth_ = 2. */
425  pc[1] = 1.;
426  rq1fac = 1.;
427 
428  /*
429  The pc array will contain the coefficients of the polynomial
430  p(x) = (x+1)*(x+2)*...*(x+nq).
431  Initially, p(x) = 1.
432  */
433  for(nq = 1; nq <= 5; nq++) {
434  fnq = (double)nq;
435  nqp1 = nq + 1;
436  /*
437  Form coefficients of p(x)*(x+nq).
438  */
439  pc[nqp1] = 0.;
440  for(i = nq + 1; i >= 2; i--)
441  pc[i] = pc[i - 1] + fnq * pc[i];
442  pc[1] *= fnq;
443  /*
444  Store coefficients in elco and tesco.
445  */
446  for(i = 1; i <= nqp1; i++)
447  elco[nq][i] = pc[i] / pc[2];
448  elco[nq][2] = 1.;
449  tesco[nq][1] = rq1fac;
450  tesco[nq][2] = ((double)nqp1) / elco[nq][1];
451  tesco[nq][3] = ((double)(nq + 2)) / elco[nq][1];
452  rq1fac /= fnq;
453  }
454  return;
455 
456 } /* end cfode */
LSODA::elco
std::array< std::array< double, 14 >, 13 > elco
Definition: lsoda.h:127
LSODA::meth_
size_t meth_
Definition: lsoda.h:135
LSODA::tesco
std::array< std::array< double, 4 >, 13 > tesco
Definition: lsoda.h:128
LSODA::nq
size_t nq
Definition: lsoda.h:137
Here is the caller graph for this function:

◆ corfailure()

void LSODA::corfailure ( double *  told,
double *  rh,
size_t *  ncf,
size_t *  corflag 
)

Definition at line 532 of file lsoda.cpp.

533 {
534  ncf++;
535  rmax = 2.;
536  tn_ = *told;
537  for(size_t j = nq; j >= 1; j--)
538  for(size_t i1 = j; i1 <= nq; i1++)
539  for(size_t i = 1; i <= n; i++)
540  yh_[i1][i] -= yh_[i1 + 1][i];
541 
542  if(fabs(h_) <= hmin * 1.00001 || *ncf == mxncf) {
543  *corflag = 2;
544  return;
545  }
546  *corflag = 1;
547  *rh = 0.25;
548  ipup = miter;
549 }
LSODA::ipup
size_t ipup
Definition: lsoda.h:146
LSODA::miter
size_t miter
Definition: lsoda.h:130
LSODA::yh_
std::vector< std::vector< double > > yh_
Definition: lsoda.h:154
LSODA::tn_
double tn_
Definition: lsoda.h:142
LSODA::hmin
double hmin
Definition: lsoda.h:142
LSODA::n
size_t n
Definition: lsoda.h:137
LSODA::h_
double h_
Definition: lsoda.h:141
LSODA::nq
size_t nq
Definition: lsoda.h:137
LSODA::mxncf
size_t mxncf
Definition: lsoda.h:130
LSODA::rmax
double rmax
Definition: lsoda.h:144
Here is the caller graph for this function:

◆ correction()

template<class Functor >
void LSODA::correction ( const size_t  neq,
std::vector< double > &  y,
Functor &  derivs,
size_t *  corflag,
double  pnorm,
double *  del,
double *  delp,
double *  told,
size_t *  ncf,
double *  rh,
size_t *  m,
void *  _data 
)

Definition at line 78 of file lsoda.tpp.

81 {
82  double rm = 0.0, rate = 0.0, dcon = 0.0;
83 
84  /*
85  Up to maxcor corrector iterations are taken. A convergence test is
86  made on the r.m.s. norm of each correction, weighted by the error
87  weight vector ewt. The sum of the corrections is accumulated in the
88  vector acor[i]. The yh_ array is not altered in the corrector loop.
89  */
90 
91  *m = 0;
92  *corflag = 0;
93  *del = 0.;
94 
95  for(size_t i = 1; i <= n; i++)
96  y[i] = yh_[1][i];
97 
98  derivs(tn_, ++y.begin(), ++savf.begin(), _data);
99 
100  nfe++;
101  /*
102  If indicated, the matrix P = I - h_ * el[1] * J is reevaluated and
103  preprocessed before starting the corrector iteration. ipup is set
104  to 0 as an indicator that this has been done.
105  */
106  while(1) {
107  if(*m == 0) {
108  if(ipup > 0) {
109  prja(neq, y, derivs, _data);
110  ipup = 0;
111  rc = 1.;
112  nslp = nst;
113  crate = 0.7;
114  if(ierpj != 0) {
115  corfailure(told, rh, ncf, corflag);
116  return;
117  }
118  }
119  for(size_t i = 1; i <= n; i++)
120  acor[i] = 0.;
121  } /* end if ( *m == 0 ) */
122  if(miter == 0) {
123  /*
124  In case of functional iteration, update y directly from
125  the result of the last function evaluation.
126  */
127  for(size_t i = 1; i <= n; i++) {
128  savf[i] = h_ * savf[i] - yh_[2][i];
129  y[i] = savf[i] - acor[i];
130  }
131  *del = vmnorm(n, y, ewt);
132  for(size_t i = 1; i <= n; i++) {
133  y[i] = yh_[1][i] + el[1] * savf[i];
134  acor[i] = savf[i];
135  }
136  }
137  /* end functional iteration */
138  /*
139  In the case of the chord method, compute the corrector error,
140  and solve the linear system with that as right-hand side and
141  P as coefficient matrix.
142  */
143  else {
144  for(size_t i = 1; i <= n; i++)
145  y[i] = h_ * savf[i] - (yh_[2][i] + acor[i]);
146 
147  solsy(y);
148  *del = vmnorm(n, y, ewt);
149 
150  for(size_t i = 1; i <= n; i++) {
151  acor[i] += y[i];
152  y[i] = yh_[1][i] + el[1] * acor[i];
153  }
154  } /* end chord method */
155  /*
156  Test for convergence. If *m > 0, an estimate of the convergence
157  rate constant is stored in crate, and this is used in the test.
158 
159  We first check for a change of iterates that is the size of
160  roundoff error. If this occurs, the iteration has converged, and a
161  new rate estimate is not formed.
162  In all other cases, force at least two iterations to estimate a
163  local Lipschitz constant estimate for Adams method.
164  On convergence, form pdest = local maximum Lipschitz constant
165  estimate. pdlast is the most recent nonzero estimate.
166  */
167  if(*del <= 100. * pnorm * ETA)
168  break;
169  if(*m != 0 || meth_ != 1) {
170  if(*m != 0) {
171  rm = 1024.0;
172  if(*del <= (1024. * *delp))
173  rm = *del / *delp;
174  rate = std::max(rate, rm);
175  crate = std::max(0.2 * crate, rm);
176  }
177  dcon = *del * std::min(1., 1.5 * crate) / (tesco[nq][2] * conit);
178  if(dcon <= 1.) {
179  pdest = std::max(pdest, rate / fabs(h_ * el[1]));
180  if(pdest != 0.)
181  pdlast = pdest;
182  break;
183  }
184  }
185  /*
186  The corrector iteration failed to converge.
187  If miter != 0 and the Jacobian is out of date, prja is called for
188  the next try. Otherwise the yh_ array is retracted to its values
189  before prediction, and h_ is reduced, if possible. If h_ cannot be
190  reduced or mxncf failures have occured, exit with corflag = 2.
191  */
192  (*m)++;
193  if(*m == maxcor || (*m >= 2 && *del > 2. * *delp)) {
194  if(miter == 0 || jcur == 1) {
195  corfailure(told, rh, ncf, corflag);
196  return;
197  }
198  ipup = miter;
199  /*
200  Restart corrector if Jacobian is recomputed.
201  */
202  *m = 0;
203  rate = 0.;
204  *del = 0.;
205  for(size_t i = 1; i <= n; i++)
206  y[i] = yh_[1][i];
207 
208  derivs(tn_, ++y.begin(), ++savf.begin(), _data);
209 
210  nfe++;
211  }
212  /*
213  Iterate corrector.
214  */
215  else {
216  *delp = *del;
217  derivs(tn_, ++y.begin(), ++savf.begin(), _data);
218  nfe++;
219  }
220  } /* end while */
221 } /* end correction */
LSODA::nst
size_t nst
Definition: lsoda.h:137
LSODA::ETA
double ETA
Definition: lsoda.h:112
LSODA::pdlast
double pdlast
Definition: lsoda.h:148
LSODA::savf
std::vector< double > savf
Definition: lsoda.h:152
LSODA::ipup
size_t ipup
Definition: lsoda.h:146
LSODA::prja
void prja(const size_t neq, std::vector< double > &y, Functor &derivs, void *_data)
Definition: lsoda.tpp:6
LSODA::crate
double crate
Definition: lsoda.h:144
LSODA::rc
double rc
Definition: lsoda.h:142
LSODA::miter
size_t miter
Definition: lsoda.h:130
LSODA::yh_
std::vector< std::vector< double > > yh_
Definition: lsoda.h:154
LSODA::jcur
size_t jcur
Definition: lsoda.h:130
LSODA::conit
double conit
Definition: lsoda.h:144
LSODA::tn_
double tn_
Definition: lsoda.h:142
LSODA::meth_
size_t meth_
Definition: lsoda.h:135
LSODA::ierpj
size_t ierpj
Definition: lsoda.h:130
LSODA::nfe
size_t nfe
Definition: lsoda.h:137
LSODA::solsy
void solsy(std::vector< double > &y)
Definition: lsoda.cpp:560
LSODA::ewt
std::vector< double > ewt
Definition: lsoda.h:151
LSODA::n
size_t n
Definition: lsoda.h:137
LSODA::vmnorm
double vmnorm(const size_t n, const std::vector< double > &v, const std::vector< double > &w)
Definition: lsoda.cpp:501
LSODA::pdest
double pdest
Definition: lsoda.h:148
LSODA::tesco
std::array< std::array< double, 4 >, 13 > tesco
Definition: lsoda.h:128
LSODA::nslp
size_t nslp
Definition: lsoda.h:147
LSODA::h_
double h_
Definition: lsoda.h:141
LSODA::nq
size_t nq
Definition: lsoda.h:137
LSODA::corfailure
void corfailure(double *told, double *rh, size_t *ncf, size_t *corflag)
Definition: lsoda.cpp:532
LSODA::maxcor
size_t maxcor
Definition: lsoda.h:130
LSODA::acor
std::vector< double > acor
Definition: lsoda.h:153
LSODA::el
std::array< double, 14 > el
Definition: lsoda.h:123
Here is the call graph for this function:
Here is the caller graph for this function:

◆ daxpy1()

void LSODA::daxpy1 ( const double  da,
const std::vector< double > &  dx,
std::vector< double > &  dy,
const size_t  n,
const size_t  offsetX = 0,
const size_t  offsetY = 0 
)

Definition at line 96 of file lsoda.cpp.

98 {
99 
100  for(size_t i = 1; i <= n; i++)
101  dy[i + offsetY] = da * dx[i + offsetX] + dy[i + offsetY];
102 }
LSODA::n
size_t n
Definition: lsoda.h:137

◆ ddot1()

double LSODA::ddot1 ( const std::vector< double > &  a,
const std::vector< double > &  b,
const size_t  n,
const size_t  offsetA = 0,
const size_t  offsetB = 0 
)

Definition at line 87 of file lsoda.cpp.

89 {
90  double sum = 0.0;
91  for(size_t i = 1; i <= n; i++)
92  sum += a[i + offsetA] * b[i + offsetB];
93  return sum;
94 }
LSODA::n
size_t n
Definition: lsoda.h:137

◆ dgefa()

void LSODA::dgefa ( std::vector< std::vector< double >> &  a,
const size_t  n,
std::vector< int > &  ipvt,
size_t *const  info 
)

Definition at line 161 of file lsoda.cpp.

163 {
164  size_t j = 0, k = 0, i = 0;
165  double t = 0.0;
166 
167  /* Gaussian elimination with partial pivoting. */
168 
169  *info = 0;
170  for(k = 1; k <= n - 1; k++) {
171  /*
172  Find j = pivot index. Note that a[k]+k-1 is the address of
173  the 0-th element of the row vector whose 1st element is a[k][k].
174  */
175  j = idamax1(a[k], n - k + 1, k - 1) + k - 1;
176  ipvt[k] = j;
177  /*
178  Zero pivot implies this row already triangularized.
179  */
180  if(a[k][j] == 0.) {
181  *info = k;
182  continue;
183  }
184  /*
185  Interchange if necessary.
186  */
187  if(j != k) {
188  t = a[k][j];
189  a[k][j] = a[k][k];
190  a[k][k] = t;
191  }
192  /*
193  Compute multipliers.
194  */
195  t = -1. / a[k][k];
196  dscal1(t, a[k], n - k, k);
197 
198  /*
199  Column elimination with row indexing.
200  */
201  for(i = k + 1; i <= n; i++) {
202  t = a[i][j];
203  if(j != k) {
204  a[i][j] = a[i][k];
205  a[i][k] = t;
206  }
207  daxpy1(t, a[k], a[i], n - k, k, k);
208  }
209  } /* end k-loop */
210 
211  ipvt[n] = n;
212  if(a[n][n] == 0.)
213  *info = n;
214 }
LSODA::n
size_t n
Definition: lsoda.h:137
LSODA::dscal1
void dscal1(const double da, std::vector< double > &dx, const size_t n, const size_t offset)
Definition: lsoda.cpp:76
LSODA::idamax1
size_t idamax1(const std::vector< double > &dx, const size_t n, const size_t offset)
Definition: lsoda.cpp:54
LSODA::daxpy1
void daxpy1(const double da, const std::vector< double > &dx, std::vector< double > &dy, const size_t n, const size_t offsetX, const size_t offsetY)
Definition: lsoda.cpp:96
Here is the caller graph for this function:

◆ dgesl()

void LSODA::dgesl ( const std::vector< std::vector< double >> &  a,
const size_t  n,
std::vector< int > &  ipvt,
std::vector< double > &  b,
const size_t  job 
)

Definition at line 105 of file lsoda.cpp.

107 {
108  size_t k, j;
109  double t;
110 
111  /*
112  Job = 0, solve a * x = b.
113  */
114  if(job == 0) {
115  /*
116  First solve L * y = b.
117  */
118  for(k = 1; k <= n; k++) {
119  t = ddot1(a[k], b, k - 1);
120  b[k] = (b[k] - t) / a[k][k];
121  }
122  /*
123  Now solve U * x = y.
124  */
125  for(k = n - 1; k >= 1; k--) {
126  b[k] = b[k] + ddot1(a[k], b, n - k, k, k);
127  j = ipvt[k];
128  if(j != k) {
129  t = b[j];
130  b[j] = b[k];
131  b[k] = t;
132  }
133  }
134  return;
135  }
136  /*
137  Job = nonzero, solve Transpose(a) * x = b.
138 
139  First solve Transpose(U) * y = b.
140  */
141  for(k = 1; k <= n - 1; k++) {
142  j = ipvt[k];
143  t = b[j];
144  if(j != k) {
145  b[j] = b[k];
146  b[k] = t;
147  }
148  daxpy1(t, a[k], b, n - k, k, k);
149  }
150  /*
151  Now solve Transpose(L) * x = y.
152  */
153  for(k = n; k >= 1; k--) {
154  b[k] = b[k] / a[k][k];
155  t = -b[k];
156  daxpy1(t, a[k], b, k - 1);
157  }
158 }
LSODA::n
size_t n
Definition: lsoda.h:137
LSODA::ddot1
double ddot1(const std::vector< double > &a, const std::vector< double > &b, const size_t n, const size_t offsetA, const size_t offsetB)
Definition: lsoda.cpp:87
LSODA::daxpy1
void daxpy1(const double da, const std::vector< double > &dx, std::vector< double > &dy, const size_t n, const size_t offsetX, const size_t offsetY)
Definition: lsoda.cpp:96

◆ dscal1()

void LSODA::dscal1 ( const double  da,
std::vector< double > &  dx,
const size_t  n,
const size_t  offset = 0 
)

Definition at line 76 of file lsoda.cpp.

78 {
79  // FIXME: n is not used here. why?
80  (void)n;
81 
82  std::transform(dx.begin() + 1 + offset, dx.end(), dx.begin() + 1 + offset,
83  [&da](double x) -> double { return da * x; });
84 }
LSODA::n
size_t n
Definition: lsoda.h:137

◆ endstoda()

void LSODA::endstoda ( void  )

Definition at line 694 of file lsoda.cpp.

695 {
696  double r = 1. / tesco[nqu][2];
697  for(size_t i = 1; i <= n; i++)
698  acor[i] *= r;
699  hold = h_;
700  jstart = 1;
701 }
LSODA::nqu
size_t nqu
Definition: lsoda.h:137
LSODA::n
size_t n
Definition: lsoda.h:137
LSODA::tesco
std::array< std::array< double, 4 >, 13 > tesco
Definition: lsoda.h:128
LSODA::jstart
int jstart
Definition: lsoda.h:132
LSODA::h_
double h_
Definition: lsoda.h:141
LSODA::acor
std::vector< double > acor
Definition: lsoda.h:153
LSODA::hold
double hold
Definition: lsoda.h:144
Here is the caller graph for this function:

◆ ewset()

void LSODA::ewset ( const std::vector< double > &  ycur)

Definition at line 241 of file lsoda.cpp.

242 {
243  switch(itol_) {
244  case 1:
245  for(size_t i = 1; i <= n; i++)
246  ewt[i] = rtol_[1] * fabs(ycur[i]) + atol_[1];
247  break;
248  case 2:
249  for(size_t i = 1; i <= n; i++)
250  ewt[i] = rtol_[1] * fabs(ycur[i]) + atol_[i];
251  break;
252  case 3:
253  for(size_t i = 1; i <= n; i++)
254  ewt[i] = rtol_[i] * fabs(ycur[i]) + atol_[1];
255  break;
256  case 4:
257  for(size_t i = 1; i <= n; i++)
258  ewt[i] = rtol_[i] * fabs(ycur[i]) + atol_[i];
259  break;
260  }
261 
262 } /* end ewset */
LSODA::ewt
std::vector< double > ewt
Definition: lsoda.h:151
LSODA::n
size_t n
Definition: lsoda.h:137
LSODA::itol_
int itol_
Definition: lsoda.h:160
LSODA::atol_
std::vector< double > atol_
Definition: lsoda.h:162
LSODA::rtol_
std::vector< double > rtol_
Definition: lsoda.h:161
Here is the caller graph for this function:

◆ fnorm()

double LSODA::fnorm ( int  n,
const std::vector< std::vector< double >> &  a,
const std::vector< double > &  w 
)

Definition at line 509 of file lsoda.cpp.

519 {
520  double an = 0, sum = 0;
521 
522  for(size_t i = 1; i <= (size_t)n; i++) {
523  sum = 0.;
524  for(size_t j = 1; j <= (size_t)n; j++)
525  sum += fabs(a[i][j]) / w[j];
526  an = max(an, sum * w[i]);
527  }
528  return an;
529 }
LSODA::n
size_t n
Definition: lsoda.h:137
Here is the caller graph for this function:

◆ get_fncalls()

int LSODA::get_fncalls ( )

Definition at line 848 of file lsoda.cpp.

848  {
849  return nfe;
850 }
LSODA::nfe
size_t nfe
Definition: lsoda.h:137
Here is the caller graph for this function:

◆ get_istate()

int LSODA::get_istate ( ) const

Definition at line 853 of file lsoda.cpp.

853  {
854  return iState;
855 }
LSODA::iState
int iState
Definition: lsoda.h:165
Here is the caller graph for this function:

◆ idamax1()

size_t LSODA::idamax1 ( const std::vector< double > &  dx,
const size_t  n,
const size_t  offset = 0 
)

Definition at line 54 of file lsoda.cpp.

55 {
56 
57  size_t v = 0, vmax = 0;
58  size_t idmax = 1;
59  for(size_t i = 1; i <= n; i++) {
60  v = abs(dx[i + offset]);
61  if(v > vmax) {
62  vmax = v;
63  idmax = i;
64  }
65  }
66  return idmax;
67 
68  // Following has failed with seg-fault. Probably issue with STL.
69  // return std::max_element( dx.begin()+1+offset, dx.begin()+1+n, LSODA::abs_compare) -
70  // dx.begin() - offset;
71 }
LSODA::n
size_t n
Definition: lsoda.h:137

◆ intdy()

void LSODA::intdy ( double  t,
int  k,
std::vector< double > &  dky,
int *  iflag 
)

Definition at line 286 of file lsoda.cpp.

287 {
288  int ic, jp1 = 0;
289  double c, r, s, tp;
290 
291  *iflag = 0;
292  if(k < 0 || k > (int)nq) {
293  fprintf(stderr, "[intdy] k = %d illegal\n", k);
294  *iflag = -1;
295  return;
296  }
297  tp = tn_ - hu - 100. * ETA * (tn_ + hu);
298  if((t - tp) * (t - tn_) > 0.) {
299  fprintf(
300  stderr, "intdy -- t = %g illegal. t not in interval tcur - hu to tcur\n", t);
301  *iflag = -2;
302  return;
303  }
304  s = (t - tn_) / h_;
305  ic = 1;
306  for(size_t jj = l - k; jj <= nq; jj++)
307  ic *= jj;
308  c = (double)ic;
309  for(size_t i = 1; i <= n; i++)
310  dky[i] = c * yh_[l][i];
311 
312  for(int j = nq - 1; j >= k; j--) {
313  jp1 = j + 1;
314  ic = 1;
315  for(int jj = jp1 - k; jj <= j; jj++)
316  ic *= jj;
317  c = (double)ic;
318 
319  for(size_t i = 1; i <= n; i++)
320  dky[i] = c * yh_[jp1][i] + s * dky[i];
321  }
322  if(k == 0)
323  return;
324  r = pow(h_, (double)(-k));
325 
326  for(size_t i = 1; i <= n; i++)
327  dky[i] *= r;
328 
329 } /* end intdy */
LSODA::hu
double hu
Definition: lsoda.h:142
LSODA::ETA
double ETA
Definition: lsoda.h:112
LSODA::yh_
std::vector< std::vector< double > > yh_
Definition: lsoda.h:154
LSODA::tn_
double tn_
Definition: lsoda.h:142
LSODA::n
size_t n
Definition: lsoda.h:137
LSODA::h_
double h_
Definition: lsoda.h:141
LSODA::nq
size_t nq
Definition: lsoda.h:137
LSODA::l
size_t l
Definition: lsoda.h:130
Here is the caller graph for this function:

◆ lsoda()

template<class Functor , class AfterStep >
void LSODA::lsoda ( Functor &  derivs,
AfterStep &  after_step,
const size_t  neq,
std::vector< double > &  y,
double *  t,
double  tout,
int  itask,
int *  istate,
int  iopt,
int  jt,
std::array< int, 7 > &  iworks,
std::array< double, 4 > &  rworks,
void *  _data 
)

Definition at line 606 of file lsoda.tpp.

609 {
610  assert(tout > *t);
611 
612  int mxstp0 = 500, mxhnl0 = 10;
613 
614  int iflag = 0, lenyh = 0, ihit = 0;
615 
616  double atoli = 0, ayi = 0, big = 0, h0 = 0, hmax = 0, hmx = 0, rh = 0, rtoli = 0,
617  tcrit = 0, tdist = 0, tnext = 0, tol = 0, tolsf = 0, tp = 0, size = 0, sum = 0,
618  w0 = 0;
619 
620  /*
621  Block a.
622  This code block is executed on every call.
623  It tests *istate and itask for legality and branches appropriately.
624  If *istate > 1 but the flag init shows that initialization has not
625  yet been done, an error return occurs.
626  If *istate = 1 and tout = t, return immediately.
627  */
628 
629  if(*istate < 1 || *istate > 3) {
630  // fprintf(stderr, "[lsoda] illegal istate = %d\n", *istate);
631  std::cerr << "[lsoda] illegal istate = " << *istate << std::endl;
632  terminate(istate);
633  return;
634  }
635  if(itask < 1 || itask > 5) {
636  fprintf(stderr, "[lsoda] illegal itask = %d\n", itask);
637  terminate(istate);
638  return;
639  }
640  if(init == 0 && (*istate == 2 || *istate == 3)) {
641  fprintf(stderr, "[lsoda] istate > 1 but lsoda not initialized\n");
642  terminate(istate);
643  return;
644  }
645 
646  /*
647  Block b.
648  The next code block is executed for the initial call ( *istate = 1 ),
649  or for a continuation call with parameter changes ( *istate = 3 ).
650  It contains checking of all inputs and various initializations.
651 
652  First check legality of the non-optional inputs neq, itol, iopt,
653  jt, ml, and mu.
654  */
655 
656  if(*istate == 1 || *istate == 3) {
657  ntrep = 0;
658  if(neq <= 0) {
659  std::cerr << "[lsoda] neq = " << neq << " is less than 1." << std::endl;
660  terminate(istate);
661  return;
662  }
663  if(*istate == 3 && neq > n) {
664  std::cerr << "[lsoda] istate = 3 and neq increased" << std::endl;
665  terminate(istate);
666  return;
667  }
668  n = neq;
669  if(itol_ < 1 || itol_ > 4) {
670  std::cerr << "[lsoda] itol = " << itol_ << " illegal" << std::endl;
671  terminate(istate);
672  return;
673  }
674  if(iopt < 0 || iopt > 1) {
675  std::cerr << "[lsoda] iopt = " << iopt << " illegal" << std::endl;
676  terminate(istate);
677  return;
678  }
679  if(jt == 3 || jt < 1 || jt > 5) {
680  std::cerr << "[lsoda] jt = " << jt << " illegal" << std::endl;
681  terminate(istate);
682  return;
683  }
684  jtyp = jt;
685  if(jt > 2) {
686  ml = iworks[0];
687  mu = iworks[1];
688  if(ml >= n) {
689  std::cerr << "[lsoda] ml = " << ml << " not between 1 and neq" << std::endl;
690  terminate(istate);
691  return;
692  }
693  if(mu >= n) {
694  std::cerr << "[lsoda] mu = " << mu << " not between 1 and neq" << std::endl;
695  terminate(istate);
696  return;
697  }
698  }
699 
700  /* Next process and check the optional inpus. */
701  /* Default options. */
702  if(iopt == 0) {
703  ixpr = 0;
704  mxstep = mxstp0;
705  mxhnil = mxhnl0;
706  hmxi = 0.;
707  hmin = 0.;
708  if(*istate == 1) {
709  h0 = 0.;
710  mxordn = mord[0];
711  mxords = mord[1];
712  }
713  }
714  /* end if ( iopt == 0 ) */
715  /* Optional inputs. */
716  else /* if ( iopt = 1 ) */
717  {
718  ixpr = iworks[2];
719  if(ixpr > 1) {
720  std::cerr << "[lsoda] ixpr = " << ixpr << " is illegal" << std::endl;
721  terminate(istate);
722  return;
723  }
724 
725  mxstep = iworks[3];
726  if(mxstep == 0)
727  mxstep = mxstp0;
728  mxhnil = iworks[4];
729 
730  if(*istate == 1) {
731  h0 = rworks[1];
732  mxordn = iworks[5];
733 
734  if(mxordn == 0)
735  mxordn = 100;
736 
737  mxordn = std::min(mxordn, mord[0]);
738  mxords = iworks[6];
739 
740  // if mxords is not given use 100.
741  if(mxords == 0)
742  mxords = 100;
743 
744  mxords = std::min(mxords, mord[1]);
745 
746  if((tout - *t) * h0 < 0.) {
747  std::cerr << "[lsoda] tout = " << tout << " behind t = " << *t
748  << ". integration direction is given by " << h0 << std::endl;
749  terminate(istate);
750  return;
751  }
752  } /* end if ( *istate == 1 ) */
753  hmax = rworks[2];
754  if(hmax < 0.) {
755  std::cerr << "[lsoda] hmax < 0." << std::endl;
756  terminate(istate);
757  return;
758  }
759  hmxi = 0.;
760  if(hmax > 0)
761  hmxi = 1. / hmax;
762 
763  hmin = rworks[3];
764  if(hmin < 0.) {
765  std::cerr << "[lsoda] hmin < 0." << std::endl;
766  terminate(istate);
767  return;
768  }
769  } /* end else */ /* end iopt = 1 */
770  } /* end if ( *istate == 1 || *istate == 3 ) */
771  /*
772  If *istate = 1, meth_ is initialized to 1.
773 
774  Also allocate memory for yh_, wm_, ewt, savf, acor, ipvt.
775  */
776  if(*istate == 1) {
777  /*
778  If memory were not freed, *istate = 3 need not reallocate memory.
779  Hence this section is not executed by *istate = 3.
780  */
781  sqrteta = sqrt(ETA);
782  meth_ = 1;
783 
784  nyh = n;
785  lenyh = 1 + std::max(mxordn, mxords);
786  resize_system(nyh, lenyh);
787 
789  //yh_.clear(); yh_.resize(lenyh + 1, std::vector<double>(nyh + 1, 0.0));
790  //wm_.clear(); wm_.resize(nyh + 1, std::vector<double>(nyh + 1, 0.0));
791  //ewt.resize(1 + nyh, 0);
792  //savf.resize(1 + nyh, 0);
793  //acor.resize(nyh + 1, 0.0);
794  //ipvt.resize(nyh + 1, 0.0);
795  }
796  /*
797  Check rtol and atol for legality.
798  */
799  if(*istate == 1 || *istate == 3) {
800  rtoli = rtol_[1];
801  atoli = atol_[1];
802  for(size_t i = 1; i <= n; i++) {
803  if(itol_ >= 3)
804  rtoli = rtol_[i];
805  if(itol_ == 2 || itol_ == 4)
806  atoli = atol_[i];
807  if(rtoli < 0.) {
808  fprintf(stderr, "[lsoda] rtol = %g is less than 0.\n", rtoli);
809  terminate(istate);
810  return;
811  }
812  if(atoli < 0.) {
813  fprintf(stderr, "[lsoda] atol = %g is less than 0.\n", atoli);
814  terminate(istate);
815  return;
816  }
817  } /* end for */
818  } /* end if ( *istate == 1 || *istate == 3 ) */
819 
820  /* If *istate = 3, set flag to signal parameter changes to stoda. */
821  if(*istate == 3) {
822  jstart = -1;
823  }
824  /*
825  Block c.
826  The next block is for the initial call only ( *istate = 1 ).
827  It contains all remaining initializations, the initial call to f,
828  and the calculation of the initial step size.
829  The error weights in ewt are inverted after being loaded.
830  */
831  if(*istate == 1) {
832  tn_ = *t;
833  tsw = *t;
834  maxord = mxordn;
835  if(itask == 4 || itask == 5) {
836  tcrit = rworks[0];
837  if((tcrit - tout) * (tout - *t) < 0.) {
838  fprintf(stderr, "[lsoda] itask = 4 or 5 and tcrit behind tout\n");
839  terminate(istate);
840  return;
841  }
842  if(h0 != 0. && (*t + h0 - tcrit) * h0 > 0.)
843  h0 = tcrit - *t;
844  }
845 
846  jstart = 0;
847  nhnil = 0;
848  nst = 0;
849  nje = 0;
850  nslast = 0;
851  hu = 0.;
852  nqu = 0;
853  mused = 0;
854  miter = 0;
855  ccmax = 0.3;
856  maxcor = 3;
857  msbp = 20;
858  mxncf = 10;
859 
860  /* Initial call to f. */
861  assert((int)yh_.size() == lenyh + 1);
862  assert(yh_[0].size() == nyh + 1);
863 
864  derivs(*t, ++y.begin(), ++yh_[2].begin(), _data);
865  nfe = 1;
866 
867  /* Load the initial value vector in yh_. */
868  for(size_t i = 1; i <= n; i++)
869  yh_[1][i] = y[i];
870 
871  /* Load and invert the ewt array. ( h_ is temporarily set to 1. ) */
872  nq = 1;
873  h_ = 1.;
874  ewset(y);
875  for(size_t i = 1; i <= n; i++) {
876  if(ewt[i] <= 0.) {
877  std::cerr << "[lsoda] ewt[" << i << "] = " << ewt[i] << " <= 0.\n" << std::endl;
878  terminate2(y, t);
879  return;
880  }
881  ewt[i] = 1. / ewt[i];
882  }
883 
884  /*
885  The coding below computes the step size, h0, to be attempted on the
886  first step, unless the user has supplied a value for this.
887  First check that tout - *t differs significantly from zero.
888  A scalar tolerance quantity tol is computed, as max(rtol[i])
889  if this is positive, or max(atol[i]/fabs(y[i])) otherwise, adjusted
890  so as to be between 100*ETA and 0.001.
891  Then the computed value h0 is given by
892 
893  h0^(-2) = 1. / ( tol * w0^2 ) + tol * ( norm(f) )^2
894 
895  where w0 = max( fabs(*t), fabs(tout) ),
896  f = the initial value of the vector f(t,y), and
897  norm() = the weighted vector norm used throughout, given by
898  the vmnorm function routine, and weighted by the
899  tolerances initially loaded into the ewt array.
900 
901  The sign of h0 is inferred from the initial values of tout and *t.
902  fabs(h0) is made < fabs(tout-*t) in any case.
903  */
904  if(h0 == 0.) {
905  tdist = fabs(tout - *t);
906  w0 = std::max(fabs(*t), fabs(tout));
907  if(tdist < 2. * ETA * w0) {
908  fprintf(stderr, "[lsoda] tout too close to t to start integration\n ");
909  terminate(istate);
910  return;
911  }
912  tol = rtol_[1];
913  if(itol_ > 2) {
914  for(size_t i = 2; i <= n; i++)
915  tol = std::max(tol, rtol_[i]);
916  }
917  if(tol <= 0.) {
918  atoli = atol_[1];
919  for(size_t i = 1; i <= n; i++) {
920  if(itol_ == 2 || itol_ == 4)
921  atoli = atol_[i];
922  ayi = fabs(y[i]);
923  if(ayi != 0.)
924  tol = std::max(tol, atoli / ayi);
925  }
926  }
927  tol = std::max(tol, 100. * ETA);
928  tol = std::min(tol, 0.001);
929  sum = vmnorm(n, yh_[2], ewt);
930  sum = 1. / (tol * w0 * w0) + tol * sum * sum;
931  h0 = 1. / sqrt(sum);
932  h0 = std::min(h0, tdist);
933  h0 = h0 * ((tout - *t >= 0.) ? 1. : -1.);
934  } /* end if ( h0 == 0. ) */
935  /*
936  Adjust h0 if necessary to meet hmax bound.
937  */
938  rh = fabs(h0) * hmxi;
939  if(rh > 1.)
940  h0 /= rh;
941 
942  /*
943  Load h_ with h0 and scale yh_[2] by h0.
944  */
945  h_ = h0;
946  for(size_t i = 1; i <= n; i++)
947  yh_[2][i] *= h0;
948  } /* if ( *istate == 1 ) */
949  /*
950  Block d.
951  The next code block is for continuation calls only ( *istate = 2 or 3 )
952  and is to check stop conditions before taking a step.
953  */
954  if(*istate == 2 || *istate == 3) {
955  nslast = nst;
956  switch(itask) {
957  case 1:
958  if((tn_ - tout) * h_ >= 0.) {
959  intdy(tout, 0, y, &iflag);
960  if(iflag != 0) {
961  fprintf(stderr,
962  "[lsoda] trouble from intdy, itask = %d, tout = %g\n", itask,
963  tout);
964  terminate(istate);
965  return;
966  }
967  *t = tout;
968  *istate = 2;
969  illin = 0;
970  return;
971  }
972  break;
973  case 2:
974  break;
975  case 3:
976  tp = tn_ - hu * (1. + 100. * ETA);
977  if((tp - tout) * h_ > 0.) {
978  fprintf(
979  stderr, "[lsoda] itask = %d and tout behind tcur - hu\n", itask);
980  terminate(istate);
981  return;
982  }
983  if((tn_ - tout) * h_ < 0.)
984  break;
985  successreturn(after_step, y, t, itask, ihit, tcrit, istate);
986  return;
987  case 4:
988  tcrit = rworks[0];
989  if((tn_ - tcrit) * h_ > 0.) {
990  fprintf(stderr, "[lsoda] itask = 4 or 5 and tcrit behind tcur\n");
991  terminate(istate);
992  return;
993  }
994  if((tcrit - tout) * h_ < 0.) {
995  fprintf(stderr, "[lsoda] itask = 4 or 5 and tcrit behind tout\n");
996  terminate(istate);
997  return;
998  }
999  if((tn_ - tout) * h_ >= 0.) {
1000  intdy(tout, 0, y, &iflag);
1001  if(iflag != 0) {
1002  fprintf(stderr,
1003  "[lsoda] trouble from intdy, itask = %d, tout = %g\n", itask,
1004  tout);
1005  terminate(istate);
1006  return;
1007  }
1008  *t = tout;
1009  *istate = 2;
1010  illin = 0;
1011  after_step(*t, ++y.begin());
1012  return;
1013  }
1014  break;
1015  case 5:
1016  if(itask == 5) {
1017  tcrit = rworks[0];
1018  if((tn_ - tcrit) * h_ > 0.) {
1019  fprintf(stderr, "[lsoda] itask = 4 or 5 and tcrit behind tcur\n");
1020  terminate(istate);
1021  return;
1022  }
1023  }
1024  hmx = fabs(tn_) + fabs(h_);
1025  ihit = fabs(tn_ - tcrit) <= (100. * ETA * hmx);
1026  if(ihit) {
1027  *t = tcrit;
1028  successreturn(after_step, y, t, itask, ihit, tcrit, istate);
1029  return;
1030  }
1031  tnext = tn_ + h_ * (1. + 4. * ETA);
1032  if((tnext - tcrit) * h_ <= 0.)
1033  break;
1034  h_ = (tcrit - tn_) * (1. - 4. * ETA);
1035  if(*istate == 2)
1036  jstart = -2;
1037  break;
1038  } /* end switch */
1039  } /* end if ( *istate == 2 || *istate == 3 ) */
1040  /*
1041  Block e.
1042  The next block is normally executed for all calls and contains
1043  the call to the one-step core integrator stoda.
1044 
1045  This is a looping point for the integration steps.
1046 
1047  First check for too many steps being taken, update ewt ( if not at
1048  start of problem). Check for too much accuracy being requested, and
1049  check for h_ below the roundoff level in *t.
1050  */
1051  while(1) {
1052  if(*istate != 1 || nst != 0) {
1053  if((nst - nslast) >= mxstep) {
1054  std::cerr << "[lsoda] " << mxstep << " steps taken before reaching tout"
1055  << std::endl;
1056  *istate = -1;
1057  terminate2(y, t);
1058  return;
1059  }
1060 
1061  ewset(yh_[1]);
1062  for(size_t i = 1; i <= n; i++) {
1063  if(ewt[i] <= 0.) {
1064  std::cerr << "[lsoda] ewt[" << i << "] = " << ewt[i] << " <= 0." << std::endl;
1065  *istate = -6;
1066  terminate2(y, t);
1067  return;
1068  }
1069  ewt[i] = 1. / ewt[i];
1070  }
1071  }
1072  tolsf = ETA * vmnorm(n, yh_[1], ewt);
1073  if(tolsf > 0.01) {
1074  tolsf = tolsf * 200.;
1075  if(nst == 0) {
1076  fprintf(stderr, "lsoda -- at start of problem, too much accuracy\n");
1077  fprintf(stderr, " requested for precision of machine,\n");
1078  fprintf(stderr, " suggested scaling factor = %g\n", tolsf);
1079  terminate(istate);
1080  return;
1081  }
1082  fprintf(stderr, "lsoda -- at t = %g, too much accuracy requested\n", *t);
1083  fprintf(stderr, " for precision of machine, suggested\n");
1084  fprintf(stderr, " scaling factor = %g\n", tolsf);
1085  *istate = -2;
1086  terminate2(y, t);
1087  return;
1088  }
1089 
1090  if((tn_ + h_) == tn_) {
1091  nhnil++;
1092  if(nhnil <= mxhnil) {
1093  fprintf(stderr, "lsoda -- warning..internal t = %g and h_ = %g are\n",
1094  tn_, h_);
1095  fprintf(stderr,
1096  " such that in the machine, t + h_ = t on the next step\n");
1097  fprintf(stderr, " solver will continue anyway.\n");
1098  if(nhnil == mxhnil) {
1099  std::cerr << "lsoda -- above warning has been issued " << nhnil
1100  << " times, " << std::endl
1101  << " it will not be issued again for this problem" << std::endl;
1102  }
1103  }
1104  }
1105 
1106  /* Call stoda */
1107  stoda(neq, y, derivs, _data);
1108  if(kflag == 0) {
1109  /*
1110  Block f.
1111  The following block handles the case of a successful return from the
1112  core integrator ( kflag = 0 ).
1113  If a method switch was just made, record tsw, reset maxord,
1114  set jstart to -1 to signal stoda to complete the switch,
1115  and do extra printing of data if ixpr = 1.
1116  Then, in any case, check for stop conditions.
1117  */
1118  if (tn_ < tout) after_step(tn_, ++y.begin());
1119 
1120  init = 1;
1121  if(meth_ != mused) {
1122  tsw = tn_;
1123  maxord = mxordn;
1124  if(meth_ == 2)
1125  maxord = mxords;
1126  jstart = -1;
1127  if(ixpr) {
1128  if(meth_ == 2)
1129  std::cerr << "[lsoda] a switch to the stiff method has occurred "
1130  << std::endl;
1131  if(meth_ == 1)
1132  std::cerr << "[lsoda] a switch to the nonstiff method has occurred"
1133  << std::endl;
1134  }
1135  } /* end if ( meth_ != mused ) */
1136  /*
1137  itask = 1.
1138  If tout has been reached, interpolate.
1139  */
1140  if(1 == itask) {
1141  if((tn_ - tout) * h_ < 0.)
1142  continue;
1143 
1144  intdy(tout, 0, y, &iflag);
1145  *t = tout;
1146  *istate = 2;
1147  illin = 0;
1148  after_step(*t, ++y.begin());
1149  return;
1150  }
1151  /*
1152  itask = 2.
1153  */
1154  if(itask == 2) {
1155  successreturn(after_step, y, t, itask, ihit, tcrit, istate);
1156  return;
1157  }
1158  /*
1159  itask = 3.
1160  Jump to exit if tout was reached.
1161  */
1162  if(itask == 3) {
1163  if((tn_ - tout) * h_ >= 0.) {
1164  successreturn(after_step, y, t, itask, ihit, tcrit, istate);
1165  return;
1166  }
1167  continue;
1168  }
1169  /*
1170  itask = 4.
1171  See if tout or tcrit was reached. Adjust h_ if necessary.
1172  */
1173  if(itask == 4) {
1174  if((tn_ - tout) * h_ >= 0.) {
1175  intdy(tout, 0, y, &iflag);
1176  *t = tout;
1177  *istate = 2;
1178  illin = 0;
1179  after_step(*t, ++y.begin());
1180  return;
1181  }
1182  else {
1183  hmx = fabs(tn_) + fabs(h_);
1184  ihit = fabs(tn_ - tcrit) <= (100. * ETA * hmx);
1185  if(ihit) {
1186  successreturn(after_step, y, t, itask, ihit, tcrit, istate);
1187  return;
1188  }
1189  tnext = tn_ + h_ * (1. + 4. * ETA);
1190  if((tnext - tcrit) * h_ <= 0.)
1191  continue;
1192  h_ = (tcrit - tn_) * (1. - 4. * ETA);
1193  jstart = -2;
1194  continue;
1195  }
1196  } /* end if ( itask == 4 ) */
1197  /*
1198  itask = 5.
1199  See if tcrit was reached and jump to exit.
1200  */
1201  if(itask == 5) {
1202  hmx = fabs(tn_) + fabs(h_);
1203  ihit = fabs(tn_ - tcrit) <= (100. * ETA * hmx);
1204  successreturn(after_step, y, t, itask, ihit, tcrit, istate);
1205  return;
1206  }
1207 
1208  } /* end if ( kflag == 0 ) */
1209  /*
1210  kflag = -1, error test failed repeatedly or with fabs(h_) = hmin.
1211  kflag = -2, convergence failed repeatedly or with fabs(h_) = hmin.
1212  */
1213  if(kflag == -1 || kflag == -2) {
1214  fprintf(stderr, "lsoda -- at t = %g and step size h_ = %g, the\n", tn_, h_);
1215  if(kflag == -1) {
1216  fprintf(stderr, " error test failed repeatedly or\n");
1217  fprintf(stderr, " with fabs(h_) = hmin\n");
1218  *istate = -4;
1219  }
1220  if(kflag == -2) {
1221  fprintf(stderr, " corrector convergence failed repeatedly or\n");
1222  fprintf(stderr, " with fabs(h_) = hmin\n");
1223  *istate = -5;
1224  }
1225  big = 0.;
1226  imxer = 1;
1227  for(size_t i = 1; i <= n; i++) {
1228  size = fabs(acor[i]) * ewt[i];
1229  if(big < size) {
1230  big = size;
1231  imxer = i;
1232  }
1233  }
1234  terminate2(y, t);
1235  return;
1236  } /* end if ( kflag == -1 || kflag == -2 ) */
1237  } /* end while */
1238 } /* end lsoda */
LSODA::nje
size_t nje
Definition: lsoda.h:137
LSODA::nyh
size_t nyh
Definition: lsoda.h:139
LSODA::mxstep
size_t mxstep
Definition: lsoda.h:138
LSODA::hu
double hu
Definition: lsoda.h:142
LSODA::nst
size_t nst
Definition: lsoda.h:137
LSODA::ETA
double ETA
Definition: lsoda.h:112
LSODA::mu
size_t mu
Definition: lsoda.h:114
LSODA::miter
size_t miter
Definition: lsoda.h:130
LSODA::terminate
void terminate(int *istate)
Definition: lsoda.cpp:217
LSODA::yh_
std::vector< std::vector< double > > yh_
Definition: lsoda.h:154
LSODA::ccmax
double ccmax
Definition: lsoda.h:141
LSODA::tn_
double tn_
Definition: lsoda.h:142
LSODA::ewset
void ewset(const std::vector< double > &ycur)
Definition: lsoda.cpp:241
LSODA::meth_
size_t meth_
Definition: lsoda.h:135
LSODA::hmin
double hmin
Definition: lsoda.h:142
LSODA::nfe
size_t nfe
Definition: lsoda.h:137
LSODA::ewt
std::vector< double > ewt
Definition: lsoda.h:151
LSODA::nqu
size_t nqu
Definition: lsoda.h:137
LSODA::n
size_t n
Definition: lsoda.h:137
LSODA::vmnorm
double vmnorm(const size_t n, const std::vector< double > &v, const std::vector< double > &w)
Definition: lsoda.cpp:501
LSODA::ml
size_t ml
Definition: lsoda.h:114
LSODA::successreturn
void successreturn(AfterStep &after_step, std::vector< double > &y, double *t, int itask, int ihit, double tcrit, int *istate)
Definition: lsoda.tpp:578
LSODA::mord
std::array< size_t, 3 > mord
Definition: lsoda.h:120
LSODA::itol_
int itol_
Definition: lsoda.h:160
LSODA::ntrep
size_t ntrep
Definition: lsoda.h:139
LSODA::mused
size_t mused
Definition: lsoda.h:134
LSODA::jstart
int jstart
Definition: lsoda.h:132
LSODA::illin
size_t illin
Definition: lsoda.h:130
LSODA::ixpr
size_t ixpr
Definition: lsoda.h:134
LSODA::kflag
int kflag
Definition: lsoda.h:132
LSODA::h_
double h_
Definition: lsoda.h:141
LSODA::nslast
size_t nslast
Definition: lsoda.h:139
LSODA::nq
size_t nq
Definition: lsoda.h:137
LSODA::terminate2
void terminate2(std::vector< double > &y, double *t)
Definition: lsoda.cpp:230
LSODA::intdy
void intdy(double t, int k, std::vector< double > &dky, int *iflag)
Definition: lsoda.cpp:286
LSODA::maxcor
size_t maxcor
Definition: lsoda.h:130
LSODA::acor
std::vector< double > acor
Definition: lsoda.h:153
LSODA::mxhnil
size_t mxhnil
Definition: lsoda.h:138
LSODA::tsw
double tsw
Definition: lsoda.h:143
LSODA::sqrteta
double sqrteta
Definition: lsoda.h:115
LSODA::hmxi
double hmxi
Definition: lsoda.h:142
LSODA::mxordn
size_t mxordn
Definition: lsoda.h:134
LSODA::init
size_t init
Definition: lsoda.h:130
LSODA::resize_system
void resize_system(int nyh, int lenyh)
Definition: lsoda.cpp:837
LSODA::nhnil
size_t nhnil
Definition: lsoda.h:139
LSODA::msbp
size_t msbp
Definition: lsoda.h:130
LSODA::atol_
std::vector< double > atol_
Definition: lsoda.h:162
LSODA::stoda
void stoda(const size_t neq, std::vector< double > &y, Functor &derivs, void *_data)
Definition: lsoda.tpp:225
LSODA::maxord
size_t maxord
Definition: lsoda.h:130
LSODA::mxncf
size_t mxncf
Definition: lsoda.h:130
LSODA::imxer
size_t imxer
Definition: lsoda.h:114
LSODA::mxords
size_t mxords
Definition: lsoda.h:134
LSODA::jtyp
size_t jtyp
Definition: lsoda.h:134
LSODA::rtol_
std::vector< double > rtol_
Definition: lsoda.h:161
Here is the call graph for this function:
Here is the caller graph for this function:

◆ lsoda_update()

template<class Functor , class AfterStep >
void LSODA::lsoda_update ( Functor &  derivs,
AfterStep &  after_step,
const size_t  neq,
std::vector< double > &  y,
double *  t,
const double  tout,
void *const  _data,
double  rtol = 1e-6,
double  atol = 1e-6 
)

Simpler interface.

f System neq, size of system. y, init values of size neq yout, results vector for size neq+1, ignore yout[0] t, start time. tout, stop time. _data rtol, relative tolerance. atol, absolute tolerance.

Definition at line 1257 of file lsoda.tpp.

1260 {
1261  std::array<int, 7> iworks = {{0}};
1262  std::array<double, 4> rworks = {{0.0}};
1263 
1264  int itask, iopt, jt;
1265 
1266  // cout << "Debug : rtol " << rtol << ". atol " << atol << std::endl;
1267 
1268  itask = 1;
1269  iopt = 0;
1270  jt = 2;
1271 
1272  yout.resize(neq + 1);
1273 
1274  // Set the tolerance. We should do it only once.
1275  rtol_.resize(neq + 1, rtol);
1276  atol_.resize(neq + 1, atol);
1277  rtol_[0] = 0;
1278  atol_[0] = 0;
1279 
1280  // Fill-in values.
1281  for(size_t i = 1; i <= neq; i++)
1282  yout[i] = y[i - 1];
1283 
1284  lsoda(derivs, after_step, neq, yout, t, tout, itask, &iState, iopt, jt, iworks, rworks, _data);
1285 
1286  for (int i=0; i<y.size(); ++i) y[i] = yout[i+1];
1287 }
LSODA::iState
int iState
Definition: lsoda.h:165
LSODA::lsoda
void lsoda(Functor &derivs, AfterStep &after_step, const size_t neq, std::vector< double > &y, double *t, double tout, int itask, int *istate, int iopt, int jt, std::array< int, 7 > &iworks, std::array< double, 4 > &rworks, void *_data)
Definition: lsoda.tpp:606
LSODA::atol_
std::vector< double > atol_
Definition: lsoda.h:162
LSODA::yout
std::vector< double > yout
Definition: lsoda.h:167
LSODA::rtol_
std::vector< double > rtol_
Definition: lsoda.h:161
Here is the call graph for this function:
Here is the caller graph for this function:

◆ methodswitch()

void LSODA::methodswitch ( double  dsm,
double  pnorm,
double *  pdh,
double *  rh 
)

Definition at line 572 of file lsoda.cpp.

573 {
574  int lm1, lm1p1, lm2, lm2p1, nqm1, nqm2;
575  double rh1, rh2, rh1it, exm2, dm2, exm1, dm1, alpha, exsm;
576 
577  /*
578  We are current using an Adams method. Consider switching to bdf.
579  If the current order is greater than 5, assume the problem is
580  not stiff, and skip this section.
581  If the Lipschitz constant and error estimate are not polluted
582  by roundoff, perform the usual test.
583  Otherwise, switch to the bdf methods if the last step was
584  restricted to insure stability ( irflag = 1 ), and stay with Adams
585  method if not. When switching to bdf with polluted error estimates,
586  in the absence of other information, double the step size.
587 
588  When the estimates are ok, we make the usual test by computing
589  the step size we could have (ideally) used on this step,
590  with the current (Adams) method, and also that for the bdf.
591  If nq > mxords, we consider changing to order mxords on switching.
592  Compare the two step sizes to decide whether to switch.
593  The step size advantage must be at least ratio = 5 to switch.
594  */
595  if(meth_ == 1) {
596  if(nq > 5)
597  return;
598  if(dsm <= (100. * pnorm * ETA) || pdest == 0.) {
599  if(irflag == 0)
600  return;
601  rh2 = 2.;
602  nqm2 = min(nq, mxords);
603  }
604  else {
605  exsm = 1. / (double)l;
606  rh1 = 1. / (1.2 * pow(dsm, exsm) + 0.0000012);
607  rh1it = 2. * rh1;
608  *pdh = pdlast * fabs(h_);
609  if((*pdh * rh1) > 0.00001)
610  rh1it = sm1[nq] / *pdh;
611  rh1 = min(rh1, rh1it);
612  if(nq > mxords) {
613  nqm2 = mxords;
614  lm2 = mxords + 1;
615  exm2 = 1. / (double)lm2;
616  lm2p1 = lm2 + 1;
617  dm2 = vmnorm(n, yh_[lm2p1], ewt) / cm2[mxords];
618  rh2 = 1. / (1.2 * pow(dm2, exm2) + 0.0000012);
619  }
620  else {
621  dm2 = dsm * (cm1[nq] / cm2[nq]);
622  rh2 = 1. / (1.2 * pow(dm2, exsm) + 0.0000012);
623  nqm2 = nq;
624  }
625  if(rh2 < ratio * rh1)
626  return;
627  }
628  /*
629  The method switch test passed. Reset relevant quantities for bdf.
630  */
631  *rh = rh2;
632  icount = 20;
633  meth_ = 2;
634  miter = jtyp;
635  pdlast = 0.;
636  nq = nqm2;
637  l = nq + 1;
638  return;
639  } /* end if ( meth_ == 1 ) */
640 
641  /*
642  We are currently using a bdf method, considering switching to Adams.
643  Compute the step size we could have (ideally) used on this step,
644  with the current (bdf) method, and also that for the Adams.
645  If nq > mxordn, we consider changing to order mxordn on switching.
646  Compare the two step sizes to decide whether to switch.
647  The step size advantage must be at least 5/ratio = 1 to switch.
648  If the step size for Adams would be so small as to cause
649  roundoff pollution, we stay with bdf.
650  */
651  exsm = 1. / (double)l;
652  if(mxordn < nq) {
653  nqm1 = mxordn;
654  lm1 = mxordn + 1;
655  exm1 = 1. / (double)lm1;
656  lm1p1 = lm1 + 1;
657  dm1 = vmnorm(n, yh_[lm1p1], ewt) / cm1[mxordn];
658  rh1 = 1. / (1.2 * pow(dm1, exm1) + 0.0000012);
659  }
660  else {
661  dm1 = dsm * (cm2[nq] / cm1[nq]);
662  rh1 = 1. / (1.2 * pow(dm1, exsm) + 0.0000012);
663  nqm1 = nq;
664  exm1 = exsm;
665  }
666  rh1it = 2. * rh1;
667  *pdh = pdnorm * fabs(h_);
668  if((*pdh * rh1) > 0.00001)
669  rh1it = sm1[nqm1] / *pdh;
670  rh1 = min(rh1, rh1it);
671  rh2 = 1. / (1.2 * pow(dsm, exsm) + 0.0000012);
672  if((rh1 * ratio) < (5. * rh2))
673  return;
674  alpha = max(0.001, rh1);
675  dm1 *= pow(alpha, exm1);
676  if(dm1 <= 1000. * ETA * pnorm)
677  return;
678  /*
679  The switch test passed. Reset relevant quantities for Adams.
680  */
681  *rh = rh1;
682  icount = 20;
683  meth_ = 1;
684  miter = 0;
685  pdlast = 0.;
686  nq = nqm1;
687  l = nq + 1;
688 } /* end methodswitch */
LSODA::ETA
double ETA
Definition: lsoda.h:112
LSODA::pdlast
double pdlast
Definition: lsoda.h:148
LSODA::ratio
double ratio
Definition: lsoda.h:148
LSODA::miter
size_t miter
Definition: lsoda.h:130
LSODA::yh_
std::vector< std::vector< double > > yh_
Definition: lsoda.h:154
LSODA::sm1
std::array< double, 13 > sm1
Definition: lsoda.h:121
LSODA::meth_
size_t meth_
Definition: lsoda.h:135
LSODA::ewt
std::vector< double > ewt
Definition: lsoda.h:151
LSODA::n
size_t n
Definition: lsoda.h:137
LSODA::vmnorm
double vmnorm(const size_t n, const std::vector< double > &v, const std::vector< double > &w)
Definition: lsoda.cpp:501
LSODA::irflag
int irflag
Definition: lsoda.h:149
LSODA::pdest
double pdest
Definition: lsoda.h:148
LSODA::icount
int icount
Definition: lsoda.h:149
LSODA::cm1
std::array< double, 13 > cm1
Definition: lsoda.h:124
LSODA::h_
double h_
Definition: lsoda.h:141
LSODA::pdnorm
double pdnorm
Definition: lsoda.h:143
LSODA::nq
size_t nq
Definition: lsoda.h:137
LSODA::mxordn
size_t mxordn
Definition: lsoda.h:134
LSODA::l
size_t l
Definition: lsoda.h:130
LSODA::mxords
size_t mxords
Definition: lsoda.h:134
LSODA::cm2
std::array< double, 6 > cm2
Definition: lsoda.h:125
LSODA::jtyp
size_t jtyp
Definition: lsoda.h:134
Here is the caller graph for this function:

◆ orderswitch()

void LSODA::orderswitch ( double *  rhup,
double  dsm,
double *  pdh,
double *  rh,
size_t *  orderflag 
)

Definition at line 716 of file lsoda.cpp.

718 {
719  size_t newq = 0;
720  double exsm, rhdn, rhsm, ddn, exdn, r;
721 
722  *orderflag = 0;
723 
724  exsm = 1. / (double)l;
725  rhsm = 1. / (1.2 * pow(dsm, exsm) + 0.0000012);
726 
727  rhdn = 0.;
728  if(nq != 1) {
729  ddn = vmnorm(n, yh_[l], ewt) / tesco[nq][1];
730  exdn = 1. / (double)nq;
731  rhdn = 1. / (1.3 * pow(ddn, exdn) + 0.0000013);
732  }
733  /*
734  If meth_ = 1, limit rh accordinfg to the stability region also.
735  */
736  if(meth_ == 1) {
737  *pdh = max(fabs(h_) * pdlast, 0.000001);
738  if(l < lmax)
739  *rhup = min(*rhup, sm1[l] / *pdh);
740  rhsm = min(rhsm, sm1[nq] / *pdh);
741  if(nq > 1)
742  rhdn = min(rhdn, sm1[nq - 1] / *pdh);
743  pdest = 0.;
744  }
745  if(rhsm >= *rhup) {
746  if(rhsm >= rhdn) {
747  newq = nq;
748  *rh = rhsm;
749  }
750  else {
751  newq = nq - 1;
752  *rh = rhdn;
753  if(kflag < 0 && *rh > 1.)
754  *rh = 1.;
755  }
756  }
757  else {
758  if(*rhup <= rhdn) {
759  newq = nq - 1;
760  *rh = rhdn;
761  if(kflag < 0 && *rh > 1.)
762  *rh = 1.;
763  }
764  else {
765  *rh = *rhup;
766  if(*rh >= 1.1) {
767  r = el[l] / (double)l;
768  nq = l;
769  l = nq + 1;
770  for(size_t i = 1; i <= n; i++)
771  yh_[l][i] = acor[i] * r;
772 
773  *orderflag = 2;
774  return;
775  }
776  else {
777  ialth = 3;
778  return;
779  }
780  }
781  }
782  /*
783  If meth_ = 1 and h_ is restricted by stability, bypass 10 percent test.
784  */
785  if(1 == meth_) {
786  if((*rh * *pdh * 1.00001) < sm1[newq])
787  if(kflag == 0 && *rh < 1.1) {
788  ialth = 3;
789  return;
790  }
791  }
792  else {
793  if(kflag == 0 && *rh < 1.1) {
794  ialth = 3;
795  return;
796  }
797  }
798  if(kflag <= -2)
799  *rh = min(*rh, 0.2);
800  /*
801  If there is a change of order, reset nq, l, and the coefficients.
802  In any case h_ is reset according to rh and the yh_ array is rescaled.
803  Then exit or redo the step.
804  */
805  if(newq == nq) {
806  *orderflag = 1;
807  return;
808  }
809  nq = newq;
810  l = nq + 1;
811  *orderflag = 2;
812 
813 } /* end orderswitch */
LSODA::pdlast
double pdlast
Definition: lsoda.h:148
LSODA::yh_
std::vector< std::vector< double > > yh_
Definition: lsoda.h:154
LSODA::sm1
std::array< double, 13 > sm1
Definition: lsoda.h:121
LSODA::meth_
size_t meth_
Definition: lsoda.h:135
LSODA::ewt
std::vector< double > ewt
Definition: lsoda.h:151
LSODA::n
size_t n
Definition: lsoda.h:137
LSODA::vmnorm
double vmnorm(const size_t n, const std::vector< double > &v, const std::vector< double > &w)
Definition: lsoda.cpp:501
LSODA::pdest
double pdest
Definition: lsoda.h:148
LSODA::ialth
size_t ialth
Definition: lsoda.h:146
LSODA::tesco
std::array< std::array< double, 4 >, 13 > tesco
Definition: lsoda.h:128
LSODA::kflag
int kflag
Definition: lsoda.h:132
LSODA::h_
double h_
Definition: lsoda.h:141
LSODA::nq
size_t nq
Definition: lsoda.h:137
LSODA::acor
std::vector< double > acor
Definition: lsoda.h:153
LSODA::lmax
size_t lmax
Definition: lsoda.h:146
LSODA::l
size_t l
Definition: lsoda.h:130
LSODA::el
std::array< double, 14 > el
Definition: lsoda.h:123
Here is the caller graph for this function:

◆ prja()

template<class Functor >
void LSODA::prja ( const size_t  neq,
std::vector< double > &  y,
Functor &  derivs,
void *  _data 
)

Definition at line 6 of file lsoda.tpp.

8 {
9  (void)neq;
10 
11  size_t i = 0, ier = 0, j = 0;
12  double fac = 0.0, hl0 = 0.0, r = 0.0, r0 = 0.0, yj = 0.0;
13  /*
14  prja is called by stoda to compute and process the matrix
15  P = I - h_ * el[1] * J, where J is an approximation to the Jacobian.
16  Here J is computed by finite differencing.
17  J, scaled by -h_ * el[1], is stored in wm_. Then the norm of J ( the
18  matrix norm consistent with the weighted max-norm on vectors given
19  by vmnorm ) is computed, and J is overwritten by P. P is then
20  subjected to LU decomposition in preparation for later solution
21  of linear systems with p as coefficient matrix. This is done
22  by dgefa if miter = 2, and by dgbfa if miter = 5.
23  */
24  nje++;
25  ierpj = 0;
26  jcur = 1;
27  hl0 = h_ * el0;
28  /*
29  If miter = 2, make n calls to f to approximate J.
30  */
31  if(miter != 2) {
32  fprintf(stderr, "[prja] miter != 2\n");
33  return;
34  }
35  if(miter == 2) {
36  fac = vmnorm(n, savf, ewt);
37  r0 = 1000. * fabs(h_) * ETA * ((double)n) * fac;
38  if(r0 == 0.)
39  r0 = 1.;
40  for(j = 1; j <= n; j++) {
41  yj = y[j];
42  r = std::max(sqrteta * fabs(yj), r0 / ewt[j]);
43  y[j] += r;
44  fac = -hl0 / r;
45  derivs(tn_, ++y.begin(), ++acor.begin(), _data);
46  for(i = 1; i <= n; i++)
47  wm_[i][j] = (acor[i] - savf[i]) * fac;
48  y[j] = yj;
49  }
50  nfe += n;
51  /*
52  Compute norm of Jacobian.
53  */
54  pdnorm = fnorm(n, wm_, ewt) / fabs(hl0);
55  /*
56  Add identity matrix.
57  */
58  for(i = 1; i <= n; i++)
59  wm_[i][i] += 1.;
60  /*
61  Do LU decomposition on P.
62  */
63  dgefa(wm_, n, ipvt, &ier);
64  if(ier != 0)
65  ierpj = 1;
66  return;
67  }
68 } /* end prja */
LSODA::nje
size_t nje
Definition: lsoda.h:137
LSODA::ETA
double ETA
Definition: lsoda.h:112
LSODA::savf
std::vector< double > savf
Definition: lsoda.h:152
LSODA::miter
size_t miter
Definition: lsoda.h:130
LSODA::jcur
size_t jcur
Definition: lsoda.h:130
LSODA::el0
double el0
Definition: lsoda.h:141
LSODA::wm_
std::vector< std::vector< double > > wm_
Definition: lsoda.h:155
LSODA::tn_
double tn_
Definition: lsoda.h:142
LSODA::ierpj
size_t ierpj
Definition: lsoda.h:130
LSODA::nfe
size_t nfe
Definition: lsoda.h:137
LSODA::ewt
std::vector< double > ewt
Definition: lsoda.h:151
LSODA::n
size_t n
Definition: lsoda.h:137
LSODA::vmnorm
double vmnorm(const size_t n, const std::vector< double > &v, const std::vector< double > &w)
Definition: lsoda.cpp:501
LSODA::fnorm
double fnorm(int n, const std::vector< std::vector< double >> &a, const std::vector< double > &w)
Definition: lsoda.cpp:509
LSODA::h_
double h_
Definition: lsoda.h:141
LSODA::pdnorm
double pdnorm
Definition: lsoda.h:143
LSODA::acor
std::vector< double > acor
Definition: lsoda.h:153
LSODA::sqrteta
double sqrteta
Definition: lsoda.h:115
LSODA::dgefa
void dgefa(std::vector< std::vector< double >> &a, const size_t n, std::vector< int > &ipvt, size_t *const info)
Definition: lsoda.cpp:161
LSODA::ipvt
std::vector< int > ipvt
Definition: lsoda.h:157
Here is the call graph for this function:
Here is the caller graph for this function:

◆ resetcoeff()

void LSODA::resetcoeff ( void  )

Definition at line 815 of file lsoda.cpp.

820 {
821  array<double, 14> ep1;
822 
823  ep1 = elco[nq];
824  for(size_t i = 1; i <= l; i++)
825  el[i] = ep1[i];
826  rc = rc * el[1] / el0;
827  el0 = el[1];
828  conit = 0.5 / (double)(nq + 2);
829 }
LSODA::elco
std::array< std::array< double, 14 >, 13 > elco
Definition: lsoda.h:127
LSODA::rc
double rc
Definition: lsoda.h:142
LSODA::el0
double el0
Definition: lsoda.h:141
LSODA::conit
double conit
Definition: lsoda.h:144
LSODA::nq
size_t nq
Definition: lsoda.h:137
LSODA::l
size_t l
Definition: lsoda.h:130
LSODA::el
std::array< double, 14 > el
Definition: lsoda.h:123
Here is the caller graph for this function:

◆ resize_system()

void LSODA::resize_system ( int  nyh,
int  lenyh 
)

Definition at line 837 of file lsoda.cpp.

837  {
838  // yh_ and wm_ need a clear() because otherwise, only additional elements will have the new length nyh+1
839  yh_.clear(); yh_.resize(lenyh + 1, std::vector<double>(nyh + 1, 0.0));
840  wm_.clear(); wm_.resize(nyh + 1, std::vector<double>(nyh + 1, 0.0));
841  ewt.resize(1 + nyh, 0);
842  savf.resize(1 + nyh, 0);
843  acor.resize(nyh + 1, 0.0);
844  ipvt.resize(nyh + 1, 0.0);
845 }
LSODA::nyh
size_t nyh
Definition: lsoda.h:139
LSODA::savf
std::vector< double > savf
Definition: lsoda.h:152
LSODA::yh_
std::vector< std::vector< double > > yh_
Definition: lsoda.h:154
LSODA::wm_
std::vector< std::vector< double > > wm_
Definition: lsoda.h:155
LSODA::ewt
std::vector< double > ewt
Definition: lsoda.h:151
LSODA::acor
std::vector< double > acor
Definition: lsoda.h:153
LSODA::ipvt
std::vector< int > ipvt
Definition: lsoda.h:157
Here is the caller graph for this function:

◆ scaleh()

void LSODA::scaleh ( double *  rh,
double *  pdh 
)

Definition at line 458 of file lsoda.cpp.

459 {
460  double r;
461  /*
462  If h_ is being changed, the h_ ratio rh is checked against rmax, hmin,
463  and hmxi, and the yh_ array is rescaled. ialth is set to l = nq + 1
464  to prevent a change of h_ for that many steps, unless forced by a
465  convergence or error test failure.
466  */
467  *rh = min(*rh, rmax);
468  *rh = *rh / max(1., fabs(h_) * hmxi * *rh);
469  /*
470  If meth_ = 1, also restrict the new step size by the stability region.
471  If this reduces h_, set irflag to 1 so that if there are roundoff
472  problems later, we can assume that is the cause of the trouble.
473  */
474  if(meth_ == 1) {
475  irflag = 0;
476  *pdh = max(fabs(h_) * pdlast, 0.000001);
477  if((*rh * *pdh * 1.00001) >= sm1[nq]) {
478  *rh = sm1[nq] / *pdh;
479  irflag = 1;
480  }
481  }
482  r = 1.;
483  for(size_t j = 2; j <= l; j++) {
484  r *= *rh;
485  for(size_t i = 1; i <= n; i++)
486  yh_[j][i] *= r;
487  }
488  h_ *= *rh;
489  rc *= *rh;
490  ialth = l;
491 
492 } /* end scaleh */
LSODA::pdlast
double pdlast
Definition: lsoda.h:148
LSODA::rc
double rc
Definition: lsoda.h:142
LSODA::yh_
std::vector< std::vector< double > > yh_
Definition: lsoda.h:154
LSODA::sm1
std::array< double, 13 > sm1
Definition: lsoda.h:121
LSODA::meth_
size_t meth_
Definition: lsoda.h:135
LSODA::n
size_t n
Definition: lsoda.h:137
LSODA::irflag
int irflag
Definition: lsoda.h:149
LSODA::ialth
size_t ialth
Definition: lsoda.h:146
LSODA::h_
double h_
Definition: lsoda.h:141
LSODA::nq
size_t nq
Definition: lsoda.h:137
LSODA::hmxi
double hmxi
Definition: lsoda.h:142
LSODA::l
size_t l
Definition: lsoda.h:130
LSODA::rmax
double rmax
Definition: lsoda.h:144
Here is the caller graph for this function:

◆ set_istate()

void LSODA::set_istate ( int  n)

Definition at line 857 of file lsoda.cpp.

857  {
858  iState = n;
859 }
LSODA::iState
int iState
Definition: lsoda.h:165
LSODA::n
size_t n
Definition: lsoda.h:137
Here is the caller graph for this function:

◆ solsy()

void LSODA::solsy ( std::vector< double > &  y)

Definition at line 560 of file lsoda.cpp.

561 {
562  iersl = 0;
563  if(miter != 2) {
564  printf("solsy -- miter != 2\n");
565  return;
566  }
567  if(miter == 2)
568  dgesl(wm_, n, ipvt, y, 0);
569  return;
570 }
LSODA::miter
size_t miter
Definition: lsoda.h:130
LSODA::wm_
std::vector< std::vector< double > > wm_
Definition: lsoda.h:155
LSODA::iersl
size_t iersl
Definition: lsoda.h:130
LSODA::n
size_t n
Definition: lsoda.h:137
LSODA::dgesl
void dgesl(const std::vector< std::vector< double >> &a, const size_t n, std::vector< int > &ipvt, std::vector< double > &b, const size_t job)
Definition: lsoda.cpp:105
LSODA::ipvt
std::vector< int > ipvt
Definition: lsoda.h:157
Here is the caller graph for this function:

◆ stoda()

template<class Functor >
void LSODA::stoda ( const size_t  neq,
std::vector< double > &  y,
Functor &  derivs,
void *  _data 
)

Definition at line 225 of file lsoda.tpp.

227 {
228  assert(neq + 1 == y.size());
229 
230  size_t corflag = 0, orderflag = 0;
231  size_t i = 0, i1 = 0, j = 0, m = 0, ncf = 0;
232  double del = 0.0, delp = 0.0, dsm = 0.0, dup = 0.0, exup = 0.0, r = 0.0, rh = 0.0,
233  rhup = 0.0, told = 0.0;
234  double pdh = 0.0, pnorm = 0.0;
235 
236  /*
237  stoda performs one step of the integration of an initial value
238  problem for a system of ordinary differential equations.
239  Note.. stoda is independent of the value of the iteration method
240  indicator miter, when this is != 0, and hence is independent
241  of the type of chord method used, or the Jacobian structure.
242  Communication with stoda is done with the following variables:
243 
244  jstart = an integer used for input only, with the following
245  values and meanings:
246 
247  0 perform the first step,
248  > 0 take a new step continuing from the last,
249  -1 take the next step with a new value of h_,
250  n, meth_, miter, and/or matrix parameters.
251  -2 take the next step with a new value of h_,
252  but with other inputs unchanged.
253 
254  kflag = a completion code with the following meanings:
255 
256  0 the step was successful,
257  -1 the requested error could not be achieved,
258  -2 corrector convergence could not be achieved,
259  -3 fatal error in prja or solsy.
260 
261  miter = corrector iteration method:
262 
263  0 functional iteration,
264  >0 a chord method corresponding to jacobian type jt.
265 
266  */
267  kflag = 0;
268  told = tn_;
269  ncf = 0;
270  ierpj = 0;
271  iersl = 0;
272  jcur = 0;
273  delp = 0.;
274 
275  /*
276  On the first call, the order is set to 1, and other variables are
277  initialized. rmax is the maximum ratio by which h_ can be increased
278  in a single step. It is initially 1.e4 to compensate for the small
279  initial h_, but then is normally equal to 10. If a filure occurs
280  (in corrector convergence or error test), rmax is set at 2 for
281  the next increase.
282  cfode is called to get the needed coefficients for both methods.
283  */
284  if(jstart == 0) {
285  lmax = maxord + 1;
286  nq = 1;
287  l = 2;
288  ialth = 2;
289  rmax = 10000.;
290  rc = 0.;
291  el0 = 1.;
292  crate = 0.7;
293  hold = h_;
294  nslp = 0;
295  ipup = miter;
296  /*
297  Initialize switching parameters. meth_ = 1 is assumed initially.
298  */
299  icount = 20;
300  irflag = 0;
301  pdest = 0.;
302  pdlast = 0.;
303  ratio = 5.;
304  cfode(2);
305  for(i = 1; i <= 5; i++)
306  cm2[i] = tesco[i][2] * elco[i][i + 1];
307  cfode(1);
308  for(i = 1; i <= 12; i++)
309  cm1[i] = tesco[i][2] * elco[i][i + 1];
310  resetcoeff();
311  } /* end if ( jstart == 0 ) */
312  /*
313  The following block handles preliminaries needed when jstart = -1.
314  ipup is set to miter to force a matrix update.
315  If an order increase is about to be considered ( ialth = 1 ),
316  ialth is reset to 2 to postpone consideration one more step.
317  If the caller has changed meth_, cfode is called to reset
318  the coefficients of the method.
319  If h_ is to be changed, yh_ must be rescaled.
320  If h_ or meth_ is being changed, ialth is reset to l = nq + 1
321  to prevent further changes in h_ for that many steps.
322  */
323  if(jstart == -1) {
324  ipup = miter;
325  lmax = maxord + 1;
326  if(ialth == 1)
327  ialth = 2;
328  if(meth_ != mused) {
329  cfode(meth_);
330  ialth = l;
331  resetcoeff();
332  }
333  if(h_ != hold) {
334  rh = h_ / hold;
335  h_ = hold;
336  scaleh(&rh, &pdh);
337  }
338  } /* if ( jstart == -1 ) */
339  if(jstart == -2) {
340  if(h_ != hold) {
341  rh = h_ / hold;
342  h_ = hold;
343  scaleh(&rh, &pdh);
344  }
345  } /* if ( jstart == -2 ) */
346 
347  /*
348  Prediction.
349  This section computes the predicted values by effectively
350  multiplying the yh_ array by the pascal triangle matrix.
351  rc is the ratio of new to old values of the coefficient h_ * el[1].
352  When rc differs from 1 by more than ccmax, ipup is set to miter
353  to force pjac to be called, if a jacobian is involved.
354  In any case, prja is called at least every msbp steps.
355  */
356  while(1) {
357  while(1) {
358  if(fabs(rc - 1.) > ccmax)
359  ipup = miter;
360  if(nst >= nslp + msbp)
361  ipup = miter;
362  tn_ += h_;
363  for(size_t j = nq; j >= 1; j--)
364  for(size_t i1 = j; i1 <= nq; i1++)
365  for(i = 1; i <= n; i++)
366  yh_[i1][i] += yh_[i1 + 1][i];
367 
368  pnorm = vmnorm(n, yh_[1], ewt);
369  correction(
370  neq, y, derivs, &corflag, pnorm, &del, &delp, &told, &ncf, &rh, &m, _data);
371  if(corflag == 0)
372  break;
373  if(corflag == 1) {
374  rh = std::max(rh, hmin / fabs(h_));
375  scaleh(&rh, &pdh);
376  continue;
377  }
378  if(corflag == 2) {
379  kflag = -2;
380  hold = h_;
381  jstart = 1;
382  return;
383  }
384  } /* end inner while ( corrector loop ) */
385 
386  /*
387  The corrector has converged. jcur is set to 0
388  to signal that the Jacobian involved may need updating later.
389  The local error test is done now.
390  */
391  jcur = 0;
392  if(m == 0)
393  dsm = del / tesco[nq][2];
394  if(m > 0)
395  dsm = vmnorm(n, acor, ewt) / tesco[nq][2];
396 
397  if(dsm <= 1.) {
398  /*
399  After a successful step, update the yh_ array.
400  Decrease icount by 1, and if it is -1, consider switching methods.
401  If a method switch is made, reset various parameters,
402  rescale the yh_ array, and exit. If there is no switch,
403  consider changing h_ if ialth = 1. Otherwise decrease ialth by 1.
404  If ialth is then 1 and nq < maxord, then acor is saved for
405  use in a possible order increase on the next step.
406  If a change in h_ is considered, an increase or decrease in order
407  by one is considered also. A change in h_ is made only if it is by
408  a factor of at least 1.1. If not, ialth is set to 3 to prevent
409  testing for that many steps.
410  */
411  kflag = 0;
412  nst++;
413  hu = h_;
414  nqu = nq;
415  mused = meth_;
416  for(size_t j = 1; j <= l; j++) {
417  r = el[j];
418  for(i = 1; i <= n; i++)
419  yh_[j][i] += r * acor[i];
420  }
421  icount--;
422  if(icount < 0) {
423  methodswitch(dsm, pnorm, &pdh, &rh);
424  if(meth_ != mused) {
425  rh = std::max(rh, hmin / fabs(h_));
426  scaleh(&rh, &pdh);
427  rmax = 10.;
428  endstoda();
429  break;
430  }
431  }
432  /*
433  No method switch is being made. Do the usual step/order selection.
434  */
435  ialth--;
436  if(ialth == 0) {
437  rhup = 0.;
438  if(l != lmax) {
439  for(i = 1; i <= n; i++)
440  savf[i] = acor[i] - yh_[lmax][i];
441  dup = vmnorm(n, savf, ewt) / tesco[nq][3];
442  exup = 1. / (double)(l + 1);
443  rhup = 1. / (1.4 * pow(dup, exup) + 0.0000014);
444  }
445 
446  orderswitch(&rhup, dsm, &pdh, &rh, &orderflag);
447 
448  /*
449  No change in h_ or nq.
450  */
451  if(orderflag == 0) {
452  endstoda();
453  break;
454  }
455  /*
456  h_ is changed, but not nq.
457  */
458  if(orderflag == 1) {
459  rh = std::max(rh, hmin / fabs(h_));
460  scaleh(&rh, &pdh);
461  rmax = 10.;
462  endstoda();
463  break;
464  }
465  /*
466  both nq and h_ are changed.
467  */
468  if(orderflag == 2) {
469  resetcoeff();
470  rh = std::max(rh, hmin / fabs(h_));
471  scaleh(&rh, &pdh);
472  rmax = 10.;
473  endstoda();
474  break;
475  }
476  } /* end if ( ialth == 0 ) */
477  if(ialth > 1 || l == lmax) {
478  endstoda();
479  break;
480  }
481 
482  for(size_t i = 1; i <= n; i++)
483  yh_[lmax][i] = acor[i];
484 
485  endstoda();
486  break;
487  }
488  /* end if ( dsm <= 1. ) */
489  /*
490  The error test failed. kflag keeps track of multiple failures.
491  Restore tn_ and the yh_ array to their previous values, and prepare
492  to try the step again. Compute the optimum step size for this or
493  one lower. After 2 or more failures, h_ is forced to decrease
494  by a factor of 0.2 or less.
495  */
496  else {
497  kflag--;
498  tn_ = told;
499  for(j = nq; j >= 1; j--) {
500  for(i1 = j; i1 <= nq; i1++)
501  for(i = 1; i <= n; i++)
502  yh_[i1][i] -= yh_[i1 + 1][i];
503  }
504  rmax = 2.;
505  if(fabs(h_) <= hmin * 1.00001) {
506  kflag = -1;
507  hold = h_;
508  jstart = 1;
509  break;
510  }
511  if(kflag > -3) {
512  rhup = 0.;
513  orderswitch(&rhup, dsm, &pdh, &rh, &orderflag);
514  if(orderflag == 1 || orderflag == 0) {
515  if(orderflag == 0)
516  rh = std::min(rh, 0.2);
517  rh = std::max(rh, hmin / fabs(h_));
518  scaleh(&rh, &pdh);
519  }
520  if(orderflag == 2) {
521  resetcoeff();
522  rh = std::max(rh, hmin / fabs(h_));
523  scaleh(&rh, &pdh);
524  }
525  continue;
526  }
527  /* if ( kflag > -3 ) */
528  /*
529  Control reaches this section if 3 or more failures have occurred.
530  If 10 failures have occurred, exit with kflag = -1.
531  It is assumed that the derivatives that have accumulated in the
532  yh_ array have errors of the wrong order. Hence the first
533  derivative is recomputed, and the order is set to 1. Then
534  h_ is reduced by a factor of 10, and the step is retried,
535  until it succeeds or h_ reaches hmin.
536  */
537  else {
538  if(kflag == -10) {
539  kflag = -1;
540  hold = h_;
541  jstart = 1;
542  break;
543  }
544  else {
545  rh = 0.1;
546  rh = std::max(hmin / fabs(h_), rh);
547  h_ *= rh;
548  for(i = 1; i <= n; i++)
549  y[i] = yh_[1][i];
550  derivs(tn_, ++y.begin(), ++savf.begin(), _data);
551  nfe++;
552  for(i = 1; i <= n; i++)
553  yh_[2][i] = h_ * savf[i];
554  ipup = miter;
555  ialth = 5;
556  if(nq == 1)
557  continue;
558  nq = 1;
559  l = 2;
560  resetcoeff();
561  continue;
562  }
563  } /* end else -- kflag <= -3 */
564  } /* end error failure handling */
565  } /* end outer while */
566 
567 } /* end stoda */
LSODA::orderswitch
void orderswitch(double *rhup, double dsm, double *pdh, double *rh, size_t *orderflag)
Definition: lsoda.cpp:716
LSODA::hu
double hu
Definition: lsoda.h:142
LSODA::nst
size_t nst
Definition: lsoda.h:137
LSODA::pdlast
double pdlast
Definition: lsoda.h:148
LSODA::ratio
double ratio
Definition: lsoda.h:148
LSODA::savf
std::vector< double > savf
Definition: lsoda.h:152
LSODA::ipup
size_t ipup
Definition: lsoda.h:146
LSODA::elco
std::array< std::array< double, 14 >, 13 > elco
Definition: lsoda.h:127
LSODA::crate
double crate
Definition: lsoda.h:144
LSODA::rc
double rc
Definition: lsoda.h:142
LSODA::miter
size_t miter
Definition: lsoda.h:130
LSODA::yh_
std::vector< std::vector< double > > yh_
Definition: lsoda.h:154
LSODA::jcur
size_t jcur
Definition: lsoda.h:130
LSODA::el0
double el0
Definition: lsoda.h:141
LSODA::ccmax
double ccmax
Definition: lsoda.h:141
LSODA::tn_
double tn_
Definition: lsoda.h:142
LSODA::scaleh
void scaleh(double *rh, double *pdh)
Definition: lsoda.cpp:458
LSODA::meth_
size_t meth_
Definition: lsoda.h:135
LSODA::hmin
double hmin
Definition: lsoda.h:142
LSODA::ierpj
size_t ierpj
Definition: lsoda.h:130
LSODA::iersl
size_t iersl
Definition: lsoda.h:130
LSODA::nfe
size_t nfe
Definition: lsoda.h:137
LSODA::correction
void correction(const size_t neq, std::vector< double > &y, Functor &derivs, size_t *corflag, double pnorm, double *del, double *delp, double *told, size_t *ncf, double *rh, size_t *m, void *_data)
Definition: lsoda.tpp:78
LSODA::ewt
std::vector< double > ewt
Definition: lsoda.h:151
LSODA::nqu
size_t nqu
Definition: lsoda.h:137
LSODA::n
size_t n
Definition: lsoda.h:137
LSODA::vmnorm
double vmnorm(const size_t n, const std::vector< double > &v, const std::vector< double > &w)
Definition: lsoda.cpp:501
LSODA::irflag
int irflag
Definition: lsoda.h:149
LSODA::pdest
double pdest
Definition: lsoda.h:148
LSODA::ialth
size_t ialth
Definition: lsoda.h:146
LSODA::icount
int icount
Definition: lsoda.h:149
LSODA::tesco
std::array< std::array< double, 4 >, 13 > tesco
Definition: lsoda.h:128
LSODA::cm1
std::array< double, 13 > cm1
Definition: lsoda.h:124
LSODA::mused
size_t mused
Definition: lsoda.h:134
LSODA::jstart
int jstart
Definition: lsoda.h:132
LSODA::nslp
size_t nslp
Definition: lsoda.h:147
LSODA::cfode
void cfode(int meth_)
Definition: lsoda.cpp:331
LSODA::kflag
int kflag
Definition: lsoda.h:132
LSODA::h_
double h_
Definition: lsoda.h:141
LSODA::nq
size_t nq
Definition: lsoda.h:137
LSODA::endstoda
void endstoda(void)
Definition: lsoda.cpp:694
LSODA::acor
std::vector< double > acor
Definition: lsoda.h:153
LSODA::resetcoeff
void resetcoeff(void)
Definition: lsoda.cpp:815
LSODA::methodswitch
void methodswitch(double dsm, double pnorm, double *pdh, double *rh)
Definition: lsoda.cpp:572
LSODA::lmax
size_t lmax
Definition: lsoda.h:146
LSODA::msbp
size_t msbp
Definition: lsoda.h:130
LSODA::hold
double hold
Definition: lsoda.h:144
LSODA::l
size_t l
Definition: lsoda.h:130
LSODA::maxord
size_t maxord
Definition: lsoda.h:130
LSODA::cm2
std::array< double, 6 > cm2
Definition: lsoda.h:125
LSODA::el
std::array< double, 14 > el
Definition: lsoda.h:123
LSODA::rmax
double rmax
Definition: lsoda.h:144
Here is the call graph for this function:
Here is the caller graph for this function:

◆ successreturn()

template<class AfterStep >
void LSODA::successreturn ( AfterStep &  after_step,
std::vector< double > &  y,
double *  t,
int  itask,
int  ihit,
double  tcrit,
int *  istate 
)

Definition at line 578 of file lsoda.tpp.

580 {
581  for(size_t i = 1; i <= n; i++)
582  y[i] = yh_[1][i];
583  *t = tn_;
584  if(itask == 4 || itask == 5)
585  if(ihit)
586  *t = tcrit;
587  *istate = 2;
588  illin = 0;
589 
590  after_step(*t, ++y.begin());
591 }
LSODA::yh_
std::vector< std::vector< double > > yh_
Definition: lsoda.h:154
LSODA::tn_
double tn_
Definition: lsoda.h:142
LSODA::n
size_t n
Definition: lsoda.h:137
LSODA::illin
size_t illin
Definition: lsoda.h:130
Here is the caller graph for this function:

◆ terminate()

void LSODA::terminate ( int *  istate)

Definition at line 217 of file lsoda.cpp.

218 {
219  if(illin == 5)
220  cerr << "[lsoda] repeated occurrence of illegal input. run aborted.. "
221  "apparent infinite loop."
222  << endl;
223  else {
224  illin++;
225  *istate = -3;
226  }
227 }
LSODA::illin
size_t illin
Definition: lsoda.h:130
Here is the caller graph for this function:

◆ terminate2()

void LSODA::terminate2 ( std::vector< double > &  y,
double *  t 
)

Definition at line 230 of file lsoda.cpp.

231 {
232  for(size_t i = 1; i <= n; i++)
233  y[i] = yh_[1][i];
234  *t = tn_;
235  illin = 0;
236  return;
237 }
LSODA::yh_
std::vector< std::vector< double > > yh_
Definition: lsoda.h:154
LSODA::tn_
double tn_
Definition: lsoda.h:142
LSODA::n
size_t n
Definition: lsoda.h:137
LSODA::illin
size_t illin
Definition: lsoda.h:130
Here is the caller graph for this function:

◆ vmnorm()

double LSODA::vmnorm ( const size_t  n,
const std::vector< double > &  v,
const std::vector< double > &  w 
)

Definition at line 501 of file lsoda.cpp.

502 {
503  double vm = 0.;
504  for(size_t i = 1; i <= n; i++)
505  vm = max(vm, fabs(v[i]) * w[i]);
506  return vm;
507 }
LSODA::n
size_t n
Definition: lsoda.h:137
Here is the caller graph for this function:

Member Data Documentation

◆ acor

std::vector<double> LSODA::acor
private

Definition at line 153 of file lsoda.h.

◆ atol_

std::vector<double> LSODA::atol_
private

Definition at line 162 of file lsoda.h.

◆ ccmax

double LSODA::ccmax
private

Definition at line 141 of file lsoda.h.

◆ cm1

std::array<double, 13> LSODA::cm1
private

Definition at line 124 of file lsoda.h.

◆ cm2

std::array<double, 6> LSODA::cm2
private

Definition at line 125 of file lsoda.h.

◆ conit

double LSODA::conit
private

Definition at line 144 of file lsoda.h.

◆ crate

double LSODA::crate
private

Definition at line 144 of file lsoda.h.

◆ el

std::array<double, 14> LSODA::el
private

Definition at line 123 of file lsoda.h.

◆ el0

double LSODA::el0
private

Definition at line 141 of file lsoda.h.

◆ elco

std::array<std::array<double, 14>, 13> LSODA::elco
private

Definition at line 127 of file lsoda.h.

◆ ETA

double LSODA::ETA = 2.2204460492503131e-16
private

Definition at line 112 of file lsoda.h.

◆ ewt

std::vector<double> LSODA::ewt
private

Definition at line 151 of file lsoda.h.

◆ h_

double LSODA::h_ = .0
private

Definition at line 141 of file lsoda.h.

◆ hmin

double LSODA::hmin
private

Definition at line 142 of file lsoda.h.

◆ hmxi

double LSODA::hmxi
private

Definition at line 142 of file lsoda.h.

◆ hold

double LSODA::hold
private

Definition at line 144 of file lsoda.h.

◆ hu

double LSODA::hu
private

Definition at line 142 of file lsoda.h.

◆ ialth

size_t LSODA::ialth
private

Definition at line 146 of file lsoda.h.

◆ icount

int LSODA::icount
private

Definition at line 149 of file lsoda.h.

◆ ierpj

size_t LSODA::ierpj
private

Definition at line 130 of file lsoda.h.

◆ iersl

size_t LSODA::iersl
private

Definition at line 130 of file lsoda.h.

◆ illin

size_t LSODA::illin
private

Definition at line 130 of file lsoda.h.

◆ imxer

size_t LSODA::imxer
private

Definition at line 114 of file lsoda.h.

◆ init

size_t LSODA::init
private

Definition at line 130 of file lsoda.h.

◆ ipup

size_t LSODA::ipup
private

Definition at line 146 of file lsoda.h.

◆ ipvt

std::vector<int> LSODA::ipvt
private

Definition at line 157 of file lsoda.h.

◆ irflag

int LSODA::irflag
private

Definition at line 149 of file lsoda.h.

◆ iState

int LSODA::iState = 1
private

Definition at line 165 of file lsoda.h.

◆ itol_

int LSODA::itol_ = 2
private

Definition at line 160 of file lsoda.h.

◆ ixpr

size_t LSODA::ixpr = 0
private

Definition at line 134 of file lsoda.h.

◆ jcur

size_t LSODA::jcur
private

Definition at line 130 of file lsoda.h.

◆ jstart

int LSODA::jstart
private

Definition at line 132 of file lsoda.h.

◆ jtyp

size_t LSODA::jtyp
private

Definition at line 134 of file lsoda.h.

◆ kflag

int LSODA::kflag
private

Definition at line 132 of file lsoda.h.

◆ l

size_t LSODA::l
private

Definition at line 130 of file lsoda.h.

◆ lmax

size_t LSODA::lmax
private

Definition at line 146 of file lsoda.h.

◆ maxcor

size_t LSODA::maxcor
private

Definition at line 130 of file lsoda.h.

◆ maxord

size_t LSODA::maxord
private

Definition at line 130 of file lsoda.h.

◆ meth_

size_t LSODA::meth_
private

Definition at line 135 of file lsoda.h.

◆ miter

size_t LSODA::miter
private

Definition at line 130 of file lsoda.h.

◆ ml

size_t LSODA::ml
private

Definition at line 114 of file lsoda.h.

◆ mord

std::array<size_t, 3> LSODA::mord
private

Definition at line 120 of file lsoda.h.

◆ msbp

size_t LSODA::msbp
private

Definition at line 130 of file lsoda.h.

◆ mu

size_t LSODA::mu
private

Definition at line 114 of file lsoda.h.

◆ mused

size_t LSODA::mused
private

Definition at line 134 of file lsoda.h.

◆ mxhnil

size_t LSODA::mxhnil
private

Definition at line 138 of file lsoda.h.

◆ mxncf

size_t LSODA::mxncf
private

Definition at line 130 of file lsoda.h.

◆ mxordn

size_t LSODA::mxordn
private

Definition at line 134 of file lsoda.h.

◆ mxords

size_t LSODA::mxords = 12
private

Definition at line 134 of file lsoda.h.

◆ mxstep

size_t LSODA::mxstep
private

Definition at line 138 of file lsoda.h.

◆ n

size_t LSODA::n
private

Definition at line 137 of file lsoda.h.

◆ nfe

size_t LSODA::nfe
private

Definition at line 137 of file lsoda.h.

◆ nhnil

size_t LSODA::nhnil
private

Definition at line 139 of file lsoda.h.

◆ nje

size_t LSODA::nje
private

Definition at line 137 of file lsoda.h.

◆ nq

size_t LSODA::nq
private

Definition at line 137 of file lsoda.h.

◆ nqu

size_t LSODA::nqu
private

Definition at line 137 of file lsoda.h.

◆ nslast

size_t LSODA::nslast
private

Definition at line 139 of file lsoda.h.

◆ nslp

size_t LSODA::nslp
private

Definition at line 147 of file lsoda.h.

◆ nst

size_t LSODA::nst
private

Definition at line 137 of file lsoda.h.

◆ ntrep

size_t LSODA::ntrep
private

Definition at line 139 of file lsoda.h.

◆ nyh

size_t LSODA::nyh
private

Definition at line 139 of file lsoda.h.

◆ param

void* LSODA::param = nullptr

Definition at line 170 of file lsoda.h.

◆ pdest

double LSODA::pdest
private

Definition at line 148 of file lsoda.h.

◆ pdlast

double LSODA::pdlast
private

Definition at line 148 of file lsoda.h.

◆ pdnorm

double LSODA::pdnorm
private

Definition at line 143 of file lsoda.h.

◆ ratio

double LSODA::ratio
private

Definition at line 148 of file lsoda.h.

◆ rc

double LSODA::rc
private

Definition at line 142 of file lsoda.h.

◆ rmax

double LSODA::rmax
private

Definition at line 144 of file lsoda.h.

◆ rtol_

std::vector<double> LSODA::rtol_
private

Definition at line 161 of file lsoda.h.

◆ savf

std::vector<double> LSODA::savf
private

Definition at line 152 of file lsoda.h.

◆ sm1

std::array<double, 13> LSODA::sm1
private

Definition at line 121 of file lsoda.h.

◆ sqrteta

double LSODA::sqrteta
private

Definition at line 115 of file lsoda.h.

◆ tesco

std::array<std::array<double, 4>, 13> LSODA::tesco
private

Definition at line 128 of file lsoda.h.

◆ tn_

double LSODA::tn_ = 0.0
private

Definition at line 142 of file lsoda.h.

◆ tsw

double LSODA::tsw
private

Definition at line 143 of file lsoda.h.

◆ wm_

std::vector<std::vector<double> > LSODA::wm_
private

Definition at line 155 of file lsoda.h.

◆ yh_

std::vector<std::vector<double> > LSODA::yh_
private

Definition at line 154 of file lsoda.h.

◆ yout

std::vector<double> LSODA::yout
private

Definition at line 167 of file lsoda.h.

The documentation for this class was generated from the following files: