From ca4cdf56d5e022b06d59f57097d2ab9cf29856f6 Mon Sep 17 00:00:00 2001 From: eginer Date: Mon, 6 Feb 2023 19:03:22 +0100 Subject: [PATCH] added non_hermit_dav --- src/non_hermit_dav/NEED | 1 + src/non_hermit_dav/biorthog.irp.f | 1156 +++++++ src/non_hermit_dav/gram_schmit.irp.f | 56 + src/non_hermit_dav/htilde_mat.irp.f | 93 + .../lapack_diag_non_hermit.irp.f | 2907 +++++++++++++++++ src/non_hermit_dav/new_routines.irp.f | 670 ++++ src/non_hermit_dav/project.irp.f | 53 + src/non_hermit_dav/utils.irp.f | 325 ++ 8 files changed, 5261 insertions(+) create mode 100644 src/non_hermit_dav/NEED create mode 100644 src/non_hermit_dav/biorthog.irp.f create mode 100644 src/non_hermit_dav/gram_schmit.irp.f create mode 100644 src/non_hermit_dav/htilde_mat.irp.f create mode 100644 src/non_hermit_dav/lapack_diag_non_hermit.irp.f create mode 100644 src/non_hermit_dav/new_routines.irp.f create mode 100644 src/non_hermit_dav/project.irp.f create mode 100644 src/non_hermit_dav/utils.irp.f diff --git a/src/non_hermit_dav/NEED b/src/non_hermit_dav/NEED new file mode 100644 index 00000000..9487075c --- /dev/null +++ b/src/non_hermit_dav/NEED @@ -0,0 +1 @@ +utils diff --git a/src/non_hermit_dav/biorthog.irp.f b/src/non_hermit_dav/biorthog.irp.f new file mode 100644 index 00000000..78fddf54 --- /dev/null +++ b/src/non_hermit_dav/biorthog.irp.f @@ -0,0 +1,1156 @@ +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 + 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(:,:) + + + ! ------------------------------------------------------------------------------------- + ! + + !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' + !do i = 1, n + ! write(*, '(1000(F16.10,X))') WR(i), WI(i) + !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 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 diff --git a/src/non_hermit_dav/gram_schmit.irp.f b/src/non_hermit_dav/gram_schmit.irp.f new file mode 100644 index 00000000..520661b8 --- /dev/null +++ b/src/non_hermit_dav/gram_schmit.irp.f @@ -0,0 +1,56 @@ +subroutine bi_ortho_gram_schmidt(wi,vi,n,ni,wk,wk_schmidt) + implicit none + BEGIN_DOC +! you enter with a set of "ni" BI-ORTHONORMAL vectors of length "n" +! +! vi(j,i) = , wi(j,i) = , = delta_{ij} S_ii, S_ii = +! +! and a vector vk(j) = +! +! you go out with a vector vk_schmidt(j) = +! +! which is Gram-Schmidt orthonormalized with respect to the "vi" +! +! = 0 +! +! |wk_schmidt> = |wk> - \sum_{i=1}^ni (/) |wi> +! +! according to Eq. (5), (6) of Computers Structures, Vol 56, No. 4, pp 605-613, 1995 +! +! https://doi.org/10.1016/0045-7949(94)00565-K + END_DOC + integer, intent(in) :: n,ni + double precision, intent(in) :: wi(n,ni),vi(n,ni),wk(n) + double precision, intent(out):: wk_schmidt(n) + double precision :: vi_wk,u_dot_v,tmp,u_dot_u + double precision, allocatable :: sii(:) + integer :: i,j + allocate( sii(ni) ) + wk_schmidt = wk + do i = 1, ni + sii(i) = u_dot_v(vi(1,i),wi(1,i),n) + enddo +! do i = 1, n +! print*,i,'wk',wk(i) +! enddo +! print*,'' +! print*,'' + do i = 1, ni +! print*,'i',i + ! Gram-Schmidt + vi_wk = u_dot_v(vi(1,i),wk,n) + vi_wk = vi_wk / sii(i) +! print*,'' + do j = 1, n +! print*,j,vi_wk,wi(j,i) + wk_schmidt(j) -= vi_wk * wi(j,i) + enddo + enddo + tmp = u_dot_u(wk_schmidt,n) + tmp = 1.d0/dsqrt(tmp) + wk_schmidt = tmp * wk_schmidt +! do j = 1, n +! print*,j,'wk_scc',wk_schmidt(j) +! enddo +! pause +end diff --git a/src/non_hermit_dav/htilde_mat.irp.f b/src/non_hermit_dav/htilde_mat.irp.f new file mode 100644 index 00000000..6d5101ac --- /dev/null +++ b/src/non_hermit_dav/htilde_mat.irp.f @@ -0,0 +1,93 @@ +BEGIN_PROVIDER [ integer, n_mat] + implicit none + n_mat = 2 +END_PROVIDER + + BEGIN_PROVIDER [ double precision, h_non_hermit, (n_mat, n_mat)] +&BEGIN_PROVIDER [ double precision, h_non_hermit_transp, (n_mat, n_mat)] +&BEGIN_PROVIDER [ double precision, reigvec_ht, (n_mat, n_mat)] +&BEGIN_PROVIDER [ double precision, leigvec_ht, (n_mat, n_mat)] +&BEGIN_PROVIDER [ double precision, eigval_ht, (n_mat)] +&BEGIN_PROVIDER [ integer, n_real_ht, (n_mat)] + implicit none + integer :: i,j + do i = 1, n_mat + read(33,*)h_non_hermit(i,1:n_mat) + enddo + print*,'' + print*,'H_mat ' + print*,'' + do i = 1, n_mat + write(*,'(1000(F16.10,X))')h_non_hermit(i,:) + enddo + do i = 1, n_mat + do j = 1, n_mat + h_non_hermit_transp(j,i) = h_non_hermit(i,j) + enddo + enddo + call non_hrmt_real_diag(n_mat,h_non_hermit,reigvec_ht,leigvec_ht,n_real_ht,eigval_ht) + + +END_PROVIDER + + +subroutine hcalc_r_tmp(v,u,N_st,sze) ! v = H u + implicit none + BEGIN_DOC + ! Template of routine for the application of H + ! + ! Here, it is done with the Hamiltonian matrix + ! + ! on the set of determinants of psi_det + ! + ! Computes $v = H | u \rangle$ + ! + END_DOC + integer, intent(in) :: N_st,sze + double precision, intent(in) :: u(sze,N_st) + double precision, intent(inout) :: v(sze,N_st) + integer :: i,j,istate + v = 0.d0 + do istate = 1, N_st + do j = 1, sze + do i = 1, sze + v(i,istate) += h_non_hermit(i,j) * u(j,istate) +! print*,i,j,h_non_hermit(i,j),u(j,istate) + enddo + enddo + enddo + print*,'HU' + do i = 1, sze + print*,v(i,1) + enddo +end + +subroutine hcalc_l_tmp(v,u,N_st,sze) ! v = H^\dagger u + implicit none + BEGIN_DOC + ! Template of routine for the application of H + ! + ! Here, it is done with the Hamiltonian matrix + ! + ! on the set of determinants of psi_det + ! + ! Computes $v = H | u \rangle$ + ! + END_DOC + integer, intent(in) :: N_st,sze + double precision, intent(in) :: u(sze,N_st) + double precision, intent(inout) :: v(sze,N_st) + integer :: i,j,istate + v = 0.d0 + do istate = 1, N_st + do j = 1, sze + do i = 1, sze + v(i,istate) += h_non_hermit_transp(i,j) * u(j,istate) + enddo + enddo + enddo + print*,'HU' + do i = 1, sze + print*,v(i,1) + enddo +end diff --git a/src/non_hermit_dav/lapack_diag_non_hermit.irp.f b/src/non_hermit_dav/lapack_diag_non_hermit.irp.f new file mode 100644 index 00000000..0d652af4 --- /dev/null +++ b/src/non_hermit_dav/lapack_diag_non_hermit.irp.f @@ -0,0 +1,2907 @@ +subroutine lapack_diag_non_sym(n, A, WR, WI, VL, VR) + + BEGIN_DOC + ! You enter with a general non hermitian matrix A(n,n) + ! + ! You get out with the real WR and imaginary part WI of the eigenvalues + ! + ! Eigvalue(n) = WR(n) + i * WI(n) + ! + ! And the left VL and right VR eigenvectors + ! + ! VL(i,j) = :: projection on the basis element |i> on the jth left eigenvector + ! + ! VR(i,j) = :: projection on the basis element |i> on the jth right eigenvector + ! + ! The real part of the matrix A can be written as A = VR D VL^T + ! + END_DOC + + implicit none + + integer, intent(in) :: n + double precision, intent(in) :: A(n,n) + double precision, intent(out) :: WR(n), WI(n), VL(n,n), VR(n,n) + + integer :: lda, ldvl, ldvr, LWORK, INFO + double precision, allocatable :: Atmp(:,:), WORK(:) + + lda = n + ldvl = n + ldvr = n + + allocate( Atmp(n,n) ) + Atmp(1:n,1:n) = A(1:n,1:n) + + allocate(WORK(1)) + LWORK = -1 ! to ask for the optimal size of WORK + call dgeev('V', 'V', n, Atmp, lda, WR, WI, VL, ldvl, VR, ldvr, WORK, LWORK, INFO) + if(INFO.gt.0)then + print*,'dgeev failed !!',INFO + stop + endif + LWORK = max(int(WORK(1)), 1) ! this is the optimal size of WORK + deallocate(WORK) + + allocate(WORK(LWORK)) + + ! Actual diagonalization + call dgeev('V', 'V', n, Atmp, lda, WR, WI, VL, ldvl, VR, ldvr, WORK, LWORK, INFO) + if(INFO.ne.0) then + print*,'dgeev failed !!', INFO + stop + endif + + deallocate(Atmp, WORK) + +end subroutine lapack_diag_non_sym + + +subroutine non_sym_diag_inv_right(n,A,leigvec,reigvec,n_real_eigv,eigval) + implicit none + 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 + integer, intent(in) :: n + double precision, intent(in) :: A(n,n) + double precision, intent(out) :: reigvec(n,n),leigvec(n,n),eigval(n) + double precision, allocatable :: Aw(:,:) + integer, intent(out) :: n_real_eigv + print*,'Computing the left/right eigenvectors ...' + character*1 :: JOBVL,JOBVR + JOBVL = "V" ! computes the left eigenvectors + JOBVR = "V" ! computes the right eigenvectors + double precision, allocatable :: WR(:),WI(:),Vl(:,:),VR(:,:),S(:,:),inv_reigvec(:,:) + integer :: i,j + integer :: n_good + integer, allocatable :: list_good(:), iorder(:) + double precision :: thr + thr = 1.d-10 + ! Eigvalue(n) = WR(n) + i * WI(n) + allocate(WR(n),WI(n),VL(n,n),VR(n,n),Aw(n,n)) + Aw = A + do i = 1, n + do j = i+1, n + if(dabs(Aw(j,j)-Aw(i,i)).lt.thr)then + Aw(j,j)+= thr + Aw(i,i)-= thr +! if(Aw(j,i) * A(i,j) .lt.0d0 )then +! if(dabs(Aw(j,i) * A(i,j)).lt.thr**(1.5d0))then +! print*,Aw(j,j),Aw(i,i) +! print*,Aw(j,i) , A(i,j) + Aw(j,i) = 0.d0 + Aw(i,j) = Aw(j,i) +! endif +! endif + endif + enddo + enddo + call lapack_diag_non_sym(n,Aw,WR,WI,VL,VR) + ! You track the real eigenvalues + n_good = 0 +! do i = 1, n +! write(*,'(100(F16.12,X))')A(:,i) +! enddo + do i = 1, n + print*,'Im part of lambda = ',dabs(WI(i)) + 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) + write(*,'(100(F10.5,X))')VR(:,i) + write(*,'(100(F10.5,X))')VR(:,i+1) + write(*,'(100(F10.5,X))')VL(:,i) + write(*,'(100(F10.5,X))')VL(:,i+1) + 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 + 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) + 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 + allocate(inv_reigvec(n_real_eigv,n_real_eigv)) +! call get_pseudo_inverse(reigvec,n_real_eigv,n_real_eigv,n_real_eigv,inv_reigvec,n_real_eigv,thr) +! do i = 1, n_real_eigv +! do j = 1, n +! leigvec(j,i) = inv_reigvec(i,j) +! enddo +! enddo + allocate( S(n_real_eigv,n_real_eigv) ) + + ! S = VL x VR + call dgemm( 'T', 'N', n_real_eigv, n_real_eigv, n_real_eigv, 1.d0 & + , leigvec, size(leigvec, 1), reigvec, size(reigvec, 1) & + , 0.d0, S, size(S, 1) ) + do i = 1,n_real_eigv + write(*,'(100(F10.5,X))')S(:,i) + enddo +! call lapack_diag_non_sym(n,S,WR,WI,VL,VR) +! print*,'Eigenvalues of S' +! do i = 1, n +! print*,WR(i),dabs(WI(i)) +! enddo + call dgemm( 'T', 'N', n_real_eigv, n_real_eigv, n_real_eigv, 1.d0 & + , leigvec, size(leigvec, 1), reigvec, size(reigvec, 1) & + , 0.d0, S, size(S, 1) ) +! call get_inv_half_svd(S, n_real_eigv, inv_reigvec) + + double precision :: accu_d,accu_nd + accu_nd = 0.d0 + accu_d = 0.d0 + do i = 1, n_real_eigv + do j = 1, n_real_eigv + if(i==j) then + accu_d += S(j,i) * S(j,i) + else + accu_nd = accu_nd + S(j,i) * S(j,i) + endif + enddo + enddo + accu_nd = dsqrt(accu_nd) + + print*,'accu_nd = ',accu_nd + if( accu_nd .lt. 1d-10 ) then + ! L x R is already bi-orthogonal + !print *, ' L & T bi-orthogonality: ok' + return + else + print*,'PB with bi-orthonormality!!' + stop + endif +end + +subroutine lapack_diag_non_sym_new(n, A, WR, WI, VL, VR) + + BEGIN_DOC + ! + ! You enter with a general non hermitian matrix A(n,n) + ! + ! You get out with the real WR and imaginary part WI of the eigenvalues + ! + ! Eigvalue(n) = WR(n) + i * WI(n) + ! + ! And the left VL and right VR eigenvectors + ! + ! VL(i,j) = :: projection on the basis element |i> on the jth left eigenvector + ! + ! VR(i,j) = :: projection on the basis element |i> on the jth right eigenvector + ! + END_DOC + + implicit none + integer, intent(in) :: n + double precision, intent(in) :: A(n,n) + double precision, intent(out) :: WR(n), WI(n), VL(n,n), VR(n,n) + + character*1 :: JOBVL,JOBVR,BALANC,SENSE + integer :: ILO, IHI + integer :: lda, ldvl, ldvr, LWORK, INFO + double precision :: ABNRM + integer, allocatable :: IWORK(:) + double precision, allocatable :: WORK(:), SCALE_array(:), RCONDE(:), RCONDV(:) + double precision, allocatable :: Atmp(:,:) + + allocate( Atmp(n,n) ) + Atmp(1:n,1:n) = A(1:n,1:n) + + JOBVL = "V" ! computes the left eigenvectors + JOBVR = "V" ! computes the right eigenvectors + BALANC = "B" ! Diagonal scaling and Permutation for optimization + SENSE = "B" + lda = n + ldvl = n + ldvr = n + allocate(WORK(1),SCALE_array(n),RCONDE(n),RCONDV(n),IWORK(2*n-2)) + LWORK = -1 ! to ask for the optimal size of WORK + call dgeevx(BALANC,JOBVL,JOBVR,SENSE,& ! CHARACTERS + n,Atmp,lda, & ! MATRIX TO DIAGONALIZE + WR,WI, & ! REAL AND IMAGINARY PART OF EIGENVALUES + VL,ldvl,VR,ldvr, & ! LEFT AND RIGHT EIGENVECTORS + ILO,IHI,SCALE_array,ABNRM,RCONDE,RCONDV, & ! OUTPUTS OF OPTIMIZATION + WORK,LWORK,IWORK,INFO) + + !if(INFO.gt.0)then + ! print*,'dgeev failed !!',INFO + if( INFO.ne.0 ) then + print *, 'dgeevx failed !!', INFO + stop + endif + + LWORK = max(int(work(1)), 1) ! this is the optimal size of WORK + deallocate(WORK) + allocate(WORK(LWORK)) + ! Actual dnon_hrmt_real_diag_newiagonalization + call dgeevx(BALANC,JOBVL,JOBVR,SENSE,& ! CHARACTERS + n,Atmp,lda, & ! MATRIX TO DIAGONALIZE + WR,WI, & ! REAL AND IMAGINARY PART OF EIGENVALUES + VL,ldvl,VR,ldvr, & ! LEFT AND RIGHT EIGENVECTORS + ILO,IHI,SCALE_array,ABNRM,RCONDE,RCONDV, & ! OUTPUTS OF OPTIMIZATION + WORK,LWORK,IWORK,INFO) + + !if(INFO.ne.0)then + ! print*,'dgeev failed !!',INFO + if( INFO.ne.0 ) then + print *, 'dgeevx failed !!', INFO + stop + endif + + deallocate( Atmp ) + deallocate( WORK, SCALE_array, RCONDE, RCONDV, IWORK ) + +end subroutine lapack_diag_non_sym_new + +! --- + +subroutine lapack_diag_non_sym_right(n, A, WR, WI, VR) + + implicit none + + integer, intent(in) :: n + double precision, intent(in) :: A(n,n) + double precision, intent(out) :: WR(n), WI(n), VR(n,n) + + integer :: i, lda, ldvl, ldvr, LWORK, INFO + double precision, allocatable :: Atmp(:,:), WORK(:), VL(:,:) + + lda = n + ldvl = 1 + ldvr = n + + allocate( Atmp(n,n), VL(1,1) ) + Atmp(1:n,1:n) = A(1:n,1:n) + + allocate(WORK(1)) + LWORK = -1 + call dgeev('N', 'V', n, Atmp, lda, WR, WI, VL, ldvl, VR, ldvr, WORK, LWORK, INFO) + if(INFO.gt.0)then + print*,'dgeev failed !!',INFO + stop + endif + + LWORK = max(int(WORK(1)), 1) ! this is the optimal size of WORK + deallocate(WORK) + + allocate(WORK(LWORK)) + + ! Actual diagonalization + call dgeev('N', 'V', n, Atmp, lda, WR, WI, VL, ldvl, VR, ldvr, WORK, LWORK, INFO) + if(INFO.ne.0) then + print*,'dgeev failed !!', INFO + stop + endif + + deallocate(Atmp, WORK, VL) + +! print *, ' JOBL = F' +! print *, ' eigenvalues' +! do i = 1, n +! write(*, '(1000(F16.10,X))') WR(i), WI(i) +! enddo +! print *, ' right eigenvect' +! do i = 1, n +! write(*, '(1000(F16.10,X))') VR(:,i) +! enddo + +end subroutine lapack_diag_non_sym_right + +! --- + +subroutine non_hrmt_real_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, n_good + double precision :: thr, threshold, accu_d, accu_nd + integer, allocatable :: list_good(:), iorder(:) + double precision, allocatable :: Aw(:,:) + double precision, allocatable :: WR(:), WI(:), Vl(:,:), VR(:,:), S(:,:), S_inv_half_tmp(:,:) + + print*, ' Computing the left/right eigenvectors with lapack ...' + + ! Eigvalue(n) = WR(n) + i * WI(n) + allocate(WR(n), WI(n), VL(n,n), VR(n,n), Aw(n,n)) + Aw = A + !print *, ' matrix to diagonalize', Aw + call lapack_diag_non_sym(n, Aw, WR, WI, VL, VR) + + ! --- + ! You track the real eigenvalues + + thr = 1d-15 + + 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 + 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) + 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 + +! print *, ' ordered eigenvalues' +! print *, ' right eigenvect' +! do i = 1, n +! print *, i, eigval(i) +! write(*, '(1000(F16.10,X))') reigvec(:,i) +! enddo + + ! --- + + allocate( S(n_real_eigv,n_real_eigv), S_inv_half_tmp(n_real_eigv,n_real_eigv) ) + + ! S = VL x VR + call dgemm( 'T', 'N', n_real_eigv, n_real_eigv, n_real_eigv, 1.d0 & + , leigvec, size(leigvec, 1), reigvec, size(reigvec, 1) & + , 0.d0, S, size(S, 1) ) + + accu_nd = 0.d0 + accu_d = 0.d0 + do i = 1, n_real_eigv + do j = 1, n_real_eigv + if(i==j) then + accu_d += S(j,i) + else + accu_nd = accu_nd + S(j,i) * S(j,i) + endif + enddo + enddo + accu_nd = dsqrt(accu_nd) + + threshold = 1.d-15 + if( (accu_nd .gt. threshold) .or. (dabs(accu_d-dble(n_real_eigv)) .gt. threshold) ) then + + print*, ' sum of off-diag S elements = ', accu_nd + print*, ' Should be zero ' + print*, ' sum of diag S elements = ', accu_d + print*, ' Should be ',n + print*, ' Not bi-orthonormal !!' + print*, ' Notice that if you are interested in ground state it is not a problem :)' + endif + +end subroutine non_hrmt_real_diag + +! --- + +subroutine lapack_diag_general_non_sym(n, A, B, WR, WI, VL, VR) + + BEGIN_DOC + ! You enter with a general non hermitian matrix A(n,n) and another B(n,n) + ! + ! You get out with the real WR and imaginary part WI of the eigenvalues + ! + ! Eigvalue(n) = (WR(n) + i * WI(n)) + ! + ! And the left VL and right VR eigenvectors + ! + ! VL(i,j) = :: projection on the basis element |i> on the jth left eigenvector + ! + ! VR(i,j) = :: projection on the basis element |i> on the jth right eigenvector + END_DOC + + implicit none + integer, intent(in) :: n + double precision, intent(in) :: A(n,n), B(n,n) + double precision, intent(out) :: WR(n), WI(n), VL(n,n), VR(n,n) + + integer :: lda, ldvl, ldvr, LWORK, INFO + integer :: n_good + double precision, allocatable :: WORK(:) + double precision, allocatable :: Atmp(:,:) + + lda = n + ldvl = n + ldvr = n + + allocate( Atmp(n,n) ) + Atmp(1:n,1:n) = A(1:n,1:n) + + allocate(WORK(1)) + LWORK = -1 + call dgeev('V', 'V', n, Atmp, lda, WR, WI, VL, ldvl, VR, ldvr, WORK, LWORK, INFO) + if(INFO.gt.0) then + print*,'dgeev failed !!',INFO + stop + endif + + LWORK = max(int(WORK(1)), 1) + deallocate(WORK) + + allocate(WORK(LWORK)) + + call dgeev('V', 'V', n, Atmp, lda, WR, WI, VL, ldvl, VR, ldvr, WORK, LWORK, INFO) + if(INFO.ne.0) then + print*,'dgeev failed !!', INFO + stop + endif + + deallocate( WORK, Atmp ) + +end subroutine lapack_diag_general_non_sym + +! --- + +subroutine non_hrmt_general_real_diag(n, A, B, reigvec, leigvec, 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) and B(n,n) + ! + ! A reigvec = eigval * B * reigvec + ! + ! (A)^\dagger leigvec = eigval * B * leigvec + ! + ! 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_good + integer, allocatable :: list_good(:), iorder(:) + double precision, allocatable :: WR(:), WI(:), Vl(:,:), VR(:,:) + double precision, allocatable :: Aw(:,:), Bw(:,:) + + print*,'Computing the left/right eigenvectors ...' + + allocate(WR(n), WI(n), VL(n,n), VR(n,n), Aw(n,n), Bw(n,n)) + Aw = A + Bw = B + + call lapack_diag_general_non_sym(n, A, B, WR, WI, VL, VR) + + ! You track the real eigenvalues + n_good = 0 + do i = 1, n + if(dabs(WI(i)) .lt. 1.d-12) 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-12)then + n_good += 1 + list_good(n_good) = i + eigval(n_good) = WR(i) + endif + enddo + 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) + print*,'n_real_eigv = ', n_real_eigv + print*,'n = ', n + do i = 1, n_real_eigv + print*,i,'eigval(i) = ', eigval(i) + do j = 1, n + reigvec(j,i) = VR(j,list_good(iorder(i))) + leigvec(j,i) = Vl(j,list_good(iorder(i))) + enddo + enddo + +end subroutine non_hrmt_general_real_diag + +! --- + +subroutine impose_biorthog_qr(m, n, thr_d, thr_nd, Vl, Vr) + + implicit none + integer, intent(in) :: m, n + double precision, intent(in) :: thr_d, thr_nd + double precision, intent(inout) :: Vl(m,n), Vr(m,n) + + integer :: i, j + integer :: LWORK, INFO + double precision :: accu_nd, accu_d + double precision, allocatable :: TAU(:), WORK(:) + double precision, allocatable :: S(:,:), R(:,:), tmp(:,:) + + ! --- + + call check_biorthog_binormalize(m, n, Vl, Vr, thr_d, thr_nd, .false.) + + ! --- + + allocate(S(n,n)) + call dgemm( 'T', 'N', n, n, m, 1.d0 & + , Vl, size(Vl, 1), Vr, size(Vr, 1) & + , 0.d0, S, size(S, 1) ) + + accu_nd = 0.d0 + accu_d = 0.d0 + do i = 1, n + do j = 1, n + if(i==j) then + accu_d += S(j,i) + else + accu_nd = accu_nd + S(j,i) * S(j,i) + endif + enddo + enddo + accu_nd = dsqrt(accu_nd) + + if((accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n))/dble(n) .lt. thr_d)) then + print *, ' bi-orthogonal vectors without QR !' + deallocate(S) + return + endif + + ! ------------------------------------------------------------------------------------- + ! QR factorization of S: S = Q x R + + + print *, ' apply QR decomposition ...' + + allocate( TAU(n), WORK(1) ) + + LWORK = -1 + call dgeqrf(n, n, S, n, TAU, WORK, LWORK, INFO) + if(INFO .ne. 0) then + print*,'dgeqrf failed !!', INFO + stop + endif + + LWORK = max(n, int(WORK(1))) + deallocate(WORK) + + allocate( WORK(LWORK) ) + call dgeqrf(n, n, S, n, TAU, WORK, LWORK, INFO) + if(INFO .ne. 0) then + print*,'dgeqrf failed !!', INFO + stop + endif + + ! save the upper triangular R + allocate( R(n,n) ) + R(:,:) = S(:,:) + + ! get Q + LWORK = -1 + call dorgqr(n, n, n, S, n, TAU, WORK, LWORK, INFO) + if(INFO .ne. 0) then + print*,'dorgqr failed !!', INFO + stop + endif + + LWORK = max(n, int(WORK(1))) + deallocate(WORK) + + allocate( WORK(LWORK) ) + call dorgqr(n, n, n, S, n, TAU, WORK, LWORK, INFO) + if(INFO .ne. 0) then + print*,'dorgqr failed !!', INFO + stop + endif + + deallocate( WORK, TAU ) + + ! + ! ------------------------------------------------------------------------------------- + + ! --- + + ! ------------------------------------------------------------------------------------- + ! get bi-orhtog left & right vectors: + ! Vr' = Vr x inv(R) + ! Vl' = inv(Q) x Vl = Q.T x Vl + + ! Q.T x Vl, where Q = S + + allocate( tmp(n,m) ) + call dgemm( 'T', 'T', n, m, n, 1.d0 & + , S, size(S, 1), Vl, size(Vl, 1) & + , 0.d0, tmp, size(tmp, 1) ) + + do i = 1, n + do j = 1, m + Vl(j,i) = tmp(i,j) + enddo + enddo + deallocate(tmp) + + ! --- + + ! inv(R) + !print *, ' inversing upper triangular matrix ...' + call dtrtri("U", "N", n, R, n, INFO) + if(INFO .ne. 0) then + print*,'dtrtri failed !!', INFO + stop + endif + !print *, ' inversing upper triangular matrix OK' + + do i = 1, n-1 + do j = i+1, n + R(j,i) = 0.d0 + enddo + enddo + + !print *, ' inv(R):' + !do i = 1, n + ! write(*, '(1000(F16.10,X))') R(i,:) + !enddo + + ! Vr x inv(R) + allocate( tmp(m,n) ) + call dgemm( 'N', 'N', m, n, n, 1.d0 & + , Vr, size(Vr, 1), R, size(R, 1) & + , 0.d0, tmp, size(tmp, 1) ) + deallocate( R ) + + do i = 1, n + do j = 1, m + Vr(j,i) = tmp(j,i) + enddo + enddo + deallocate(tmp) + + return +end subroutine impose_biorthog_qr + +! --- + +subroutine impose_biorthog_lu(m, n, Vl, Vr, S) + + implicit none + integer, intent(in) :: m, n + double precision, intent(inout) :: Vl(m,n), Vr(m,n), S(n,n) + + integer :: i, j + integer :: INFO + double precision :: nrm + integer, allocatable :: IPIV(:) + double precision, allocatable :: L(:,:), tmp(:,:), vectmp(:) + !double precision, allocatable :: T(:,:), ll(:,:), rr(:,:), tt(:,:) + + !allocate( T(n,n) ) + !T(:,:) = S(:,:) + + print *, ' apply LU decomposition ...' + + ! ------------------------------------------------------------------------------------- + ! LU factorization of S: S = P x L x U + + allocate( IPIV(n) ) + + call dgetrf(n, n, S, n, IPIV, INFO) + if(INFO .ne. 0) then + print*, 'dgetrf failed !!', INFO + stop + endif + + ! check | S - P x L x U | + !allocate( ll(n,n), rr(n,n), tmp(n,n) ) + !ll = S + !rr = S + !do i = 1, n-1 + ! ll(i,i) = 1.d0 + ! do j = i+1, n + ! ll(i,j) = 0.d0 + ! rr(j,i) = 0.d0 + ! enddo + !enddo + !ll(n,n) = 1.d0 + !call dgemm( 'N', 'N', n, n, n, 1.d0 & + ! , ll, size(ll, 1), rr, size(rr, 1) & + ! , 0.d0, tmp, size(tmp, 1) ) + ! deallocate(ll, rr) + !allocate( vectmp(n) ) + !do j = n-1, 1, -1 + ! i = IPIV(j) + ! if(i.ne.j) then + ! print *, j, i + ! vectmp(:) = tmp(i,:) + ! tmp(i,:) = tmp(j,:) + ! tmp(j,:) = vectmp(:) + ! endif + !enddo + !deallocate( vectmp ) + !nrm = 0.d0 + !do i = 1, n + ! do j = 1, n + ! nrm += dabs(tmp(j,i) - T(j,i)) + ! enddo + !enddo + !deallocate( tmp ) + !print*, '|L.T x R - S| =', nrm + !stop + + ! ------ + ! inv(L) + ! ------ + + allocate( L(n,n) ) + L(:,:) = S(:,:) + + call dtrtri("L", "U", n, L, n, INFO) + if(INFO .ne. 0) then + print*, 'dtrtri failed !!', INFO + stop + endif + do i = 1, n-1 + L(i,i) = 1.d0 + do j = i+1, n + L(i,j) = 0.d0 + enddo + enddo + L(n,n) = 1.d0 + + ! ------ + ! inv(U) + ! ------ + + call dtrtri("U", "N", n, S, n, INFO) + if(INFO .ne. 0) then + print*, 'dtrtri failed !!', INFO + stop + endif + + do i = 1, n-1 + do j = i+1, n + S(j,i) = 0.d0 + enddo + enddo + + ! + ! ------------------------------------------------------------------------------------- + + ! --- + + ! ------------------------------------------------------------------------------------- + ! get bi-orhtog left & right vectors: + ! Vr' = Vr x inv(U) + ! Vl' = inv(L) x inv(P) x Vl + + ! inv(P) x Vl + allocate( vectmp(n) ) + do j = n-1, 1, -1 + i = IPIV(j) + if(i.ne.j) then + vectmp(:) = L(:,j) + L(:,j) = L(:,i) + L(:,i) = vectmp(:) + endif + enddo + deallocate( vectmp ) + + ! Vl' + allocate( tmp(m,n) ) + call dgemm( 'N', 'T', m, n, n, 1.d0 & + , Vl, size(Vl, 1), L, size(L, 1) & + , 0.d0, tmp, size(tmp, 1) ) + deallocate(L) + + Vl = tmp + deallocate(tmp) + + ! --- + + ! Vr x inv(U) + allocate( tmp(m,n) ) + call dgemm( 'N', 'N', m, n, n, 1.d0 & + , Vr, size(Vr, 1), S, size(S, 1) & + , 0.d0, tmp, size(tmp, 1) ) + Vr = tmp + deallocate(tmp) + + !allocate( tmp(n,n) ) + !call dgemm( 'T', 'N', n, n, m, 1.d0 & + ! , Vl, size(Vl, 1), Vr, size(Vr, 1) & + ! , 0.d0, tmp, size(tmp, 1) ) + !nrm = 0.d0 + !do i = 1, n + ! do j = 1, n + ! nrm += dabs(tmp(j,i)) + ! enddo + !enddo + !deallocate( tmp ) + !print*, '|L.T x R| =', nrm + !stop + + return +end subroutine impose_biorthog_lu + +! --- + +subroutine check_EIGVEC(n, m, A, eigval, leigvec, reigvec, thr_diag, thr_norm, stop_ifnot) + + implicit none + integer, intent(in) :: n, m + logical, intent(in) :: stop_ifnot + double precision, intent(in) :: A(n,n), eigval(m), leigvec(n,m), reigvec(n,m), thr_diag, thr_norm + + integer :: i, j + double precision :: tmp, tmp_abs, tmp_nrm, tmp_rel, tmp_dif + double precision :: V_nrm, U_nrm + double precision, allocatable :: Mtmp(:,:) + + allocate( Mtmp(n,m) ) + + ! --- + + Mtmp = 0.d0 + call dgemm( 'N', 'N', n, m, n, 1.d0 & + , A, size(A, 1), reigvec, size(reigvec, 1) & + , 0.d0, Mtmp, size(Mtmp, 1) ) + + V_nrm = 0.d0 + tmp_nrm = 0.d0 + tmp_abs = 0.d0 + do j = 1, m + + tmp = 0.d0 + U_nrm = 0.d0 + do i = 1, n + tmp = tmp + dabs(Mtmp(i,j) - eigval(j) * reigvec(i,j)) + tmp_nrm = tmp_nrm + dabs(Mtmp(i,j)) + U_nrm = U_nrm + reigvec(i,j) * reigvec(i,j) + enddo + + tmp_abs = tmp_abs + tmp + V_nrm = V_nrm + U_nrm + !write(*,'(I4,X,(100(F25.16,X)))') j,eigval(j), tmp, U_nrm + + enddo + + if(tmp_abs.lt.10.d-10)then + tmp_rel = thr_diag/10.d0 + else + tmp_rel = tmp_abs / tmp_nrm + endif + tmp_dif = dabs(V_nrm - dble(m)) + + if( stop_ifnot .and. ((tmp_rel .gt. thr_diag) .or. (tmp_dif .gt. thr_norm)) ) then + print *, ' error in right-eigenvectors' + print *, ' err tol = ',thr_diag, thr_norm + print *, '(tmp_rel .gt. thr_diag) = ',(tmp_rel .gt. thr_diag) + print *, '(tmp_dif .gt. thr_norm) = ',(tmp_dif .gt. thr_norm) + print *, ' err estim = ', tmp_abs, tmp_rel + print *, ' CR norm = ', V_nrm + stop + endif + + ! --- + + Mtmp = 0.d0 + call dgemm( 'T', 'N', n, m, n, 1.d0 & + , A, size(A, 1), leigvec, size(leigvec, 1) & + , 0.d0, Mtmp, size(Mtmp, 1) ) + + V_nrm = 0.d0 + tmp_nrm = 0.d0 + tmp_abs = 0.d0 + do j = 1, m + + tmp = 0.d0 + U_nrm = 0.d0 + do i = 1, n + tmp = tmp + dabs(Mtmp(i,j) - eigval(j) * leigvec(i,j)) + tmp_nrm = tmp_nrm + dabs(Mtmp(i,j)) + U_nrm = U_nrm + leigvec(i,j) * leigvec(i,j) + enddo + + tmp_abs = tmp_abs + tmp + V_nrm = V_nrm + U_nrm + !write(*,'(I4,X,(100(F25.16,X)))') j,eigval(j), tmp, U_nrm + + enddo + + if(tmp_abs.lt.10.d-10)then + tmp_rel = thr_diag/10.d0 + else + tmp_rel = tmp_abs / tmp_nrm + endif + if( stop_ifnot .and. ((tmp_rel .gt. thr_diag) .or. (tmp_dif .gt. thr_norm)) ) then + print *, ' error in left-eigenvectors' + print *, ' err tol = ',thr_diag, thr_norm + print *, '(tmp_rel .gt. thr_diag) = ',(tmp_rel .gt. thr_diag) + print *, '(tmp_dif .gt. thr_norm) = ',(tmp_dif .gt. thr_norm) + print *, ' err estim = ', tmp_abs, tmp_rel + print *, ' CR norm = ', V_nrm + stop + endif + + ! --- + + deallocate( Mtmp ) + +end subroutine check_EIGVEC + +! --- + +subroutine check_degen(n, m, eigval, leigvec, reigvec) + + implicit none + integer, intent(in) :: n, m + double precision, intent(in) :: eigval(m) + double precision, intent(inout) :: leigvec(n,m), reigvec(n,m) + + integer :: i, j + double precision :: ei, ej, de, de_thr, accu_nd + double precision, allocatable :: S(:,:) + + de_thr = 1d-6 + + do i = 1, m-1 + ei = eigval(i) + + do j = i+1, m + ej = eigval(j) + de = dabs(ei - ej) + + if(de .lt. de_thr) then + + leigvec(:,i) = 0.d0 + leigvec(:,j) = 0.d0 + leigvec(i,i) = 1.d0 + leigvec(j,j) = 1.d0 + + reigvec(:,i) = 0.d0 + reigvec(:,j) = 0.d0 + reigvec(i,i) = 1.d0 + reigvec(j,j) = 1.d0 + + endif + + enddo + enddo + + ! --- + + allocate( S(m,m) ) + + ! S = VL x VR + call dgemm( 'T', 'N', m, m, n, 1.d0 & + , leigvec, size(leigvec, 1), reigvec, size(reigvec, 1) & + , 0.d0, S, size(S, 1) ) + + accu_nd = 0.d0 + do i = 1, m + do j = 1, m + if(i==j) cycle + accu_nd = accu_nd + S(j,i) * S(j,i) + enddo + enddo + accu_nd = dsqrt(accu_nd) + + deallocate( S ) + + print *, ' check_degen: L & T bi-orthogonality: ok' + print *, ' accu_nd = ', accu_nd + + if( accu_nd .lt. 1d-8 ) then + return + else + stop + endif + +end subroutine check_degen + +! --- + +subroutine impose_weighted_orthog_svd(n, m, W, C) + + implicit none + + integer, intent(in) :: n, m + double precision, intent(inout) :: C(n,m), W(n,n) + + integer :: i, j, num_linear_dependencies + double precision :: threshold + double precision, allocatable :: S(:,:), tmp(:,:) + double precision, allocatable :: U(:,:), Vt(:,:), D(:) + + !print *, ' apply SVD to orthogonalize & normalize weighted vectors' + + ! --- + + ! C.T x W x C + allocate(S(m,m)) + allocate(tmp(m,n)) + call dgemm( 'T', 'N', m, n, n, 1.d0 & + , C, size(C, 1), W, size(W, 1) & + , 0.d0, tmp, size(tmp, 1) ) + call dgemm( 'N', 'N', m, m, n, 1.d0 & + , tmp, size(tmp, 1), C, size(C, 1) & + , 0.d0, S, size(S, 1) ) + deallocate(tmp) + + !print *, ' overlap bef SVD: ' + !do i = 1, m + ! write(*, '(1000(F16.10,X))') S(i,:) + !enddo + + ! --- + + allocate(U(m,m), Vt(m,m), D(m)) + + call svd(S, m, U, m, D, Vt, m, m, m) + + deallocate(S) + + threshold = 1.d-6 + num_linear_dependencies = 0 + do i = 1, m + if(abs(D(i)) <= threshold) then + D(i) = 0.d0 + num_linear_dependencies += 1 + else + ASSERT (D(i) > 0.d0) + D(i) = 1.d0 / dsqrt(D(i)) + endif + enddo + if(num_linear_dependencies > 0) then + write(*,*) ' linear dependencies = ', num_linear_dependencies + write(*,*) ' m = ', m + stop + endif + + ! --- + + allocate(tmp(n,m)) + + ! tmp <-- C x U + call dgemm( 'N', 'N', n, m, m, 1.d0 & + , C, size(C, 1), U, size(U, 1) & + , 0.d0, tmp, size(tmp, 1) ) + + deallocate(U, Vt) + + ! C <-- tmp x sigma^-0.5 + do j = 1, m + do i = 1, n + C(i,j) = tmp(i,j) * D(j) + enddo + enddo + + deallocate(D, tmp) + + ! --- + + ! C.T x W x C + allocate(S(m,m)) + allocate(tmp(m,n)) + call dgemm( 'T', 'N', m, n, n, 1.d0 & + , C, size(C, 1), W, size(W, 1) & + , 0.d0, tmp, size(tmp, 1) ) + call dgemm( 'N', 'N', m, m, n, 1.d0 & + , tmp, size(tmp, 1), C, size(C, 1) & + , 0.d0, S, size(S, 1) ) + deallocate(tmp) + + !print *, ' overlap aft SVD: ' + !do i = 1, m + ! write(*, '(1000(F16.10,X))') S(i,:) + !enddo + + deallocate(S) + + ! --- + +end subroutine impose_weighted_orthog_svd + +! --- + +subroutine impose_orthog_svd(n, m, C) + + implicit none + + integer, intent(in) :: n, m + double precision, intent(inout) :: C(n,m) + + integer :: i, j, num_linear_dependencies + double precision :: threshold + double precision, allocatable :: S(:,:), tmp(:,:) + double precision, allocatable :: U(:,:), Vt(:,:), D(:) + + !print *, ' apply SVD to orthogonalize & normalize vectors' + + ! --- + + allocate(S(m,m)) + + ! S = C.T x C + call dgemm( 'T', 'N', m, m, n, 1.d0 & + , C, size(C, 1), C, size(C, 1) & + , 0.d0, S, size(S, 1) ) + + !print *, ' eigenvec overlap bef SVD: ' + !do i = 1, m + ! write(*, '(1000(F16.10,X))') S(i,:) + !enddo + + ! --- + + allocate(U(m,m), Vt(m,m), D(m)) + + call svd(S, m, U, m, D, Vt, m, m, m) + + deallocate(S) + + threshold = 1.d-6 + num_linear_dependencies = 0 + do i = 1, m + if(abs(D(i)) <= threshold) then + write(*,*) ' D(i) = ', D(i) + D(i) = 0.d0 + num_linear_dependencies += 1 + else + ASSERT (D(i) > 0.d0) + D(i) = 1.d0 / dsqrt(D(i)) + endif + enddo + if(num_linear_dependencies > 0) then + write(*,*) ' linear dependencies = ', num_linear_dependencies + write(*,*) ' m = ', m + write(*,*) ' try with Graham-Schmidt' + stop + endif + + ! --- + + allocate(tmp(n,m)) + + ! tmp <-- C x U + call dgemm( 'N', 'N', n, m, m, 1.d0 & + , C, size(C, 1), U, size(U, 1) & + , 0.d0, tmp, size(tmp, 1) ) + + deallocate(U, Vt) + + ! C <-- tmp x sigma^-0.5 + do j = 1, m + do i = 1, n + C(i,j) = tmp(i,j) * D(j) + enddo + enddo + + deallocate(D, tmp) + + ! --- + + allocate(S(m,m)) + + ! S = C.T x C + call dgemm( 'T', 'N', m, m, n, 1.d0 & + , C, size(C, 1), C, size(C, 1) & + , 0.d0, S, size(S, 1) ) + + !print *, ' eigenvec overlap aft SVD: ' + !do i = 1, m + ! write(*, '(1000(F16.10,X))') S(i,:) + !enddo + + deallocate(S) + + ! --- + +end subroutine impose_orthog_svd + +! --- + +subroutine impose_orthog_svd_overlap(n, m, C, overlap) + + implicit none + + integer, intent(in) :: n, m + double precision, intent(in ) :: overlap(n,n) + double precision, intent(inout) :: C(n,m) + + integer :: i, j, num_linear_dependencies + double precision :: threshold + double precision, allocatable :: S(:,:), tmp(:,:), Stmp(:,:) + double precision, allocatable :: U(:,:), Vt(:,:), D(:) + + print *, ' apply SVD to orthogonalize vectors' + + ! --- + + ! S = C.T x overlap x C + allocate(S(m,m), Stmp(n,m)) + call dgemm( 'N', 'N', n, m, n, 1.d0 & + , overlap, size(overlap, 1), C, size(C, 1) & + , 0.d0, Stmp, size(Stmp, 1) ) + call dgemm( 'T', 'N', m, m, n, 1.d0 & + , C, size(C, 1), Stmp, size(Stmp, 1) & + , 0.d0, S, size(S, 1) ) + deallocate(Stmp) + + !print *, ' eigenvec overlap bef SVD: ' + !do i = 1, m + ! write(*, '(1000(F16.10,X))') S(i,:) + !enddo + + ! --- + + allocate(U(m,m), Vt(m,m), D(m)) + + call svd(S, m, U, m, D, Vt, m, m, m) + + deallocate(S) + + threshold = 1.d-6 + num_linear_dependencies = 0 + do i = 1, m + if(abs(D(i)) <= threshold) then + D(i) = 0.d0 + num_linear_dependencies += 1 + else + ASSERT (D(i) > 0.d0) + D(i) = 1.d0 / dsqrt(D(i)) + endif + enddo + if(num_linear_dependencies > 0) then + write(*,*) ' linear dependencies = ', num_linear_dependencies + write(*,*) ' m = ', m + stop + endif + + ! --- + + allocate(tmp(n,m)) + + ! tmp <-- C x U + call dgemm( 'N', 'N', n, m, m, 1.d0 & + , C, size(C, 1), U, size(U, 1) & + , 0.d0, tmp, size(tmp, 1) ) + + deallocate(U, Vt) + + ! C <-- tmp x sigma^-0.5 + do j = 1, m + do i = 1, n + C(i,j) = tmp(i,j) * D(j) + enddo + enddo + + deallocate(D, tmp) + + ! --- + + ! S = C.T x overlap x C + allocate(S(m,m), Stmp(n,m)) + call dgemm( 'N', 'N', n, m, n, 1.d0 & + , overlap, size(overlap, 1), C, size(C, 1) & + , 0.d0, Stmp, size(Stmp, 1) ) + call dgemm( 'T', 'N', m, m, n, 1.d0 & + , C, size(C, 1), Stmp, size(Stmp, 1) & + , 0.d0, S, size(S, 1) ) + deallocate(Stmp) + + !print *, ' eigenvec overlap aft SVD: ' + !do i = 1, m + ! write(*, '(1000(F16.10,X))') S(i,:) + !enddo + deallocate(S) + +end subroutine impose_orthog_svd_overlap + +! --- + +subroutine impose_orthog_GramSchmidt(n, m, C) + + implicit none + + integer, intent(in) :: n, m + double precision, intent(inout) :: C(n,m) + + integer :: i, j, k + double precision :: Ojk, Ojj, fact_ct + double precision, allocatable :: S(:,:) + + print *, '' + print *, ' apply Gram-Schmidt to orthogonalize & normalize vectors' + print *, '' + + ! --- + + allocate(S(m,m)) + call dgemm( 'T', 'N', m, m, n, 1.d0 & + , C, size(C, 1), C, size(C, 1) & + , 0.d0, S, size(S, 1) ) + + print *, ' eigenvec overlap bef Gram-Schmidt: ' + do i = 1, m + write(*, '(1000(F16.10,X))') S(i,:) + enddo + + ! --- + + do k = 2, m + do j = 1, k-1 + + Ojk = 0.d0 + Ojj = 0.d0 + do i = 1, n + Ojk = Ojk + C(i,j) * C(i,k) + Ojj = Ojj + C(i,j) * C(i,j) + enddo + fact_ct = Ojk / Ojj + + do i = 1, n + C(i,k) = C(i,k) - fact_ct * C(i,j) + enddo + + enddo + enddo + + do k = 1, m + fact_ct = 0.d0 + do i = 1, n + fact_ct = fact_ct + C(i,k) * C(i,k) + enddo + fact_ct = dsqrt(fact_ct) + do i = 1, n + C(i,k) = C(i,k) / fact_ct + enddo + enddo + + ! --- + + call dgemm( 'T', 'N', m, m, n, 1.d0 & + , C, size(C, 1), C, size(C, 1) & + , 0.d0, S, size(S, 1) ) + + print *, ' eigenvec overlap aft Gram-Schmidt: ' + do i = 1, m + write(*, '(1000(F16.10,X))') S(i,:) + enddo + + deallocate(S) + + ! --- + +end subroutine impose_orthog_GramSchmidt + +! --- + +subroutine impose_orthog_ones(n, deg_num, C) + + + implicit none + + integer, intent(in) :: n + integer, intent(in) :: deg_num(n) + double precision, intent(inout) :: C(n,n) + + integer :: i, j, ii, di, dj + + print *, '' + print *, ' orthogonalize vectors by hand' + print *, '' + + do i = 1, n-1 + di = deg_num(i) + + if(di .gt. 1) then + + do ii = 1, di + C(: ,i+ii-1) = 0.d0 + C(i+ii-1,i+ii-1) = 1.d0 + enddo + + do j = i+di+1, n + dj = deg_num(j) + if(dj .eq. di) then + do ii = 1, dj + C(:, j+ii-1) = 0.d0 + C(j+ii-1,j+ii-1) = 1.d0 + enddo + endif + enddo + + endif + enddo + +end subroutine impose_orthog_ones + +! --- + +subroutine impose_orthog_degen_eigvec(n, e0, C0) + + implicit none + + integer, intent(in) :: n + double precision, intent(in) :: e0(n) + double precision, intent(inout) :: C0(n,n) + + integer :: i, j, k, m + double precision :: ei, ej, de, de_thr + integer, allocatable :: deg_num(:) + double precision, allocatable :: C(:,:) + + ! --- + + allocate( deg_num(n) ) + do i = 1, n + deg_num(i) = 1 + enddo + + de_thr = thr_degen_tc + + do i = 1, n-1 + ei = e0(i) + + ! already considered in degen vectors + if(deg_num(i).eq.0) cycle + + do j = i+1, n + ej = e0(j) + de = dabs(ei - ej) + + if(de .lt. de_thr) then + deg_num(i) = deg_num(i) + 1 + deg_num(j) = 0 + endif + + enddo + enddo + + + !do i = 1, n + ! if(deg_num(i) .gt. 1) then + ! print *, ' degen on', i, deg_num(i) + ! endif + !enddo + + ! --- + +! call impose_orthog_ones(n, deg_num, C0) + + do i = 1, n + m = deg_num(i) + + if(m .gt. 1) then + !if(m.eq.3) then + + allocate(C(n,m)) + do j = 1, m + C(1:n,j) = C0(1:n,i+j-1) + enddo + + ! --- + + ! C <= C U sigma^-0.5 + call impose_orthog_svd(n, m, C) + + ! --- + + ! C = I + !C = 0.d0 + !do j = 1, m + ! C(i+j-1,j) = 1.d0 + !enddo + + ! --- + +! call impose_orthog_GramSchmidt(n, m, C) + + ! --- + + do j = 1, m + C0(1:n,i+j-1) = C(1:n,j) + enddo + deallocate(C) + + endif + enddo + +end subroutine impose_orthog_degen_eigvec + +! --- + +subroutine get_halfinv_svd(n, S) + + implicit none + + integer, intent(in) :: n + double precision, intent(inout) :: S(n,n) + + integer :: num_linear_dependencies + integer :: i, j, k + double precision :: accu_d, accu_nd, thresh + double precision, parameter :: threshold = 1.d-6 + double precision, allocatable :: U(:,:), Vt(:,:), D(:) + double precision, allocatable :: S0(:,:), Stmp(:,:), Stmp2(:,:) + + allocate( S0(n,n) ) + S0(1:n,1:n) = S(1:n,1:n) + + allocate(U(n,n), Vt(n,n), D(n)) + call svd(S, n, U, n, D, Vt, n, n, n) + + num_linear_dependencies = 0 + do i = 1, n + if(abs(D(i)) <= threshold) then + D(i) = 0.d0 + num_linear_dependencies += 1 + else + ASSERT (D(i) > 0.d0) + D(i) = 1.d0 / dsqrt(D(i)) + endif + enddo + write(*,*) ' linear dependencies', num_linear_dependencies + + S(:,:) = 0.d0 + do k = 1, n + if(D(k) /= 0.d0) then + do j = 1, n + do i = 1, n + S(i,j) = S(i,j) + U(i,k) * D(k) * Vt(k,j) + enddo + enddo + endif + enddo + deallocate(U, D, Vt) + + allocate( Stmp(n,n), Stmp2(n,n) ) + Stmp = 0.d0 + Stmp2 = 0.d0 + ! S^-1/2 x S + call dgemm( 'N', 'N', n, n, n, 1.d0 & + , S, size(S, 1), S0, size(S0, 1) & + , 0.d0, Stmp, size(Stmp, 1) ) + ! ( S^-1/2 x S ) x S^-1/2 + call dgemm( 'N', 'N', n, n, n, 1.d0 & + , Stmp, size(Stmp, 1), S, size(S, 1) & + , 0.d0, Stmp2, size(Stmp2, 1) ) + + accu_nd = 0.d0 + accu_d = 0.d0 + thresh = 1.d-10 + do i = 1, n + do j = 1, n + if(i==j) then + accu_d += Stmp2(j,i) + else + accu_nd = accu_nd + Stmp2(j,i) * Stmp2(j,i) + endif + enddo + enddo + accu_nd = dsqrt(accu_nd) + if( accu_nd.gt.thresh .or. dabs(accu_d-dble(n)).gt.thresh) then + print*, ' after S^-1/2: sum of off-diag S elements = ', accu_nd + print*, ' after S^-1/2: sum of diag S elements = ', accu_d + do i = 1, n + write(*,'(1000(F16.10,X))') Stmp2(i,:) + enddo + stop + endif + + deallocate(S0, Stmp, Stmp2) + +end subroutine get_halfinv_svd + +! --- + +subroutine check_biorthog_binormalize(n, m, Vl, Vr, thr_d, thr_nd, stop_ifnot) + + implicit none + + integer, intent(in) :: n, m + logical, intent(in) :: stop_ifnot + double precision, intent(in) :: thr_d, thr_nd + double precision, intent(inout) :: Vl(n,m), Vr(n,m) + + integer :: i, j + double precision :: accu_d, accu_nd, s_tmp + double precision, allocatable :: S(:,:) + + !print *, ' check bi-orthonormality' + + ! --- + + allocate(S(m,m)) + call dgemm( 'T', 'N', m, m, n, 1.d0 & + , Vl, size(Vl, 1), Vr, size(Vr, 1) & + , 0.d0, S, size(S, 1) ) + !print *, ' overlap matrix before:' + !do i = 1, m + ! write(*,'(1000(F16.10,X))') S(i,:) + !enddo + + ! S(i,i) = -1 + do i = 1, m + if(S(i,i) .lt. 0.d0) then + !if( (S(i,i) + 1.d0) .lt. thr_d ) then + do j = 1, n + Vl(j,i) = -1.d0 * Vl(j,i) + enddo + !S(i,i) = 1.d0 + S(i,i) = -S(i,i) + endif + enddo + + accu_d = 0.d0 + accu_nd = 0.d0 + do i = 1, m + do j = 1, m + if(i==j) then + accu_d = accu_d + S(i,i) + else + accu_nd = accu_nd + S(j,i) * S(j,i) + endif + enddo + enddo + accu_nd = dsqrt(accu_nd) / dble(m) + !print*, ' diag acc bef = ', accu_d + !print*, ' nondiag acc bef = ', accu_nd + + ! --- + + if( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(m))/dble(m) .gt. thr_d) ) then + + do i = 1, m + if(S(i,i) <= 0.d0) then + print *, ' overap negative' + print *, i, S(i,i) + exit + endif + if(dabs(S(i,i) - 1.d0) .gt. thr_d) then + s_tmp = 1.d0 / dsqrt(S(i,i)) + do j = 1, n + Vl(j,i) = Vl(j,i) * s_tmp + Vr(j,i) = Vr(j,i) * s_tmp + enddo + endif + + enddo + + endif + + ! --- + + call dgemm( 'T', 'N', m, m, n, 1.d0 & + , Vl, size(Vl, 1), Vr, size(Vr, 1) & + , 0.d0, S, size(S, 1) ) + !print *, ' overlap matrix after:' + !do i = 1, m + ! write(*,'(1000(F16.10,X))') S(i,:) + !enddo + + accu_d = 0.d0 + accu_nd = 0.d0 + do i = 1, m + do j = 1, m + if(i==j) then + accu_d = accu_d + S(i,i) + else + accu_nd = accu_nd + S(j,i) * S(j,i) + endif + enddo + enddo + accu_nd = dsqrt(accu_nd) / dble(m) + !print *, ' diag acc aft = ', accu_d + !print *, ' nondiag acc aft = ', accu_nd + + deallocate(S) + + ! --- + + if( stop_ifnot .and. ((accu_nd .gt. thr_nd) .or. (dabs(accu_d-dble(m))/dble(m) .gt. thr_d)) ) then + print *, accu_nd, thr_nd + print *, dabs(accu_d-dble(m))/dble(m), thr_d + print *, ' biorthog_binormalize failed !' + stop + endif + +end subroutine check_biorthog_binormalize + +! --- + +subroutine check_weighted_biorthog(n, m, W, Vl, Vr, thr_d, thr_nd, accu_d, accu_nd, S, stop_ifnot) + + implicit none + + integer, intent(in) :: n, m + double precision, intent(in) :: Vl(n,m), Vr(n,m), W(n,n) + double precision, intent(in) :: thr_d, thr_nd + logical, intent(in) :: stop_ifnot + double precision, intent(out) :: accu_d, accu_nd, S(m,m) + + integer :: i, j + double precision, allocatable :: SS(:,:), tmp(:,:) + + print *, ' check weighted bi-orthogonality' + + ! --- + + allocate(tmp(m,n)) + call dgemm( 'T', 'N', m, n, n, 1.d0 & + , Vl, size(Vl, 1), W, size(W, 1) & + , 0.d0, tmp, size(tmp, 1) ) + call dgemm( 'N', 'N', m, m, n, 1.d0 & + , tmp, size(tmp, 1), Vr, size(Vr, 1) & + , 0.d0, S, size(S, 1) ) + deallocate(tmp) + + !print *, ' overlap matrix:' + !do i = 1, m + ! write(*,'(1000(F16.10,X))') S(i,:) + !enddo + + accu_d = 0.d0 + accu_nd = 0.d0 + do i = 1, m + do j = 1, m + if(i==j) then + accu_d = accu_d + dabs(S(i,i)) + else + accu_nd = accu_nd + S(j,i) * S(j,i) + endif + enddo + enddo + accu_nd = dsqrt(accu_nd) + + print *, ' accu_nd = ', accu_nd + print *, ' accu_d = ', dabs(accu_d-dble(m))/dble(m) + + ! --- + + if( stop_ifnot .and. ((accu_nd .gt. thr_nd) .or. dabs(accu_d-dble(m))/dble(m) .gt. thr_d) ) then + print *, ' non bi-orthogonal vectors !' + print *, ' accu_nd = ', accu_nd + print *, ' accu_d = ', dabs(accu_d-dble(m))/dble(m) + !print *, ' overlap matrix:' + !do i = 1, m + ! write(*,'(1000(F16.10,X))') S(i,:) + !enddo + stop + endif + +end subroutine check_weighted_biorthog + +! --- + +subroutine check_biorthog(n, m, Vl, Vr, accu_d, accu_nd, S, thr_d, thr_nd, stop_ifnot) + + implicit none + + integer, intent(in) :: n, m + double precision, intent(in) :: Vl(n,m), Vr(n,m) + logical, intent(in) :: stop_ifnot + double precision, intent(in) :: thr_d, thr_nd + double precision, intent(out) :: accu_d, accu_nd, S(m,m) + + integer :: i, j + double precision, allocatable :: SS(:,:) + + !print *, ' check bi-orthogonality' + + ! --- + + call dgemm( 'T', 'N', m, m, n, 1.d0 & + , Vl, size(Vl, 1), Vr, size(Vr, 1) & + , 0.d0, S, size(S, 1) ) + + !print *, ' overlap matrix:' + !do i = 1, m + ! write(*,'(1000(F16.10,X))') S(i,:) + !enddo + + accu_d = 0.d0 + accu_nd = 0.d0 + do i = 1, m + do j = 1, m + if(i==j) then + accu_d = accu_d + dabs(S(i,i)) + else + accu_nd = accu_nd + S(j,i) * S(j,i) + endif + enddo + enddo + accu_nd = dsqrt(accu_nd) / dble(m) + + !print *, ' accu_nd = ', accu_nd + !print *, ' accu_d = ', dabs(accu_d-dble(m))/dble(m) + + ! --- + + if(stop_ifnot .and. ((accu_nd .gt. thr_nd) .or. dabs(accu_d-dble(m))/dble(m) .gt. thr_d)) then + print *, ' non bi-orthogonal vectors !' + print *, ' accu_nd = ', accu_nd + print *, ' accu_d = ', dabs(accu_d-dble(m))/dble(m) + !print *, ' overlap matrix:' + !do i = 1, m + ! write(*,'(1000(F16.10,X))') S(i,:) + !enddo + stop + endif + +end subroutine check_biorthog + +! --- + +subroutine check_orthog(n, m, V, accu_d, accu_nd, S) + + implicit none + + integer, intent(in) :: n, m + double precision, intent(in) :: V(n,m) + double precision, intent(out) :: accu_d, accu_nd, S(m,m) + + integer :: i, j + + S = 0.d0 + call dgemm( 'T', 'N', m, m, n, 1.d0 & + , V, size(V, 1), V, size(V, 1) & + , 0.d0, S, size(S, 1) ) + + !print *, '' + !print *, ' overlap matrix:' + !do i = 1, m + ! write(*,'(1000(F16.10,X))') S(i,:) + !enddo + !print *, '' + + accu_d = 0.d0 + accu_nd = 0.d0 + do i = 1, m + do j = 1, m + if(i==j) then + accu_d = accu_d + dabs(S(i,i)) + else + accu_nd = accu_nd + S(j,i) * S(j,i) + endif + enddo + enddo + accu_nd = dsqrt(accu_nd) + + !print*, ' diag acc: ', accu_d + !print*, ' nondiag acc: ', accu_nd + +end subroutine check_orthog + +! --- + +subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0) + + implicit none + + integer, intent(in) :: n + double precision, intent(in) :: e0(n) + double precision, intent(inout) :: L0(n,n), R0(n,n) + + logical :: complex_root + integer :: i, j, k, m + double precision :: ei, ej, de, de_thr + double precision :: accu_d, accu_nd + integer, allocatable :: deg_num(:) + double precision, allocatable :: L(:,:), R(:,:), S(:,:), S_inv_half(:,:) + + ! --- + + allocate( deg_num(n) ) + do i = 1, n + deg_num(i) = 1 + enddo + + de_thr = thr_degen_tc + + do i = 1, n-1 + ei = e0(i) + + ! already considered in degen vectors + if(deg_num(i).eq.0) cycle + + do j = i+1, n + ej = e0(j) + de = dabs(ei - ej) + + if(de .lt. de_thr) then + deg_num(i) = deg_num(i) + 1 + deg_num(j) = 0 + endif + + enddo + enddo + + !do i = 1, n + ! if(deg_num(i) .gt. 1) then + ! print *, ' degen on', i, deg_num(i), e0(i) + ! endif + !enddo + + ! --- + + do i = 1, n + m = deg_num(i) + + if(m .gt. 1) then + + allocate(L(n,m)) + allocate(R(n,m)) + + do j = 1, m + L(1:n,j) = L0(1:n,i+j-1) + R(1:n,j) = R0(1:n,i+j-1) + enddo + + ! --- + + call impose_orthog_svd(n, m, L) + call impose_orthog_svd(n, m, R) + !call impose_orthog_GramSchmidt(n, m, L) + !call impose_orthog_GramSchmidt(n, m, R) + + ! --- + + !allocate(S(m,m)) + !call dgemm( 'T', 'N', m, m, n, 1.d0 & + ! , L, size(L, 1), R, size(R, 1) & + ! , 0.d0, S, size(S, 1) ) + !allocate(S_inv_half(m,m)) + !call get_inv_half_nonsymmat_diago(S, m, S_inv_half, complex_root) + !if(complex_root) then + ! print*, ' complex roots in inv_half !!! ' + ! stop + !endif + !call bi_ortho_s_inv_half(m, L, R, S_inv_half) + !deallocate(S, S_inv_half) + + call impose_biorthog_svd(n, m, L, R) + + !call impose_biorthog_qr(n, m, thr_d, thr_nd, L, R) + + ! --- + + do j = 1, m + L0(1:n,i+j-1) = L(1:n,j) + R0(1:n,i+j-1) = R(1:n,j) + enddo + + deallocate(L, R) + + endif + enddo + +end subroutine impose_biorthog_degen_eigvec + +! --- + +subroutine impose_orthog_biorthog_degen_eigvec(n, thr_d, thr_nd, e0, L0, R0) + + implicit none + + integer, intent(in) :: n + double precision, intent(in) :: thr_d, thr_nd + double precision, intent(in) :: e0(n) + double precision, intent(inout) :: L0(n,n), R0(n,n) + + integer :: i, j, k, m + double precision :: ei, ej, de, de_thr + double precision :: accu_d, accu_nd + integer, allocatable :: deg_num(:) + double precision, allocatable :: L(:,:), R(:,:), S(:,:) + + ! --- + + allocate( deg_num(n) ) + do i = 1, n + deg_num(i) = 1 + enddo + + de_thr = thr_degen_tc + + do i = 1, n-1 + ei = e0(i) + + ! already considered in degen vectors + if(deg_num(i).eq.0) cycle + + do j = i+1, n + ej = e0(j) + de = dabs(ei - ej) + + if(de .lt. de_thr) then + deg_num(i) = deg_num(i) + 1 + deg_num(j) = 0 + endif + + enddo + enddo + + do i = 1, n + if(deg_num(i).gt.1) then + print *, ' degen on', i, deg_num(i) + endif + enddo + + ! --- + + do i = 1, n + m = deg_num(i) + + if(m .gt. 1) then + + allocate(L(n,m)) + allocate(R(n,m)) + + do j = 1, m + L(1:n,j) = L0(1:n,i+j-1) + R(1:n,j) = R0(1:n,i+j-1) + enddo + + ! --- + + call impose_orthog_svd(n, m, L) + call impose_orthog_svd(n, m, R) + + ! --- + + call impose_biorthog_qr(n, m, thr_d, thr_nd, L, R) + + allocate(S(m,m)) + call check_biorthog(n, m, L, R, accu_d, accu_nd, S, thr_d, thr_nd, .true.) + !call check_biorthog(n, m, L, L, accu_d, accu_nd, S, thr_d, thr_nd, .true.) + !call check_biorthog(n, m, R, R, accu_d, accu_nd, S, thr_d, thr_nd, .false.) + deallocate(S) + + ! --- + + do j = 1, m + L0(1:n,i+j-1) = L(1:n,j) + R0(1:n,i+j-1) = R(1:n,j) + enddo + + deallocate(L, R) + + endif + enddo + +end subroutine impose_orthog_biorthog_degen_eigvec + +! --- + +subroutine impose_unique_biorthog_degen_eigvec(n, thr_d, thr_nd, e0, C0, W0, L0, R0) + + implicit none + + integer, intent(in) :: n + double precision, intent(in) :: thr_d, thr_nd + double precision, intent(in) :: e0(n), W0(n,n), C0(n,n) + double precision, intent(inout) :: L0(n,n), R0(n,n) + + logical :: complex_root + integer :: i, j, k, m + double precision :: ei, ej, de, de_thr + integer, allocatable :: deg_num(:) + double precision, allocatable :: L(:,:), R(:,:), C(:,:) + double precision, allocatable :: S(:,:), S_inv_half(:,:), tmp(:,:) + + ! --- + + allocate( deg_num(n) ) + do i = 1, n + deg_num(i) = 1 + enddo + + de_thr = thr_degen_tc + + do i = 1, n-1 + ei = e0(i) + + ! already considered in degen vectors + if(deg_num(i).eq.0) cycle + + do j = i+1, n + ej = e0(j) + de = dabs(ei - ej) + + if(de .lt. de_thr) then + deg_num(i) = deg_num(i) + 1 + deg_num(j) = 0 + endif + + enddo + enddo + + !do i = 1, n + ! if(deg_num(i) .gt. 1) then + ! print *, ' degen on', i, deg_num(i) + ! endif + !enddo + + ! --- + + do i = 1, n + m = deg_num(i) + + if(m .gt. 1) then + + allocate(L(n,m)) + allocate(R(n,m)) + allocate(C(n,m)) + + do j = 1, m + L(1:n,j) = L0(1:n,i+j-1) + R(1:n,j) = R0(1:n,i+j-1) + C(1:n,j) = C0(1:n,i+j-1) + enddo + + ! --- + + call impose_orthog_svd(n, m, L) + call impose_orthog_svd(n, m, R) + + ! --- + + + ! TODO: + ! select C correctly via overlap + ! or via selecting degen in HF + + !call max_overlap_qr(n, m, C, L) + !call max_overlap_qr(n, m, C, R) + + + allocate(tmp(m,n)) + allocate(S(m,m)) + + call dgemm( 'T', 'N', m, n, n, 1.d0 & + , L, size(L, 1), W0, size(W0, 1) & + , 0.d0, tmp, size(tmp, 1) ) + call dgemm( 'N', 'N', m, m, n, 1.d0 & + , tmp, size(tmp, 1), C, size(C, 1) & + , 0.d0, S, size(S, 1) ) + + call max_overlap_qr(n, m, S, L) + !call max_overlap_invprod(n, m, S, L) + + call dgemm( 'T', 'N', m, n, n, 1.d0 & + , C, size(C, 1), W0, size(W0, 1) & + , 0.d0, tmp, size(tmp, 1) ) + call dgemm( 'N', 'N', m, m, n, 1.d0 & + , tmp, size(tmp, 1), R, size(R, 1) & + , 0.d0, S, size(S, 1) ) + + call max_overlap_qr(n, m, S, R) + !call max_overlap_invprod(n, m, S, R) + + deallocate(S, tmp) + + ! --- + + allocate(S(m,m), S_inv_half(m,m)) + call dgemm( 'T', 'N', m, m, n, 1.d0 & + , L, size(L, 1), R, size(R, 1) & + , 0.d0, S, size(S, 1) ) + call get_inv_half_nonsymmat_diago(S, m, S_inv_half, complex_root) + if(complex_root)then + call impose_biorthog_svd(n, m, L, R) + !call impose_biorthog_qr(n, m, thr_d, thr_nd, L, R) + else + call bi_ortho_s_inv_half(m, L, R, S_inv_half) + endif + deallocate(S, S_inv_half) + + ! --- + + do j = 1, m + L0(1:n,i+j-1) = L(1:n,j) + R0(1:n,i+j-1) = R(1:n,j) + enddo + + deallocate(L, R, C) + + endif + enddo + +end subroutine impose_unique_biorthog_degen_eigvec + +! --- + +subroutine max_overlap_qr(m, n, S0, V) + + implicit none + integer, intent(in) :: m, n + double precision, intent(in) :: S0(n,n) + double precision, intent(inout) :: V(m,n) + + integer :: i, j + integer :: LWORK, INFO + double precision, allocatable :: TAU(:), WORK(:) + double precision, allocatable :: S(:,:), tmp(:,:) + + allocate(S(n,n)) + S = S0 + + ! --- + + allocate( TAU(n), WORK(1) ) + + LWORK = -1 + call dgeqrf(n, n, S, n, TAU, WORK, LWORK, INFO) + if(INFO .ne. 0) then + print*,'dgeqrf failed !!', INFO + stop + endif + + LWORK = max(n, int(WORK(1))) + deallocate(WORK) + + allocate( WORK(LWORK) ) + call dgeqrf(n, n, S, n, TAU, WORK, LWORK, INFO) + if(INFO .ne. 0) then + print*,'dgeqrf failed !!', INFO + stop + endif + + ! get Q in S matrix + LWORK = -1 + call dorgqr(n, n, n, S, n, TAU, WORK, LWORK, INFO) + if(INFO .ne. 0) then + print*,'dorgqr failed !!', INFO + stop + endif + + LWORK = max(n, int(WORK(1))) + deallocate(WORK) + + allocate( WORK(LWORK) ) + call dorgqr(n, n, n, S, n, TAU, WORK, LWORK, INFO) + if(INFO .ne. 0) then + print*,'dorgqr failed !!', INFO + stop + endif + + deallocate( WORK, TAU ) + + ! --- + + ! V0.T <-- Q.T x V0.T, where Q = S + + allocate( tmp(n,m) ) + + call dgemm( 'T', 'T', n, m, n, 1.d0 & + , S, size(S, 1), V, size(V, 1) & + , 0.d0, tmp, size(tmp, 1) ) + + deallocate(S) + + do i = 1, n + do j = 1, m + V(j,i) = tmp(i,j) + enddo + enddo + + deallocate(tmp) + + ! --- + + return +end subroutine max_overlap_qr + +! --- + +subroutine max_overlap_invprod(n, m, S, V) + + implicit none + integer, intent(in) :: m, n + double precision, intent(in) :: S(m,m) + double precision, intent(inout) :: V(n,m) + + integer :: i + double precision, allocatable :: invS(:,:), tmp(:,:) + + allocate(invS(m,m)) + call get_inverse(S, size(S, 1), m, invS, size(invS, 1)) + print *, ' overlap ' + do i = 1, m + write(*, '(1000(F16.10,X))') S(i,:) + enddo + print *, ' inv overlap ' + do i = 1, m + write(*, '(1000(F16.10,X))') invS(i,:) + enddo + + allocate(tmp(n,m)) + tmp = V + + call dgemm( 'N', 'N', n, m, m, 1.d0 & + , tmp, size(tmp, 1), invS, size(invS, 1) & + , 0.d0, V, size(V, 1) ) + + deallocate(tmp, invS) + + return +end subroutine max_overlap_invprod + +! --- + +subroutine impose_biorthog_svd(n, m, L, R) + + implicit none + + integer, intent(in) :: n, m + double precision, intent(inout) :: L(n,m), R(n,m) + + integer :: i, j, num_linear_dependencies + double precision :: threshold + double precision, allocatable :: S(:,:), tmp(:,:) + double precision, allocatable :: U(:,:), V(:,:), Vt(:,:), D(:) + + ! --- + + allocate(S(m,m)) + + call dgemm( 'T', 'N', m, m, n, 1.d0 & + , L, size(L, 1), R, size(R, 1) & + , 0.d0, S, size(S, 1) ) + + !print *, ' overlap bef SVD: ' + !do i = 1, m + ! write(*, '(1000(F16.10,X))') S(i,:) + !enddo + + ! --- + + allocate(U(m,m), Vt(m,m), D(m)) + + call svd(S, m, U, m, D, Vt, m, m, m) + + deallocate(S) + + threshold = 1.d-6 + num_linear_dependencies = 0 + do i = 1, m + if(abs(D(i)) <= threshold) then + D(i) = 0.d0 + num_linear_dependencies += 1 + else + ASSERT (D(i) > 0.d0) + D(i) = 1.d0 / dsqrt(D(i)) + endif + enddo + if(num_linear_dependencies > 0) then + write(*,*) ' linear dependencies = ', num_linear_dependencies + write(*,*) ' m = ', m + stop + endif + + allocate(V(m,m)) + do i = 1, m + do j = 1, m + V(j,i) = Vt(i,j) + enddo + enddo + deallocate(Vt) + + ! --- + + allocate(tmp(n,m)) + + ! tmp <-- R x V + call dgemm( 'N', 'N', n, m, m, 1.d0 & + , R, size(R, 1), V, size(V, 1) & + , 0.d0, tmp, size(tmp, 1) ) + deallocate(V) + ! R <-- tmp x sigma^-0.5 + do j = 1, m + do i = 1, n + R(i,j) = tmp(i,j) * D(j) + enddo + enddo + + ! tmp <-- L x U + call dgemm( 'N', 'N', n, m, m, 1.d0 & + , L, size(L, 1), U, size(U, 1) & + , 0.d0, tmp, size(tmp, 1) ) + deallocate(U) + ! L <-- tmp x sigma^-0.5 + do j = 1, m + do i = 1, n + L(i,j) = tmp(i,j) * D(j) + enddo + enddo + + deallocate(D, tmp) + + ! --- + + allocate(S(m,m)) + call dgemm( 'T', 'N', m, m, n, 1.d0 & + , L, size(L, 1), R, size(R, 1) & + , 0.d0, S, size(S, 1) ) + + !print *, ' overlap aft SVD: ' + !do i = 1, m + ! write(*, '(1000(F16.10,X))') S(i,:) + !enddo + + deallocate(S) + + ! --- + +end subroutine impose_biorthog_svd + +! --- + +subroutine impose_weighted_biorthog_qr(m, n, thr_d, thr_nd, Vl, W, Vr) + + implicit none + integer, intent(in) :: m, n + double precision, intent(in) :: thr_d, thr_nd + double precision, intent(inout) :: Vl(m,n), W(m,m), Vr(m,n) + + integer :: i, j + integer :: LWORK, INFO + double precision :: accu_nd, accu_d + double precision, allocatable :: TAU(:), WORK(:) + double precision, allocatable :: S(:,:), R(:,:), tmp(:,:), Stmp(:,:) + + + call check_weighted_biorthog_binormalize(m, n, Vl, W, Vr, thr_d, thr_nd, .false.) + + ! --- + + allocate(Stmp(n,m), S(n,n)) + call dgemm( 'T', 'N', n, m, m, 1.d0 & + , Vl, size(Vl, 1), W, size(W, 1) & + , 0.d0, Stmp, size(Stmp, 1) ) + call dgemm( 'N', 'N', n, n, m, 1.d0 & + , Stmp, size(Stmp, 1), Vr, size(Vr, 1) & + , 0.d0, S, size(S, 1) ) + deallocate(Stmp) + + accu_nd = 0.d0 + accu_d = 0.d0 + do i = 1, n + do j = 1, n + if(i==j) then + accu_d += S(j,i) + else + accu_nd = accu_nd + S(j,i) * S(j,i) + endif + enddo + enddo + accu_nd = dsqrt(accu_nd) + + if((accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n))/dble(n) .lt. thr_d)) then + print *, ' bi-orthogonal vectors without QR !' + deallocate(S) + return + endif + + ! ------------------------------------------------------------------------------------- + ! QR factorization of S: S = Q x R + + + print *, ' apply QR decomposition ...' + + allocate( TAU(n), WORK(1) ) + + LWORK = -1 + call dgeqrf(n, n, S, n, TAU, WORK, LWORK, INFO) + if(INFO .ne. 0) then + print*,'dgeqrf failed !!', INFO + stop + endif + + LWORK = max(n, int(WORK(1))) + deallocate(WORK) + + allocate( WORK(LWORK) ) + call dgeqrf(n, n, S, n, TAU, WORK, LWORK, INFO) + if(INFO .ne. 0) then + print*,'dgeqrf failed !!', INFO + stop + endif + + ! save the upper triangular R + allocate( R(n,n) ) + R(:,:) = S(:,:) + + ! get Q + LWORK = -1 + call dorgqr(n, n, n, S, n, TAU, WORK, LWORK, INFO) + if(INFO .ne. 0) then + print*,'dorgqr failed !!', INFO + stop + endif + + LWORK = max(n, int(WORK(1))) + deallocate(WORK) + + allocate( WORK(LWORK) ) + call dorgqr(n, n, n, S, n, TAU, WORK, LWORK, INFO) + if(INFO .ne. 0) then + print*,'dorgqr failed !!', INFO + stop + endif + + deallocate( WORK, TAU ) + + ! + ! ------------------------------------------------------------------------------------- + + ! --- + + ! ------------------------------------------------------------------------------------- + ! get bi-orhtog left & right vectors: + ! Vr' = Vr x inv(R) + ! Vl' = inv(Q) x Vl = Q.T x Vl + + ! Q.T x Vl, where Q = S + + allocate( tmp(n,m) ) + call dgemm( 'T', 'T', n, m, n, 1.d0 & + , S, size(S, 1), Vl, size(Vl, 1) & + , 0.d0, tmp, size(tmp, 1) ) + + do i = 1, n + do j = 1, m + Vl(j,i) = tmp(i,j) + enddo + enddo + deallocate(tmp) + + ! --- + + ! inv(R) + !print *, ' inversing upper triangular matrix ...' + call dtrtri("U", "N", n, R, n, INFO) + if(INFO .ne. 0) then + print*,'dtrtri failed !!', INFO + stop + endif + !print *, ' inversing upper triangular matrix OK' + + do i = 1, n-1 + do j = i+1, n + R(j,i) = 0.d0 + enddo + enddo + + !print *, ' inv(R):' + !do i = 1, n + ! write(*, '(1000(F16.10,X))') R(i,:) + !enddo + + ! Vr x inv(R) + allocate( tmp(m,n) ) + call dgemm( 'N', 'N', m, n, n, 1.d0 & + , Vr, size(Vr, 1), R, size(R, 1) & + , 0.d0, tmp, size(tmp, 1) ) + deallocate( R ) + + do i = 1, n + do j = 1, m + Vr(j,i) = tmp(j,i) + enddo + enddo + deallocate(tmp) + + call check_weighted_biorthog_binormalize(m, n, Vl, W, Vr, thr_d, thr_nd, .false.) + + return +end subroutine impose_weighted_biorthog_qr + +! --- + +subroutine check_weighted_biorthog_binormalize(n, m, Vl, W, Vr, thr_d, thr_nd, stop_ifnot) + + implicit none + + integer, intent(in) :: n, m + logical, intent(in) :: stop_ifnot + double precision, intent(in) :: thr_d, thr_nd + double precision, intent(inout) :: Vl(n,m), W(n,n), Vr(n,m) + + integer :: i, j + double precision :: accu_d, accu_nd, s_tmp + double precision, allocatable :: S(:,:), Stmp(:,:) + + print *, ' check weighted bi-orthonormality' + + ! --- + + allocate(Stmp(m,n), S(m,m)) + call dgemm( 'T', 'N', m, n, n, 1.d0 & + , Vl, size(Vl, 1), W, size(W, 1) & + , 0.d0, Stmp, size(Stmp, 1) ) + call dgemm( 'N', 'N', m, m, n, 1.d0 & + , Stmp, size(Stmp, 1), Vr, size(Vr, 1) & + , 0.d0, S, size(S, 1) ) + deallocate(Stmp) + !print *, ' overlap matrix before:' + !do i = 1, m + ! write(*,'(1000(F16.10,X))') S(i,:) + !enddo + + ! S(i,i) = -1 + do i = 1, m + if( (S(i,i) + 1.d0) .lt. thr_d ) then + do j = 1, n + Vl(j,i) = -1.d0 * Vl(j,i) + enddo + S(i,i) = 1.d0 + endif + enddo + + accu_d = 0.d0 + accu_nd = 0.d0 + do i = 1, m + do j = 1, m + if(i==j) then + accu_d = accu_d + S(i,i) + else + accu_nd = accu_nd + S(j,i) * S(j,i) + endif + enddo + enddo + accu_nd = dsqrt(accu_nd) / dble(m) + print*, ' diag acc: ', accu_d + print*, ' nondiag acc: ', accu_nd + + ! --- + + if( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(m))/dble(m) .gt. thr_d) ) then + + do i = 1, m + print *, i, S(i,i) + if(dabs(S(i,i) - 1.d0) .gt. thr_d) then + s_tmp = 1.d0 / dsqrt(S(i,i)) + do j = 1, n + Vl(j,i) = Vl(j,i) * s_tmp + Vr(j,i) = Vr(j,i) * s_tmp + enddo + endif + enddo + + endif + + ! --- + + allocate(Stmp(m,n)) + call dgemm( 'T', 'N', m, n, n, 1.d0 & + , Vl, size(Vl, 1), W, size(W, 1) & + , 0.d0, Stmp, size(Stmp, 1) ) + call dgemm( 'N', 'N', m, m, n, 1.d0 & + , Stmp, size(Stmp, 1), Vr, size(Vr, 1) & + , 0.d0, S, size(S, 1) ) + deallocate(Stmp) + !print *, ' overlap matrix after:' + !do i = 1, m + ! write(*,'(1000(F16.10,X))') S(i,:) + !enddo + + accu_d = 0.d0 + accu_nd = 0.d0 + do i = 1, m + do j = 1, m + if(i==j) then + accu_d = accu_d + S(i,i) + else + accu_nd = accu_nd + S(j,i) * S(j,i) + endif + enddo + enddo + accu_nd = dsqrt(accu_nd) / dble(m) + print *, ' diag acc: ', accu_d + print *, ' nondiag acc: ', accu_nd + + deallocate(S) + + ! --- + + if( stop_ifnot .and. ((accu_nd .gt. thr_nd) .or. (dabs(accu_d-dble(m))/dble(m) .gt. thr_d)) ) then + print *, accu_nd, thr_nd + print *, dabs(accu_d-dble(m))/dble(m), thr_d + print *, ' weighted biorthog_binormalize failed !' + stop + endif + +end subroutine check_weighted_biorthog_binormalize + +! --- + +subroutine impose_weighted_biorthog_svd(n, m, overlap, L, R) + + implicit none + + integer, intent(in) :: n, m + double precision, intent(in) :: overlap(n,n) + double precision, intent(inout) :: L(n,m), R(n,m) + + integer :: i, j, num_linear_dependencies + double precision :: threshold + double precision, allocatable :: S(:,:), tmp(:,:),Stmp(:,:) + double precision, allocatable :: U(:,:), V(:,:), Vt(:,:), D(:) + + ! --- + + allocate(S(m,m),Stmp(n,m)) + + ! S = C.T x overlap x C + call dgemm( 'N', 'N', n, m, n, 1.d0 & + , overlap, size(overlap, 1), R, size(R, 1) & + , 0.d0, Stmp, size(Stmp, 1) ) + call dgemm( 'T', 'N', m, m, n, 1.d0 & + , L, size(L, 1), Stmp, size(Stmp, 1) & + , 0.d0, S, size(S, 1) ) + deallocate(Stmp) + + !print *, ' overlap bef SVD: ' + !do i = 1, m + ! write(*, '(1000(F25.16,X))') S(i,:) + !enddo + + ! --- + + allocate(U(m,m), Vt(m,m), D(m)) + + call svd(S, m, U, m, D, Vt, m, m, m) + + deallocate(S) + + threshold = 1.d-6 + num_linear_dependencies = 0 + do i = 1, m + if(abs(D(i)) <= threshold) then + D(i) = 0.d0 + num_linear_dependencies += 1 + else + ASSERT (D(i) > 0.d0) + D(i) = 1.d0 / dsqrt(D(i)) + endif + enddo + if(num_linear_dependencies > 0) then + write(*,*) ' linear dependencies = ', num_linear_dependencies + write(*,*) ' m = ', m + stop + endif + + allocate(V(m,m)) + do i = 1, m + do j = 1, m + V(j,i) = Vt(i,j) + enddo + enddo + deallocate(Vt) + + ! --- + + allocate(tmp(n,m)) + + ! tmp <-- R x V + call dgemm( 'N', 'N', n, m, m, 1.d0 & + , R, size(R, 1), V, size(V, 1) & + , 0.d0, tmp, size(tmp, 1) ) + deallocate(V) + ! R <-- tmp x sigma^-0.5 + do j = 1, m + do i = 1, n + R(i,j) = tmp(i,j) * D(j) + enddo + enddo + + ! tmp <-- L x U + call dgemm( 'N', 'N', n, m, m, 1.d0 & + , L, size(L, 1), U, size(U, 1) & + , 0.d0, tmp, size(tmp, 1) ) + deallocate(U) + ! L <-- tmp x sigma^-0.5 + do j = 1, m + do i = 1, n + L(i,j) = tmp(i,j) * D(j) + enddo + enddo + + deallocate(D, tmp) + + ! --- + + allocate(S(m,m),Stmp(n,m)) + ! S = C.T x overlap x C + call dgemm( 'N', 'N', n, m, n, 1.d0 & + , overlap, size(overlap, 1), R, size(R, 1) & + , 0.d0, Stmp, size(Stmp, 1) ) + call dgemm( 'T', 'N', m, m, n, 1.d0 & + , L, size(L, 1), Stmp, size(Stmp, 1) & + , 0.d0, S, size(S, 1) ) + deallocate(Stmp) + + !print *, ' overlap aft SVD with overlap: ' + !do i = 1, m + ! write(*, '(1000(F16.10,X))') S(i,:) + !enddo + + deallocate(S) + + return +end subroutine impose_weighted_biorthog_svd + +! --- + diff --git a/src/non_hermit_dav/new_routines.irp.f b/src/non_hermit_dav/new_routines.irp.f new file mode 100644 index 00000000..8db044d3 --- /dev/null +++ b/src/non_hermit_dav/new_routines.irp.f @@ -0,0 +1,670 @@ +subroutine non_hrmt_diag_split_degen_bi_orthog(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 + double precision :: shift_current + double precision :: r,thr,accu_d, accu_nd + integer, allocatable :: 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(:) + double precision :: accu,thr_cut, thr_norm=1d0 + + + thr_cut = 1.d-15 + print*,'Computing the left/right eigenvectors ...' + print*,'Using the degeneracy splitting algorithm' + ! initialization + shift_current = 1.d-15 + iteration = 0 + print*,'***** iteration = ',iteration + + + ! 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)) +! print*,'Aw' + do i = 1, n + iorder_origin(i) = i + diag_elem(i) = A(i,i) +! write(*,'(100(F16.10,X))')A(:,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 + + allocate(WR(n), WI(n), VL(n,n), VR(n,n), Aw(n,n)) + allocate(im_part(n),iorder(n)) + allocate( S(n,n) ) + + + Aw = A_save + call cancel_small_elmts(aw,n,thr_cut) + call lapack_diag_non_sym(n,Aw,WR,WI,VL,VR) + do i = 1, n + im_part(i) = -dabs(WI(i)) + iorder(i) = i + enddo + call dsort(im_part, iorder, n) + n_real_eigv = 0 + do i = 1, n + if(dabs(WI(i)).lt.1.d-20)then + n_real_eigv += 1 + else +! print*,'Found an imaginary component to eigenvalue' +! print*,'Re(i) + Im(i)',WR(i),WI(i) + endif + enddo + if(n_real_eigv.ne.n)then + shift_current = max(10.d0 * dabs(im_part(1)),shift_current*10.d0) + print*,'Largest imaginary part found in eigenvalues = ',im_part(1) + print*,'Splitting the degeneracies by ',shift_current + else + print*,'All eigenvalues are real !' + endif + + + do while(n_real_eigv.ne.n) + iteration += 1 + print*,'***** iteration = ',iteration + if(shift_current.gt.1.d-3)then + print*,'shift_current > 1.d-3 !!' + print*,'Your matrix intrinsically contains complex eigenvalues' + stop + endif + Aw = A_save + call cancel_small_elmts(Aw,n,thr_cut) + call split_matrix_degen(Aw,n,shift_current) + call lapack_diag_non_sym(n,Aw,WR,WI,VL,VR) + n_real_eigv = 0 + do i = 1, n + if(dabs(WI(i)).lt.1.d-20)then + n_real_eigv+= 1 + else +! print*,'Found an imaginary component to eigenvalue' +! print*,'Re(i) + Im(i)',WR(i),WI(i) + endif + enddo + if(n_real_eigv.ne.n)then + 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_current*10.d0) + print*,'Largest imaginary part found in eigenvalues = ',im_part(1) + print*,'Splitting the degeneracies by ',shift_current + else + print*,'All eigenvalues are real !' + endif + enddo + !!!!!!!!!!!!!!!! SORTING THE EIGENVALUES + do i = 1, n + eigval(i) = WR(i) + iorder(i) = i + enddo + call dsort(eigval,iorder,n) + do i = 1, n +! print*,'eigval(i) = ',eigval(i) + reigvec_tmp(:,i) = VR(:,iorder(i)) + leigvec_tmp(:,i) = Vl(:,iorder(i)) + enddo + +!!! ONCE ALL EIGENVALUES ARE REAL ::: CHECK BI-ORTHONORMALITY + ! check bi-orthogonality + call check_biorthog(n, n, leigvec_tmp, reigvec_tmp, accu_d, accu_nd, S, thresh_biorthog_diag, thresh_biorthog_nondiag, .false.) + print *, ' accu_nd bi-orthog = ', accu_nd + if(accu_nd .lt. thresh_biorthog_nondiag) then + print *, ' bi-orthogonality: ok' + else + print *, ' ' + print *, ' bi-orthogonality: not imposed yet' + print *, ' ' + print *, ' ' + print *, ' orthog between degen eigenvect' + print *, ' ' + double precision, allocatable :: S_nh_inv_half(:,:) + allocate(S_nh_inv_half(n,n)) + logical :: complex_root + deallocate(S_nh_inv_half) + call impose_orthog_degen_eigvec(n, eigval, reigvec_tmp) + call impose_orthog_degen_eigvec(n, eigval, leigvec_tmp) + call check_biorthog(n, n, leigvec_tmp, reigvec_tmp, accu_d, accu_nd, S, thresh_biorthog_diag, thresh_biorthog_nondiag, .false.) + if(accu_nd .lt. thresh_biorthog_nondiag) then + print *, ' bi-orthogonality: ok' + else + print*,'New vectors not bi-orthonormals at ',accu_nd + call impose_biorthog_qr(n, n, leigvec_tmp, reigvec_tmp, S) + call check_biorthog(n, n, leigvec_tmp, reigvec_tmp, accu_d, accu_nd, S, thresh_biorthog_diag, thresh_biorthog_nondiag, .false.) + if(accu_nd .lt. thresh_biorthog_nondiag) then + print *, ' bi-orthogonality: ok' + else + print*,'New vectors not bi-orthonormals at ',accu_nd + print*,'Must be a deep problem ...' + stop + endif + endif + endif + + !! EIGENVECTORS SORTED AND BI-ORTHONORMAL + do i = 1, n + do j = 1, n + VR(iorder_origin(j),i) = reigvec_tmp(j,i) + VL(iorder_origin(j),i) = leigvec_tmp(j,i) + enddo + enddo + + !! RECOMPUTING THE EIGENVALUES + eigval = 0.d0 + do i = 1, n + iorder(i) = i + accu = 0.d0 + do j = 1, n + accu += VL(j,i) * VR(j,i) + do k = 1, n + eigval(i) += VL(j,i) * A(j,k) * VR(k,i) + enddo + enddo + eigval(i) *= 1.d0/accu +! print*,'eigval(i) = ',eigval(i) + enddo + !! RESORT JUST TO BE SURE + call dsort(eigval, iorder, n) + do i = 1, n + do j = 1, n + reigvec(j,i) = VR(j,iorder(i)) + leigvec(j,i) = VL(j,iorder(i)) + enddo + enddo + print*,'Checking for final reigvec/leigvec' + shift_current = max(1.d-10,shift_current) + print*,'Thr for eigenvectors = ',shift_current + call check_EIGVEC(n, n, A, eigval, leigvec, reigvec, shift_current, thr_norm, .false.) + call check_biorthog(n, n, leigvec, reigvec, accu_d, accu_nd, S, thresh_biorthog_diag, thresh_biorthog_nondiag, .false.) + print *, ' accu_nd bi-orthog = ', accu_nd + + if(accu_nd .lt. thresh_biorthog_nondiag) then + print *, ' bi-orthogonality: ok' + else + print*,'Something went wrong in non_hrmt_diag_split_degen_bi_orthog' + print*,'Eigenvectors are not bi orthonormal ..' + print*,'accu_nd = ',accu_nd + stop + endif + +end + + + +subroutine non_hrmt_diag_split_degen_s_inv_half(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 + double precision :: shift_current + double precision :: r,thr,accu_d, accu_nd + integer, allocatable :: 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(:) + double precision :: accu,thr_cut, thr_norm=1.d0 + double precision, allocatable :: S_nh_inv_half(:,:) + logical :: complex_root + + + thr_cut = 1.d-15 + print*,'Computing the left/right eigenvectors ...' + print*,'Using the degeneracy splitting algorithm' + ! initialization + shift_current = 1.d-15 + iteration = 0 + print*,'***** iteration = ',iteration + + + ! 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)) +! print*,'Aw' + do i = 1, n + iorder_origin(i) = i + diag_elem(i) = A(i,i) +! write(*,'(100(F16.10,X))')A(:,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 + + allocate(WR(n), WI(n), VL(n,n), VR(n,n), Aw(n,n)) + allocate(im_part(n),iorder(n)) + allocate( S(n,n) ) + allocate(S_nh_inv_half(n,n)) + + + Aw = A_save + call cancel_small_elmts(aw,n,thr_cut) + call lapack_diag_non_sym(n,Aw,WR,WI,VL,VR) + do i = 1, n + im_part(i) = -dabs(WI(i)) + iorder(i) = i + enddo + call dsort(im_part, iorder, n) + n_real_eigv = 0 + do i = 1, n + if(dabs(WI(i)).lt.1.d-20)then + n_real_eigv += 1 + else +! print*,'Found an imaginary component to eigenvalue' +! print*,'Re(i) + Im(i)',WR(i),WI(i) + endif + enddo + if(n_real_eigv.ne.n)then + shift_current = max(10.d0 * dabs(im_part(1)),shift_current*10.d0) + print*,'Largest imaginary part found in eigenvalues = ',im_part(1) + print*,'Splitting the degeneracies by ',shift_current + else + print*,'All eigenvalues are real !' + endif + + + do while(n_real_eigv.ne.n) + iteration += 1 + print*,'***** iteration = ',iteration + if(shift_current.gt.1.d-3)then + print*,'shift_current > 1.d-3 !!' + print*,'Your matrix intrinsically contains complex eigenvalues' + stop + endif + Aw = A_save +! thr_cut = shift_current + call cancel_small_elmts(Aw,n,thr_cut) + call split_matrix_degen(Aw,n,shift_current) + call lapack_diag_non_sym(n,Aw,WR,WI,VL,VR) + n_real_eigv = 0 + do i = 1, n + if(dabs(WI(i)).lt.1.d-20)then + n_real_eigv+= 1 + else +! print*,'Found an imaginary component to eigenvalue' +! print*,'Re(i) + Im(i)',WR(i),WI(i) + endif + enddo + if(n_real_eigv.ne.n)then + 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_current*10.d0) + print*,'Largest imaginary part found in eigenvalues = ',im_part(1) + print*,'Splitting the degeneracies by ',shift_current + else + print*,'All eigenvalues are real !' + endif + enddo + !!!!!!!!!!!!!!!! SORTING THE EIGENVALUES + do i = 1, n + eigval(i) = WR(i) + iorder(i) = i + enddo + call dsort(eigval,iorder,n) + do i = 1, n +! print*,'eigval(i) = ',eigval(i) + reigvec_tmp(:,i) = VR(:,iorder(i)) + leigvec_tmp(:,i) = Vl(:,iorder(i)) + enddo + +!!! ONCE ALL EIGENVALUES ARE REAL ::: CHECK BI-ORTHONORMALITY + ! check bi-orthogonality + call check_biorthog(n, n, leigvec_tmp, reigvec_tmp, accu_d, accu_nd, S, thresh_biorthog_diag, thresh_biorthog_nondiag, .false.) + print *, ' accu_nd bi-orthog = ', accu_nd + if(accu_nd .lt. thresh_biorthog_nondiag) then + print *, ' bi-orthogonality: ok' + else + print *, ' ' + print *, ' bi-orthogonality: not imposed yet' + if(complex_root) then + print *, ' ' + print *, ' ' + print *, ' orthog between degen eigenvect' + print *, ' ' + ! bi-orthonormalization using orthogonalization of left, right and then QR between left and right + call impose_orthog_degen_eigvec(n, eigval, reigvec_tmp) ! orthogonalization of reigvec + call impose_orthog_degen_eigvec(n, eigval, leigvec_tmp) ! orthogonalization of leigvec + call check_biorthog(n, n, leigvec_tmp, reigvec_tmp, accu_d, accu_nd, S, thresh_biorthog_diag, thresh_biorthog_nondiag, .false.) + + if(accu_nd .lt. thresh_biorthog_nondiag) then + print *, ' bi-orthogonality: ok' + else + print*,'New vectors not bi-orthonormals at ', accu_nd + call get_inv_half_nonsymmat_diago(S, n, S_nh_inv_half, complex_root) + if(complex_root)then + call impose_biorthog_qr(n, n, leigvec_tmp, reigvec_tmp, S) ! bi-orthonormalization using QR + else + print*,'S^{-1/2} exists !!' + call bi_ortho_s_inv_half(n,leigvec_tmp,reigvec_tmp,S_nh_inv_half) ! use of S^{-1/2} bi-orthonormalization + endif + endif + else ! the matrix S^{-1/2} exists + print*,'S^{-1/2} exists !!' + call bi_ortho_s_inv_half(n,leigvec_tmp,reigvec_tmp,S_nh_inv_half) ! use of S^{-1/2} bi-orthonormalization + endif + call check_biorthog(n, n, leigvec_tmp, reigvec_tmp, accu_d, accu_nd, S, thresh_biorthog_diag, thresh_biorthog_nondiag, .false.) + if(accu_nd .lt. thresh_biorthog_nondiag) then + print *, ' bi-orthogonality: ok' + else + print*,'New vectors not bi-orthonormals at ',accu_nd + print*,'Must be a deep problem ...' + stop + endif + endif + + !! EIGENVECTORS SORTED AND BI-ORTHONORMAL + do i = 1, n + do j = 1, n + VR(iorder_origin(j),i) = reigvec_tmp(j,i) + VL(iorder_origin(j),i) = leigvec_tmp(j,i) + enddo + enddo + + !! RECOMPUTING THE EIGENVALUES + eigval = 0.d0 + do i = 1, n + iorder(i) = i + accu = 0.d0 + do j = 1, n + accu += VL(j,i) * VR(j,i) + do k = 1, n + eigval(i) += VL(j,i) * A(j,k) * VR(k,i) + enddo + enddo + eigval(i) *= 1.d0/accu +! print*,'eigval(i) = ',eigval(i) + enddo + !! RESORT JUST TO BE SURE + call dsort(eigval, iorder, n) + do i = 1, n + do j = 1, n + reigvec(j,i) = VR(j,iorder(i)) + leigvec(j,i) = VL(j,iorder(i)) + enddo + enddo + print*,'Checking for final reigvec/leigvec' + shift_current = max(1.d-10,shift_current) + print*,'Thr for eigenvectors = ',shift_current + call check_EIGVEC(n, n, A, eigval, leigvec, reigvec, shift_current, thr_norm, .false.) + call check_biorthog(n, n, leigvec, reigvec, accu_d, accu_nd, S, thresh_biorthog_diag, thresh_biorthog_nondiag, .false.) + print *, ' accu_nd bi-orthog = ', accu_nd + + if(accu_nd .lt. thresh_biorthog_nondiag) then + print *, ' bi-orthogonality: ok' + else + print*,'Something went wrong in non_hrmt_diag_split_degen_bi_orthog' + print*,'Eigenvectors are not bi orthonormal ..' + print*,'accu_nd = ',accu_nd + stop + endif + +end + + +subroutine non_hrmt_fock_mat(n, A, leigvec, reigvec, n_real_eigv, eigval) + + BEGIN_DOC + ! + ! routine returning the eigenvalues and left/right eigenvectors of the TC fock matrix + ! + 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 + double precision :: shift_current + double precision :: r,thr,accu_d, accu_nd + integer, allocatable :: 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(:) + double precision :: accu,thr_cut + double precision, allocatable :: S_nh_inv_half(:,:) + logical :: complex_root + double precision :: thr_norm=1d0 + + + thr_cut = 1.d-15 + print*,'Computing the left/right eigenvectors ...' + print*,'Using the degeneracy splitting algorithm' + ! initialization + shift_current = 1.d-15 + iteration = 0 + print*,'***** iteration = ',iteration + + + ! 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)) +! print*,'Aw' + do i = 1, n + iorder_origin(i) = i + diag_elem(i) = A(i,i) +! write(*,'(100(F16.10,X))')A(:,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 + + allocate(WR(n), WI(n), VL(n,n), VR(n,n), Aw(n,n)) + allocate(im_part(n),iorder(n)) + allocate( S(n,n) ) + allocate(S_nh_inv_half(n,n)) + + + Aw = A_save + call cancel_small_elmts(aw,n,thr_cut) + call lapack_diag_non_sym(n,Aw,WR,WI,VL,VR) + do i = 1, n + im_part(i) = -dabs(WI(i)) + iorder(i) = i + enddo + call dsort(im_part, iorder, n) + n_real_eigv = 0 + do i = 1, n + if(dabs(WI(i)).lt.1.d-20)then + n_real_eigv += 1 + else +! print*,'Found an imaginary component to eigenvalue' +! print*,'Re(i) + Im(i)',WR(i),WI(i) + endif + enddo + if(n_real_eigv.ne.n)then + shift_current = max(10.d0 * dabs(im_part(1)),shift_current*10.d0) + print*,'Largest imaginary part found in eigenvalues = ',im_part(1) + print*,'Splitting the degeneracies by ',shift_current + else + print*,'All eigenvalues are real !' + endif + + + do while(n_real_eigv.ne.n) + iteration += 1 + print*,'***** iteration = ',iteration + if(shift_current.gt.1.d-3)then + print*,'shift_current > 1.d-3 !!' + print*,'Your matrix intrinsically contains complex eigenvalues' + stop + endif + Aw = A_save +! thr_cut = shift_current + call cancel_small_elmts(Aw,n,thr_cut) + call split_matrix_degen(Aw,n,shift_current) + call lapack_diag_non_sym(n,Aw,WR,WI,VL,VR) + n_real_eigv = 0 + do i = 1, n + if(dabs(WI(i)).lt.1.d-20)then + n_real_eigv+= 1 + else +! print*,'Found an imaginary component to eigenvalue' +! print*,'Re(i) + Im(i)',WR(i),WI(i) + endif + enddo + if(n_real_eigv.ne.n)then + 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_current*10.d0) + print*,'Largest imaginary part found in eigenvalues = ',im_part(1) + print*,'Splitting the degeneracies by ',shift_current + else + print*,'All eigenvalues are real !' + endif + enddo + !!!!!!!!!!!!!!!! SORTING THE EIGENVALUES + do i = 1, n + eigval(i) = WR(i) + iorder(i) = i + enddo + call dsort(eigval,iorder,n) + do i = 1, n +! print*,'eigval(i) = ',eigval(i) + reigvec_tmp(:,i) = VR(:,iorder(i)) + leigvec_tmp(:,i) = Vl(:,iorder(i)) + enddo + +!!! ONCE ALL EIGENVALUES ARE REAL ::: CHECK BI-ORTHONORMALITY + ! check bi-orthogonality + call check_biorthog(n, n, leigvec_tmp, reigvec_tmp, accu_d, accu_nd, S, thresh_biorthog_diag, thresh_biorthog_nondiag, .false.) + print *, ' accu_nd bi-orthog = ', accu_nd + if(accu_nd .lt. thresh_biorthog_nondiag) then + print *, ' bi-orthogonality: ok' + else + print *, ' ' + print *, ' bi-orthogonality: not imposed yet' + print *, ' ' + print *, ' ' + print *, ' Using impose_unique_biorthog_degen_eigvec' + print *, ' ' + ! bi-orthonormalization using orthogonalization of left, right and then QR between left and right + call impose_unique_biorthog_degen_eigvec(n, eigval, mo_coef, leigvec_tmp, reigvec_tmp) + call check_biorthog(n, n, leigvec_tmp, reigvec_tmp, accu_d, accu_nd, S, thresh_biorthog_diag, thresh_biorthog_nondiag, .false.) + print*,'accu_nd = ',accu_nd + if(accu_nd .lt. thresh_biorthog_nondiag) then + print *, ' bi-orthogonality: ok' + else + print*,'New vectors not bi-orthonormals at ',accu_nd + call get_inv_half_nonsymmat_diago(S, n, S_nh_inv_half,complex_root) + if(complex_root)then + print*,'S^{-1/2} does not exits, using QR bi-orthogonalization' + call impose_biorthog_qr(n, n, leigvec_tmp, reigvec_tmp, S) ! bi-orthonormalization using QR + else + print*,'S^{-1/2} exists !!' + call bi_ortho_s_inv_half(n,leigvec_tmp,reigvec_tmp,S_nh_inv_half) ! use of S^{-1/2} bi-orthonormalization + endif + endif + call check_biorthog(n, n, leigvec_tmp, reigvec_tmp, accu_d, accu_nd, S, thresh_biorthog_diag, thresh_biorthog_nondiag, .false.) + if(accu_nd .lt. thresh_biorthog_nondiag) then + print *, ' bi-orthogonality: ok' + else + print*,'New vectors not bi-orthonormals at ',accu_nd + print*,'Must be a deep problem ...' + stop + endif + endif + + !! EIGENVECTORS SORTED AND BI-ORTHONORMAL + do i = 1, n + do j = 1, n + VR(iorder_origin(j),i) = reigvec_tmp(j,i) + VL(iorder_origin(j),i) = leigvec_tmp(j,i) + enddo + enddo + + !! RECOMPUTING THE EIGENVALUES + eigval = 0.d0 + do i = 1, n + iorder(i) = i + accu = 0.d0 + do j = 1, n + accu += VL(j,i) * VR(j,i) + do k = 1, n + eigval(i) += VL(j,i) * A(j,k) * VR(k,i) + enddo + enddo + eigval(i) *= 1.d0/accu +! print*,'eigval(i) = ',eigval(i) + enddo + !! RESORT JUST TO BE SURE + call dsort(eigval, iorder, n) + do i = 1, n + do j = 1, n + reigvec(j,i) = VR(j,iorder(i)) + leigvec(j,i) = VL(j,iorder(i)) + enddo + enddo + print*,'Checking for final reigvec/leigvec' + shift_current = max(1.d-10,shift_current) + print*,'Thr for eigenvectors = ',shift_current + call check_EIGVEC(n, n, A, eigval, leigvec, reigvec, shift_current, thr_norm, .false.) + call check_biorthog(n, n, leigvec, reigvec, accu_d, accu_nd, S, thresh_biorthog_diag, thresh_biorthog_nondiag, .false.) + print *, ' accu_nd bi-orthog = ', accu_nd + + if(accu_nd .lt. thresh_biorthog_nondiag) then + print *, ' bi-orthogonality: ok' + else + print*,'Something went wrong in non_hrmt_diag_split_degen_bi_orthog' + print*,'Eigenvectors are not bi orthonormal ..' + print*,'accu_nd = ',accu_nd + stop + endif + +end + + diff --git a/src/non_hermit_dav/project.irp.f b/src/non_hermit_dav/project.irp.f new file mode 100644 index 00000000..c04719ac --- /dev/null +++ b/src/non_hermit_dav/project.irp.f @@ -0,0 +1,53 @@ +subroutine h_non_hermite(v,u,Hmat,a,N_st,sze) + implicit none + BEGIN_DOC + ! Template of routine for the application of H + ! + ! Here, it is done with the Hamiltonian matrix + ! + ! on the set of determinants of psi_det + ! + ! Computes $v = a * H | u \rangle$ + ! + END_DOC + integer, intent(in) :: N_st,sze + double precision, intent(in) :: u(sze,N_st), Hmat(sze,sze), a + double precision, intent(inout) :: v(sze,N_st) + integer :: i,j,k + do k = 1, N_st + do j = 1, sze + do i = 1, sze + v(i,k) += a * u(j,k) * Hmat(i,j) + enddo + enddo + enddo +end + + +subroutine exp_tau_H(u,v,hmat,tau,et,N_st,sze) + implicit none + BEGIN_DOC +! realises v = (1 - tau (H - et)) u + END_DOC + integer, intent(in) :: N_st,sze + double precision, intent(in) :: hmat(sze,sze), u(sze,N_st), tau, et + double precision, intent(out):: v(sze,N_st) + double precision :: a + integer :: i,j + v = (1.d0 + tau * et) * u + a = -1.d0 * tau + call h_non_hermite(v,u,Hmat,a,N_st,sze) +end + +double precision function project_phi0(u,Hmat0,N_st,sze) + implicit none + integer, intent(in) :: N_st,sze + double precision, intent(in) :: u(sze,N_st), Hmat0(sze) + integer :: j + project_phi0 = 0.d0 + do j = 1, sze + project_phi0 += u(j,1) * Hmat0(j) + enddo + project_phi0 *= 1.d0 / u(1,1) +end + diff --git a/src/non_hermit_dav/utils.irp.f b/src/non_hermit_dav/utils.irp.f new file mode 100644 index 00000000..7f331a6b --- /dev/null +++ b/src/non_hermit_dav/utils.irp.f @@ -0,0 +1,325 @@ + +subroutine get_inv_half_svd(matrix, n, matrix_inv_half) + + BEGIN_DOC + ! :math:`X = S^{-1/2}` obtained by SVD + END_DOC + + implicit none + + integer, intent(in) :: n + double precision, intent(in) :: matrix(n,n) + double precision, intent(out) :: matrix_inv_half(n,n) + + integer :: num_linear_dependencies + integer :: LDA, LDC + integer :: info, i, j, k + double precision, parameter :: threshold = 1.d-6 + double precision, allocatable :: U(:,:),Vt(:,:), D(:),matrix_half(:,:),D_half(:) + + double precision :: accu_d,accu_nd + + LDA = size(matrix, 1) + LDC = size(matrix_inv_half, 1) + if(LDA .ne. LDC) then + print*, ' LDA != LDC' + stop + endif + + print*, ' n = ', n + print*, ' LDA = ', LDA + print*, ' LDC = ', LDC + + double precision,allocatable :: WR(:),WI(:),VL(:,:),VR(:,:) + allocate(WR(n),WI(n),VL(n,n),VR(n,n)) + call lapack_diag_non_sym(n,matrix,WR,WI,VL,VR) + do i = 1, n + print*,'WR,WI',WR(i),WI(i) + enddo + + + allocate(U(LDC,n), Vt(LDA,n), D(n)) + + call svd(matrix, LDA, U, LDC, D, Vt, LDA, n, n) + double precision, allocatable :: tmp1(:,:),tmp2(:,:),D_mat(:,:) + allocate(tmp1(n,n),tmp2(n,n),D_mat(n,n),matrix_half(n,n),D_half(n)) + D_mat = 0.d0 + do i = 1,n + D_mat(i,i) = D(i) + enddo + ! matrix = U D Vt + ! tmp1 = U D + tmp1 = 0.d0 + call dgemm( 'N', 'N', n, n, n, 1.d0 & + , U, size(U, 1), D_mat, size(D_mat, 1) & + , 0.d0, tmp1, size(tmp1, 1) ) + ! tmp2 = tmp1 X Vt = matrix + tmp2 = 0.d0 + call dgemm( 'N', 'N', n, n, n, 1.d0 & + , tmp1, size(tmp1, 1), Vt, size(Vt, 1) & + , 0.d0, tmp2, size(tmp2, 1) ) + print*,'Checking the recomposition of the matrix' + accu_nd = 0.d0 + accu_d = 0.d0 + do i = 1, n + accu_d += dabs(tmp2(i,i) - matrix(i,i)) + do j = 1, n + if(i==j)cycle + accu_nd += dabs(tmp2(j,i) - matrix(j,i)) + enddo + enddo + print*,'accu_d =',accu_d + print*,'accu_nd =',accu_nd + print*,'passed the recomposition' + + num_linear_dependencies = 0 + do i = 1, n + if(abs(D(i)) <= threshold) then + D(i) = 0.d0 + num_linear_dependencies += 1 + else + ASSERT (D(i) > 0.d0) + D_half(i) = dsqrt(D(i)) + D(i) = 1.d0 / dsqrt(D(i)) + endif + enddo + write(*,*) ' linear dependencies', num_linear_dependencies + + matrix_inv_half = 0.d0 + matrix_half = 0.d0 + do k = 1, n + if(D(k) /= 0.d0) then + do j = 1, n + do i = 1, n +! matrix_inv_half(i,j) = matrix_inv_half(i,j) + U(i,k) * D(k) * Vt(k,j) + matrix_inv_half(i,j) = matrix_inv_half(i,j) + U(i,k) * D(k) * Vt(j,k) + matrix_half(i,j) = matrix_half(i,j) + U(i,k) * D_half(k) * Vt(j,k) + enddo + enddo + endif + enddo + print*,'testing S^1/2 * S^1/2= S' + ! tmp1 = S^1/2 X S^1/2 + tmp1 = 0.d0 + call dgemm( 'N', 'N', n, n, n, 1.d0 & + , matrix_half, size(matrix_half, 1), matrix_half, size(matrix_half, 1) & + , 0.d0, tmp1, size(tmp1, 1) ) + accu_nd = 0.d0 + accu_d = 0.d0 + do i = 1, n + accu_d += dabs(tmp1(i,i) - matrix(i,i)) + do j = 1, n + if(i==j)cycle + accu_nd += dabs(tmp1(j,i) - matrix(j,i)) + enddo + enddo + print*,'accu_d =',accu_d + print*,'accu_nd =',accu_nd + +! print*,'S inv half' +! do i = 1, n +! write(*, '(1000(F16.10,X))') matrix_inv_half(i,:) +! enddo + + double precision, allocatable :: pseudo_inverse(:,:),identity(:,:) + allocate( pseudo_inverse(n,n),identity(n,n)) + call get_pseudo_inverse(matrix,n,n,n,pseudo_inverse,n,threshold) + + ! S^-1 X S = 1 +! identity = 0.d0 +! call dgemm( 'N', 'N', n, n, n, 1.d0 & +! , matrix, size(matrix, 1), pseudo_inverse, size(pseudo_inverse, 1) & +! , 0.d0, identity, size(identity, 1) ) + print*,'Checking S^-1/2 X S^-1/2 = S^-1 ?' + ! S^-1/2 X S^-1/2 = S^-1 ? + tmp1 = 0.d0 + call dgemm( 'N', 'N', n, n, n, 1.d0 & + ,matrix_inv_half, size(matrix_inv_half, 1), matrix_inv_half, size(matrix_inv_half, 1) & + , 0.d0, tmp1, size(tmp1, 1) ) + accu_nd = 0.d0 + accu_d = 0.d0 + do i = 1, n + accu_d += dabs(1.d0 - pseudo_inverse(i,i)) + do j = 1, n + if(i==j)cycle + accu_nd += dabs(tmp1(j,i) - pseudo_inverse(j,i)) + enddo + enddo + print*,'accu_d =',accu_d + print*,'accu_nd =',accu_nd + + stop +! +! ! ( S^-1/2 x S ) x S^-1/2 +! Stmp2 = 0.d0 +! call dgemm( 'N', 'N', n, n, n, 1.d0 & +! , Stmp, size(Stmp, 1), matrix_inv_half, size(matrix_inv_half, 1) & +! , 0.d0, Stmp2, size(Stmp2, 1) ) + + ! S^-1/2 x ( S^-1/2 x S ) +! Stmp2 = 0.d0 +! call dgemm( 'N', 'N', n, n, n, 1.d0 & +! , matrix_inv_half, size(matrix_inv_half, 1), Stmp, size(Stmp, 1) & +! , 0.d0, Stmp2, size(Stmp2, 1) ) + +! do i = 1, n +! do j = 1, n +! if(i==j) then +! accu_d += Stmp2(j,i) +! else +! accu_nd = accu_nd + Stmp2(j,i) * Stmp2(j,i) +! endif +! enddo +! enddo +! accu_nd = dsqrt(accu_nd) +! print*, ' after S^-1/2: sum of off-diag S elements = ', accu_nd +! print*, ' after S^-1/2: sum of diag S elements = ', accu_d +! do i = 1, n +! write(*,'(1000(F16.10,X))') Stmp2(i,:) +! enddo + + !double precision :: thresh + !thresh = 1.d-10 + !if( accu_nd.gt.thresh .or. dabs(accu_d-dble(n)).gt.thresh) then + ! stop + !endif + +end subroutine get_inv_half_svd + +! --- + +subroutine get_inv_half_nonsymmat_diago(matrix, n, matrix_inv_half, complex_root) + + BEGIN_DOC + ! input: S = matrix + ! output: S^{-1/2} = matrix_inv_half obtained by diagonalization + ! + ! S = VR D VL^T + ! = VR D^{1/2} D^{1/2} VL^T + ! = VR D^{1/2} VL^T VR D^{1/2} VL^T + ! = S^{1/2} S^{1/2} with S = VR D^{1/2} VL^T + ! + ! == > S^{-1/2} = VR D^{-1/2} VL^T + ! + END_DOC + + implicit none + + integer, intent(in) :: n + double precision, intent(in) :: matrix(n,n) + logical, intent(out) :: complex_root + double precision, intent(out) :: matrix_inv_half(n,n) + + integer :: i, j + double precision :: accu_d, accu_nd + double precision, allocatable :: WR(:), WI(:), VL(:,:), VR(:,:), S(:,:), S_diag(:) + double precision, allocatable :: tmp1(:,:), D_mat(:,:) + + complex_root = .False. + + matrix_inv_half = 0.D0 + print*,'Computing S^{-1/2}' + + allocate(WR(n), WI(n), VL(n,n), VR(n,n)) + call lapack_diag_non_sym(n, matrix, WR, WI, VL, VR) + + allocate(S(n,n)) + call check_biorthog(n, n, VL, VR, accu_d, accu_nd, S) + print*,'accu_nd S^{-1/2}',accu_nd + if(accu_nd.gt.1.d-10) then + complex_root = .True. ! if vectors are not bi-orthogonal return + print*,'Eigenvectors of S are not bi-orthonormal, skipping S^{-1/2}' + return + endif + + allocate(S_diag(n)) + do i = 1, n + S_diag(i) = 1.d0/dsqrt(S(i,i)) + if(dabs(WI(i)).gt.1.d-20.or.WR(i).lt.0.d0)then ! check that eigenvalues are real and positive + complex_root = .True. + print*,'Eigenvalues of S have imaginary part ' + print*,'WR(i),WI(i)',WR(i), WR(i) + print*,'Skipping S^{-1/2}' + return + endif + enddo + deallocate(S) + + if(complex_root) return + + ! normalization of vectors + do i = 1, n + if(S_diag(i).eq.1.d0) cycle + do j = 1,n + VL(j,i) *= S_diag(i) + VR(j,i) *= S_diag(i) + enddo + enddo + deallocate(S_diag) + + allocate(tmp1(n,n), D_mat(n,n)) + + D_mat = 0.d0 + do i = 1, n + D_mat(i,i) = 1.d0/dsqrt(WR(i)) + enddo + deallocate(WR, WI) + + ! tmp1 = VR D^{-1/2} + tmp1 = 0.d0 + call dgemm( 'N', 'N', n, n, n, 1.d0 & + , VR, size(VR, 1), D_mat, size(D_mat, 1) & + , 0.d0, tmp1, size(tmp1, 1) ) + deallocate(VR, D_mat) + + ! S^{-1/2} = tmp1 X VL^T + matrix_inv_half = 0.d0 + call dgemm( 'N', 'T', n, n, n, 1.d0 & + , tmp1, size(tmp1, 1), VL, size(VL, 1) & + , 0.d0, matrix_inv_half, size(matrix_inv_half, 1) ) + deallocate(tmp1, VL) + +end + +! --- + +subroutine bi_ortho_s_inv_half(n,leigvec,reigvec,S_nh_inv_half) + implicit none + integer, intent(in) :: n + double precision, intent(in) :: S_nh_inv_half(n,n) + double precision, intent(inout) :: leigvec(n,n),reigvec(n,n) + BEGIN_DOC + ! bi-orthonormalization of left and right vectors + ! + ! S = VL^T VR + ! + ! S^{-1/2} S S^{-1/2} = 1 = S^{-1/2} VL^T VR S^{-1/2} = VL_new^T VR_new + ! + ! VL_new = VL (S^{-1/2})^T + ! + ! VR_new = VR S^{^{-1/2}} + END_DOC + double precision,allocatable :: vl_tmp(:,:),vr_tmp(:,:) + print*,'Bi-orthonormalization using S^{-1/2}' + allocate(vl_tmp(n,n),vr_tmp(n,n)) + vl_tmp = leigvec + vr_tmp = reigvec + ! VL_new = VL (S^{-1/2})^T + call dgemm( 'N', 'T', n, n, n, 1.d0 & + , vl_tmp, size(vl_tmp, 1), S_nh_inv_half, size(S_nh_inv_half, 1) & + , 0.d0, leigvec, size(leigvec, 1) ) + ! VR_new = VR S^{^{-1/2}} + call dgemm( 'N', 'N', n, n, n, 1.d0 & + , vr_tmp, size(vr_tmp, 1), S_nh_inv_half, size(S_nh_inv_half, 1) & + , 0.d0, reigvec, size(reigvec, 1) ) + double precision :: accu_d, accu_nd + double precision,allocatable :: S(:,:) + allocate(S(n,n)) + call check_biorthog(n, n, leigvec, reigvec, accu_d, accu_nd, S) + if(dabs(accu_d - n).gt.1.d-10 .or. accu_nd .gt.1.d-8 )then + print*,'Pb in bi_ortho_s_inv_half !!' + print*,'accu_d =',accu_d + print*,'accu_nd =',accu_nd + stop + endif +end