ddastp.f

Go to the documentation of this file.

      SUBROUTINE ddastp (X, Y, YPRIME, NEQ, RES, JAC, H, WT, JSTART,
     +   IDID, RPAR, IPAR, PHI, DELTA, E, WM, IWM, ALPHA, BETA, GAMMA,
     +   PSI, SIGMA, CJ, CJOLD, HOLD, S, HMIN, UROUND, IPHASE, JCALC,
     +   K, KOLD, NS, NONNEG, NTEMP)
C***BEGIN PROLOGUE  DDASTP
C***SUBSIDIARY
C***PURPOSE  Perform one step of the DDASSL integration.
C***LIBRARY   SLATEC (DASSL)
C***TYPE      DOUBLE PRECISION (SDASTP-S, DDASTP-D)
C***AUTHOR  PETZOLD, LINDA R., (LLNL)
C***DESCRIPTION
C-----------------------------------------------------------------------
C     DDASTP SOLVES A SYSTEM OF DIFFERENTIAL/
C     ALGEBRAIC EQUATIONS OF THE FORM
C     G(X,Y,YPRIME) = 0,  FOR ONE STEP (NORMALLY
C     FROM X TO X+H).
C
C     THE METHODS USED ARE MODIFIED DIVIDED
C     DIFFERENCE,FIXED LEADING COEFFICIENT
C     FORMS OF BACKWARD DIFFERENTIATION
C     FORMULAS. THE CODE ADJUSTS THE STEPSIZE
C     AND ORDER TO CONTROL THE LOCAL ERROR PER
C     STEP.
C
C
C     THE PARAMETERS REPRESENT
C     X  --        INDEPENDENT VARIABLE
C     Y  --        SOLUTION VECTOR AT X
C     YPRIME --    DERIVATIVE OF SOLUTION VECTOR
C                  AFTER SUCCESSFUL STEP
C     NEQ --       NUMBER OF EQUATIONS TO BE INTEGRATED
C     RES --       EXTERNAL USER-SUPPLIED SUBROUTINE
C                  TO EVALUATE THE RESIDUAL.  THE CALL IS
C                  CALL RES(X,Y,YPRIME,DELTA,IRES,RPAR,IPAR)
C                  X,Y,YPRIME ARE INPUT.  DELTA IS OUTPUT.
C                  ON INPUT, IRES=0.  RES SHOULD ALTER IRES ONLY
C                  IF IT ENCOUNTERS AN ILLEGAL VALUE OF Y OR A
C                  STOP CONDITION.  SET IRES=-1 IF AN INPUT VALUE
C                  OF Y IS ILLEGAL, AND DDASTP WILL TRY TO SOLVE
C                  THE PROBLEM WITHOUT GETTING IRES = -1.  IF
C                  IRES=-2, DDASTP RETURNS CONTROL TO THE CALLING
C                  PROGRAM WITH IDID = -11.
C     JAC --       EXTERNAL USER-SUPPLIED ROUTINE TO EVALUATE
C                  THE ITERATION MATRIX (THIS IS OPTIONAL)
C                  THE CALL IS OF THE FORM
C                  CALL JAC(X,Y,YPRIME,PD,CJ,RPAR,IPAR)
C                  PD IS THE MATRIX OF PARTIAL DERIVATIVES,
C                  PD=DG/DY+CJ*DG/DYPRIME
C     H --         APPROPRIATE STEP SIZE FOR NEXT STEP.
C                  NORMALLY DETERMINED BY THE CODE
C     WT --        VECTOR OF WEIGHTS FOR ERROR CRITERION.
C     JSTART --    INTEGER VARIABLE SET 0 FOR
C                  FIRST STEP, 1 OTHERWISE.
C     IDID --      COMPLETION CODE WITH THE FOLLOWING MEANINGS:
C                  IDID= 1 -- THE STEP WAS COMPLETED SUCCESSFULLY
C                  IDID=-6 -- THE ERROR TEST FAILED REPEATEDLY
C                  IDID=-7 -- THE CORRECTOR COULD NOT CONVERGE
C                  IDID=-8 -- THE ITERATION MATRIX IS SINGULAR
C                  IDID=-9 -- THE CORRECTOR COULD NOT CONVERGE.
C                             THERE WERE REPEATED ERROR TEST
C                             FAILURES ON THIS STEP.
C                  IDID=-10-- THE CORRECTOR COULD NOT CONVERGE
C                             BECAUSE IRES WAS EQUAL TO MINUS ONE
C                  IDID=-11-- IRES EQUAL TO -2 WAS ENCOUNTERED,
C                             AND CONTROL IS BEING RETURNED TO
C                             THE CALLING PROGRAM
C     RPAR,IPAR -- REAL AND INTEGER PARAMETER ARRAYS THAT
C                  ARE USED FOR COMMUNICATION BETWEEN THE
C                  CALLING PROGRAM AND EXTERNAL USER ROUTINES
C                  THEY ARE NOT ALTERED BY DDASTP
C     PHI --       ARRAY OF DIVIDED DIFFERENCES USED BY
C                  DDASTP. THE LENGTH IS NEQ*(K+1),WHERE
C                  K IS THE MAXIMUM ORDER
C     DELTA,E --   WORK VECTORS FOR DDASTP OF LENGTH NEQ
C     WM,IWM --    REAL AND INTEGER ARRAYS STORING
C                  MATRIX INFORMATION SUCH AS THE MATRIX
C                  OF PARTIAL DERIVATIVES,PERMUTATION
C                  VECTOR,AND VARIOUS OTHER INFORMATION.
C
C     THE OTHER PARAMETERS ARE INFORMATION
C     WHICH IS NEEDED INTERNALLY BY DDASTP TO
C     CONTINUE FROM STEP TO STEP.
C
C-----------------------------------------------------------------------
C***ROUTINES CALLED  DDAJAC, DDANRM, DDASLV, DDATRP
C***REVISION HISTORY  (YYMMDD)
C   830315  DATE WRITTEN
C   901009  Finished conversion to SLATEC 4.0 format (F.N.Fritsch)
C   901019  Merged changes made by C. Ulrich with SLATEC 4.0 format.
C   901026  Added explicit declarations for all variables and minor
C           cosmetic changes to prologue.  (FNF)
C***END PROLOGUE  DDASTP
C
      INTEGER  NEQ, JSTART, IDID, IPAR(*), IWM(*), IPHASE, JCALC, K,
     *   KOLD, NS, NONNEG, NTEMP
      DOUBLE PRECISION
     *   x, y(*), yprime(*), h, wt(*), rpar(*), phi(neq,*), delta(*),
     *   e(*), wm(*), alpha(*), beta(*), gamma(*), psi(*), sigma(*), cj,
     *   cjold, hold, s, hmin, uround
      EXTERNAL  res, jac
