ddaini.f

Go to the documentation of this file.

      SUBROUTINE ddaini (X, Y, YPRIME, NEQ, RES, JAC, H, WT, IDID, RPAR,
     +   ipar, phi, delta, e, wm, iwm, hmin, uround, nonneg, ntemp)
C***BEGIN PROLOGUE  DDAINI
C***SUBSIDIARY
C***PURPOSE  Initialization routine for DDASSL.
C***LIBRARY   SLATEC (DASSL)
C***TYPE      DOUBLE PRECISION (SDAINI-S, DDAINI-D)
C***AUTHOR  PETZOLD, LINDA R., (LLNL)
C***DESCRIPTION
C-----------------------------------------------------------------
C     DDAINI TAKES ONE STEP OF SIZE H OR SMALLER
C     WITH THE BACKWARD EULER METHOD, TO
C     FIND YPRIME.  X AND Y ARE UPDATED TO BE CONSISTENT WITH THE
C     NEW STEP.  A MODIFIED DAMPED NEWTON ITERATION IS USED TO
C     SOLVE THE CORRECTOR ITERATION.
C
C     THE INITIAL GUESS FOR YPRIME IS USED IN THE
C     PREDICTION, AND IN FORMING THE ITERATION
C     MATRIX, BUT IS NOT INVOLVED IN THE
C     ERROR TEST. THIS MAY HAVE TROUBLE
C     CONVERGING IF THE INITIAL GUESS IS NO
C     GOOD, OR IF G(X,Y,YPRIME) DEPENDS
C     NONLINEARLY ON YPRIME.
C
C     THE PARAMETERS REPRESENT:
C     X --         INDEPENDENT VARIABLE
C     Y --         SOLUTION VECTOR AT X
C     YPRIME --    DERIVATIVE OF SOLUTION VECTOR
C     NEQ --       NUMBER OF EQUATIONS
C     H --         STEPSIZE. IMDER MAY USE A STEPSIZE
C                  SMALLER THAN H.
C     WT --        VECTOR OF WEIGHTS FOR ERROR
C                  CRITERION
C     IDID --      COMPLETION CODE WITH THE FOLLOWING MEANINGS
C                  IDID= 1 -- YPRIME WAS FOUND SUCCESSFULLY
C                  IDID=-12 -- DDAINI FAILED TO FIND YPRIME
C     RPAR,IPAR -- REAL AND INTEGER PARAMETER ARRAYS
C                  THAT ARE NOT ALTERED BY DDAINI
C     PHI --       WORK SPACE FOR DDAINI
C     DELTA,E --   WORK SPACE FOR DDAINI
C     WM,IWM --    REAL AND INTEGER ARRAYS STORING
C                  MATRIX INFORMATION
C
C-----------------------------------------------------------------
C***ROUTINES CALLED  DDAJAC, DDANRM, DDASLV
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   901030  Minor corrections to declarations.  (FNF)
C***END PROLOGUE  DDAINI
C
      INTEGER  neq, idid, ipar(*), iwm(*), nonneg, ntemp
      double precision
     *   x, y(*), yprime(*), h, wt(*), rpar(*), phi(neq,*), delta(*),
     *   e(*), wm(*), hmin, uround
      EXTERNAL  res, jac
C
      EXTERNAL  ddajac, ddanrm, ddaslv
      DOUBLE PRECISION  ddanrm
C
      INTEGER  i, ier, ires, jcalc, lnje, lnre, m, maxit, mjac, ncf,
     *   nef, nsf
      double precision
     *   cj, damp, delnrm, err, oldnrm, r, rate, s, xold, ynorm
      LOGICAL  convgd
C
      parameter(lnre=12)
      parameter(lnje=13)
C
      DATA maxit/10/,mjac/5/
      DATA damp/0.75d0/
C
C
C---------------------------------------------------
C     BLOCK 1.
C     INITIALIZATIONS.
C---------------------------------------------------
C
C***FIRST EXECUTABLE STATEMENT  DDAINI
      idid=1
      nef=0
      ncf=0
      nsf=0
      xold=x
      ynorm=ddanrm(neq,y,wt,rpar,ipar)
C
C     SAVE Y AND YPRIME IN PHI
      DO 100 i=1,neq
         phi(i,1)=y(i)
100      phi(i,2)=yprime(i)
C
C
C----------------------------------------------------
C     BLOCK 2.
C     DO ONE BACKWARD EULER STEP.
C----------------------------------------------------
C
C     SET UP FOR START OF CORRECTOR ITERATION
200   cj=1.0d0/h
      x=x+h
C
C     PREDICT SOLUTION AND DERIVATIVE
      DO 250 i=1,neq
