mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-11-08 22:43:38 +01:00
2955 lines
67 KiB
Fortran
2955 lines
67 KiB
Fortran
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) = <i|Psi_left(j)> :: projection on the basis element |i> on the jth left eigenvector
|
|
!
|
|
! VR(i,j) = <i|Psi_right(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 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) = <i|Psi_left(j)> :: projection on the basis element |i> on the jth left eigenvector
|
|
!
|
|
! VR(i,j) = <i|Psi_right(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 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 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) = <i|Psi_left(j)> :: projection on the basis element |i> on the jth left eigenvector
|
|
!
|
|
! VR(i,j) = <i|Psi_right(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 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 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_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 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_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 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_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_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_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_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_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 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 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(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 S s'il y a besoin
|
|
!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))
|
|
!print*, i, S(i,i)
|
|
else
|
|
accu_nd = accu_nd + S(j,i) * S(j,i)
|
|
endif
|
|
enddo
|
|
enddo
|
|
!accu_nd = dsqrt(accu_nd) / dble(m*m)
|
|
accu_nd = dsqrt(accu_nd) / dble(m)
|
|
|
|
if((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)
|
|
else
|
|
print *, ' vectors are bi-orthogonals'
|
|
endif
|
|
|
|
! ---
|
|
|
|
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_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 reorder_degen_eigvec(n, deg_num, e0, L0, R0)
|
|
|
|
implicit none
|
|
|
|
integer, intent(in) :: n
|
|
double precision, intent(inout) :: e0(n), L0(n,n), R0(n,n)
|
|
integer, intent(out) :: deg_num(n)
|
|
|
|
logical :: complex_root
|
|
integer :: i, j, k, m, ii, j_tmp
|
|
double precision :: ei, ej, de, de_thr
|
|
double precision :: accu_d, accu_nd
|
|
double precision :: e0_tmp, L0_tmp(n), R0_tmp(n)
|
|
double precision, allocatable :: L(:,:), R(:,:), S(:,:), S_inv_half(:,:)
|
|
|
|
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
|
|
|
|
ii = 0
|
|
do j = i+1, n
|
|
ej = e0(j)
|
|
de = dabs(ei - ej)
|
|
|
|
if(de .lt. de_thr) then
|
|
ii = ii + 1
|
|
|
|
j_tmp = i + ii
|
|
|
|
deg_num(j_tmp) = 0
|
|
|
|
e0_tmp = e0(j_tmp)
|
|
e0(j_tmp) = e0(j)
|
|
e0(j) = e0_tmp
|
|
|
|
L0_tmp(1:n) = L0(1:n,j_tmp)
|
|
L0(1:n,j_tmp) = L0(1:n,j)
|
|
L0(1:n,j) = L0_tmp(1:n)
|
|
|
|
R0_tmp(1:n) = R0(1:n,j_tmp)
|
|
R0(1:n,j_tmp) = R0(1:n,j)
|
|
R0(1:n,j) = R0_tmp(1:n)
|
|
endif
|
|
enddo
|
|
|
|
deg_num(i) = ii + 1
|
|
enddo
|
|
|
|
ii = 0
|
|
do i = 1, n
|
|
if(deg_num(i) .gt. 1) then
|
|
!print *, ' degen on', i, deg_num(i), e0(i)
|
|
ii = ii + 1
|
|
endif
|
|
enddo
|
|
|
|
if(ii .eq. 0) then
|
|
print*, ' WARNING: bi-orthogonality is lost but there is no degeneracies'
|
|
print*, ' rotations may change energy'
|
|
stop
|
|
endif
|
|
|
|
print *, ii, ' type of degeneracies'
|
|
|
|
! ---
|
|
|
|
! do i = 1, n
|
|
! m = deg_num(i)
|
|
!
|
|
! if(m .gt. 1) then
|
|
!
|
|
! allocate(L(n,m))
|
|
! allocate(R(n,m),S(m,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 dgemm( 'T', 'N', m, m, n, 1.d0 &
|
|
! ! , L, size(L, 1), R, size(R, 1) &
|
|
! ! , 0.d0, S, size(S, 1) )
|
|
! !print*, 'Overlap matrix '
|
|
! !accu_nd = 0.d0
|
|
! !do j = 1, m
|
|
! ! write(*,'(100(F16.10,X))') S(1:m,j)
|
|
! ! do k = 1, m
|
|
! ! if(j==k) cycle
|
|
! ! accu_nd += dabs(S(j,k))
|
|
! ! enddo
|
|
! !enddo
|
|
! !print*,'accu_nd = ',accu_nd
|
|
!! if(accu_nd .gt.1.d-10) then
|
|
!! stop
|
|
!! endif
|
|
!
|
|
! 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, S)
|
|
!
|
|
! endif
|
|
! enddo
|
|
!
|
|
end
|
|
|
|
! ---
|
|
|
|
subroutine impose_biorthog_degen_eigvec(n, deg_num, e0, L0, R0)
|
|
|
|
implicit none
|
|
|
|
integer, intent(in) :: n, deg_num(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
|
|
double precision, allocatable :: L(:,:), R(:,:), S(:,:), S_inv_half(:,:)
|
|
|
|
!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), R(n,m), S(m,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
|
|
|
|
! ---
|
|
|
|
!print*, 'Overlap matrix before'
|
|
call dgemm( 'T', 'N', m, m, n, 1.d0 &
|
|
, L, size(L, 1), R, size(R, 1) &
|
|
, 0.d0, S, size(S, 1) )
|
|
|
|
accu_nd = 0.d0
|
|
do j = 1, m
|
|
!write(*,'(100(F16.10,X))') S(1:m,j)
|
|
do k = 1, m
|
|
if(j==k) cycle
|
|
accu_nd += dabs(S(j,k))
|
|
enddo
|
|
enddo
|
|
|
|
if(accu_nd .lt. 1d-12) then
|
|
deallocate(S, L, R)
|
|
cycle
|
|
endif
|
|
|
|
!print*, ' accu_nd before = ', accu_nd
|
|
|
|
call impose_biorthog_svd(n, m, L, R)
|
|
|
|
!print*, 'Overlap matrix after'
|
|
call dgemm( 'T', 'N', m, m, n, 1.d0 &
|
|
, L, size(L, 1), R, size(R, 1) &
|
|
, 0.d0, S, size(S, 1) )
|
|
accu_nd = 0.d0
|
|
do j = 1, m
|
|
!write(*,'(100(F16.10,X))') S(1:m,j)
|
|
do k = 1, m
|
|
if(j==k) cycle
|
|
accu_nd += dabs(S(j,k))
|
|
enddo
|
|
enddo
|
|
!print*,' accu_nd after = ', accu_nd
|
|
if(accu_nd .gt. 1d-12) then
|
|
print*, ' accu_nd =', accu_nd
|
|
print*, ' your strategy for degenerates orbitals failed !'
|
|
print*, m, 'deg on', i
|
|
stop
|
|
endif
|
|
|
|
deallocate(S)
|
|
|
|
! ---
|
|
|
|
!call impose_orthog_svd(n, m, L)
|
|
!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_inverse(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_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_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 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_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 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)
|
|
|
|
! ---
|
|
|
|
! R <-- R x V x D^{-0.5}
|
|
! L <-- L x U x D^{-0.5}
|
|
|
|
do i = 1, m
|
|
do j = 1, m
|
|
V(j,i) = V(j,i) * D(i)
|
|
U(j,i) = U(j,i) * D(i)
|
|
enddo
|
|
enddo
|
|
|
|
allocate(tmp(n,m))
|
|
tmp(:,:) = R(:,:)
|
|
call dgemm( 'N', 'N', n, m, m, 1.d0 &
|
|
, tmp, size(tmp, 1), V, size(V, 1) &
|
|
, 0.d0, R, size(R, 1))
|
|
|
|
tmp(:,:) = L(:,:)
|
|
call dgemm( 'N', 'N', n, m, m, 1.d0 &
|
|
, tmp, size(tmp, 1), U, size(U, 1) &
|
|
, 0.d0, L, size(L, 1))
|
|
|
|
deallocate(tmp, U, V, D)
|
|
|
|
end
|
|
|
|
! ---
|
|
|
|
subroutine impose_biorthog_inverse(n, m, L, R)
|
|
|
|
implicit none
|
|
|
|
integer, intent(in) :: n, m
|
|
double precision, intent(inout) :: L(n,m)
|
|
double precision, intent(in) :: R(n,m)
|
|
double precision, allocatable :: Lt(:,:),S(:,:)
|
|
integer :: i,j
|
|
allocate(Lt(m,n))
|
|
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
|
|
|
|
call get_pseudo_inverse(R,n,n,m,Lt,m,1.d-6)
|
|
do i = 1, m
|
|
do j = 1, n
|
|
L(j,i) = Lt(i,j)
|
|
enddo
|
|
enddo
|
|
! ---
|
|
|
|
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,Lt)
|
|
|
|
|
|
end
|
|
|
|
! ---
|
|
|
|
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 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 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
|
|
|
|
! ---
|
|
|