C
      EXTERNAL  ddajac, ddanrm, ddaslv, ddatrp
      DOUBLE PRECISION  DDANRM
C
      INTEGER  I, IER, IRES, J, J1, KDIFF, KM1, KNEW, KP1, KP2, LCTF,
     *   letf, lmxord, lnje, lnre, lnst, m, maxit, ncf, nef, nsf, nsp1
      double precision
     *   alpha0, alphas, cjlast, ck, delnrm, enorm, erk, erkm1,
     *   erkm2, erkp1, err, est, hnew, oldnrm, pnorm, r, rate, temp1,
     *   temp2, terk, terkm1, terkm2, terkp1, xold, xrate
      LOGICAL  CONVGD
C
      PARAMETER (LMXORD=3)
      parameter(lnst=11)
      parameter(lnre=12)
      parameter(lnje=13)
      parameter(letf=14)
      parameter(lctf=15)
C
      DATA maxit/4/
      DATA xrate/0.25d0/
C
C
C
C
C
C-----------------------------------------------------------------------
C     BLOCK 1.
C     INITIALIZE. ON THE FIRST CALL,SET
C     THE ORDER TO 1 AND INITIALIZE
C     OTHER VARIABLES.
C-----------------------------------------------------------------------
C
C     INITIALIZATIONS FOR ALL CALLS
C***FIRST EXECUTABLE STATEMENT  DDASTP
      idid=1
      xold=x
      ncf=0
      nsf=0
      nef=0
      IF(jstart .NE. 0) GO TO 120
C
C     IF THIS IS THE FIRST STEP,PERFORM
C     OTHER INITIALIZATIONS
      iwm(letf) = 0
      iwm(lctf) = 0
      k=1
      kold=0
      hold=0.0d0
      jstart=1
      psi(1)=h
      cjold = 1.0d0/h
      cj = cjold
      s = 100.d0
      jcalc = -1
      delnrm=1.0d0
      iphase = 0
      ns=0
