% FIGURE 1.10.1
      function X = DCG(A,IROW,JCOL,NZ,B,N,D)
%
%  FUNCTION DCG SOLVES THE SYMMETRIC LINEAR SYSTEM A*X=B, USING THE
%     CONJUGATE GRADIENT ITERATIVE METHOD.  THE NON-ZEROS OF A ARE STORED
%     IN SPARSE FORMAT.  
%
%  ARGUMENTS
%
%             ON INPUT                          ON OUTPUT
%             --------                          ---------
%
%    A      - A(IZ) IS THE MATRIX ELEMENT IN
%             ROW IROW(IZ), COLUMN JCOL(IZ),
%             FOR IZ=1,...,NZ. 
%
%    IROW   - (SEE A).
%
%    JCOL   - (SEE A).
%
%    NZ     - NUMBER OF NONZEROS.
%
%    X      -                                   AN N-VECTOR CONTAINING
%                                               THE SOLUTION.
%
%    B      - THE RIGHT HAND SIDE N-VECTOR.
%
%    N      - SIZE OF MATRIX A.
%
%    D      - VECTOR HOLDING A DIAGONAL 
%             PRECONDITIONING MATRIX.
%             D = DIAGONAL(A) IS RECOMMENDED.
%             D(I) = 1 FOR NO PRECONDITIONING
%
%-----------------------------------------------------------------------
%                                  X0 = 0
%                                  R0 = D^(-1)*B
%                                  P0 = R0
      R0MAX = 0;
      for I=1:N
         X(I) = 0;
         R(I) = B(I)/D(I);
         R0MAX = max(R0MAX,abs(R(I)));
         P(I) = R(I);
      end
%                                  NITER = MAX NUMBER OF ITERATIONS
      NITER = 3*N;
      for ITER=1:NITER
%                                  AP = A*P
         for I=1:N
            AP(I) = 0;
         end
         for IZ=1:NZ
            I = IROW(IZ);
            J = JCOL(IZ);
            AP(I) = AP(I) + A(IZ)*P(J);
         end
%                                  PAP = (P,AP)
%                                  RP = (R,D*P)
         PAP = 0.0;
         RP = 0.0;
         for I=1:N
            PAP = PAP + P(I)*AP(I);
            RP = RP + R(I)*D(I)*P(I);
         end
%                                  LAMBDA = (R,D*P)/(P,AP)
         LAMBDA = RP/PAP;
%                                  X = X + LAMBDA*P
%                                  R = R - LAMBDA*D^(-1)*AP
         for I=1:N
            X(I) = X(I) + LAMBDA*P(I);
            R(I) = R(I) - LAMBDA*AP(I)/D(I);
         end
%                                  RAP = (R,AP)
         RAP = 0.0;
         for I=1:N
            RAP = RAP + R(I)*AP(I);
         end
%                                  ALPHA = -(R,AP)/(P,AP)
         ALPHA = -RAP/PAP;
%                                  P = R + ALPHA*P
         for I=1:N
            P(I) = R(I) + ALPHA*P(I);
         end
%                                  RMAX = MAX OF RESIDUAL (R)
         RMAX = 0;
         for I=1:N
            RMAX = max(RMAX,abs(R(I)));
         end
%                                  CHECK FOR CONVERGENCE
         if (RMAX <= 1.e-10*R0MAX) 
            disp ([' Number of iterations = ',num2str(ITER)])
            return
         end
      end
%                                  DCG DOES NOT CONVERGE
      error ('***** DCG DOES NOT CONVERGE *****')
      end
