pure subroutine dggev_rel(matrix, eig_vecr, eig_val, eig_vecl)
real(rk), dimension(:,:), intent(in) :: matrix
real(rk), dimension(:), allocatable, intent(out) :: eig_val
real(rk), dimension(:,:), allocatable, intent(out) :: eig_vecr
real(rk), dimension(:,:), allocatable, intent(out), optional :: eig_vecl
real(rk), dimension(size(matrix,1)) :: alphar, alphai, beta
real(rk), dimension(size(matrix,1),size(matrix,1)) :: vl, vr
real(rk), dimension(:), allocatable :: work
integer :: m, lwork, info
real(rk) :: work1(1)
real(rk), dimension(size(matrix,1),size(matrix,1)) :: A
interface
pure subroutine dggev(fjobvl, fjobvr, fn, fa, flda, fb, ldb, &
falphar, falphai, fbeta, fvl, fldvl, fvr, fldvr, fwork, flwork, finfo)
import :: rk
character, intent(in) :: fjobvl, fjobvr
integer, intent(in) :: fn, flda, ldb, fldvl, fldvr, flwork, finfo
real(rk), intent(inout) :: fa(flda, *), fb(ldb, *)
real(rk), intent(out) :: falphar(fn), falphai(fn), fbeta(fn)
real(rk), intent(inout) :: fvl(fldvl, *), fvr(fldvr, *)
real(rk), intent(out) :: fwork(flwork)
end subroutine
end interface
m = size(matrix, 1)
allocate(eig_val(m), eig_vecr(m, m))
eig_vecr = matrix
A = matrix
if (present(eig_vecl)) then
eig_vecl = matrix
call dggev('V', 'V', m, eig_vecl, m,&
eig_vecr, m, alphar, alphai, beta, vl, m, vr, m, work1, -1, info)
else
call dggev('N', 'V', m, eig_vecr, m,&
A, m, alphar, alphai, beta, vl, m, vr, m, work1, -1, info)
end if
lwork = nint(work1(1))
allocate(work(lwork))
if (present(eig_vecl)) then
call dggev('V', 'V', m, eig_vecl, m, &
eig_vecr, m, alphar, alphai, beta, vl, m, vr, m, work, lwork, info)
else
call dggev('N', 'V', m, eig_vecr, m,&
A, m, alphar, alphai, beta, vl, m, vr, m, work, lwork, info)
end if
eig_val = alphar / beta
deallocate(work)
end subroutine dggev_rel