module gauss

  implicit none
  integer,                            private :: n
  real, dimension (:,:), allocatable, private :: d
  real, dimension (:),   allocatable, private :: drow

contains

  subroutine gauss_init (size)
    integer, intent (in) :: size
    n = size
    allocate (d (n,n+1), drow (n+1))
  end subroutine gauss_init

  function gauss_solve (a,b,x) result (code)
    logical                                :: code
    real,    dimension (:,:), intent (in)  :: a
    real,    dimension (:),   intent (in)  :: b
    real,    dimension (:),   intent (out) :: x

    integer, dimension (1)                 :: mloc
    integer                                :: k

    d (:,1:n) = a
    d (:,n+1) = b

    ! triangularization phase
    do k = 1, n
       ! pivoting
       mloc = maxloc (abs(d(k:n,k)))
       drow (k:n+1) = d (k,k:n+1)
       d (k, k:n+1) = d (k-1+mloc(1), k:n+1)
       d (k-1+mloc(1), k:n+1) = drow (k:n+1)

       if (abs(d(k,k)) < epsilon (d(1,1))) then
          code = .false.
          return
       end if
       d (k,k:n+1)       = d (k,k:n+1) / d (k,k)

       d (k+1:n,k+1:n+1) = &
       d (k+1:n,k+1:n+1) - spread (d (k,k+1:n+1), 1, n-k  ) &
       *                   spread (d (k+1:n,k  ), 2, n-k+1)
    end do

    ! back substitution phase
    do k = n, 1, -1
       x (k) = d (k, n+1) - sum (d(k,k+1:n)*x(k+1:n))
    end do

    code = .true.
  end function gauss_solve

  subroutine gauss_close ()
    deallocate (d, drow)
  end subroutine gauss_close

end module gauss
