!> Program for calculating linear equations via Gauss elimination technique.
!!
!! The program first transforms the provided coefficient matrix to triangluar
!! form, and calculates the result by backsubstitution. If during the
!! transformation to triangular form a leading zero is obtained (equations
!! linearly dependent), the program stops.
program GaussElimination
  implicit none

  !! precision
  integer, parameter :: dp = selected_real_kind(14, 300)

  !! Minimal leading element
  real(dp), parameter :: minpivot = 1e-10_dp
  
  real(dp), allocatable :: aa(:,:), bb(:), xx(:)
  integer :: nvar
  real(dp) :: tmp
  integer :: ii, jj, imax

  write(*,"(A)", advance="no") "Enter number of variables: "
  read (*,*) nvar

  allocate(aa(nvar, nvar))
  allocate(bb(nvar))
  allocate(xx(nvar))

  ! Reading in coefficients
  ! Matrix coefficients are stored row wise.
  write(*,"(/,A)") "Enter coefficients (A):"
  do ii = 1, nvar
    do jj = 1, nvar
      write(*,"(A,I0,A,I0,A)", advance="no") " A(", ii, ",", jj, ")= "
      read (*, *) aa(ii, jj)
    end do
  end do

  ! Reading in RHS.
  write(*,"(/,A)") "Enter right hand side (b):"
  do ii = 1, nvar
    write(*,"(A,I0,A)", advance="no") " B(", ii, ")= "
    read (*,*) bb(ii)
  end do

  !! Elimination: aa to lower triangle form
  lpElim: do ii = 1, nvar
    !! Look for biggest element
    imax = maxloc(abs(aa(ii:,ii)), dim=1) + ii - 1
    if (abs(aa(imax,ii)) < minpivot) then
      write(*,*) "Sorry, couldn't proceed, equations probably linearly &
          &dependent!"
      stop
    end if

    if (imax /= ii) then
      !! Exchange rows to have the one with the biggest leading coefficient
      !! as the first
      xx = aa(ii,:)
      aa(ii,:) = aa(imax,:)
      aa(imax,:) = xx
      !! Do the same exchange in b as well.
      tmp = bb(ii)
      bb(ii) = bb(imax)
      bb(imax) = tmp
    end if
      
    !! Subtract actual column from all subsequent columns to set leading
    !! elements to zero.
    do jj = ii + 1, nvar
      tmp = aa(jj,ii) / aa(ii,ii)
      aa(jj,ii:) = aa(jj,ii:) - tmp * aa(ii,ii:)
      bb(jj) = bb(jj) - tmp * bb(ii)
    end do
  end do lpElim

  !! Back substitution
  do ii = nvar, 1, -1
    if (ii == nvar) then
      tmp = bb(ii)
    else
      tmp = bb(ii) - dot_product(aa(ii,ii+1:), xx(ii+1:))
    end if
    xx(ii) = tmp / aa(ii, ii)
  end do

  write(*,"(/,A)") "The components of the solution vector (x):"
  do ii = 1, nvar
    write(*,"(A,I0,A,E14.6)") " x(", ii, ")= ", xx(ii)
  end do

  deallocate(aa, bb, xx)

end program GaussElimination