120   CONTINUE
C
C
C
C
C
C-----------------------------------------------------------------------
C     BLOCK 2
C     COMPUTE COEFFICIENTS OF FORMULAS FOR
C     THIS STEP.
C-----------------------------------------------------------------------
200   CONTINUE
      kp1=k+1
      kp2=k+2
      km1=k-1
      xold=x
      IF(h.NE.hold.OR.k .NE. kold) ns = 0
      ns=min(ns+1,kold+2)
      nsp1=ns+1
      IF(kp1 .LT. ns)GO TO 230
C
      beta(1)=1.0d0
      alpha(1)=1.0d0
      temp1=h
      gamma(1)=0.0d0
      sigma(1)=1.0d0
      DO 210 i=2,kp1
         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
210      CONTINUE
      psi(kp1)=temp1
230   CONTINUE
C
C     COMPUTE ALPHAS, ALPHA0
      alphas = 0.0d0
      alpha0 = 0.0d0
      DO 240 i = 1,k
        alphas = alphas - 1.0d0/i
        alpha0 = alpha0 - alpha(i)
240     CONTINUE
C
C     COMPUTE LEADING COEFFICIENT CJ
      cjlast = cj
      cj = -alphas/h
C
C     COMPUTE VARIABLE STEPSIZE ERROR COEFFICIENT CK
      ck = abs(alpha(kp1) + alphas - alpha0)
      ck = max(ck,alpha(kp1))
C
C     DECIDE WHETHER NEW JACOBIAN IS NEEDED
      temp1 = (1.0d0 - xrate)/(1.0d0 + xrate)
      temp2 = 1.0d0/temp1
      IF (cj/cjold .LT. temp1 .OR. cj/cjold .GT. temp2) jcalc = -1
      IF (cj .NE. cjlast) s = 100.d0
C
C     CHANGE PHI TO PHI STAR
      IF(kp1 .LT. nsp1) GO TO 280
      DO 270 j=nsp1,kp1
         DO 260 i=1,neq
260         phi(i,j)=beta(j)*phi(i,j)
270      CONTINUE
280   CONTINUE
C
C     UPDATE TIME
      x=x+h
C
C
C
C
C
C-----------------------------------------------------------------------
C     BLOCK 3
C     PREDICT THE SOLUTION AND DERIVATIVE,
C     AND SOLVE THE CORRECTOR EQUATION
C-----------------------------------------------------------------------
C
C     FIRST,PREDICT THE SOLUTION AND DERIVATIVE
300   CONTINUE
      DO 310 i=1,neq
         y(i)=phi(i,1)
310      yprime(i)=0.0d0
      DO 330 j=2,kp1
         DO 320 i=1,neq
            y(i)=y(i)+phi(i,j)
320         yprime(i)=yprime(i)+gamma(j)*phi(i,j)
330   CONTINUE
      pnorm = ddanrm(neq,y,wt,rpar,ipar)
C
C
C
C     SOLVE THE CORRECTOR EQUATION USING A
C     MODIFIED NEWTON SCHEME.
      convgd= .true.
      m=0
      iwm(lnre)=iwm(lnre)+1
      ires = 0
      CALL res(x,y,yprime,delta,ires,rpar,ipar)
      IF (ires .LT. 0) GO TO 380
C
C
C     IF INDICATED,REEVALUATE THE
C     ITERATION MATRIX PD = DG/DY + CJ*DG/DYPRIME
C     (WHERE G(X,Y,YPRIME)=0). SET
C     JCALC TO 0 AS AN INDICATOR THAT
C     THIS HAS BEEN DONE.
      IF(jcalc .NE. -1)GO TO 340
      iwm(lnje)=iwm(lnje)+1
      jcalc=0
      CALL ddajac(neq,x,y,yprime,delta,cj,h,
     * ier,wt,e,wm,iwm,res,ires,uround,jac,rpar,
     * ipar,ntemp)
      cjold=cj
      s = 100.d0
      IF (ires .LT. 0) GO TO 380
      IF(ier .NE. 0)GO TO 380
      nsf=0
C
C
C     INITIALIZE THE ERROR ACCUMULATION VECTOR E.
340   CONTINUE
      DO 345 i=1,neq
345      e(i)=0.0d0
C
C
C     CORRECTOR LOOP.
350   CONTINUE
C
C     MULTIPLY RESIDUAL BY TEMP1 TO ACCELERATE CONVERGENCE
      temp1 = 2.0d0/(1.0d0 + cj/cjold)
      DO 355 i = 1,neq
