dslvk.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 dslvk (NEQ, Y, TN, YPRIME, SAVR, X, EWT, WM, IWM,
     *   res, ires, psol, iersl, cj, eplin, sqrtn, rsqrtn, rhok,
     *   rpar, ipar)
C
C***BEGIN PROLOGUE  DSLVK
C***REFER TO  DDASPK
C***DATE WRITTEN   890101   (YYMMDD)
C***REVISION DATE  900926   (YYMMDD)
C***REVISION DATE  940928   Removed MNEWT and added RHOK in call list.
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C DSLVK uses a restart algorithm and interfaces to DSPIGM for
C the solution of the linear system arising from a Newton iteration.
C
C In addition to variables described elsewhere,
C communication with DSLVK uses the following variables..
C WM    = Real work space containing data for the algorithm
C         (Krylov basis vectors, Hessenberg matrix, etc.).
C IWM   = Integer work space containing data for the algorithm.
C X     = The right-hand side vector on input, and the solution vector
C         on output, of length NEQ.
C IRES  = Error flag from RES.
C IERSL = Output flag ..
C         IERSL =  0 means no trouble occurred (or user RES routine
C                    returned IRES < 0)
C         IERSL =  1 means the iterative method failed to converge
C                    (DSPIGM returned IFLAG > 0.)
C         IERSL = -1 means there was a nonrecoverable error in the
C                    iterative solver, and an error exit will occur.
C-----------------------------------------------------------------------
C***ROUTINES CALLED
C   DSCAL, DCOPY, DSPIGM
C
C***END PROLOGUE  DSLVK
C
      INTEGER NEQ, IWM, IRES, IERSL, IPAR
      DOUBLE PRECISION Y, TN, YPRIME, SAVR, X, EWT, WM, CJ, EPLIN,
     1   sqrtn, rsqrtn, rhok, rpar
      dimension y(*), yprime(*), savr(*), x(*), ewt(*),
     1  wm(*), iwm(*), rpar(*), ipar(*)
C
      INTEGER IFLAG, IRST, NRSTS, NRMAX, LR, LDL, LHES, LGMR, LQ, LV,
     1        lwk, lz, maxlp1, npsl
      INTEGER NLI, NPS, NCFL, NRE, MAXL, KMP, MITER
      EXTERNAL  res, psol
C
      parameter(lnre=12, lncfl=16, lnli=20, lnps=21)
      parameter(llocwp=29, llciwp=30)
      parameter(lmiter=23, lmaxl=24, lkmp=25, lnrmax=26)
C
C-----------------------------------------------------------------------
C IRST is set to 1, to indicate restarting is in effect.
C NRMAX is the maximum number of restarts.
C-----------------------------------------------------------------------
      DATA irst/1/
C
      liwp = iwm(llciwp)
      nli = iwm(lnli)
      nps = iwm(lnps)
      ncfl = iwm(lncfl)
      nre = iwm(lnre)
      lwp = iwm(llocwp)
      maxl = iwm(lmaxl)
      kmp = iwm(lkmp)
      nrmax = iwm(lnrmax)
      miter = iwm(lmiter)
      iersl = 0
      ires = 0
C-----------------------------------------------------------------------
C Use a restarting strategy to solve the linear system
C P*X = -F.  Parse the work vector, and perform initializations.
C Note that zero is the initial guess for X.
C-----------------------------------------------------------------------
      maxlp1 = maxl + 1
      lv = 1
      lr = lv + neq*maxl
      lhes = lr + neq + 1
      lq = lhes + maxl*maxlp1
      lwk = lq + 2*maxl
      ldl = lwk + min0(1,maxl-kmp)*neq
      lz = ldl + neq
      CALL dscal(neq, rsqrtn, ewt, 1)
      CALL dcopy(neq, x, 1, wm(lr), 1)
      DO 110 i = 1,neq
 110     x(i) = 0.d0
C-----------------------------------------------------------------------
C Top of loop for the restart algorithm.  Initial pass approximates
C X and sets up a transformed system to perform subsequent restarts
C to update X.  NRSTS is initialized to -1, because restarting
C does not occur until after the first pass.
C Update NRSTS; conditionally copy DL to R; call the DSPIGM
C algorithm to solve A*Z = R;  updated counters;  update X with
C the residual solution.
C Note:  if convergence is not achieved after NRMAX restarts,
C then the linear solver is considered to have failed.
C-----------------------------------------------------------------------
      nrsts = -1
 115  CONTINUE
      nrsts = nrsts + 1
      IF (nrsts .GT. 0) CALL dcopy(neq, wm(ldl), 1, wm(lr),1)
      CALL dspigm(neq, tn, y, yprime, savr, wm(lr), ewt, maxl, maxlp1,
     1   kmp, eplin, cj, res, ires, nres, psol, npsl, wm(lz), wm(lv),
     2   wm(lhes), wm(lq), lgmr, wm(lwp), iwm(liwp), wm(lwk),
     3   wm(ldl), rhok, iflag, irst, nrsts, rpar, ipar)
      nli = nli + lgmr
      nps = nps + npsl
      nre = nre + nres
      DO 120 i = 1,neq
 120     x(i) = x(i) + wm(lz+i-1)
      IF ((iflag .EQ. 1) .AND. (nrsts .LT. nrmax) .AND. (ires .EQ. 0))
     1   go to 115
C-----------------------------------------------------------------------
C The restart scheme is finished.  Test IRES and IFLAG to see if
C convergence was not achieved, and set flags accordingly.
C-----------------------------------------------------------------------
      IF (ires .LT. 0) THEN
         ncfl = ncfl + 1
      ELSE IF (iflag .NE. 0) THEN
         ncfl = ncfl + 1
         IF (iflag .GT. 0) iersl = 1
         IF (iflag .LT. 0) iersl = -1
      ENDIF
C-----------------------------------------------------------------------
C Update IWM with counters, rescale EWT, and return.
C-----------------------------------------------------------------------
      iwm(lnli)  = nli
      iwm(lnps)  = nps
      iwm(lncfl) = ncfl
      iwm(lnre)  = nre
      CALL dscal(neq, sqrtn, ewt, 1)
      RETURN
C
C------END OF SUBROUTINE DSLVK------------------------------------------
      END