/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:33:09 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv pf=,p_ddasl2 s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include #include #include "p_ddasl2.h" #include "drddasl2.h" /* PARAMETER translations */ #define LDC NEQ #define LIW (30 + NEQ) #define LRW (45 + (5 + 2*MAXCON + 4)*NEQ + SQ(NEQ)) #define LTD NEQ #define MAXCON 0 #define NDIG 8 #define NEPOCH 10 #define NEQ 4 #define TOL (powi(10.e0,-NDIG)) /* end of PARAMETER translations */ /* COMMON translations */ struct t_counts { long int kf2, ks2, kf3, ks3; } counts; /* end of COMMON translations */ int main( ) { LOGICAL32 match; long int _l0, i, idid, info[16], iwork[LIW], j; double atol[NEQ], c[LTD][LDC], ftol, rnktol, rtol[NEQ], rwork[LRW], t, tout[NEPOCH], y[NEQ], y1[NEPOCH][NEQ], y2[NEPOCH][NEQ], yprime[NEQ]; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Atol = &atol[0] - 1; long *const Info = &info[0] - 1; long *const Iwork = &iwork[0] - 1; double *const Rtol = &rtol[0] - 1; double *const Rwork = &rwork[0] - 1; double *const Tout = &tout[0] - 1; double *const Y = &y[0] - 1; double *const Yprime = &yprime[0] - 1; /* end of OFFSET VECTORS */ /*>> 2009-10-27 DRDDASL2 Krogh Declared DUM(4,*) for NAG compiler. *>> 2009-10-19 DRDDASL2 Krogh Changed def. of d(4,4) so 0 doesn't abort. *>> 2008-10-26 DRDDASL2 Krogh Moved Format statements up for C conv. *>> 2008-10-24 DRDDASL2 Krogh Removed in INCLUDE statement & cDEC$... *>> 2008-09-04 DRDDASL2 Hanson added starting computation of y' *>> 2008-08-26 DRDDASL2 Hanson added row dimensions to evaluators *>> 2001-10-11 DRDDASL2 R. J. Hanson Example 2 Code, with Download */ /* Solve Enright and Pryce stiff test problem E5. * The equation is presented as y'= F(t,y), y_0 given. * This is solved by defining the residual function * f(t,y,y')=F(t,y)-y'. Integration is done twice. */ /* The first time analytic partials are provided. */ /* The second time, difference quotients are used for * partials, with the integrator requesting only values * of f(t,y,y'). The two solutions are equivalent but * not exactly equal. */ /*--D replaces "?": DR?DASL2, ?DASLS, ?DASLX, ?DASSF2, ?DASSF3, ?COPY *-- & ?ROTG */ /* Set number of equations: */ /* Set number of constraints. */ /* Work space sizes: */ /*++S Default NDIG = 3 *++ Default NDIG = 8 *++ Substitute for NDIG below */ /* Tolerances: */ for (i = 1; i <= NEQ; i++) { Atol[i] = TOL; Rtol[i] = TOL; } /* Setup options: */ for (i = 1; i <= 16; i++) { Info[i] = 0; } /* Compute, factor and solve the partial derivative matrix in routine * DDASSF2. This illustrates how to completely control the linear * algebra. Using analytic partial derivatives. */ Info[5] = 5; /* Compute the initial value of YPRIME(*): */ t = 0; ftol = TOL; rnktol = TOL; /* Give starting values to y and y' before Newton method: */ for (i = 1; i <= NEQ; i++) { Y[i] = 0.e0; Yprime[i] = 0.e0; } Y[1] = 1.76e-3; ddasls( ddassf2, NEQ, &t, y, yprime, info, ftol, rnktol, (double*)c, ADR(_l0,LDC), LTD, &idid, rwork, LRW, iwork, LIW ); /* Set output points for integration: */ Tout[1] = 0.e0; Tout[2] = 1.e-3; for (i = 3; i <= NEPOCH; i++) { Tout[i] = 10*Tout[i - 1]; } for (i = 2; i <= NEPOCH; i++) { t = Tout[i - 1]; ddaslx( ddassf2, NEQ, &t, y, yprime, Tout[i], info, rtol, atol, &idid, rwork, LRW, iwork, LIW ); dcopy( NEQ, y, 1, &y1[i - 1][0], 1 ); } /* Compute the initial value of YPRIME(*): */ t = 0; ftol = TOL; rnktol = TOL; /* Give starting values to y and y' before Newton method: */ for (i = 1; i <= NEQ; i++) { Y[i] = 0.e0; Yprime[i] = 0.e0; } Y[1] = 1.76e-3; ddasls( ddassf3, NEQ, &t, y, yprime, info, ftol, rnktol, (double*)c, ADR(_l0,LDC), LTD, &idid, rwork, LRW, iwork, LIW ); /* Restart the integration, but do not use analytic derivatives. */ Info[1] = 0; /* Use divided differences instead of user-provided partials. * And the default linear solver provided. */ Info[5] = 1; for (i = 2; i <= NEPOCH; i++) { t = Tout[i - 1]; ddaslx( ddassf3, NEQ, &t, y, yprime, Tout[i], info, rtol, atol, &idid, rwork, LRW, iwork, LIW ); dcopy( NEQ, y, 1, &y2[i - 1][0], 1 ); } /* See if the solutions "match" the requested tolerances. */ match = TRUE; for (j = 2; j <= NEPOCH; j++) { for (i = 1; i <= NEQ; i++) { y2[j - 1][i - 1] = (y2[j - 1][i - 1] - y1[j - 1][i - 1])/ (fabs( y1[j - 1][i - 1] )*Rtol[i] + Atol[i]); match = match && (fabs( y2[j - 1][i - 1] ) <= 100.e0); } } /* Output the relative errors and a summary. */ printf("\n Example Results for a Stiff ODE Problem, E5\n"); printf("\n T RErr in y_1 RErr in y_2 RErr in y_3 RErr in y_4\n"); for (j = 2; j <= NEPOCH; j++) { printf("%14.4e", Tout[j]); for (i = 1; i <= NEQ; i++) { printf("%14.4e", y2[j - 1][i - 1]); } printf("\n"); } if (match) { printf("\n No. of Residual Evaluations No. of User Solves\n Partials Provided- %4ld %4ld\n", counts.kf2, counts.ks2); printf(" Divided Differences- %4ld %4ld\n", counts.kf3, counts.ks3); } else { printf(" Integration Values do not Match\n"); } exit(0); } /* end of function */ void /*FUNCTION*/ ddassf2( double *t, double y[], double yprime[], double delta[], double *dum, long *ldd, double *cj, long *ires, double rwork[], long iwork[]) { #define DUM(I_,J_) (*(dum+(I_)*(*ldd)+(J_))) long int i; double u; static double b0, b1, b2, b3, b4, d[4][4], dc[3], ds[3]; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Dc = &dc[0] - 1; double *const Delta = &delta[0] - 1; double *const Ds = &ds[0] - 1; long *const Iwork = &iwork[0] - 1; double *const Rwork = &rwork[0] - 1; double *const Y = &y[0] - 1; double *const Yprime = &yprime[0] - 1; /* end of OFFSET VECTORS */ /* Routine for the Enright and Pryce problem E5. * DUM(*,*) is not used in this example. */ /* This is the setup call. */ if (*ires == 0) { b0 = 1.76e-03; b1 = 7.89e-10; b2 = 1.1e07; b3 = 1.13e09; b4 = 1.13e03; /* Set initial conditions. */ Y[1] = 0.e0; dcopy( 4, y, 0, y, 1 ); Y[1] = b0; /* Count evaluations in KF2 and KS2. */ counts.kf2 = 0; counts.ks2 = 0; return; } /* The sytem residual value: */ if (*ires == 1) { Delta[1] = -b1*Y[1] - b2*Y[1]*Y[3] - Yprime[1]; Delta[2] = b1*Y[1] - b3*Y[2]*Y[3] - Yprime[2]; Delta[3] = b1*Y[1] - b2*Y[1]*Y[3] + b4*Y[4] - b3*Y[2]*Y[3] - Yprime[3]; Delta[4] = b2*Y[1]*Y[3] - b4*Y[4] - Yprime[4]; counts.kf2 += 1; return; } /* The partial derivative matrix: */ if (*ires == 2) { d[0][3] = b2*Y[3]; d[0][0] = -b1 - d[0][3] - *cj; d[0][1] = b1; d[0][2] = b1 - d[0][3]; d[1][0] = 0.e0; d[1][2] = -b3*Y[3]; d[1][1] = d[1][2] - *cj; d[1][3] = 0.e0; d[2][0] = -b2*Y[1]; d[2][1] = -b3*Y[2]; d[2][2] = d[2][0] + d[2][1] - *cj; d[2][3] = -d[2][0]; d[3][0] = 0.e0; d[3][1] = 0.e0; d[3][2] = b4; d[3][3] = -b4 - *cj; /* This matrix is right factored to lower triangular for with three * plane rotations. Use planes 1 and 3: */ /* Tell integrator system is non-singular. */ *ires = 0; drotg( &d[0][0], &d[2][0], &Dc[1], &Ds[1] ); u = Dc[1]*d[0][1] + Ds[1]*d[2][1]; d[2][1] = -Ds[1]*d[0][1] + Dc[1]*d[2][1]; d[0][1] = u; u = Dc[1]*d[0][2] + Ds[1]*d[2][2]; d[2][2] = -Ds[1]*d[0][2] + Dc[1]*d[2][2]; d[0][2] = u; u = Dc[1]*d[0][3] + Ds[1]*d[2][3]; d[2][3] = -Ds[1]*d[0][3] + Dc[1]*d[2][3]; d[0][3] = u; d[0][0] = 1.e0/d[0][0]; /* Eliminate in planes 2 and 3. */ drotg( &d[1][1], &d[2][1], &Dc[2], &Ds[2] ); u = Dc[2]*d[1][2] + Ds[2]*d[2][2]; d[2][2] = -Ds[2]*d[1][2] + Dc[2]*d[2][2]; d[1][2] = u; u = Dc[2]*d[1][3] + Ds[2]*d[2][3]; d[2][3] = -Ds[2]*d[1][3] + Dc[2]*d[2][3]; d[1][3] = u; d[1][1] = 1.e0/d[1][1]; /* Eliminate in planes 3 and 4. */ drotg( &d[2][2], &d[3][2], &Dc[3], &Ds[3] ); u = Dc[3]*d[2][3] + Ds[3]*d[3][3]; d[3][3] = -Ds[3]*d[2][3] + Dc[3]*d[3][3]; if (d[3][3] == 0.e0) { *ires = 1; return; } d[3][3] = 1.e0/d[3][3]; d[2][3] = u; d[2][2] = 1.e0/d[2][2]; return; } /* Solve the corrector equation. */ if (*ires == 4) { /* Count number of solves. */ counts.ks2 += 1; /* Forward substitute: */ Delta[1] *= d[0][0]; Delta[4] += -Delta[1]*d[0][3]; Delta[3] += -Delta[1]*d[0][2]; Delta[2] += -Delta[1]*d[0][1]; Delta[2] *= d[1][1]; Delta[4] += -Delta[2]*d[1][3]; Delta[3] += -Delta[2]*d[1][2]; Delta[3] *= d[2][2]; Delta[4] += -Delta[3]*d[2][3]; Delta[4] *= d[3][3]; /* Apply three plane rotations,(3,4),(2,3),(1,3). */ u = Dc[3]*Delta[3] - Ds[3]*Delta[4]; Delta[4] = Ds[3]*Delta[3] + Dc[3]*Delta[4]; Delta[3] = u; u = Dc[2]*Delta[2] - Ds[2]*Delta[3]; Delta[3] = Ds[2]*Delta[2] + Dc[2]*Delta[3]; Delta[2] = u; u = Dc[1]*Delta[1] - Ds[1]*Delta[3]; Delta[3] = Ds[1]*Delta[1] + Dc[1]*Delta[3]; Delta[1] = u; return; } /* This is df/dy' for the starting procedure. This * problem is a linear ODE. */ if (*ires == 7) { for (i = 1; i <= 4; i++) { d[i - 1][i - 1] = -1.e0; } return; } return; #undef DUM } /* end of function */ /* PARAMETER translations */ #define B1 7.89e-10 #define B2 1.1e07 #define B3 1.13e09 #define B4 1.13e03 /* end of PARAMETER translations */ void /*FUNCTION*/ ddassf3( double *t, double y[], double yprime[], double delta[], double *d, long *ldd, double *cj, long *ires, double rwork[], long iwork[]) { #define D(I_,J_) (*(d+(I_)*(*ldd)+(J_))) long int i; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Delta = &delta[0] - 1; long *const Iwork = &iwork[0] - 1; double *const Rwork = &rwork[0] - 1; double *const Y = &y[0] - 1; double *const Yprime = &yprime[0] - 1; /* end of OFFSET VECTORS */ /* Routine for the Enright and Pryce problem E5. */ /* This is the setup call. */ if (*ires == 0) { /* Count evaluations in KF3 and KS3. */ counts.kf3 = 0; counts.ks3 = 0; } /* The sytem residual value: */ if (*ires == 1) { Delta[1] = -B1*Y[1] - B2*Y[1]*Y[3] - Yprime[1]; Delta[2] = B1*Y[1] - B3*Y[2]*Y[3] - Yprime[2]; Delta[3] = B1*Y[1] - B2*Y[1]*Y[3] + B4*Y[4] - B3*Y[2]*Y[3] - Yprime[3]; Delta[4] = B2*Y[1]*Y[3] - B4*Y[4] - Yprime[4]; counts.kf3 += 1; } /* The partial derivative matrix: */ if (*ires == 2) { D(0,3) = B2*Y[3]; D(0,0) = -B1 - D(0,3) - *cj; D(0,1) = B1; D(0,2) = B1 - D(0,3); D(1,1) = -B3*Y[3] - *cj; D(1,2) = D(1,1); D(2,0) = -B2*Y[1]; D(2,1) = -B3*Y[2]; D(2,2) = D(2,0) + D(2,1) - *cj; D(2,3) = -D(2,0); D(3,2) = B4; D(3,3) = -B4 - *cj; } if (*ires == 3) { printf("Should be no call to DDASSF3 with IRES = 3.\n"); } if (*ires == 4) { printf("Should be no call to DDASSF3 with IRES = 4.\n"); } /* This is df/dy' for the starting procedure. This * problem is a linear ODE of index 0. So only IRES=1,6 occur. */ if (*ires == 7) { for (i = 1; i <= 4; i++) { D(i - 1,i - 1) = -1.e0; } return; } return; #undef D } /* end of function */