# matrix_exponential Function

## public function matrix_exponential(matrix) result(exp_matrix)

### Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in) :: matrix(:,:)

## Source Code

    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