dnedd.f

Go to the documentation of this file.

C Work performed under the auspices of the U.S. Department of Energy
C by Lawrence Livermore National Laboratory under contract number 
C W-7405-Eng-48.
C
      SUBROUTINE dnedd(X,Y,YPRIME,NEQ,RES,JACD,PDUM,H,WT,
     *   jstart,idid,rpar,ipar,phi,gamma,dumsvr,delta,e,
     *   wm,iwm,cj,cjold,cjlast,s,uround,dume,dums,dumr,
     *   epcon,jcalc,jfdum,kp1,nonneg,ntype,iernls)
C
C***BEGIN PROLOGUE  DNEDD
C***REFER TO  DDASPK
C***DATE WRITTEN   891219   (YYMMDD)
C***REVISION DATE  900926   (YYMMDD)
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C     DNEDD solves a nonlinear system of
C     algebraic equations of the form
C     G(X,Y,YPRIME) = 0 for the unknown Y.
C
C     The method used is a modified Newton scheme.
C
C     The parameters represent
C
C     X         -- Independent variable.
C     Y         -- Solution vector.
C     YPRIME    -- Derivative of solution vector.
C     NEQ       -- Number of unknowns.
C     RES       -- External user-supplied subroutine
C                  to evaluate the residual.  See RES description
C                  in DDASPK prologue.
C     JACD      -- External user-supplied routine to evaluate the
C                  Jacobian.  See JAC description for the case
C                  INFO(12) = 0 in the DDASPK prologue.
C     PDUM      -- Dummy argument.
C     H         -- Appropriate step size for next step.
C     WT        -- Vector of weights for error criterion.
C     JSTART    -- Indicates first call to this routine.
C                  If JSTART = 0, then this is the first call,
C                  otherwise it is not.
C     IDID      -- Completion flag, output by DNEDD.
C                  See IDID description in DDASPK prologue.
C     RPAR,IPAR -- Real and integer arrays used for communication
C                  between the calling program and external user
C                  routines.  They are not altered within DASPK.
C     PHI       -- Array of divided differences used by
C                  DNEDD.  The length is NEQ*(K+1),where
C                  K is the maximum order.
C     GAMMA     -- Array used to predict Y and YPRIME.  The length
C                  is MAXORD+1 where MAXORD is the maximum order.
C     DUMSVR    -- Dummy argument.
C     DELTA     -- Work vector for NLS of length NEQ.
C     E         -- Error accumulation vector for NLS 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     CJ        -- Parameter always proportional to 1/H.
C     CJOLD     -- Saves the value of CJ as of the last call to DMATD.
C                  Accounts for changes in CJ needed to
C                  decide whether to call DMATD.
C     CJLAST    -- Previous value of CJ.
C     S         -- A scalar determined by the approximate rate
C                  of convergence of the Newton iteration and used
C                  in the convergence test for the Newton iteration.
C
C                  If RATE is defined to be an estimate of the
C                  rate of convergence of the Newton iteration,
C                  then S = RATE/(1.D0-RATE).
C
C                  The closer RATE is to 0., the faster the Newton
C                  iteration is converging; the closer RATE is to 1.,
C                  the slower the Newton iteration is converging.
C
C                  On the first Newton iteration with an up-dated
C                  preconditioner S = 100.D0, Thus the initial
C                  RATE of convergence is approximately 1.
C
C                  S is preserved from call to call so that the rate
C                  estimate from a previous step can be applied to
C                  the current step.
C     UROUND    -- Unit roundoff.
C     DUME      -- Dummy argument.
C     DUMS      -- Dummy argument.
C     DUMR      -- Dummy argument.
C     EPCON     -- Tolerance to test for convergence of the Newton
C                  iteration.
C     JCALC     -- Flag used to determine when to update
C                  the Jacobian matrix.  In general:
C
C                  JCALC = -1 ==> Call the DMATD routine to update
C                                 the Jacobian matrix.
C                  JCALC =  0 ==> Jacobian matrix is up-to-date.
C                  JCALC =  1 ==> Jacobian matrix is out-dated,
C                                 but DMATD will not be called unless
C                                 JCALC is set to -1.
C     JFDUM     -- Dummy argument.
C     KP1       -- The current order(K) + 1;  updated across calls.
C     NONNEG    -- Flag to determine nonnegativity constraints.
C     NTYPE     -- Identification code for the NLS routine.
C                   0  ==> modified Newton; direct solver.
C     IERNLS    -- Error flag for nonlinear solver.
C                   0  ==> nonlinear solver converged.
C                   1  ==> recoverable error inside nonlinear solver.
C                  -1  ==> unrecoverable error inside nonlinear solver.
C
C     All variables with "DUM" in their names are dummy variables
C     which are not used in this routine.
C
C     Following is a list and description of local variables which
C     may not have an obvious usage.  They are listed in roughly the
C     order they occur in this subroutine.
C
C     The following group of variables are passed as arguments to
C     the Newton iteration solver.  They are explained in greater detail
C     in DNSD:
C        TOLNEW, MULDEL, MAXIT, IERNEW
C
C     IERTYP -- Flag which tells whether this subroutine is correct.
C               0 ==> correct subroutine.
C               1 ==> incorrect subroutine.
C 
C-----------------------------------------------------------------------
C***ROUTINES CALLED
C   DDWNRM, RES, DMATD, DNSD
C
C***END PROLOGUE  DNEDD
C
C
      IMPLICIT DOUBLE PRECISION(a-h,o-z)
      dimension y(*),yprime(*),wt(*)
      dimension delta(*),e(*)
      dimension wm(*),iwm(*), rpar(*),ipar(*)
      dimension phi(neq,*),gamma(*)
      EXTERNAL  res, jacd
