!    GAUSSIAN ELIMINATION DEMONSTRATION - SIMPLE VERSION
!
     PROGRAM P141
     IMPLICIT NONE
     REAL :: M1(5,6),SOL(5)
     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)
     IMPLICIT NONE
     REAL ,INTENT(IN OUT) :: M1(:,:),SOL(:)
     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
!
     PRINT *,'  This is Program P141 - Gaussian elimination'
     PRINT *,'PROGRAM IS READING DATA INTO ARRAYS'
     CALL MATSIN(M1,N)
     PRINT *,'SOLVING SYSTEM OF EQUATIONS'
     CALL GAUSS(M1,SOL,N)
     PRINT *,'SOLUTION:'
     PRINT 90,SOL
90   FORMAT(' | ',F8.3,' |')
     STOP
     END PROGRAM P141


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='P141.DAT')      ! UNIT=5 is the default input
!
!     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('  ',F8.3))
!
      RETURN
      END SUBROUTINE PRNMAT
!
SUBROUTINE GAUSS(M1,SOL,N)
     IMPLICIT NONE
     REAL ,INTENT(IN OUT) :: M1(:,:),SOL(:)
     INTEGER ,INTENT(IN OUT) :: N
! 
!  THIS ROUTINE PERFORMS GAUSSIAN ELIMINATION AND BACKSUBSTITUTION.
!  WE HAVE SIMPLIFIED THE PROBLEM BY NOT WORRYING ABOUT MATRICES
!  WITHOUT A SOLUTION. IF WE WANTED TO CONSIDER THAT CASE, ALL THE 
!  CODE THAT FOLLOWS WOULD STILL BE VALID, BUT WE WOULD HAVE TO ADD
!  MORE TESTS.
!
     REAL :: M2(:,:),TEMP
     INTEGER :: I,J,K
!
!     
!  INSTEAD OF MODIFYING THE ORIGINAL ARRAY, WE WILL PRODUCE A WORKING COPY
!  OF IT
!
     M2 = M1
!
L1:  DO I=1,N
!        
     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