PROGRAM Linear_Systems
!-----------------------------------------------------------------------
! Program to solve a linear system using Gaussian elimination.
! Variables used are:
!     N        : number of equations and unknowns
!     I, J     : subscripts
!     LinSys   : matrix for the linear system
!     X        : solution
!     Singular : indicates if system is (nearly) singular
! 
! Input:   The number of equations, the coefficients, and the
!          constants of the linear system 
! Output:  The solution of the linear system or a message indicating
!          that the system is (nearly) singular
!-----------------------------------------------------------------------

  IMPLICIT NONE
  REAL, DIMENSION(:, :), ALLOCATABLE :: LinSys
  REAL, DIMENSION(:), ALLOCATABLE :: X
  INTEGER :: N, I, J
  LOGICAL :: Singular

  ! Read N and allocate N x (N + 1) matrix LinSys 
  ! and one-dimensional array X of size N

  WRITE (*, '(1X, A)', ADVANCE = "NO")  "Enter number of equations: "
  READ *, N
  ALLOCATE(LinSys(N, N+1), X(N))

  ! Read coefficients and constants

  DO I = 1, N
     PRINT *, "Enter coefficients and constant of equation", I, ":"
     READ *, (LinSys(I,J), J = 1, N + 1)
  END DO

! Use subroutine Gaussian_Elimination to find the solution,
! and then display the solution

  CALL Gaussian_Elimination(LinSys, N, X, Singular)
  IF (.NOT. Singular) THEN
     PRINT *, "Solution is:"
     DO I = 1, N
        PRINT '(1X, "X(", I2, ") =", F8.3)', I, X(I)
     END DO
  ELSE
     PRINT *, "Matrix is (nearly) singular"
  END IF

  DEALLOCATE(LinSys, X)

CONTAINS

  !-Gaussian_Elimination------------------------------------------------
  ! Subroutine to find solution of a linear system of N equations in N 
  ! unknowns using Gaussian elimination, provided a unique solution 
  ! exists.  The coefficients and constants of the linear system are
  ! stored in the matrix LinSys.  If the system is singular, Singular 
  ! is returned as true, and the solution X is undefined.  Local 
  ! variables used are:
  !   I, J       : subscripts
  !   Multiplier : multiplier used to eliminate an unknown 
  !   AbsPivot   : absolute value of pivot element 
  !   PivotRow   : row containing pivot element
  !   Epsilon    : a small positive real value ("almost zero")
  !   Temp       : used to interchange rows of matrix
  ! 
  ! Accepts: Two-dimensional array LinSys and integer N 
  ! Returns: One-dimensional array X and logical value Singular 
  !--------------------------------------------------------------------

    SUBROUTINE Gaussian_Elimination(LinSys, N, X, Singular)
  
    INTEGER, INTENT(IN) :: N
    REAL, DIMENSION(N, N+1), INTENT(INOUT) :: LinSys
    REAL, DIMENSION(N), INTENT(OUT) :: X
    LOGICAL, INTENT(OUT) :: Singular
    REAL, DIMENSION(N+1) :: Temp
    REAL :: AbsPivot, Multiplier
    REAL, PARAMETER :: Epsilon = 1E-7
    INTEGER :: PivotRow
    
    Singular = .FALSE.
    DO I = 1, N
  
       ! Locate pivot element
 
       AbsPivot = ABS(LinSys(I, I))
       PivotRow = I
       DO J = I + 1, N
          IF (ABS(LinSys(J,I)) > AbsPivot) THEN
             AbsPivot = ABS(LinSys(J,I))
             PivotRow = J
          END IF
       END DO
  
       ! Check if matrix is (nearly) singular
  		
       IF (AbsPivot < Epsilon) THEN
          Singular = .TRUE.
          RETURN
       END IF
  
       ! It isn't, so interchange rows PivotRow and I if necessary
  
       IF (PivotRow /= I) THEN
          Temp(I:N+1) = LinSys(I, I:N+1)
          LinSys(I, I:N+1) = LinSys(PivotRow, I:N+1)
          LinSys(PivotRow, I:N+1) = Temp(I:N+1)
       END IF
  
       ! Eliminate Ith unknown from equations I + 1, ..., N
  
       DO J = I + 1, N
          Multiplier = -LinSys(J,I) / LinSys(I,I)
          LinSys(J, I:N+1) = LinSys(J, I:N+1)  &
                            + Multiplier * LinSys(I, I:N+1)
       END DO
  
    END DO
  
    ! Find the solutions by back substitution
  
    X(N) = LinSys(N, N + 1) / LinSys(N,N)
    DO J = N - 1, 1, -1
       X(J) = (LinSys(J, N + 1) -   &
              SUM(LinSys(J, J+1:N) * X(J+1:N))) / LinSys(J,J)
    END DO
  
  END SUBROUTINE Gaussian_Elimination
  
END PROGRAM Linear_Systems