C/ FIGURE 1.10.1
      SUBROUTINE DCG(A,IROW,JCOL,NZ,X,B,N,D)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C                                  DECLARATIONS FOR ARGUMENTS
      DOUBLE PRECISION A(NZ),B(N),X(N),D(N)
      INTEGER IROW(NZ),JCOL(NZ)
C                                  DECLARATIONS FOR LOCAL VARIABLES
      DOUBLE PRECISION R(N),P(N),AP(N),LAMBDA
C
C  SUBROUTINE DCG SOLVES THE SYMMETRIC LINEAR SYSTEM A*X=B, USING THE
C     CONJUGATE GRADIENT ITERATIVE METHOD.  THE NON-ZEROS OF A ARE STORED
C     IN SPARSE FORMAT.  
C
C  ARGUMENTS
C
C             ON INPUT                          ON OUTPUT
C             --------                          ---------
C
C    A      - A(IZ) IS THE MATRIX ELEMENT IN
C             ROW IROW(IZ), COLUMN JCOL(IZ),
C             FOR IZ=1,...,NZ. 
C
C    IROW   - (SEE A).
C
C    JCOL   - (SEE A).
C
C    NZ     - NUMBER OF NONZEROS.
C
C    X      -                                   AN N-VECTOR CONTAINING
C                                               THE SOLUTION.
C
C    B      - THE RIGHT HAND SIDE N-VECTOR.
C
C    N      - SIZE OF MATRIX A.
C
C    D      - VECTOR HOLDING A DIAGONAL 
C             PRECONDITIONING MATRIX.
C             D = DIAGONAL(A) IS RECOMMENDED.
C             D(I) = 1 FOR NO PRECONDITIONING
C
C-----------------------------------------------------------------------
C                                  X0 = 0
C                                  R0 = D**(-1)*B
C                                  P0 = R0
      R0MAX = 0
      DO 10 I=1,N
         X(I) = 0
         R(I) = B(I)/D(I)
         R0MAX = MAX(R0MAX,ABS(R(I)))
         P(I) = R(I)
   10 CONTINUE
C                                  NITER = MAX NUMBER OF ITERATIONS
      NITER = 3*N 
      DO 90 ITER=1,NITER
C                                  AP = A*P
         DO 20 I=1,N
            AP(I) = 0
   20    CONTINUE
         DO 30 IZ=1,NZ
            I = IROW(IZ)
            J = JCOL(IZ)
            AP(I) = AP(I) + A(IZ)*P(J)
   30    CONTINUE
C                                  PAP = (P,AP)
C                                  RP = (R,D*P)
         PAP = 0.0
         RP = 0.0
         DO 40 I=1,N
            PAP = PAP + P(I)*AP(I)
            RP = RP + R(I)*D(I)*P(I)
   40    CONTINUE
C                                  LAMBDA = (R,D*P)/(P,AP)
         LAMBDA = RP/PAP
C                                  X = X + LAMBDA*P
C                                  R = R - LAMBDA*D**(-1)*AP
         DO 50 I=1,N
            X(I) = X(I) + LAMBDA*P(I)
            R(I) = R(I) - LAMBDA*AP(I)/D(I)
   50    CONTINUE
C                                  RAP = (R,AP)
         RAP = 0.0
         DO 60 I=1,N
            RAP = RAP + R(I)*AP(I)
   60    CONTINUE
C                                  ALPHA = -(R,AP)/(P,AP)
         ALPHA = -RAP/PAP
C                                  P = R + ALPHA*P
         DO 70 I=1,N
            P(I) = R(I) + ALPHA*P(I)
   70    CONTINUE
C                                  RMAX = MAX OF RESIDUAL (R)
         RMAX = 0
         DO 80 I=1,N
            RMAX = MAX(RMAX,ABS(R(I)))
   80    CONTINUE
C                                  CHECK IF CONVERGED
         IF (RMAX.LE.1.D-10*R0MAX) THEN
            PRINT *, ' Number of iterations = ',ITER
            RETURN
         ENDIF
   90 CONTINUE
C                                  DCG DOES NOT CONVERGE
      PRINT 100
  100 FORMAT('***** DCG DOES NOT CONVERGE *****')
      RETURN
      END
