      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 PCG(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)
      include 'mpif.h'
      ALLOCATABLE ABAND(:,:)
      NLDI = 0
      NUDI = 0
      DO 10 IZ=1,NZ
         NLDI =  MAX(IR(IZ)-IC(IZ),NLDI)
         NUDI =  MAX(IC(IZ)-IR(IZ),NUDI)
   10 CONTINUE
C          NLD = NUMBER OF LOWER DIAGONALS, NUD = NUMBER OF UPPER DIAGONALS
      CALL MPI_ALLREDUCE(NLDI,NLD,1,MPI_INTEGER,MPI_MAX,ICOMM,IERR)
      CALL MPI_ALLREDUCE(NUDI,NUD,1,MPI_INTEGER,MPI_MAX,ICOMM,IERR)
      ALLOCATE (ABAND(-NUD-NLD:NLD,(N-1)/NPES+1))
C          COPY LOCAL PROCESSOR'S COLUMNS OF SPARSE MATRIX A TO 
C          LOCAL PROCESSOR'S COLUMNS OF BAND MATRIX ABAND
      ABAND = 0
      DO 20 IZ=1,NZ
         I = IR(IZ)
         J = IC(IZ)
         ABAND(I-J,(J-1)/NPES+1) = A(IZ)
   20 CONTINUE
      CALL PBAND(ABAND,N,NLD,NUD,X,B)
C          DTD3N EXPECTS THE SOLUTION TO BE RETURNED IN B.  
      B = X
      RETURN
      END


      SUBROUTINE PCG(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),PI(N),API(N),AP(N),DI(N),D(N),LAMBDA
      include 'mpif.h'
C
C  SUBROUTINE PCG 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.  THE COLUMNS OF A ARE DISTRIBUTED 
C     CYCLICALLY OVER THE AVAILABLE PROCESSORS.
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. ELEMENTS WITH
C             MOD(JCOL(IZ)-1,NPES)=ITASK
C             ARE STORED ON PROCESSOR ITASK.
C
C    IROW   - (SEE A).
C
C    JCOL   - (SEE A).
C
C    NZ     - NUMBER OF NONZEROS STORED ON
C             THE LOCAL PROCESSOR.
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                                  NPES = NUMBER OF PROCESSORS
      CALL MPI_COMM_SIZE (MPI_COMM_WORLD,NPES,IERR)
C                                  ITASK = MY PROCESSOR NUMBER
      CALL MPI_COMM_RANK (MPI_COMM_WORLD,ITASK,IERR)
C
      DI(1:N) = 0
      DO 5 I=1,NZ
C                                  PCG REQUIRES THE DIAGONAL OF A, FOR 
C                                  "JACOBI" PRECONDITIONING 
         IF (IROW(I).EQ.JCOL(I)) DI(IROW(I)) = A(I)
    5 CONTINUE 
C                                  CALL MPI_ALLREDUCE TO COLLECT DIAGONAL
      CALL MPI_ALLREDUCE(DI,D,N,MPI_DOUBLE_PRECISION,
     & MPI_SUM,MPI_COMM_WORLD,IERR)
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
            API(I) = 0
   20    CONTINUE
         DO 30 IZ=1,NZ
            I = IROW(IZ)
            J = JCOL(IZ)
            API(I) = API(I) + A(IZ)*P(J)
            IF (I.NE.J .AND. SYMM) API(J) = API(J) + A(IZ)*P(I)
   30    CONTINUE
C                                  MPI_ALLREDUCE COLLECTS THE VECTORS API
C                                  (API = LOCAL(A)*P) FROM ALL PROCESSORS
C                                  AND ADDS THEM TOGETHER, THEN SENDS
C                                  THE RESULT, AP, BACK TO ALL PROCESSORS.
         CALL MPI_ALLREDUCE(API,AP,N,MPI_DOUBLE_PRECISION,
     &    MPI_SUM,MPI_COMM_WORLD,IERR)
C                                  PAP = (P,AP)
C                                  RP = (R,D*P)
         PAPI = 0.0
         RPI = 0.0
         DO 40 I=ITASK+1,N,NPES
            PAPI = PAPI + P(I)*AP(I)
            RPI = RPI + R(I)*D(I)*P(I)
   40    CONTINUE
         CALL MPI_ALLREDUCE(PAPI,PAP,1,MPI_DOUBLE_PRECISION,
     &    MPI_SUM,MPI_COMM_WORLD,IERR)
         CALL MPI_ALLREDUCE(RPI,RP,1,MPI_DOUBLE_PRECISION,
     &    MPI_SUM,MPI_COMM_WORLD,IERR)
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=ITASK+1,N,NPES
            X(I) = X(I) + LAMBDA*P(I)
            R(I) = R(I) - LAMBDA*AP(I)/D(I)
   50    CONTINUE
