/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:45 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "sdasj.h" #include /* 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*/ sdasj( long neq, long *ldd, float *x, float y[], float yprime[], float delta[], float h, float wt[], float e[], float wm[], long iwork[], float rwork[], void (*sdasf)(float*,float[],float[],float[],float[],long*,float*,long*,float[],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 float del, delinv, squr, ypsave, ysave; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Delta = &delta[0] - 1; float *const E = &e[0] - 1; long *const Info = &info[0] - 1; long *const Iwork = &iwork[0] - 1; float *const Rwork = &rwork[0] - 1; float *const Wm = &wm[0] - 1; float *const Wt = &wt[0] - 1; float *const Y = &y[0] - 1; float *const Yprime = &yprime[0] - 1; /* end of OFFSET VECTORS */ /* Copyright (c) 2006, Math a la Carte, Inc. *>> 2010-08-26 sdasj Krogh Changed declaration of info to info(*). *>> 2008-08-26 sdasj Hanson add argument of leading dimension to sdasf *>> 2001-12-12 sdasj Krogh Changed code for reverse communication *>> 2001-11-23 sdasj Krogh Changed many names per library conventions. *>> 2001-11-04 sdasj Krogh Fixes for F77 and conversion to single & C *>> 2001-11-01 sdasj Hanson Provide code to Math a la Carte. *--S replaces "?": ?dasj, ?dasf, ?daslx, ?dasdb, ?swap, ?gbfa, ?gefa ****BEGIN PROLOGUE SDASJ ****SUBSIDIARY ****PURPOSE Compute the iteration matrix for SDASLX and form the * LU-decomposition. ****LIBRARY SLATEC (SDASLX) ****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 SDASF 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 * SDASF = 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 SDASF, AND LESS THAN ZERO OTHERWISE. (IF IRES * IS LESS THAN ZERO, THE MATRIX WAS NOT COMPLETED) *----------------------------------------------------------------------- ****ROUTINES CALLED SGBFA, SGEFA ****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 SDASJ * */ /* POINTERS INTO IWORK */ /* POINTERS INTO RWORK */ /* POINTERS INTO INFO */ /****FIRST EXECUTABLE STATEMENT SDASJ */ 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) sdasdb( 2, neq, *x, y, yprime, info, rwork, iwork, *ires, rwork, rwork ); if (Iwork[LMAT] >= 0) { (*sdasf)( 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 = sqrtf( 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*fmaxf( fabsf( Y[i] ), fmaxf( fabsf( h*Yprime[i] ), fabsf( Wt[i] ) ) ); del = signf( 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) sdasdb( 2, neq, *x, y, yprime, info, rwork, iwork, *ires, rwork, rwork ); if (Iwork[LMAT] >= 0) { (*sdasf)( 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(*). */ sswap( neq, e, 1, delta, 1 ); goto L_240; } /* REVERSE ENTRY 2: */ L_50: ; sswap( neq, e, 1, delta, 1 ); L_60: ; if (Info[IDB] != 0) sdasdb( 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) sdasdb( 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 { sgefa( 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) sdasdb( 2, neq, *x, y, yprime, info, rwork, iwork, *ires, rwork, rwork ); if (Iwork[LMAT] < 0) { Iwork[REVLOC] = 3; Iwork[LIRES] = *ires; goto L_240; } (*sdasf)( 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 = sqrtf( 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*fmaxf( fabsf( Y[n] ), fmaxf( fabsf( h*Yprime[n] ), fabsf( Wt[n] ) ) ); del = signf( del, h*Yprime[n] ); del = (Y[n] + del) - Y[n]; Y[n] += del; Yprime[n] += Rwork[LCJ]*del; } if (Info[IDB] != 0) sdasdb( 2, neq, *x, y, yprime, info, rwork, iwork, *ires, rwork, rwork ); if (Iwork[LMAT] >= 0) { (*sdasf)( 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(*). */ sswap( neq, e, 1, delta, 1 ); goto L_240; } /* REVERSE ENTRY 4: */ L_170: ; sswap( neq, e, 1, delta, 1 ); L_180: ; if (Info[IDB] != 0) sdasdb( 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*fmaxf( fabsf( Y[n] ), fmaxf( fabsf( h*Yprime[n] ), fabsf( Wt[n] ) ) ); del = signf( 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) sdasdb( 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 { sgbfa( 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) sdasdb( 2, neq, *x, y, yprime, info, rwork, iwork, *ires, rwork, rwork ); (*sdasf)( x, y, yprime, e, wm, ldd, &Rwork[LCJ], ires, rwork, iwork ); if (Info[IDB] != 0) sdasdb( 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 SDASJ ------ */ } /* end of function */