      SUBROUTINE DTD3M(N,NZ,IR,IC,A,B,JOB,SPD)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION IR(NZ),IC(NZ),A(NZ),B(N)
      LOGICAL SPD
      COMMON /DTDP27/ ITASK,NPES,ICOMM
C
C  DTD3M IS A LOCALLY-WRITTEN CODE WHICH SOLVES A SPARSE SYMMETRIC
C  LINEAR SYSTEM A*X=B.  IT WILL BE CALLED WHEN ISOLVE=5 FOR A SYMMETRIC
C  2D GALERKIN PROBLEM, AND WHEN ISOLVE=4 FOR A 2D OR 3D COLLOCATION PROBLEM.  
C  IF YOU ACTIVATE DTD3M,DTD3N, YOU SHOULD ALSO INCREASE ISMX2D AND ISMX3D TO 
C  AT LEAST 5 IN THE FILE 'pde2d.f'.
C
C  N  - NUMBER OF EQUATIONS AND UNKNOWNS (INPUT)
C  NZ - NUMBER OF NONZEROS IN THE UPPER TRIANGLE OF A (INPUT)
C  IR - ROW NUMBERS OF THE NONZEROS IN THE UPPER TRIANGLE OF A (INPUT)
C  IC - COLUMN NUMBERS OF THE NONZEROS IN THE UPPER TRIANGLE OF A (INPUT)
C  A  - NONZERO ELEMENTS OF THE UPPER TRIANGLE OF A (INPUT).  A(J)
C         CONTAINS ELEMENT (IR(J),IC(J)), J=1,...,NZ, OF THE MATRIX,
C         WHERE IC(J).GE.IR(J).
C  B  - ON INPUT, B WILL CONTAIN THE RIGHT HAND SIDE OF THE LINEAR
C         SYSTEM.  ON OUTPUT, B SHOULD CONTAIN THE SOLUTION, X.
C  JOB- JOB PARAMETER (INPUT).  IF JOB=2, THIS MEANS THAT THE MATRIX
C         A HAS CHANGED SINCE THE LAST CALL TO DTD3M, WHILE JOB=3 MEANS
C         A HAS NOT CHANGED.  THUS, IF YOU WISH, YOU CAN COMPUTE AN LU
C         DECOMPOSITION OF A WHEN JOB=2 AND SAVE IT, AND USE THIS
C         DECOMPOSITION TO SOLVE THE SYSTEM MORE RAPIDLY, WHEN JOB=3.
C  SPD- .TRUE. IF LINEAR SYSTEM IS POSITIVE DEFINITE. (INPUT)
C
C  IF MORE THAN ONE PROCESSOR IS USED (NPES > 1), THE MATRIX WILL BE
C  DISTRIBUTED OVER THE PROCESSORS, AND ONLY THOSE NONZEROS WITH
C  MOD(IC(J)-1,NPES) = ITASK WILL BE PASSED TO DTD3M ON PROCESSOR ITASK.
C
C      CALL DTDPS (3,
C     & ' Local solver not available.  Try another value for ISOLVE.',
C     & 0.0D0,0)
      DIMENSION X(N)
      CALL DCG(A,IR,IC,NZ,X,B,N,.TRUE.) 
C          DTD3M EXPECTS THE SOLUTION TO BE RETURNED IN B.  
      B = X 
      RETURN
      END


      SUBROUTINE DTD3N(N,NZ,IR,IC,A,B,JOB)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION IR(NZ),IC(NZ),A(NZ),B(N)
      COMMON /DTDP27/ ITASK,NPES,ICOMM
C
C  DTD3N IS A LOCALLY-WRITTEN CODE WHICH SOLVES A SPARSE NONSYMMETRIC
C  LINEAR SYSTEM A*X=B.  IT WILL BE CALLED WHEN ISOLVE=5 FOR A NONSYMMETRIC 
C  2D GALERKIN PROBLEM, AND WHEN ISOLVE=5 FOR A 2D OR 3D COLLOCATION PROBLEM.  
C  IF YOU ACTIVATE DTD3M,DTD3N, YOU SHOULD ALSO INCREASE ISMX2D AND ISMX3D
C  TO AT LEAST 5 IN THE FILE 'pde2d.f'.
C
C  N  - NUMBER OF EQUATIONS AND UNKNOWNS (INPUT)
C  NZ - NUMBER OF NONZEROS IN A (INPUT)
C  IR - ROW NUMBERS OF THE NONZEROS OF A (INPUT)
C  IC - COLUMN NUMBERS OF THE NONZEROS OF A (INPUT)
C  A  - NONZERO ELEMENTS OF A (INPUT).  A(J) CONTAINS ELEMENT
C         (IR(J),IC(J)), J=1,...,NZ, OF THE MATRIX.
C  B  - ON INPUT, B WILL CONTAIN THE RIGHT HAND SIDE OF THE LINEAR
C         SYSTEM.  ON OUTPUT, B SHOULD CONTAIN THE SOLUTION, X.
C  JOB- JOB PARAMETER (INPUT).  IF JOB=2, THIS MEANS THAT THE MATRIX
C         A HAS CHANGED SINCE THE LAST CALL TO DTD3N, WHILE JOB=3 MEANS
C         A HAS NOT CHANGED.  THUS, IF YOU WISH, YOU CAN COMPUTE AN LU
C         DECOMPOSITION OF A WHEN JOB=2 AND SAVE IT, AND USE THIS
C         DECOMPOSITION TO SOLVE THE SYSTEM MORE RAPIDLY, WHEN JOB=3.
C
C  IF MORE THAN ONE PROCESSOR IS USED (NPES > 1), THE MATRIX WILL BE
C  DISTRIBUTED OVER THE PROCESSORS, AND ONLY THOSE NONZEROS WITH
C  MOD(IC(J)-1,NPES) = ITASK WILL BE PASSED TO DTD3N ON PROCESSOR ITASK.
C
C     CALL DTDPS (3,
C    & ' Local solver not available.  Try another value for ISOLVE.',
C    & 0.0D0,0)
      DIMENSION X(N)
      CALL DCG(A,IR,IC,NZ,X,B,N,.FALSE.)
C          DTD3N EXPECTS THE SOLUTION TO BE RETURNED IN B.  
      B = X
      RETURN
      END


      SUBROUTINE DCG(A,IROW,JCOL,NZ,X,B,N,SYMM)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C                                  DECLARATIONS FOR ARGUMENTS
      DOUBLE PRECISION A(NZ),B(N),X(N)
      INTEGER IROW(NZ),JCOL(NZ)
      LOGICAL SYMM
C                                  DECLARATIONS FOR LOCAL VARIABLES
      DOUBLE PRECISION R(N),P(N),AP(N),D(N),LAMBDA
C
C  SUBROUTINE DCG SOLVES THE LINEAR SYSTEM A*X=B, USING THE JACOBI
C     CONJUGATE GRADIENT ITERATIVE METHOD.  THE NON-ZEROS OF A ARE 
C     STORED 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    SYMM   - .TRUE. IF A IS SYMMETRIC, IN
C             WHICH CASE ONLY THE NONZEROS
C             IN THE UPPER TRIANGLE ARE INPUT
C-----------------------------------------------------------------------
C
      DO 5 I=1,NZ
C                                  DCG REQUIRES THE DIAGONAL OF A, FOR 
C                                  "JACOBI" PRECONDITIONING 
         IF (IROW(I).EQ.JCOL(I)) D(IROW(I)) = A(I)
    5 CONTINUE 
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)
            IF (I.NE.J .AND. SYMM) AP(J) = AP(J) + A(IZ)*P(I)
   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