C                                  RAP = (R,AP)
         RAPI = 0.0
         DO 60 I=ITASK+1,N,NPES
            RAPI = RAPI + R(I)*AP(I)
   60    CONTINUE
         CALL MPI_ALLREDUCE(RAPI,RAP,1,MPI_DOUBLE_PRECISION,
     &    MPI_SUM,MPI_COMM_WORLD,IERR)
C                                  ALPHA = -(R,AP)/(P,AP)
         ALPHA = -RAP/PAP
C                                  P = R + ALPHA*P
               IF (SYMM) THEN
         PI(1:N) = 0
         DO 65 I=ITASK+1,N,NPES
            PI(I) = R(I) + ALPHA*P(I)
   65    CONTINUE
         CALL MPI_ALLREDUCE(PI,P,N,MPI_DOUBLE_PRECISION,
     &    MPI_SUM,MPI_COMM_WORLD,IERR)
               ELSE
         DO 70 I=ITASK+1,N,NPES
            P(I) = R(I) + ALPHA*P(I)
   70    CONTINUE
               ENDIF
C                                  RMAX = MAX OF RESIDUAL (R)
         RMAXI = 0
         DO 80 I=ITASK+1,N,NPES
            RMAXI = MAX(RMAXI,ABS(R(I)))
   80    CONTINUE
         CALL MPI_ALLREDUCE(RMAXI,RMAX,1,MPI_DOUBLE_PRECISION,
     &    MPI_MAX,MPI_COMM_WORLD,IERR)
C                                  IF CONVERGED, MERGE PORTIONS OF X
C                                  STORED ON DIFFERENT PROCESSORS
         IF (RMAX.LE.1.D-10*R0MAX) THEN
            IF (ITASK.EQ.0) PRINT *, ' Number of iterations = ',ITER
            CALL MPI_ALLREDUCE(X,R,N,MPI_DOUBLE_PRECISION,
     &       MPI_SUM,MPI_COMM_WORLD,IERR)
            X(1:N) = R(1:N)
            RETURN
         ENDIF
   90 CONTINUE
C                                  PCG DOES NOT CONVERGE
      IF (ITASK.EQ.0) PRINT 100
  100 FORMAT(' ***** PCG does not converge *****')
      RETURN
      END


      SUBROUTINE PBAND(A,N,NLD,NUD,X,B)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C                                 DECLARATIONS FOR ARGUMENTS
      DOUBLE PRECISION A(-NUD-NLD:NLD,*),X(N),B(N)
C                                 DECLARATIONS FOR LOCAL VARIABLES
      DOUBLE PRECISION B_(N),LJI,COLUMNI(N)
      INCLUDE 'mpif.h'
C
C  SUBROUTINE PBAND SOLVES THE LINEAR SYSTEM A*X=B, WHERE A IS A
C    BAND MATRIX, WITH COLUMNS STORED CYCLICALLY ON THE AVAILABLE
C    PROCESSORS.
C
C  ARGUMENTS
C
C             ON INPUT                          ON OUTPUT
C             --------                          ---------
C
C    A      - THE BAND MATRIX, DIMENSIONED      DESTROYED.
C               A(-NUD-NLD:NLD,(N-1)/NPES+1)
C             IN THE CALLING PROGRAM.  ELEMENT
C             (I,J) OF THE MATRIX IS STORED IN
C               A(I-J,(J-1)/NPES+1) 
C             FOR J=ITASK+1+(K-1)*NPES,K=1,2...
C
C    N      - THE SIZE OF MATRIX A.
C
C    NLD    - NUMBER OF NONZERO LOWER DIAGONALS
C             IN A, I.E., NUMBER OF DIAGONALS
C             BELOW THE MAIN DIAGONAL.
C
C    NUD    - NUMBER OF NONZERO UPPER DIAGONALS
C             IN A, I.E., NUMBER OF DIAGONALS
C             ABOVE THE MAIN DIAGONAL.
C
C    X      -                                   AN N-VECTOR CONTAINING
C                                               THE SOLUTION.
C
C    B      - THE RIGHT HAND SIDE N-VECTOR.     
C
C
C-----------------------------------------------------------------------
C                              NPES = NUMBER OF PROCESSORS
      CALL MPI_COMM_SIZE (MPI_COMM_WORLD,NPES,IERR)
C                              ITASK = MY PROCESSOR NUMBER (0,1,...,NPES-1).
C                              I WILL NEVER TOUCH ANY COLUMNS OF A EXCEPT
C                              MY COLUMNS, ITASK+1+ K*NPES, K=0,1,2,...
      CALL MPI_COMM_RANK (MPI_COMM_WORLD,ITASK,IERR)
C                              COPY B TO B_, SO B WILL NOT BE ALTERED
      B_(1:N) = B(1:N) 
