subroutine eig_real(matrix, e_values, e_vectors, t_left_ev)
! for very restricted matrices do a diag routine here!
real(dp), intent(in) :: matrix(:,:)
real(dp), intent(out) :: e_values(size(matrix,1))
real(dp), intent(out), optional :: e_vectors(size(matrix,1),size(matrix,1))
logical, intent(in), optional :: t_left_ev
! get the specifics for the eigenvectors still..
! i think i need a bigger work, and maybe also a flag for how many
! eigenvectors i want.. maybe also the number of eigenvalues..
integer :: n, info
real(dp) :: work(4*size(matrix,1)), tmp_matrix(size(matrix,1),size(matrix,2))
real(dp) :: left_ev(size(matrix,1),size(matrix,1)), dummy_eval(size(matrix,1))
real(dp) :: right_ev(size(matrix,1),size(matrix,1))
integer :: sort_ind(size(matrix,1))
character :: left, right
! and convention is: we only want the right eigenvectors!!
! and always assume real-only eigenvalues
if (present(e_vectors)) then
if (present(t_left_ev)) then
if (t_left_ev) then
left = 'V'
right = 'N'
else
left = 'N'
right = 'V'
end if
else
left = 'N'
right = 'V'
end if
n = size(matrix,1)
tmp_matrix = matrix
left = 'V'
right = 'V'
call dgeev(&
left, &
right, &
n, &
tmp_matrix, &
n, &
e_values, &
dummy_eval, &
left_ev, &
n, &
right_ev, &
n, &
work, &
4*n, &
info)
sort_ind = [(n, n = 1, size(matrix,1))]
call sort(e_values, sort_ind)
if (present(t_left_ev)) then
if (t_left_ev) then
e_vectors = left_ev(:,sort_ind)
else
e_vectors = right_ev(:,sort_ind)
end if
else
e_vectors = right_ev(:,sort_ind)
end if
else
e_values = calc_eigenvalues(matrix)
end if
end subroutine eig_real