dspigm.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 dspigm (NEQ, TN, Y, YPRIME, SAVR, R, WGHT, MAXL,
     *   maxlp1, kmp, eplin, cj, res, ires, nre, psol, npsl, z, v,
     *   hes, q, lgmr, wp, iwp, wk, dl, rhok, iflag, irst, nrsts,
     *   rpar, ipar)
C
C***BEGIN PROLOGUE  DSPIGM
C***DATE WRITTEN   890101   (YYMMDD)
C***REVISION DATE  900926   (YYMMDD)
C***REVISION DATE  940927   Removed MNEWT and added RHOK in call list.
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C This routine solves the linear system A * Z = R using a scaled
C preconditioned version of the generalized minimum residual method.
C An initial guess of Z = 0 is assumed.
C
C      On entry
C
C          NEQ = Problem size, passed to PSOL.
C
C           TN = Current Value of T.
C
C            Y = Array Containing current dependent variable vector.
C
C       YPRIME = Array Containing current first derivative of Y.
C
C         SAVR = Array containing current value of G(T,Y,YPRIME).
C
C            R = The right hand side of the system A*Z = R.
C                R is also used as work space when computing
C                the final approximation and will therefore be
C                destroyed.
C                (R is the same as V(*,MAXL+1) in the call to DSPIGM.)
C
C         WGHT = The vector of length NEQ containing the nonzero
C                elements of the diagonal scaling matrix.
C
C         MAXL = The maximum allowable order of the matrix H.
C
C       MAXLP1 = MAXL + 1, used for dynamic dimensioning of HES.
C
C          KMP = The number of previous vectors the new vector, VNEW,
C                must be made orthogonal to.  (KMP .LE. MAXL.)
C
C        EPLIN = Tolerance on residuals R-A*Z in weighted rms norm.
C
C           CJ = Scalar proportional to current value of
C                1/(step size H).
C
C           WK = Real work array used by routine DATV and PSOL.
C
C           DL = Real work array used for calculation of the residual
C                norm RHO when the method is incomplete (KMP.LT.MAXL)
C                and/or when using restarting.
C
C           WP = Real work array used by preconditioner PSOL.
C
C          IWP = Integer work array used by preconditioner PSOL.
C
C         IRST = Method flag indicating if restarting is being
C                performed.  IRST .GT. 0 means restarting is active,
C                while IRST = 0 means restarting is not being used.
C
C        NRSTS = Counter for the number of restarts on the current
C                call to DSPIGM.  If NRSTS .GT. 0, then the residual
C                R is already scaled, and so scaling of R is not
C                necessary.
C
C
C      On Return
C
C         Z    = The final computed approximation to the solution
C                of the system A*Z = R.
C
C         LGMR = The number of iterations performed and
C                the current order of the upper Hessenberg
C                matrix HES.
C
C         NRE  = The number of calls to RES (i.e. DATV)
C
C         NPSL = The number of calls to PSOL.
C
C         V    = The neq by (LGMR+1) array containing the LGMR
C                orthogonal vectors V(*,1) to V(*,LGMR).
C
C         HES  = The upper triangular factor of the QR decomposition
C                of the (LGMR+1) by LGMR upper Hessenberg matrix whose
C                entries are the scaled inner-products of A*V(*,I)
C                and V(*,K).
C
C         Q    = Real array of length 2*MAXL containing the components
C                of the givens rotations used in the QR decomposition
C                of HES.  It is loaded in DHEQR and used in DHELS.
C
C         IRES = Error flag from RES.
C
C           DL = Scaled preconditioned residual,
C                (D-inverse)*(P-inverse)*(R-A*Z). Only loaded when
C                performing restarts of the Krylov iteration.
C
C         RHOK = Weighted norm of final preconditioned residual.
C
C        IFLAG = Integer error flag..
C                0 Means convergence in LGMR iterations, LGMR.LE.MAXL.
C                1 Means the convergence test did not pass in MAXL
C                  iterations, but the new residual norm (RHO) is
C                  .LT. the old residual norm (RNRM), and so Z is
C                  computed.
C                2 Means the convergence test did not pass in MAXL
C                  iterations, new residual norm (RHO) .GE. old residual
C                  norm (RNRM), and the initial guess, Z = 0, is
C                  returned.
C                3 Means there was a recoverable error in PSOL
C                  caused by the preconditioner being out of date.
C               -1 Means there was an unrecoverable error in PSOL.
C
C-----------------------------------------------------------------------
C***ROUTINES CALLED
C   PSOL, DNRM2, DSCAL, DATV, DORTH, DHEQR, DCOPY, DHELS, DAXPY
C
C***END PROLOGUE  DSPIGM
C
      INTEGER NEQ,MAXL,MAXLP1,KMP,IRES,NRE,NPSL,LGMR,IWP,
     1   iflag,irst,nrsts,ipar
      DOUBLE PRECISION TN,Y,YPRIME,SAVR,R,WGHT,EPLIN,CJ,Z,V,HES,Q,WP,WK,
     1   dl,rhok,rpar
      dimension y(*), yprime(*), savr(*), r(*), wght(*), z(*),
     1   v(neq,*), hes(maxlp1,*), q(*), wp(*), iwp(*), wk(*), dl(*),
     2   rpar(*), ipar(*)
      INTEGER I, IER, INFO, IP1, I2, J, K, LL, LLP1
      DOUBLE PRECISION RNRM,C,DLNRM,PROD,RHO,S,SNORMW,DNRM2,TEM
      EXTERNAL  res, psol
