matrix_exponential Function

    function matrix_exponential(matrix) result(exp_matrix)
        ! calculate the matrix exponential of a real, symmetric 2-D matrix with lapack
        ! routines
        ! i need A = UDU^-1
        ! e^A = Ue^DU^-1
        HElement_t(dp), intent(in) :: matrix(:,:)
        HElement_t(dp) :: exp_matrix(size(matrix,1),size(matrix,2))

        ! maybe i need to allocate this stuff:
        HElement_t(dp) :: vectors(size(matrix,1),size(matrix,2))
        real(dp) :: values(size(matrix,1))
        HElement_t(dp) :: work(3*size(matrix,1)-1)
        HElement_t(dp) :: inverse(size(matrix,1),size(matrix,2))
        HElement_t(dp) :: exp_diag(size(matrix,1),size(matrix,2))
        integer :: info, n

        n = size(matrix,1)

        ! first i need to diagonalise the matrix and calculate the
        ! eigenvectors
        vectors = matrix
#ifdef CMPLX_
        block
            real(dp), allocatable :: rwork(:)
            allocate(rwork(max(1, 3*n - 2)))
            call zheev('V', 'U', n, vectors, n, values, work, 3*n-1, rwork, info)
            deallocate(rwork)
        end block
#else
        call dsyev('V', 'U', n, vectors, n, values, work, 3*n-1,info)
#endif
        ! now i have the eigenvectors, which i need the inverse of
        ! it is rotation only or? so i would just need a transpose or?
        inverse = transpose(vectors)

        ! i need to construct exp(eigenvalues) as a diagonal matrix!
        exp_diag = matrix_diag(exp(values))

        exp_matrix = blas_matmul(blas_matmul(vectors,exp_diag), inverse)

    end function matrix_exponential