matrix_diag Function

    function matrix_diag(vector) result(diag)
        ! constructs a diagonal matrix with the vector on the diagonal
        real(dp), intent(in) :: vector(:)
        HElement_t(dp) :: diag(size(vector),size(vector))

        integer :: i

        diag = 0.0_dp

        do i = 1, size(vector)
            diag(i,i) = vector(i)
        end do

    end function matrix_diag