C
      ier = 0
      iflag = 0
      lgmr = 0
      npsl = 0
      nre = 0
C-----------------------------------------------------------------------
C The initial guess for Z is 0.  The initial residual is therefore
C the vector R.  Initialize Z to 0.
C-----------------------------------------------------------------------
      DO 10 i = 1,neq
 10     z(i) = 0.0d0
C-----------------------------------------------------------------------
C Apply inverse of left preconditioner to vector R if NRSTS .EQ. 0.
C Form V(*,1), the scaled preconditioned right hand side.
C-----------------------------------------------------------------------
      IF (nrsts .EQ. 0) THEN
         CALL psol(neq, tn, y, yprime, savr, wk, cj, wght, wp, iwp,
     1      r, eplin, ier, rpar, ipar)
         npsl = 1
         IF (ier .NE. 0) go to 300
         DO 30 i = 1,neq
 30         v(i,1) = r(i)*wght(i)
      ELSE
         DO 35 i = 1,neq
 35         v(i,1) = r(i)
      ENDIF
C-----------------------------------------------------------------------
C Calculate norm of scaled vector V(*,1) and normalize it
C If, however, the norm of V(*,1) (i.e. the norm of the preconditioned
C residual) is .le. EPLIN, then return with Z=0.
C-----------------------------------------------------------------------
      rnrm = dnrm2(neq, v, 1)
      IF (rnrm .LE. eplin) THEN
        rhok = rnrm
        RETURN
        ENDIF
      tem = 1.0d0/rnrm
      CALL dscal(neq, tem, v(1,1), 1)
C-----------------------------------------------------------------------
C Zero out the HES array.
C-----------------------------------------------------------------------
      DO 65 j = 1,maxl
        DO 60 i = 1,maxlp1
 60       hes(i,j) = 0.0d0
 65     CONTINUE
C-----------------------------------------------------------------------
C Main loop to compute the vectors V(*,2) to V(*,MAXL).
C The running product PROD is needed for the convergence test.
C-----------------------------------------------------------------------
      prod = 1.0d0
      DO 90 ll = 1,maxl
        lgmr = ll
C-----------------------------------------------------------------------
C Call routine DATV to compute VNEW = ABAR*V(LL), where ABAR is
C the matrix A with scaling and inverse preconditioner factors applied.
C Call routine DORTH to orthogonalize the new vector VNEW = V(*,LL+1).
C call routine DHEQR to update the factors of HES.
C-----------------------------------------------------------------------
        CALL datv(neq, y, tn, yprime, savr, v(1,ll), wght, z,
     1     res, ires, psol, v(1,ll+1), wk, wp, iwp, cj, eplin,
     1     ier, nre, npsl, rpar, ipar)
        IF (ires .LT. 0) RETURN
        IF (ier .NE. 0) go to 300
        CALL dorth(v(1,ll+1), v, hes, neq, ll, maxlp1, kmp, snormw)
        hes(ll+1,ll) = snormw
        CALL dheqr(hes, maxlp1, ll, q, info, ll)
        IF (info .EQ. ll) go to 120