355     delta(i) = delta(i) * temp1
C
C     COMPUTE A NEW ITERATE (BACK-SUBSTITUTION).
C     STORE THE CORRECTION IN DELTA.
      CALL ddaslv(neq,delta,wm,iwm)
C
C     UPDATE Y,E,AND YPRIME
      DO 360 i=1,neq
         y(i)=y(i)-delta(i)
         e(i)=e(i)-delta(i)
360      yprime(i)=yprime(i)-cj*delta(i)
C
C     TEST FOR CONVERGENCE OF THE ITERATION
      delnrm=ddanrm(neq,delta,wt,rpar,ipar)
      IF (delnrm .LE. 100.d0*uround*pnorm) GO TO 375
      IF (m .GT. 0) GO TO 365
         oldnrm = delnrm
         GO TO 367
365   rate = (delnrm/oldnrm)**(1.0d0/m)
      IF (rate .GT. 0.90d0) GO TO 370
      s = rate/(1.0d0 - rate)
367   IF (s*delnrm .LE. 0.33d0) GO TO 375
C
C     THE CORRECTOR HAS NOT YET CONVERGED.
C     UPDATE M AND TEST WHETHER THE
C     MAXIMUM NUMBER OF ITERATIONS HAVE
C     BEEN TRIED.
      m=m+1
      IF(m.GE.maxit)GO TO 370
C
C     EVALUATE THE RESIDUAL
C     AND GO BACK TO DO ANOTHER ITERATION
      iwm(lnre)=iwm(lnre)+1
      ires = 0
      CALL res(x,y,yprime,delta,ires,
     *  rpar,ipar)
      IF (ires .LT. 0) GO TO 380
      GO TO 350
C
C
C     THE CORRECTOR FAILED TO CONVERGE IN MAXIT
C     ITERATIONS. IF THE ITERATION MATRIX
C     IS NOT CURRENT,RE-DO THE STEP WITH
C     A NEW ITERATION MATRIX.
370   CONTINUE
      IF(jcalc.EQ.0)GO TO 380
      jcalc=-1
      GO TO 300
C
C
C     THE ITERATION HAS CONVERGED.  IF NONNEGATIVITY OF SOLUTION IS
C     REQUIRED, SET THE SOLUTION NONNEGATIVE, IF THE PERTURBATION
C     TO DO IT IS SMALL ENOUGH.  IF THE CHANGE IS TOO LARGE, THEN
C     CONSIDER THE CORRECTOR ITERATION TO HAVE FAILED.
375   IF(nonneg .EQ. 0) GO TO 390
      DO 377 i = 1,neq
377      delta(i) = min(y(i),0.0d0)
      delnrm = ddanrm(neq,delta,wt,rpar,ipar)
      IF(delnrm .GT. 0.33d0) GO TO 380
      DO 378 i = 1,neq
378      e(i) = e(i) - delta(i)
      GO TO 390
C
C
C     EXITS FROM BLOCK 3
C     NO CONVERGENCE WITH CURRENT ITERATION
C     MATRIX,OR SINGULAR ITERATION MATRIX
380   convgd= .false.
390   jcalc = 1
      IF(.NOT.convgd)GO TO 600
C
C
C
C
C
C-----------------------------------------------------------------------
C     BLOCK 4
C     ESTIMATE THE ERRORS AT ORDERS K,K-1,K-2
C     AS IF CONSTANT STEPSIZE WAS USED. ESTIMATE
C     THE LOCAL ERROR AT ORDER K AND TEST
C     WHETHER THE CURRENT STEP IS SUCCESSFUL.
C-----------------------------------------------------------------------
C
C     ESTIMATE ERRORS AT ORDERS K,K-1,K-2
      enorm = ddanrm(neq,e,wt,rpar,ipar)
      erk = sigma(k+1)*enorm
      terk = (k+1)*erk
      est = erk
      knew=k
      IF(k .EQ. 1)GO TO 430
      DO 405 i = 1,neq
405     delta(i) = phi(i,kp1) + e(i)
      erkm1=sigma(k)*ddanrm(neq,delta,wt,rpar,ipar)
      terkm1 = k*erkm1
      IF(k .GT. 2)GO TO 410
      IF(terkm1 .LE. 0.5d0*terk)GO TO 420
      GO TO 430