250     y(i)=y(i)+h*yprime(i)
C
      jcalc=-1
      m=0
      convgd=.true.
C
C
C     CORRECTOR LOOP.
300   iwm(lnre)=iwm(lnre)+1
      ires=0
C
      CALL res(x,y,yprime,delta,ires,rpar,ipar)
      IF (ires.LT.0) go to 430
C
C
C     EVALUATE THE ITERATION MATRIX
      IF (jcalc.NE.-1) go to 310
      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)
C
      s=1000000.d0
      IF (ires.LT.0) go to 430
      IF (ier.NE.0) go to 430
      nsf=0
C
C
C
C     MULTIPLY RESIDUAL BY DAMPING FACTOR
310   CONTINUE
      DO 320 i=1,neq
320      delta(i)=delta(i)*damp
C
C     COMPUTE A NEW ITERATE (BACK SUBSTITUTION)
C     STORE THE CORRECTION IN DELTA
C
      CALL ddaslv(neq,delta,wm,iwm)
C
C     UPDATE Y AND YPRIME
      DO 330 i=1,neq
         y(i)=y(i)-delta(i)
330      yprime(i)=yprime(i)-cj*delta(i)
C
C     TEST FOR CONVERGENCE OF THE ITERATION.
C
      delnrm=ddanrm(neq,delta,wt,rpar,ipar)
      IF (delnrm.LE.100.d0*uround*ynorm)
     *   go to 400
C
      IF (m.GT.0) go to 340
         oldnrm=delnrm
         go to 350
C
340   rate=(delnrm/oldnrm)**(1.0d0/m)
      IF (rate.GT.0.90d0) go to 430
      s=rate/(1.0d0-rate)
C
350   IF (s*delnrm .LE. 0.33d0) go to 400
C
C
C     THE CORRECTOR HAS NOT YET CONVERGED. UPDATE
C     M AND AND TEST WHETHER THE MAXIMUM
C     NUMBER OF ITERATIONS HAVE BEEN TRIED.
C     EVERY MJAC ITERATIONS, GET A NEW
C     ITERATION MATRIX.
C
      m=m+1
      IF (m.GE.maxit) go to 430
C
      IF ((m/mjac)*mjac.EQ.m) jcalc=-1
      go to 300
C
C
C     THE ITERATION HAS CONVERGED.
C     CHECK NONNEGATIVITY CONSTRAINTS
400   IF (nonneg.EQ.0) go to 450
      DO 410 i=1,neq
410      delta(i)=min(y(i),0.0d0)
C
      delnrm=ddanrm(neq,delta,wt,rpar,ipar)
      IF (delnrm.GT.0.33d0) go to 430
C
      DO 420 i=1,neq
         y(i)=y(i)-delta(i)
420      yprime(i)=yprime(i)-cj*delta(i)
      go to 450
C
C
C     EXITS FROM CORRECTOR LOOP.
430   convgd=.false.
450   IF (.NOT.convgd) go to 600
C
C
C
C-----------------------------------------------------
C     BLOCK 3.
C     THE CORRECTOR ITERATION CONVERGED.
C     DO ERROR TEST.
C-----------------------------------------------------
C
      DO 510 i=1,neq
510      e(i)=y(i)-phi(i,1)
      err=ddanrm(neq,e,wt,rpar,ipar)
C
      IF (err.LE.1.0d0) RETURN
C
C
C
C--------------------------------------------------------
C     BLOCK 4.
C     THE BACKWARD EULER STEP FAILED. RESTORE X, Y
C     AND YPRIME TO THEIR ORIGINAL VALUES.
C     REDUCE STEPSIZE AND TRY AGAIN, IF
C     POSSIBLE.
C---------------------------------------------------------
C
600   CONTINUE
      x = xold
      DO 610 i=1,neq
         y(i)=phi(i,1)
610      yprime(i)=phi(i,2)
C
      IF (convgd) go to 640
      IF (ier.EQ.0) go to 620
         nsf=nsf+1
         h=h*0.25d0
         IF (nsf.LT.3.AND.abs(h).GE.hmin) go to 690
         idid=-12
         RETURN
620   IF (ires.GT.-2) go to 630
         idid=-12
         RETURN
630   ncf=ncf+1
      h=h*0.25d0
      IF (ncf.LT.10.AND.abs(h).GE.hmin) go to 690
         idid=-12
         RETURN
C
640   nef=nef+1
      r=0.90d0/(2.0d0*err+0.0001d0)
      r=max(0.1d0,min(0.5d0,r))
      h=h*r
      IF (abs(h).GE.hmin.AND.nef.LT.10) go to 690
         idid=-12
         RETURN
690      go to 200
C
C-------------END OF SUBROUTINE DDAINI----------------------
      END