/*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_ddasl1 s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include #include #include "p_ddasl1.h" #include "drddasl1.h" /* PARAMETER translations */ #define LDC (2*NEQ) #define LIW (30 + NEQ) #define LRW (45 + (5 + 2*MAXCON + 4)*NEQ + SQ(NEQ)) #define LTD NEQ #define MAXCON 2 #define NDIG 8 #define NEQ 2 #define NTIMES 10 #define TOL (powi(10.e0,-NDIG)) /* end of PARAMETER translations */ /* COMMON translations */ struct t_counts { long int kount0, kount1; } counts; /* end of COMMON translations */ int main( ) { long int _l0, i, idid, info[16], iwork[LIW]; double atol[NEQ], c[LTD][LDC], ftol, rnktol, rtol[NEQ], rwork[LRW], soln0[NEQ][NTIMES], soln1[NEQ][NTIMES], solp0[NEQ][NTIMES], solp1[NEQ][NTIMES], t, tout, y[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 Y = &y[0] - 1; double *const Yprime = &yprime[0] - 1; /* end of OFFSET VECTORS */ /*>> 2009-10-21 DRDDASL1 Hanson/Krogh Fixed initialization. *>> 2008-10-26 DRDDASL1 Krogh Moved Format statements up for C conv. *>> 2008-10-24 DRDDASL1 Krogh Removed in INCLUDE statement & cDEC$... *>> 2008-09-02 DRDDASL1 Hanson added starting computation of y' *>> 2008-08-26 DRDDASL1 Hanson added row dimensions to evaluators *>> 2006-04-10 DRDDASL1 Krogh Removed declaration of unused E. *>> 2002-05-29 DRDDASL1 Krogh Changes for C conversion problem. *>> 2001-10-11 DRDDASL1 R. J. Hanson Document Example Code, */ /* Solve a C. W. Gear index=2 problem. * Reduce it to an index=1, solve it. * Reduce it further to index=0, solve it. * Compare results, which are equivalent but * not exactly equal. */ /*--D replaces "?": DR?DASL1, ?DASLX, ?DASLS, ?DASSF1, ?DASSF0 */ /*++S Default NDIG = 4 *++ Default NDIG = 8 *++ Substitute for NDIG below */ /* Set number of equations: */ /* Set maximum number of constraints. */ /* Work space sizes: */ /* Tolerances: */ for (i = 1; i <= NEQ; i++) { Atol[i] = TOL; Rtol[i] = TOL; } /* Setup options: */ for (i = 1; i <= 16; i++) { Info[i] = 0; } /* Use partial derivatives provided in evaluator routines: */ Info[5] = 2; /* Constrain solution with 1 constraint: */ Info[10] = 1; /* Compute the initial value of YPRIME(*): */ t = 0; ftol = TOL; rnktol = TOL; /* Assign initial y, and guess for y', then get initial y'. */ Y[1] = 1.e0; Y[2] = 0.e0; Yprime[1] = 0.e0; Yprime[2] = 0.e0; ddasls( ddassf1, NEQ, &t, y, yprime, info, ftol, rnktol, (double*)c, ADR(_l0,LDC), LTD, &idid, rwork, LRW, iwork, LIW ); printf(" Example Results for a Transformed Index-2 DAE Problem\n\n T Y1/Y2 Y1P/Y2P\n"); printf(" Initial conditions for y,y' at t=0, index 1 system\n"); printf("%15.6e%15.6e%15.6e\n %15.6e%15.6e\n \n", t, Y[1], Yprime[1], Y[2], Yprime[2]); /* Allow up to 5000 steps: */ Info[12] = 5000; for (i = 1; i <= NTIMES; i++) { /* Integrate from T=I-1 to TOUT=T+1. Final TOUT=10. */ t = i - 1; tout = t + 1; ddaslx( ddassf1, NEQ, &t, y, yprime, tout, info, rtol, atol, &idid, rwork, LRW, iwork, LIW ); printf("%15.6e%15.6e%15.6e\n %15.6e%15.6e\n \n", tout, Y[1], Yprime[1], Y[2], Yprime[2]); /* Save solution and derivative for comparison to index 0 values. */ soln1[0][i - 1] = Y[1]; soln1[1][i - 1] = Y[2]; solp1[0][i - 1] = Yprime[1]; solp1[1][i - 1] = Yprime[2]; } /* Start integration again. */ Info[1] = 0; for (i = 1; i <= NEQ; i++) { Atol[i] = TOL; Rtol[i] = TOL; } /* Switch from 1 to 2 constraints, and use the index 0 system. */ Info[10] = 2; printf(" Differences, (Index-1) - (Index-0) Values\n\n T Y1/Y2 Y1P/Y2P\n"); t = 0.e0; /* Assign initial y, and guess for y', then get initial y'. */ Y[1] = 1.e0; Y[2] = 0.e0; Yprime[1] = 0.e0; Yprime[2] = 0.e0; ddasls( ddassf0, NEQ, &t, y, yprime, info, ftol, rnktol, (double*)c, ADR(_l0,LDC), LTD, &idid, rwork, LRW, iwork, LIW ); printf(" Initial conditions for y,y' at t=0, index 0 system\n"); printf("%15.6e%15.6e%15.6e\n %15.6e%15.6e\n \n", t, Y[1], Yprime[1], Y[2], Yprime[2]); for (i = 1; i <= NTIMES; i++) { /* Integrate from T=I-1 to TOUT=T+1. Final TOUT=10. */ t = i - 1; tout = t + 1; ddaslx( ddassf0, NEQ, &t, y, yprime, tout, info, rtol, atol, &idid, rwork, LRW, iwork, LIW ); /* Use solution and derivative differences for comparison * with index 1 values. */ soln0[0][i - 1] = soln1[0][i - 1] - Y[1]; soln0[1][i - 1] = soln1[1][i - 1] - Y[2]; solp0[0][i - 1] = solp1[0][i - 1] - Yprime[1]; solp0[1][i - 1] = solp1[1][i - 1] - Yprime[2]; printf("%15.6e%15.6e%15.6e\n %15.6e%15.6e\n \n", tout, soln0[0][i - 1], solp0[0][i - 1], soln0[1][i - 1], solp0[1][i - 1]); } /* Print number of residual evaluations for index 1 and index 0 * problems. */ printf(" Index-1 and Index-0 residual evalutions:%5ld%5ld\n", counts.kount1, counts.kount0); exit(0); } /* end of function */ void /*FUNCTION*/ ddassf1( double *t, double y[], double yprime[], double delta[], double *d1, long *ldd, double *cj, long *ires, double rwork[], long iwork[]) { #define D1(I_,J_) (*(d1+(I_)*(*ldd)+(J_))) double eta, one, two, zero; /* 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 Gear index=2 problem. * One equation is differentiated to reduce it to index 1, * with a constraint on the index 2 equation. */ one = 1.e0; two = 2.e0; zero = 0.e0; eta = 10.e0; /* This is the setup call. */ if (*ires == 0) { counts.kount1 = 0; } /* The system residual value. */ if (*ires == 1) { Delta[1] = Yprime[1] + eta**t*Yprime[2] + (one + eta)*Y[2] - sin( *t ); /* This second equation comes from differentiating the second equation in * section C.1, and subtracting the result from the first equation. */ Delta[2] = Y[2] - two*sin( *t ); /* Count function evaluations. */ counts.kount1 += 1; } /* The partial of the iteration matrix with respect to y. This is an * index 1 system. d1 is set to 0 prior to all calls here. Partials * are based on equations above. Note that \partial y'_i / y_i is c_j. */ if (*ires == 2) { D1(0,0) = *cj; D1(1,0) = (one + eta) + *cj*eta**t; D1(1,1) = one; } /* The constraining equation after the corrector has converged. * Both partials and residuals are required. This is for the second * equation in C.1. */ if (*ires == 5) { D1(0,2) = one; D1(1,2) = eta**t; Delta[3] = Y[1] + eta**t*Y[2] - cos( *t ); } /* The values of IRES=6,7 and 8 occur for the starting procedure * that solves for y'. First the partials with respect to t of what * is computedb when ires is 1. Here and below we are computing * partials of f not of the iteration matrix. */ if (*ires == 6) { Delta[1] = eta*Yprime[2] - cos( *t ); Delta[2] = -two*cos( *t ); } /* Compute the partial w.r.t y' of the equations defined when IRES=1 */ if (*ires == 7) { D1(0,0) = one; D1(1,0) = eta**t; } /* Compute the partial w.r.t y of the equations defined when IRES=1 */ if (*ires == 8) { D1(1,0) = one + eta; D1(1,1) = one; } return; #undef D1 } /* end of function */ void /*FUNCTION*/ ddassf0( double *t, double y[], double yprime[], double delta[], double *d0, long *ldd, double *cj, long *ires, double rwork[], long iwork[]) { #define D0(I_,J_) (*(d0+(I_)*(*ldd)+(J_))) double eta, one, two, zero; /* 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 Gear index=2 problem. * One equation is differentiated twice to reduce it to index 0. * This gives constraints on the index 2 and index 1 equations. */ one = 1.e0; two = 2.e0; zero = 0.e0; eta = 10.e0; /* This is the setup call. */ if (*ires == 0) { counts.kount0 = 0; } /* The system residual value. */ if (*ires == 1) { Delta[1] = Yprime[1] + eta**t*Yprime[2] + (one + eta)*Y[2] - sin( *t ); Delta[2] = Yprime[2] - two*cos( *t ); /* Count function evaluations. */ counts.kount0 += 1; } /* The mixed partial derivative matrix. * This is an index 0 system. */ if (*ires == 2) { D0(0,0) = *cj; D0(1,0) = (one + eta) + *cj*eta**t; D0(1,1) = *cj; } /* The constraining equations after the corrector has converged. * Both partials and residuals are required. */ if (*ires == 5) { D0(0,2) = one; D0(1,2) = eta**t; D0(0,3) = zero; D0(1,3) = one; Delta[3] = Y[1] + eta**t*Y[2] - cos( *t ); Delta[4] = Y[2] - two*sin( *t ); } /* The partial w.r.t y' for the starting procedure * Since this is an index 0 system the cases IRES=6,8 will not occur */ if (*ires == 7) { D0(0,0) = one; D0(1,0) = eta**t; D0(1,1) = one; } return; #undef D0 } /* end of function */