mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-08-29 15:23:43 +02:00
1216 lines
32 KiB
Fortran
1216 lines
32 KiB
Fortran
subroutine non_hrmt_diag_split_degen(n, A, leigvec, reigvec, n_real_eigv, eigval)
|
|
|
|
BEGIN_DOC
|
|
!
|
|
! routine which returns the sorted REAL EIGENVALUES ONLY and corresponding LEFT/RIGHT eigenvetors
|
|
!
|
|
! of a non hermitian matrix A(n,n)
|
|
!
|
|
! n_real_eigv is the number of real eigenvalues, which might be smaller than the dimension "n"
|
|
!
|
|
END_DOC
|
|
|
|
implicit none
|
|
|
|
integer, intent(in) :: n
|
|
double precision, intent(in) :: A(n,n)
|
|
integer, intent(out) :: n_real_eigv
|
|
double precision, intent(out) :: reigvec(n,n), leigvec(n,n), eigval(n)
|
|
double precision, allocatable :: reigvec_tmp(:,:), leigvec_tmp(:,:)
|
|
|
|
integer :: i, j, n_degen,k , iteration
|
|
integer :: n_good
|
|
double precision :: shift,shift_current
|
|
double precision :: r,thr
|
|
integer, allocatable :: list_good(:), iorder_origin(:),iorder(:)
|
|
double precision, allocatable :: WR(:), WI(:), Vl(:,:), VR(:,:),S(:,:)
|
|
double precision, allocatable :: Aw(:,:),diag_elem(:),A_save(:,:)
|
|
double precision, allocatable :: im_part(:),re_part(:)
|
|
|
|
|
|
print*,'Computing the left/right eigenvectors ...'
|
|
print*,'Using the degeneracy splitting algorithm'
|
|
|
|
|
|
! pre-processing the matrix :: sorting by diagonal elements
|
|
allocate(reigvec_tmp(n,n), leigvec_tmp(n,n))
|
|
allocate(diag_elem(n),iorder_origin(n),A_save(n,n))
|
|
do i = 1, n
|
|
iorder_origin(i) = i
|
|
diag_elem(i) = A(i,i)
|
|
enddo
|
|
call dsort(diag_elem, iorder_origin, n)
|
|
do i = 1, n
|
|
do j = 1, n
|
|
A_save(j,i) = A(iorder_origin(j),iorder_origin(i))
|
|
enddo
|
|
enddo
|
|
|
|
shift = 1.d-15
|
|
shift_current = shift
|
|
iteration = 1
|
|
logical :: good_ortho
|
|
good_ortho = .False.
|
|
do while(n_real_eigv.ne.n.or. .not.good_ortho)
|
|
if(shift.gt.1.d-3)then
|
|
print*,'shift > 1.d-3 !!'
|
|
print*,'Your matrix intrinsically contains complex eigenvalues'
|
|
stop
|
|
endif
|
|
print*,'***** iteration = ',iteration
|
|
print*,'shift = ',shift
|
|
allocate(WR(n), WI(n), VL(n,n), VR(n,n), Aw(n,n))
|
|
Aw = A_save
|
|
do i = 1, n
|
|
do j = 1, n
|
|
if(dabs(Aw(j,i)).lt.shift)then
|
|
Aw(j,i) = 0.d0
|
|
endif
|
|
enddo
|
|
enddo
|
|
call lapack_diag_non_sym(n,Aw,WR,WI,VL,VR)
|
|
allocate(im_part(n),iorder(n))
|
|
do i = 1, n
|
|
im_part(i) = -dabs(WI(i))
|
|
iorder(i) = i
|
|
enddo
|
|
call dsort(im_part, iorder, n)
|
|
|
|
shift_current = max(10.d0 * dabs(im_part(1)),shift)
|
|
print*,'Largest imaginary part found in eigenvalues = ',im_part(1)
|
|
print*,'Splitting the degeneracies by ',shift_current
|
|
Aw = A_save
|
|
call split_matrix_degen(Aw,n,shift_current)
|
|
deallocate( im_part, iorder )
|
|
call lapack_diag_non_sym(n,Aw,WR,WI,VL,VR)
|
|
! You track the real eigenvalues
|
|
n_good = 0
|
|
do i = 1, n
|
|
if(dabs(WI(i)).lt.1.d-20)then
|
|
n_good += 1
|
|
else
|
|
print*,'Found an imaginary component to eigenvalue'
|
|
print*,'Re(i) + Im(i)',WR(i),WI(i)
|
|
endif
|
|
enddo
|
|
allocate( list_good(n_good), iorder(n_good) )
|
|
n_good = 0
|
|
do i = 1, n
|
|
if(dabs(WI(i)).lt.1.d-20)then
|
|
n_good += 1
|
|
list_good(n_good) = i
|
|
eigval(n_good) = WR(i)
|
|
endif
|
|
enddo
|
|
deallocate( WR, WI )
|
|
|
|
n_real_eigv = n_good
|
|
do i = 1, n_good
|
|
iorder(i) = i
|
|
enddo
|
|
|
|
! You sort the real eigenvalues
|
|
call dsort(eigval, iorder, n_good)
|
|
|
|
reigvec(:,:) = 0.d0
|
|
leigvec(:,:) = 0.d0
|
|
do i = 1, n_real_eigv
|
|
do j = 1, n
|
|
reigvec_tmp(j,i) = VR(j,list_good(iorder(i)))
|
|
leigvec_tmp(j,i) = Vl(j,list_good(iorder(i)))
|
|
enddo
|
|
enddo
|
|
|
|
if(n_real_eigv == n)then
|
|
allocate(S(n,n))
|
|
call check_bi_ortho(reigvec_tmp,leigvec_tmp,n,S,accu_nd)
|
|
print*,'accu_nd = ',accu_nd
|
|
double precision :: accu_nd
|
|
good_ortho = accu_nd .lt. 1.d-10
|
|
deallocate(S)
|
|
endif
|
|
|
|
deallocate( list_good, iorder )
|
|
deallocate( VL, VR, Aw)
|
|
shift *= 10.d0
|
|
iteration += 1
|
|
enddo
|
|
do i = 1, n
|
|
do j = 1, n
|
|
reigvec(iorder_origin(j),i) = reigvec_tmp(j,i)
|
|
leigvec(iorder_origin(j),i) = leigvec_tmp(j,i)
|
|
enddo
|
|
enddo
|
|
|
|
end subroutine non_hrmt_diag_split_degen
|
|
|
|
! ---
|
|
|
|
subroutine non_hrmt_real_diag_new(n, A, leigvec, reigvec, n_real_eigv, eigval)
|
|
|
|
BEGIN_DOC
|
|
!
|
|
! routine which returns the sorted REAL EIGENVALUES ONLY and corresponding LEFT/RIGHT eigenvetors
|
|
!
|
|
! of a non hermitian matrix A(n,n)
|
|
!
|
|
! n_real_eigv is the number of real eigenvalues, which might be smaller than the dimension "n"
|
|
!
|
|
END_DOC
|
|
|
|
implicit none
|
|
|
|
integer, intent(in) :: n
|
|
double precision, intent(in) :: A(n,n)
|
|
integer, intent(out) :: n_real_eigv
|
|
double precision, intent(out) :: reigvec(n,n), leigvec(n,n), eigval(n)
|
|
|
|
integer :: i, j
|
|
integer :: n_good
|
|
double precision :: shift,shift_current
|
|
double precision :: r,thr
|
|
integer, allocatable :: list_good(:), iorder(:)
|
|
double precision, allocatable :: WR(:), WI(:), Vl(:,:), VR(:,:)
|
|
double precision, allocatable :: Aw(:,:)
|
|
double precision, allocatable :: im_part(:)
|
|
|
|
|
|
print*,'Computing the left/right eigenvectors ...'
|
|
|
|
! Eigvalue(n) = WR(n) + i * WI(n)
|
|
shift = 1.d-10
|
|
do while(n_real_eigv.ne.n.or.shift.gt.1.d-3)
|
|
allocate(WR(n), WI(n), VL(n,n), VR(n,n), Aw(n,n))
|
|
Aw = A
|
|
call lapack_diag_non_sym(n,Aw,WR,WI,VL,VR)
|
|
allocate(im_part(n), iorder(n))
|
|
do i = 1, n
|
|
im_part(i) = -dabs(WI(i))
|
|
iorder(i) = i
|
|
enddo
|
|
shift_current = max(10.d0 * dabs(im_part(1)),shift)
|
|
print*,'adding random number of magnitude ',shift_current
|
|
Aw = A
|
|
do i = 1, n
|
|
call RANDOM_NUMBER(r)
|
|
Aw(i,i) += shift_current * r
|
|
enddo
|
|
deallocate( im_part, iorder )
|
|
call lapack_diag_non_sym(n,Aw,WR,WI,VL,VR)
|
|
|
|
! You track the real eigenvalues
|
|
thr = 1.d-10
|
|
n_good = 0
|
|
do i = 1, n
|
|
if(dabs(WI(i)).lt.thr)then
|
|
n_good += 1
|
|
else
|
|
print*,'Found an imaginary component to eigenvalue'
|
|
print*,'Re(i) + Im(i)',WR(i),WI(i)
|
|
endif
|
|
enddo
|
|
|
|
allocate( list_good(n_good), iorder(n_good) )
|
|
n_good = 0
|
|
do i = 1, n
|
|
if(dabs(WI(i)).lt.thr)then
|
|
n_good += 1
|
|
list_good(n_good) = i
|
|
eigval(n_good) = WR(i)
|
|
endif
|
|
enddo
|
|
|
|
deallocate( WR, WI )
|
|
|
|
n_real_eigv = n_good
|
|
do i = 1, n_good
|
|
iorder(i) = i
|
|
enddo
|
|
|
|
! You sort the real eigenvalues
|
|
call dsort(eigval, iorder, n_good)
|
|
|
|
reigvec(:,:) = 0.d0
|
|
leigvec(:,:) = 0.d0
|
|
do i = 1, n_real_eigv
|
|
do j = 1, n
|
|
reigvec(j,i) = VR(j,list_good(iorder(i)))
|
|
leigvec(j,i) = Vl(j,list_good(iorder(i)))
|
|
enddo
|
|
enddo
|
|
|
|
deallocate( list_good, iorder )
|
|
deallocate( VL, VR, Aw)
|
|
shift *= 10.d0
|
|
enddo
|
|
if(shift.gt.1.d-3)then
|
|
print*,'shift > 1.d-3 !!'
|
|
print*,'Your matrix intrinsically contains complex eigenvalues'
|
|
endif
|
|
|
|
end subroutine non_hrmt_real_diag_new
|
|
|
|
! ---
|
|
|
|
subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, eigval)
|
|
|
|
BEGIN_DOC
|
|
!
|
|
! routine which returns the sorted REAL EIGENVALUES ONLY and corresponding LEFT/RIGHT eigenvetors
|
|
! of a non hermitian matrix A(n,n)
|
|
!
|
|
! n_real_eigv is the number of real eigenvalues, which might be smaller than the dimension "n"
|
|
!
|
|
END_DOC
|
|
|
|
implicit none
|
|
integer, intent(in) :: n
|
|
double precision, intent(in) :: A(n,n)
|
|
double precision, intent(in) :: thr_d, thr_nd
|
|
integer, intent(out) :: n_real_eigv
|
|
double precision, intent(out) :: reigvec(n,n), leigvec(n,n), eigval(n)
|
|
|
|
integer :: i, j,k
|
|
integer :: n_good
|
|
double precision :: thr, thr_cut, thr_diag, thr_norm
|
|
double precision :: accu_d, accu_nd
|
|
|
|
integer, allocatable :: list_good(:), iorder(:), deg_num(:)
|
|
double precision, allocatable :: WR(:), WI(:), VL(:,:), VR(:,:)
|
|
double precision, allocatable :: S(:,:)
|
|
double precision, allocatable :: phi_1_tilde(:),phi_2_tilde(:),chi_1_tilde(:),chi_2_tilde(:)
|
|
|
|
allocate(phi_1_tilde(n),phi_2_tilde(n),chi_1_tilde(n),chi_2_tilde(n))
|
|
|
|
|
|
! -------------------------------------------------------------------------------------
|
|
!
|
|
|
|
!print *, ' '
|
|
!print *, ' Computing the left/right eigenvectors ...'
|
|
!print *, ' '
|
|
|
|
allocate(WR(n), WI(n), VL(n,n), VR(n,n))
|
|
|
|
!print *, ' fock matrix'
|
|
!do i = 1, n
|
|
! write(*, '(1000(F16.10,X))') A(i,:)
|
|
!enddo
|
|
|
|
!thr_cut = 1.d-15
|
|
!call cancel_small_elmts(A, n, thr_cut)
|
|
|
|
!call lapack_diag_non_sym_right(n, A, WR, WI, VR)
|
|
call lapack_diag_non_sym(n, A, WR, WI, VL, VR)
|
|
!call lapack_diag_non_sym_new(n, A, WR, WI, VL, VR)
|
|
|
|
|
|
|
|
!print *, ' '
|
|
!print *, ' eigenvalues'
|
|
i = 1
|
|
do while(i .le. n)
|
|
!write(*, '(I3,X,1000(F16.10,X))')i, WR(i), WI(i)
|
|
if(.false.)then
|
|
if(WI(i).ne.0.d0)then
|
|
print*,'*****************'
|
|
print*,'WARNING ! IMAGINARY EIGENVALUES !!!'
|
|
write(*, '(1000(F16.10,X))') WR(i), WI(i+1)
|
|
! phi = VR(:,i), psi = VR(:,i+1), |Phi_i> = phi + j psi , |Phi_i+1> = phi - j psi
|
|
! chi = VL(:,i), xhi = VL(:,i+1), |Chi_i> = chi + j xhi , |Chi_i+1> = chi - j xhi
|
|
!
|
|
accu_chi_phi = 0.d0
|
|
accu_xhi_psi = 0.d0
|
|
accu_chi_psi = 0.d0
|
|
accu_xhi_phi = 0.d0
|
|
double precision :: accu_chi_phi, accu_xhi_psi, accu_chi_psi, accu_xhi_phi
|
|
double precision :: mat_ovlp(2,2),eigval_tmp(2),eigvec(2,2),mat_ovlp_orig(2,2)
|
|
do j = 1, n
|
|
accu_chi_phi += VL(j,i) * VR(j,i)
|
|
accu_xhi_psi += VL(j,i+1) * VR(j,i+1)
|
|
accu_chi_psi += VL(j,i) * VR(j,i+1)
|
|
accu_xhi_phi += VL(j,i+1) * VR(j,i)
|
|
enddo
|
|
mat_ovlp_orig(1,1) = accu_chi_phi
|
|
mat_ovlp_orig(2,1) = accu_xhi_phi
|
|
mat_ovlp_orig(1,2) = accu_chi_psi
|
|
mat_ovlp_orig(2,2) = accu_xhi_psi
|
|
print*,'old overlap matrix '
|
|
write(*,'(100(F16.10,X))')mat_ovlp_orig(1:2,1)
|
|
write(*,'(100(F16.10,X))')mat_ovlp_orig(1:2,2)
|
|
|
|
|
|
mat_ovlp(1,1) = accu_xhi_phi
|
|
mat_ovlp(2,1) = accu_chi_phi
|
|
mat_ovlp(1,2) = accu_xhi_psi
|
|
mat_ovlp(2,2) = accu_chi_psi
|
|
!print*,'accu_chi_phi = ',accu_chi_phi
|
|
!print*,'accu_xhi_psi = ',accu_xhi_psi
|
|
!print*,'accu_chi_psi = ',accu_chi_psi
|
|
!print*,'accu_xhi_phi = ',accu_xhi_phi
|
|
print*,'new overlap matrix '
|
|
write(*,'(100(F16.10,X))')mat_ovlp(1:2,1)
|
|
write(*,'(100(F16.10,X))')mat_ovlp(1:2,2)
|
|
call lapack_diag(eigval_tmp,eigvec,mat_ovlp,2,2)
|
|
print*,'eigval_tmp(1) = ',eigval_tmp(1)
|
|
print*,'eigvec(1) = ',eigvec(1:2,1)
|
|
print*,'eigval_tmp(2) = ',eigval_tmp(2)
|
|
print*,'eigvec(2) = ',eigvec(1:2,2)
|
|
print*,'*****************'
|
|
phi_1_tilde = 0.d0
|
|
phi_2_tilde = 0.d0
|
|
chi_1_tilde = 0.d0
|
|
chi_2_tilde = 0.d0
|
|
do j = 1, n
|
|
phi_1_tilde(j) += VR(j,i) * eigvec(1,1) + VR(j,i+1) * eigvec(2,1)
|
|
phi_2_tilde(j) += VR(j,i) * eigvec(1,2) + VR(j,i+1) * eigvec(2,2)
|
|
chi_1_tilde(j) += VL(j,i+1) * eigvec(1,1) + VL(j,i) * eigvec(2,1)
|
|
chi_2_tilde(j) += VL(j,i+1) * eigvec(1,2) + VL(j,i) * eigvec(2,2)
|
|
enddo
|
|
VR(1:n,i) = phi_1_tilde(1:n)
|
|
VR(1:n,i+1) = phi_2_tilde(1:n)
|
|
! Vl(1:n,i) = -chi_1_tilde(1:n)
|
|
! Vl(1:n,i+1) = chi_2_tilde(1:n)
|
|
i+=1
|
|
endif
|
|
endif
|
|
i+=1
|
|
enddo
|
|
!print *, ' right eigenvect bef'
|
|
!do i = 1, n
|
|
! write(*, '(1000(F16.10,X))') VR(:,i)
|
|
!enddo
|
|
!print *, ' left eigenvect bef'
|
|
!do i = 1, n
|
|
! write(*, '(1000(F16.10,X))') VL(:,i)
|
|
!enddo
|
|
|
|
thr_diag = 1d-06
|
|
thr_norm = 1d+10
|
|
!call check_EIGVEC(n, n, A, WR, VL, VR, thr_diag, thr_norm, .false.)
|
|
|
|
!
|
|
! -------------------------------------------------------------------------------------
|
|
|
|
! ---
|
|
|
|
! -------------------------------------------------------------------------------------
|
|
! track & sort the real eigenvalues
|
|
|
|
n_good = 0
|
|
!thr = 100d0
|
|
thr = Im_thresh_tcscf
|
|
do i = 1, n
|
|
!print*, 'Re(i) + Im(i)', WR(i), WI(i)
|
|
if(dabs(WI(i)) .lt. thr) then
|
|
n_good += 1
|
|
else
|
|
print*, 'Found an imaginary component to eigenvalue on i = ', i
|
|
print*, 'Re(i) + Im(i)', WR(i), WI(i)
|
|
endif
|
|
enddo
|
|
|
|
if(n_good.ne.n)then
|
|
print*,'there are some imaginary eigenvalues '
|
|
thr_diag = 1d-03
|
|
n_good = n
|
|
endif
|
|
allocate(list_good(n_good), iorder(n_good))
|
|
|
|
n_good = 0
|
|
do i = 1, n
|
|
n_good += 1
|
|
list_good(n_good) = i
|
|
eigval(n_good) = WR(i)
|
|
enddo
|
|
|
|
deallocate( WR, WI )
|
|
|
|
n_real_eigv = n_good
|
|
do i = 1, n_good
|
|
iorder(i) = i
|
|
enddo
|
|
call dsort(eigval, iorder, n_good)
|
|
|
|
reigvec(:,:) = 0.d0
|
|
leigvec(:,:) = 0.d0
|
|
do i = 1, n_real_eigv
|
|
do j = 1, n
|
|
reigvec(j,i) = VR(j,list_good(iorder(i)))
|
|
leigvec(j,i) = VL(j,list_good(iorder(i)))
|
|
enddo
|
|
enddo
|
|
|
|
deallocate( list_good, iorder )
|
|
deallocate( VL, VR )
|
|
|
|
ASSERT(n==n_real_eigv)
|
|
|
|
!print *, ' eigenvalues'
|
|
!do i = 1, n
|
|
! write(*, '(1000(F16.10,X))') eigval(i)
|
|
!enddo
|
|
!print *, ' right eigenvect aft ord'
|
|
!do i = 1, n
|
|
! write(*, '(1000(F16.10,X))') reigvec(:,i)
|
|
!enddo
|
|
!print *, ' left eigenvect aft ord'
|
|
!do i = 1, n
|
|
! write(*, '(1000(F16.10,X))') leigvec(:,i)
|
|
!enddo
|
|
|
|
!
|
|
! -------------------------------------------------------------------------------------
|
|
|
|
! ---
|
|
|
|
! -------------------------------------------------------------------------------------
|
|
! check bi-orthogonality
|
|
|
|
thr_diag = 10.d0
|
|
thr_norm = 1d+10
|
|
|
|
allocate( S(n_real_eigv,n_real_eigv) )
|
|
call check_biorthog(n, n_real_eigv, leigvec, reigvec, accu_d, accu_nd, S, thr_d, thr_nd, .false.)
|
|
|
|
if( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n_real_eigv))/dble(n_real_eigv) .lt. thr_d) ) then
|
|
|
|
print *, ' lapack vectors are normalized and bi-orthogonalized'
|
|
deallocate(S)
|
|
return
|
|
|
|
! accu_nd is modified after adding the normalization
|
|
elseif( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n_real_eigv))/dble(n_real_eigv) .gt. thr_d) ) then
|
|
|
|
print *, ' lapack vectors are not normalized but bi-orthogonalized'
|
|
call check_biorthog_binormalize(n, n_real_eigv, leigvec, reigvec, thr_d, thr_nd, .true.)
|
|
|
|
call check_biorthog(n, n_real_eigv, leigvec, reigvec, accu_d, accu_nd, S, thr_d, thr_nd, .true.)
|
|
call check_EIGVEC(n, n, A, eigval, leigvec, reigvec, thr_diag, thr_norm, .true.)
|
|
|
|
deallocate(S)
|
|
return
|
|
|
|
else
|
|
|
|
print *, ' lapack vectors are not normalized neither bi-orthogonalized'
|
|
|
|
! ---
|
|
|
|
allocate(deg_num(n))
|
|
call reorder_degen_eigvec(n, deg_num, eigval, leigvec, reigvec)
|
|
call impose_biorthog_degen_eigvec(n, deg_num, eigval, leigvec, reigvec)
|
|
deallocate(deg_num)
|
|
|
|
call check_biorthog(n, n_real_eigv, leigvec, reigvec, accu_d, accu_nd, S, thr_d, thr_nd, .false.)
|
|
if( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n_real_eigv)) .gt. thr_d) ) then
|
|
call check_biorthog_binormalize(n, n_real_eigv, leigvec, reigvec, thr_d, thr_nd, .true.)
|
|
endif
|
|
call check_biorthog(n, n_real_eigv, leigvec, reigvec, accu_d, accu_nd, S, thr_d, thr_nd, .true.)
|
|
|
|
!call check_EIGVEC(n, n, A, eigval, leigvec, reigvec, thr_diag, thr_norm, .true.)
|
|
|
|
deallocate(S)
|
|
|
|
endif
|
|
|
|
!
|
|
! -------------------------------------------------------------------------------------
|
|
|
|
return
|
|
|
|
end subroutine non_hrmt_bieig
|
|
|
|
! ---
|
|
|
|
subroutine non_hrmt_bieig_random_diag(n, A, leigvec, reigvec, n_real_eigv, eigval)
|
|
|
|
BEGIN_DOC
|
|
!
|
|
! routine which returns the sorted REAL EIGENVALUES ONLY and corresponding LEFT/RIGHT eigenvetors
|
|
! of a non hermitian matrix A(n,n)
|
|
!
|
|
! n_real_eigv is the number of real eigenvalues, which might be smaller than the dimension "n"
|
|
!
|
|
END_DOC
|
|
|
|
implicit none
|
|
integer, intent(in) :: n
|
|
double precision, intent(in) :: A(n,n)
|
|
integer, intent(out) :: n_real_eigv
|
|
double precision, intent(out) :: reigvec(n,n), leigvec(n,n), eigval(n)
|
|
|
|
integer :: i, j
|
|
integer :: n_good
|
|
double precision :: thr
|
|
double precision :: accu_nd
|
|
|
|
integer, allocatable :: list_good(:), iorder(:)
|
|
double precision, allocatable :: Aw(:,:)
|
|
double precision, allocatable :: WR(:), WI(:), VL(:,:), VR(:,:)
|
|
double precision, allocatable :: S(:,:)
|
|
double precision :: r
|
|
|
|
|
|
! -------------------------------------------------------------------------------------
|
|
!
|
|
|
|
print *, 'Computing the left/right eigenvectors ...'
|
|
allocate( WR(n), WI(n), VL(n,n), VR(n,n), Aw(n,n) )
|
|
|
|
Aw(:,:) = A(:,:)
|
|
call lapack_diag_non_sym_new(n, Aw, WR, WI, VL, VR)
|
|
|
|
thr = 1.d-12
|
|
double precision, allocatable :: im_part(:)
|
|
n_good = 0
|
|
do i = 1, n
|
|
if( dabs(WI(i)).lt.thr ) then
|
|
n_good += 1
|
|
else
|
|
print*, 'Found an imaginary component to eigenvalue on i = ', i
|
|
print*, 'Re(i) + Im(i)', WR(i), WI(i)
|
|
endif
|
|
enddo
|
|
print*,'n_good = ',n_good
|
|
if(n_good .lt. n)then
|
|
print*,'Removing degeneracies to remove imaginary parts'
|
|
allocate(im_part(n),iorder(n))
|
|
r = 0.d0
|
|
do i = 1, n
|
|
im_part(i) = -dabs(WI(i))
|
|
iorder(i) = i
|
|
enddo
|
|
call dsort(im_part,iorder,n)
|
|
thr = 10.d0 * dabs(im_part(1))
|
|
print*,'adding random numbers on the diagonal of magnitude ',thr
|
|
Aw(:,:) = A(:,:)
|
|
do i = 1, n
|
|
call RANDOM_NUMBER(r)
|
|
print*,'r = ',r*thr
|
|
Aw(i,i) += thr * r
|
|
enddo
|
|
print*,'Rediagonalizing the matrix with random numbers'
|
|
call lapack_diag_non_sym_new(n, Aw, WR, WI, VL, VR)
|
|
deallocate(im_part,iorder)
|
|
endif
|
|
deallocate( Aw )
|
|
|
|
!
|
|
! -------------------------------------------------------------------------------------
|
|
|
|
! ---
|
|
|
|
! -------------------------------------------------------------------------------------
|
|
! track & sort the real eigenvalues
|
|
|
|
n_good = 0
|
|
thr = 1.d-5
|
|
do i = 1, n
|
|
if( dabs(WI(i)).lt.thr ) then
|
|
n_good += 1
|
|
else
|
|
print*, 'Found an imaginary component to eigenvalue on i = ', i
|
|
print*, 'Re(i) + Im(i)', WR(i), WI(i)
|
|
endif
|
|
enddo
|
|
print*,'n_good = ',n_good
|
|
allocate( list_good(n_good), iorder(n_good) )
|
|
|
|
n_good = 0
|
|
do i = 1, n
|
|
if( dabs(WI(i)).lt.thr ) then
|
|
n_good += 1
|
|
list_good(n_good) = i
|
|
eigval(n_good) = WR(i)
|
|
endif
|
|
enddo
|
|
|
|
deallocate( WR, WI )
|
|
|
|
n_real_eigv = n_good
|
|
do i = 1, n_good
|
|
iorder(i) = i
|
|
enddo
|
|
call dsort(eigval, iorder, n_good)
|
|
|
|
reigvec(:,:) = 0.d0
|
|
leigvec(:,:) = 0.d0
|
|
do i = 1, n_real_eigv
|
|
do j = 1, n
|
|
reigvec(j,i) = VR(j,list_good(iorder(i)))
|
|
leigvec(j,i) = VL(j,list_good(iorder(i)))
|
|
enddo
|
|
enddo
|
|
|
|
deallocate( list_good, iorder )
|
|
deallocate( VL, VR )
|
|
|
|
!
|
|
! -------------------------------------------------------------------------------------
|
|
|
|
! ---
|
|
|
|
! -------------------------------------------------------------------------------------
|
|
! check bi-orthogonality
|
|
|
|
allocate( S(n_real_eigv,n_real_eigv) )
|
|
|
|
! S = VL x VR
|
|
call dgemm( 'T', 'N', n_real_eigv, n_real_eigv, n, 1.d0 &
|
|
, leigvec, size(leigvec, 1), reigvec, size(reigvec, 1) &
|
|
, 0.d0, S, size(S, 1) )
|
|
|
|
accu_nd = 0.d0
|
|
do i = 1, n_real_eigv
|
|
do j = 1, n_real_eigv
|
|
if(i==j) cycle
|
|
accu_nd = accu_nd + S(j,i) * S(j,i)
|
|
enddo
|
|
enddo
|
|
accu_nd = dsqrt(accu_nd)
|
|
|
|
if(accu_nd .lt. thresh_biorthog_nondiag) then
|
|
! L x R is already bi-orthogonal
|
|
|
|
print *, ' L & T bi-orthogonality: ok'
|
|
deallocate( S )
|
|
return
|
|
|
|
else
|
|
! impose bi-orthogonality
|
|
|
|
print *, ' L & T bi-orthogonality: not imposed yet'
|
|
print *, ' accu_nd = ', accu_nd
|
|
call impose_biorthog_qr(n, n_real_eigv, thresh_biorthog_diag, thresh_biorthog_nondiag, leigvec, reigvec)
|
|
deallocate( S )
|
|
|
|
endif
|
|
|
|
!
|
|
! -------------------------------------------------------------------------------------
|
|
|
|
return
|
|
|
|
end subroutine non_hrmt_bieig_random_diag
|
|
|
|
! ---
|
|
|
|
subroutine non_hrmt_real_im(n, A, leigvec, reigvec, n_real_eigv, eigval)
|
|
|
|
BEGIN_DOC
|
|
!
|
|
! routine which returns the EIGENVALUES sorted the REAL part and corresponding LEFT/RIGHT eigenvetors
|
|
! of a non hermitian matrix A(n,n)
|
|
!
|
|
! n_real_eigv is the number of real eigenvalues, which might be smaller than the dimension "n"
|
|
!
|
|
END_DOC
|
|
|
|
implicit none
|
|
integer, intent(in) :: n
|
|
double precision, intent(in) :: A(n,n)
|
|
integer, intent(out) :: n_real_eigv
|
|
double precision, intent(out) :: reigvec(n,n), leigvec(n,n), eigval(n)
|
|
|
|
integer :: i, j
|
|
integer :: n_bad
|
|
double precision :: thr
|
|
double precision :: accu_nd
|
|
|
|
integer, allocatable :: iorder(:)
|
|
double precision, allocatable :: Aw(:,:)
|
|
double precision, allocatable :: WR(:), WI(:), VL(:,:), VR(:,:)
|
|
double precision, allocatable :: S(:,:)
|
|
double precision :: r
|
|
|
|
! -------------------------------------------------------------------------------------
|
|
!
|
|
|
|
print *, 'Computing the left/right eigenvectors ...'
|
|
allocate( WR(n), WI(n), VL(n,n), VR(n,n), Aw(n,n), iorder(n))
|
|
|
|
Aw(:,:) = A(:,:)
|
|
do i = 1, n
|
|
call RANDOM_NUMBER(r)
|
|
Aw(i,i) += 10.d-10* r
|
|
enddo
|
|
call lapack_diag_non_sym(n, Aw, WR, WI, VL, VR)
|
|
|
|
! -------------------------------------------------------------------------------------
|
|
! track & sort the real eigenvalues
|
|
|
|
i = 1
|
|
thr = 1.d-15
|
|
n_real_eigv = 0
|
|
do while (i.le.n)
|
|
! print*,i,dabs(WI(i))
|
|
if( dabs(WI(i)).gt.thr ) then
|
|
print*, 'Found an imaginary component to eigenvalue on i = ', i
|
|
print*, 'Re(i) , Im(i) ', WR(i), WI(i)
|
|
iorder(i) = i
|
|
eigval(i) = WR(i)
|
|
i+=1
|
|
print*, 'Re(i+1),Im(i+1)',WR(i), WI(i)
|
|
iorder(i) = i
|
|
eigval(i) = WR(i)
|
|
i+=1
|
|
else
|
|
n_real_eigv += 1
|
|
iorder(i) = i
|
|
eigval(i) = WR(i)
|
|
i+=1
|
|
endif
|
|
enddo
|
|
call dsort(eigval, iorder, n)
|
|
reigvec(:,:) = 0.d0
|
|
leigvec(:,:) = 0.d0
|
|
do i = 1, n
|
|
do j = 1, n
|
|
reigvec(j,i) = VR(j,iorder(i))
|
|
leigvec(j,i) = VL(j,iorder(i))
|
|
enddo
|
|
enddo
|
|
|
|
deallocate( iorder )
|
|
deallocate( VL, VR )
|
|
|
|
!
|
|
! -------------------------------------------------------------------------------------
|
|
|
|
! ---
|
|
|
|
! -------------------------------------------------------------------------------------
|
|
! check bi-orthogonality
|
|
|
|
allocate( S(n,n) )
|
|
|
|
! S = VL x VR
|
|
call dgemm( 'T', 'N', n, n, n, 1.d0 &
|
|
, leigvec, size(leigvec, 1), reigvec, size(reigvec, 1) &
|
|
, 0.d0, S, size(S, 1) )
|
|
|
|
accu_nd = 0.d0
|
|
do i = 1, n
|
|
do j = 1, n
|
|
if(i==j) cycle
|
|
accu_nd = accu_nd + S(j,i) * S(j,i)
|
|
enddo
|
|
enddo
|
|
accu_nd = dsqrt(accu_nd)
|
|
|
|
deallocate( S )
|
|
|
|
end subroutine non_hrmt_real_im
|
|
|
|
! ---
|
|
|
|
subroutine non_hrmt_generalized_real_im(n, A, B, leigvec, reigvec, n_real_eigv, eigval)
|
|
|
|
BEGIN_DOC
|
|
!
|
|
! routine which returns the EIGENVALUES sorted the REAL part and corresponding LEFT/RIGHT eigenvetors
|
|
! for A R = lambda B R and A^\dagger L = lambda B^\dagger L
|
|
!
|
|
! n_real_eigv is the number of real eigenvalues, which might be smaller than the dimension "n"
|
|
!
|
|
END_DOC
|
|
|
|
implicit none
|
|
integer, intent(in) :: n
|
|
double precision, intent(in) :: A(n,n),B(n,n)
|
|
integer, intent(out) :: n_real_eigv
|
|
double precision, intent(out) :: reigvec(n,n), leigvec(n,n), eigval(n)
|
|
|
|
integer :: i, j
|
|
integer :: n_bad
|
|
double precision :: thr
|
|
double precision :: accu_nd
|
|
|
|
integer, allocatable :: iorder(:)
|
|
double precision, allocatable :: Aw(:,:),Bw(:,:)
|
|
double precision, allocatable :: WR(:), WI(:), VL(:,:), VR(:,:), beta(:)
|
|
double precision, allocatable :: S(:,:)
|
|
double precision :: r
|
|
|
|
! -------------------------------------------------------------------------------------
|
|
!
|
|
|
|
print *, 'Computing the left/right eigenvectors ...'
|
|
allocate( WR(n), WI(n), VL(n,n), VR(n,n), Aw(n,n), Bw(n,n),iorder(n),beta(n))
|
|
|
|
Aw(:,:) = A(:,:)
|
|
Bw(:,:) = B(:,:)
|
|
call lapack_diag_general_non_sym(n,Aw,Bw,WR,beta,WI,VL,VR)
|
|
|
|
! -------------------------------------------------------------------------------------
|
|
! track & sort the real eigenvalues
|
|
|
|
i = 1
|
|
thr = 1.d-10
|
|
n_real_eigv = 0
|
|
do while (i.le.n)
|
|
if( dabs(WI(i)).gt.thr ) then
|
|
print*, 'Found an imaginary component to eigenvalue on i = ', i
|
|
print*, 'Re(i) , Im(i) ', WR(i), WI(i)
|
|
iorder(i) = i
|
|
eigval(i) = WR(i)/(beta(i) + 1.d-10)
|
|
i+=1
|
|
print*, 'Re(i+1),Im(i+1)',WR(i), WI(i)
|
|
iorder(i) = i
|
|
eigval(i) = WR(i)/(beta(i) + 1.d-10)
|
|
i+=1
|
|
else
|
|
n_real_eigv += 1
|
|
iorder(i) = i
|
|
eigval(i) = WR(i)/(beta(i) + 1.d-10)
|
|
i+=1
|
|
endif
|
|
enddo
|
|
call dsort(eigval, iorder, n)
|
|
reigvec(:,:) = 0.d0
|
|
leigvec(:,:) = 0.d0
|
|
do i = 1, n
|
|
do j = 1, n
|
|
reigvec(j,i) = VR(j,iorder(i))
|
|
leigvec(j,i) = VL(j,iorder(i))
|
|
enddo
|
|
enddo
|
|
|
|
deallocate( iorder )
|
|
deallocate( VL, VR )
|
|
|
|
!
|
|
! -------------------------------------------------------------------------------------
|
|
|
|
! ---
|
|
|
|
! -------------------------------------------------------------------------------------
|
|
! check bi-orthogonality
|
|
|
|
allocate( S(n,n) )
|
|
|
|
! S = VL x VR
|
|
call dgemm( 'T', 'N', n, n, n, 1.d0 &
|
|
, leigvec, size(leigvec, 1), reigvec, size(reigvec, 1) &
|
|
, 0.d0, S, size(S, 1) )
|
|
|
|
accu_nd = 0.d0
|
|
do i = 1, n
|
|
do j = 1, n
|
|
if(i==j) cycle
|
|
accu_nd = accu_nd + S(j,i) * S(j,i)
|
|
enddo
|
|
enddo
|
|
accu_nd = dsqrt(accu_nd)
|
|
|
|
deallocate( S )
|
|
|
|
end subroutine non_hrmt_generalized_real_im
|
|
|
|
! ---
|
|
|
|
subroutine non_hrmt_bieig_fullvect(n, A, leigvec, reigvec, n_real_eigv, eigval)
|
|
|
|
BEGIN_DOC
|
|
!
|
|
! routine which returns the sorted REAL EIGENVALUES ONLY and corresponding LEFT/RIGHT eigenvetors
|
|
! of a non hermitian matrix A(n,n)
|
|
!
|
|
! n_real_eigv is the number of real eigenvalues, which might be smaller than the dimension "n"
|
|
!
|
|
END_DOC
|
|
|
|
implicit none
|
|
integer, intent(in) :: n
|
|
double precision, intent(in) :: A(n,n)
|
|
integer, intent(out) :: n_real_eigv
|
|
double precision, intent(out) :: reigvec(n,n), leigvec(n,n), eigval(n)
|
|
|
|
integer :: i, j
|
|
integer :: n_good
|
|
double precision :: thr
|
|
double precision :: accu_nd
|
|
|
|
integer, allocatable :: iorder(:)
|
|
double precision, allocatable :: Aw(:,:)
|
|
double precision, allocatable :: WR(:), WI(:), VL(:,:), VR(:,:)
|
|
double precision, allocatable :: S(:,:)
|
|
double precision, allocatable :: eigval_sorted(:)
|
|
|
|
|
|
! -------------------------------------------------------------------------------------
|
|
!
|
|
|
|
print *, 'Computing the left/right eigenvectors ...'
|
|
|
|
allocate( WR(n), WI(n), VL(n,n), VR(n,n), Aw(n,n) )
|
|
Aw(:,:) = A(:,:)
|
|
|
|
call lapack_diag_non_sym_new(n, Aw, WR, WI, VL, VR)
|
|
|
|
deallocate( Aw )
|
|
|
|
!
|
|
! -------------------------------------------------------------------------------------
|
|
|
|
! ---
|
|
|
|
! -------------------------------------------------------------------------------------
|
|
! track & sort the real eigenvalues
|
|
|
|
allocate( eigval_sorted(n), iorder(n) )
|
|
|
|
n_good = 0
|
|
thr = 1.d-10
|
|
|
|
do i = 1, n
|
|
|
|
iorder(i) = i
|
|
eigval_sorted(i) = WR(i)
|
|
|
|
if(dabs(WI(i)) .gt. thr) then
|
|
print*, ' Found an imaginary component to eigenvalue on i = ', i
|
|
print*, ' Re(i) + Im(i)', WR(i), WI(i)
|
|
else
|
|
n_good += 1
|
|
endif
|
|
|
|
enddo
|
|
|
|
n_real_eigv = n_good
|
|
|
|
call dsort(eigval_sorted, iorder, n)
|
|
|
|
reigvec(:,:) = 0.d0
|
|
leigvec(:,:) = 0.d0
|
|
do i = 1, n
|
|
eigval(i) = WR(i)
|
|
do j = 1, n
|
|
reigvec(j,i) = VR(j,iorder(i))
|
|
leigvec(j,i) = VL(j,iorder(i))
|
|
enddo
|
|
enddo
|
|
|
|
deallocate( eigval_sorted, iorder )
|
|
deallocate( WR, WI )
|
|
deallocate( VL, VR )
|
|
|
|
!
|
|
! -------------------------------------------------------------------------------------
|
|
|
|
! ---
|
|
|
|
! -------------------------------------------------------------------------------------
|
|
! check bi-orthogonality
|
|
|
|
allocate( S(n,n) )
|
|
|
|
! S = VL x VR
|
|
call dgemm( 'T', 'N', n, n, n, 1.d0 &
|
|
, leigvec, size(leigvec, 1), reigvec, size(reigvec, 1) &
|
|
, 0.d0, S, size(S, 1) )
|
|
|
|
accu_nd = 0.d0
|
|
do i = 1, n
|
|
do j = 1, n
|
|
if(i==j) cycle
|
|
accu_nd = accu_nd + S(j,i) * S(j,i)
|
|
enddo
|
|
enddo
|
|
accu_nd = dsqrt(accu_nd)
|
|
|
|
if(accu_nd .lt. thresh_biorthog_nondiag) then
|
|
! L x R is already bi-orthogonal
|
|
|
|
!print *, ' L & T bi-orthogonality: ok'
|
|
deallocate( S )
|
|
return
|
|
|
|
else
|
|
! impose bi-orthogonality
|
|
|
|
!print *, ' L & T bi-orthogonality: not imposed yet'
|
|
!print *, ' accu_nd = ', accu_nd
|
|
call impose_biorthog_qr(n, n, thresh_biorthog_diag, thresh_biorthog_nondiag, leigvec, reigvec)
|
|
deallocate( S )
|
|
|
|
endif
|
|
|
|
!
|
|
! -------------------------------------------------------------------------------------
|
|
|
|
return
|
|
|
|
end subroutine non_hrmt_bieig_fullvect
|
|
|
|
! ---
|
|
|
|
|
|
subroutine split_matrix_degen(aw,n,shift)
|
|
implicit none
|
|
BEGIN_DOC
|
|
! subroutines that splits the degeneracies of a matrix by adding a splitting of magnitude thr * n_degen/2
|
|
!
|
|
! WARNING !! THE MATRIX IS ASSUMED TO BE PASSED WITH INCREASING DIAGONAL ELEMENTS
|
|
END_DOC
|
|
double precision,intent(inout) :: Aw(n,n)
|
|
double precision,intent(in) :: shift
|
|
integer, intent(in) :: n
|
|
integer :: i,j,n_degen
|
|
logical :: keep_on
|
|
i=1
|
|
do while(i.lt.n)
|
|
if(dabs(Aw(i,i)-Aw(i+1,i+1)).lt.shift)then
|
|
j=1
|
|
keep_on = .True.
|
|
do while(keep_on)
|
|
if(i+j.gt.n)then
|
|
keep_on = .False.
|
|
exit
|
|
endif
|
|
if(dabs(Aw(i,i)-Aw(i+j,i+j)).lt.shift)then
|
|
j+=1
|
|
else
|
|
keep_on=.False.
|
|
exit
|
|
endif
|
|
enddo
|
|
n_degen = j
|
|
j=0
|
|
keep_on = .True.
|
|
do while(keep_on)
|
|
if(i+j+1.gt.n)then
|
|
keep_on = .False.
|
|
exit
|
|
endif
|
|
if(dabs(Aw(i+j,i+j)-Aw(i+j+1,i+j+1)).lt.shift)then
|
|
Aw(i+j,i+j) += (j-n_degen/2) * shift
|
|
j+=1
|
|
else
|
|
keep_on = .False.
|
|
exit
|
|
endif
|
|
enddo
|
|
Aw(i+n_degen-1,i+n_degen-1) += (n_degen-1-n_degen/2) * shift
|
|
i+=n_degen
|
|
else
|
|
i+=1
|
|
endif
|
|
enddo
|
|
|
|
end
|
|
|
|
subroutine give_degen(a,n,shift,list_degen,n_degen_list)
|
|
implicit none
|
|
BEGIN_DOC
|
|
! returns n_degen_list :: the number of degenerated SET of elements (i.e. with |A(i)-A(i+1)| below shift)
|
|
!
|
|
! for each of these sets, list_degen(1,i) = first degenerate element of the set i,
|
|
!
|
|
! list_degen(2,i) = last degenerate element of the set i.
|
|
END_DOC
|
|
double precision,intent(in) :: A(n)
|
|
double precision,intent(in) :: shift
|
|
integer, intent(in) :: n
|
|
integer, intent(out) :: list_degen(2,n),n_degen_list
|
|
integer :: i,j,n_degen,k
|
|
logical :: keep_on
|
|
double precision,allocatable :: Aw(:)
|
|
list_degen = -1
|
|
allocate(Aw(n))
|
|
Aw = A
|
|
i=1
|
|
k = 0
|
|
do while(i.lt.n)
|
|
if(dabs(Aw(i)-Aw(i+1)).lt.shift)then
|
|
k+=1
|
|
j=1
|
|
list_degen(1,k) = i
|
|
keep_on = .True.
|
|
do while(keep_on)
|
|
if(i+j.gt.n)then
|
|
keep_on = .False.
|
|
exit
|
|
endif
|
|
if(dabs(Aw(i)-Aw(i+j)).lt.shift)then
|
|
j+=1
|
|
else
|
|
keep_on=.False.
|
|
exit
|
|
endif
|
|
enddo
|
|
n_degen = j
|
|
list_degen(2,k) = list_degen(1,k)-1 + n_degen
|
|
j=0
|
|
keep_on = .True.
|
|
do while(keep_on)
|
|
if(i+j+1.gt.n)then
|
|
keep_on = .False.
|
|
exit
|
|
endif
|
|
if(dabs(Aw(i+j)-Aw(i+j+1)).lt.shift)then
|
|
Aw(i+j) += (j-n_degen/2) * shift
|
|
j+=1
|
|
else
|
|
keep_on = .False.
|
|
exit
|
|
endif
|
|
enddo
|
|
Aw(i+n_degen-1) += (n_degen-1-n_degen/2) * shift
|
|
i+=n_degen
|
|
else
|
|
i+=1
|
|
endif
|
|
enddo
|
|
n_degen_list = k
|
|
|
|
end
|
|
|
|
subroutine cancel_small_elmts(aw,n,shift)
|
|
implicit none
|
|
BEGIN_DOC
|
|
! subroutines that splits the degeneracies of a matrix by adding a splitting of magnitude thr * n_degen/2
|
|
!
|
|
! WARNING !! THE MATRIX IS ASSUMED TO BE PASSED WITH INCREASING DIAGONAL ELEMENTS
|
|
END_DOC
|
|
double precision,intent(inout) :: Aw(n,n)
|
|
double precision,intent(in) :: shift
|
|
integer, intent(in) :: n
|
|
integer :: i,j
|
|
do i = 1, n
|
|
do j = 1, n
|
|
if(dabs(Aw(j,i)).lt.shift)then
|
|
Aw(j,i) = 0.d0
|
|
endif
|
|
enddo
|
|
enddo
|
|
end
|
|
|
|
subroutine check_bi_ortho(reigvec,leigvec,n,S,accu_nd)
|
|
implicit none
|
|
integer, intent(in) :: n
|
|
double precision,intent(in) :: reigvec(n,n),leigvec(n,n)
|
|
double precision, intent(out) :: S(n,n),accu_nd
|
|
BEGIN_DOC
|
|
! retunrs the overlap matrix S = Leigvec^T Reigvec
|
|
!
|
|
! and the square root of the sum of the squared off-diagonal elements of S
|
|
END_DOC
|
|
integer :: i,j
|
|
! S = VL x VR
|
|
call dgemm( 'T', 'N', n, n, n, 1.d0 &
|
|
, leigvec, size(leigvec, 1), reigvec, size(reigvec, 1) &
|
|
, 0.d0, S, size(S, 1) )
|
|
accu_nd = 0.d0
|
|
do i = 1, n
|
|
do j = 1, n
|
|
if(i.ne.j) then
|
|
accu_nd = accu_nd + S(j,i) * S(j,i)
|
|
endif
|
|
enddo
|
|
enddo
|
|
accu_nd = dsqrt(accu_nd)
|
|
|
|
end
|