!! leq - linear equations; Gauss elimination method
!! and backward substitution

module linear
integer, parameter, public :: kind=selected_real_kind(12)
public :: solve
contains
  subroutine solve(a, piv_tol, b, ok)
    ! arguments
    real(kind), intent(inout), dimension(:,:) :: a
    ! The matrix a.
    real(kind), intent(in) :: piv_tol
    ! Smallest acceptable pivot.
    real(kind), intent(inout), dimension(:) :: b
    ! The right-hand side vector on
    ! entry. Overwritten by the solution.
    logical, intent(out) :: ok
    ! True after a successful entry and false otherwise.
    ! Local variables
    integer :: i ! Row index.
    integer :: j ! Column index.
    integer :: n ! Matrix order.
    real(kind), dimension(size(b)) :: row
    ! Automatic array needed for workspace;
    real(kind) :: element ! Workspace variable.
    n = size(b)
    ok = size(a, 1) == n .and. size(a, 2) == n
    if (.not.ok) then
       return
    end if
    do j = 1, n
       ! Update elements in column j.
       do i = 1, j - 1
          a(i+1:n, j) = a(i+1:n, j) - a(i,j) * a(i+1:n, i)
       end do
       ! Find pivot and check its size
       i = maxloc(abs(a(j:n, j)), dim=1) + j - 1
       if (abs(a(i, j)) < piv_tol) then
          ok = .false.
          return
       end if
       ! If necessary, apply row interchange
       if (i/=j) then
          row = a(j, :); a(j, :) = a(i, :); a(i, :) = row
          element = b(j); b(j) = b(i); b(i) = element
       end if
       ! Compute elements j+1 : n of j-th column.
       a(j+1:n, j) = a(j+1:n, j)/a(j, j)
    end do
    ! Forward substitution
    do i = 1, n-1
       b(i+1:n) = b(i+1:n) - b(i)*a(i+1:n, i)
    end do
    ! Back-substitution
    do j = n, 1, -1
       b(j) = b(j)/a(j, j)
       b(1:j-1) = b(1:j-1) - b(j)*a(1:j-1, j)
    end do
  end subroutine solve
end module linear

program main
  use linear
  integer :: i, n
  real(kind), allocatable :: a(:, :), b(:), c(:), d(:,:)
  logical :: ok
  print *, ' Matrix order?'
  read *, n
  allocate ( a(n, n), b(n), c(n), d(n,n) )
!!$  do i = 1, n
!!$     write (*, '(a, i2, a)') ' Elements of row ', i, ' of a?'
!!$     read *, a(i,:)
!!$     write (*, '(a, i2, a)') ' Component ', i, ' of b?'
!!$     read *, b(i)
!!$  end do
  call random_number(a) ! destroyed
  call random_number(b) ! destroyed
  c = b
  d = a
  call solve(a, maxval(abs(a))*1.0e-10, b, ok)
  if (ok) then
     write (*, '(/,a,/,(5f12.4))') ' Solution is', b
     ! check
     h = norm2(c - matmul(d,b))
     write (*, '(/,a,f12.4)') ' Norm |ax-b|:', h
  else
     print *, ' The matrix is singular'
  end if
  
end program main
