75 lines
1.9 KiB
Forth
75 lines
1.9 KiB
Forth
subroutine lapack_g_non_sym_real(n, H, LDH, S, LDS, beta, &
|
|
n_real, vl, LDVL, vr, LDVR)
|
|
implicit none
|
|
integer, intent(in) :: n, LDH, LDS, LDVL, LDVR
|
|
double precision, intent(in) :: H(LDH,n), S(LDS,n)
|
|
double precision, intent(out) :: VL(LDVL,n), VR(LDVR,n), beta(n)
|
|
integer, intent(out) :: n_real
|
|
|
|
integer :: lwork, info, i,j
|
|
double precision, allocatable :: work(:), Htmp(:,:), Stmp(:,:)
|
|
double precision, allocatable :: alphar(:), alphai(:), vltmp(:,:), vrtmp(:,:)
|
|
integer, allocatable :: iorder(:), list_good(:)
|
|
|
|
|
|
lwork = -1
|
|
allocate(work(1), alphar(n), alphai(n), vltmp(n,n), vrtmp(n,n), &
|
|
Htmp(n,n), Stmp(n,n), list_good(n))
|
|
|
|
Htmp(1:n,1:n) = H(1:n,1:n)
|
|
Stmp(1:n,1:n) = S(1:n,1:n)
|
|
call dggev('V', 'V', n, Htmp, size(Htmp,1), Stmp, size(Stmp,1), alphar, alphai, beta, &
|
|
vltmp, size(vltmp,1), vrtmp, size(vrtmp,1), work, lwork, info)
|
|
|
|
lwork = int(work(1))
|
|
deallocate(work)
|
|
allocate(work(lwork))
|
|
call dggev('V', 'V', n, Htmp, size(Htmp,1), Stmp, size(Stmp,1), alphar, alphai, beta, &
|
|
vltmp, size(vltmp,1), vrtmp, size(vrtmp,1), work, lwork, info)
|
|
|
|
deallocate(work)
|
|
if (info /= 0) then
|
|
stop 'DGGEV Diagonalization failed'
|
|
endif
|
|
|
|
allocate(iorder(n))
|
|
n_real = 0
|
|
do i=1,n
|
|
if (dabs(alphai(i)) < 1.d-10) then
|
|
n_real += 1
|
|
list_good(n_real) = i
|
|
else
|
|
alphar(i) = dble(huge(1.0))
|
|
endif
|
|
enddo
|
|
|
|
do i=1,n
|
|
if (dabs(beta(i)/(alphar(i)+1.d-12)) < 1.d-10) then
|
|
beta(i) = 0.d0
|
|
else
|
|
beta(i) = alphar(i)/beta(i)
|
|
endif
|
|
enddo
|
|
|
|
do i=1,n_real
|
|
iorder(i) = i
|
|
do j=1,n
|
|
vrtmp(j,i) = vrtmp(j,list_good(i))
|
|
vltmp(j,i) = vltmp(j,list_good(i))
|
|
end do
|
|
end do
|
|
|
|
call dsort(beta, iorder, n_real)
|
|
|
|
do i=1,n_real
|
|
do j=1,n
|
|
vr(j,i) = vrtmp(j,iorder(i))
|
|
vl(j,i) = vltmp(j,iorder(i))
|
|
end do
|
|
end do
|
|
|
|
deallocate(vrtmp, vltmp, iorder, Htmp, Stmp)
|
|
|
|
end
|
|
|