/*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_ddasl3 s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include #include #include #include "p_ddasl3.h" #include "drddasl3.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 0 #define NDIG 9 #define NEQ 5 #define TOL (powi(10.e0,-NDIG)) /* end of PARAMETER translations */ /* COMMON translations */ struct t_counts { long int kr, kf, ks; } counts; /* end of COMMON translations */ int main( ) { LOGICAL32 constraint; long int _l0, _n, i, idid, info[16], iwork[LIW]; double atol[NEQ], c[LTD][LDC], drl, drv, ftol, length, rnktol, rtol[NEQ], rwork[LRW], 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-28 DRDDASL3 Remove implicit none. *>> 2009-10-21 DRDDASL3 Hanson/Krogh Fixed initialization. *>> 2008-09-04 DRDDASL3 Hanson added starting computation of y' *>> 2008-08-26 DRDDASL3 Hanson added row dimensions to evaluators *>> 2001-10-11 DRDDASL3 R. J. Hanson Example 3 Code, with Download */ /* Solve a planar pendulum problem in rectangular coordinates. * The equation is transformed from "Index 3" to "Index 1" * by differentiating the length constraint twice. But * this results in a (non-physical) drift away from the constraint * for the length. This example shows this drift at t=100. The prob * is fixed in Example DRDDASL5 by providing information that * enables the integrator to move back onto the constraints. */ /*--D replaces "?": DR?DASL3, ?DASLS, ?DASLX, ?DASSF3, ?COPY, ?ROTG * Version 6, Intel 10.1 */ /* Set number of equations: */ /* Set number of constraints. */ /* Work space sizes: */ /*++S Default NDIG = 5 *++ Default NDIG = 9 *++ 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 the initial value of YPRIME(*): * Base computation of the initial y' on the * report by Krogh, Hanson, (2008), "Solving * Constrained Differential-Algebraic Systems * Using Projections," (8/2008). */ t = 0; /* Tolerances for a small F and a rank tolerance: */ ftol = pow(DBL_EPSILON,2./3.); rnktol = ftol; /* Provide partials in dense format. */ Info[5] = 2; /* Assign initial y, and guess for y', then get initial y'. */ for (i = 1; i <= NEQ; i++) { Y[i] = 0.e0; Yprime[i] = 0.e0; } /* (Actual initial values for y(1) is set in ddassf3.) */ ddasls( ddassf3, NEQ, &t, y, yprime, info, ftol, rnktol, (double*)c, ADR(_l0,LDC), NEQ, &idid, rwork, LRW, iwork, LIW ); /* Allow more than nominal number of steps. */ Info[12] = 50000; length = 1.1e0; /* Use partial derivatives provided in DDASSF3, and do our own * linear algebra. */ Info[5] = 8; for (i = 1; i <= 10; i++) { /* Integrate from T=I-1 to TOUT=T+1. Final TOUT=10. */ t = i - 1; tout = t + 1; ddaslx( ddassf3, NEQ, &t, y, yprime, tout, info, rtol, atol, &idid, rwork, LRW, iwork, LIW ); /* Compute residuals on length and its variation. They * should be smaller than the tolerances used. Relatively * soon the constraints are not satisfied. */ drl = (sqrt( SQ(Y[1]) + SQ(Y[2]) ) - length)/(Rtol[1]*length + Atol[1]); drv = (Y[1]*Y[3] + Y[2]*Y[4])/(Rtol[1]*(fabs( Y[1]*Y[3] ) + fabs( Y[1]*Y[4] )) + Atol[1]); constraint = fabs( drl ) <= 1.e0 && fabs( drv ) <= 1.e0; if (!constraint) goto L_40; } L_40: ; printf(" Example Results for a Simple Pendulum Problem\n"); printf("\n T y_1 y_2 y_3 y_4 y_5\n%14.4e", tout); for(_n=0L; _n < sizeof(y)/sizeof(double); _n++) printf("%14.4e", y[_n]); printf("\n"); printf(" No. Residual Evaluations Factorizations No. Projection Solves-\n %6ld %6ld %6ld\n\n", counts.kr, counts.kf, counts.ks); printf("At the time value%6.2f the integration was stopped.\n", tout); if (!constraint) { printf("With no constraints, the length is expected to drift. The error in the\nlength and its derivative / (error tolerance) are: %8.2f%8.2f\n", drl, drv); } else { printf(" The pendulum length and its variation have small errors.\n"); } exit(0); } /* end of function */ /* PARAMETER translations */ #define GRAVITY 9.806650e0 #define LENGTH 1.1e0 #define MASS 1.e0 #define ONE 1.e0 #define TWO 2.e0 /* 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_))) static LOGICAL32 factor; double u, v; static double dc[4], ds[4], lsq, mg, tj, yt[5]; /* 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; double *const Yt = &yt[0] - 1; /* end of OFFSET VECTORS */ /* Routine for the swinging simple pendulum problem. * Conservation of length is not enforced. */ /* This is the setup call. */ if (*ires == 0) { mg = MASS*GRAVITY; lsq = SQ(LENGTH); Y[1] = LENGTH; counts.kr = 0; counts.kf = 0; counts.ks = 0; } /* The sytem residual value. */ if (*ires == 1) { Delta[1] = Y[3] - Yprime[1]; Delta[2] = Y[4] - Yprime[2]; Delta[3] = -Y[1]*Y[5] - MASS*Yprime[3]; Delta[4] = -Y[2]*Y[5] - MASS*Yprime[4] - mg; Delta[5] = MASS*(SQ(Y[3]) + SQ(Y[4])) - mg*Y[2] - lsq*Y[5]; /* Count residual evaluations. */ counts.kr += 1; } /* The partial derivative matrix. */ if (*ires == 2) { /* Save value of CJ and Y(*), and YPRIME(*), if needed. These are * local variables used in the solving step. */ tj = *cj; dcopy( 5, y, 1, yt, 1 ); } /* This matrix is factored with five plane rotations. * Rotations 1 and 2 are identical but in planes * (1,3) and (2,4). This example illustrates using * a special solver for the system. */ if (*ires == 3) { factor = TRUE; /* This passes the error code for singularity back to the integrator. * It is an issue when the user solves their own linear systems. */ *ires = 0; } /* Solve the corrector equation. */ if (*ires == 4) { if (factor) { /* Count number of factorization steps. */ counts.kf += 1; D(0,0) = -tj; D(0,2) = -Yt[5]; D(1,1) = -tj; D(1,3) = -Yt[5]; D(1,4) = -mg; D(2,0) = ONE; D(2,2) = -MASS*tj; D(2,4) = TWO*MASS*Yt[3]; D(3,1) = ONE; D(3,3) = D(2,2); D(3,4) = TWO*MASS*Yt[4]; D(4,2) = -Yt[1]; D(4,3) = -Yt[2]; D(4,4) = -lsq; drotg( &D(0,0), &D(0,2), &Dc[1], &Ds[1] ); D(2,0) = Dc[1] + Ds[1]*D(2,2); D(2,2) = -Ds[1] + Dc[1]*D(2,2); D(3,1) = D(2,0); D(3,3) = D(2,2); D(1,1) = D(0,0); D(4,0) = Ds[1]*D(4,2); D(4,2) *= Dc[1]; D(4,1) = Ds[1]*D(4,3); D(4,3) *= Dc[1]; D(0,0) = ONE/D(0,0); /* Eliminate now in planes 2 and 5, column 2. */ drotg( &D(1,1), &D(1,4), &Dc[2], &Ds[2] ); u = Dc[2]*D(2,1) + Ds[2]*D(2,4); D(2,4) = -Ds[2]*D(2,1) + Dc[2]*D(2,4); D(2,1) = u; v = Dc[2]*D(3,1) + Ds[2]*D(3,4); D(3,4) = -Ds[2]*D(3,1) + Dc[2]*D(3,4); D(3,1) = v; u = Dc[2]*D(4,1) + Ds[2]*D(4,4); D(4,4) = -Ds[2]*D(4,1) + Dc[2]*D(4,4); D(4,1) = u; D(1,1) = ONE/D(1,1); /* Eliminate in planes 3 and 5, column 3. */ drotg( &D(2,2), &D(2,4), &Dc[3], &Ds[3] ); u = Dc[3]*D(3,2) + Ds[3]*D(3,4); D(3,4) = -Ds[3]*D(3,2) + Dc[3]*D(3,4); D(3,2) = u; v = Dc[3]*D(4,2) + Ds[3]*D(4,4); D(4,4) = -Ds[3]*D(4,2) + Dc[3]*D(4,4); D(4,2) = v; D(2,2) = ONE/D(2,2); /* Eliminate in planes 4 and 5, column 4. */ drotg( &D(3,3), &D(3,4), &Dc[4], &Ds[4] ); u = Dc[4]*D(4,3) + Ds[4]*D(4,4); D(3,3) = ONE/D(3,3); D(4,4) = ONE/(-Ds[4]*D(4,3) + Dc[4]*D(4,4)); D(4,3) = u; /* It does not matter if FACTOR is TRUE or FALSE here. * Things work either way, but with less work when FACTOR = FALSE., * since the factorization step is done only when the integrator * requests it and not each solve step. */ factor = FALSE; } /* Count number of solve steps. */ counts.ks += 1; /* Apply rotations to right-hand side. */ u = Dc[1]*Delta[1] + Ds[1]*Delta[3]; v = Dc[1]*Delta[2] + Ds[1]*Delta[4]; Delta[3] = -Ds[1]*Delta[1] + Dc[1]*Delta[3]; Delta[4] = -Ds[1]*Delta[2] + Dc[1]*Delta[4]; Delta[1] = u; Delta[2] = v; u = Dc[2]*Delta[2] + Ds[2]*Delta[5]; Delta[5] = -Ds[2]*Delta[2] + Dc[2]*Delta[5]; Delta[2] = u; u = Dc[3]*Delta[3] + Ds[3]*Delta[5]; Delta[5] = -Ds[3]*Delta[3] + Dc[3]*Delta[5]; Delta[3] = u; u = Dc[4]*Delta[4] + Ds[4]*Delta[5]; Delta[5] = -Ds[4]*Delta[4] + Dc[4]*Delta[5]; Delta[4] = u; /* Back substitute: */ Delta[5] *= D(4,4); Delta[4] += -Delta[5]*D(4,3); Delta[3] += -Delta[5]*D(4,2); Delta[2] += -Delta[5]*D(4,1); Delta[1] += -Delta[5]*D(4,0); Delta[4] *= D(3,3); Delta[3] += -Delta[4]*D(3,2); Delta[2] += -Delta[4]*D(3,1); Delta[1] += -Delta[4]*D(3,0); Delta[3] *= D(2,2); Delta[2] += -Delta[3]*D(2,1); Delta[1] += -Delta[3]*D(2,0); Delta[2] *= D(1,1); Delta[1] += -Delta[2]*D(1,0); Delta[1] *= D(0,0); } if (*ires == 6) { /* This is the partial of F, the residual function, w.r.t. T. * There is no explicit time dependence so this is zero. * The contents of DELTA(:) are set zero by the calling program * so there is nothing to do. * DELTA(1:6)=ZERO */ } /* This is the partial of F w.r.t. y'. * Values not assigned are set zero by the calling program. */ if (*ires == 7) { D(0,0) = -ONE; D(1,1) = -ONE; D(2,2) = -MASS; D(3,3) = -MASS; counts.kf += 1; } /* This is the partial of F w.r.t. y. * Values not assigned are set zero by the calling program. */ if (*ires == 8) { D(0,2) = -Y[5]; D(1,3) = -Y[5]; D(1,4) = -mg; D(2,0) = ONE; D(2,4) = TWO*MASS*Y[3]; D(3,1) = ONE; D(3,4) = TWO*MASS*Y[4]; D(4,2) = -Y[1]; D(4,3) = -Y[2]; D(4,4) = -lsq; counts.kf += 1; } return; #undef D } /* end of function */