determine_eigenspaces Subroutine

    pure subroutine determine_eigenspaces(lambda,dimensions)
        real(dp), intent(inout) :: lambda(:)
        integer, allocatable, intent(out) :: dimensions(:)
        integer, allocatable :: d_buffer(:)
        integer :: i, low
        integer :: n_spaces
        debug_function_name("determine_eigenspaces")

        ASSERT(all(lambda(: size(lambda) - 1) <= lambda(2 : )))

        allocate(d_buffer(size(lambda)), source=0)

        low = 1
        n_spaces = 1
        do i=1, size(lambda)
            d_buffer(n_spaces) = d_buffer(n_spaces)+1
            if (i + 1 <= size(lambda)) then
                if (.not. isclose(lambda(i), lambda(i + 1), &
                                   epsilon(lambda) * 1.0e3_dp, 1.0e-8_dp)) then
                    n_spaces = n_spaces + 1
                    lambda(low : i) = mean(lambda(low : i))
                    low = i+1
                end if
            else
                lambda(low:i) = mean(lambda(low:i))
            end if
        end do

        dimensions = d_buffer(:n_spaces)

        contains

            real(dp) pure function mean(X)
                real(dp), intent(in) :: X(:)
                mean = sum(X) / size(X)
            end function mean
    end subroutine determine_eigenspaces