/*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_ddasl4 s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include <stdio.h>
#include <stdlib.h>
#include "p_ddasl4.h"
#include "drddasl4.h"
/* PARAMETER translations */
#define LIW (30 + NEQ)
#define LRW (45 + (5 + 2*MAXCON + 4)*NEQ)
#define MAXCON 0
#define NDIG 8
#define NEPOCH 10
#define NEQ (NSTATE + NVAREQ*NSTATE)
#define NSTATE 4
#define NVAREQ 5
#define TOL (powi(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;
double 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 */
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 Ival = &ival[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 */
/*>> 2008-10-29 DRDDASL4 Krogh Fixed print output for C.
*>> 2008-10-26 DRDDASL4 Krogh Changed Fortran 90 type statement to F77.
*>> 2008-10-26 DRDDASL4 Krogh Moved Format statements up for C conv.
*>> 2008-10-24 DRDDASL4 Krogh Removed in INCLUDE statement & cDEC$...
*>> 2008-09-17 DRDDASL4 Hanson added starting values for y'
*>> 2008-08-26 DRDDASL4 Hanson added row dimension to evaluator
*>> 2001-10-11 DRDDASL4 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. */
/*--D 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;
dcopy( NEQ, y, 0, y, 1 );
dcopy( 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:
ddaslx( ddassf4, 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;
dcopy( SQ(NSTATE), (double*)gy, 0, (double*)gy, 1 );
dcopy( NSTATE*NVAREQ, (double*)gp, 0, (double*)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: */
dcopy( NSTATE*NVAREQ, &Y[NSTATE + 1], 1, (double*)yp,
1 );
/* Extract (DY/DP)': */
dcopy( NSTATE*NVAREQ, &Yprime[NSTATE + 1], 1, (double*)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] += ddot( NSTATE, &gy[0][i - 1],
NSTATE, &yp[j - 1][0], 1 ) - ypp[j - 1][i - 1];
}
}
/* Move the results where the integrator uses them: */
dcopy( NSTATE, delta, 1, &Rwork[Iwork[4]], 1 );
dcopy( NSTATE*NVAREQ, (double*)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. */
dcopy( SQ(NSTATE), (double*)gy, 1, (double*)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. */
dgefa( (double*)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++)
{
dgesl( (double*)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(double); _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*/ ddassf4(
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_)))
/* 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 */
/* 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 */