module precision_module
  implicit none
  integer, parameter :: long = selected_real_kind(15,307)
end module precision_module
module ge_module
  implicit none

contains
  subroutine gaussian_elimination(a,n,b,x,singular)
! routine to solve a system ax=b
! using gaussian elimination
! with partial pivoting
! the code is based on the linpack routines
! sgefa and sgesl
! and operates on columns rather than rows!
    use precision_module
    implicit none
! matrix a and vector b are over-written
! arguments
    integer, intent (in) :: n
    real (long), intent (inout) :: a(:,:), b(:)
    real (long), intent (out) :: x(:)
    logical, intent (out) :: singular
! local variables
    integer :: i, j, k, pivot_row
    real (long) :: pivot, sum, element
    real (long), parameter :: eps = 1.e-13_long
! work through the matrix column by column
    do k = 1, n - 1
! find largest element in column k for pivot
      pivot_row = maxval(maxloc(abs(a(k:n,k)))) + k - 1
! test to see if a is singular
! if so return to main program
      if (abs(a(pivot_row,k))<=eps) then
        singular = .true.
        return
      else
        singular = .false.
      end if
! exchange elements in column k if largest is
! not on the diagonal
      if (pivot_row/=k) then
        element = a(pivot_row,k)
        a(pivot_row,k) = a(k,k)
        a(k,k) = element
        element = b(pivot_row)
        b(pivot_row) = b(k)
        b(k) = element
      end if
! compute multipliers
! elements of column k below diagonal
! are set to these multipliers for use
! in elimination later on
      a(k+1:n,k) = a(k+1:n,k)/a(k,k)
! row elimination performed by columns for efficiency
      do j = k + 1, n
        pivot = a(pivot_row,j)
        if (pivot_row/=k) then
! swap if pivot row is not k
          a(pivot_row,j) = a(k,j)
          a(k,j) = pivot
        end if
        a(k+1:n,j) = a(k+1:n,j) - pivot*a(k+1:n,k)
      end do
! apply same operations to b
      b(k+1:n) = b(k+1:n) - a(k+1:n,k)*b(k)
    end do
! backward substitution
    do i = n, 1, -1
      sum = 0.0
      do j = i + 1, n
        sum = sum + a(i,j)*x(j)
      end do
      x(i) = (b(i)-sum)/a(i,i)
    end do
  end subroutine gaussian_elimination
end module ge_module
program ch2504
  use precision_module
  use ge_module
  implicit none
  integer :: i, n
  real (long), allocatable :: a(:,:), b(:), x(:)
  logical :: singular

  print *, 'number of equations?'
  read *, n
  allocate (a(1:n,1:n),b(1:n),x(1:n))
  do i = 1, n
    print *, 'input elements of row ', i, ' of a'
    read *, a(i,1:n)
    print *, 'input element ', i, ' of b'
    read *, b(i)
  end do
  call gaussian_elimination(a,n,b,x,singular)
  if (singular) then
    print *, 'matrix is singular'
  else
    print *, 'solution x:'
    print *, x(1:n)
  end if
end program ch2504