C-----------------------------------------------------------------------
C Update RHO, the estimate of the norm of the residual R - A*ZL.
C If KMP .LT. MAXL, then the vectors V(*,1),...,V(*,LL+1) are not
C necessarily orthogonal for LL .GT. KMP.  The vector DL must then
C be computed, and its norm used in the calculation of RHO.
C-----------------------------------------------------------------------
        prod = prod*q(2*ll)
        rho = abs(prod*rnrm)
        IF ((ll.GT.kmp) .AND. (kmp.LT.maxl)) THEN
          IF (ll .EQ. kmp+1) THEN
            CALL dcopy(neq, v(1,1), 1, dl, 1)
            DO 75 i = 1,kmp
              ip1 = i + 1
              i2 = i*2
              s = q(i2)
              c = q(i2-1)
              DO 70 k = 1,neq
 70             dl(k) = s*dl(k) + c*v(k,ip1)
 75           CONTINUE
            ENDIF
          s = q(2*ll)
          c = q(2*ll-1)/snormw
          llp1 = ll + 1
          DO 80 k = 1,neq
 80         dl(k) = s*dl(k) + c*v(k,llp1)
          dlnrm = dnrm2(neq, dl, 1)
          rho = rho*dlnrm
          ENDIF
C-----------------------------------------------------------------------
C Test for convergence.  If passed, compute approximation ZL.
C If failed and LL .LT. MAXL, then continue iterating.
C-----------------------------------------------------------------------
        IF (rho .LE. eplin) go to 200
        IF (ll .EQ. maxl) go to 100
C-----------------------------------------------------------------------
C Rescale so that the norm of V(1,LL+1) is one.
C-----------------------------------------------------------------------
        tem = 1.0d0/snormw
        CALL dscal(neq, tem, v(1,ll+1), 1)
 90     CONTINUE
 100  CONTINUE
      IF (rho .LT. rnrm) go to 150
 120  CONTINUE
      iflag = 2
      DO 130 i = 1,neq
 130     z(i) = 0.d0
      RETURN
 150  iflag = 1
C-----------------------------------------------------------------------
C The tolerance was not met, but the residual norm was reduced.
C If performing restarting (IRST .gt. 0) calculate the residual vector
C RL and store it in the DL array.  If the incomplete version is
C being used (KMP .lt. MAXL) then DL has already been calculated.
C-----------------------------------------------------------------------
      IF (irst .GT. 0) THEN
         IF (kmp .EQ. maxl) THEN
C
C           Calculate DL from the V(I)'s.
C
            CALL dcopy(neq, v(1,1), 1, dl, 1)
            maxlm1 = maxl - 1
            DO 175 i = 1,maxlm1
               ip1 = i + 1
               i2 = i*2
               s = q(i2)
               c = q(i2-1)
               DO 170 k = 1,neq
 170              dl(k) = s*dl(k) + c*v(k,ip1)
 175        CONTINUE
            s = q(2*maxl)
            c = q(2*maxl-1)/snormw
            DO 180 k = 1,neq
 180           dl(k) = s*dl(k) + c*v(k,maxlp1)
         ENDIF
C
C        Scale DL by RNRM*PROD to obtain the residual RL.
C
         tem = rnrm*prod
         CALL dscal(neq, tem, dl, 1)
      ENDIF
C-----------------------------------------------------------------------
C Compute the approximation ZL to the solution.
C Since the vector Z was used as work space, and the initial guess
C of the Newton correction is zero, Z must be reset to zero.
C-----------------------------------------------------------------------
 200  CONTINUE
      ll = lgmr
      llp1 = ll + 1
      DO 210 k = 1,llp1
 210    r(k) = 0.0d0
      r(1) = rnrm
      CALL dhels(hes, maxlp1, ll, q, r)
      DO 220 k = 1,neq
 220    z(k) = 0.0d0
      DO 230 i = 1,ll
        CALL daxpy(neq, r(i), v(1,i), 1, z, 1)
 230    CONTINUE
      DO 240 i = 1,neq
 240    z(i) = z(i)/wght(i)
C Load RHO into RHOK.
      rhok = rho
      RETURN
C-----------------------------------------------------------------------
C This block handles error returns forced by routine PSOL.
C-----------------------------------------------------------------------
 300  CONTINUE
      IF (ier .LT. 0) iflag = -1
      IF (ier .GT. 0) iflag = 3
C
      RETURN
C
C------END OF SUBROUTINE DSPIGM-----------------------------------------
      END