/*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_ddasl7 s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include #include #include #include "p_ddasl7.h" #include "drddasl7.h" /* PARAMETER translations */ #define LIW (NEQ + 30 + 2) #define LRW (45 + (MAXORD + MAXCON + 4)*NEQ + (2*ML + MU + 1)*NEQ + 2*(NEQ/(ML + MU + 1) + 1)) #define MAXCON 0 #define MAXORD 5 #define ML 5 #define MU 0 #define NDIG 8 #define NEQ 25 #define TOL (powi(10.e0,-NDIG)) /* end of PARAMETER translations */ int main( ) { long int i, idid, info[16], iout, ipvt[NEQ], ires, iwork[LIW], lda, nerr, nqu; double atol[1], delta[NEQ], er, erm, ero, hu, rtol[1], rwork[LRW], t, tout, y[NEQ], yprime[NEQ]; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Atol = &atol[0] - 1; double *const Delta = &delta[0] - 1; long *const Info = &info[0] - 1; long *const Ipvt = &ipvt[0] - 1; long *const Iwork = &iwork[0] - 1; double *const Rtol = &rtol[0] - 1; double *const Rwork = &rwork[0] - 1; double *const Y = &y[0] - 1; double *const Yprime = &yprime[0] - 1; /* end of OFFSET VECTORS */ /*>> 2008-10-26 DRDDASL7 Krogh Moved Format statements up for C conv. *>> 2008-10-24 DRDDASL7 Krogh Removed in INCLUDE statement & cDEC$... *>> 2008-08-27 DRDDASL7 compute initial y' value, change user codes *>> 2006-04-26 DRDDASL7, Krogh Dimensioned ATOL and RTOL. *>> 2006-04-24 DRDDASL7, Krogh Moved count initialization up. *>> 2002-01-18 DRDDASL7, R. J. Hanson Example Code */ /* DDASLX is used to solve an ODE problem, with a banded Jacobian. * Two runs with analytic partials are used. Both forward and * reverse communication usage examples are illustrated. */ /*--D replaces "?": DR?DASL7, ?DASLX, ?DASSF7, ?EDIT2, ?GBFA, ?GBSL */ /*++S Default NDIG = 4 *++ Default NDIG = 8 *++ Substitute for NDIG below */ /* The work space LRW has the banded matrix size * because with INFO(5)=4, the Jacobian matrix is stored in * banded form. (Only MAXCON = 0 is allowed with the band solver.) */ /* A constant coefficient, banded matrix: * */ for (i = 1; i <= 16; i++) { Info[i] = 0; } /* Banded matrix with user providing derivatives. */ Info[5] = 4; Iwork[1] = ML; Iwork[2] = MU; lda = 2*ML + MU + 1; /* Tolerance: */ Atol[1] = TOL; Rtol[1] = 0.e0; printf(" Demo Program for DDASLX\n\n\n Problem 7: y' = A * y , where A is a banded lower triangular matrix\n NEQ =%3d ML =%2d MU =%2d( size, bandwidths)\n RTOL =%10.1e ATOL =%10.1e( rel, abs tolerance)\n", NEQ, ML, MU, Rtol[1], Atol[1]); t = 0.0e0; for (i = 1; i <= NEQ; i++) { Y[i] = 0.0e0; Delta[i] = 0.0e0; } Y[1] = 1.0e0; /* These are function and Jacobian evaluation counters. */ Iwork[LIW - 1] = 0; Iwork[LIW] = 0; /* The first call initializes internal data values, and the second * call gives a consistent value of YPRIME. Note the reversed * positions of DELTA, YPRIME. This usage computes a consistent value * of YPRIME in the case of Index 0 systems. */ for (i = 0; i <= 1; i++) { ires = i; ddassf7( &t, y, delta, yprime, rwork, &lda, &Rwork[1], &ires, rwork, iwork ); } tout = 0.01e0; ero = 0.0e0; nerr = 0; printf("\n Example Results for an Index-0 Banded ODE Problem, Soln in Forward Comm\n\n T Max Err BDF Order Last Step=H\n"); /* This shows passing data from evaluation routine ?DASF = ?DASF7 * to the calling program. */ for (iout = 1; iout <= 6; iout++) { ddaslx( ddassf7, NEQ, &t, y, yprime, tout, info, rtol, atol, &idid, rwork, LRW, iwork, LIW ); dedit2( y, t, &erm ); hu = Rwork[7]; nqu = Iwork[8]; printf(" %15.5e%14.3e%6ld%14.3e\n", t, erm, nqu, hu); if (idid < 0) goto L_50; er = erm/Atol[1]; ero = fmax( ero, er ); if (er > 1.0e0) { nerr += 1; } /* Advance to the next output point. */ tout *= 10.0e0; } L_50: ; /* Start over but solve banded linear system in reverse communication mod */ t = 0.0e0; for (i = 1; i <= NEQ; i++) { Y[i] = 0.0e0; Delta[i] = 0.0e0; } Y[1] = 1.0e0; /* The first call initializes internal data values, and the second * call gives a consistent value of YPRIME. Note the reversed * positions of DELTA, YPRIME. This computes a consistent value * of YPRIME in the case of Index 0 systems. */ for (i = 0; i <= 1; i++) { ires = i; ddassf7( &t, y, delta, yprime, rwork, &lda, &Rwork[1], &ires, rwork, iwork ); } tout = 0.01e0; ero = 0.0e0; nerr = 0; Info[1] = 0; printf("\n Example Results for an Index-0 Banded ODE Problem, Soln in Reverse Comm\n\n T Max Err BDF Order Last Step=H\n"); Iwork[LIW - 1] = 0; Iwork[LIW] = 0; /* All work done using reverse communication: */ Info[5] = -14; for (iout = 1; iout <= 6; iout++) { L_80: ddaslx( ddassf7, NEQ, &t, y, yprime, tout, 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_90; ires = Iwork[3]; /* Evaluate residuals or partials: */ if (ires <= 2) { ddassf7( &t, y, yprime, &Rwork[Iwork[4]], &Rwork[Iwork[5]], &lda, &Rwork[1], &ires, rwork, iwork ); } /* Factor the banded matrix: */ if (ires == 3) { /* The matrix is contained in the RWORK array. * With this request, IWORK(5) points to the matrix, and IWORK(3) == IRES * Note that IWORK(3) gets the error code (==0 means non-singular matrix) */ dgbfa( &Rwork[Iwork[5]], lda, NEQ, ML, MU, ipvt, &Iwork[3] ); } /* Solve with the banded matrix: */ if (ires == 4) { /* The right-hand side is contained in the RWORK array. * With this request, IWORK(4) points to the right-hand side. */ dgbsl( &Rwork[Iwork[5]], lda, NEQ, ML, MU, ipvt, &Rwork[Iwork[4]], 0 ); } goto L_80; L_90: ; dedit2( y, t, &erm ); hu = Rwork[7]; nqu = Iwork[8]; printf(" %15.5e%14.3e%6ld%14.3e\n", t, erm, nqu, hu); if (idid < 0) goto L_110; er = erm/Atol[1]; ero = fmax( ero, er ); if (er > 1.0e0) { nerr += 1; } /* Advance to the next output point. */ tout *= 10.0e0; } L_110: ; exit(0); } /* end of function */ void /*FUNCTION*/ ddassf7( double *t, double y[], double yprime[], double delta[], double *pd, long *ldp, double *cj, long *ires, double rwork[], long iwork[]) { #define PD(I_,J_) (*(pd+(I_)*(*ldp)+(J_))) long int _d_l, _d_m, _do0, _do1, i, j, k, mband; static long int ml, mu, ng; double d; static double alph1, alph2; /* 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 */ /* Example from SLATEC distribution. * This exercises the banded matrix solver. */ /* The last two value of IWORK() pass evaluation counters * back to the main program. */ /* This is the setup call. */ if (*ires == 0) { alph1 = 1.0e0; alph2 = 1.0e0; ng = 5; ml = 5; mu = 0; } /* The system residual value. */ if (*ires == 1) { for (j = 1; j <= ng; j++) { for (i = 1; i <= ng; i++) { k = i + (j - 1)*ng; d = -2.0e0*Y[k]; if (i != 1) d += Y[k - 1]*alph1; if (j != 1) d += Y[k - ng]*alph2; Delta[k] = d - Yprime[k]; } } /* Count residual evaluations using end of integer work array. */ Iwork[LIW - 1] += 1; } /* The partial derivative matrix. */ if (*ires == 2) { mband = ml + mu + 1; for (j = 1; j <= NEQ; j++) { PD(j - 1,mband - 1) = -2.0e0 - *cj; PD(j - 1,mband) = alph1; PD(j - 1,mband + 1) = 0.0e0; PD(j - 1,mband + 2) = 0.0e0; PD(j - 1,mband + 3) = 0.0e0; PD(j - 1,mband + 4) = alph2; } for (j = 1, _do0=DOCNT(j,NEQ,_do1 = ng); _do0 > 0; j += _do1, _do0--) { PD(j - 1,mband) = 0.0e0; } /* Count partial evaluations using end of integer work array. */ Iwork[LIW] += 1; } return; #undef PD } /* end of function */ void /*FUNCTION*/ dedit2( double y[], double t, double *erm) { long int _l0, i, j, k; double a1, a2, big, er, ex, ri, rj, yt; static double alph1 = 1.0e0; static double alph2 = 1.0e0; static long ng = 5; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /****BEGIN PROLOGUE DEDIT2 ****SUBSIDIARY ****LIBRARY SLATEC (DASSL) ****TYPE DOUBLE PRECISION (EDIT2-S, DEDIT2-D) ****AUTHOR PETZOLD, LINDA R., (LLNL) ****ROUTINES CALLED (NONE) ****END PROLOGUE DEDIT2 */ /****FIRST EXECUTABLE STATEMENT DEDIT2 */ *erm = 0.0e0; if (t == 0.0e0) return; ex = 0.0e0; big = .5e0*log( DBL_MAX ); /* Compute matrix exponential * initial data vector, * for this particular lower triangular matrix. */ if (t < big) ex = exp( -2.0e0*t ); a2 = 1.0e0; for (j = 1; j <= ng; j++) { a1 = 1.0e0; for (i = 1; i <= ng; i++) { k = i + (j - 1)*ng; yt = powi(t,i + j - 2)*ex*a1*a2; er = fabs( Y[k] - yt ); *erm = fmax( *erm, er ); ri = i; a1 = a1*alph1/ri; } rj = j; a2 = a2*alph2/rj; } return; } /* end of function */