COPASI API  4.16.103
drcheck.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 #include <algorithm>
21 
22 #include "copasi.h"
23 
24 #include "CInternalSolver.h"
25 
26 #include "lapack/blaswrap.h"
27 #include "lapack/lapackwrap.h"
28 
29 #define dls001_1 (mdls001_._1)
30 #define dls001_2 (mdls001_._2)
31 #define dls001_3 (mdls001_._3)
32 
33 #define dlsr01_1 (mdlsr01_._1)
34 #define dlsr01_2 (mdlsr01_._2)
35 #define dlsr01_3 (mdlsr01_._3)
36 
37 double d_sign(const double & a, const double & b);
38 
39 static double c_b3 = 1.0;
40 
41 static const C_INT c__0 = 0;
42 static C_INT c__1 = 1;
43 
44 /* DECK DRCHEK */
45 /* Subroutine */
46 C_INT CInternalSolver::drchek_(const C_INT *job, evalG g, C_INT *neq, double *
47  y, double *yh, C_INT *nyh, double *g0, double *g1,
48  double *gx, C_INT *jroot, C_INT *irt)
49 {
50  /* System generated locals */
51  C_INT yh_dim1, yh_offset, i__1;
52  double d__1;
53 
54  /* Local variables */
55  C_INT i__;
56  double x, t1, temp1, temp2;
57  C_INT iflag, jflag;
58  double hming;
59  bool zroot;
60 
61  /* ----------------------------------------------------------------------- */
62  /* This routine checks for the presence of a root in the vicinity of */
63  /* the current T, in a manner depending on the input flag JOB. It calls */
64  /* Subroutine DROOTS to locate the root as precisely as possible. */
65 
66  /* In addition to variables described previously, DRCHEK */
67  /* uses the following for communication: */
68  /* JOB = integer flag indicating type of call: */
69  /* JOB = 1 means the problem is being initialized, and DRCHEK */
70  /* is to look for a root at or very near the initial T. */
71  /* JOB = 2 means a continuation call to the solver was just */
72  /* made, and DRCHEK is to check for a root in the */
73  /* relevant part of the step last taken. */
74  /* JOB = 3 means a successful step was just taken, and DRCHEK */
75  /* is to look for a root in the interval of the step. */
76  /* G0 = array of length NG, containing the value of g at T = T0. */
77  /* G0 is input for JOB .ge. 2, and output in all cases. */
78  /* G1,GX = arrays of length NG for work space. */
79  /* IRT = completion flag: */
80  /* IRT = 0 means no root was found. */
81  /* IRT = -1 means JOB = 1 and a root was found too near to T. */
82  /* IRT = 1 means a legitimate root was found (JOB = 2 or 3). */
83  /* On return, T0 is the root location, and Y is the */
84  /* corresponding solution vector. */
85  /* T0 = value of T at one endpoint of interval of interest. Only */
86  /* roots beyond T0 in the direction of integration are sought. */
87  /* T0 is input if JOB .ge. 2, and output in all cases. */
88  /* T0 is updated by DRCHEK, whether a root is found or not. */
89  /* TLAST = last value of T returned by the solver (input only). */
90  /* TOUTC = copy of TOUT (input only). */
91  /* IRFND = input flag showing whether the last step taken had a root. */
92  /* IRFND = 1 if it did, = 0 if not. */
93  /* ITASKC = copy of ITASK (input only). */
94  /* NGC = copy of NG (input only). */
95  /* ----------------------------------------------------------------------- */
96  /* Parameter adjustments */
97  --neq;
98  --y;
99  yh_dim1 = *nyh;
100  yh_offset = 1 + yh_dim1;
101  yh -= yh_offset;
102  --g0;
103  --g1;
104  --gx;
105  --jroot;
106 
107  /* Function Body */
108  *irt = 0;
109  i__1 = dlsr01_1.ngc;
110 
111  for (i__ = 1; i__ <= i__1; ++i__)
112  {
113  /* L10: */
114  jroot[i__] = 0;
115  }
116 
117  hming = (fabs(dls001_1.tn) + fabs(dls001_1.h__)) * dls001_1.uround * 100.;
118 
119  switch (*job)
120  {
121  case 1: goto L100;
122  case 2: goto L200;
123  case 3: goto L300;
124  }
125 
126  /* Evaluate g at initial T, and check for zero values. ------------------ */
127 L100:
128  dlsr01_1.t0 = dls001_1.tn;
129  (*g)(&neq[1], &dlsr01_1.t0, &y[1], &dlsr01_1.ngc, &g0[1]);
130  dlsr01_1.nge = 1;
131  zroot = false;
132  i__1 = dlsr01_1.ngc;
133 
134  for (i__ = 1; i__ <= i__1; ++i__)
135  {
136  /* L110: */
137  if ((d__1 = g0[i__], fabs(d__1)) <= 0.)
138  {
139  zroot = true;
140  }
141  }
142 
143  if (! zroot)
144  {
145  goto L190;
146  }
147 
148  /* g has a zero at T. Look at g at T + (small increment). -------------- */
149  /* Computing MAX */
150  d__1 = hming / fabs(dls001_1.h__);
151  temp2 = std::max(d__1, .1);
152  temp1 = temp2 * dls001_1.h__;
153  dlsr01_1.t0 += temp1;
154  i__1 = dls001_1.n;
155 
156  for (i__ = 1; i__ <= i__1; ++i__)
157  {
158  /* L120: */
159  y[i__] += temp2 * yh[i__ + (yh_dim1 << 1)];
160  }
161 
162  (*g)(&neq[1], &dlsr01_1.t0, &y[1], &dlsr01_1.ngc, &g0[1]);
163  ++dlsr01_1.nge;
164  zroot = false;
165  i__1 = dlsr01_1.ngc;
166 
167  for (i__ = 1; i__ <= i__1; ++i__)
168  {
169  /* L130: */
170  if ((d__1 = g0[i__], fabs(d__1)) <= 0.)
171  {
172  zroot = true;
173  }
174  }
175 
176  if (! zroot)
177  {
178  goto L190;
179  }
180 
181  /* g has a zero at T and also close to T. Take error return. ----------- */
182  *irt = -1;
183  return 0;
184 
185 L190:
186  return 0;
187 
188 L200:
189 
190  if (dlsr01_1.irfnd == 0)
191  {
192  goto L260;
193  }
194 
195  /* If a root was found on the previous step, evaluate G0 = g(T0). ------- */
196  dintdy_(&dlsr01_1.t0, &c__0, &yh[yh_offset], nyh, &y[1], &iflag);
197  (*g)(&neq[1], &dlsr01_1.t0, &y[1], &dlsr01_1.ngc, &g0[1]);
198  ++dlsr01_1.nge;
199  zroot = false;
200  i__1 = dlsr01_1.ngc;
201 
202  for (i__ = 1; i__ <= i__1; ++i__)
203  {
204  /* L210: */
205  if ((d__1 = g0[i__], fabs(d__1)) <= 0.)
206  {
207  zroot = true;
208  }
209  }
210 
211  if (! zroot)
212  {
213  goto L260;
214  }
215 
216  /* g has a zero at T0. Look at g at T + (small increment). ------------- */
217  temp1 = d_sign(hming, dls001_1.h__);
218  dlsr01_1.t0 += temp1;
219 
220  if ((dlsr01_1.t0 - dls001_1.tn) * dls001_1.h__ < 0.)
221  {
222  goto L230;
223  }
224 
225  temp2 = temp1 / dls001_1.h__;
226  i__1 = dls001_1.n;
227 
228  for (i__ = 1; i__ <= i__1; ++i__)
229  {
230  /* L220: */
231  y[i__] += temp2 * yh[i__ + (yh_dim1 << 1)];
232  }
233 
234  goto L240;
235 L230:
236  dintdy_(&dlsr01_1.t0, &c__0, &yh[yh_offset], nyh, &y[1], &iflag);
237 L240:
238  (*g)(&neq[1], &dlsr01_1.t0, &y[1], &dlsr01_1.ngc, &g0[1]);
239  ++dlsr01_1.nge;
240  zroot = false;
241  i__1 = dlsr01_1.ngc;
242 
243  for (i__ = 1; i__ <= i__1; ++i__)
244  {
245  if ((d__1 = g0[i__], fabs(d__1)) > 0.)
246  {
247  goto L250;
248  }
249 
250  jroot[i__] = 1;
251  zroot = true;
252 L250:
253  ;
254  }
255 
256  if (! zroot)
257  {
258  goto L260;
259  }
260 
261  /* g has a zero at T0 and also close to T0. Return root. --------------- */
262  *irt = 1;
263  return 0;
264  /* G0 has no zero components. Proceed to check relevant interval. ------ */
265 L260:
266 
267  if (dls001_1.tn == dlsr01_1.tlast)
268  {
269  goto L390;
270  }
271 
272 L300:
273 
274  /* Set T1 to TN or TOUTC, whichever comes first, and get g at T1. ------- */
275  if (dlsr01_1.itaskc == 2 || dlsr01_1.itaskc == 3 || dlsr01_1.itaskc == 5)
276  {
277  goto L310;
278  }
279 
280  if ((dlsr01_1.toutc - dls001_1.tn) * dls001_1.h__ >= 0.)
281  {
282  goto L310;
283  }
284 
285  t1 = dlsr01_1.toutc;
286 
287  if ((t1 - dlsr01_1.t0) * dls001_1.h__ <= 0.)
288  {
289  goto L390;
290  }
291 
292  dintdy_(&t1, &c__0, &yh[yh_offset], nyh, &y[1], &iflag);
293  goto L330;
294 L310:
295  t1 = dls001_1.tn;
296  i__1 = dls001_1.n;
297 
298  for (i__ = 1; i__ <= i__1; ++i__)
299  {
300  /* L320: */
301  y[i__] = yh[i__ + yh_dim1];
302  }
303 
304 L330:
305  (*g)(&neq[1], &t1, &y[1], &dlsr01_1.ngc, &g1[1]);
306  ++dlsr01_1.nge;
307  /* Call DROOTS to search for root in interval from T0 to T1. ------------ */
308  jflag = 0;
309 L350:
310  droots_(&dlsr01_1.ngc, &hming, &jflag, &dlsr01_1.t0, &t1, &g0[1], &g1[1],
311  &gx[1], &x, &jroot[1]);
312 
313  if (jflag > 1)
314  {
315  goto L360;
316  }
317 
318  dintdy_(&x, &c__0, &yh[yh_offset], nyh, &y[1], &iflag);
319  (*g)(&neq[1], &x, &y[1], &dlsr01_1.ngc, &gx[1]);
320  ++dlsr01_1.nge;
321  goto L350;
322 L360:
323  dlsr01_1.t0 = x;
324  dcopy_(&dlsr01_1.ngc, &gx[1], &c__1, &g0[1], &c__1);
325 
326  if (jflag == 4)
327  {
328  goto L390;
329  }
330 
331  /* Found a root. Interpolate to X and return. -------------------------- */
332  dintdy_(&x, &c__0, &yh[yh_offset], nyh, &y[1], &iflag);
333  *irt = 1;
334  return 0;
335 
336 L390:
337  return 0;
338  /* ----------------------- End of Subroutine DRCHEK ---------------------- */
339 } /* drchek_ */
340 
341 /* DECK DROOTS */
342 /* Subroutine */
343 C_INT CInternalSolver::droots_(C_INT *ng, double *hmin, C_INT *jflag,
344  double *x0, double *x1, double *g0, double *g1,
345  double *gx, double *x, C_INT *jroot)
346 {
347  /* Initialized data */
348 
349  static double zero = 0.;
350  static double half = .5;
351  static double tenth = .1;
352  static double five = 5.;
353 
354  /* System generated locals */
355  C_INT i__1;
356  double d__1;
357 
358  /* Local variables */
359  C_INT i__;
360  double t2, tmax;
361  bool xroot, zroot, sgnchg;
362  C_INT imxold, nxlast;
363  double fracsub, fracint;
364 
365  /* ----------------------------------------------------------------------- */
366  /* This subroutine finds the leftmost root of a set of arbitrary */
367  /* functions gi(x) (i = 1,...,NG) in an interval (X0,X1). Only roots */
368  /* of odd multiplicity (i.e. changes of sign of the gi) are found. */
369  /* Here the sign of X1 - X0 is arbitrary, but is constant for a given */
370  /* problem, and -leftmost- means nearest to X0. */
371  /* The values of the vector-valued function g(x) = (gi, i=1...NG) */
372  /* are communicated through the call sequence of DROOTS. */
373  /* The method used is the Illinois algorithm. */
374 
375  /* Reference: */
376  /* Kathie L. Hiebert and Lawrence F. Shampine, Implicitly Defined */
377  /* Output Points for Solutions of ODEs, Sandia Report SAND80-0180, */
378  /* February 1980. */
379 
380  /* Description of parameters. */
381 
382  /* NG = number of functions gi, or the number of components of */
383  /* the vector valued function g(x). Input only. */
384 
385  /* HMIN = resolution parameter in X. Input only. When a root is */
386  /* found, it is located only to within an error of HMIN in X. */
387  /* Typically, HMIN should be set to something on the order of */
388  /* 100 * UROUND * MAX(ABS(X0),ABS(X1)), */
389  /* where UROUND is the unit roundoff of the machine. */
390 
391  /* JFLAG = integer flag for input and output communication. */
392 
393  /* On input, set JFLAG = 0 on the first call for the problem, */
394  /* and leave it unchanged until the problem is completed. */
395  /* (The problem is completed when JFLAG .ge. 2 on return.) */
396 
397  /* On output, JFLAG has the following values and meanings: */
398  /* JFLAG = 1 means DROOTS needs a value of g(x). Set GX = g(X) */
399  /* and call DROOTS again. */
400  /* JFLAG = 2 means a root has been found. The root is */
401  /* at X, and GX contains g(X). (Actually, X is the */
402  /* rightmost approximation to the root on an interval */
403  /* (X0,X1) of size HMIN or less.) */
404  /* JFLAG = 3 means X = X1 is a root, with one or more of the gi */
405  /* being zero at X1 and no sign changes in (X0,X1). */
406  /* GX contains g(X) on output. */
407  /* JFLAG = 4 means no roots (of odd multiplicity) were */
408  /* found in (X0,X1) (no sign changes). */
409 
410  /* X0,X1 = endpoints of the interval where roots are sought. */
411  /* X1 and X0 are input when JFLAG = 0 (first call), and */
412  /* must be left unchanged between calls until the problem is */
413  /* completed. X0 and X1 must be distinct, but X1 - X0 may be */
414  /* of either sign. However, the notion of -left- and -right- */
415  /* will be used to mean nearer to X0 or X1, respectively. */
416  /* When JFLAG .ge. 2 on return, X0 and X1 are output, and */
417  /* are the endpoints of the relevant interval. */
418 
419  /* G0,G1 = arrays of length NG containing the vectors g(X0) and g(X1), */
420  /* respectively. When JFLAG = 0, G0 and G1 are input and */
421  /* none of the G0(i) should be zero. */
422  /* When JFLAG .ge. 2 on return, G0 and G1 are output. */
423 
424  /* GX = array of length NG containing g(X). GX is input */
425  /* when JFLAG = 1, and output when JFLAG .ge. 2. */
426 
427  /* X = independent variable value. Output only. */
428  /* When JFLAG = 1 on output, X is the point at which g(x) */
429  /* is to be evaluated and loaded into GX. */
430  /* When JFLAG = 2 or 3, X is the root. */
431  /* When JFLAG = 4, X is the right endpoint of the interval, X1. */
432 
433  /* JROOT = integer array of length NG. Output only. */
434  /* When JFLAG = 2 or 3, JROOT indicates which components */
435  /* of g(x) have a root at X. JROOT(i) is 1 if the i-th */
436  /* component has a root, and JROOT(i) = 0 otherwise. */
437  /* ----------------------------------------------------------------------- */
438  /* Parameter adjustments */
439  --jroot;
440  --gx;
441  --g1;
442  --g0;
443 
444  /* Function Body */
445 
446  if (*jflag == 1)
447  {
448  goto L200;
449  }
450 
451  /* JFLAG .ne. 1. Check for change in sign of g or zero at X1. ---------- */
452  dlsr01_2.imax = 0;
453  tmax = zero;
454  zroot = false;
455  i__1 = *ng;
456 
457  for (i__ = 1; i__ <= i__1; ++i__)
458  {
459  if ((d__1 = g1[i__], fabs(d__1)) > zero)
460  {
461  goto L110;
462  }
463 
464  zroot = true;
465  goto L120;
466  /* At this point, G0(i) has been checked and cannot be zero. ------------ */
467 L110:
468 
469  if (d_sign(c_b3, g0[i__]) == d_sign(c_b3, g1[i__]))
470  {
471  goto L120;
472  }
473 
474  t2 = (d__1 = g1[i__] / (g1[i__] - g0[i__]), fabs(d__1));
475 
476  if (t2 <= tmax)
477  {
478  goto L120;
479  }
480 
481  tmax = t2;
482  dlsr01_2.imax = i__;
483 L120:
484  ;
485  }
486 
487  if (dlsr01_2.imax > 0)
488  {
489  goto L130;
490  }
491 
492  sgnchg = false;
493  goto L140;
494 L130:
495  sgnchg = true;
496 L140:
497 
498  if (! sgnchg)
499  {
500  goto L400;
501  }
502 
503  /* There is a sign change. Find the first root in the interval. -------- */
504  xroot = false;
505  nxlast = 0;
506  dlsr01_2.last = 1;
507 
508  /* Repeat until the first root in the interval is found. Loop point. --- */
509 L150:
510 
511  if (xroot)
512  {
513  goto L300;
514  }
515 
516  if (nxlast == dlsr01_2.last)
517  {
518  goto L160;
519  }
520 
521  dlsr01_2.alpha = 1.;
522  goto L180;
523 L160:
524 
525  if (dlsr01_2.last == 0)
526  {
527  goto L170;
528  }
529 
530  dlsr01_2.alpha *= .5;
531  goto L180;
532 L170:
533  dlsr01_2.alpha *= 2.;
534 L180:
535  dlsr01_2.x2 = *x1 - (*x1 - *x0) * g1[dlsr01_2.imax] / (g1[dlsr01_2.imax]
536  - dlsr01_2.alpha * g0[dlsr01_2.imax]);
537 
538  /* If X2 is too close to X0 or X1, adjust it inward, by a fractional ---- */
539 
540  /* distance that is between 0.1 and 0.5. -------------------------------- */
541  if ((d__1 = dlsr01_2.x2 - *x0, fabs(d__1)) < half * *hmin)
542  {
543  fracint = (d__1 = *x1 - *x0, fabs(d__1)) / *hmin;
544  fracsub = tenth;
545 
546  if (fracint <= five)
547  {
548  fracsub = half / fracint;
549  }
550 
551  dlsr01_2.x2 = *x0 + fracsub * (*x1 - *x0);
552  }
553 
554  if ((d__1 = *x1 - dlsr01_2.x2, fabs(d__1)) < half * *hmin)
555  {
556  fracint = (d__1 = *x1 - *x0, fabs(d__1)) / *hmin;
557  fracsub = tenth;
558 
559  if (fracint <= five)
560  {
561  fracsub = half / fracint;
562  }
563 
564  dlsr01_2.x2 = *x1 - fracsub * (*x1 - *x0);
565  }
566 
567  *jflag = 1;
568  *x = dlsr01_2.x2;
569  /* Return to the calling routine to get a value of GX = g(X). ----------- */
570  return 0;
571  /* Check to see in which interval g changes sign. ----------------------- */
572 L200:
573  imxold = dlsr01_2.imax;
574  dlsr01_2.imax = 0;
575  tmax = zero;
576  zroot = false;
577  i__1 = *ng;
578 
579  for (i__ = 1; i__ <= i__1; ++i__)
580  {
581  if ((d__1 = gx[i__], fabs(d__1)) > zero)
582  {
583  goto L210;
584  }
585 
586  zroot = true;
587  goto L220;
588  /* Neither G0(i) nor GX(i) can be zero at this point. ------------------- */
589 L210:
590 
591  if (d_sign(c_b3, g0[i__]) == d_sign(c_b3, gx[i__]))
592  {
593  goto L220;
594  }
595 
596  t2 = (d__1 = gx[i__] / (gx[i__] - g0[i__]), fabs(d__1));
597 
598  if (t2 <= tmax)
599  {
600  goto L220;
601  }
602 
603  tmax = t2;
604  dlsr01_2.imax = i__;
605 L220:
606  ;
607  }
608 
609  if (dlsr01_2.imax > 0)
610  {
611  goto L230;
612  }
613 
614  sgnchg = false;
615  dlsr01_2.imax = imxold;
616  goto L240;
617 L230:
618  sgnchg = true;
619 L240:
620  nxlast = dlsr01_2.last;
621 
622  if (! sgnchg)
623  {
624  goto L250;
625  }
626 
627  /* Sign change between X0 and X2, so replace X1 with X2. ---------------- */
628  *x1 = dlsr01_2.x2;
629  dcopy_(ng, &gx[1], &c__1, &g1[1], &c__1);
630  dlsr01_2.last = 1;
631  xroot = false;
632  goto L270;
633 L250:
634 
635  if (! zroot)
636  {
637  goto L260;
638  }
639 
640  /* Zero value at X2 and no sign change in (X0,X2), so X2 is a root. ----- */
641  *x1 = dlsr01_2.x2;
642  dcopy_(ng, &gx[1], &c__1, &g1[1], &c__1);
643  xroot = true;
644  goto L270;
645  /* No sign change between X0 and X2. Replace X0 with X2. --------------- */
646 L260:
647  dcopy_(ng, &gx[1], &c__1, &g0[1], &c__1);
648  *x0 = dlsr01_2.x2;
649  dlsr01_2.last = 0;
650  xroot = false;
651 L270:
652 
653  if ((d__1 = *x1 - *x0, fabs(d__1)) <= *hmin)
654  {
655  xroot = true;
656  }
657 
658  goto L150;
659 
660  /* Return with X1 as the root. Set JROOT. Set X = X1 and GX = G1. ----- */
661 L300:
662  *jflag = 2;
663  *x = *x1;
664  dcopy_(ng, &g1[1], &c__1, &gx[1], &c__1);
665  i__1 = *ng;
666 
667  for (i__ = 1; i__ <= i__1; ++i__)
668  {
669  jroot[i__] = 0;
670 
671  if ((d__1 = g1[i__], fabs(d__1)) > zero)
672  {
673  goto L310;
674  }
675 
676  jroot[i__] = 1;
677  goto L320;
678 L310:
679 
680  if (d_sign(c_b3, g0[i__]) != d_sign(c_b3, g1[i__]))
681  {
682  jroot[i__] = 1;
683  }
684 
685 L320:
686  ;
687  }
688 
689  return 0;
690 
691  /* No sign change in the interval. Check for zero at right endpoint. --- */
692 L400:
693 
694  if (! zroot)
695  {
696  goto L420;
697  }
698 
699  /* Zero value at X1 and no sign change in (X0,X1). Return JFLAG = 3. --- */
700  *x = *x1;
701  dcopy_(ng, &g1[1], &c__1, &gx[1], &c__1);
702  i__1 = *ng;
703 
704  for (i__ = 1; i__ <= i__1; ++i__)
705  {
706  jroot[i__] = 0;
707 
708  if ((d__1 = g1[i__], fabs(d__1)) <= zero)
709  {
710  jroot[i__] = 1;
711  }
712 
713  /* L410: */
714  }
715 
716  *jflag = 3;
717  return 0;
718 
719  /* No sign changes in this interval. Set X = X1, return JFLAG = 4. ----- */
720 L420:
721  dcopy_(ng, &g1[1], &c__1, &gx[1], &c__1);
722  *x = *x1;
723  *jflag = 4;
724  return 0;
725  /* ----------------------- End of Subroutine DROOTS ---------------------- */
726 } /* droots_ */
#define C_INT
Definition: copasi.h:115
void(* evalG)(const C_INT *, const double *, const double *, const C_INT *, double *)
Definition: common.h:32
int dcopy_(integer *n, doublereal *dx, integer *incx, doublereal *dy, integer *incy)
C_INT droots_(C_INT *ng, double *hmin, C_INT *jflag, double *x0, double *x1, double *g0, double *g1, double *gx, double *x, C_INT *jroot)
Definition: drcheck.cpp:343
static C_INT c__1
Definition: drcheck.cpp:42
C_INT dintdy_(double *t, const C_INT *k, double *yh, C_INT *nyh, double *dky, C_INT *iflag)
Definition: dintdy.cpp:53
#define dlsr01_1
Definition: drcheck.cpp:33
#define dlsr01_2
Definition: drcheck.cpp:34
switch(yytype)
static double c_b3
Definition: drcheck.cpp:39
#define dls001_1
Definition: drcheck.cpp:29
double d_sign(const double &a, const double &b)
Definition: CLSODA.cpp:2471
C_INT drchek_(const C_INT *job, evalG g, C_INT *neq, double *y, double *yh, C_INT *nyh, double *g0, double *g1, double *gx, C_INT *jroot, C_INT *irt)
Definition: drcheck.cpp:46
static const C_INT c__0
Definition: drcheck.cpp:41
#define max(a, b)
Definition: f2c.h:176