33 double d_sign(
const double & a,
const double & b);
38 #define dls001_1 (mdls001_.lsoda)
39 #define dlsa01_1 (mdlsa01_.lsoda)
40 #define dlsr01_1 (mdlsr01_.lsodar)
42 static const double c_b76 = 0.0;
125 *t,
double *tout,
C_INT *itol,
double *rtol,
double *
158 C_INT len1c, len1n, len1s, iflag;
160 C_INT leniw, lenwm = 0, lenyh, imxer;
163 double rtoli, tdist, tolsf, tnext;
1310 if (*istate < 1 || *istate > 3)
1315 if (*itask < 1 || *itask > 5)
1375 if (*itol < 1 || *itol > 4)
1380 if (*iopt < 0 || *iopt > 1)
1385 if (*jt == 3 || *jt < 1 || *jt > 5)
1516 if ((*tout - *t) * h0 < 0.)
1587 if (*lrw < lyhnew - 1 + lenyh)
1613 lenwm = ((ml << 1) + mu + 1) *
dls001_1.n + 2;
1626 lenrw = len1 + len2;
1627 lenrwc = len1c + len2;
1640 if (*istate == 1 && *lrw < lenrwc)
1645 if (*istate == 1 && *liw < leniwc)
1650 if (*istate == 3 && *lrw < lenrwc)
1655 if (*istate == 3 && *liw < leniwc)
1670 msg =
"DLSODAR- Warning.. RWORK length is sufficient for now, but ";
1673 msg =
" may not be later. Integration will proceed anyway. ";
1676 msg =
" Length needed is LENRW = I1, while LRW = I2.";
1690 msg =
"DLSODAR- Warning.. IWORK length is sufficient for now, but ";
1693 msg =
" may not be later. Integration will proceed anyway. ";
1696 msg =
" Length needed is LENIW = I1, while LIW = I2.";
1705 for (i__ = 1; i__ <= i__1; ++i__)
1712 if (*itol == 2 || *itol == 4)
1754 for (i__ = i1; i__ <= i__1; ++i__)
1769 dls001_1.uround = std::numeric_limits< C_FLOAT64 >::epsilon();
1774 if (*itask != 4 && *itask != 5)
1781 if ((tcrit - *tout) * (*tout - *t) < 0.)
1786 if (h0 != 0. && (*t + h0 - tcrit) * h0 > 0.)
1807 (*f)(&neq[1], t, &y[1], &rwork[lf0]);
1812 for (i__ = 1; i__ <= i__1; ++i__)
1815 rwork[i__ +
dls001_1.lyh - 1] = y[i__];
1825 for (i__ = 1; i__ <= i__1; ++i__)
1827 if (rwork[i__ +
dls001_1.lewt - 1] <= 0.)
1867 tdist = (d__1 = *tout - *t, fabs(d__1));
1869 d__1 = fabs(*t), d__2 = fabs(*tout);
1872 if (tdist <
dls001_1.uround * 2. * w0)
1886 for (i__ = 1; i__ <= i__1; ++i__)
1890 d__1 = tol, d__2 = rtol[i__];
1904 for (i__ = 1; i__ <= i__1; ++i__)
1906 if (*itol == 2 || *itol == 4)
1911 ayi = (d__1 = y[i__], fabs(d__1));
1916 d__1 = tol, d__2 = atoli / ayi;
1925 d__1 = tol, d__2 =
dls001_1.uround * 100.;
1931 sum = 1. / (tol * w0 * w0) + tol * (d__1 * d__1);
1932 h0 = 1. / sqrt(sum);
1949 for (i__ = 1; i__ <= i__1; ++i__)
1952 rwork[i__ + lf0 - 1] = h0 * rwork[i__ + lf0 - 1];
1993 if (*itask == 1 || *itask == 4)
2047 if ((tp - *tout) *
dls001_1.h__ > 0.)
2067 if ((tcrit - *tout) *
dls001_1.h__ < 0.)
2096 ihit = (d__1 =
dls001_1.tn - tcrit, fabs(d__1)) <=
dls001_1.uround * 100. *
2116 if ((tnext - tcrit) *
dls001_1.h__ <= 0.)
2123 if (*istate == 2 &&
dls001_1.jstart >= 0)
2167 for (i__ = 1; i__ <= i__1; ++i__)
2169 if (rwork[i__ +
dls001_1.lewt - 1] <= 0.)
2209 msg =
"DLSODAR- Warning..Internal T(=R1) and H(=R2) are ";
2212 msg =
" such that in the machine, T + H = T on the next step ";
2215 msg =
" (H = step size). Solver will continue anyway.";
2224 msg =
"DLSODAR- Above warning has been issued I1 times. ";
2227 msg =
" It will not be issued again for this problem.";
2289 msg =
"DLSODAR- A switch to the BDF (stiff) method has occurred";
2296 msg =
"DLSODAR- A switch to the Adams (nonstiff) method occurred";
2301 msg =
" at T = R1, tentative step size H = R2, step NST = I1";
2368 ihit = (d__1 =
dls001_1.tn - tcrit, fabs(d__1)) <=
dls001_1.uround * 100. *
2378 if ((tnext - tcrit) *
dls001_1.h__ <= 0.)
2394 ihit = (d__1 =
dls001_1.tn - tcrit, fabs(d__1)) <=
dls001_1.uround * 100. *
2406 for (i__ = 1; i__ <= i__1; ++i__)
2409 y[i__] = rwork[i__ +
dls001_1.lyh - 1];
2414 if (*itask != 4 && *itask != 5)
2451 msg =
"DLSODAR- At current T (=R1), MXSTEP (=I1) steps ";
2454 msg =
" taken on this call before reaching TOUT ";
2461 ewti = rwork[
dls001_1.lewt + i__ - 1];
2462 msg =
"DLSODAR- At T(=R1), EWT(I1) has become R2 .le. 0.";
2469 msg =
"DLSODAR- At T (=R1), too much accuracy requested ";
2472 msg =
" for precision of machine.. See TOLSF (=R2) ";
2480 msg =
"DLSODAR- At T(=R1), step size H(=R2), the error ";
2483 msg =
" test failed repeatedly or with ABS(H) = HMIN";
2490 msg =
"DLSODAR- At T (=R1) and step size H (=R2), the ";
2493 msg =
" corrector convergence failed repeatedly ";
2496 msg =
" or with ABS(H) = HMIN ";
2503 msg =
"DLSODAR- At current T(=R1), RWORK length too small";
2506 msg =
" to proceed. The integration was otherwise successful.";
2513 msg =
"DLSODAR- At current T(=R1), IWORK length too small";
2516 msg =
" to proceed. The integration was otherwise successful.";
2527 for (i__ = 1; i__ <= i__1; ++i__)
2529 size = (d__1 = rwork[i__ +
dls001_1.lacor - 1] * rwork[i__ +
2548 for (i__ = 1; i__ <= i__1; ++i__)
2551 y[i__] = rwork[i__ +
dls001_1.lyh - 1];
2577 msg =
"DLSODAR- ISTATE(=I1) illegal.";
2588 msg =
"DLSODAR- ITASK (=I1) illegal.";
2593 msg =
"DLSODAR- ISTATE.gt.1 but DLSODAR not initialized.";
2598 msg =
"DLSODAR- NEQ (=I1) .lt. 1 ";
2603 msg =
"DLSODAR- ISTATE = 3 and NEQ increased (I1 to I2).";
2608 msg =
"DLSODAR- ITOL (=I1) illegal. ";
2613 msg =
"DLSODAR- IOPT (=I1) illegal. ";
2618 msg =
"DLSODAR- JT (=I1) illegal. ";
2623 msg =
"DLSODAR- ML (=I1) illegal: .lt.0 or .ge.NEQ (=I2)";
2628 msg =
"DLSODAR- MU (=I1) illegal: .lt.0 or .ge.NEQ (=I2)";
2633 msg =
"DLSODAR- IXPR (=I1) illegal. ";
2638 msg =
"DLSODAR- MXSTEP (=I1) .lt. 0 ";
2643 msg =
"DLSODAR- MXHNIL (=I1) .lt. 0 ";
2648 msg =
"DLSODAR- TOUT (=R1) behind T (=R2) ";
2649 mxerrwd(msg, &
c__40, &
c__14, &
c__0, &
c__0, &
c__0, &
c__0, &
c__2, tout, t, (
2651 msg =
" Integration direction is given by H0 (=R1) ";
2656 msg =
"DLSODAR- HMAX (=R1) .lt. 0.0 ";
2661 msg =
"DLSODAR- HMIN (=R1) .lt. 0.0 ";
2666 msg =
"DLSODAR- RWORK length needed, LENRW(=I1), exceeds LRW(=I2) ";
2671 msg =
"DLSODAR- IWORK length needed, LENIW(=I1), exceeds LIW(=I2) ";
2676 msg =
"DLSODAR- RTOL(I1) is R1 .lt. 0.0 ";
2681 msg =
"DLSODAR- ATOL(I1) is R1 .lt. 0.0 ";
2686 ewti = rwork[
dls001_1.lewt + i__ - 1];
2687 msg =
"DLSODAR- EWT(I1) is R1 .le. 0.0 ";
2692 msg =
"DLSODAR- TOUT(=R1) too close to T(=R2) to start integration.";
2693 mxerrwd(msg, &
c__60, &
c__22, &
c__0, &
c__0, &
c__0, &
c__0, &
c__2, tout, t, (
2697 msg =
"DLSODAR- ITASK = I1 and TOUT (=R1) behind TCUR - HU (= R2) ";
2698 mxerrwd(msg, &
c__60, &
c__23, &
c__0, &
c__1, itask, &
c__0, &
c__2, tout, &tp,
2702 msg =
"DLSODAR- ITASK = 4 or 5 and TCRIT (=R1) behind TCUR (=R2) ";
2707 msg =
"DLSODAR- ITASK = 4 or 5 and TCRIT (=R1) behind TOUT (=R2) ";
2712 msg =
"DLSODAR- At start of problem, too much accuracy ";
2715 msg =
" requested for precision of machine.. See TOLSF (=R1) ";
2721 msg =
"DLSODAR- Trouble in DINTDY. ITASK = I1, TOUT = R1";
2726 msg =
"DLSODAR- MXORDN (=I1) .lt. 0 ";
2731 msg =
"DLSODAR- MXORDS (=I1) .lt. 0 ";
2736 msg =
"DLSODAR- NG (=I1) .lt. 0 ";
2741 msg =
"DLSODAR- NG changed (from I1 to I2) illegally, ";
2744 msg =
" i.e. not immediately after a root was found.";
2749 msg =
"DLSODAR- One or more components of g has a root ";
2752 msg =
" too near to the initial point. ";
2763 msg =
"DLSODAR- Run aborted.. apparent infinite loop. ";
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(* evalG)(const C_INT *, const double *, const double *, const C_INT *, double *)
static const C_INT mxstp0
static const C_INT c__105
int dcopy_(integer *n, doublereal *dx, integer *incx, doublereal *dy, integer *incy)
static const C_INT c__202
static const C_INT c__103
static const C_INT c__303
C_INT dintdy_(double *t, const C_INT *k, double *yh, C_INT *nyh, double *dky, C_INT *iflag)
static const C_INT c__205
void(* evalF)(const C_INT *, const double *, const double *, double *)
C_INT operator()(evalF f, C_INT *neq, double *y, double *t, double *tout, C_INT *itol, double *rtol, double *atol, C_INT *itask, C_INT *istate, C_INT *iopt, double *rwork, C_INT *lrw, C_INT *iwork, C_INT *liw, evalJ jac, C_INT *jt, evalG g, C_INT *ng, C_INT *jroot)
static const C_INT c__106
C_INT dewset_(C_INT *n, C_INT *itol, double *rtol, double *atol, double *ycur, double *ewt)
double d_sign(const double &a, const double &b)
void(* evalJ)(const C_INT *, const double *, const double *, const C_INT *, const C_INT *, double *, const C_INT *)
static const C_INT c__204
C_INT dsolsy_(double *wm, C_INT *iwm, double *x, double *tem)
static const C_INT c__101
static const C_INT c__102
static const C_INT c__207
static const C_INT c__206
static const C_INT c__107
static const double c_b76
static const C_INT mxhnl0
C_INT dprja_(C_INT *neq, double *y, double *yh, C_INT *nyh, double *ewt, double *ftem, double *savf, double *wm, C_INT *iwm, evalF f, evalJ jac)
double dmnorm_(C_INT *n, double *v, double *w)
static const C_INT c__201
static const C_INT c__203
static const C_INT mord[2]
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)
static const C_INT c__104