C
      parameter(lnre=12, lnje=13)
C
      SAVE muldel, maxit, xrate
      DATA muldel/1/, maxit/4/, xrate/0.25d0/
C
C     Verify that this is the correct subroutine.
C
      iertyp = 0
      IF (ntype .NE. 0) THEN
         iertyp = 1
         go to 380
         ENDIF
C
C     If this is the first step, perform initializations.
C
      IF (jstart .EQ. 0) THEN
         cjold = cj
         jcalc = -1
         ENDIF
C
C     Perform all other initializations.
C
      iernls = 0
C
C     Decide whether new Jacobian is needed.
C
      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-----------------------------------------------------------------------
C     Entry point for updating the Jacobian with current
C     stepsize.
C-----------------------------------------------------------------------
300   CONTINUE
C
C     Initialize all error flags to zero.
C
      ierj = 0
      ires = 0
      iernew = 0
C
C     Predict the solution and derivative and compute the tolerance
C     for the Newton iteration.
C
      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 = ddwnrm(neq,y,wt,rpar,ipar)
      tolnew = 100.d0*uround*pnorm
C     
C     Call RES to initialize DELTA.
C
      iwm(lnre)=iwm(lnre)+1
      CALL res(x,y,yprime,cj,delta,ires,rpar,ipar)
      IF (ires .LT. 0) go to 380
C
C     If indicated, reevaluate the iteration matrix 
C     J = dG/dY + CJ*dG/dYPRIME (where G(X,Y,YPRIME)=0).
C     Set JCALC to 0 as an indicator that this has been done.
C
      IF(jcalc .EQ. -1) THEN
         iwm(lnje)=iwm(lnje)+1
         jcalc=0
         CALL dmatd(neq,x,y,yprime,delta,cj,h,ierj,wt,e,wm,iwm,
     *              res,ires,uround,jacd,rpar,ipar)
         cjold=cj
         s = 100.d0
         IF (ires .LT. 0) go to 380
         IF(ierj .NE. 0)go to 380
      ENDIF
C
C     Call the nonlinear Newton solver.
C
      temp1 = 2.0d0/(1.0d0 + cj/cjold)
      CALL dnsd(x,y,yprime,neq,res,pdum,wt,rpar,ipar,dumsvr,
     *          delta,e,wm,iwm,cj,dums,dumr,dume,epcon,s,temp1,
     *          tolnew,muldel,maxit,ires,idum,iernew)
C
      IF (iernew .GT. 0 .AND. jcalc .NE. 0) THEN
C
C        The Newton iteration had a recoverable failure with an old
C        iteration matrix.  Retry the step with a new iteration matrix.
C
         jcalc = -1
         go to 300
      ENDIF
C
      IF (iernew .NE. 0) go to 380
C
C     The Newton iteration has converged.  If nonnegativity of
C     solution is required, set the solution nonnegative, if the
C     perturbation to do it is small enough.  If the change is too
C     large, then consider the corrector iteration to have failed.
C
375   IF(nonneg .EQ. 0) go to 390
      DO 377 i = 1,neq
377      delta(i) = min(y(i),0.0d0)
      delnrm = ddwnrm(neq,delta,wt,rpar,ipar)
      IF(delnrm .GT. epcon) go to 380
      DO 378 i = 1,neq
378      e(i) = e(i) - delta(i)
      go to 390
C
C
C     Exits from nonlinear solver.
C     No convergence with current iteration
C     matrix, or singular iteration matrix.
C     Compute IERNLS and IDID accordingly.
C
380   CONTINUE
      IF (ires .LE. -2 .OR. iertyp .NE. 0) THEN
         iernls = -1
         IF (ires .LE. -2) idid = -11
         IF (iertyp .NE. 0) idid = -15
      ELSE
         iernls = 1
         IF (ires .LT. 0) idid = -10
         IF (ierj .NE. 0) idid = -8
      ENDIF
C
390   jcalc = 1
      RETURN
C
C------END OF SUBROUTINE DNEDD------------------------------------------
      END