share/odepack/fortran/dpcg.f

   1 *DECK DPCG
   2       SUBROUTINE DPCG (NEQ, TN, Y, SAVF, R, WGHT, N, MAXL, DELTA, HL0,
   3      1 JPRE, MNEWT, F, PSOL, NPSL, X, P, W, Z, LPCG, WP, IWP, WK, IFLAG)
   4       EXTERNAL F, PSOL
   5       INTEGER NEQ, N, MAXL, JPRE, MNEWT, NPSL, LPCG, IWP, IFLAG
   6       DOUBLE PRECISION TN,Y,SAVF,R,WGHT,DELTA,HL0,X,P,W,Z,WP,WK
   7       DIMENSION NEQ(*), Y(*), SAVF(*), R(*), WGHT(*), X(*), P(*), W(*),
   8      1   Z(*), WP(*), IWP(*), WK(*)
   9 C-----------------------------------------------------------------------
  10 C This routine computes the solution to the system A*x = b using a
  11 C preconditioned version of the Conjugate Gradient algorithm.
  12 C It is assumed here that the matrix A and the preconditioner
  13 C matrix M are symmetric positive definite or nearly so.
  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            R = the right hand side of the system A*x = b.
  27 C
  28 C         WGHT = array of length N containing scale factors.
  29 C                1/WGHT(i) are the diagonal elements of the diagonal
  30 C                scaling matrix D.
  31 C
  32 C            N = the order of the matrix A, and the lengths
  33 C                of the vectors Y, SAVF, R, WGHT, P, W, Z, WK, and X.
  34 C
  35 C         MAXL = the maximum allowable number of iterates.
  36 C
  37 C        DELTA = tolerance on residuals b - A*x in weighted RMS-norm.
  38 C
  39 C          HL0 = current value of (step size h) * (coefficient l0).
  40 C
  41 C         JPRE = preconditioner type flag.
  42 C
  43 C        MNEWT = Newton iteration counter (.ge. 0).
  44 C
  45 C           WK = real work array used by routine DATP.
  46 C
  47 C           WP = real work array used by preconditioner PSOL.
  48 C
  49 C          IWP = integer work array used by preconditioner PSOL.
  50 C
  51 C      On return
  52 C
  53 C         X    = the final computed approximation to the solution
  54 C                of the system A*x = b.
  55 C
  56 C         LPCG = the number of iterations performed, and current
  57 C                order of the upper Hessenberg matrix HES.
  58 C
  59 C         NPSL = the number of calls to PSOL.
  60 C
  61 C        IFLAG = integer error flag:
  62 C                0 means convergence in LPCG iterations, LPCG .le. MAXL.
  63 C                1 means the convergence test did not pass in MAXL
  64 C                  iterations, but the residual norm is .lt. 1,
  65 C                  or .lt. norm(b) if MNEWT = 0, and so X is computed.
  66 C                2 means the convergence test did not pass in MAXL
  67 C                  iterations, residual .gt. 1, and X is undefined.
  68 C                3 means there was a recoverable error in PSOL
  69 C                  caused by the preconditioner being out of date.
  70 C                4 means there was a zero denominator in the algorithm.
  71 C                  The system matrix or preconditioner matrix is not
  72 C                  sufficiently close to being symmetric pos. definite.
  73 C               -1 means there was a nonrecoverable error in PSOL.
  74 C
  75 C-----------------------------------------------------------------------
  76       INTEGER I, IER
  77       DOUBLE PRECISION ALPHA,BETA,BNRM,PTW,RNRM,DDOT,DVNORM,ZTR,ZTR0
  78 C
  79       IFLAG = 0
  80       NPSL = 0
  81       LPCG = 0
  82       DO 10 I = 1,N
  83  10     X(I) = 0.0D0
  84       BNRM = DVNORM (N, R, WGHT)
  85 C Test for immediate return with X = 0 or X = b. -----------------------
  86       IF (BNRM .GT. DELTA) GO TO 20
  87       IF (MNEWT .GT. 0) RETURN
  88       CALL DCOPY (N, R, 1, X, 1)
  89       RETURN
  90 C
  91  20   ZTR = 0.0D0
  92 C Loop point for PCG iterations. ---------------------------------------
  93  30   CONTINUE
  94       LPCG = LPCG + 1
  95       CALL DCOPY (N, R, 1, Z, 1)
  96       IER = 0
  97       IF (JPRE .EQ. 0) GO TO 40
  98       CALL PSOL (NEQ, TN, Y, SAVF, WK, HL0, WP, IWP, Z, 3, IER)
  99       NPSL = NPSL + 1
 100       IF (IER .NE. 0) GO TO 100
 101  40   CONTINUE
 102       ZTR0 = ZTR
 103       ZTR = DDOT (N, Z, 1, R, 1)
 104       IF (LPCG .NE. 1) GO TO 50
 105       CALL DCOPY (N, Z, 1, P, 1)
 106       GO TO 70
 107  50   CONTINUE
 108       IF (ZTR0 .EQ. 0.0D0) GO TO 200
 109       BETA = ZTR/ZTR0
 110       DO 60 I = 1,N
 111  60     P(I) = Z(I) + BETA*P(I)
 112  70   CONTINUE
 113 C-----------------------------------------------------------------------
 114 C  Call DATP to compute A*p and return the answer in W.
 115 C-----------------------------------------------------------------------
 116       CALL DATP (NEQ, Y, SAVF, P, WGHT, HL0, WK, F, W)
 117 C
 118       PTW = DDOT (N, P, 1, W, 1)
 119       IF (PTW .EQ. 0.0D0) GO TO 200
 120       ALPHA = ZTR/PTW
 121       CALL DAXPY (N, ALPHA, P, 1, X, 1)
 122       ALPHA = -ALPHA
 123       CALL DAXPY (N, ALPHA, W, 1, R, 1)
 124       RNRM = DVNORM (N, R, WGHT)
 125       IF (RNRM .LE. DELTA) RETURN
 126       IF (LPCG .LT. MAXL) GO TO 30
 127       IFLAG = 2
 128       IF (RNRM .LE. 1.0D0) IFLAG = 1
 129       IF (RNRM .LE. BNRM .AND. MNEWT .EQ. 0) IFLAG = 1
 130       RETURN
 131 C-----------------------------------------------------------------------
 132 C This block handles error returns from PSOL.
 133 C-----------------------------------------------------------------------
 134  100  CONTINUE
 135       IF (IER .LT. 0) IFLAG = -1
 136       IF (IER .GT. 0) IFLAG = 3
 137       RETURN
 138 C-----------------------------------------------------------------------
 139 C This block handles division by zero errors.
 140 C-----------------------------------------------------------------------
 141  200  CONTINUE
 142       IFLAG = 4
 143       RETURN
 144 C----------------------- End of Subroutine DPCG ------------------------
 145       END