/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:55 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "ddas1.h" #include /* PARAMETER translations */ #define IDB 6 #define LCJ 1 #define LCNSTR 10 #define LHMIN 10 #define LIRES 3 #define LMAT 9 #define LMXORD 6 #define LNJE 15 #define LNRE 14 #define LROUND 9 #define REVLOC 21 /* end of PARAMETER translations */ void /*FUNCTION*/ ddas1( double *x, double y[], double yprime[], long neq, long *ldd, void (*ddasf)(double*,double[],double[],double[],double[],long*,double*,long*,double[],long[]), long info[], double *h, double wt[], long *idid, double *phi, double delta[], double e[], double wm[], long iwork[], double rwork[]) { #define PHI(I_,J_) (*(phi+(I_)*(neq)+(J_))) static LOGICAL32 convgd; static long int i, ires, j, jcalc, locate, m, ncf, nef, nsf; static double delnrm, err, oldnrm, r, rate, s, xold, ynorm; static long maxit = 10; static long mjac = 5; static double damp = 0.75e0; /* 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. *>> 2009-09-30 ddas1 Krogh Removed unused variable sc. *>> 2008-08-26 ddas1 Hanson add argument of leading dimension to ddasf *>> 2006-05-18 ddas1 Hanson added checks for inconsistent constraints *>> 2006-04-16 ddas1 Krogh declared j. *>> 2004-04-13 ddas1 Hanson cleared phi array at outset. *>> 2003-03-06 ddas1 Hanson changed weights to use reciprocals *>> 2001-12-12 ddas1 Krogh Changed code for reverse communication *>> 2001-11-23 ddas1 Krogh Changed many names per library conventions. *>> 2001-11-01 ddas1 Hanson Provide code to Math a la Carte. *--D replaces "?": ?das1, ?dasf, ?dasj, ?daslv, ?dasco, *--& ?daslx, ?dasnm, ?dasdb ****BEGIN PROLOGUE DDAS1 ****SUBSIDIARY ****PURPOSE Initialization routine for DDASLX. ****LIBRARY SLATEC (DDASLX) ****TYPE DOUBLE PRECISION (SDAS1-S, DDAS1-D) ****AUTHOR Petzold, Linda R., (LLNL) ****DESCRIPTION *----------------------------------------------------------------- * DDAS1 TAKES ONE STEP OF SIZE H OR SMALLER WITH THE BACKWARD EULER * METHOD, TO FIND YPRIME. X AND Y ARE UPDATED TO BE CONSISTENT WITH * THE NEW STEP. A MODIFIED DAMPED NEWTON ITERATION IS USED TO SOLVE * THE CORRECTOR ITERATION. * * THE INITIAL GUESS FOR YPRIME IS USED IN THE PREDICTION, AND IN * FORMING THE ITERATION MATRIX, BUT IS NOT INVOLVED IN THE ERROR * TEST. THIS MAY HAVE TROUBLE CONVERGING IF THE INITIAL GUESS IS NO * GOOD, OR IF G(X,Y,YPRIME) DEPENDS NONLINEARLY ON YPRIME. * * THE PARAMETERS REPRESENT: * X -- INDEPENDENT VARIABLE * Y -- SOLUTION VECTOR AT X * YPRIME -- DERIVATIVE OF SOLUTION VECTOR * NEQ -- NUMBER OF EQUATIONS * LDD -- FIRST (ROW) DIMENSION FOR MATIX. * DDASF -- THE USER DEFINED ROUTINE. * INFO -- THE USER SUPPLIED ARRAY DEFINING OPTIONS. * H -- STEPSIZE. IMDER MAY USE A STEPSIZE SMALLER THAN H. * WT -- VECTOR OF WEIGHTS FOR ERROR CRITERION * IDID -- COMPLETION CODE WITH THE FOLLOWING MEANINGS * IDID= 1 -- YPRIME WAS FOUND SUCCESSFULLY * IDID=-10 -- DDAS1 FAILED TO FIND YPRIME * IDID=-28 -- IRES NOT RESET WHEN IRES=5 WITH USER DEFINED MATRIX. * PHI -- WORK SPACE FOR DDAS1 * DELTA,E -- WORK SPACE FOR DDAS1 * WM,IWORK -- REAL AND INTEGER ARRAYS STORING * MATRIX INFORMATION * RWORK -- THE USUAL WORK ARRAY. * *----------------------------------------------------------------- ****ROUTINES CALLED DDASF, DDASNM, DDASLV ****REVISION HISTORY (YYMMDD) * 830315 DATE WRITTEN * 901009 Finished conversion to SLATEC 4.0 format (F.N.Fritsch) * 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) * 901030 Minor corrections to declarations. (FNF) * 981118 Use one external user routine, RJH. ****END PROLOGUE DDAS1 * */ /* POINTERS INTO IWORK */ /* POINTERS INTO RWORK */ /* POINTERS INTO INFO */ /*--------------------------------------------------- * BLOCK 1. * INITIALIZATIONS. *--------------------------------------------------- * ****FIRST EXECUTABLE STATEMENT DDAS1 */ if (Iwork[REVLOC] != 0) { ires = Iwork[LIRES]; locate = Iwork[REVLOC]%8; Iwork[REVLOC] /= 8; /* 1 2 3 4 5 6 */ switch (locate) { case 1: goto L_50; case 2: goto L_70; case 3: goto L_100; case 4: goto L_320; case 5: goto L_75; case 6: goto L_360; } } /* Control drops through here on first call: */ *idid = 1; nef = 0; ncf = 0; nsf = 0; xold = *x; ynorm = ddasnm( neq, y, wt, rwork, iwork ); /* SAVE Y AND YPRIME IN PHI */ for (i = 1; i <= neq; i++) { PHI(0,i - 1) = Y[i]; PHI(1,i - 1) = Yprime[i]; } /* Clear rest of phi(*,*) array to define differences. */ for (j = 3; j <= (Iwork[LMXORD] + 1); j++) { for (i = 1; i <= neq; i++) { PHI(j - 1,i - 1) = 0.e0; } } /*---------------------------------------------------- * BLOCK 2. * DO ONE BACKWARD EULER STEP. *---------------------------------------------------- * * SET UP FOR START OF CORRECTOR ITERATION */ L_20: Rwork[LCJ] = 1.0e0/ *h; *x += *h; /* PREDICT SOLUTION AND DERIVATIVE */ for (i = 1; i <= neq; i++) { Y[i] += *h*Yprime[i]; } jcalc = -1; m = 0; convgd = TRUE; /* CORRECTOR LOOP. */ L_40: Iwork[LNRE] += 1; ires = 1; if (Info[IDB] != 0) ddasdb( 2, neq, *x, y, yprime, info, rwork, iwork, ires, rwork, rwork ); if (Iwork[LMAT] >= 0) { (*ddasf)( x, y, yprime, delta, wm, ldd, &Rwork[LCJ], &ires, rwork, iwork ); } else { Iwork[REVLOC] = 8*Iwork[REVLOC] + 1; Iwork[LIRES] = ires; return; } /* REVERSE ENTRY 1: */ L_50: ; /* Test signal for invalid request. */ if (Info[IDB] != 0) ddasdb( 3, neq, *x, y, yprime, info, rwork, iwork, ires, rwork, rwork ); if (ires < 0) { if (ires == -1) goto L_190; if (ires == -2) goto L_240; Info[IDB] = -ires; ires = 0; } /* EVALUATE THE ITERATION MATRIX */ if (jcalc != -1) goto L_80; Iwork[LNJE] += 1; jcalc = 0; /* REVERSE ENTRY 2: */ L_70: ; ddasj( neq, ldd, x, y, yprime, delta, *h, wt, e, wm, iwork, rwork, ddasf, info, &ires ); /* See if reverse communication needed: */ if (Iwork[REVLOC] != 0) { if (Iwork[REVLOC] < 0) { Iwork[REVLOC] = 5; } else { Iwork[REVLOC] = 8*Iwork[REVLOC] + 2; } return; } /* REVERSE ENTRY 5: After user computes and factors the matrix. */ L_75: ; s = 1000000.e0; if (ires < 0) { if (ires == -1) goto L_190; if (ires == -2) goto L_240; Info[IDB] = -ires; ires = 0; } if (ires != 0) goto L_190; nsf = 0; /* MULTIPLY RESIDUAL BY DAMPING FACTOR */ L_80: ; for (i = 1; i <= neq; i++) { Delta[i] *= damp; } /* COMPUTE A NEW ITERATE (BACK SUBSTITUTION) * STORE THE CORRECTION IN DELTA * */ ddaslv( neq, ldd, x, y, yprime, delta, ddasf, info, iwork, rwork ); if (Iwork[REVLOC] != 0) { Iwork[REVLOC] = 3; return; } /* REVERSE ENTRY 3: */ L_100: ; /* UPDATE Y AND YPRIME */ for (i = 1; i <= neq; i++) { Y[i] -= Delta[i]; Yprime[i] += -Rwork[LCJ]*Delta[i]; } /* TEST FOR CONVERGENCE OF THE ITERATION. * */ delnrm = ddasnm( neq, delta, wt, rwork, iwork ); if (delnrm <= 100.e0*Rwork[LROUND]*ynorm) goto L_200; if (m > 0) goto L_120; oldnrm = delnrm; goto L_130; L_120: rate = pow(delnrm/oldnrm,1.0e0/m); if (rate > 0.90e0) goto L_190; s = rate/(1.0e0 - rate); L_130: if (s*delnrm <= 0.33e0) goto L_200; /* THE CORRECTOR HAS NOT YET CONVERGED. UPDATE * M AND AND TEST WHETHER THE MAXIMUM * NUMBER OF ITERATIONS HAVE BEEN TRIED. * EVERY MJAC ITERATIONS, GET A NEW * ITERATION MATRIX. * */ m += 1; if (m >= maxit) goto L_190; if ((m/mjac)*mjac == m) jcalc = -1; goto L_40; /* THE ITERATION HAS CONVERGED. */ /* EXITS FROM CORRECTOR LOOP. */ L_190: convgd = FALSE; goto L_220; L_200: ; /*----------------------------------------------------- * BLOCK 3. * THE CORRECTOR ITERATION CONVERGED. * DO ERROR TEST. *----------------------------------------------------- * */ for (i = 1; i <= neq; i++) { E[i] = Y[i] - PHI(0,i - 1); } err = ddasnm( neq, e, wt, rwork, iwork ); if (err <= 1.0e0) goto L_280; /*-------------------------------------------------------- * BLOCK 4. * THE BACKWARD EULER STEP FAILED. RESTORE X, Y * AND YPRIME TO THEIR ORIGINAL VALUES. * REDUCE STEPSIZE AND TRY AGAIN, IF * POSSIBLE. *--------------------------------------------------------- * */ L_220: ; *x = xold; for (i = 1; i <= neq; i++) { Y[i] = PHI(0,i - 1); Yprime[i] = PHI(1,i - 1); } if (convgd) goto L_260; if (ires != 0) { nsf += 1; *h *= 0.25e0; if ((nsf < 3) && (fabs( *h ) >= Rwork[LHMIN])) goto L_270; *idid = -10; goto L_280; } if (ires != -2) goto L_250; L_240: *idid = -10; goto L_280; L_250: ncf += 1; *h *= 0.25e0; if ((ncf < 10) && (fabs( *h ) >= Rwork[LHMIN])) goto L_270; *idid = -10; goto L_280; L_260: nef += 1; r = 0.90e0/(2.0e0*err + 0.0001e0); r = fmax( 0.1e0, fmin( 0.5e0, r ) ); *h *= r; if ((fabs( *h ) >= Rwork[LHMIN]) && (nef < 10)) goto L_270; *idid = -10; goto L_280; L_270: goto L_20; L_280: ; /* CHECK for CONSTRAINTS */ if (Iwork[LCNSTR] == 0) goto L_400; ires = 5; if (Iwork[LCNSTR] >= 0) { if (Info[IDB] != 0) ddasdb( 2, neq, *x, y, yprime, info, rwork, iwork, ires, rwork, rwork ); (*ddasf)( x, y, yprime, delta, wm, ldd, &Rwork[LCJ], &ires, rwork, iwork ); } else { Iwork[REVLOC] = 4; Iwork[LIRES] = ires; return; } /* REVERSE ENTRY 4: */ L_320: ; if (Info[IDB] != 0) ddasdb( 3, neq, *x, y, yprime, info, rwork, iwork, ires, rwork, rwork ); if (ires != 5) { if (ires < 0) { if (ires == -1) goto L_190; if (ires == -2) goto L_240; Info[IDB] = -ires; ires = 0; } } else { i = labs( Iwork[LMAT] ); if ((i == 5) || (i == 6)) { *idid = -28; goto L_400; } /* If IRES=5, compute a weighted least-distance Newton step * back to the constraints. Put move into delta(:). */ ddasco( &Wm[neq + 1], *ldd, neq, wt, delta ); } /* Check that the moves back to the constraints are smaller than * sizes of the error tolerances. If they are then apply the * moves. Otherwise evaluate the constraints and return an error flag. */ for (i = 1; i <= neq; i++) { /* Recall that tol_i = 1/ wt_i. */ if (fabs( Delta[i]*Wt[i] ) > 1) goto L_350; } /* If moving onto the constraints was reasonable, make the move. */ for (i = 1; i <= neq; i++) { Y[i] -= Delta[i]; Yprime[i] += -Rwork[LCJ]*Delta[i]; } goto L_400; L_350: ; /* Here the projection back to the constraints was * .gt. the requested tolerance in some component. * Evaluate the constraints at the current values. */ if (Iwork[LMAT] >= 0) { (*ddasf)( x, y, yprime, delta, wm, ldd, &Rwork[LCJ], &ires, rwork, iwork ); } else { Iwork[REVLOC] = 6; Iwork[LIRES] = ires; return; } /* REVERSE ENTRY 6: */ L_360: ; /* This error flag will result in printing a message * and residuals on the constraints. */ *idid = -29; /*-------------END OF SUBROUTINE DDAS1---------------------- */ L_400: return; #undef PHI } /* end of function */