LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 03:23:29 +01:00
Code Releases Activity
b25489e14c
qp2/plugins/local/non_hermit_dav/biorthog.irp.f

1227 lines
32 KiB
Fortran
Raw Blame History

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(:)
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_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'
! ---
! call impose_orthog_degen_eigvec(n, eigval, reigvec)
! call impose_orthog_degen_eigvec(n, eigval, leigvec)
call reorder_degen_eigvec(n, eigval, leigvec, reigvec)
call impose_biorthog_degen_eigvec(n, eigval, leigvec, reigvec)
!call impose_orthog_biorthog_degen_eigvec(n, thr_d, thr_nd, eigval, leigvec, reigvec)
!call impose_unique_biorthog_degen_eigvec(n, eigval, mo_coef, ao_overlap, leigvec, reigvec)
! ---
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 impose_biorthog_qr(n, n_real_eigv, thr_d, thr_nd, leigvec, reigvec)
!call impose_biorthog_lu(n, n_real_eigv, thr_d, thr_nd, leigvec, reigvec)
! ---
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
View Git Blame Copy Permalink