/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:40 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include "ddasj.h"
#include <stdlib.h>
/* PARAMETER translations */
#define IDB 6
#define LCJ 1
#define LIPVT 31
#define LIRES 3
#define LMAT 9
#define LML 1
#define LMU 2
#define LNPD 18
#define LROUND 9
#define NTEMP 26
#define REVLOC 21
/* end of PARAMETER translations */
void /*FUNCTION*/ ddasj(
long neq,
long *ldd,
double *x,
double y[],
double yprime[],
double delta[],
double h,
double wt[],
double e[],
double wm[],
long iwork[],
double rwork[],
void (*ddasf)(double*,double[],double[],double[],double[],long*,double*,long*,double[],long[]),
long info[],
long *ires)
{
static long int _d_l, _d_m, _do0, _do1, _do2, _do3, i, i1, i2,
ier, ii, ipsave, isave, j, k, l, lmata, locate, mba, mband, meb1,
meband, msave, n, nrow;
static double del, delinv, squr, ypsave, ysave;
/* OFFSET Vectors w/subscript range: 1 to dimension */
double *const Delta = &delta[0] - 1;
double *const E = &e[0] - 1;
long *const Info = &info[0] - 1;
long *const Iwork = &iwork[0] - 1;
double *const Rwork = &rwork[0] - 1;
double *const Wm = &wm[0] - 1;
double *const Wt = &wt[0] - 1;
double *const Y = &y[0] - 1;
double *const Yprime = &yprime[0] - 1;
/* end of OFFSET VECTORS */
/* Copyright (c) 2006, Math a la Carte, Inc.
*>> 2010-08-26 ddasj Krogh Changed declaration of info to info(*).
*>> 2008-08-26 ddasj Hanson add argument of leading dimension to ddasf
*>> 2001-12-12 ddasj Krogh Changed code for reverse communication
*>> 2001-11-23 ddasj Krogh Changed many names per library conventions.
*>> 2001-11-04 ddasj Krogh Fixes for F77 and conversion to single & C
*>> 2001-11-01 ddasj Hanson Provide code to Math a la Carte.
*--D replaces "?": ?dasj, ?dasf, ?daslx, ?dasdb, ?swap, ?gbfa, ?gefa
****BEGIN PROLOGUE DDASJ
****SUBSIDIARY
****PURPOSE Compute the iteration matrix for DDASLX and form the
* LU-decomposition.
****LIBRARY SLATEC (DDASLX)
****AUTHOR Petzold, Linda R., (LLNL)
****DESCRIPTION
*----------------------------------------------------------------------- */
/* THIS ROUTINE COMPUTES THE ITERATION MATRIX PD=DG/DY+CJ*DG/DYPRIME
* (WHERE G(X,Y,YPRIME)=0, AND CJ IS CONTAINED IN RWORK(LCJ)).
* HERE PD IS COMPUTED BY THE USER-SUPPLIED ROUTINE DDASF IF |INFO(5)|
* IS 2, 4, 9, OR 13 AND IS COMPUTED BY FINITE DIFFERENCES IF |INFO(5)|
* IS 1, 3, 8, OR 12. IF INFO(5) < 0, THEN COMPUTATIONS OF J ARE DONE
* BY USING REVERSE COMMUNICTION IF THESE COMPUTATIONS ARE DONE BY THE
* THE USER. WHEN |INFO(5) = 5 OR 6, THE USER IS DOING ALL
* COMPUTATIONS ASSOCIATED WITH J AND THE ASSOCIATED LINEAR ALGEBRA.
* THE PARAMETERS HAVE THE FOLLOWING MEANINGS.
* NEQ = NUMBER OF EQUATIONS
* LDD = FIRST (ROW) DIMENSION OF THE MATRIX.
* X = CURRENT VALUE OF THE INDEPENDENT VARIABLE.
* Y = ARRAY CONTAINING PREDICTED VALUES
* YPRIME = ARRAY CONTAINING PREDICTED DERIVATIVES
* DELTA = RESIDUAL EVALUATED AT (X,Y,YPRIME)
* (USED ONLY FOR BAND MATRICES AND FOR REVERSE
* COMMUNICATION.)
* H = CURRENT STEPSIZE IN INTEGRATION
* WT = VECTOR OF WEIGHTS FOR COMPUTING NORMS
* E = WORK SPACE (TEMPORARY) OF LENGTH NEQ
* WM = REAL WORK SPACE FOR MATRICES. ON
* OUTPUT IT CONTAINS THE LU DECOMPOSITION
* OF THE ITERATION MATRIX.
* IWORK = INTEGER WORK SPACE CONTAINING MATRIX INFORMATION
* DDASF = NAME OF THE EXTERNAL USER-SUPPLIED ROUTINE
* TO EVALUATE THE RESIDUAL FUNCTION G(X,Y,YPRIME)
* IRES = FLAG WHICH IS EQUAL TO ZERO IF NO ILLEGAL VALUES
* IN DDASF, AND LESS THAN ZERO OTHERWISE. (IF IRES
* IS LESS THAN ZERO, THE MATRIX WAS NOT COMPLETED)
*-----------------------------------------------------------------------
****ROUTINES CALLED DGBFA, DGEFA
****REVISION HISTORY (YYMMDD)
* 830315 DATE WRITTEN
* 901009 Finished conversion to SLATEC 4.0 format (F.N.Fritsch)
* 901010 Modified three MAX calls to be all on one line. (FNF)
* 901019 Merged changes made by C. Ulrich with SLATEC 4.0 format.
* 901026 Added explicit declarations for all variables and minor
* cosmetic changes to prologue. (FNF)
* 901101 Corrected PURPOSE. (FNF)
****END PROLOGUE DDASJ
* */
/* POINTERS INTO IWORK */
/* POINTERS INTO RWORK */
/* POINTERS INTO INFO */
/****FIRST EXECUTABLE STATEMENT DDASJ */
if (Iwork[REVLOC] != 0)
{
*ires = Iwork[LIRES];
locate = Iwork[REVLOC]%8;
Iwork[REVLOC] /= 8;
switch (locate)
{
case 1: goto L_90;
case 2: goto L_50;
case 3: goto L_210;
case 4: goto L_170;
}
}
/* The first entry drops through here. */
lmata = labs( Iwork[LMAT] );
/* 1 2 3 4 5 6 7 8 9 10 11 12 13 14 */
switch (lmata)
{
case 1: goto L_30;
case 2: goto L_10;
case 3: goto L_140;
case 4: goto L_120;
case 5: goto L_250;
case 6: goto L_280;
case 7: goto L_30;
case 8: goto L_10;
case 9: goto L_140;
case 10: goto L_120;
case 11: goto L_30;
case 12: goto L_10;
case 13: goto L_140;
case 14: goto L_120;
}
goto L_30;
/* Dense full user-supplied matrix */
L_10:
;
for (i = 1; i <= Iwork[LNPD]; i++)
{
Wm[i] = 0.0e0;
}
*ires = 2;
if (Info[IDB] != 0)
ddasdb( 2, neq, *x, y, yprime, info, rwork, iwork, *ires,
rwork, rwork );
if (Iwork[LMAT] >= 0)
{
(*ddasf)( x, y, yprime, e, wm, ldd, &Rwork[LCJ], ires, rwork,
iwork );
goto L_90;
}
/* Reverse communication to compute the partials. */
Iwork[REVLOC] = 1;
Iwork[LIRES] = *ires;
/* This location will be the same as WM when the user responds. */
goto L_240;
/* DENSE FINITE-DIFFERENCE-GENERATED MATRIX */
L_30:
;
*ires = 1;
nrow = 0;
squr = sqrt( Rwork[LROUND] );
i = 0;
/* Loop over the columns to generate the approximate Jacobian. */
L_40:
i += 1;
if (i > neq)
goto L_80;
del = squr*fmax( fabs( Y[i] ), fmax( fabs( h*Yprime[i] ), fabs( Wt[i] ) ) );
del = sign( del, h*Yprime[i] );
del = (Y[i] + del) - Y[i];
ysave = Y[i];
ypsave = Yprime[i];
Y[i] += del;
Yprime[i] += Rwork[LCJ]*del;
if (Info[IDB] != 0)
ddasdb( 2, neq, *x, y, yprime, info, rwork, iwork, *ires,
rwork, rwork );
if (Iwork[LMAT] >= 0)
{
(*ddasf)( x, y, yprime, e, wm, ldd, &Rwork[LCJ], ires, rwork,
iwork );
goto L_60;
}
else
{
Iwork[REVLOC] = 2;
Iwork[LIRES] = *ires;
/* This is placed in the start of DELTA(*), in units of RWORK(*). */
dswap( neq, e, 1, delta, 1 );
goto L_240;
}
/* REVERSE ENTRY 2: */
L_50:
;
dswap( neq, e, 1, delta, 1 );
L_60:
;
if (Info[IDB] != 0)
ddasdb( 3, neq, *x, y, yprime, info, rwork, iwork, *ires,
rwork, rwork );
if (*ires < 0)
{
if (*ires >= -2)
goto L_240;
Info[IDB] = -*ires;
*ires = 1;
}
delinv = 1.0e0/del;
for (l = 1; l <= neq; l++)
{
Wm[nrow + l] = (E[l] - Delta[l])*delinv;
}
nrow += neq;
Y[i] = ysave;
Yprime[i] = ypsave;
goto L_40;
L_80:
;
/* DO DENSE-MATRIX LU DECOMPOSITION ON PD
* REVERSE ENTRY 1: */
L_90:
;
if (Info[IDB] != 0)
ddasdb( 3, neq, *x, y, yprime, info, rwork, iwork, *ires,
rwork, rwork );
if (*ires < 0)
{
if (*ires >= -2)
goto L_240;
Info[IDB] = -*ires;
*ires = 0;
}
if (lmata > 4)
{
*ires = 3;
if (lmata <= 10)
goto L_260;
goto L_290;
}
else
{
dgefa( wm, *ldd, neq, &Iwork[LIPVT], &ier );
if (ier == 0)
*ires = 0;
}
goto L_240;
/* BANDED USER-SUPPLIED MATRIX */
L_120:
for (i = 1; i <= Iwork[LNPD]; i++)
{
Wm[i] = 0.0e0;
}
*ires = 2;
if (Info[IDB] != 0)
ddasdb( 2, neq, *x, y, yprime, info, rwork, iwork, *ires,
rwork, rwork );
if (Iwork[LMAT] < 0)
{
Iwork[REVLOC] = 3;
Iwork[LIRES] = *ires;
goto L_240;
}
(*ddasf)( x, y, yprime, e, wm, ldd, &Rwork[LCJ], ires, rwork,
iwork );
meband = 2*Iwork[LML] + Iwork[LMU] + 1;
goto L_210;
/* BANDED FINITE-DIFFERENCE-GENERATED MATRIX */
L_140:
mband = Iwork[LML] + Iwork[LMU] + 1;
mba = min( mband, neq );
meband = mband + Iwork[LML];
meb1 = meband - 1;
msave = (neq/mband) + 1;
isave = Iwork[NTEMP] - 1;
ipsave = isave + msave;
*ires = 1;
squr = sqrt( Rwork[LROUND] );
j = 0;
L_150:
;
j += 1;
if (j > mba)
goto L_220;
for (n = j, _do0=DOCNT(n,neq,_do1 = mband); _do0 > 0; n += _do1, _do0--)
{
k = (n - j)/mband + 1;
Wm[isave + k] = Y[n];
Wm[ipsave + k] = Yprime[n];
del = squr*fmax( fabs( Y[n] ), fmax( fabs( h*Yprime[n] ),
fabs( Wt[n] ) ) );
del = sign( del, h*Yprime[n] );
del = (Y[n] + del) - Y[n];
Y[n] += del;
Yprime[n] += Rwork[LCJ]*del;
}
if (Info[IDB] != 0)
ddasdb( 2, neq, *x, y, yprime, info, rwork, iwork, *ires,
rwork, rwork );
if (Iwork[LMAT] >= 0)
{
(*ddasf)( x, y, yprime, e, wm, ldd, &Rwork[LCJ], ires, rwork,
iwork );
if (*ires < 0)
{
if (*ires >= -2)
goto L_240;
Info[IDB] = -*ires;
*ires = 0;
}
goto L_180;
}
else
{
Iwork[REVLOC] = 4;
*ires = 1;
Iwork[LIRES] = *ires;
/* This is placed in the start of DELTA(*), in units of RWORK(*). */
dswap( neq, e, 1, delta, 1 );
goto L_240;
}
/* REVERSE ENTRY 4: */
L_170:
;
dswap( neq, e, 1, delta, 1 );
L_180:
;
if (Info[IDB] != 0)
ddasdb( 3, neq, *x, y, yprime, info, rwork, iwork, *ires,
rwork, rwork );
for (n = j, _do2=DOCNT(n,neq,_do3 = mband); _do2 > 0; n += _do3, _do2--)
{
k = (n - j)/mband + 1;
Y[n] = Wm[isave + k];
Yprime[n] = Wm[ipsave + k];
del = squr*fmax( fabs( Y[n] ), fmax( fabs( h*Yprime[n] ),
fabs( Wt[n] ) ) );
del = sign( del, h*Yprime[n] );
del = (Y[n] + del) - Y[n];
delinv = 1.0e0/del;
i1 = max( 1, n - Iwork[LMU] );
i2 = min( neq, n + Iwork[LML] );
ii = n*meb1 - Iwork[LML];
for (i = i1; i <= i2; i++)
{
Wm[ii + i] = (E[i] - Delta[i])*delinv;
}
}
goto L_150;
/* DO LU DECOMPOSITION OF BANDED PD
* REVERSE ENTRY 3 */
L_210:
if (Info[IDB] != 0)
ddasdb( 3, neq, *x, y, yprime, info, rwork, iwork, *ires,
rwork, rwork );
L_220:
;
if (*ires < 0)
{
if (*ires >= -2)
goto L_240;
Info[IDB] = -*ires;
*ires = 0;
}
if (lmata > 4)
{
*ires = 3;
if (lmata <= 10)
goto L_260;
goto L_290;
}
else
{
dgbfa( wm, meband, neq, Iwork[LML], Iwork[LMU], &Iwork[LIPVT],
&ier );
if (ier == 0)
*ires = 0;
}
L_240:
;
return;
/* User defined matrix, Get Jacobian and factor in one step */
L_250:
;
*ires = 2;
/* Enters here if already have matrix and want to factor */
L_260:
;
if (Info[IDB] != 0)
ddasdb( 2, neq, *x, y, yprime, info, rwork, iwork, *ires,
rwork, rwork );
(*ddasf)( x, y, yprime, e, wm, ldd, &Rwork[LCJ], ires, rwork,
iwork );
if (Info[IDB] != 0)
ddasdb( 3, neq, *x, y, yprime, info, rwork, iwork, *ires,
rwork, rwork );
if (*ires < -2)
{
Info[IDB] = -*ires;
*ires = 0;
}
goto L_240;
/* User defined but with reverse communication */
L_280:
*ires = 2;
/* Enter here if matrix is computed. */
L_290:
;
Iwork[LIRES] = *ires;
Iwork[REVLOC] = -1;
goto L_240;
/*------ END OF SUBROUTINE DDASJ ------ */
} /* end of function */