! Figure 6.9/10
module linear
   public :: solve
   integer, parameter, public :: kind10=selected_real_kind(10)
 
contains
   subroutine solve(a, piv_tol, b, ok)
   ! arguments
      real(kind=kind10), intent(inout), dimension(:,:) :: a
                       ! The matrix a.
      real(kind=kind10), intent(in) :: piv_tol 
                       ! Smallest acceptable pivot.
      real(kind=kind10), 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=kind10), dimension(size(b)) :: row 
                       ! Automatic array needed for workspace; 
                       ! size is described in Section 8.11.2.
      real(kind=kind10) :: 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 (using maxval just to 
!     obtain a scalar).
         i = maxval(maxloc(abs(a(j:n, j)))) + j - 1
            ! maxval and maxloc are in Sections 8.10.1 and 8.13.
         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=kind10), allocatable, dimension(:,:) :: a
   real(kind=kind10), allocatable, dimension(:) ::  b
   logical :: ok
   print *, " Matrix order?"
   read *,  n
   allocate ( a(n, n), b(n) )
   do i = 1, n
      write(unit=*,fmt="(a, i2, a)") " Elements of row ", i, " of a?"
           ! Edit descriptors are described in Section 9.11
      read *, a(i,:)
      write(unit=*,fmt="(a, i2, a)") " Component ", i, " of b?"
      read *, b(i)
   end do
   call solve(a, maxval(abs(a))*1.0e-10, b, ok)
   if (ok) then
      write(unit=*,fmt="(/,a,/,(5f12.4))") " Solution is", b
   else
      print *, " The matrix is singular"
   end if
end program main
