/*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 "ddastp.h" #include /* PARAMETER translations */ #define IDB 6 #define IFIRST 16 #define ISMOOT 13 #define LCJ 1 #define LCJOLD 6 #define LCNSTR 10 #define LCTF 17 #define LETF 16 #define LHMIN 10 #define LHOLD 7 #define LIRES 3 #define LJCALC 19 #define LK 7 #define LKOLD 8 #define LMAT 9 #define LMXORD 6 #define LNJAC 8 #define LNJE 15 #define LNRE 14 #define LNS 11 #define LNST 13 #define LPHASE 20 #define LROUND 9 #define REVLOC 21 /* end of PARAMETER translations */ void /*FUNCTION*/ ddastp( 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[], double alpha[], double beta[], double gamma[], double psi[], double sigma[], long *k) { #define PHI(I_,J_) (*(phi+(I_)*(neq)+(J_))) static LOGICAL32 convgd; static long int i, ires, j, j1, kdiff, km1, knew, kp1, kp2, locate, m, ncf, nef, nsf, nsp1; static double alpha0, alphas, cjlast, ck, conrate, delnrm, enorm, erk, erkm1, erkm2, erkp1, err, est, estold, hnew, oldnrm, pnorm, r, rate, ratio, sc, temp1, temp2, terk, terkm1, terkm2, terkp1, xold; static long maxit = 4; static double xrate = 0.25e0; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Alpha = &alpha[0] - 1; double *const Beta = &beta[0] - 1; double *const Delta = &delta[0] - 1; double *const E = &e[0] - 1; double *const Gamma = &gamma[0] - 1; long *const Info = &info[0] - 1; long *const Iwork = &iwork[0] - 1; double *const Psi = &psi[0] - 1; double *const Rwork = &rwork[0] - 1; double *const Sigma = &sigma[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. *>> 2008-10-24 ddastp Krogh Declared dnrm2 *>> 2008-08-26 ddastp Hanson add argument of leading dimension to ddasf *>> 2006-05-18 ddastp Hanson Install test for inconsistent constraints *>> 2006-04-14 ddastp Krogh Zero high differences on order increase. *>> 2003-03-05 ddastp Hanson Install Soderlind stepsize code. *>> 2002-06-26 ddastp Krogh Insured iwork(lk) has current value. *>> 2001-12-12 ddastp Krogh Changed code for reverse communication *>> 2001-11-23 ddastp Krogh Changed many names per library conventions. *>> 2001-11-04 ddastp Krogh Fixes for F77 and conversion to single *>> 2001-11-01 ddastp Hanson Provide code to Math a la Carte. *--D replaces "?": ?dastp, ?dasj, ?dasco, ?dasin, ?dasf, *-- & ?dasnm, ?daslv, ?daslx, ?dasdb, ?dasgh, ?nrm2, ?copy */ /****BEGIN PROLOGUE DDASTP ****SUBSIDIARY ****PURPOSE Perform one step of the DDASLX integration. ****LIBRARY SLATEC (DDASLX) ****TYPE DOUBLE PRECISION (SDASTP-S, DDASTP-D) ****AUTHOR Petzold, Linda R., (LLNL) ****DESCRIPTION *----------------------------------------------------------------------- * DDASTP SOLVES A SYSTEM OF DIFFERENTIAL/ * ALGEBRAIC EQUATIONS OF THE FORM * G(X,Y,YPRIME) = 0, FOR ONE STEP (NORMALLY * FROM X TO X+H). * * THE METHODS USED ARE MODIFIED DIVIDED * DIFFERENCE,FIXED LEADING COEFFICIENT * FORMS OF BACKWARD DIFFERENTIATION * FORMULAS. THE CODE ADJUSTS THE STEPSIZE * AND ORDER TO CONTROL THE LOCAL ERROR PER * STEP. * * * THE PARAMETERS REPRESENT * X -- INDEPENDENT VARIABLE * Y -- SOLUTION VECTOR AT X * YPRIME -- DERIVATIVE OF SOLUTION VECTOR AFTER SUCCESSFUL STEP * NEQ -- NUMBER OF EQUATIONS TO BE INTEGRATED * DDASF -- EXTERNAL USER-SUPPLIED SUBROUTINE * TO EVALUATE THE RESIDUAL. THE CALL IS * CALL DDASF(X,Y,YPRIME,DELTA,D,LDD,C,IRES,RWORK,IWORK) * X,Y,YPRIME ARE INPUT. DELTA IS OUTPUT. * ON INPUT, IRES=0. DDASF SHOULD ALTER IRES ONLY * IF IT ENCOUNTERS AN ILLEGAL VALUE OF Y OR A * STOP CONDITION. SET IRES=-1 IF AN INPUT VALUE * OF Y IS ILLEGAL, AND DDASTP WILL TRY TO SOLVE * THE PROBLEM WITHOUT GETTING IRES = -1. IF * IRES=-2, DDASTP RETURNS CONTROL TO THE CALLING * PROGRAM WITH IDID = -1. IRES CAN ALSO BE SET TO * LARGE NEGATIVE VALUES TO SET DEBUGGING PRINT. * H -- APPROPRIATE STEP SIZE FOR NEXT STEP. * NORMALLY DETERMINED BY THE CODE * WT -- VECTOR OF WEIGHTS FOR ERROR CRITERION. * IDID -- COMPLETION CODE WITH THE FOLLOWING MEANINGS: * IDID= 1 -- THE STEP WAS COMPLETED SUCCESSFULLY * IDID=-1 -- IRES EQUAL TO -2 WAS ENCOUNTERED, * AND CONTROL IS BEING RETURNED TO * THE CALLING PROGRAM * IDID=-4 -- THE CORRECTOR COULD NOT CONVERGE * BECAUSE IRES WAS EQUAL TO MINUS ONE * IDID=-6 -- THE ERROR TEST FAILED REPEATEDLY * IDID=-7 -- THE CORRECTOR COULD NOT CONVERGE * IDID=-8 -- THE ITERATION MATRIX IS SINGULAR * IDID=-9 -- THE CORRECTOR COULD NOT CONVERGE. * THERE WERE REPEATED ERROR TEST * FAILURES ON THIS STEP. * PHI -- ARRAY OF DIVIDED DIFFERENCES USED BY * DDASTP. THE LENGTH IS NEQ*(K+1),WHERE * K IS THE MAXIMUM ORDER * DELTA,E -- WORK VECTORS FOR DDASTP OF LENGTH NEQ * WM,IWORK -- REAL AND INTEGER ARRAYS STORING * MATRIX INFORMATION SUCH AS THE MATRIX * OF PARTIAL DERIVATIVES,PERMUTATION * VECTOR, AND VARIOUS OTHER INFORMATION. * RWORK -- THE USUAL WORK ARRAY. * ALPHA, BETA, GAMMA, PSI, SIGMA -- USED FOR INTEGRATION * COEFFICIENTS. * K -- THE CURRENT INTEGRATION ORDER. * * THE OTHER PARAMETERS ARE INFORMATION * WHICH IS NEEDED INTERNALLY BY DDASTP TO * CONTINUE FROM STEP TO STEP. * *----------------------------------------------------------------------- ****ROUTINES CALLED DDASJ, DDASNM, DDASLV, DDASIN ****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) * 981119 Replaced RES, JAC by DDASF, RJH. ****END PROLOGUE DDASTP * */ /* DDASGH -- computes the new stepsize. There are three controllers * one can choose: H211b(b=4), PI.4.2 and standard * control. * CTRNM -- the name of controller: * CTRNM = H211B4: H221B(B=4) controller <= only one used * = PI42: PI.4.2 controller * = STAND: standard controller * default controller is H211b4 <= only one used */ /* POINTERS INTO IWORK */ /* POINTERS INTO RWORK */ /* POINTERS INTO INFO */ /*----------------------------------------------------------------------- * BLOCK 1. * INITIALIZE. ON THE FIRST CALL,SET * THE ORDER TO 1 AND INITIALIZE * OTHER VARIABLES. *----------------------------------------------------------------------- * * INITIALIZATIONS FOR ALL CALLS ****FIRST EXECUTABLE STATEMENT DDASTP */ locate = Iwork[REVLOC]%8; if (locate > 0) { ires = Iwork[LIRES]; Iwork[REVLOC] /= 8; switch (locate) { case 1: goto L_130; case 2: goto L_140; case 3: goto L_145; case 4: goto L_190; case 5: goto L_230; case 6: goto L_260; case 7: goto L_286; } } /* No reverse communication active */ if (locate < 0) { /* FIRST STEP INITIALIZATIONS */ Iwork[LETF] = 0; Iwork[LCTF] = 0; Rwork[LHOLD] = 0.0e0; Psi[1] = *h; Rwork[LCJOLD] = 1.0e0/ *h; Rwork[LCJ] = Rwork[LCJOLD]; Rwork[LNJAC] = 100.e0; Iwork[LJCALC] = -1; delnrm = 1.0e0; Iwork[LPHASE] = 0; Iwork[LNS] = 0; Iwork[REVLOC] = 0; ratio = 1.e0; } /* Not doing reverse communication */ *idid = 1; xold = *x; ncf = 0; nsf = 0; nef = 0; /*----------------------------------------------------------------------- * BLOCK 2 * COMPUTE COEFFICIENTS OF FORMULAS FOR * THIS STEP. *----------------------------------------------------------------------- */ L_20: ; *k = Iwork[LK]; kp1 = *k + 1; kp2 = *k + 2; km1 = *k - 1; xold = *x; if ((*h != Rwork[LHOLD]) || (*k != Iwork[LKOLD])) Iwork[LNS] = 0; Iwork[LNS] = min( Iwork[LNS] + 1, Iwork[LKOLD] + 2 ); nsp1 = Iwork[LNS] + 1; if (kp1 < Iwork[LNS]) goto L_40; Beta[1] = 1.0e0; Alpha[1] = 1.0e0; temp1 = *h; Gamma[1] = 0.0e0; Sigma[1] = 1.0e0; for (i = 2; i <= kp1; i++) { temp2 = Psi[i - 1]; Psi[i - 1] = temp1; Beta[i] = Beta[i - 1]*Psi[i - 1]/temp2; temp1 = temp2 + *h; Alpha[i] = *h/temp1; Sigma[i] = (i - 1)*Sigma[i - 1]*Alpha[i]; Gamma[i] = Gamma[i - 1] + Alpha[i - 1]/ *h; } Psi[kp1] = temp1; L_40: ; /* COMPUTE ALPHAS, ALPHA0 */ alphas = 0.0e0; alpha0 = 0.0e0; for (i = 1; i <= *k; i++) { alphas += -1.0e0/i; alpha0 -= Alpha[i]; } /* COMPUTE LEADING COEFFICIENT RWORK(LCJ) */ cjlast = Rwork[LCJ]; Rwork[LCJ] = -alphas/ *h; /* COMPUTE VARIABLE STEPSIZE ERROR COEFFICIENT CK */ ck = fabs( Alpha[kp1] + alphas - alpha0 ); ck = fmax( ck, Alpha[kp1] ); /* DECIDE WHETHER NEW JACOBIAN IS NEEDED */ temp1 = (1.0e0 - xrate)/(1.0e0 + xrate); temp2 = 1.0e0/temp1; if ((Rwork[LCJ]/Rwork[LCJOLD] < temp1) || (Rwork[LCJ]/Rwork[LCJOLD] > temp2)) Iwork[LJCALC] = -1; if (Rwork[LCJ] != cjlast) Rwork[LNJAC] = 100.e0; /* CHANGE PHI TO PHI STAR */ if (kp1 < nsp1) goto L_80; for (j = nsp1; j <= kp1; j++) { for (i = 1; i <= neq; i++) { PHI(j - 1,i - 1) *= Beta[j]; } } L_80: ; /* UPDATE TIME */ *x += *h; /*----------------------------------------------------------------------- * BLOCK 3 * PREDICT THE SOLUTION AND DERIVATIVE, * AND SOLVE THE CORRECTOR EQUATION *----------------------------------------------------------------------- * * stepping past TOUT. Y(:)} is obtained by interpolation. * YPRIME(:) is obtained by interpolation. * 4 The integration has paused for reverse communication. Respond * based on the values of IRES * Task Interupted * -1 IRES set to -2 by the user. * -2 Accuracy requested exceeds machine precision. RTOL and ATOL * have been increased. * -3 There have been too many steps between output points. * Quit or Restart Integration * * FIRST,PREDICT THE SOLUTION AND DERIVATIVE */ L_90: ; for (i = 1; i <= neq; i++) { Y[i] = PHI(0,i - 1); Yprime[i] = 0.0e0; } for (j = 2; j <= kp1; j++) { for (i = 1; i <= neq; i++) { Y[i] += PHI(j - 1,i - 1); Yprime[i] += Gamma[j]*PHI(j - 1,i - 1); } } pnorm = ddasnm( neq, y, wt, rwork, iwork ); /* SOLVE THE CORRECTOR EQUATION USING A * MODIFIED NEWTON SCHEME. */ convgd = TRUE; m = 0; 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[LIRES] = ires; /* Need to put locations [from E(*)] to DELTA(*). */ Iwork[REVLOC] = 8*Iwork[REVLOC] + 1; return; } /* REVERSE ENTRY 1: */ L_130: ; if (Info[IDB] != 0) ddasdb( 3, neq, *x, y, yprime, info, rwork, iwork, ires, rwork, rwork ); if (ires < 0) { if (ires >= -2) goto L_290; Info[IDB] = -ires; ires = 0; } /* IF INDICATED,REEVALUATE THE * ITERATION MATRIX PD = DG/DY + RWORK(LCJ)*DG/DYPRIME * (WHERE G(X,Y,YPRIME)=0). SET * IWORK(LJCALC) TO 0 AS AN INDICATOR THAT * THIS HAS BEEN DONE. */ if (Iwork[LJCALC] != -1) goto L_150; Iwork[LNJE] += 1; Iwork[LJCALC] = 0; /* REVERSE ENTRY 2: */ L_140: ; ddasj( neq, ldd, x, y, yprime, delta, *h, wt, e, wm, iwork, rwork, ddasf, info, &ires ); if (Iwork[REVLOC] != 0) { if (Iwork[REVLOC] < 0) { Iwork[REVLOC] = 3; } else { Iwork[REVLOC] = 8*Iwork[REVLOC] + 2; } return; } /* REVERSE ENTRY 3: */ L_145: ; Rwork[LCJOLD] = Rwork[LCJ]; Rwork[LNJAC] = 100.e0; if (ires < 0) { if (ires >= -2) goto L_290; Info[IDB] = -ires; ires = 0; } if (ires != 0) goto L_290; nsf = 0; /* INITIALIZE THE ERROR ACCUMULATION VECTOR E. */ L_150: ; for (i = 1; i <= neq; i++) { E[i] = 0.0e0; } /* CORRECTOR LOOP. */ L_170: ; /* MULTIPLY RESIDUAL BY TEMP1 TO ACCELERATE CONVERGENCE */ temp1 = 2.0e0/(1.0e0 + Rwork[LCJ]/Rwork[LCJOLD]); for (i = 1; i <= neq; i++) { Delta[i] *= temp1; } /* 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] = 4; return; } /* REVERSE ENTRY 4: */ L_190: ; if (Info[IDB] != 0) ddasdb( 4, neq, *x, y, yprime, info, rwork, iwork, ires, wt, rwork ); /* UPDATE Y, E, AND YPRIME */ for (i = 1; i <= neq; i++) { Y[i] -= Delta[i]; E[i] -= Delta[i]; Yprime[i] += -Rwork[LCJ]*Delta[i]; } /* TEST FOR CONVERGENCE OF THE ITERATION */ delnrm = ddasnm( neq, delta, wt, rwork, iwork ); /* Set convergence rate allowed as Soderlind proposes: */ conrate = 0.125e0; if (Info[ISMOOT] != 0) conrate = .33e0; if (delnrm <= 100.e0*Rwork[LROUND]*pnorm) goto L_250; if (m > 0) goto L_210; oldnrm = delnrm; goto L_220; L_210: rate = pow(delnrm/oldnrm,1.0e0/m); if (rate > 0.90e0) goto L_240; Rwork[LNJAC] = rate/(1.0e0 - rate); L_220: ; if ((Rwork[LNJAC]*delnrm) <= conrate) goto L_250; /* THE CORRECTOR HAS NOT YET CONVERGED. * UPDATE M AND TEST WHETHER THE * MAXIMUM NUMBER OF ITERATIONS HAVE * BEEN TRIED. */ m += 1; if (m >= maxit) goto L_240; /* EVALUATE THE RESIDUAL * AND GO BACK TO DO ANOTHER ITERATION */ 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[LIRES] = ires; Iwork[REVLOC] = 5; return; } /* REVERSE ENTRY 5: */ L_230: ; if (Info[IDB] != 0) ddasdb( 3, neq, *x, y, yprime, info, rwork, iwork, ires, rwork, rwork ); if (ires < 0) { if (ires >= -2) goto L_290; Info[IDB] = -ires; ires = 0; } goto L_170; /* THE CORRECTOR FAILED TO CONVERGE IN MAXIT * ITERATIONS. IF THE ITERATION MATRIX * IS NOT CURRENT,RE-DO THE STEP WITH * A NEW ITERATION MATRIX. */ L_240: ; if (Iwork[LJCALC] == 0) goto L_290; Iwork[LJCALC] = -1; goto L_90; /* THE ITERATION HAS CONVERGED. IF CONSTRAINTS on SOLUTION are * REQUIRED, SET THE SOLUTION, IF THE PERTURBATION * TO DO IT IS SMALL ENOUGH. IF THE CHANGE IS TOO LARGE, THEN * CONSIDER THE CORRECTOR ITERATION TO HAVE FAILED. * Also have an option right here to project back onto constraints. * This is primarily for problems that started with index > 1. */ L_250: if (Iwork[LCNSTR] == 0) goto L_300; /* Define weights to use for projecting back to constraints. */ for (i = 1; i <= neq; i++) { /* DELTA(I)=ONE/WT(I) */ Delta[i] = Wt[i]; /* DELTA(I)=sqrt(WT(I)) */ } ires = 5; 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[LIRES] = ires; /* Need to put constraint change in DELTA(*). */ Iwork[REVLOC] = 6; return; } /* REVERSE ENTRY 6: */ L_260: ; if (Info[IDB] != 0) ddasdb( 3, neq, *x, y, yprime, info, rwork, iwork, ires, rwork, rwork ); if (ires < 0) { if (ires >= -2) goto L_290; Info[IDB] = -ires; ires = 0; } /* The user has computed info for constraints. Compute a single weighted * least-distance Newton step back to the constraints. Put move into * delta(:). */ if (labs( Iwork[LMAT] ) <= 4) ddasco( &Wm[neq + 1], *ldd, neq, wt, delta ); /* If ldd .eq. neq there are no constraints. Nothing needs * to be done. */ if (*ldd == neq) goto L_300; delnrm = ddasnm( neq, delta, wt, rwork, iwork ); sc = ddasnm( neq, e, wt, rwork, iwork ); if (delnrm <= sc) { for (i = 1; i <= neq; i++) { E[i] -= Delta[i]; } } else { for (i = 1; i <= neq; i++) { /* Scale change so it is no larger that the current size of e(:). */ E[i] += -(sc/delnrm)*Delta[i]; } } /* Check if the residual norm for the linear system was positive. * If so then compare this value with what a perturbation in y implies * for a perturbation in the constraints. */ if (Delta[neq + 1] == 0) goto L_300; /* Save a copy of the solution */ dcopy( neq, y, 1, delta, 1 ); /* Perturb the solution values in units of the largest amount of * error allowed on each component. A random number could be used * here instead of 1, -1, 1, .... */ sc = 1; for (i = 1; i <= neq; i++) { Y[i] += sc/Wt[i]; sc = -sc; } ires = 5; sc = Delta[neq + 1]; if (Iwork[LMAT] > 0) { (*ddasf)( x, y, yprime, delta, wm, ldd, &Rwork[LCJ], &ires, rwork, iwork ); } else { Iwork[LIRES] = ires; Iwork[REVLOC] = 7; return; } /* REVERSE ENTRY 7: */ L_286: ; /* Compute the residual norm on the constraints with the * perturbation. */ delnrm = dnrm2( *ldd - neq, &Delta[neq + 1], 1 ); /* If the solution residual norm is almost as large * as the perturbation norm then there is a serious problem. */ if (sc >= 0.5*delnrm) { /* This is the error flag assigned to this condition. */ *idid = -30; /* These values are needed to respond to the error condition. */ Delta[1] = sc; Delta[2] = delnrm; return; } /* Restore the solution if the system appears consistent. */ dcopy( neq, delta, 1, y, 1 ); goto L_300; /* EXITS FROM BLOCK 3 * NO CONVERGENCE WITH CURRENT ITERATION * MATRIX,OR SINGULAR ITERATION MATRIX */ L_290: convgd = FALSE; if (ires == -2) { *idid = -1; return; } L_300: Iwork[LJCALC] = 1; if (!convgd) goto L_490; /*----------------------------------------------------------------------- * BLOCK 4 * ESTIMATE THE ERRORS AT ORDERS K,K-1,K-2 * AS IF CONSTANT STEPSIZE WAS USED. ESTIMATE * THE LOCAL ERROR AT ORDER K AND TEST * stepping past TOUT. Y(:)} is obtained by interpolation. * YPRIME(:) is obtained by interpolation. * 4 The integration has paused for reverse communication. Respond * based on the values of IRES. * Task Interupted * -1 IRES set to -2 by the user. * -2 Accuracy requested exceeds machine precision. RTOL and ATOL * have been increased. * -3 There have been too many steps between output points. * Quit or Restart Integration * * WHETHER THE CURRENT STEP IS SUCCESSFUL. *----------------------------------------------------------------------- * * ESTIMATE ERRORS AT ORDERS K,K-1,K-2 */ if (Info[IFIRST] != 0) estold = est; enorm = ddasnm( neq, e, wt, rwork, iwork ); erk = Sigma[*k + 1]*enorm; terk = (*k + 1)*erk; est = erk; knew = *k; if (*k == 1) goto L_350; /* To check if order should be decreased */ for (i = 1; i <= neq; i++) { Delta[i] = PHI(kp1 - 1,i - 1) + E[i]; } erkm1 = Sigma[*k]*ddasnm( neq, delta, wt, rwork, iwork ); terkm1 = *k*erkm1; if (*k > 2) goto L_320; if (terkm1 <= 0.5e0*terk) goto L_340; goto L_350; L_320: ; for (i = 1; i <= neq; i++) { Delta[i] += PHI(*k - 1,i - 1); } erkm2 = Sigma[*k - 1]*ddasnm( neq, delta, wt, rwork, iwork ); terkm2 = (*k - 1)*erkm2; if (fmax( terkm1, terkm2 ) > terk) goto L_350; /* LOWER THE ORDER */ L_340: ; knew = *k - 1; est = erkm1; /* CALCULATE THE LOCAL ERROR FOR THE CURRENT STEP * TO SEE IF THE STEP WAS SUCCESSFUL */ L_350: ; err = ck*enorm; if (err > 1.0e0) goto L_490; /*----------------------------------------------------------------------- * BLOCK 5 * THE STEP IS SUCCESSFUL. DETERMINE * THE BEST ORDER AND STEPSIZE FOR * THE NEXT STEP. UPDATE THE DIFFERENCES * FOR THE NEXT STEP. *----------------------------------------------------------------------- */ *idid = 1; Iwork[LNST] += 1; kdiff = *k - Iwork[LKOLD]; Iwork[LKOLD] = *k; Rwork[LHOLD] = *h; /* ESTIMATE THE ERROR AT ORDER K+1 UNLESS: * ALREADY DECIDED TO LOWER ORDER, OR * ALREADY USING MAXIMUM ORDER, OR * STEPSIZE NOT CONSTANT, OR * ORDER RAISED IN PREVIOUS STEP */ if ((knew == km1) || (*k == Iwork[LMXORD])) Iwork[LPHASE] = 1; if (Iwork[LPHASE] == 0) goto L_400; if (knew == km1) goto L_390; if (*k == Iwork[LMXORD]) goto L_410; /* Free the order and let the error sequence decide if the order * gets increased or otherwise changed. */ if (Info[ISMOOT] != 0 && ((kp1 >= Iwork[LNS]) || (kdiff == 1))) goto L_410; for (i = 1; i <= neq; i++) { /* stepping past TOUT. Y(:)} is obtained by interpolation. * YPRIME(:) is obtained by interpolation. * 4 The integration has paused for reverse communication. Respond * based on the values of IRES. * Task Interupted * -1 IRES set to -2 by the user. * -2 Accuracy requested exceeds machine precision. RTOL and ATOL * have been increased. * -3 There have been too many steps between output points. * Quit or Restart Integration * */ Delta[i] = E[i] - PHI(kp2 - 1,i - 1); } erkp1 = (1.0e0/(*k + 2))*ddasnm( neq, delta, wt, rwork, iwork ); terkp1 = (*k + 2)*erkp1; if (*k > 1) goto L_370; if (terkp1 >= 0.5e0*terk) goto L_410; goto L_380; L_370: if (terkm1 <= fmin( terk, terkp1 )) goto L_390; if (terkp1 >= terk || *k == Iwork[LMXORD]) goto L_410; /* RAISE ORDER */ L_380: *k = kp1; est = erkp1; if (*k < Iwork[LMXORD]) { for (i = 1; i <= neq; i++) { PHI(*k + 1,i - 1) = 0.e0; } } goto L_410; /* LOWER ORDER */ L_390: *k = km1; est = erkm1; goto L_410; /* IF IWORK(LPHASE) = 0, INCREASE ORDER BY 1 AND MULTIPLY STEPSIZE BY * FACTOR 2 */ L_400: *k = kp1; if (Info[IFIRST] == 0) estold = est; /* Selectively use Soderlind's smoothing: */ if (Info[ISMOOT] != 0) { hnew = *h*2.0e0; *h = hnew; goto L_440; } /* DETERMINE THE APPROPRIATE STEPSIZE FOR * THE NEXT STEP. */ L_410: ; if (Info[ISMOOT] != 0) { hnew = *h; temp2 = *k + 1; r = pow(2.0e0*est + 0.0001e0,-1.0e0/temp2); if (r < 2.0e0) goto L_420; hnew = 2.0e0**h; goto L_430; L_420: if (r > 1.0e0) goto L_430; r = fmax( 0.5e0, fmin( 0.9e0, r ) ); hnew = *h*r; L_430: *h = hnew; goto L_440; } else { /* New Soderlind logic goes here: */ temp2 = *k + 1; ddasgh( est, estold, &ratio, *h, &hnew, temp2 ); *h = hnew; } /* UPDATE DIFFERENCES FOR NEXT STEP */ L_440: ; if (Iwork[LKOLD] == Iwork[LMXORD]) goto L_460; for (i = 1; i <= neq; i++) { PHI(kp2 - 1,i - 1) = E[i]; } L_460: ; for (i = 1; i <= neq; i++) { PHI(kp1 - 1,i - 1) += E[i]; } for (j1 = 2; j1 <= kp1; j1++) { j = kp1 - j1 + 1; for (i = 1; i <= neq; i++) { PHI(j - 1,i - 1) += PHI(j,i - 1); } } Iwork[LK] = *k; Info[IFIRST] = 1; return; /*----------------------------------------------------------------------- * BLOCK 6 * THE STEP IS UNSUCCESSFUL. RESTORE X,PSI,PHI * DETERMINE APPROPRIATE STEPSIZE FOR * CONTINUING THE INTEGRATION, OR EXIT WITH * AN ERROR FLAG IF THERE HAVE BEEN MANY * FAILURES. *----------------------------------------------------------------------- */ L_490: Iwork[LPHASE] = 1; /* RESTORE X,PHI,PSI */ *x = xold; if (kp1 < nsp1) goto L_520; for (j = nsp1; j <= kp1; j++) { temp1 = 1.0e0/Beta[j]; for (i = 1; i <= neq; i++) { PHI(j - 1,i - 1) *= temp1; } } L_520: ; for (i = 2; i <= kp1; i++) { Psi[i - 1] = Psi[i] - *h; } /* TEST WHETHER FAILURE IS DUE TO CORRECTOR ITERATION * OR ERROR TEST */ if (convgd) goto L_560; Iwork[LCTF] += 1; /* THE NEWTON ITERATION FAILED TO CONVERGE WITH * A CURRENT ITERATION MATRIX. DETERMINE THE CAUSE * OF THE FAILURE AND TAKE APPROPRIATE ACTION. */ if (ires == 0) goto L_540; /* THE ITERATION MATRIX IS SINGULAR. REDUCE * THE STEPSIZE BY A FACTOR OF 4. IF * THIS HAPPENS THREE TIMES IN A ROW ON * THE SAME STEP, RETURN WITH AN ERROR FLAG */ nsf += 1; r = 0.25e0; *h *= r; if ((nsf < 3) && (fabs( *h ) >= Rwork[LHMIN])) goto L_600; *idid = -8; goto L_590; /* THE NEWTON ITERATION FAILED TO CONVERGE FOR A REASON * OTHER THAN A SINGULAR ITERATION MATRIX. IF IRES = -2, THEN * RETURN. OTHERWISE, REDUCE THE STEPSIZE AND TRY AGAIN, UNLESS * TOO MANY FAILURES HAVE OCCURRED. */ L_540: ; if (ires == -2) { *idid = -1; } else { ncf += 1; r = 0.25e0; *h *= r; if ((ncf < 10) && (fabs( *h ) >= Rwork[LHMIN])) goto L_600; *idid = -7; if (ires == -1) *idid = -4; if (nef >= 3) *idid = -9; } goto L_590; /* THE NEWTON SCHEME CONVERGED, AND THE CAUSE * OF THE FAILURE WAS THE ERROR ESTIMATE * EXCEEDING THE TOLERANCE. */ L_560: nef += 1; Iwork[LETF] += 1; if (nef > 1) goto L_570; /* ON FIRST ERROR TEST FAILURE, KEEP CURRENT ORDER OR LOWER * ORDER BY ONE. COMPUTE NEW STEPSIZE BASED ON DIFFERENCES * OF THE SOLUTION. */ *k = knew; temp2 = *k + 1; if (Info[ISMOOT] != 0) { r = 0.90e0*pow(2.0e0*est + 0.0001e0,-1.0e0/temp2); r = fmax( 0.25e0, fmin( 0.9e0, r ) ); *h *= r; } else { ddasgh( est, estold, &ratio, *h, &hnew, temp2 ); *h = hnew; } if (fabs( *h ) >= Rwork[LHMIN]) goto L_600; *idid = -6; goto L_590; /* ON SECOND ERROR TEST FAILURE, USE THE CURRENT ORDER OR * DECREASE ORDER BY ONE. REDUCE THE STEPSIZE BY A FACTOR OF * FOUR. */ L_570: if (nef > 2) goto L_580; Iwork[LK] = knew; *h *= 0.25e0; if (fabs( *h ) >= Rwork[LHMIN]) goto L_600; *idid = -6; goto L_590; /* ON THIRD AND SUBSEQUENT ERROR TEST FAILURES, SET THE ORDER TO * ONE AND REDUCE THE STEPSIZE BY A FACTOR OF FOUR. */ L_580: Iwork[LK] = 1; *h *= 0.25e0; if (fabs( *h ) >= Rwork[LHMIN]) goto L_600; *idid = -6; goto L_590; /* FOR ALL CRASHES, RESTORE Y TO ITS LAST VALUE, * INTERPOLATE TO FIND YPRIME AT LAST X, AND RETURN */ L_590: ; ddasin( *x, *x, y, yprime, neq, *k, phi, psi ); *k = *k; return; /* GO BACK AND TRY THIS STEP AGAIN */ L_600: ; goto L_20; /*------END OF SUBROUTINE DDASTP------ */ #undef PHI } /* end of function */