d5/d46/dspigm_8f_source.html

 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

datv
subroutine datv(NEQ, Y, TN, YPRIME, SAVR, V, WGHT, YPTEM, RES, IRES, PSOL, Z, VTEM, WP, IWP, CJ, EPLIN, IER, NRE, NPSL, RPAR, IPAR)
Definition: datv.f:8

dhels
subroutine dhels(A, LDA, N, Q, B)
Definition: dhels.f:6

dheqr
subroutine dheqr(A, LDA, N, Q, INFO, IJOB)
Definition: dheqr.f:6

dorth
subroutine dorth(VNEW, V, HES, N, LL, LDHES, KMP, SNORMW)
Definition: dorth.f:6

dspigm
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)
Definition: dspigm.f:9