share/odepack/fortran/dspiom.f

   1 *DECK DSPIOM
   2       SUBROUTINE DSPIOM (NEQ, TN, Y, SAVF, B, WGHT, N, MAXL, KMP, DELTA,
   3      1            HL0, JPRE, MNEWT, F, PSOL, NPSL, X, V, HES, IPVT,
   4      2            LIOM, WP, IWP, WK, IFLAG)
   5       EXTERNAL F, PSOL
   6       INTEGER NEQ,N,MAXL,KMP,JPRE,MNEWT,NPSL,IPVT,LIOM,IWP,IFLAG
   7       DOUBLE PRECISION TN,Y,SAVF,B,WGHT,DELTA,HL0,X,V,HES,WP,WK
   8       DIMENSION NEQ(*), Y(*), SAVF(*), B(*), WGHT(*), X(*), V(N,*),
   9      1   HES(MAXL,MAXL), IPVT(*), WP(*), IWP(*), WK(*)
  10 C-----------------------------------------------------------------------
  11 C This routine solves the linear system A * x = b using a scaled
  12 C preconditioned version of the Incomplete Orthogonalization 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 the
  28 C                final approximation.
  29 C                (B is the same as V(*,MAXL+1) in the call to DSPIOM.)
  30 C
  31 C         WGHT = array of length N containing scale factors.
  32 C                1/WGHT(i) are the diagonal elements of the diagonal
  33 C                scaling matrix D.
  34 C
  35 C         N    = the order of the matrix A, and the lengths
  36 C                of the vectors Y, SAVF, B, WGHT, and X.
  37 C
  38 C         MAXL = the maximum allowable order of the matrix HES.
  39 C
  40 C          KMP = the number of previous vectors the new vector VNEW
  41 C                must be made orthogonal to.  KMP .le. MAXL.
  42 C
  43 C        DELTA = tolerance on residuals b - A*x in weighted RMS-norm.
  44 C
  45 C          HL0 = current value of (step size h) * (coefficient l0).
  46 C
  47 C         JPRE = preconditioner type flag.
  48 C
  49 C        MNEWT = Newton iteration counter (.ge. 0).
  50 C
  51 C           WK = real work array of length N used by DATV and PSOL.
  52 C
  53 C           WP = real work array used by preconditioner PSOL.
  54 C
  55 C          IWP = integer work array used by preconditioner PSOL.
  56 C
  57 C      On return
  58 C
  59 C         X    = the final computed approximation to the solution
  60 C                of the system A*x = b.
  61 C
  62 C         V    = the N by (LIOM+1) array containing the LIOM
  63 C                orthogonal vectors V(*,1) to V(*,LIOM).
  64 C
  65 C         HES  = the LU factorization of the LIOM by LIOM upper
  66 C                Hessenberg matrix whose entries are the
  67 C                scaled inner products of A*V(*,k) and V(*,i).
  68 C
  69 C         IPVT = an integer array containg pivoting information.
  70 C                It is loaded in DHEFA and used in DHESL.
  71 C
  72 C         LIOM = the number of iterations performed, and current
  73 C                order of the upper Hessenberg matrix HES.
  74 C
  75 C         NPSL = the number of calls to PSOL.
  76 C
  77 C        IFLAG = integer error flag:
  78 C                0 means convergence in LIOM iterations, LIOM.le.MAXL.
  79 C                1 means the convergence test did not pass in MAXL
  80 C                  iterations, but the residual norm is .lt. 1,
  81 C                  or .lt. norm(b) if MNEWT = 0, and so X is computed.
  82 C                2 means the convergence test did not pass in MAXL
  83 C                  iterations, residual .gt. 1, and X is undefined.
  84 C                3 means there was a recoverable error in PSOL
  85 C                  caused by the preconditioner being out of date.
  86 C               -1 means there was a nonrecoverable error in PSOL.
  87 C
  88 C-----------------------------------------------------------------------
  89       INTEGER I, IER, INFO, J, K, LL, LM1
  90       DOUBLE PRECISION BNRM, BNRM0, PROD, RHO, SNORMW, DNRM2, TEM
  91 C
  92       IFLAG = 0
  93       LIOM = 0
  94       NPSL = 0
  95 C-----------------------------------------------------------------------
  96 C The initial residual is the vector b.  Apply scaling to b, and test
  97 C for an immediate return with X = 0 or X = b.
  98 C-----------------------------------------------------------------------
  99       DO 10 I = 1,N
 100  10     V(I,1) = B(I)*WGHT(I)
 101       BNRM0 = DNRM2 (N, V, 1)
 102       BNRM = BNRM0
 103       IF (BNRM0 .GT. DELTA) GO TO 30
 104       IF (MNEWT .GT. 0) GO TO 20
 105       CALL DCOPY (N, B, 1, X, 1)
 106       RETURN
 107  20   DO 25 I = 1,N
 108  25     X(I) = 0.0D0
 109       RETURN
 110  30   CONTINUE
 111 C Apply inverse of left preconditioner to vector b. --------------------
 112       IER = 0
 113       IF (JPRE .EQ. 0 .OR. JPRE .EQ. 2) GO TO 55
 114       CALL PSOL (NEQ, TN, Y, SAVF, WK, HL0, WP, IWP, B, 1, IER)
 115       NPSL = 1
 116       IF (IER .NE. 0) GO TO 300
 117 C Calculate norm of scaled vector V(*,1) and normalize it. -------------
 118       DO 50 I = 1,N
 119  50     V(I,1) = B(I)*WGHT(I)
 120       BNRM = DNRM2(N, V, 1)
 121       DELTA = DELTA*(BNRM/BNRM0)
 122  55   TEM = 1.0D0/BNRM
 123       CALL DSCAL (N, TEM, V(1,1), 1)
 124 C Zero out the HES array. ----------------------------------------------
 125       DO 65 J = 1,MAXL
 126         DO 60 I = 1,MAXL
 127  60       HES(I,J) = 0.0D0
 128  65     CONTINUE
 129 C-----------------------------------------------------------------------
 130 C Main loop on LL = l to compute the vectors V(*,2) to V(*,MAXL).
 131 C The running product PROD is needed for the convergence test.
 132 C-----------------------------------------------------------------------
 133       PROD = 1.0D0
 134       DO 90 LL = 1,MAXL
 135         LIOM = LL
 136 C-----------------------------------------------------------------------
 137 C Call routine DATV to compute VNEW = Abar*v(l), where Abar is
 138 C the matrix A with scaling and inverse preconditioner factors applied.
 139 C Call routine DORTHOG to orthogonalize the new vector vnew = V(*,l+1).
 140 C Call routine DHEFA to update the factors of HES.
 141 C-----------------------------------------------------------------------
 142         CALL DATV (NEQ, Y, SAVF, V(1,LL), WGHT, X, F, PSOL, V(1,LL+1),
 143      1        WK, WP, IWP, HL0, JPRE, IER, NPSL)
 144         IF (IER .NE. 0) GO TO 300
 145         CALL DORTHOG (V(1,LL+1), V, HES, N, LL, MAXL, KMP, SNORMW)
 146         CALL DHEFA (HES, MAXL, LL, IPVT, INFO, LL)
 147         LM1 = LL - 1
 148         IF (LL .GT. 1 .AND. IPVT(LM1) .EQ. LM1) PROD = PROD*HES(LL,LM1)
 149         IF (INFO .NE. LL) GO TO 70
 150 C-----------------------------------------------------------------------
 151 C The last pivot in HES was found to be zero.
 152 C If vnew = 0 or l = MAXL, take an error return with IFLAG = 2.
 153 C otherwise, continue the iteration without a convergence test.
 154 C-----------------------------------------------------------------------
 155         IF (SNORMW .EQ. 0.0D0) GO TO 120
 156         IF (LL .EQ. MAXL) GO TO 120
 157         GO TO 80
 158 C-----------------------------------------------------------------------
 159 C Update RHO, the estimate of the norm of the residual b - A*x(l).
 160 C test for convergence.  If passed, compute approximation x(l).
 161 C If failed and l .lt. MAXL, then continue iterating.
 162 C-----------------------------------------------------------------------
 163  70     CONTINUE
 164         RHO = BNRM*SNORMW*ABS(PROD/HES(LL,LL))
 165         IF (RHO .LE. DELTA) GO TO 200
 166         IF (LL .EQ. MAXL) GO TO 100
 167 C If l .lt. MAXL, store HES(l+1,l) and normalize the vector v(*,l+1).
 168  80     CONTINUE
 169         HES(LL+1,LL) = SNORMW
 170         TEM = 1.0D0/SNORMW
 171         CALL DSCAL (N, TEM, V(1,LL+1), 1)
 172  90     CONTINUE
 173 C-----------------------------------------------------------------------
 174 C l has reached MAXL without passing the convergence test:
 175 C If RHO is not too large, compute a solution anyway and return with
 176 C IFLAG = 1.  Otherwise return with IFLAG = 2.
 177 C-----------------------------------------------------------------------
 178  100  CONTINUE
 179       IF (RHO .LE. 1.0D0) GO TO 150
 180       IF (RHO .LE. BNRM .AND. MNEWT .EQ. 0) GO TO 150
 181  120  CONTINUE
 182       IFLAG = 2
 183       RETURN
 184  150  IFLAG = 1
 185 C-----------------------------------------------------------------------
 186 C Compute the approximation x(l) to the solution.
 187 C Since the vector X was used as work space, and the initial guess
 188 C of the Newton correction is zero, X must be reset to zero.
 189 C-----------------------------------------------------------------------
 190  200  CONTINUE
 191       LL = LIOM
 192       DO 210 K = 1,LL
 193  210    B(K) = 0.0D0
 194       B(1) = BNRM
 195       CALL DHESL (HES, MAXL, LL, IPVT, B)
 196       DO 220 K = 1,N
 197  220    X(K) = 0.0D0
 198       DO 230 I = 1,LL
 199         CALL DAXPY (N, B(I), V(1,I), 1, X, 1)
 200  230    CONTINUE
 201       DO 240 I = 1,N
 202  240    X(I) = X(I)/WGHT(I)
 203       IF (JPRE .LE. 1) RETURN
 204       CALL PSOL (NEQ, TN, Y, SAVF, WK, HL0, WP, IWP, X, 2, IER)
 205       NPSL = NPSL + 1
 206       IF (IER .NE. 0) GO TO 300
 207       RETURN
 208 C-----------------------------------------------------------------------
 209 C This block handles error returns forced by routine PSOL.
 210 C-----------------------------------------------------------------------
 211  300  CONTINUE
 212       IF (IER .LT. 0) IFLAG = -1
 213       IF (IER .GT. 0) IFLAG = 3
 214       RETURN
 215 C----------------------- End of Subroutine DSPIOM ----------------------
 216       END