share/odepack/fortran/dspigmr.f

   1 *DECK DSPIGMR
   2       SUBROUTINE DSPIGMR (NEQ, TN, Y, SAVF, B, WGHT, N, MAXL, MAXLP1,
   3      1  KMP, DELTA, HL0, JPRE, MNEWT, F, PSOL, NPSL, X, V, HES, Q,
   4      2  LGMR, WP, IWP, WK, DL, IFLAG)
   5       EXTERNAL F, PSOL
   6       INTEGER NEQ,N,MAXL,MAXLP1,KMP,JPRE,MNEWT,NPSL,LGMR,IWP,IFLAG
   7       DOUBLE PRECISION TN,Y,SAVF,B,WGHT,DELTA,HL0,X,V,HES,Q,WP,WK,DL
   8       DIMENSION NEQ(*), Y(*), SAVF(*), B(*), WGHT(*), X(*), V(N,*),
   9      1    HES(MAXLP1,*), Q(*), WP(*), IWP(*), WK(*), DL(*)
  10 C-----------------------------------------------------------------------
  11 C This routine solves the linear system A * x = b using a scaled
  12 C preconditioned version of the Generalized Minimal Residual method.
  13 C An initial guess of x = 0 is assumed.
  14 C-----------------------------------------------------------------------
  15 C
  16 C      On entry
  17 C
  18 C          NEQ = problem size, passed to F and PSOL (NEQ(1) = N).
  19 C
  20 C           TN = current value of t.
  21 C
  22 C            Y = array containing current dependent variable vector.
  23 C
  24 C         SAVF = array containing current value of f(t,y).
  25 C
  26 C            B = the right hand side of the system A*x = b.
  27 C                B is also used as work space when computing
  28 C                the final approximation.
  29 C                (B is the same as V(*,MAXL+1) in the call to DSPIGMR.)
  30 C
  31 C         WGHT = the vector of length N containing the nonzero
  32 C                elements of the diagonal scaling matrix.
  33 C
  34 C            N = the order of the matrix A, and the lengths
  35 C                of the vectors WGHT, B and X.
  36 C
  37 C         MAXL = the maximum allowable order of the matrix HES.
  38 C
  39 C       MAXLP1 = MAXL + 1, used for dynamic dimensioning of HES.
  40 C
  41 C          KMP = the number of previous vectors the new vector VNEW
  42 C                must be made orthogonal to.  KMP .le. MAXL.
  43 C
  44 C        DELTA = tolerance on residuals b - A*x in weighted RMS-norm.
  45 C
  46 C          HL0 = current value of (step size h) * (coefficient l0).
  47 C
  48 C         JPRE = preconditioner type flag.
  49 C
  50 C        MNEWT = Newton iteration counter (.ge. 0).
  51 C
  52 C           WK = real work array used by routine DATV and PSOL.
  53 C
  54 C           DL = real work array used for calculation of the residual
  55 C                norm RHO when the method is incomplete (KMP .lt. MAXL).
  56 C                Not needed or referenced in complete case (KMP = MAXL).
  57 C
  58 C           WP = real work array used by preconditioner PSOL.
  59 C
  60 C          IWP = integer work array used by preconditioner PSOL.
  61 C
  62 C      On return
  63 C
  64 C         X    = the final computed approximation to the solution
  65 C                of the system A*x = b.
  66 C
  67 C         LGMR = the number of iterations performed and
  68 C                the current order of the upper Hessenberg
  69 C                matrix HES.
  70 C
  71 C         NPSL = the number of calls to PSOL.
  72 C
  73 C         V    = the N by (LGMR+1) array containing the LGMR
  74 C                orthogonal vectors V(*,1) to V(*,LGMR).
  75 C
  76 C         HES  = the upper triangular factor of the QR decomposition
  77 C                of the (LGMR+1) by lgmr upper Hessenberg matrix whose
  78 C                entries are the scaled inner-products of A*V(*,i)
  79 C                and V(*,k).
  80 C
  81 C         Q    = real array of length 2*MAXL containing the components
  82 C                of the Givens rotations used in the QR decomposition
  83 C                of HES.  It is loaded in DHEQR and used in DHELS.
  84 C
  85 C        IFLAG = integer error flag:
  86 C                0 means convergence in LGMR iterations, LGMR .le. MAXL.
  87 C                1 means the convergence test did not pass in MAXL
  88 C                  iterations, but the residual norm is .lt. 1,
  89 C                  or .lt. norm(b) if MNEWT = 0, and so x is computed.
  90 C                2 means the convergence test did not pass in MAXL
  91 C                  iterations, residual .gt. 1, and X is undefined.
  92 C                3 means there was a recoverable error in PSOL
  93 C                  caused by the preconditioner being out of date.
  94 C               -1 means there was a nonrecoverable error in PSOL.
  95 C
  96 C-----------------------------------------------------------------------
  97       INTEGER I, IER, INFO, IP1, I2, J, K, LL, LLP1
  98       DOUBLE PRECISION BNRM,BNRM0,C,DLNRM,PROD,RHO,S,SNORMW,DNRM2,TEM
  99 C
 100       IFLAG = 0
 101       LGMR = 0
 102       NPSL = 0
 103 C-----------------------------------------------------------------------
 104 C The initial residual is the vector b.  Apply scaling to b, and test
 105 C for an immediate return with X = 0 or X = b.
 106 C-----------------------------------------------------------------------
 107       DO 10 I = 1,N
 108  10     V(I,1) = B(I)*WGHT(I)
 109       BNRM0 = DNRM2 (N, V, 1)
 110       BNRM = BNRM0
 111       IF (BNRM0 .GT. DELTA) GO TO 30
 112       IF (MNEWT .GT. 0) GO TO 20
 113       CALL DCOPY (N, B, 1, X, 1)
 114       RETURN
 115  20   DO 25 I = 1,N
 116  25     X(I) = 0.0D0
 117       RETURN
 118  30   CONTINUE
 119 C Apply inverse of left preconditioner to vector b. --------------------
 120       IER = 0
 121       IF (JPRE .EQ. 0 .OR. JPRE .EQ. 2) GO TO 55
 122       CALL PSOL (NEQ, TN, Y, SAVF, WK, HL0, WP, IWP, B, 1, IER)
 123       NPSL = 1
 124       IF (IER .NE. 0) GO TO 300
 125 C Calculate norm of scaled vector V(*,1) and normalize it. -------------
 126       DO 50 I = 1,N
 127  50     V(I,1) = B(I)*WGHT(I)
 128       BNRM = DNRM2 (N, V, 1)
 129       DELTA = DELTA*(BNRM/BNRM0)
 130  55   TEM = 1.0D0/BNRM
 131       CALL DSCAL (N, TEM, V(1,1), 1)
 132 C Zero out the HES array. ----------------------------------------------
 133       DO 65 J = 1,MAXL
 134         DO 60 I = 1,MAXLP1
 135  60       HES(I,J) = 0.0D0
 136  65     CONTINUE
 137 C-----------------------------------------------------------------------
 138 C Main loop to compute the vectors V(*,2) to V(*,MAXL).
 139 C The running product PROD is needed for the convergence test.
 140 C-----------------------------------------------------------------------
 141       PROD = 1.0D0
 142       DO 90 LL = 1,MAXL
 143         LGMR = LL
 144 C-----------------------------------------------------------------------
 145 C Call routine DATV to compute VNEW = Abar*v(ll), where Abar is
 146 C the matrix A with scaling and inverse preconditioner factors applied.
 147 C Call routine DORTHOG to orthogonalize the new vector VNEW = V(*,LL+1).
 148 C Call routine DHEQR to update the factors of HES.
 149 C-----------------------------------------------------------------------
 150         CALL DATV (NEQ, Y, SAVF, V(1,LL), WGHT, X, F, PSOL, V(1,LL+1),
 151      1        WK, WP, IWP, HL0, JPRE, IER, NPSL)
 152         IF (IER .NE. 0) GO TO 300
 153         CALL DORTHOG (V(1,LL+1), V, HES, N, LL, MAXLP1, KMP, SNORMW)
 154         HES(LL+1,LL) = SNORMW
 155         CALL DHEQR (HES, MAXLP1, LL, Q, INFO, LL)
 156         IF (INFO .EQ. LL) GO TO 120
 157 C-----------------------------------------------------------------------
 158 C Update RHO, the estimate of the norm of the residual b - A*xl.
 159 C If KMP .lt. MAXL, then the vectors V(*,1),...,V(*,LL+1) are not
 160 C necessarily orthogonal for LL .gt. KMP.  The vector DL must then
 161 C be computed, and its norm used in the calculation of RHO.
 162 C-----------------------------------------------------------------------
 163         PROD = PROD*Q(2*LL)
 164         RHO = ABS(PROD*BNRM)
 165         IF ((LL.GT.KMP) .AND. (KMP.LT.MAXL)) THEN
 166           IF (LL .EQ. KMP+1) THEN
 167             CALL DCOPY (N, V(1,1), 1, DL, 1)
 168             DO 75 I = 1,KMP
 169               IP1 = I + 1
 170               I2 = I*2
 171               S = Q(I2)
 172               C = Q(I2-1)
 173               DO 70 K = 1,N
 174  70             DL(K) = S*DL(K) + C*V(K,IP1)
 175  75           CONTINUE
 176             ENDIF
 177           S = Q(2*LL)
 178           C = Q(2*LL-1)/SNORMW
 179           LLP1 = LL + 1
 180           DO 80 K = 1,N
 181  80         DL(K) = S*DL(K) + C*V(K,LLP1)
 182           DLNRM = DNRM2 (N, DL, 1)
 183           RHO = RHO*DLNRM
 184           ENDIF
 185 C-----------------------------------------------------------------------
 186 C Test for convergence.  If passed, compute approximation xl.
 187 C if failed and LL .lt. MAXL, then continue iterating.
 188 C-----------------------------------------------------------------------
 189         IF (RHO .LE. DELTA) GO TO 200
 190         IF (LL .EQ. MAXL) GO TO 100
 191 C-----------------------------------------------------------------------
 192 C Rescale so that the norm of V(1,LL+1) is one.
 193 C-----------------------------------------------------------------------
 194         TEM = 1.0D0/SNORMW
 195         CALL DSCAL (N, TEM, V(1,LL+1), 1)
 196  90     CONTINUE
 197  100  CONTINUE
 198       IF (RHO .LE. 1.0D0) GO TO 150
 199       IF (RHO .LE. BNRM .AND. MNEWT .EQ. 0) GO TO 150
 200  120  CONTINUE
 201       IFLAG = 2
 202       RETURN
 203  150  IFLAG = 1
 204 C-----------------------------------------------------------------------
 205 C Compute the approximation xl to the solution.
 206 C Since the vector X was used as work space, and the initial guess
 207 C of the Newton correction is zero, X must be reset to zero.
 208 C-----------------------------------------------------------------------
 209  200  CONTINUE
 210       LL = LGMR
 211       LLP1 = LL + 1
 212       DO 210 K = 1,LLP1
 213  210    B(K) = 0.0D0
 214       B(1) = BNRM
 215       CALL DHELS (HES, MAXLP1, LL, Q, B)
 216       DO 220 K = 1,N
 217  220    X(K) = 0.0D0
 218       DO 230 I = 1,LL
 219         CALL DAXPY (N, B(I), V(1,I), 1, X, 1)
 220  230    CONTINUE
 221       DO 240 I = 1,N
 222  240    X(I) = X(I)/WGHT(I)
 223       IF (JPRE .LE. 1) RETURN
 224       CALL PSOL (NEQ, TN, Y, SAVF, WK, HL0, WP, IWP, X, 2, IER)
 225       NPSL = NPSL + 1
 226       IF (IER .NE. 0) GO TO 300
 227       RETURN
 228 C-----------------------------------------------------------------------
 229 C This block handles error returns forced by routine PSOL.
 230 C-----------------------------------------------------------------------
 231  300  CONTINUE
 232       IF (IER .LT. 0) IFLAG = -1
 233       IF (IER .GT. 0) IFLAG = 3
 234 C
 235       RETURN
 236 C----------------------- End of Subroutine DSPIGMR ---------------------
 237       END