mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-06-01 10:15:18 +02:00
326 lines
9.2 KiB
Fortran
326 lines
9.2 KiB
Fortran
|
|
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
|