410   CONTINUE
      DO 415 i = 1,neq
415     delta(i) = phi(i,k) + delta(i)
      erkm2=sigma(k-1)*ddanrm(neq,delta,wt,rpar,ipar)
      terkm2 = (k-1)*erkm2
      IF(max(terkm1,terkm2).GT.terk)GO TO 430
C     LOWER THE ORDER
420   CONTINUE
      knew=k-1
      est = erkm1
C
C
C     CALCULATE THE LOCAL ERROR FOR THE CURRENT STEP
C     TO SEE IF THE STEP WAS SUCCESSFUL
430   CONTINUE
      err = ck * enorm
      IF(err .GT. 1.0d0)GO TO 600
C
C
C
C
C
C-----------------------------------------------------------------------
C     BLOCK 5
C     THE STEP IS SUCCESSFUL. DETERMINE
C     THE BEST ORDER AND STEPSIZE FOR
C     THE NEXT STEP. UPDATE THE DIFFERENCES
C     FOR THE NEXT STEP.
C-----------------------------------------------------------------------
      idid=1
      iwm(lnst)=iwm(lnst)+1
      kdiff=k-kold
      kold=k
      hold=h
C
C
C     ESTIMATE THE ERROR AT ORDER K+1 UNLESS:
C        ALREADY DECIDED TO LOWER ORDER, OR
C        ALREADY USING MAXIMUM ORDER, OR
C        STEPSIZE NOT CONSTANT, OR
C        ORDER RAISED IN PREVIOUS STEP
      IF(knew.EQ.km1.OR.k.EQ.iwm(lmxord))iphase=1
      IF(iphase .EQ. 0)GO TO 545
      IF(knew.EQ.km1)GO TO 540
      IF(k.EQ.iwm(lmxord)) GO TO 550
      IF(kp1.GE.ns.OR.kdiff.EQ.1)GO TO 550
      DO 510 i=1,neq
510      delta(i)=e(i)-phi(i,kp2)
      erkp1 = (1.0d0/(k+2))*ddanrm(neq,delta,wt,rpar,ipar)
      terkp1 = (k+2)*erkp1
      IF(k.GT.1)GO TO 520
      IF(terkp1.GE.0.5d0*terk)GO TO 550
      GO TO 530
520   IF(terkm1.LE.min(terk,terkp1))GO TO 540
      IF(terkp1.GE.terk.OR.k.EQ.iwm(lmxord))GO TO 550
C
C     RAISE ORDER
530   k=kp1
      est = erkp1
      GO TO 550
C
C     LOWER ORDER
540   k=km1
      est = erkm1
      GO TO 550
C
C     IF IPHASE = 0, INCREASE ORDER BY ONE AND MULTIPLY STEPSIZE BY
C     FACTOR TWO
545   k = kp1
      hnew = h*2.0d0
      h = hnew
      GO TO 575
C
C
C     DETERMINE THE APPROPRIATE STEPSIZE FOR
C     THE NEXT STEP.
550   hnew=h
      temp2=k+1
      r=(2.0d0*est+0.0001d0)**(-1.0d0/temp2)
      IF(r .LT. 2.0d0) GO TO 555
      hnew = 2.0d0*h
      GO TO 560
555   IF(r .GT. 1.0d0) GO TO 560
      r = max(0.5d0,min(0.9d0,r))
      hnew = h*r
560   h=hnew
C
C
C     UPDATE DIFFERENCES FOR NEXT STEP
575   CONTINUE
      IF(kold.EQ.iwm(lmxord))GO TO 585
      DO 580 i=1,neq
580      phi(i,kp2)=e(i)
585   CONTINUE
      DO 590 i=1,neq
590      phi(i,kp1)=phi(i,kp1)+e(i)
      DO 595 j1=2,kp1
         j=kp1-j1+1
         DO 595 i=1,neq
595      phi(i,j)=phi(i,j)+phi(i,j+1)
      RETURN
C
C
C
C
C
C-----------------------------------------------------------------------
C     BLOCK 6
C     THE STEP IS UNSUCCESSFUL. RESTORE X,PSI,PHI
C     DETERMINE APPROPRIATE STEPSIZE FOR
C     CONTINUING THE INTEGRATION, OR EXIT WITH
C     AN ERROR FLAG IF THERE HAVE BEEN MANY
C     FAILURES.
C-----------------------------------------------------------------------
600   iphase = 1
C
C     RESTORE X,PHI,PSI
      x=xold
      IF(kp1.LT.nsp1)GO TO 630
      DO 620 j=nsp1,kp1
         temp1=1.0d0/beta(j)
         DO 610 i=1,neq
