PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-11-14 18:13:51 +01:00
Code Releases Activity
528bf20e1e
QuantumPackage/src/non_hermit_dav/utils.irp.f
eginer ca4cdf56d5 added non_hermit_dav
2023-02-06 19:03:22 +01:00

326 lines
9.2 KiB
Fortran
Raw Blame History

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
View Git Blame Copy Permalink