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