27 #define dls001_1 (mdls001_._1)
28 #define dls001_2 (mdls001_._2)
29 #define dls001_3 (mdls001_._3)
31 #define dlsa01_1 (mdlsa01_._1)
32 #define dlsa01_2 (mdlsa01_._2)
33 #define dlsa01_3 (mdlsa01_._3)
38 #define pow_dd(__x, __y) pow(*__x, *__y)
46 C_INT *nyh,
double *yh1,
double *ewt,
double *savf,
47 double *acor,
double *wm,
C_INT *iwm,
52 static const double sm1[12] =
54 .5, .575, .55, .45, .35, .25, .2, .15, .1, .075, .05,
59 C_INT yh_dim1, yh_offset, i__1, i__2;
60 double d__1, d__2, d__3;
66 double rh = 0.0, rm, dm1, dm2;
68 double rh1, rh2, del, ddn;
70 double pdh = 0.0, dsm, dup, exm1, exm2;
78 double exsm, rhup, rate, exup, rh1it, alpha;
86 yh_offset = 1 + yh_dim1;
229 for (i__ = 1; i__ <= 5; ++i__)
233 i__ + 1 + i__ * 13 - 14];
238 for (i__ = 1; i__ <= 12; ++i__)
242 i__ + 1 + i__ * 13 - 14];
281 for (i__ = 1; i__ <= i__1; ++i__)
345 if (rh * pdh * 1.00001 < sm1[
dls001_3.nq - 1])
356 for (j = 2; j <= i__1; ++j)
361 for (i__ = 1; i__ <= i__2; ++i__)
364 yh[i__ + j * yh_dim1] *= r__;
401 for (jb = 1; jb <= i__2; ++jb)
407 for (i__ = i1; i__ <= i__1; ++i__)
410 yh1[i__] += yh1[i__ + *nyh];
429 for (i__ = 1; i__ <= i__2; ++i__)
432 y[i__] = yh[i__ + yh_dim1];
435 f(&neq[1], &
dls001_3.tn, &y[1], &savf[1]);
448 (*pjac)(&neq[1], &y[1], &yh[yh_offset], nyh, &ewt[1], &acor[1], &savf[1],
449 &wm[1], &iwm[1], f, jac);
463 for (i__ = 1; i__ <= i__2; ++i__)
482 for (i__ = 1; i__ <= i__2; ++i__)
484 savf[i__] =
dls001_3.h__ * savf[i__] - yh[i__ + (yh_dim1 << 1)];
486 y[i__] = savf[i__] - acor[i__];
492 for (i__ = 1; i__ <= i__2; ++i__)
494 y[i__] = yh[i__ + yh_dim1] +
dls001_3.el[0] * savf[i__];
496 acor[i__] = savf[i__];
508 for (i__ = 1; i__ <= i__2; ++i__)
511 y[i__] =
dls001_3.h__ * savf[i__] - (yh[i__ + (yh_dim1 << 1)] + acor[
515 (*slvs)(&wm[1], &iwm[1], &y[1], &savf[1]);
530 for (i__ = 1; i__ <= i__2; ++i__)
534 y[i__] = yh[i__ + yh_dim1] +
dls001_3.el[0] * acor[i__];
551 if (del <= pnorm * 100. *
dls001_3.uround)
568 if (del <= delp * 1024.)
579 d__1 = 1., d__2 =
dls001_3.crate * 1.5;
607 if (m >= 2 && del > delp * 2.)
613 f(&neq[1], &
dls001_3.tn, &y[1], &savf[1]);
641 for (jb = 1; jb <= i__2; ++jb)
647 for (i__ = i1; i__ <= i__1; ++i__)
650 yh1[i__] -= yh1[i__ + *nyh];
721 for (j = 1; j <= i__2; ++j)
725 for (i__ = 1; i__ <= i__1; ++i__)
728 yh[i__ + j * yh_dim1] +=
dls001_3.el[j - 1] * acor[i__];
788 rh1 = 1. / (
pow_dd(&dsm, &exsm) * 1.2 + 1.2e-6);
792 if (pdh * rh1 > 1e-5)
810 rh2 = 1. / (
pow_dd(&dm2, &exm2) * 1.2 + 1.2e-6);
815 rh2 = 1. / (
pow_dd(&dm2, &exsm) * 1.2 + 1.2e-6);
858 rh1 = 1. / (
pow_dd(&dm1, &exm1) * 1.2 + 1.2e-6);
863 rh1 = 1. / (
pow_dd(&dm1, &exsm) * 1.2 + 1.2e-6);
870 if (pdh * rh1 > 1e-5)
872 rh1it = sm1[nqm1 - 1] / pdh;
876 rh2 = 1. / (
pow_dd(&dsm, &exsm) * 1.2 + 1.2e-6);
878 if (rh1 *
dlsa01_1.ratio < rh2 * 5.)
884 dm1 =
pow_dd(&alpha, &exm1) * dm1;
886 if (dm1 <=
dls001_3.uround * 1e3 * pnorm)
922 for (i__ = 1; i__ <= i__1; ++i__)
925 yh[i__ +
dls001_3.lmax * yh_dim1] = acor[i__];
942 for (jb = 1; jb <= i__1; ++jb)
948 for (i__ = i1; i__ <= i__2; ++i__)
951 yh1[i__] -= yh1[i__ + *nyh];
991 for (i__ = 1; i__ <= i__1; ++i__)
994 savf[i__] = acor[i__] - yh[i__ +
dls001_3.lmax * yh_dim1];
1000 rhup = 1. / (
pow_dd(&dup, &exup) * 1.4 + 1.4e-6);
1003 rhsm = 1. / (
pow_dd(&dsm, &exsm) * 1.2 + 1.2e-6);
1014 rhdn = 1. / (
pow_dd(&ddn, &exdn) * 1.3 + 1.3e-6);
1030 d__1 = rhup, d__2 = sm1[
dls001_3.l - 1] / pdh;
1035 d__1 = rhsm, d__2 = sm1[
dls001_3.nq - 1] / pdh;
1041 d__1 = rhdn, d__2 = sm1[
dls001_3.nq - 2] / pdh;
1091 for (i__ = 1; i__ <= i__1; ++i__)
1094 yh[i__ + (newq + 1) * yh_dim1] = acor[i__] * r__;
1109 if (rh * pdh * 1.00001 >= sm1[newq - 1])
1116 if (
dls001_3.kflag == 0 && rh < 1.1)
1170 for (i__ = 1; i__ <= i__1; ++i__)
1173 y[i__] = yh[i__ + yh_dim1];
1176 f(&neq[1], &
dls001_3.tn, &y[1], &savf[1]);
1180 for (i__ = 1; i__ <= i__1; ++i__)
1183 yh[i__ + (yh_dim1 << 1)] =
dls001_3.h__ * savf[i__];
1217 for (i__ = 1; i__ <= i__1; ++i__)
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)
void(* evalF)(const C_INT *, const double *, const double *, double *)
void(* evalJ)(const C_INT *, const double *, const double *, const C_INT *, const C_INT *, double *, const C_INT *)
C_INT dcfode_(C_INT *meth, double *elco, double *tesco)
double dmnorm_(C_INT *n, double *v, double *w)