COPASI API  4.16.103
dstoda.cpp
Go to the documentation of this file.
1 // Copyright (C) 2010 - 2014 by Pedro Mendes, Virginia Tech Intellectual
2 // Properties, Inc., University of Heidelberg, and The University
3 // of Manchester.
4 // All rights reserved.
5 
6 // Copyright (C) 2008 - 2009 by Pedro Mendes, Virginia Tech Intellectual
7 // Properties, Inc., EML Research, gGmbH, University of Heidelberg,
8 // and The University of Manchester.
9 // All rights reserved.
10 
11 // Copyright (C) 2006 - 2007 by Pedro Mendes, Virginia Tech Intellectual
12 // Properties, Inc. and EML Research, gGmbH.
13 // All rights reserved.
14 
15 //
16 // This C++ code is based on an f2c conversion of the Fortran
17 // library ODEPACK available at: http://www.netlib.org/odepack/
18 
19 #include <cmath>
20 
21 #include <algorithm>
22 
23 #include "copasi.h"
24 
25 #include "CInternalSolver.h"
26 
27 #define dls001_1 (mdls001_._1)
28 #define dls001_2 (mdls001_._2)
29 #define dls001_3 (mdls001_._3)
30 
31 #define dlsa01_1 (mdlsa01_._1)
32 #define dlsa01_2 (mdlsa01_._2)
33 #define dlsa01_3 (mdlsa01_._3)
34 
35 static C_INT c__1 = 1;
36 static C_INT c__2 = 2;
37 
38 #define pow_dd(__x, __y) pow(*__x, *__y)
39 
40 #include "dcfode.h"
41 #include "dmnorm.h"
42 
43 /* DECK DSTODA */
44 /* Subroutine */
45 C_INT CInternalSolver::dstoda_(C_INT *neq, double *y, double *yh,
46  C_INT *nyh, double *yh1, double *ewt, double *savf,
47  double *acor, double *wm, C_INT *iwm,
48  evalF f, evalJ jac, PJAC * pjac, SLVS * slvs)
49 {
50  /* Initialized data */
51 
52  static const double sm1[12] =
53  {
54  .5, .575, .55, .45, .35, .25, .2, .15, .1, .075, .05,
55  .025
56  };
57 
58  /* System generated locals */
59  C_INT yh_dim1, yh_offset, i__1, i__2;
60  double d__1, d__2, d__3;
61 
62  /* Local variables */
63  C_INT i__, j, m;
64  double r__;
65  C_INT i1, jb;
66  double rh = 0.0, rm, dm1, dm2;
67  C_INT lm1, lm2;
68  double rh1, rh2, del, ddn;
69  C_INT ncf;
70  double pdh = 0.0, dsm, dup, exm1, exm2;
71  C_INT nqm1, nqm2;
72  double dcon, delp;
73  C_INT lm1p1, lm2p1;
74  double exdn, rhdn;
75  C_INT iret;
76  double told, rhsm;
77  C_INT newq;
78  double exsm, rhup, rate, exup, rh1it, alpha;
79  C_INT iredo = 0;
80  double pnorm;
81 
82  /* Parameter adjustments */
83  --neq;
84  --y;
85  yh_dim1 = *nyh;
86  yh_offset = 1 + yh_dim1;
87  yh -= yh_offset;
88  --yh1;
89  --ewt;
90  --savf;
91  --acor;
92  --wm;
93  --iwm;
94 
95  /* Function Body */
96  /* ----------------------------------------------------------------------- */
97  /* DSTODA performs one step of the integration of an initial value */
98  /* problem for a system of ordinary differential equations. */
99  /* Note: DSTODA is independent of the value of the iteration method */
100  /* indicator MITER, when this is .ne. 0, and hence is independent */
101  /* of the type of chord method used, or the Jacobian structure. */
102  /* Communication with DSTODA is done with the following variables: */
103 
104  /* Y = an array of length .ge. N used as the Y argument in */
105  /* all calls to F and JAC. */
106  /* NEQ = integer array containing problem size in NEQ(1), and */
107  /* passed as the NEQ argument in all calls to F and JAC. */
108  /* YH = an NYH by LMAX array containing the dependent variables */
109  /* and their approximate scaled derivatives, where */
110  /* LMAX = MAXORD + 1. YH(i,j+1) contains the approximate */
111  /* j-th derivative of y(i), scaled by H**j/factorial(j) */
112  /* (j = 0,1,...,NQ). On entry for the first step, the first */
113  /* two columns of YH must be set from the initial values. */
114  /* NYH = a constant integer .ge. N, the first dimension of YH. */
115  /* YH1 = a one-dimensional array occupying the same space as YH. */
116  /* EWT = an array of length N containing multiplicative weights */
117  /* for local error measurements. Local errors in y(i) are */
118  /* compared to 1.0/EWT(i) in various error tests. */
119  /* SAVF = an array of working storage, of length N. */
120  /* ACOR = a work array of length N, used for the accumulated */
121  /* corrections. On a successful return, ACOR(i) contains */
122  /* the estimated one-step local error in y(i). */
123  /* WM,IWM = real and integer work arrays associated with matrix */
124  /* operations in chord iteration (MITER .ne. 0). */
125  /* PJAC = name of routine to evaluate and preprocess Jacobian matrix */
126  /* and P = I - H*EL0*Jac, if a chord method is being used. */
127  /* It also returns an estimate of norm(Jac) in PDNORM. */
128  /* SLVS = name of routine to solve linear system in chord iteration. */
129  /* CCMAX = maximum relative change in H*EL0 before PJAC is called. */
130  /* H = the step size to be attempted on the next step. */
131  /* H is altered by the error control algorithm during the */
132  /* problem. H can be either positive or negative, but its */
133  /* sign must remain constant throughout the problem. */
134  /* HMIN = the minimum absolute value of the step size H to be used. */
135  /* HMXI = inverse of the maximum absolute value of H to be used. */
136  /* HMXI = 0.0 is allowed and corresponds to an infinite HMAX. */
137  /* HMIN and HMXI may be changed at any time, but will not */
138  /* take effect until the next change of H is considered. */
139  /* TN = the independent variable. TN is updated on each step taken. */
140  /* JSTART = an integer used for input only, with the following */
141  /* values and meanings: */
142  /* 0 perform the first step. */
143  /* .gt.0 take a new step continuing from the last. */
144  /* -1 take the next step with a new value of H, */
145  /* N, METH, MITER, and/or matrix parameters. */
146  /* -2 take the next step with a new value of H, */
147  /* but with other inputs unchanged. */
148  /* On return, JSTART is set to 1 to facilitate continuation. */
149  /* KFLAG = a completion code with the following meanings: */
150  /* 0 the step was succesful. */
151  /* -1 the requested error could not be achieved. */
152  /* -2 corrector convergence could not be achieved. */
153  /* -3 fatal error in PJAC or SLVS. */
154  /* A return with KFLAG = -1 or -2 means either */
155  /* ABS(H) = HMIN or 10 consecutive failures occurred. */
156  /* On a return with KFLAG negative, the values of TN and */
157  /* the YH array are as of the beginning of the last */
158  /* step, and H is the last step size attempted. */
159  /* MAXORD = the maximum order of integration method to be allowed. */
160  /* MAXCOR = the maximum number of corrector iterations allowed. */
161  /* MSBP = maximum number of steps between PJAC calls (MITER .gt. 0). */
162  /* MXNCF = maximum number of convergence failures allowed. */
163  /* METH = current method. */
164  /* METH = 1 means Adams method (nonstiff) */
165  /* METH = 2 means BDF method (stiff) */
166  /* METH may be reset by DSTODA. */
167  /* MITER = corrector iteration method. */
168  /* MITER = 0 means functional iteration. */
169  /* MITER = JT .gt. 0 means a chord iteration corresponding */
170  /* to Jacobian type JT. (The DLSODA/DLSODAR argument JT is */
171  /* communicated here as JTYP, but is not used in DSTODA */
172  /* except to load MITER following a method switch.) */
173  /* MITER may be reset by DSTODA. */
174  /* N = the number of first-order differential equations. */
175  /* ----------------------------------------------------------------------- */
176  dls001_3.kflag = 0;
177  told = dls001_3.tn;
178  ncf = 0;
179  dls001_3.ierpj = 0;
180  dls001_3.iersl = 0;
181  dls001_3.jcur = 0;
182  dls001_3.icf = 0;
183  delp = 0.;
184 
185  if (dls001_3.jstart > 0)
186  {
187  goto L200;
188  }
189 
190  if (dls001_3.jstart == -1)
191  {
192  goto L100;
193  }
194 
195  if (dls001_3.jstart == -2)
196  {
197  goto L160;
198  }
199 
200  /* ----------------------------------------------------------------------- */
201  /* On the first call, the order is set to 1, and other variables are */
202  /* initialized. RMAX is the maximum ratio by which H can be increased */
203  /* in a single step. It is initially 1.E4 to compensate for the small */
204  /* initial H, but then is normally equal to 10. If a failure */
205  /* occurs (in corrector convergence or error test), RMAX is set at 2 */
206  /* for the next increase. */
207  /* DCFODE is called to get the needed coefficients for both methods. */
208  /* ----------------------------------------------------------------------- */
209  dls001_3.lmax = dls001_3.maxord + 1;
210  dls001_3.nq = 1;
211  dls001_3.l = 2;
212  dls001_3.ialth = 2;
213  dls001_3.rmax = 1e4;
214  dls001_3.rc = 0.;
215  dls001_3.el0 = 1.;
216  dls001_3.crate = .7;
217  dls001_3.hold = dls001_3.h__;
218  dls001_3.nslp = 0;
219  dls001_3.ipup = dls001_3.miter;
220  iret = 3;
221  /* Initialize switching parameters. METH = 1 is assumed initially. ----- */
222  dlsa01_1.icount = 20;
223  dlsa01_1.irflag = 0;
224  dlsa01_1.pdest = 0.;
225  dlsa01_1.pdlast = 0.;
226  dlsa01_1.ratio = 5.;
227  dcfode_(&c__2, dls001_3.elco, dls001_3.tesco);
228 
229  for (i__ = 1; i__ <= 5; ++i__)
230  {
231  /* L10: */
232  dlsa01_1.cm2[i__ - 1] = dls001_3.tesco[i__ * 3 - 2] * dls001_3.elco[
233  i__ + 1 + i__ * 13 - 14];
234  }
235 
236  dcfode_(&c__1, dls001_3.elco, dls001_3.tesco);
237 
238  for (i__ = 1; i__ <= 12; ++i__)
239  {
240  /* L20: */
241  dlsa01_1.cm1[i__ - 1] = dls001_3.tesco[i__ * 3 - 2] * dls001_3.elco[
242  i__ + 1 + i__ * 13 - 14];
243  }
244 
245  goto L150;
246  /* ----------------------------------------------------------------------- */
247  /* The following block handles preliminaries needed when JSTART = -1. */
248  /* IPUP is set to MITER to force a matrix update. */
249  /* If an order increase is about to be considered (IALTH = 1), */
250  /* IALTH is reset to 2 to postpone consideration one more step. */
251  /* If the caller has changed METH, DCFODE is called to reset */
252  /* the coefficients of the method. */
253  /* If H is to be changed, YH must be rescaled. */
254  /* If H or METH is being changed, IALTH is reset to L = NQ + 1 */
255  /* to prevent further changes in H for that many steps. */
256  /* ----------------------------------------------------------------------- */
257 L100:
258  dls001_3.ipup = dls001_3.miter;
259  dls001_3.lmax = dls001_3.maxord + 1;
260 
261  if (dls001_3.ialth == 1)
262  {
263  dls001_3.ialth = 2;
264  }
265 
266  if (dls001_3.meth == dlsa01_1.mused)
267  {
268  goto L160;
269  }
270 
271  dcfode_(&dls001_3.meth, dls001_3.elco, dls001_3.tesco);
272  dls001_3.ialth = dls001_3.l;
273  iret = 1;
274  /* ----------------------------------------------------------------------- */
275  /* The el vector and related constants are reset */
276  /* whenever the order NQ is changed, or at the start of the problem. */
277  /* ----------------------------------------------------------------------- */
278 L150:
279  i__1 = dls001_3.l;
280 
281  for (i__ = 1; i__ <= i__1; ++i__)
282  {
283  /* L155: */
284  dls001_3.el[i__ - 1] = dls001_3.elco[i__ + dls001_3.nq * 13 - 14];
285  }
286 
287  dls001_3.nqnyh = dls001_3.nq * *nyh;
288  dls001_3.rc = dls001_3.rc * dls001_3.el[0] / dls001_3.el0;
289  dls001_3.el0 = dls001_3.el[0];
290  dls001_3.conit = .5 / (dls001_3.nq + 2);
291 
292  switch (iret)
293  {
294  case 1: goto L160;
295  case 2: goto L170;
296  case 3: goto L200;
297  }
298 
299  /* ----------------------------------------------------------------------- */
300  /* If H is being changed, the H ratio RH is checked against */
301  /* RMAX, HMIN, and HMXI, and the YH array rescaled. IALTH is set to */
302  /* L = NQ + 1 to prevent a change of H for that many steps, unless */
303  /* forced by a convergence or error test failure. */
304  /* ----------------------------------------------------------------------- */
305 L160:
306 
307  if (dls001_3.h__ == dls001_3.hold)
308  {
309  goto L200;
310  }
311 
312  rh = dls001_3.h__ / dls001_3.hold;
313  dls001_3.h__ = dls001_3.hold;
314  iredo = 3;
315  goto L175;
316 L170:
317  /* Computing MAX */
318  d__1 = rh, d__2 = dls001_3.hmin / fabs(dls001_3.h__);
319  rh = std::max(d__1, d__2);
320 L175:
321  rh = std::min(rh, dls001_3.rmax);
322  /* Computing MAX */
323  d__1 = 1., d__2 = fabs(dls001_3.h__) * dls001_3.hmxi * rh;
324  rh /= std::max(d__1, d__2);
325 
326  /* ----------------------------------------------------------------------- */
327 
328  /* If METH = 1, also restrict the new step size by the stability region. */
329 
330  /* If this reduces H, set IRFLAG to 1 so that if there are roundoff */
331 
332  /* problems later, we can assume that is the cause of the trouble. */
333 
334  /* ----------------------------------------------------------------------- */
335  if (dls001_3.meth == 2)
336  {
337  goto L178;
338  }
339 
340  dlsa01_1.irflag = 0;
341  /* Computing MAX */
342  d__1 = fabs(dls001_3.h__) * dlsa01_1.pdlast;
343  pdh = std::max(d__1, 1e-6);
344 
345  if (rh * pdh * 1.00001 < sm1[dls001_3.nq - 1])
346  {
347  goto L178;
348  }
349 
350  rh = sm1[dls001_3.nq - 1] / pdh;
351  dlsa01_1.irflag = 1;
352 L178:
353  r__ = 1.;
354  i__1 = dls001_3.l;
355 
356  for (j = 2; j <= i__1; ++j)
357  {
358  r__ *= rh;
359  i__2 = dls001_3.n;
360 
361  for (i__ = 1; i__ <= i__2; ++i__)
362  {
363  /* L180: */
364  yh[i__ + j * yh_dim1] *= r__;
365  }
366  }
367 
368  dls001_3.h__ *= rh;
369  dls001_3.rc *= rh;
370  dls001_3.ialth = dls001_3.l;
371 
372  if (iredo == 0)
373  {
374  goto L690;
375  }
376 
377  /* ----------------------------------------------------------------------- */
378  /* This section computes the predicted values by effectively */
379  /* multiplying the YH array by the Pascal triangle matrix. */
380  /* RC is the ratio of new to old values of the coefficient H*EL(1). */
381  /* When RC differs from 1 by more than CCMAX, IPUP is set to MITER */
382  /* to force PJAC to be called, if a Jacobian is involved. */
383  /* In any case, PJAC is called at least every MSBP steps. */
384  /* ----------------------------------------------------------------------- */
385 L200:
386 
387  if ((d__1 = dls001_3.rc - 1., fabs(d__1)) > dls001_3.ccmax)
388  {
389  dls001_3.ipup = dls001_3.miter;
390  }
391 
392  if (dls001_3.nst >= dls001_3.nslp + dls001_3.msbp)
393  {
394  dls001_3.ipup = dls001_3.miter;
395  }
396 
397  dls001_3.tn += dls001_3.h__;
398  i1 = dls001_3.nqnyh + 1;
399  i__2 = dls001_3.nq;
400 
401  for (jb = 1; jb <= i__2; ++jb)
402  {
403  i1 -= *nyh;
404  /* DIR$ IVDEP */
405  i__1 = dls001_3.nqnyh;
406 
407  for (i__ = i1; i__ <= i__1; ++i__)
408  {
409  /* L210: */
410  yh1[i__] += yh1[i__ + *nyh];
411  }
412 
413  /* L215: */
414  }
415 
416  pnorm = dmnorm_(&dls001_3.n, &yh1[1], &ewt[1]);
417  /* ----------------------------------------------------------------------- */
418  /* Up to MAXCOR corrector iterations are taken. A convergence test is */
419  /* made on the RMS-norm of each correction, weighted by the error */
420  /* weight vector EWT. The sum of the corrections is accumulated in the */
421  /* vector ACOR(i). The YH array is not altered in the corrector loop. */
422  /* ----------------------------------------------------------------------- */
423 L220:
424  m = 0;
425  rate = 0.;
426  del = 0.;
427  i__2 = dls001_3.n;
428 
429  for (i__ = 1; i__ <= i__2; ++i__)
430  {
431  /* L230: */
432  y[i__] = yh[i__ + yh_dim1];
433  }
434 
435  f(&neq[1], &dls001_3.tn, &y[1], &savf[1]);
436  ++dls001_3.nfe;
437 
438  if (dls001_3.ipup <= 0)
439  {
440  goto L250;
441  }
442 
443  /* ----------------------------------------------------------------------- */
444  /* If indicated, the matrix P = I - H*EL(1)*J is reevaluated and */
445  /* preprocessed before starting the corrector iteration. IPUP is set */
446  /* to 0 as an indicator that this has been done. */
447  /* ----------------------------------------------------------------------- */
448  (*pjac)(&neq[1], &y[1], &yh[yh_offset], nyh, &ewt[1], &acor[1], &savf[1],
449  &wm[1], &iwm[1], f, jac);
450  dls001_3.ipup = 0;
451  dls001_3.rc = 1.;
452  dls001_3.nslp = dls001_3.nst;
453  dls001_3.crate = .7;
454 
455  if (dls001_3.ierpj != 0)
456  {
457  goto L430;
458  }
459 
460 L250:
461  i__2 = dls001_3.n;
462 
463  for (i__ = 1; i__ <= i__2; ++i__)
464  {
465  /* L260: */
466  acor[i__] = 0.;
467  }
468 
469 L270:
470 
471  if (dls001_3.miter != 0)
472  {
473  goto L350;
474  }
475 
476  /* ----------------------------------------------------------------------- */
477  /* In the case of functional iteration, update Y directly from */
478  /* the result of the last function evaluation. */
479  /* ----------------------------------------------------------------------- */
480  i__2 = dls001_3.n;
481 
482  for (i__ = 1; i__ <= i__2; ++i__)
483  {
484  savf[i__] = dls001_3.h__ * savf[i__] - yh[i__ + (yh_dim1 << 1)];
485  /* L290: */
486  y[i__] = savf[i__] - acor[i__];
487  }
488 
489  del = dmnorm_(&dls001_3.n, &y[1], &ewt[1]);
490  i__2 = dls001_3.n;
491 
492  for (i__ = 1; i__ <= i__2; ++i__)
493  {
494  y[i__] = yh[i__ + yh_dim1] + dls001_3.el[0] * savf[i__];
495  /* L300: */
496  acor[i__] = savf[i__];
497  }
498 
499  goto L400;
500  /* ----------------------------------------------------------------------- */
501  /* In the case of the chord method, compute the corrector error, */
502  /* and solve the linear system with that as right-hand side and */
503  /* P as coefficient matrix. */
504  /* ----------------------------------------------------------------------- */
505 L350:
506  i__2 = dls001_3.n;
507 
508  for (i__ = 1; i__ <= i__2; ++i__)
509  {
510  /* L360: */
511  y[i__] = dls001_3.h__ * savf[i__] - (yh[i__ + (yh_dim1 << 1)] + acor[
512  i__]);
513  }
514 
515  (*slvs)(&wm[1], &iwm[1], &y[1], &savf[1]);
516 
517  if (dls001_3.iersl < 0)
518  {
519  goto L430;
520  }
521 
522  if (dls001_3.iersl > 0)
523  {
524  goto L410;
525  }
526 
527  del = dmnorm_(&dls001_3.n, &y[1], &ewt[1]);
528  i__2 = dls001_3.n;
529 
530  for (i__ = 1; i__ <= i__2; ++i__)
531  {
532  acor[i__] += y[i__];
533  /* L380: */
534  y[i__] = yh[i__ + yh_dim1] + dls001_3.el[0] * acor[i__];
535  }
536 
537  /* ----------------------------------------------------------------------- */
538  /* Test for convergence. If M .gt. 0, an estimate of the convergence */
539  /* rate constant is stored in CRATE, and this is used in the test. */
540 
541  /* We first check for a change of iterates that is the size of */
542  /* roundoff error. If this occurs, the iteration has converged, and a */
543  /* new rate estimate is not formed. */
544  /* In all other cases, force at least two iterations to estimate a */
545  /* local Lipschitz constant estimate for Adams methods. */
546  /* On convergence, form PDEST = local maximum Lipschitz constant */
547  /* estimate. PDLAST is the most recent nonzero estimate. */
548  /* ----------------------------------------------------------------------- */
549 L400:
550 
551  if (del <= pnorm * 100. * dls001_3.uround)
552  {
553  goto L450;
554  }
555 
556  if (m == 0 && dls001_3.meth == 1)
557  {
558  goto L405;
559  }
560 
561  if (m == 0)
562  {
563  goto L402;
564  }
565 
566  rm = 1024.;
567 
568  if (del <= delp * 1024.)
569  {
570  rm = del / delp;
571  }
572 
573  rate = std::max(rate, rm);
574  /* Computing MAX */
575  d__1 = dls001_3.crate * .2;
576  dls001_3.crate = std::max(d__1, rm);
577 L402:
578  /* Computing MIN */
579  d__1 = 1., d__2 = dls001_3.crate * 1.5;
580  dcon = del * std::min(d__1, d__2) / (dls001_3.tesco[dls001_3.nq * 3 - 2] *
581  dls001_3.conit);
582 
583  if (dcon > 1.)
584  {
585  goto L405;
586  }
587 
588  /* Computing MAX */
589  d__2 = dlsa01_1.pdest, d__3 = rate / (d__1 = dls001_3.h__ * dls001_3.el[0]
590  , fabs(d__1));
591  dlsa01_1.pdest = std::max(d__2, d__3);
592 
593  if (dlsa01_1.pdest != 0.)
594  {
595  dlsa01_1.pdlast = dlsa01_1.pdest;
596  }
597 
598  goto L450;
599 L405:
600  ++m;
601 
602  if (m == dls001_3.maxcor)
603  {
604  goto L410;
605  }
606 
607  if (m >= 2 && del > delp * 2.)
608  {
609  goto L410;
610  }
611 
612  delp = del;
613  f(&neq[1], &dls001_3.tn, &y[1], &savf[1]);
614  ++dls001_3.nfe;
615  goto L270;
616  /* ----------------------------------------------------------------------- */
617  /* The corrector iteration failed to converge. */
618  /* If MITER .ne. 0 and the Jacobian is out of date, PJAC is called for */
619  /* the next try. Otherwise the YH array is retracted to its values */
620  /* before prediction, and H is reduced, if possible. If H cannot be */
621  /* reduced or MXNCF failures have occurred, exit with KFLAG = -2. */
622  /* ----------------------------------------------------------------------- */
623 L410:
624 
625  if (dls001_3.miter == 0 || dls001_3.jcur == 1)
626  {
627  goto L430;
628  }
629 
630  dls001_3.icf = 1;
631  dls001_3.ipup = dls001_3.miter;
632  goto L220;
633 L430:
634  dls001_3.icf = 2;
635  ++ncf;
636  dls001_3.rmax = 2.;
637  dls001_3.tn = told;
638  i1 = dls001_3.nqnyh + 1;
639  i__2 = dls001_3.nq;
640 
641  for (jb = 1; jb <= i__2; ++jb)
642  {
643  i1 -= *nyh;
644  /* DIR$ IVDEP */
645  i__1 = dls001_3.nqnyh;
646 
647  for (i__ = i1; i__ <= i__1; ++i__)
648  {
649  /* L440: */
650  yh1[i__] -= yh1[i__ + *nyh];
651  }
652 
653  /* L445: */
654  }
655 
656  if (dls001_3.ierpj < 0 || dls001_3.iersl < 0)
657  {
658  goto L680;
659  }
660 
661  if (fabs(dls001_3.h__) <= dls001_3.hmin * 1.00001)
662  {
663  goto L670;
664  }
665 
666  if (ncf == dls001_3.mxncf)
667  {
668  goto L670;
669  }
670 
671  rh = .25;
672  dls001_3.ipup = dls001_3.miter;
673  iredo = 1;
674  goto L170;
675  /* ----------------------------------------------------------------------- */
676  /* The corrector has converged. JCUR is set to 0 */
677  /* to signal that the Jacobian involved may need updating later. */
678  /* The local error test is made and control passes to statement 500 */
679  /* if it fails. */
680  /* ----------------------------------------------------------------------- */
681 L450:
682  dls001_3.jcur = 0;
683 
684  if (m == 0)
685  {
686  dsm = del / dls001_3.tesco[dls001_3.nq * 3 - 2];
687  }
688 
689  if (m > 0)
690  {
691  dsm = dmnorm_(&dls001_3.n, &acor[1], &ewt[1]) / dls001_3.tesco[
692  dls001_3.nq * 3 - 2];
693  }
694 
695  if (dsm > 1.)
696  {
697  goto L500;
698  }
699 
700  /* ----------------------------------------------------------------------- */
701  /* After a successful step, update the YH array. */
702  /* Decrease ICOUNT by 1, and if it is -1, consider switching methods. */
703  /* If a method switch is made, reset various parameters, */
704  /* rescale the YH array, and exit. If there is no switch, */
705  /* consider changing H if IALTH = 1. Otherwise decrease IALTH by 1. */
706  /* If IALTH is then 1 and NQ .lt. MAXORD, then ACOR is saved for */
707  /* use in a possible order increase on the next step. */
708  /* If a change in H is considered, an increase or decrease in order */
709  /* by one is considered also. A change in H is made only if it is by a */
710  /* factor of at least 1.1. If not, IALTH is set to 3 to prevent */
711  /* testing for that many steps. */
712  /* ----------------------------------------------------------------------- */
713  dls001_3.kflag = 0;
714  iredo = 0;
715  ++dls001_3.nst;
716  dls001_3.hu = dls001_3.h__;
717  dls001_3.nqu = dls001_3.nq;
718  dlsa01_1.mused = dls001_3.meth;
719  i__2 = dls001_3.l;
720 
721  for (j = 1; j <= i__2; ++j)
722  {
723  i__1 = dls001_3.n;
724 
725  for (i__ = 1; i__ <= i__1; ++i__)
726  {
727  /* L460: */
728  yh[i__ + j * yh_dim1] += dls001_3.el[j - 1] * acor[i__];
729  }
730  }
731 
732  --dlsa01_1.icount;
733 
734  if (dlsa01_1.icount >= 0)
735  {
736  goto L488;
737  }
738 
739  if (dls001_3.meth == 2)
740  {
741  goto L480;
742  }
743 
744  /* ----------------------------------------------------------------------- */
745  /* We are currently using an Adams method. Consider switching to BDF. */
746  /* If the current order is greater than 5, assume the problem is */
747  /* not stiff, and skip this section. */
748  /* If the Lipschitz constant and error estimate are not polluted */
749  /* by roundoff, go to 470 and perform the usual test. */
750  /* Otherwise, switch to the BDF methods if the last step was */
751  /* restricted to insure stability (irflag = 1), and stay with Adams */
752  /* method if not. When switching to BDF with polluted error estimates, */
753  /* in the absence of other information, double the step size. */
754 
755  /* When the estimates are OK, we make the usual test by computing */
756 
757  /* the step size we could have (ideally) used on this step, */
758 
759  /* with the current (Adams) method, and also that for the BDF. */
760 
761  /* If NQ .gt. MXORDS, we consider changing to order MXORDS on switching. */
762 
763  /* Compare the two step sizes to decide whether to switch. */
764 
765  /* The step size advantage must be at least RATIO = 5 to switch. */
766 
767  /* ----------------------------------------------------------------------- */
768  if (dls001_3.nq > 5)
769  {
770  goto L488;
771  }
772 
773  if (dsm > pnorm * 100. * dls001_3.uround && dlsa01_1.pdest != 0.)
774  {
775  goto L470;
776  }
777 
778  if (dlsa01_1.irflag == 0)
779  {
780  goto L488;
781  }
782 
783  rh2 = 2.;
784  nqm2 = std::min(dls001_3.nq, dlsa01_1.mxords);
785  goto L478;
786 L470:
787  exsm = 1. / dls001_3.l;
788  rh1 = 1. / (pow_dd(&dsm, &exsm) * 1.2 + 1.2e-6);
789  rh1it = rh1 * 2.;
790  pdh = dlsa01_1.pdlast * fabs(dls001_3.h__);
791 
792  if (pdh * rh1 > 1e-5)
793  {
794  rh1it = sm1[dls001_3.nq - 1] / pdh;
795  }
796 
797  rh1 = std::min(rh1, rh1it);
798 
799  if (dls001_3.nq <= dlsa01_1.mxords)
800  {
801  goto L474;
802  }
803 
804  nqm2 = dlsa01_1.mxords;
805  lm2 = dlsa01_1.mxords + 1;
806  exm2 = 1. / lm2;
807  lm2p1 = lm2 + 1;
808  dm2 = dmnorm_(&dls001_3.n, &yh[lm2p1 * yh_dim1 + 1], &ewt[1]) /
809  dlsa01_1.cm2[dlsa01_1.mxords - 1];
810  rh2 = 1. / (pow_dd(&dm2, &exm2) * 1.2 + 1.2e-6);
811  goto L476;
812 L474:
813  dm2 = dsm * (dlsa01_1.cm1[dls001_3.nq - 1] / dlsa01_1.cm2[dls001_3.nq - 1]
814  );
815  rh2 = 1. / (pow_dd(&dm2, &exsm) * 1.2 + 1.2e-6);
816  nqm2 = dls001_3.nq;
817 L476:
818 
819  if (rh2 < dlsa01_1.ratio * rh1)
820  {
821  goto L488;
822  }
823 
824  /* THE SWITCH TEST PASSED. RESET RELEVANT QUANTITIES FOR BDF. ---------- */
825 L478:
826  rh = rh2;
827  dlsa01_1.icount = 20;
828  dls001_3.meth = 2;
829  dls001_3.miter = dlsa01_1.jtyp;
830  dlsa01_1.pdlast = 0.;
831  dls001_3.nq = nqm2;
832  dls001_3.l = dls001_3.nq + 1;
833  goto L170;
834  /* ----------------------------------------------------------------------- */
835  /* We are currently using a BDF method. Consider switching to Adams. */
836  /* Compute the step size we could have (ideally) used on this step, */
837  /* with the current (BDF) method, and also that for the Adams. */
838  /* If NQ .gt. MXORDN, we consider changing to order MXORDN on switching. */
839  /* Compare the two step sizes to decide whether to switch. */
840  /* The step size advantage must be at least 5/RATIO = 1 to switch. */
841  /* If the step size for Adams would be so small as to cause */
842  /* roundoff pollution, we stay with BDF. */
843  /* ----------------------------------------------------------------------- */
844 L480:
845  exsm = 1. / dls001_3.l;
846 
847  if (dlsa01_1.mxordn >= dls001_3.nq)
848  {
849  goto L484;
850  }
851 
852  nqm1 = dlsa01_1.mxordn;
853  lm1 = dlsa01_1.mxordn + 1;
854  exm1 = 1. / lm1;
855  lm1p1 = lm1 + 1;
856  dm1 = dmnorm_(&dls001_3.n, &yh[lm1p1 * yh_dim1 + 1], &ewt[1]) /
857  dlsa01_1.cm1[dlsa01_1.mxordn - 1];
858  rh1 = 1. / (pow_dd(&dm1, &exm1) * 1.2 + 1.2e-6);
859  goto L486;
860 L484:
861  dm1 = dsm * (dlsa01_1.cm2[dls001_3.nq - 1] / dlsa01_1.cm1[dls001_3.nq - 1]
862  );
863  rh1 = 1. / (pow_dd(&dm1, &exsm) * 1.2 + 1.2e-6);
864  nqm1 = dls001_3.nq;
865  exm1 = exsm;
866 L486:
867  rh1it = rh1 * 2.;
868  pdh = dlsa01_1.pdnorm * fabs(dls001_3.h__);
869 
870  if (pdh * rh1 > 1e-5)
871  {
872  rh1it = sm1[nqm1 - 1] / pdh;
873  }
874 
875  rh1 = std::min(rh1, rh1it);
876  rh2 = 1. / (pow_dd(&dsm, &exsm) * 1.2 + 1.2e-6);
877 
878  if (rh1 * dlsa01_1.ratio < rh2 * 5.)
879  {
880  goto L488;
881  }
882 
883  alpha = std::max(.001, rh1);
884  dm1 = pow_dd(&alpha, &exm1) * dm1;
885 
886  if (dm1 <= dls001_3.uround * 1e3 * pnorm)
887  {
888  goto L488;
889  }
890 
891  /* The switch test passed. Reset relevant quantities for Adams. -------- */
892  rh = rh1;
893  dlsa01_1.icount = 20;
894  dls001_3.meth = 1;
895  dls001_3.miter = 0;
896  dlsa01_1.pdlast = 0.;
897  dls001_3.nq = nqm1;
898  dls001_3.l = dls001_3.nq + 1;
899  goto L170;
900 
901  /* No method switch is being made. Do the usual step/order selection. -- */
902 L488:
903  --dls001_3.ialth;
904 
905  if (dls001_3.ialth == 0)
906  {
907  goto L520;
908  }
909 
910  if (dls001_3.ialth > 1)
911  {
912  goto L700;
913  }
914 
915  if (dls001_3.l == dls001_3.lmax)
916  {
917  goto L700;
918  }
919 
920  i__1 = dls001_3.n;
921 
922  for (i__ = 1; i__ <= i__1; ++i__)
923  {
924  /* L490: */
925  yh[i__ + dls001_3.lmax * yh_dim1] = acor[i__];
926  }
927 
928  goto L700;
929  /* ----------------------------------------------------------------------- */
930  /* The error test failed. KFLAG keeps track of multiple failures. */
931  /* Restore TN and the YH array to their previous values, and prepare */
932  /* to try the step again. Compute the optimum step size for this or */
933  /* one lower order. After 2 or more failures, H is forced to decrease */
934  /* by a factor of 0.2 or less. */
935  /* ----------------------------------------------------------------------- */
936 L500:
937  --dls001_3.kflag;
938  dls001_3.tn = told;
939  i1 = dls001_3.nqnyh + 1;
940  i__1 = dls001_3.nq;
941 
942  for (jb = 1; jb <= i__1; ++jb)
943  {
944  i1 -= *nyh;
945  /* DIR$ IVDEP */
946  i__2 = dls001_3.nqnyh;
947 
948  for (i__ = i1; i__ <= i__2; ++i__)
949  {
950  /* L510: */
951  yh1[i__] -= yh1[i__ + *nyh];
952  }
953 
954  /* L515: */
955  }
956 
957  dls001_3.rmax = 2.;
958 
959  if (fabs(dls001_3.h__) <= dls001_3.hmin * 1.00001)
960  {
961  goto L660;
962  }
963 
964  if (dls001_3.kflag <= -3)
965  {
966  goto L640;
967  }
968 
969  iredo = 2;
970  rhup = 0.;
971  goto L540;
972  /* ----------------------------------------------------------------------- */
973  /* Regardless of the success or failure of the step, factors */
974  /* RHDN, RHSM, and RHUP are computed, by which H could be multiplied */
975  /* at order NQ - 1, order NQ, or order NQ + 1, respectively. */
976  /* In the case of failure, RHUP = 0.0 to avoid an order increase. */
977  /* The largest of these is determined and the new order chosen */
978  /* accordingly. If the order is to be increased, we compute one */
979  /* additional scaled derivative. */
980  /* ----------------------------------------------------------------------- */
981 L520:
982  rhup = 0.;
983 
984  if (dls001_3.l == dls001_3.lmax)
985  {
986  goto L540;
987  }
988 
989  i__1 = dls001_3.n;
990 
991  for (i__ = 1; i__ <= i__1; ++i__)
992  {
993  /* L530: */
994  savf[i__] = acor[i__] - yh[i__ + dls001_3.lmax * yh_dim1];
995  }
996 
997  dup = dmnorm_(&dls001_3.n, &savf[1], &ewt[1]) / dls001_3.tesco[
998  dls001_3.nq * 3 - 1];
999  exup = 1. / (dls001_3.l + 1);
1000  rhup = 1. / (pow_dd(&dup, &exup) * 1.4 + 1.4e-6);
1001 L540:
1002  exsm = 1. / dls001_3.l;
1003  rhsm = 1. / (pow_dd(&dsm, &exsm) * 1.2 + 1.2e-6);
1004  rhdn = 0.;
1005 
1006  if (dls001_3.nq == 1)
1007  {
1008  goto L550;
1009  }
1010 
1011  ddn = dmnorm_(&dls001_3.n, &yh[dls001_3.l * yh_dim1 + 1], &ewt[1]) /
1012  dls001_3.tesco[dls001_3.nq * 3 - 3];
1013  exdn = 1. / dls001_3.nq;
1014  rhdn = 1. / (pow_dd(&ddn, &exdn) * 1.3 + 1.3e-6);
1015  /* If METH = 1, limit RH according to the stability region also. -------- */
1016 L550:
1017 
1018  if (dls001_3.meth == 2)
1019  {
1020  goto L560;
1021  }
1022 
1023  /* Computing MAX */
1024  d__1 = fabs(dls001_3.h__) * dlsa01_1.pdlast;
1025  pdh = std::max(d__1, 1e-6);
1026 
1027  if (dls001_3.l < dls001_3.lmax)
1028  {
1029  /* Computing MIN */
1030  d__1 = rhup, d__2 = sm1[dls001_3.l - 1] / pdh;
1031  rhup = std::min(d__1, d__2);
1032  }
1033 
1034  /* Computing MIN */
1035  d__1 = rhsm, d__2 = sm1[dls001_3.nq - 1] / pdh;
1036  rhsm = std::min(d__1, d__2);
1037 
1038  if (dls001_3.nq > 1)
1039  {
1040  /* Computing MIN */
1041  d__1 = rhdn, d__2 = sm1[dls001_3.nq - 2] / pdh;
1042  rhdn = std::min(d__1, d__2);
1043  }
1044 
1045  dlsa01_1.pdest = 0.;
1046 L560:
1047 
1048  if (rhsm >= rhup)
1049  {
1050  goto L570;
1051  }
1052 
1053  if (rhup > rhdn)
1054  {
1055  goto L590;
1056  }
1057 
1058  goto L580;
1059 L570:
1060 
1061  if (rhsm < rhdn)
1062  {
1063  goto L580;
1064  }
1065 
1066  newq = dls001_3.nq;
1067  rh = rhsm;
1068  goto L620;
1069 L580:
1070  newq = dls001_3.nq - 1;
1071  rh = rhdn;
1072 
1073  if (dls001_3.kflag < 0 && rh > 1.)
1074  {
1075  rh = 1.;
1076  }
1077 
1078  goto L620;
1079 L590:
1080  newq = dls001_3.l;
1081  rh = rhup;
1082 
1083  if (rh < 1.1)
1084  {
1085  goto L610;
1086  }
1087 
1088  r__ = dls001_3.el[dls001_3.l - 1] / dls001_3.l;
1089  i__1 = dls001_3.n;
1090 
1091  for (i__ = 1; i__ <= i__1; ++i__)
1092  {
1093  /* L600: */
1094  yh[i__ + (newq + 1) * yh_dim1] = acor[i__] * r__;
1095  }
1096 
1097  goto L630;
1098 L610:
1099  dls001_3.ialth = 3;
1100  goto L700;
1101  /* If METH = 1 and H is restricted by stability, bypass 10 percent test. */
1102 L620:
1103 
1104  if (dls001_3.meth == 2)
1105  {
1106  goto L622;
1107  }
1108 
1109  if (rh * pdh * 1.00001 >= sm1[newq - 1])
1110  {
1111  goto L625;
1112  }
1113 
1114 L622:
1115 
1116  if (dls001_3.kflag == 0 && rh < 1.1)
1117  {
1118  goto L610;
1119  }
1120 
1121 L625:
1122 
1123  if (dls001_3.kflag <= -2)
1124  {
1125  rh = std::min(rh, .2);
1126  }
1127 
1128  /* ----------------------------------------------------------------------- */
1129 
1130  /* If there is a change of order, reset NQ, L, and the coefficients. */
1131 
1132  /* In any case H is reset according to RH and the YH array is rescaled. */
1133 
1134  /* Then exit from 690 if the step was OK, or redo the step otherwise. */
1135 
1136  /* ----------------------------------------------------------------------- */
1137  if (newq == dls001_3.nq)
1138  {
1139  goto L170;
1140  }
1141 
1142 L630:
1143  dls001_3.nq = newq;
1144  dls001_3.l = dls001_3.nq + 1;
1145  iret = 2;
1146  goto L150;
1147  /* ----------------------------------------------------------------------- */
1148  /* Control reaches this section if 3 or more failures have occured. */
1149  /* If 10 failures have occurred, exit with KFLAG = -1. */
1150  /* It is assumed that the derivatives that have accumulated in the */
1151  /* YH array have errors of the wrong order. Hence the first */
1152  /* derivative is recomputed, and the order is set to 1. Then */
1153  /* H is reduced by a factor of 10, and the step is retried, */
1154  /* until it succeeds or H reaches HMIN. */
1155  /* ----------------------------------------------------------------------- */
1156 L640:
1157 
1158  if (dls001_3.kflag == -10)
1159  {
1160  goto L660;
1161  }
1162 
1163  rh = .1;
1164  /* Computing MAX */
1165  d__1 = dls001_3.hmin / fabs(dls001_3.h__);
1166  rh = std::max(d__1, rh);
1167  dls001_3.h__ *= rh;
1168  i__1 = dls001_3.n;
1169 
1170  for (i__ = 1; i__ <= i__1; ++i__)
1171  {
1172  /* L645: */
1173  y[i__] = yh[i__ + yh_dim1];
1174  }
1175 
1176  f(&neq[1], &dls001_3.tn, &y[1], &savf[1]);
1177  ++dls001_3.nfe;
1178  i__1 = dls001_3.n;
1179 
1180  for (i__ = 1; i__ <= i__1; ++i__)
1181  {
1182  /* L650: */
1183  yh[i__ + (yh_dim1 << 1)] = dls001_3.h__ * savf[i__];
1184  }
1185 
1186  dls001_3.ipup = dls001_3.miter;
1187  dls001_3.ialth = 5;
1188 
1189  if (dls001_3.nq == 1)
1190  {
1191  goto L200;
1192  }
1193 
1194  dls001_3.nq = 1;
1195  dls001_3.l = 2;
1196  iret = 3;
1197  goto L150;
1198  /* ----------------------------------------------------------------------- */
1199  /* All returns are made through this section. H is saved in HOLD */
1200  /* to allow the caller to change H on the next step. */
1201  /* ----------------------------------------------------------------------- */
1202 L660:
1203  dls001_3.kflag = -1;
1204  goto L720;
1205 L670:
1206  dls001_3.kflag = -2;
1207  goto L720;
1208 L680:
1209  dls001_3.kflag = -3;
1210  goto L720;
1211 L690:
1212  dls001_3.rmax = 10.;
1213 L700:
1214  r__ = 1. / dls001_3.tesco[dls001_3.nqu * 3 - 2];
1215  i__1 = dls001_3.n;
1216 
1217  for (i__ = 1; i__ <= i__1; ++i__)
1218  {
1219  /* L710: */
1220  acor[i__] *= r__;
1221  }
1222 
1223 L720:
1224  dls001_3.hold = dls001_3.h__;
1225  dls001_3.jstart = 1;
1226  return 0;
1227  /* ----------------------- End of Subroutine DSTODA ---------------------- */
1228 } /* dstoda_ */
C_INT dstoda_(C_INT *neq, double *y, double *yh, C_INT *nyh, double *yh1, double *ewt, double *savf, double *acor, double *wm, C_INT *iwm, evalF f, evalJ jac, PJAC *pjac, SLVS *slvs)
Definition: dstoda.cpp:45
#define C_INT
Definition: copasi.h:115
Definition: common.h:190
#define pow_dd(__x, __y)
Definition: dstoda.cpp:38
void(* evalF)(const C_INT *, const double *, const double *, double *)
Definition: common.h:29
static C_INT c__2
Definition: dstoda.cpp:36
void(* evalJ)(const C_INT *, const double *, const double *, const C_INT *, const C_INT *, double *, const C_INT *)
Definition: common.h:30
C_INT dcfode_(C_INT *meth, double *elco, double *tesco)
Definition: dcfode.cpp:22
static C_INT c__1
Definition: dstoda.cpp:35
#define dls001_3
Definition: dstoda.cpp:29
#define dlsa01_1
Definition: dstoda.cpp:31
double dmnorm_(C_INT *n, double *v, double *w)
Definition: dmnorm.cpp:35
#define min(a, b)
Definition: f2c.h:175
Definition: common.h:130
#define max(a, b)
Definition: f2c.h:176