C                              BEGIN FORWARD ELIMINATION
      DO 35 I=1,N
C                              JTASK IS PROCESSOR THAT OWNS ACTIVE COLUMN
         JTASK = MOD(I-1,NPES)
         IF (ITASK.EQ.JTASK) THEN
C                              IF JTASK IS ME, SAVE ACTIVE COLUMN IN
C                              VECTOR COLUMNI
            DO 5 J=I,MIN(I+NLD,N)
               IME = (I-1)/NPES+1
               COLUMNI(J) = A(J-I,IME)
    5       CONTINUE
         ENDIF
C                              RECEIVE COLUMNI FROM PROCESSOR JTASK
         CALL MPI_BCAST(COLUMNI(I),MIN(I+NLD,N)-I+1,MPI_DOUBLE_PRECISION
     &    ,JTASK,MPI_COMM_WORLD,IERR)
C                              SEARCH FROM A(I,I) ON DOWN FOR LARGEST
C                              POTENTIAL PIVOT, A(L,I)
         BIG = ABS(COLUMNI(I))
         L = I
         DO 10 J=I+1,MIN(I+NLD,N)
            IF (ABS(COLUMNI(J)).GT.BIG) THEN
               BIG = ABS(COLUMNI(J))
               L = J
            ENDIF
   10    CONTINUE
C                              IF LARGEST POTENTIAL PIVOT IS ZERO,
C                              MATRIX IS SINGULAR
         IF (BIG.EQ.0.0) GO TO 50
C                              I0 IS FIRST COLUMN >= I THAT BELONGS TO ME
         L0 = (I-1+NPES-(ITASK+1))/NPES
         I0 = ITASK+1+L0*NPES
C                              SWITCH ROW I WITH ROW L, TO BRING UP
C                              LARGEST PIVOT; BUT ONLY IN MY COLUMNS
         DO 15 K=I0,MIN(I+NUD+NLD,N),NPES
            KME = (K-1)/NPES+1
            TEMP = A(L-K,KME)
            A(L-K,KME) = A(I-K,KME)
            A(I-K,KME) = TEMP
   15    CONTINUE
         TEMP = COLUMNI(L)
         COLUMNI(L) = COLUMNI(I)
         COLUMNI(I) = TEMP
C                              SWITCH B_(I) AND B_(L)
         TEMP = B_(L)
         B_(L) = B_(I)
         B_(I) = TEMP
         DO 25 K=I0,MIN(I+NLD+NUD,N),NPES
            KME = (K-1)/NPES+1
            LJI = A(I-K,KME)/COLUMNI(I)
            IF (LJI.NE.0.0) THEN
C                              SUBTRACT MULTIPLE OF ROW I FROM ROW J;
C                              BUT ONLY IN MY COLUMNS
               DO 20 J=I+1,MIN(I+NLD,N)
                  A(J-K,KME) = A(J-K,KME) - LJI*COLUMNI(J)
   20          CONTINUE
            ENDIF
   25    CONTINUE         
         LJI = B_(I)/COLUMNI(I)
         DO 30 J=I+1,MIN(I+NLD,N)
C                              SUBTRACT MULTIPLE OF B_(I) FROM B_(J)
             B_(J) = B_(J) - LJI*COLUMNI(J)
   30    CONTINUE
   35 CONTINUE
C                              SOLVE U*X=B_ USING BACK SUBSTITUTION.
      DO 45 I=N,1,-1
C                              I0 IS FIRST COLUMN >= I+1 THAT BELONGS
C                              TO ME
         L0 = (I+NPES-(ITASK+1))/NPES
         I0 = ITASK+1+L0*NPES
         SUMI = 0.0
         DO 40 J=I0,MIN(I+NLD+NUD,N),NPES
            JME = (J-1)/NPES+1
            SUMI = SUMI + A(I-J,JME)*X(J)
   40    CONTINUE 
         CALL MPI_ALLREDUCE(SUMI,SUM,1,MPI_DOUBLE_PRECISION,
     &    MPI_SUM,MPI_COMM_WORLD,IERR)
C                              JTASK IS PROCESSOR THAT OWNS A(I,I)
         JTASK = MOD(I-1,NPES)
C                              IF JTASK IS ME, CALCULATE X(I)
         IF (ITASK.EQ.JTASK) THEN 
            IME = (I-1)/NPES+1
            X(I) = (B_(I)-SUM)/A(0,IME)
         ENDIF
C                              RECEIVE X(I) FROM PROCESSOR JTASK
         CALL MPI_BCAST(X(I),1,MPI_DOUBLE_PRECISION,JTASK,
     &   MPI_COMM_WORLD,IERR)
   45 CONTINUE
      RETURN
   50 IF (ITASK.EQ.0) PRINT 55
   55 FORMAT (' ***** The matrix is singular *****')
      RETURN
      END