610         phi(i,j)=temp1*phi(i,j)
620      CONTINUE
630   CONTINUE
      DO 640 i=2,kp1
640      psi(i-1)=psi(i)-h
C
C
C     TEST WHETHER FAILURE IS DUE TO CORRECTOR ITERATION
C     OR ERROR TEST
      IF(convgd)GO TO 660
      iwm(lctf)=iwm(lctf)+1
C
C
C     THE NEWTON ITERATION FAILED TO CONVERGE WITH
C     A CURRENT ITERATION MATRIX.  DETERMINE THE CAUSE
C     OF THE FAILURE AND TAKE APPROPRIATE ACTION.
      IF(ier.EQ.0)GO TO 650
C
C     THE ITERATION MATRIX IS SINGULAR. REDUCE
C     THE STEPSIZE BY A FACTOR OF 4. IF
C     THIS HAPPENS THREE TIMES IN A ROW ON
C     THE SAME STEP, RETURN WITH AN ERROR FLAG
      nsf=nsf+1
      r = 0.25d0
      h=h*r
      IF (nsf .LT. 3 .AND. abs(h) .GE. hmin) GO TO 690
      idid=-8
      GO TO 675
C
C
C     THE NEWTON ITERATION FAILED TO CONVERGE FOR A REASON
C     OTHER THAN A SINGULAR ITERATION MATRIX.  IF IRES = -2, THEN
C     RETURN.  OTHERWISE, REDUCE THE STEPSIZE AND TRY AGAIN, UNLESS
C     TOO MANY FAILURES HAVE OCCURRED.
650   CONTINUE
      IF (ires .GT. -2) GO TO 655
      idid = -11
      GO TO 675
655   ncf = ncf + 1
      r = 0.25d0
      h = h*r
      IF (ncf .LT. 10 .AND. abs(h) .GE. hmin) GO TO 690
      idid = -7
      IF (ires .LT. 0) idid = -10
      IF (nef .GE. 3) idid = -9
      GO TO 675
C
C
C     THE NEWTON SCHEME CONVERGED,AND THE CAUSE
C     OF THE FAILURE WAS THE ERROR ESTIMATE
C     EXCEEDING THE TOLERANCE.
660   nef=nef+1
      iwm(letf)=iwm(letf)+1
      IF (nef .GT. 1) GO TO 665
C
C     ON FIRST ERROR TEST FAILURE, KEEP CURRENT ORDER OR LOWER
C     ORDER BY ONE.  COMPUTE NEW STEPSIZE BASED ON DIFFERENCES
C     OF THE SOLUTION.
      k = knew
      temp2 = k + 1
      r = 0.90d0*(2.0d0*est+0.0001d0)**(-1.0d0/temp2)
      r = max(0.25d0,min(0.9d0,r))
      h = h*r
      IF (abs(h) .GE. hmin) GO TO 690
      idid = -6
      GO TO 675
C
C     ON SECOND ERROR TEST FAILURE, USE THE CURRENT ORDER OR
C     DECREASE ORDER BY ONE.  REDUCE THE STEPSIZE BY A FACTOR OF
C     FOUR.
665   IF (nef .GT. 2) GO TO 670
      k = knew
      h = 0.25d0*h
      IF (abs(h) .GE. hmin) GO TO 690
      idid = -6
      GO TO 675
C
C     ON THIRD AND SUBSEQUENT ERROR TEST FAILURES, SET THE ORDER TO
C     ONE AND REDUCE THE STEPSIZE BY A FACTOR OF FOUR.
670   k = 1
      h = 0.25d0*h
      IF (abs(h) .GE. hmin) GO TO 690
      idid = -6
      GO TO 675
C
C
C
C
C     FOR ALL CRASHES, RESTORE Y TO ITS LAST VALUE,
C     INTERPOLATE TO FIND YPRIME AT LAST X, AND RETURN
675   CONTINUE
      CALL ddatrp(x,x,y,yprime,neq,k,phi,psi)
      RETURN
C
C
C     GO BACK AND TRY THIS STEP AGAIN
690   GO TO 200
C
C------END OF SUBROUTINE DDASTP------
      END