dnsd.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 dnsd(X,Y,YPRIME,NEQ,RES,PDUM,WT,RPAR,IPAR,
     *   dumsvr,delta,e,wm,iwm,cj,dums,dumr,dume,epcon,
     *   s,confac,tolnew,muldel,maxit,ires,idum,iernew)
C
C***BEGIN PROLOGUE  DNSD
C***REFER TO  DDASPK
C***DATE WRITTEN   891219   (YYMMDD)
C***REVISION DATE  900926   (YYMMDD)
C***REVISION DATE  950126   (YYMMDD)
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C     DNSD 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     PDUM      -- Dummy argument.
C     WT        -- Vector of weights for error criterion.
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     DUMSVR    -- Dummy argument.
C     DELTA     -- Work vector for DNSD of length NEQ.
C     E         -- Error accumulation vector for DNSD 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 (step size).
C     DUMS      -- Dummy argument.
C     DUMR      -- Dummy argument.
C     DUME      -- Dummy argument.
C     EPCON     -- Tolerance to test for convergence of the Newton
C                  iteration.
C     S         -- Used for error convergence tests.
C                  In the Newton iteration: S = RATE/(1 - RATE),
C                  where RATE is the estimated rate of convergence
C                  of the Newton iteration.
C                  The calling routine passes the initial value
C                  of S to the Newton iteration.
C     CONFAC    -- A residual scale factor to improve convergence.
C     TOLNEW    -- Tolerance on the norm of Newton correction in
C                  alternative Newton convergence test.
C     MULDEL    -- A flag indicating whether or not to multiply
C                  DELTA by CONFAC.
C                  0  ==> do not scale DELTA by CONFAC.
C                  1  ==> scale DELTA by CONFAC.
C     MAXIT     -- Maximum allowed number of Newton iterations.
C     IRES      -- Error flag returned from RES.  See RES description
C                  in DDASPK prologue.  If IRES = -1, then IERNEW
C                  will be set to 1.
C                  If IRES < -1, then IERNEW will be set to -1.
C     IDUM      -- Dummy argument.
C     IERNEW    -- Error flag for Newton iteration.
C                   0  ==> Newton iteration converged.
C                   1  ==> recoverable error inside Newton iteration.
C                  -1  ==> unrecoverable error inside Newton iteration.
C
C     All arguments with "DUM" in their names are dummy arguments
C     which are not used in this routine.
C-----------------------------------------------------------------------
C
C***ROUTINES CALLED
C   DSLVD, DDWNRM, RES
C
C***END PROLOGUE  DNSD
C
C
      IMPLICIT DOUBLE PRECISION(a-h,o-z)
      dimension y(*),yprime(*),wt(*),delta(*),e(*)
      dimension wm(*),iwm(*), rpar(*),ipar(*)
      EXTERNAL  res
C
      parameter(lnre=12, lnni=19)
C
C     Initialize Newton counter M and accumulation vector E. 
C
      m = 0
      DO 100 i=1,neq
100     e(i)=0.0d0
C
C     Corrector loop.
C
300   CONTINUE
      iwm(lnni) = iwm(lnni) + 1
C
C     If necessary, multiply residual by convergence factor.
C
      IF (muldel .EQ. 1) THEN
         DO 320 i = 1,neq
320        delta(i) = delta(i) * confac
        ENDIF
C
C     Compute a new iterate (back-substitution).
C     Store the correction in DELTA.
C
      CALL dslvd(neq,delta,wm,iwm)
C
C     Update Y, E, and YPRIME.
C
      DO 340 i=1,neq
         y(i)=y(i)-delta(i)
         e(i)=e(i)-delta(i)
340      yprime(i)=yprime(i)-cj*delta(i)
C
C     Test for convergence of the iteration.
C
      delnrm=ddwnrm(neq,delta,wt,rpar,ipar)
      IF (delnrm .LE. tolnew) go to 370
      IF (m .EQ. 0) THEN
        oldnrm = delnrm
      ELSE
        rate = (delnrm/oldnrm)**(1.0d0/m)
        IF (rate .GT. 0.9d0) go to 380
        s = rate/(1.0d0 - rate)
      ENDIF
      IF (s*delnrm .LE. epcon) go to 370
C
C     The corrector has not yet converged.
C     Update M and test whether the
C     maximum number of iterations have
C     been tried.
C
      m=m+1
      IF(m.GE.maxit) go to 380
C
C     Evaluate the residual,
C     and go back to do another iteration.
C
      iwm(lnre)=iwm(lnre)+1
      CALL res(x,y,yprime,cj,delta,ires,rpar,ipar)
      IF (ires .LT. 0) go to 380
      go to 300
C
C     The iteration has converged.
C
370   RETURN
C
C     The iteration has not converged.  Set IERNEW appropriately.
C
380   CONTINUE
      IF (ires .LE. -2 ) THEN
         iernew = -1
      ELSE
         iernew = 1
      ENDIF
      RETURN
C
C
C------END OF SUBROUTINE DNSD-------------------------------------------
      END