29 #define dls001_1 (mdls001_._1)
30 #define dls001_2 (mdls001_._2)
31 #define dls001_3 (mdls001_._3)
33 #define dlsr01_1 (mdlsr01_._1)
34 #define dlsr01_2 (mdlsr01_._2)
35 #define dlsr01_3 (mdlsr01_._3)
37 double d_sign(
const double & a,
const double & b);
47 y,
double *yh,
C_INT *nyh,
double *g0,
double *g1,
51 C_INT yh_dim1, yh_offset, i__1;
56 double x, t1, temp1, temp2;
100 yh_offset = 1 + yh_dim1;
111 for (i__ = 1; i__ <= i__1; ++i__)
134 for (i__ = 1; i__ <= i__1; ++i__)
137 if ((d__1 = g0[i__], fabs(d__1)) <= 0.)
156 for (i__ = 1; i__ <= i__1; ++i__)
159 y[i__] += temp2 * yh[i__ + (yh_dim1 << 1)];
167 for (i__ = 1; i__ <= i__1; ++i__)
170 if ((d__1 = g0[i__], fabs(d__1)) <= 0.)
202 for (i__ = 1; i__ <= i__1; ++i__)
205 if ((d__1 = g0[i__], fabs(d__1)) <= 0.)
228 for (i__ = 1; i__ <= i__1; ++i__)
231 y[i__] += temp2 * yh[i__ + (yh_dim1 << 1)];
243 for (i__ = 1; i__ <= i__1; ++i__)
245 if ((d__1 = g0[i__], fabs(d__1)) > 0.)
292 dintdy_(&t1, &
c__0, &yh[yh_offset], nyh, &y[1], &iflag);
298 for (i__ = 1; i__ <= i__1; ++i__)
301 y[i__] = yh[i__ + yh_dim1];
305 (*g)(&neq[1], &t1, &y[1], &
dlsr01_1.ngc, &g1[1]);
311 &gx[1], &x, &jroot[1]);
318 dintdy_(&x, &
c__0, &yh[yh_offset], nyh, &y[1], &iflag);
319 (*g)(&neq[1], &x, &y[1], &
dlsr01_1.ngc, &gx[1]);
332 dintdy_(&x, &
c__0, &yh[yh_offset], nyh, &y[1], &iflag);
344 double *x0,
double *x1,
double *g0,
double *g1,
345 double *gx,
double *x,
C_INT *jroot)
349 static double zero = 0.;
350 static double half = .5;
351 static double tenth = .1;
352 static double five = 5.;
361 bool xroot, zroot, sgnchg;
362 C_INT imxold, nxlast;
363 double fracsub, fracint;
457 for (i__ = 1; i__ <= i__1; ++i__)
459 if ((d__1 = g1[i__], fabs(d__1)) > zero)
474 t2 = (d__1 = g1[i__] / (g1[i__] - g0[i__]), fabs(d__1));
541 if ((d__1 =
dlsr01_2.x2 - *x0, fabs(d__1)) < half * *hmin)
543 fracint = (d__1 = *x1 - *x0, fabs(d__1)) / *hmin;
548 fracsub = half / fracint;
551 dlsr01_2.x2 = *x0 + fracsub * (*x1 - *x0);
554 if ((d__1 = *x1 -
dlsr01_2.x2, fabs(d__1)) < half * *hmin)
556 fracint = (d__1 = *x1 - *x0, fabs(d__1)) / *hmin;
561 fracsub = half / fracint;
564 dlsr01_2.x2 = *x1 - fracsub * (*x1 - *x0);
579 for (i__ = 1; i__ <= i__1; ++i__)
581 if ((d__1 = gx[i__], fabs(d__1)) > zero)
596 t2 = (d__1 = gx[i__] / (gx[i__] - g0[i__]), fabs(d__1));
653 if ((d__1 = *x1 - *x0, fabs(d__1)) <= *hmin)
667 for (i__ = 1; i__ <= i__1; ++i__)
671 if ((d__1 = g1[i__], fabs(d__1)) > zero)
704 for (i__ = 1; i__ <= i__1; ++i__)
708 if ((d__1 = g1[i__], fabs(d__1)) <= zero)
void(* evalG)(const C_INT *, const double *, const double *, const C_INT *, double *)
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)
C_INT dintdy_(double *t, const C_INT *k, double *yh, C_INT *nyh, double *dky, C_INT *iflag)
double d_sign(const double &a, const double &b)
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)