SUBROUTINE DDTSL3(N,C,D,E,LDB,NSYS,B) C C*********************************************************************** C C THIS PROGRAM SOLVES THE LINEAR SYSTEMS C C A*X(J) = B(J), J = 1,...,NSYS, C C WHERE THE MATRIX A IS STRICTLY OR IRREDUCIBLY DIAGONALLY C DOMINANT. THE METHOD IMPLEMENTED IS A VARIANT OF THE "MILLAY" C ALGORITHM FOR SOLVING POSITIVE DEFINITE TRIDIAGONAL SYSTEMS. C SEE J.J. DONGARRA, J.R. BUNCH, C.B. MOLER AND G.W. STEWART (1979), C LINPACK USERS' GUIDE, SIAM PUBLICATIONS, PHILADELPHIA, PA., AND C AND J.J. DONGARRA AND G.W. STEWART (1984), LINPACK - A PACKAGE C FOR SOLVING LINEAR SYSTEMS, SOURCES AND DEVELOPMENT OF MATHE- C EMATICAL SOFTWARE, (W.R. COWELL, EDITOR), PRENTICE-HALL INC., NEW C JERSEY, PP. 20-48. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C ***** PARAMETERS ***** C C *** ON ENTRY ... C C N - INTEGER C THE ORDER OF THE TRIDIAGONAL MATRIX A. C C C - DOUBLE PRECISION(N) C THE SUBDIAGONAL OF THE TRIDIAGONAL MATRIX A. C C(2) THROUGH C(N) SHOULD CONTAIN THE SUBDIAGONAL C ELEMENTS OF A. C C D - DOUBLE PRECISION(N) C THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL C MATRIX A. C C E - DOUBLE PRECISION(N) C THE SUPERDIAGONAL OF THE TRIDIAGONAL MATRIX A. C E(1) THROUGH E(N-1) SHOULD CONTAIN THE SUPER- C DIAGONAL ELEMENTS OF A. C C LDB - INTEGER C THE LEADING DIMENSION OF B IN THE CALLING C PROGRAM. C LDB MUST BE GREATER THAN OR EQUAL TO N. C C NSYS - INTEGER C THE NUMBER OF RIGHT HAND SIDE VECTORS. C IF NSYS = 1, LDB MAY BE SET EQUAL TO N. C C B - DOUBLE PRECISION(LDB,NSYS) C THE SET OF RIGHT HAND SIDE VECTORS. C THE J-TH COLUMN OF B SHOULD CONTAIN THE C RIGHT HAND SIDE OF THE J-TH SYSTEM OF C EQUATIONS, J = 1,...,NSYS. C IN THE CALLING PROGRAM, B IS DIMENSIONED AS C B(LDB,NCOLS), FOR SOME NCOLS GREATER THAN OR C EQUAL TO NSYS. C C *** ON RETURN ... C C C,D,E - VECTORS CONTAINING THE DESIRED DECOMPOSITION OF C THE MATRIX A. C C B - DOUBLE PRECISION(LDB,NSYS) C THE SET OF SOLUTION VECTORS. C THE J-TH COLUMN OF B CONTAINS THE SOLUTION OF C OF THE J-TH SYSTEM OF EQUATIONS, J = 1,...,NSYS. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C ***** AUXILIARY PROGRAMS ***** C C DDTFAC(N,C,D,E) C - DECOMPOSES THE MATRIX A USING A VARIANT C OF THE MILLAY ALGORITHM FOR DIAGONALLY C DOMINANT TRIDIAGONAL MATRICES AND IS USED FOR C THIS PURPOSE IN D D T S L 3. C THE ARGUMENTS ARE ALL AS IN D D T S L 3. C C DDTSYS(N,C,D,E,LDB,NSYS,B) C - SOLVES THE SYSTEMS C C A*X(J) = B(J), J = 1,...,NSYS, C C ONCE THE MATRIX A IS DECOMPOSED. C THE ARGUMENTS ARE ALL AS IN D D T S L 3. C C*********************************************************************** C DOUBLE PRECISION C,D,E,B INTEGER N,NSYS,LDB DIMENSION C(N),D(N),E(N),B(LDB,NSYS) CALL DDTFAC(N,C,D,E) CALL DDTSYS(N,C,D,E,LDB,NSYS,B) RETURN END SUBROUTINE DDTFAC(N,C,D,E) C C*********************************************************************** C C D D T F A C DECOMPOSES THE DIAGONALLY DOMINANT TRIDIAGONAL MATRIX C A USING A VARIANT OF THE "MILLAY" ALGORITHM. C THE MATRIX A IS STORED IN THE VECTORS C, D, E. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C ***** PARAMETERS ***** C C *** ON ENTRY ... C C N - INTEGER C THE ORDER OF THE TRIDIAGONAL MATRIX A. C C C - DOUBLE PRECISION(N) C THE SUBDIAGONAL OF THE TRIDIAGONAL MATRIX A. C C(2) THROUGH C(N) SHOULD CONTAIN THE SUBDIAGONAL C ELEMENTS OF A. C C D - DOUBLE PRECISION(N) C THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL C MATRIX A. C C E - DOUBLE PRECISION(N) C THE SUPERDIAGONAL OF THE TRIDIAGONAL MATRIX A. C E(1) THROUGH E(N-1) SHOULD CONTAIN THE SUPER- C DIAGONAL ELEMENTS OF A. C C *** ON RETURN ... C C C,D,E - VECTORS CONTAINING THE DESIRED DECOMPOSITION OF C THE MATRIX A. C C*********************************************************************** C DOUBLE PRECISION C,D,E,ONE INTEGER N,NPLS1,M,MPLS1,MPLS2,I,I2,NMINM DIMENSION C(N),D(N),E(N) DATA ONE/1.0D0/ C C*********************************************************************** C C DEFINE CONSTANTS USED THROUGHOUT .... C C*********************************************************************** C NPLS1 = N+1 M = (N-1)/2 NMINM = N - M MPLS1 = M+1 MPLS2 = M+2 C C*********************************************************************** C C CHECK FOR 1 X 1 CASE. IF N = 1, RETURN .... C C*********************************************************************** C D(1) = ONE/D(1) IF(N .EQ. 1) RETURN C C*********************************************************************** C C IF N IS GREATER THAN 2, SIMULTANEOUSLY ... C ZERO TOP HALF OF SUPERDIAGONAL AND THE BOTTOM HALF OF C SUBDIAGONAL USING COLUMN ELIMINATIONS .... C C*********************************************************************** C IF(N .GT. 2) THEN D(N) = ONE/D(N) DO 10 I = 1,M-1 E(I) = E(I)*D(I) D(I+1) = ONE/(D(I+1)-C(I+1)*E(I)) I2 = NPLS1-I C(I2) = C(I2)*D(I2) D(I2-1) = ONE/(D(I2-1)-E(I2-1)*C(I2)) 10 CONTINUE E(M) = E(M)*D(M) D(MPLS1) = D(MPLS1) - C(MPLS1)*E(M) C(NPLS1-M) = C(NPLS1-M)*D(NPLS1-M) D(NMINM) = D(NMINM) - E(NMINM)*C(NPLS1-M) D(MPLS1) = ONE/D(MPLS1) ENDIF IF (MOD(N,2) .EQ. 0) THEN C ************************************************************************ C C CLEAN UP POSSIBLE 2 X 2 BLOCK AT CENTER .... C ************************************************************************ C E(MPLS1) = E(MPLS1)*D(MPLS1) D(MPLS2) = ONE/(D(MPLS2)-C(MPLS2)*E(MPLS1)) ENDIF C RETURN END SUBROUTINE DDTSYS(N,C,D,E,LDB,NSYS,B) C C*********************************************************************** C C D D T S Y S SOLVES THE LINEAR SYSTEMS WITH COEFFICIENT MATRIX A C AND RIGHT HAND SIDES B(J), J = 1,...,NSYS, USING THE DECOMPOSITION C OF A GENERATED IN D D T F A C. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C ***** PARAMETERS ***** C C *** ON ENTRY ... C C N - INTEGER C THE ORDER OF THE TRIDIAGONAL MATRIX A. C C C,D,E - DOUBLE PRECISION C OUTPUT FROM DDTFAC. C C LDB - INTEGER C THE LEADING DIMENSION OF B IN THE CALLING C PROGRAM. C LDB MUST BE GREATER THAN OR EQUAL TO N. C C NSYS - INTEGER C THE NUMBER OF RIGHT HAND SIDE VECTORS. C IF NSYS = 1, LDB MAY BE SET EQUAL TO N. C C B - DOUBLE PRECISION(LDB,NSYS) C THE SET OF RIGHT HAND SIDE VECTORS. C THE J-TH COLUMN OF B SHOULD CONTAIN THE C RIGHT HAND SIDE OF THE J-TH SYSTEM OF C EQUATIONS, J = 1,...,NSYS. C C *** ON RETURN ... C C B - THE COLUMNS OF B ARE THE SOLUTION VECTORS. C C*********************************************************************** C DOUBLE PRECISION C,D,E,B INTEGER N,LDB,NSYS INTEGER I,M,MPLS1,MPLS2,I1,I2,NPLS1,NMINM,ICOL DIMENSION C(N),D(N),E(N),B(LDB,NSYS) C C*********************************************************************** C C DEFINE THE CONSTANTS USED THROUGHOUT .... C C*********************************************************************** C NPLS1 = N + 1 M = (N-1)/2 NMINM = N - M MPLS1 = M + 1 MPLS2 = M + 2 C C********************************************************************** C C FOR EACH SYSTEM OF EQUATIONS ..... C C*********************************************************************** C DO 30 ICOL = 1,NSYS C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C CHECK FOR 1 BY 1 CASE. IF N = 1, EXIT .... C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C B(1,ICOL) = B(1,ICOL)*D(1) IF(N.EQ.1)GO TO 30 C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C IF N IS GREATER THAN OR EQUAL TO 2, PERFORM FORWARD C ELIMINATION STARTING FROM THE TOP AND BOTTOM SIMUL- C TANEOUSLY, AND WORKING INWARD TOWARDS THE CENTER .... C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IF(N.GT.2) THEN B(N,ICOL) = B(N,ICOL)*D(N) DO 10 I = 1,M-1 B(I+1,ICOL) = (B(I+1,ICOL)-C(I+1)*B(I,ICOL))*D(I+1) I2 = NPLS1-I B(I2-1,ICOL) = (B(I2-1,ICOL)-E(I2-1)*B(I2,ICOL))*D(I2-1) 10 CONTINUE B(MPLS1,ICOL) = B(MPLS1,ICOL) - C(MPLS1)*B(M,ICOL) B(NMINM,ICOL) = B(NMINM,ICOL) - E(NMINM)*B(NPLS1-M,ICOL) B(MPLS1,ICOL) = B(MPLS1,ICOL)*D(MPLS1) ENDIF IF(MOD(N,2) .EQ. 0)THEN C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C CLEAN UP POSSIBLE 2 X 2 BLOCK AT CENTER, C AND IF N = 2, EXIT ..... C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C B(MPLS2,ICOL) = (B(MPLS2,ICOL)-C(MPLS2)*B(MPLS1,ICOL))* * D(MPLS2) B(MPLS1,ICOL) = B(MPLS1,ICOL)-E(MPLS1)*B(MPLS2,ICOL) ENDIF IF(N .EQ. 2) GO TO 30 C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C BACK SOLVE STARTING AT THE CENTER, MOVING TOWARDS THE TOP C AND BOTTOM SIMULTANEOUSLY..... C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DO 20 I = 1,M I1 = MPLS1-I B(I1,ICOL) = B(I1,ICOL)-E(I1)*B(I1+1,ICOL) I2 = NMINM+I B(I2,ICOL) = B(I2,ICOL)-C(I2)*B(I2-1,ICOL) 20 CONTINUE 30 CONTINUE RETURN END