Example 20.2

FORTRAN Version

!    GAUSSIAN ELIMINATION DEMONSTRATION - BETTER VERSION
!
     PROGRAM P142
     IMPLICIT NONE
     REAL :: M1(5,6),SOL(5),TOL
     INTEGER :: N
     INTERFACE
     SUBROUTINE MATSIN(M1,N)
     IMPLICIT NONE
     REAL ,INTENT(IN OUT) :: M1(:,:)
     INTEGER ,INTENT(IN OUT) :: N
     END SUBROUTINE MATSIN  
     SUBROUTINE GAUSS(M1,SOL,N,TOL)
     IMPLICIT NONE
     REAL ,INTENT(IN OUT) :: M1(:,:),SOL(:),TOL
     INTEGER ,INTENT(IN OUT) :: N
     SUBROUTINE PRNMAT(M3,N)
     IMPLICIT NONE
     REAL ,INTENT(IN OUT) :: M3(:,:)
     INTEGER ,INTENT(IN OUT) :: N
     END SUBROUTINE PRNMAT 
     END SUBROUTINE GAUSS
     END INTERFACE
!
     N=5  ! NUMBER OF EQUATIONS
     TOL=.001
!
     PRINT *,'  This is Program P142 - Gaussian elimination'
     PRINT *,'PROGRAM IS READING DATA INTO ARRAYS'
     CALL MATSIN(M1,N)
     PRINT *,'SOLVING SYSTEM OF EQUATIONS'
     CALL GAUSS(M1,SOL,N,TOL)
     IF(TOL<0) THEN
        PRINT *,'SYSTEM OF EQUATIONS HAS NO SINGLE SOLUTION'
        STOP
     ENDIF
     PRINT *,'SOLUTION:'
     PRINT 90,SOL
90   FORMAT(' | ',F8.3,' |')
     STOP
     END PROGRAM P142


SUBROUTINE MATSIN(M1,N)
     IMPLICIT NONE
     REAL ,INTENT(IN OUT) :: M1(:,:)
     INTEGER ,INTENT(IN OUT) :: N
     INTEGER :: I,J
!
!     Tell program where data for READ is coming from
      OPEN(UNIT=5, FILE='P142.DAT') 
!
!     READ IN M1
!     ONE ROW PER CARD
!
L1:   DO I=1,N
         READ 27,(M1(I,J),J=1,N+1)
      END DO L1
27   FORMAT(10(F5.2))
     RETURN
     END SUBROUTINE MATSIN
!
SUBROUTINE PRNMAT(M3,N)
     IMPLICIT NONE
     REAL ,INTENT(IN OUT) :: M3(:,:)
     INTEGER ,INTENT(IN OUT) :: N
     INTEGER :: I,J
!
!
L3:   DO I=1,N
         PRINT 202, (M3(I,J),J=1,N+1)
      END DO L3
202   FORMAT(10('  ',F7.3))
!
      RETURN
      END SUBROUTINE PRNMAT
!
SUBROUTINE GAUSS(M1,SOL,N,TOL)
     IMPLICIT NONE
     REAL ,INTENT(IN OUT) :: M1(:,:),SOL(:),TOL
     INTEGER ,INTENT(IN OUT) :: N
! 
!  THIS ROUTINE PERFORMS GAUSSIAN ELIMINATION AND BACKSUBSTITUTION.
!  IN THIS VERSION, WE CONSIDER THE CASE WHERE THE SYSTEM OF EQUATIONS
!  HAS NO SINGLE SOLUTION ( INFINITELY MANY OR NONE ).  TO MAKE IT EASY
!  TO CHECK FOR THIS POSSIBILITY, THE COMPUTER CHOOSES WHICH ROW HAS THE
!  GREATEST LEADING CO-EFFICIENT, AND USES THIS ROW IN THE ELIMINATION
!  PROCESS. IF IT CANNOT FIND A NON-ZERO ROW (ZERO WITHIN A TOLERANCE
!  SET BY THE ROUTINE WHICH CALLS THIS SUBROUTINE), IT RETURNS TO THE
!  CALLING ROUTINE WITH A TOL=-90 FLAG, AND THE USER IS TOLD THAT THE
!  SYSTEM HAS NO SOLUTION.
!
     REAL :: M2(:,:),TEMP,MAX
     INTEGER :: SWAP,I,J,K
!
!     
!  INSTEAD OF MODIFYING THE ORIGINAL ARRAY, WE WILL PRODUCE A WORKING COPY
!  OF IT
!
     M2 = M1
!
L1:  DO I=1,N
        MAX=-3.0
        DO K=I,N
           IF(ABS(M2(K,I))>MAX) THEN
              MAX=ABS(M2(K,I))
              SWAP=K
           ENDIF
        END DO
        IF(MAXI) THEN
           DO M=I,N+1
              TEMP=M2(I,M)
              M2(I,M)=M2(SWAP,M)
              M2(SWAP,M)=TEMP
           ENDDO
        ENDIF
!        
     L2: DO J=I+1,N
           TEMP=M2(J,I)/M2(I,I)   
       L3: DO K=I,N+1
              M2(J,K)= M2(J,K) - TEMP*M2(I,K)
           END DO L3
        END DO L2
     END DO L1
!
     PRINT *,'TRIANGULARIZED MATRIX'
     CALL PRNMAT(M2,5)
!
!  MATRIX IS NOW TRIANGULAR. USE BACKSUBSTITUTION TO SOLVE
!
     SOL(N)=M2(N,N+1)/M2(N,N)
L4:  DO I=N-1,1,-1     
        TEMP=0.0
    L5: DO K=N,I+1,-1
           TEMP=TEMP+M2(I,K)*SOL(K)
        END DO L5
        SOL(I)=(M2(I,N+1) - TEMP)/M2(I,I)
     END DO L4
     RETURN
     END SUBROUTINE GAUSS

45.3 23.0 63.7 2.1 54.6 3.4
64.6 3.3 75.7 25.3 45.8 74.3
45.7 9.0 2.3 67.2 34.9 23.4
54.2 2.9 25.9 21.3 4.2 9.1
5.8 64.3 91.3 43.4 21.6 43.9

[FTN90 Version 1.12 Copyright (c)SALFORD SOFTWARE LTD 1992  &  ]
[                   (c)THE NUMERICAL ALGORITHMS GROUP 1991,1992]
    NO ERRORS  [FTN90]
Program entered
   This is Program P142 - Gaussian elimination
 PROGRAM IS READING DATA INTO ARRAYS
 SOLVING SYSTEM OF EQUATIONS
 TRIANGULARIZED MATRIX
   64.600    3.370    5.720    5.340    5.870    4.300  
    0.000   20.637   59.689   -1.595    0.513   -2.615  
    0.000    0.000  -22.873   -3.031   -3.395   -1.801  
    0.000    0.000    0.000    4.708    2.637    4.893  
    0.000    0.000    0.000    0.000   -2.945   -0.279  
 SOLUTION:
 |   -0.024 |
 |    0.138 |
 |   -0.066 |
 |    0.986 |
 |    0.094 |
Fortran-90 STOP

Last modified: 25/07/97