PROGRAM BACSUBS C NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994 C To accompany the text: C NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992 C Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A. C This free software is complements of the author. C C Algorithm 3.1 (Back Substitution). C Section 3.3, Upper-Triangular Linear Systems, Page 145 C PARAMETER(MaxR=10) INTEGER InRC,N REAL A,B,DET,X CHARACTER Ans*1,Ach*1,Bch*1,Xch*1 DIMENSION A(1:MaxR,1:MaxR),B(1:MaxR),X(1:MaxR) 10 CALL INPUTS(A,B,N,InRC,DET) CALL BACKSUB(A,B,N,X,DET) CALL RESULTS(A,B,N,X,DET) WRITE(9,*)' ' WRITE(9,*)'WANT TO SOLVE ANOTHER UPPER-TRIANGULAR SYSTEM ? ' READ(9,'(A)') Ans IF (Ans.EQ.'Y' .OR. Ans.EQ.'y') GOTO 10 STOP END SUBROUTINE BACKSUB(A,B,N,X,DET) PARAMETER(MaxR=10) INTEGER J,K,N,R REAL A,B,DET,SUM,X DIMENSION A(1:MaxR,1:MaxR),B(1:MaxR),X(1:MaxR) DO K=1,N IF (A(K,K).EQ.0) THEN DET=0 RETURN ENDIF ENDDO DET=A(N,N) X(N)=B(N)/A(N,N) DO R=(N-1),1,-1 DET=DET*A(R,R) SUM= 0 DO J=(R+1),N SUM=SUM+A(R,J)*X(J) ENDDO X(R)= (B(R)-SUM)/A(R,R) ENDDO RETURN END SUBROUTINE XBACKSUB(A,B,N,X,DET) C This subroutine uses simulated WHILE loop(s). PARAMETER(MaxR=10) INTEGER J,K,N,R REAL A,B,DET,SUM,X DIMENSION A(1:MaxR,1:MaxR),B(1:MaxR),X(1:MaxR) DO 10 K=1,N IF (A(K,K).EQ.0) THEN DET=0 RETURN ENDIF 10 CONTINUE DET=A(N,N) X(N)=B(N)/A(N,N) DO 30 R=(N-1),1,-1 DET=DET*A(R,R) SUM= 0 DO 20 J=(R+1),N SUM=SUM+A(R,J)*X(J) 20 CONTINUE X(R)= (B(R)-SUM)/A(R,R) 30 CONTINUE RETURN END SUBROUTINE INPUTS(A,B,N,InRC,DET) PARAMETER(MaxR=10) INTEGER C,I,InRC,N,R REAL A,B,DET CHARACTER Ans*1,Ach*1 DIMENSION A(1:MaxR,1:MaxR),B(1:MaxR) DO 10 I=1,15 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)'SOLUTION OF AN UPPER-TRIANGULAR LINEAR SYSTEM:' WRITE(9,*)' ' WRITE(9,*)' A*X = B' WRITE(9,*)' ' WRITE(9,*)'A IS AN N BY N UPPER-TRIANGULAR MATRIX.' WRITE(9,*)' ' WRITE(9,*)'B IS AN N DIMENSIONAL VECTOR OF CONSTANTS.' WRITE(9,*)' ' WRITE(9,*)'X IS THE N DIMENSIONAL SOLUTION VECTOR.' WRITE(9,*)' ' WRITE(9,*)'ENTER THE NUMBER OF EQUATIONS N = ' READ(9,*) N InRC=0 WRITE(9,*)'WANT TO INPUT COLUMNS ? ' READ(9,'(A)') Ans IF (Ans.EQ.'Y' .OR. Ans.EQ.'y') THEN InRC=1 ENDIF WRITE(9,*)' ' WRITE(9,*)'ENTER THE UPPER-TRIANGULAR MATRIX A(I,J):' CALL MATIN('A',A,N,InRC) WRITE(9,*)' ' WRITE(9,*)'ENTER THE COLUMN VECTOR B(J):' CALL VECIN('B',B,N) RETURN END SUBROUTINE MATIN(Ach,A,N,InRC) PARAMETER(MaxR=10) INTEGER Col,InRC,N,Row REAL A CHARACTER Ach*1 DIMENSION A(1:MaxR,1:MaXR) DO 20 Row=1,N DO 10 Col=1,N A(Row,Col)=0 10 CONTINUE 20 CONTINUE WRITE(9,*)'GIVE THE ELEMENTS OF THE MATRIX ',Ach WRITE(9,*)' ' IF (InRC.EQ.0) THEN DO 40 Row=1,N WRITE(9,*)'GIVE THE ELEMENTS OF ROW ',Row WRITE(9,*)' ' DO 30 Col=Row,N WRITE(9,*) Ach,'(',Row,',',Col,') = ' READ(9,*) A(Row,Col) 30 CONTINUE 40 CONTINUE ELSE DO 60 Col=1,N WRITE(9,*)'GIVE THE ELEMENTS OF COLUMN ',Col WRITE(9,*)' ' DO 50 Row=1,Col WRITE(9,*) Ach,'(',Row,',',Col,') = ' READ(9,*) A(Row,Col) 50 CONTINUE 60 CONTINUE ENDIF RETURN END SUBROUTINE VECIN(Bch,B,N) PARAMETER(MaxR=10) INTEGER N,Row REAL B CHARACTER Bch*1 DIMENSION B(1:MaxR) WRITE(9,*)'GIVE THE ELEMENTS OF THE COLUMN VECTOR ',Bch WRITE(9,*)' ' DO 10 Row=1,N WRITE(9,*) Bch,'(',Row,') = ' READ(9,*) B(Row) 10 CONTINUE RETURN END SUBROUTINE MATPRINT(A,N) PARAMETER(MaxR=10) INTEGER Col,N,Row REAL A DIMENSION A(1:MaxR,1:MaxR) DO 10 Row=1,N WRITE(9,999) ( A(Row,Col), Col=1,N ) 10 CONTINUE 999 FORMAT(1X,5F15.7,4X) RETURN END SUBROUTINE VECPRINT(Bch,B,Xch,X,N) PARAMETER(MaxR=10) INTEGER N,Row REAL B,X CHARACTER Bch*1,Xch*1 DIMENSION B(1:MaxR),X(1:MaxR) WRITE(9,*) Bch,' COEFFICIENT VECTOR ',Xch, +' SOLUTION VECTOR' WRITE(9,*) DO 10 Row=1,N WRITE(9,*) Bch,'(',Row,') =',B(Row),' ', +Xch,'(',Row,') =',X(Row) 10 CONTINUE RETURN END SUBROUTINE RESULTS(A,B,N,X,DET) PARAMETER(MaxR=10) INTEGER Col,I,N,Row REAL A,B,DET,X CHARACTER Bch*1,Xch*1 DIMENSION A(1:MaxR,1:MaxR),B(1:MaxR),X(1:MaxR) DO 10 I=1,15 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)'THE UPPER-TRIANGULAR MATRIX A(I,J) IS:' WRITE(9,*)' ' CALL MATPRINT(A,N) IF (DET.EQ.0) THEN WRITE(9,*)'THE MATRIX IS SINGULAR.' WRITE(9,*) WRITE(9,*)'THE BACK-SUBSTITUTION ALGORITHM DOES NOT APPLY.' ELSE WRITE(9,*) CALL VECPRINT('B',B,'X',X,N) ENDIF WRITE(9,*)' ' WRITE(9,*)'THE DETERMINANT`S VALUE IS DET A = ',DET RETURN END PROGRAM GAUSJORD C NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994 C To accompany the text: C NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992 C Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A. C This free software is complements of the author. C C Algorithm 3.2 (Upper-Triangularization Followed by Back Substitution). C Section 3.4, Gaussian Elimination and Pivoting, Page 156 C PARAMETER(MaxR=10,MaxC=11) INTEGER InRC,N REAL A,A1,DET,X CHARACTER Ans*1,Ach*1,Bch*1,Xch*1 DIMENSION A(1:MaxR,1:MaxC),A1(1:MaxR,1:MaxC),X(1:MaxR) 10 CALL INPUTS(A,A1,N,InRC) CALL GAUSSJO(A,X,N,DET) CALL RESULTS(A1,X,N,DET) WRITE(9,*)' ' WRITE(9,*)'WANT TO SOLVE ANOTHER LINEAR SYSTEM ? ' READ(9,'(A)') Ans IF (Ans.EQ.'Y' .OR. Ans.EQ.'y') GOTO 10 STOP END SUBROUTINE GAUSSJO(A,X,N,DET) PARAMETER(MaxR=10,MaxC=11) INTEGER C,J,K,N,P,R,Row,T REAL A,DET,M,SUM,X DIMENSION A(1:MaxR,1:MaxC),Row(1:MaxR),X(1:MaxR) DO J=1,N Row(J)=J ENDDO DO P=1,(N-1) DO K=(P+1),N IF (ABS(A(Row(K),P)).GT.ABS(A(Row(P),P))) THEN T = Row(P) Row(P)=Row(K) Row(K)=T ENDIF ENDDO IF (A(Row(P),P).EQ.0) THEN DET=0 RETURN ENDIF DO K=(P+1),N M = A(Row(K),P)/A(Row(P),P) DO C=(P+1),(N+1) A(Row(K),C)=A(Row(K),C)-M*A(Row(P),C) ENDDO ENDDO ENDDO DET=A(Row(N),N) IF (DET.NE.0) THEN X(N)=A(Row(N),(N+1))/A(Row(N),N) DO K=(N-1),1,-1 DET=DET*A(Row(K),K) SUM=0 DO C=(K+1),N SUM=SUM+A(Row(K),C)*X(C) ENDDO X(K)=(A(Row(K),(N+1))-SUM)/A(Row(K),K) ENDDO ENDIF RETURN END SUBROUTINE XGAUSSJO(A,X,N,DET) C This subroutine uses simulated WHILE loop(s). PARAMETER(MaxR=10,MaxC=11) INTEGER C,J,K,N,P,R,Row,T REAL A,DET,M,SUM,X DIMENSION A(1:MaxR,1:MaxC),Row(1:MaxR),X(1:MaxR) DO 10 J=1,N Row(J)=J 10 CONTINUE DO 50 P=1,(N-1) DO 20 K=(P+1),N IF (ABS(A(Row(K),P)).GT.ABS(A(Row(P),P))) THEN T = Row(P) Row(P)=Row(K) Row(K)=T ENDIF 20 CONTINUE IF (A(Row(P),P).EQ.0) THEN DET=0 RETURN ENDIF DO 40 K=(P+1),N M = A(Row(K),P)/A(Row(P),P) DO 30 C=(P+1),(N+1) A(Row(K),C)=A(Row(K),C)-M*A(Row(P),C) 30 CONTINUE 40 CONTINUE 50 CONTINUE DET=A(Row(N),N) IF (DET.NE.0) THEN X(N)=A(Row(N),(N+1))/A(Row(N),N) DO 70 K=(N-1),1,-1 DET=DET*A(Row(K),K) SUM=0 DO 60 C=(K+1),N SUM=SUM+A(Row(K),C)*X(C) 60 CONTINUE X(K)=(A(Row(K),(N+1))-SUM)/A(Row(K),K) 70 CONTINUE ENDIF RETURN END SUBROUTINE INPUTS(A,A1,N,InRC) PARAMETER(MaxR=10,MaxC=11) INTEGER Col,I,InRC,N,Row REAL A CHARACTER Ans*1,Ach*1 DIMENSION A(1:MaxR,1:MaxC),A1(1:MaxR,1:MaxC) DO 10 I=1,15 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)'SOLUTION OF A LINEAR SYSTEM:' WRITE(9,*)' ' WRITE(9,*)' A*X = B' WRITE(9,*)' ' WRITE(9,*)'GAUSSIAN ELIMINATION IS USED TO UPPER-TRIANGULARIZE,' WRITE(9,*)' ' WRITE(9,*)'AND THEN BACK-SUBSTITUTION IS USED TO SOLVE FOR X.' WRITE(9,*)' ' WRITE(9,*)'A IS AN N BY N NONSINGULAR MATRIX.' WRITE(9,*)' ' WRITE(9,*)'B IS AN N DIMENSIONAL VECTOR OF CONSTANTS.' WRITE(9,*)' ' WRITE(9,*)'X IS THE N DIMENSIONAL SOLUTION VECTOR.' WRITE(9,*)' ' WRITE(9,*)'ENTER THE NUMBER OF EQUATIONS N = ' READ(9,*) N InRC=0 WRITE(9,*)'WANT TO INPUT COLUMNS ? ' READ(9,'(A)') Ans IF (Ans.EQ.'Y' .OR. Ans.EQ.'y') THEN InRC=1 ENDIF WRITE(9,*)' ' WRITE(9,*)'ENTER THE COEFFICIENT MATRIX A:' CALL MATIN('A',A,N,InRC) WRITE(9,*)' ' WRITE(9,*)'ENTER THE COLUMN VECTOR B:' CALL VECIN('B',A,N,InRC) DO 30 Row=1,N DO 20 Col=1,N+1 A1(Row,Col)=A(Row,Col) 20 CONTINUE 30 CONTINUE RETURN END SUBROUTINE MATIN(Ach,A,N,InRC) PARAMETER(MaxR=10,MaxC=11) INTEGER Col,InRC,N,Row REAL A CHARACTER Ach*1 DIMENSION A(1:MaxR,1:MaxC) DO 20 Row=1,N DO 10 Col=1,(N+1) A(Row,Col)=0 10 CONTINUE 20 CONTINUE WRITE(9,*)'GIVE THE ELEMENTS OF THE MATRIX ',Ach WRITE(9,*)' ' IF (InRC.EQ.0) THEN DO 40 Row=1,N WRITE(9,*)'GIVE THE ELEMENTS OF ROW ',Row WRITE(9,*)' ' DO 30 Col=1,N WRITE(9,*) Ach,'(',Row,',',Col,') = ' READ(9,*) A(Row,Col) 30 CONTINUE 40 CONTINUE ELSE DO 60 Col=1,N WRITE(9,*)'GIVE THE ELEMENTS OF COLUMN ',Col WRITE(9,*)' ' DO 50 Row=1,N WRITE(9,*) Ach,'(',Row,',',Col,') = ' READ(9,*) A(Row,Col) 50 CONTINUE 60 CONTINUE ENDIF RETURN END SUBROUTINE VECIN(Ach,A,N,InRC) PARAMETER(MaxR=10,MaxC=11) INTEGER Col,InRC,N,Row REAL A CHARACTER Ach*1 DIMENSION A(1:MaxR,1:MaxC) WRITE(9,*)'GIVE THE ELEMENTS OF THE COLUMN VECTOR ',Ach WRITE(9,*)' ' Col=N+1 DO 10 Row=1,N WRITE(9,*) Ach,'(',Row,') = ' READ(9,*) A(Row,Col) 10 CONTINUE RETURN END SUBROUTINE MATPRINT(A,N) PARAMETER(MaxR=10,MaxC=11) INTEGER Col,N,Row REAL A DIMENSION A(1:MaxR,1:MaxC) DO 10 Row=1,N WRITE(9,999) ( A(Row,Col), Col=1,N ) 10 CONTINUE 999 FORMAT(1X,5F15.7,4X) RETURN END SUBROUTINE VECPRINT(Bch,Xch,A,X,N) PARAMETER(MaxR=10,MaxC=11) INTEGER Col,N,Row REAL A,X CHARACTER Bch*1,Xch*1 DIMENSION A(1:MaxR,MaxC),X(1:MaxR) Col=N+1 WRITE(9,*) Bch,' COEFFICIENT VECTOR ',Xch, +' SOLUTION VECTOR' WRITE(9,*) DO 10 Row=1,N WRITE(9,*) Bch,'(',Row,') =',A(Row,Col),' ', +Xch,'(',Row,') =',X(Row) 10 CONTINUE RETURN END SUBROUTINE RESULTS(A,X,N,DET) PARAMETER(MaxR=10,MaxC=11) INTEGER I,N REAL A,DET,X CHARACTER Bch*1,Xch*1 DIMENSION A(1:MaxR,1:MaxC),X(1:MaxR) DO 10 I=1,15 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)'COMPUTATION OF THE SOLUTION FOR THE LINEAR SYSTEM A*X = + B.' WRITE(9,*)' ' WRITE(9,*)'GAUSSIAN ELIMINATION WAS USED TO UPPER-TRIANGULARIZE,' WRITE(9,*)' ' WRITE(9,*)'AND THEN BACK-SUBSTITUTION IS USED TO SOLVE FOR X.' WRITE(9,*)' ' WRITE(9,*)'THE COEFFICIENT MATRIX A IS:' WRITE(9,*)' ' CALL MATPRINT(A,N) IF (DET.EQ.0) THEN WRITE(9,*) WRITE(9,*)'THE MATRIX IS SINGULAR.' WRITE(9,*) WRITE(9,*)'A ZERO PIVOT ELEMENT WAS ENCOUNTERED.' WRITE(9,*)' ' WRITE(9,*)'THE METHOD DOES NOT APPLY.' ELSE WRITE(9,*) CALL VECPRINT('B','X',A,X,N) ENDIF WRITE(9,*)' ' WRITE(9,*)'THE DETERMINANT`S VALUE IS DET A = ',DET RETURN END PROGRAM LUFACTOR C NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994 C To accompany the text: C NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992 C Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A. C This free software is complements of the author. C C Algorithm 3.3 (PA = LU Factorization with Pivoting). C Section 3.6, Triangular Factorizaton, Page 175 C PARAMETER(MaxR=10) INTEGER InRC,N,Row REAL A,A1,B,DET,X CHARACTER Ans*1,Ach*1,Bch*1,Xch*1 DIMENSION A(1:MaxR,1:MaxR),A1(1:MaxR,1:MaxR) DIMENSION B(1:MaxR),Row(1:MaxR),X(1:MaxR) 10 CALL INPUTS(A,A1,B,N,InRC) CALL FACTOR(A,Row,N,DET) 20 WRITE(9,*)' ' WRITE(9,*)'ENTER THE COLUMN VECTOR B:' CALL VECIN('B',B,N) IF (DET.NE.0) THEN CALL SOLVE(A,Row,B,X,N) ENDIF CALL RESULTS(A1,B,X,N,DET) IF (DET.NE.0) THEN WRITE(9,*)'WANT TO SOLVE A*X=B WITH A NEW VECTOR B ? ' READ(9,'(A)') Ans IF (Ans.EQ.'Y' .OR. Ans.EQ.'y') GOTO 20 ENDIF WRITE(9,*)'WANT TO SOLVE ANOTHER LINEAR SYSTEM ? ' READ(9,'(A)') Ans IF (Ans.EQ.'Y' .OR. Ans.EQ.'y') GOTO 10 END SUBROUTINE FACTOR(A,Row,N,DET) PARAMETER(MaxR=10) INTEGER C,J,InRC,K,N,P,Row,RowK,RowP,T REAL A,A1,B,DET,X CHARACTER Ans*1,Ach*1,Bch*1,Xch*1 DIMENSION A(1:MaxR,1:MaxR),Row(1:MaxR) DET=1 DO J=1,N Row(J)=J ENDDO DO P=1,(N-1) DO K=(P+1),N IF (ABS(A(Row(K),P)).GT.ABS(A(Row(P),P))) THEN T=Row(P) Row(P)=Row(K) Row(K)=T DET=-DET ENDIF ENDDO DET=DET*A(Row(P),P) IF (DET.EQ.0) RETURN DO K=(P+1),N RowK=Row(K) RowP=Row(P) A(RowK,P)=A(RowK,P)/A(RowP,P) DO C=(P+1),N A(RowK,C)=A(RowK,C)-A(RowK,P)*A(RowP,C) ENDDO ENDDO ENDDO DET=DET*A(Row(N),N) RETURN END SUBROUTINE SOLVE(A,Row,B,X,N) PARAMETER(MaxR=10) INTEGER C,K,N,Row,RowK REAL A,B,Sum,X DIMENSION A(1:MaxR,1:MaxR),B(1:MaxR),Row(1:MaxR),X(1:MaxR) DO K=1,N IF (A(Row(K),K).EQ.0) RETURN ENDDO X(1)=B(Row(1)) DO K=2,N Sum=0 RowK=Row(K) DO C=1,(K-1) Sum=Sum+A(RowK,C)*X(C) ENDDO X(K)=B(RowK)-Sum ENDDO X(N)=X(N)/A(Row(N),N) DO K=(N-1),1,-1 Sum=0 RowK=Row(K) DO C=(K+1),N Sum=Sum+A(RowK,C)*X(C) ENDDO X(K)=(X(K)-Sum)/A(RowK,K) ENDDO RETURN END SUBROUTINE XFACTOR(A,Row,N,DET) C This subroutine uses simulated WHILE loop(s). PARAMETER(MaxR=10) INTEGER C,J,InRC,K,N,P,Row,RowK,RowP,T REAL A,A1,B,DET,X CHARACTER Ans*1,Ach*1,Bch*1,Xch*1 DIMENSION A(1:MaxR,1:MaxR),Row(1:MaxR) DET=1 DO 10 J=1,N Row(J)=J 10 CONTINUE DO 50 P=1,(N-1) DO 20 K=(P+1),N IF (ABS(A(Row(K),P)).GT.ABS(A(Row(P),P))) THEN T=Row(P) Row(P)=Row(K) Row(K)=T DET=-DET ENDIF 20 CONTINUE DET=DET*A(Row(P),P) IF (DET.EQ.0) RETURN DO 40 K=(P+1),N RowK=Row(K) RowP=Row(P) A(RowK,P)=A(RowK,P)/A(RowP,P) DO 30 C=(P+1),N A(RowK,C)=A(RowK,C)-A(RowK,P)*A(RowP,C) 30 CONTINUE 40 CONTINUE 50 CONTINUE DET=DET*A(Row(N),N) RETURN END SUBROUTINE XSOLVE(A,Row,B,X,N) C This subroutine uses simulated WHILE loop(s). PARAMETER(MaxR=10) INTEGER C,K,N,Row,RowK REAL A,B,Sum,X DIMENSION A(1:MaxR,1:MaxR),B(1:MaxR),Row(1:MaxR),X(1:MaxR) DO 10 K=1,N IF (A(Row(K),K).EQ.0) RETURN 10 CONTINUE X(1)=B(Row(1)) DO 30 K=2,N Sum=0 RowK=Row(K) DO 20 C=1,(K-1) Sum=Sum+A(RowK,C)*X(C) 20 CONTINUE X(K)=B(RowK)-Sum 30 CONTINUE X(N)=X(N)/A(Row(N),N) DO 50 K=(N-1),1,-1 Sum=0 RowK=Row(K) DO 40 C=(K+1),N Sum=Sum+A(RowK,C)*X(C) 40 CONTINUE X(K)=(X(K)-Sum)/A(RowK,K) 50 CONTINUE RETURN END SUBROUTINE INPUTS(A,A1,B,N,InRC) PARAMETER(MaxR=10) INTEGER Col,I,InRC,N,Row REAL A,B CHARACTER Ans*1,Ach*1 DIMENSION A(1:MaxR,1:MaxR),A1(1:MaxR,1:MaxR),B(1:MaxR) DO 10 I=1,12 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)'SOLUTION OF A LINEAR SYSTEM:' WRITE(9,*)' ' WRITE(9,*)' A*X = B' WRITE(9,*)' ' WRITE(9,*)'THE TRIANGULAR FACTORIZATION L*U = P*A IS CONSTRUCTED' WRITE(9,*)' ' WRITE(9,*)'FIRST, THE SOLUTION Y TO L*Y = P*B IS FOUND,' WRITE(9,*)' ' WRITE(9,*)'SECOND, THE SOLUTION X TO U*X = Y IS FOUND.' WRITE(9,*)' ' WRITE(9,*)'A IS AN N BY N NONSINGULAR MATRIX.' WRITE(9,*)' ' WRITE(9,*)'B IS AN N DIMENSIONAL VECTOR OF CONSTANTS.' WRITE(9,*)' ' WRITE(9,*)'X IS THE N DIMENSIONAL SOLUTION VECTOR.' WRITE(9,*)' ' WRITE(9,*)'ENTER THE NUMBER OF EQUATIONS N = ' READ(9,*) N InRC=0 WRITE(9,*)'WANT TO INPUT COLUMNS ? ' READ(9,'(A)') Ans IF (Ans.EQ.'Y' .OR. Ans.EQ.'y') THEN InRC=1 ENDIF WRITE(9,*)' ' WRITE(9,*)'ENTER THE COEFFICIENT MATRIX A:' CALL MATIN('A',A,N,InRC) DO 30 Row=1,N DO 20 Col=1,N A1(Row,Col)=A(Row,Col) 20 CONTINUE 30 CONTINUE RETURN END SUBROUTINE MATIN(Ach,A,N,InRC) PARAMETER(MaxR=10) INTEGER Col,InRC,N,Row REAL A CHARACTER Ach*1 DIMENSION A(1:MaxR,1:MaxR) DO 20 Row=1,N DO 10 Col=1,N A(Row,Col)=0 10 CONTINUE 20 CONTINUE WRITE(9,*)'GIVE THE ELEMENTS OF THE MATRIX ',Ach WRITE(9,*)' ' IF (InRC.EQ.0) THEN DO 40 Row=1,N WRITE(9,*)'GIVE THE ELEMENTS OF ROW ',Row WRITE(9,*)' ' DO 30 Col=1,N WRITE(9,*) Ach,'(',Row,',',Col,') = ' READ(9,*) A(Row,Col) 30 CONTINUE 40 CONTINUE ELSE DO 60 Col=1,N WRITE(9,*)'GIVE THE ELEMENTS OF COLUMN ',Col WRITE(9,*)' ' DO 50 Row=1,N WRITE(9,*) Ach,'(',Row,',',Col,') = ' READ(9,*) A(Row,Col) 50 CONTINUE 60 CONTINUE ENDIF RETURN END SUBROUTINE VECIN(Bch,B,N) PARAMETER(MaxR=10) INTEGER N,Row REAL A CHARACTER Bch*1 DIMENSION B(1:MaxR) WRITE(9,*)'GIVE THE ELEMENTS OF THE COLUMN VECTOR ',Bch WRITE(9,*)' ' DO 10 Row=1,N WRITE(9,*) Bch,'(',Row,') = ' READ(9,*) B(Row) 10 CONTINUE RETURN END SUBROUTINE MATPRINT(A,N) PARAMETER(MaxR=10) INTEGER Col,N,Row REAL A DIMENSION A(1:MaxR,1:MaxR) DO 10 Row=1,N WRITE(9,999) ( A(Row,Col), Col=1,N ) 10 CONTINUE 999 FORMAT(1X,5F15.7,4X) RETURN END SUBROUTINE VECPRINT(Bch,Xch,B,X,N) PARAMETER(MaxR=10) INTEGER N,Row REAL B,X CHARACTER Bch*1,Xch*1 DIMENSION B(1:MaxR),X(1:MaxR) WRITE(9,*) Bch,' COEFFICIENT VECTOR ',Xch, +' SOLUTION VECTOR' WRITE(9,*) DO 10 Row=1,N WRITE(9,*) Bch,'(',Row,') =',B(Row),' ', +Xch,'(',Row,') =',X(Row) 10 CONTINUE RETURN END SUBROUTINE RESULTS(A,B,X,N,DET) PARAMETER(MaxR=10) INTEGER I,N REAL A,B,DET,X CHARACTER Bch*1,Xch*1 DIMENSION A(1:MaxR,1:MaxR),B(1:MaxR),X(1:MaxR) DO 10 I=1,15 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)'COMPUTATION OF THE SOLUTION FOR THE LINEAR SYSTEM A*X = + B.' IF (DET.NE.0) THEN WRITE(9,*)'THE TRIANGULAR FACTORIZATION L*U = P*A WAS CONSTRUCTED' WRITE(9,*)'FIRST, THE SOLUTION Y TO L*Y = P*B WAS FOUND,' WRITE(9,*)'SECOND, THE SOLUTION X TO U*X = Y WAS FOUND.' ENDIF WRITE(9,*) CALL MATPRINT(A,N) IF (DET.EQ.0) THEN WRITE(9,*) WRITE(9,*)'THE MATRIX IS SINGULAR.' WRITE(9,*) WRITE(9,*)'A ZERO PIVOT ELEMENT WAS ENCOUNTERED.' WRITE(9,*)' ' WRITE(9,*)'THE MATRIX DOES NOT HAVE TRIANGULAR FACTORIZATION.' WRITE(9,*)' ' WRITE(9,*)'THE METHOD DOES NOT APPLY.' ELSE WRITE(9,*) CALL VECPRINT('B','X',A,X,N) ENDIF WRITE(9,*)' ' WRITE(9,*)'THE DETERMINANT`S VALUE IS DET A = ',DET RETURN END PROGRAM JACOBI C NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994 C To accompany the text: C NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992 C Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A. C This free software is complements of the author. C C Algorithm 3.4 (Jacobi Iteration). C Section 3.7, Iterative Methods for Linear Systems, Page 186 C PARAMETER(MaxR=25,Max=99,Tol=1E-6) INTEGER Cond,Count,InRC,N,Row REAL A,A1,B,DET,Sep,X CHARACTER Ans*1,Ach*1,Bch*1,Xch*1 DIMENSION A(1:MaxR,1:MaxR),A1(1:MaxR,1:MaxR) DIMENSION B(1:MaxR),Row(1:MaxR),X(1:MaxR) 10 CALL INPUTS(A,A1,B,N,InRC) WRITE(9,*)'ENTER the starting vector.' Xch='X' CALL VECIN(Xch,X,N) CALL ITERATE(A,B,X,N,Max,Tol,Sep,Cond,Count) CALL RESULTS(A,B,X,N,Tol,Sep,Cond,Count) WRITE(9,*)' ' WRITE(9,*)'Want to solve another linear system ? ' READ(9,'(A)') Ans IF (Ans.EQ.'Y' .OR. Ans.EQ.'y') GOTO 10 END SUBROUTINE ITERATE(A,B,X,N,Max,Tol,Sep,Cond,Count) PARAMETER(MaxR=25) INTEGER C,Cond,Count,J,K,N,Max,R REAL A,A1,B,DET,X,Xold,Row,Sep,Sum,T,Tol DIMENSION A(1:MaxR,1:MaxR) DIMENSION B(1:MaxR),X(1:MaxR),Xold(1:MaxR) Sep=1 K=0 Cond=1 DO J=1,N Xold(J)=X(J) ENDDO DO R=1,N Row=0 DO C=1,N Row=Row+ABS(A(R,C)) ENDDO IF (Row.GE.2*ABS(A(R,R))) Cond=0 ENDDO IF (Cond.EQ.0) RETURN WRITE(9,1000) K, ( Xold(J), J=1,N ) WHILE ((K.LT.Max).AND.(Sep.GT.Tol)) DO R=1,N Sum=B(R) DO C=1,N IF (C.NE.R) Sum=Sum-A(R,C)*Xold(C) ENDDO X(R)=Sum/A(R,R) ENDDO Sep=0 DO J=1,N Sep=Sep+ABS(X(J)-Xold(J)) Xold(J)=X(J) ENDDO K=K+1 WRITE(9,1000) K, ( Xold(J), J=1,N ) REPEAT Count=K 1000 FORMAT(I2,4X,5F15.7) PAUSE RETURN END SUBROUTINE XITERATE(A,B,X,N,Max,Tol,Sep,Cond,Count) C This subroutine uses simulated WHILE loop(s). PARAMETER(MaxR=25) INTEGER C,Cond,Count,J,K,N,Max,R REAL A,A1,B,DET,X,Xold,Row,Sep,Sum,T,Tol CHARACTER Ans*1,Ach*1,Bch*1,Xch*1 DIMENSION A(1:MaxR,1:MaxR) DIMENSION B(1:MaxR),X(1:MaxR),Xold(1:MaxR) Sep=1 K=0 Cond=1 WRITE(9,*)'ENTER the starting vector.' Xch='X' CALL VECIN(Xch,X,N) DO 10 J=1,N Xold(J)=X(J) 10 CONTINUE DO 30 R=1,N Row=0 DO 20 C=1,N Row=Row+ABS(A(R,C)) 20 CONTINUE IF (Row.GE.2*ABS(A(R,R))) Cond=0 30 CONTINUE IF (Cond.EQ.0) RETURN WRITE(9,1000) K, ( Xold(J), J=1,N ) 50 IF ((K.LT.Max).AND.(Sep.GT.Tol)) THEN DO 70 R=1,N Sum=B(R) DO 60 C=1,N IF (C.NE.R) Sum=Sum-A(R,C)*Xold(C) 60 CONTINUE X(R)=Sum/A(R,R) 70 CONTINUE Sep=0 DO 80 J=1,N Sep=Sep+ABS(X(J)-Xold(J)) 80 CONTINUE DO 90 J=1,N Xold(J)=X(J) 90 CONTINUE K=K+1 WRITE(9,1000) K, ( Xold(J), J=1,N ) GOTO 50 ENDIF Count=K 1000 FORMAT(I2,4X,5F15.7) PAUSE RETURN END SUBROUTINE INPUTS(A,A1,B,N,InRC) PARAMETER(MaxR=25) INTEGER Col,I,InRC,N,Row REAL A,B CHARACTER Ans*1,Ach*1,Bch*1 DIMENSION A(1:MaxR,1:MaxR),A1(1:MaxR,1:MaxR),B(1:MaxR) DO 10 I=1,12 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)' ' WRITE(9,*)' Jacobi iteration is is performed to find t +he solution to the linear' WRITE(9,*)' ' WRITE(9,*)' system A*X = B. The matrix must be diagonally d +ominant. Starting with' WRITE(9,*)' ' WRITE(9,*)' (k) (k) + (k) ' WRITE(9,*)' X = (0,0,...,0), a sequence { X = (x ,x ,. +..,x ) } is generated' WRITE(9,*)' 0 k 1 2 + N ' WRITE(9,*)' ' WRITE(9,*)' which converges to the solution.' WRITE(9,*)' ' WRITE(9,*)' (k+1)' WRITE(9,*)' the j-th coordinate x of X is computed + using the formula:' WRITE(9,*)' j k+1' WRITE(9,*)' ' WRITE(9,*)' (k+1) (k+1) (k+1)' WRITE(9,*)' x = [b - a x -...- a x ' WRITE(9,*)' j j j,1 1 j,j-1 j-1 ' WRITE(9,*)' ' WRITE(9,*)' (k) (k) ' WRITE(9,*)' - a x -...- a x ]/a + for j=1,2,...,N. ' WRITE(9,*)' j,j+1 j+1 j,N N j,j' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)' give the dimension N = ' READ(9,*) N InRC=0 WRITE(9,*)'Want to input columns ? ' READ(9,'(A)') Ans IF (Ans.EQ.'Y' .OR. Ans.EQ.'y') THEN InRC=1 ENDIF WRITE(9,*)' ' WRITE(9,*)'ENTER the coefficient matrix A:' Ach='A' CALL MATIN(Ach,A,N,InRC) WRITE(9,*)'The matrix A is:' CALL MATPRINT(A,N) Bch='B' CALL VECIN(Bch,B,N) RETURN END SUBROUTINE MATIN(Ach,A,N,InRC) PARAMETER(MaxR=25) INTEGER Col,InRC,N,Row REAL A CHARACTER Ach*1 DIMENSION A(1:MaxR,1:MaxR) DO 20 Row=1,N DO 10 Col=1,N A(Row,Col)=0 10 CONTINUE 20 CONTINUE WRITE(9,*)'Give the elements of the matrix ',Ach WRITE(9,*)' ' IF (InRC.EQ.0) THEN DO 40 Row=1,N WRITE(9,*)'ENTER the elements of row ',Row WRITE(9,*)' ' DO 30 Col=1,N WRITE(9,*) Ach,'(',Row,',',Col,') = ' READ(9,*) A(Row,Col) 30 CONTINUE 40 CONTINUE ELSE DO 60 Col=1,N WRITE(9,*)'ENTER the elements of column ',Col WRITE(9,*)' ' DO 50 Row=1,N WRITE(9,*) Ach,'(',Row,',',Col,') = ' READ(9,*) A(Row,Col) 50 CONTINUE 60 CONTINUE ENDIF RETURN END SUBROUTINE VECIN(Bch,B,N) PARAMETER(MaxR=25) INTEGER N,Row REAL B CHARACTER Bch*1 DIMENSION B(1:MaxR) WRITE(9,*)'ENTER the elements of the column vector ',Bch WRITE(9,*)' ' DO 10 Row=1,N WRITE(9,*) Bch,'(',Row,') = ' READ(9,*) B(Row) 10 CONTINUE RETURN END SUBROUTINE MATPRINT(A,N) PARAMETER(MaxR=25) INTEGER Col,N,Row REAL A DIMENSION A(1:MaxR,1:MaxR) DO 10 Row=1,N WRITE(9,999) ( A(Row,Col), Col=1,N ) 10 CONTINUE 999 FORMAT(1X,5F15.7,4X) RETURN END SUBROUTINE VECPRINT(Bch,Xch,B,X,N) PARAMETER(MaxR=25) INTEGER N,Row REAL B,X CHARACTER Bch*1,Xch*1 DIMENSION B(1:MaxR),X(1:MaxR) WRITE(9,*) Bch,' COEFFICIENT VECTOR ',Xch, +' SOLUTION VECTOR' WRITE(9,*) DO 10 Row=1,N WRITE(9,*) Bch,'(',Row,') =',B(Row),' ',Xch,'(',Row,') =' +,X(Row) 10 CONTINUE RETURN END SUBROUTINE RESULTS(A,B,X,N,Tol,Sep,Cond,Count) PARAMETER(MaxR=25) INTEGER C,Cond,Count,K,N,R REAL A,B,X,Tol,Sep CHARACTER Bch*1,Xch*1 DIMENSION A(1:MaxR,1:MaxR),B(1:MaxR),X(1:MaxR) DO 10 I=1,15 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)' ' WRITE(9,*)'The matrix A(I,J) is:' IF (N.LT.6) THEN WRITE(9,*)' ' ENDIF CALL MATPRINT(A,N) WRITE(9,*)' ' WRITE(9,*)'The vector B is:' IF (N.LT.6) THEN WRITE(9,*)' ' ENDIF DO 20 K=1,N WRITE(9,*)' B(',K,') = ',B(K) 20 CONTINUE WRITE(9,*)' ' IF (Cond.EQ.1) THEN IF (Sep.LT.Tol) THEN WRITE(9,*)' Jacobi iteration solved A*X = B.' WRITE(9,*)'It took',Count,' iterations to find the solution:' ELSE WRITE(9,*)' Jacobi iteration did not converge.' WRITE(9,*)'The status after ',Count,' iterations is' ENDIF IF (N.LT.6) THEN WRITE(9,*)' ' ENDIF WRITE(9,*)' ' DO 30 K=1,N WRITE(9,*)' X(',K,') = ',X(K) 30 CONTINUE WRITE(9,*)' ' WRITE(9,*)' +-',Sep,' is the estimated accuracy.' ELSE WRITE(9,*)' ' WRITE(9,*)'The matrix is not diagonally dominant.' WRITE(9,*)' ' WRITE(9,*)'Jacobi iteration cannot be used.' ENDIF RETURN END PROGRAM GAUSSIT C NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994 C To accompany the text: C NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992 C Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A. C This free software is complements of the author. C C Algorithm 3.5 (Gauss-Seidel Iteration). C Section 3.7, Iterative Methods for Linear Systems, Page 187 C PARAMETER(MaxR=25,Max=99,Tol=1E-7) INTEGER Cond,Count,InRC,N,Row REAL A,A1,B,DET,Sep,X CHARACTER Ans*1,Ach*1,Bch*1,Xch*1 DIMENSION A(1:MaxR,1:MaxR),A1(1:MaxR,1:MaxR) DIMENSION B(1:MaxR),Row(1:MaxR),X(1:MaxR) 10 CALL INPUTS(A,A1,B,N,InRC) WRITE(9,*)'ENTER the starting vector.' Xch='X' CALL VECIN(Xch,X,N) CALL ITERATE(A,B,X,N,Max,Tol,Sep,Cond,Count) CALL RESULTS(A,B,X,N,Tol,Sep,Cond,Count) WRITE(9,*)' ' WRITE(9,*)'Want to solve another linear system ? ' READ(9,'(A)') Ans IF (Ans.EQ.'Y' .OR. Ans.EQ.'y') GOTO 10 END SUBROUTINE ITERATE(A,B,X,N,Max,Tol,Sep,Cond,Count) PARAMETER(MaxR=25) INTEGER C,Cond,Count,J,K,N,Max,R REAL A,A1,B,DET,X,Xold,Row,Sep,Sum,T,Tol CHARACTER Ans*1,Ach*1,Bch*1,Xch*1 DIMENSION A(1:MaxR,1:MaxR) DIMENSION B(1:MaxR),X(1:MaxR),Xold(1:MaxR) Sep=1 K=0 Cond=1 DO J=1,N Xold(J)=X(J) ENDDO DO R=1,N Row=0 DO C=1,N Row=Row+ABS(A(R,C)) ENDDO IF (Row.GE.2*ABS(A(R,R))) Cond=0 ENDDO IF (Cond.EQ.0) RETURN WRITE(9,1000) K, ( Xold(J), J=1,N ) WHILE ((K.LT.Max).AND.(Sep.GT.Tol)) DO R=1,N Sum=B(R) DO C=1,N IF (C.NE.R) Sum=Sum-A(R,C)*X(C) ENDDO X(R)=Sum/A(R,R) ENDDO Sep=0 DO J=1,N Sep=Sep+ABS(X(J)-Xold(J)) Xold(J)=X(J) ENDDO K=K+1 WRITE(9,1000) K, ( Xold(J), J=1,N ) REPEAT Count=K 1000 FORMAT(I2,4X,5F15.7) PAUSE RETURN END SUBROUTINE XITERATE(A,B,X,N,Max,Tol,Sep,Cond,Count) C This subroutine uses simulated WHILE loop(s). PARAMETER(MaxR=25) INTEGER C,Cond,Count,J,K,N,Max,R REAL A,A1,B,DET,X,Xold,Row,Sep,Sum,T,Tol CHARACTER Ans*1,Ach*1,Bch*1,Xch*1 DIMENSION A(1:MaxR,1:MaxR) DIMENSION B(1:MaxR),X(1:MaxR),Xold(1:MaxR) Sep=1 K=0 Cond=1 DO 10 J=1,N Xold(J)=X(J) 10 CONTINUE DO 30 R=1,N Row=0 DO 20 C=1,N Row=Row+ABS(A(R,C)) 20 CONTINUE IF (Row.GE.2*ABS(A(R,R))) Cond=0 30 CONTINUE IF (Cond.EQ.0) RETURN DO 40 J=1,10 WRITE(9,*) ' ' 40 CONTINUE WRITE(9,1000) K, ( Xold(J), J=1,N ) 50 IF ((K.LT.Max).AND.(Sep.GT.Tol)) THEN DO 70 R=1,N Sum=B(R) DO 60 C=1,N IF (C.NE.R) Sum=Sum-A(R,C)*X(C) 60 CONTINUE X(R)=Sum/A(R,R) 70 CONTINUE Sep=0 DO 80 J=1,N Sep=Sep+ABS(X(J)-Xold(J)) Xold(J)=X(J) 80 CONTINUE K=K+1 WRITE(9,1000) K, ( Xold(J), J=1,N ) GOTO 50 ENDIF Count=K 1000 FORMAT(I2,4X,5F15.7) PAUSE RETURN END SUBROUTINE INPUTS(A,A1,B,N,InRC) PARAMETER(MaxR=25) INTEGER Col,I,InRC,N,Row REAL A,B CHARACTER Ans*1,Ach*1,Bch*1 DIMENSION A(1:MaxR,1:MaxR),A1(1:MaxR,1:MaxR),B(1:MaxR) DO 10 I=1,12 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)' ' WRITE(9,*)' Gauss-Seidel iteration is is performed to find t +he solution to the linear' WRITE(9,*)' ' WRITE(9,*)' system A*X = B. The matrix must be diagonally d +ominant. Starting with' WRITE(9,*)' ' WRITE(9,*)' (k) (k) + (k) ' WRITE(9,*)' X = (0,0,...,0), a sequence { X = (x ,x ,. +..,x ) } is generated' WRITE(9,*)' 0 k 1 2 + N ' WRITE(9,*)' ' WRITE(9,*)' which converges to the solution.' WRITE(9,*)' ' WRITE(9,*)' (k+1)' WRITE(9,*)' the j-th coordinate x of X is computed + using the formula:' WRITE(9,*)' j k+1' WRITE(9,*)' ' WRITE(9,*)' (k+1) (k+1) (k+1)' WRITE(9,*)' x = [b - a x -...- a x ' WRITE(9,*)' j j j,1 1 j,j-1 j-1 ' WRITE(9,*)' ' WRITE(9,*)' (k) (k) ' WRITE(9,*)' - a x -...- a x ]/a + for j=1,2,...,N. ' WRITE(9,*)' j,j+1 j+1 j,N N j,j' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)' give the dimension N = ' READ(9,*) N InRC=0 WRITE(9,*)'Want to input columns ? ' READ(9,'(A)') Ans IF (Ans.EQ.'Y' .OR. Ans.EQ.'y') THEN InRC=1 ENDIF WRITE(9,*)' ' WRITE(9,*)'ENTER the coefficient matrix A:' Ach='A' CALL MATIN(Ach,A,N,InRC) WRITE(9,*)'The matrix A is:' CALL MATPRINT(A,N) Bch='B' CALL VECIN(Bch,B,N) RETURN END SUBROUTINE MATIN(Ach,A,N,InRC) PARAMETER(MaxR=25) INTEGER Col,InRC,N,Row REAL A CHARACTER Ach*1 DIMENSION A(1:MaxR,1:MaxR) DO 20 Row=1,N DO 10 Col=1,N A(Row,Col)=0 10 CONTINUE 20 CONTINUE WRITE(9,*)'Give the elements of the matrix ',Ach WRITE(9,*)' ' IF (InRC.EQ.0) THEN DO 40 Row=1,N WRITE(9,*)'ENTER the elements of row ',Row WRITE(9,*)' ' DO 30 Col=1,N WRITE(9,*) Ach,'(',Row,',',Col,') = ' READ(9,*) A(Row,Col) 30 CONTINUE 40 CONTINUE ELSE DO 60 Col=1,N WRITE(9,*)'ENTER the elements of column ',Col WRITE(9,*)' ' DO 50 Row=1,N WRITE(9,*) Ach,'(',Row,',',Col,') = ' READ(9,*) A(Row,Col) 50 CONTINUE 60 CONTINUE ENDIF RETURN END SUBROUTINE VECIN(Bch,B,N) PARAMETER(MaxR=25) INTEGER N,Row REAL B CHARACTER Bch*1 DIMENSION B(1:MaxR) WRITE(9,*)'ENTER the elements of the column vector ',Bch WRITE(9,*)' ' DO 10 Row=1,N WRITE(9,*) Bch,'(',Row,') = ' READ(9,*) B(Row) 10 CONTINUE RETURN END SUBROUTINE MATPRINT(A,N) PARAMETER(MaxR=25) INTEGER Col,N,Row REAL A DIMENSION A(1:MaxR,1:MaxR) DO 10 Row=1,N WRITE(9,999) ( A(Row,Col), Col=1,N ) 10 CONTINUE 999 FORMAT(1X,5F15.7,4X) RETURN END SUBROUTINE VECPRINT(Bch,Xch,B,X,N) PARAMETER(MaxR=25) INTEGER N,Row REAL B,X CHARACTER Bch*1,Xch*1 DIMENSION B(1:MaxR),X(1:MaxR) WRITE(9,*) Bch,' COEFFICIENT VECTOR ',Xch, +' SOLUTION VECTOR' WRITE(9,*) DO 10 Row=1,N WRITE(9,*) Bch,'(',Row,') =',B(Row),' ',Xch,'(',Row,') =' +,X(Row) 10 CONTINUE RETURN END SUBROUTINE RESULTS(A,B,X,N,Tol,Sep,Cond,Count) PARAMETER(MaxR=25) INTEGER C,Cond,Count,K,N,R REAL A,B,X,Tol,Sep CHARACTER Bch*1,Xch*1 DIMENSION A(1:MaxR,1:MaxR),B(1:MaxR),X(1:MaxR) DO 10 I=1,15 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)' ' WRITE(9,*)'The matrix A(I,J) is:' IF (N.LT.6) THEN WRITE(9,*)' ' ENDIF CALL MATPRINT(A,N) WRITE(9,*)' ' WRITE(9,*)'The vector B is:' IF (N.LT.6) THEN WRITE(9,*)' ' ENDIF DO 20 K=1,N WRITE(9,*)' B(',K,') = ',B(K) 20 CONTINUE WRITE(9,*)' ' IF (Cond.EQ.1) THEN IF (Sep.LT.Tol) THEN WRITE(9,*)'Gauss-Seidel iteration solved A*X = B.' WRITE(9,*)'It took',Count,' iterations to find the solution:' ELSE WRITE(9,*)'Gauss-Seidel iteration did not converge.' WRITE(9,*)'The status after ',Count,' iterations is' ENDIF IF (N.LT.6) THEN WRITE(9,*)' ' ENDIF WRITE(9,*)' ' DO 30 K=1,N WRITE(9,*)' X(',K,') = ',X(K) 30 CONTINUE WRITE(9,*)' ' WRITE(9,*)' +-',Sep,' is the estimated accuracy.' ELSE WRITE(9,*)' ' WRITE(9,*)'The matrix is not diagonally dominant.' WRITE(9,*)' ' WRITE(9,*)'Gauss-Seidel iteration cannot be used.' ENDIF RETURN END