/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:33:15 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv pf=,p_sdasl4 s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include #include #include "p_sdasl4.h" #include "drsdasl4.h" /* PARAMETER translations */ #define LIW (30 + NEQ) #define LRW (45 + (5 + 2*MAXCON + 4)*NEQ) #define MAXCON 0 #define NDIG 3 #define NEPOCH 10 #define NEQ (NSTATE + NVAREQ*NSTATE) #define NSTATE 4 #define NVAREQ 5 #define TOL (powif(10.e0,-NDIG)) /* end of PARAMETER translations */ int main( ) { long int _n, i, idid, ier, info[16], ipvt[NSTATE], ires, ival[3], iwork[LIW], j, k, l; float atol[NEQ], b0, b1, b2, b3, b4, cj, d[NSTATE][NSTATE], delta[NSTATE], gp[NVAREQ][NSTATE], gy[NSTATE][NSTATE], rtol[NEQ], rwork[LRW], t, tout[NEPOCH], y[NEQ], yp[NVAREQ][NSTATE], ypp[NVAREQ][NSTATE], yprime[NEQ]; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Atol = &atol[0] - 1; float *const Delta = &delta[0] - 1; long *const Info = &info[0] - 1; long *const Ipvt = &ipvt[0] - 1; long *const Ival = &ival[0] - 1; long *const Iwork = &iwork[0] - 1; float *const Rtol = &rtol[0] - 1; float *const Rwork = &rwork[0] - 1; float *const Tout = &tout[0] - 1; float *const Y = &y[0] - 1; float *const Yprime = &yprime[0] - 1; /* end of OFFSET VECTORS */ /*>> 2008-10-29 DRSDASL4 Krogh Fixed print output for C. *>> 2008-10-26 DRSDASL4 Krogh Changed Fortran 90 type statement to F77. *>> 2008-10-26 DRSDASL4 Krogh Moved Format statements up for C conv. *>> 2008-10-24 DRSDASL4 Krogh Removed in INCLUDE statement & cDEC$... *>> 2008-09-17 DRSDASL4 Hanson added starting values for y' *>> 2008-08-26 DRSDASL4 Hanson added row dimension to evaluator *>> 2001-10-11 DRSDASL4 R. J. Hanson Example 4 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'. Variational equations, with respect * to model parameters b0,b1,b2,b3, and b4, are integrated together * with state equations. The block structure of the linear algebra * solve step is used. The first NSTATE are state equations, * and the remaining blocks are for the variational equations. */ /* Reverse communication is used for the evaluation and * linear solve steps. */ /*--S replaces "?": DR?DASL4, ?DASLX, ?DASSF4, ?COPY, ?GEFA, ?GESL, ?DOT */ /* Set number of equations: */ /* Set number of constraints. */ /* Work space sizes: */ /*.....The work space does not need the NEQ**2-size additional *.....array piece because linear algebra is done in user work space. * This is flagged by INFO(5)=-6. */ /*++S Default NDIG = 3 *++ Default NDIG = 8 *++ Substitute for NDIG below */ /* Tolerances: */ for (i = 1; i <= NSTATE; i++) { Atol[i] = TOL; Rtol[i] = TOL; } /* Put more emphasis on relative error for the variational equations. * (How to set these tolerances depends on the problem.) */ for (i = NSTATE + 1; i <= NEQ; i++) { Atol[i] = 0; Rtol[i] = TOL; } /* Setup options: */ for (i = 1; i <= 16; i++) { Info[i] = 0; } /* Use partial derivatives and reverse communication * Do the linear algebra ourselves */ Info[5] = -6; /* 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]; } /* Have reverse communiation for f, and partials, * solver, and storage of partials. */ b0 = 1.76e-03; b1 = 7.89e-10; b2 = 1.1e07; b3 = 1.13e09; b4 = 1.13e03; /* Set initial conditions for state and variational equations: */ Y[1] = 0.e0; Yprime[1] = 0.e0; scopy( NEQ, y, 0, y, 1 ); scopy( NEQ, yprime, 0, yprime, 1 ); Y[1] = b0; Y[NSTATE + 1] = 1.e0; /* Set the initial value of YPRIME(*): */ Yprime[1] = -b1*b0; Yprime[2] = b1*b0; Yprime[3] = b1*b0; Yprime[NSTATE + 1] = -b1; Yprime[NSTATE + 2] = b1; Yprime[NSTATE + 3] = b1; Yprime[2*NSTATE + 1] = -b0; Yprime[2*NSTATE + 2] = b0; Yprime[2*NSTATE + 3] = b0; /* Count residuals, factorizations and solves. */ Ival[1] = 0; Ival[2] = 0; Ival[3] = 0; /* Print heading. */ printf(" Example Results for a Stiff ODE Problem, E5.\n Five additional variational equations are integrated.\n\n T y_1 y_2 y_3 y_4\n"); /* Allow more than nominal number of steps. */ Info[12] = 5000; for (l = 2; l <= NEPOCH; l++) { t = Tout[l - 1]; /* Branch here to re-enter routine. */ L_50: sdaslx( sdassf4, NEQ, &t, y, yprime, Tout[l], info, rtol, atol, &idid, rwork, LRW, iwork, LIW ); /* When IDID == 4 the code returned for reverse communication. * Otherwise the integration has finished up to this output point. */ if (idid != 4) goto L_100; /* Get flag for user action: */ ires = Iwork[3]; /* Will need partials for IRES=1,2: */ if (ires <= 2) { gy[0][0] = 0.e0; gp[0][0] = 0.e0; scopy( SQ(NSTATE), (float*)gy, 0, (float*)gy, 1 ); scopy( NSTATE*NVAREQ, (float*)gp, 0, (float*)gp, 1 ); /* Evaluate DG/DY matrix: */ gy[0][3] = b2*Y[3]; gy[0][0] = -b1 - gy[0][3]; gy[0][1] = b1; gy[0][2] = b1 - gy[0][3]; gy[1][0] = 0.e0; gy[1][2] = -b3*Y[3]; gy[1][1] = gy[1][2]; gy[1][3] = 0.e0; gy[2][0] = -b2*Y[1]; gy[2][1] = -b3*Y[2]; gy[2][2] = gy[2][0] + gy[2][1]; gy[2][3] = -gy[2][0]; gy[3][0] = 0.e0; gy[3][1] = 0.e0; gy[3][2] = b4; gy[3][3] = -b4; /* Evaluate DG/DP matrix: */ gp[1][0] = -Y[1]; gp[1][1] = Y[1]; gp[1][2] = Y[1]; gp[2][0] = -Y[1]*Y[3]; gp[2][2] = -Y[1]*Y[3]; gp[3][1] = -Y[2]*Y[3]; gp[3][2] = Y[2]*Y[3]; gp[4][2] = Y[4]; gp[4][3] = -Y[4]; /* Extract DY/DP: */ scopy( NSTATE*NVAREQ, &Y[NSTATE + 1], 1, (float*)yp, 1 ); /* Extract (DY/DP)': */ scopy( NSTATE*NVAREQ, &Yprime[NSTATE + 1], 1, (float*)ypp, 1 ); } /* Evaluate state and variational equations. */ if (ires == 1) { /* Count residual evaluations. */ Ival[1] += 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]; /* Compute matrix arithmetic, DG/DY * DY/DP + DG/DP: * Update DG/DP = DG/DP - D (Y')/DP: */ for (j = 1; j <= NVAREQ; j++) { for (i = 1; i <= NSTATE; i++) { gp[j - 1][i - 1] += sdot( NSTATE, &gy[0][i - 1], NSTATE, &yp[j - 1][0], 1 ) - ypp[j - 1][i - 1]; } } /* Move the results where the integrator uses them: */ scopy( NSTATE, delta, 1, &Rwork[Iwork[4]], 1 ); scopy( NSTATE*NVAREQ, (float*)gp, 1, &Rwork[Iwork[4] + NSTATE], 1 ); goto L_50; } if (ires == 2) { /* Evaluate partial matrix, DG/DY-C*I_4. * Each of the NSTATE vectors has this * matrix as a diagonal block. The C value * is located in RWORK(1). * There are NVAREQ+1 such vectors. * Taking advantage of this fact allows the * corrector equation to be factored and solved * using a system of size equal to the number of * state vector components. */ scopy( SQ(NSTATE), (float*)gy, 1, (float*)d, 1 ); cj = Rwork[1]; for (i = 1; i <= NSTATE; i++) { d[i - 1][i - 1] -= cj; } /* Factor the 4 by 4 matrix using the Linpack code. * Count factorizations. */ Ival[2] += 1; /* This is where savings occur, compared to solving a full system. */ sgefa( (float*)d, NSTATE, NSTATE, ipvt, &ier ); /* This passes the error code, after factorization, * back to the integrator. */ Iwork[3] = ier; goto L_50; } /* Solve with NVAREQ+1 right-hand sides using the Linpack code. * This is where savings occur, compared to solving a full system. */ if (ires == 4) { /* Count solves. */ Ival[3] += 1; /* Index to place solutions of corrector equations. */ k = Iwork[4]; for (j = 1; j <= (NVAREQ + 1); j++) { sgesl( (float*)d, NSTATE, NSTATE, ipvt, &Rwork[k], 0 ); k += NSTATE; } goto L_50; } L_100: ; /* Print time, state and variational equation results. */ printf("%14.4e%14.4e%14.4e%14.4e%14.4e\n", Tout[l], y[0], y[1], y[2], y[3]); for(_n=4L; _n < sizeof(y)/sizeof(float ); _n+=4) printf(" %14.4e%14.4e%14.4e%14.4e\n", y[_n], y[_n+1], y[_n+2], y[_n+3]); } /* WRITE (*,'(1P,5D14.4/(14x,4D14.4))') tout(L),y */ printf(" No. of Residual Evaluations Factorizations No. of User Solves\nReverse Communication-"); for(_n=0L; _n < sizeof(ival)/sizeof(long); _n++) printf("%9ld", ival[_n]); printf(" \n"); exit(0); } /* end of function */ void /*FUNCTION*/ sdassf4( float *t, float y[], float yprime[], float delta[], float *d, long *ldd, float *cj, long *ires, float rwork[], long iwork[]) { #define D(I_,J_) (*(d+(I_)*(*ldd)+(J_))) /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Delta = &delta[0] - 1; long *const Iwork = &iwork[0] - 1; float *const Rwork = &rwork[0] - 1; float *const Y = &y[0] - 1; float *const Yprime = &yprime[0] - 1; /* end of OFFSET VECTORS */ /* Dummy routine for reverse communication: * This is a double precision version. */ /* This is the setup call. It can be ignored since * all problem information is provided using reverse communication. */ if (*ires == 0) { return; } /* Other values of IRES will not occur. */ return; #undef D } /* end of function */