dnsk.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 dnsk(X,Y,YPRIME,NEQ,RES,PSOL,WT,RPAR,IPAR,
     *   SAVR,DELTA,E,WM,IWM,CJ,SQRTN,RSQRTN,EPLIN,EPCON,
     *   S,CONFAC,TOLNEW,MULDEL,MAXIT,IRES,IERSL,IERNEW)
C
C***BEGIN PROLOGUE  DNSK
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     DNSK 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     PSOL      -- External user-supplied routine to solve
C                  a linear system using preconditioning.
C                  See explanation inside DDASPK.
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     SAVR      -- Work vector for DNSK of length NEQ.
C     DELTA     -- Work vector for DNSK of length NEQ.
C     E         -- Error accumulation vector for DNSK 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     SQRTN     -- Square root of NEQ.
C     RSQRTN    -- reciprical of square root of NEQ.
C     EPLIN     -- Tolerance for linear system solver.
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.D0-RATE),
C                  where RATE is the estimated rate of convergence
C                  of the Newton iteration.
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                  The calling routine sends 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     IERSL     -- Error flag for linear system solver.
C                  See IERSL description in subroutine DSLVK.
C                  If IERSL = 1, then IERNEW will be set to 1.
C                  If IERSL < 0, then IERNEW will be set to -1.
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
C***ROUTINES CALLED
C   RES, DSLVK, DDWNRM
C
C***END PROLOGUE  DNSK
C
C
      IMPLICIT DOUBLE PRECISION(a-h,o-z)
      dimension y(*),yprime(*),wt(*),delta(*),e(*),savr(*)
      dimension wm(*),iwm(*), rpar(*),ipar(*)
      EXTERNAL  res, psol
C
      parameter(lnni=19, lnre=12)
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     Save residual in SAVR.
C
      DO 340 i = 1,neq
340     savr(i) = delta(i)
C
C     Compute a new iterate.  Store the correction in DELTA.
C
      CALL dslvk (neq, y, x, yprime, savr, delta, wt, wm, iwm,
     *   res, ires, psol, iersl, cj, eplin, sqrtn, rsqrtn, rhok,
     *   rpar, ipar)
      IF (ires .NE. 0 .OR. iersl .NE. 0) GO TO 380
C
C     Update Y, E, and YPRIME.
C
      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.
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.  Update M and test whether
C     the maximum number of iterations have been tried.
C
      m = m + 1
      IF (m .GE. maxit) GO TO 380
C
C     Evaluate the residual, 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 .OR. iersl .LT. 0) THEN
         iernew = -1
      ELSE
         iernew = 1
      ENDIF
      RETURN
C
C
C------END OF SUBROUTINE DNSK-------------------------------------------
      END