mirror of
https://gitlab.com/scemama/qp_plugins_scemama.git
synced 2024-12-22 12:23:37 +01:00
Added Abdallah's file
This commit is contained in:
parent
b177795477
commit
7039831efc
167
devel/svdwf/Evar_TruncSVD.irp.f
Normal file
167
devel/svdwf/Evar_TruncSVD.irp.f
Normal file
@ -0,0 +1,167 @@
|
||||
program Evar_TruncSVD
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! study energy variation with truncated SVD
|
||||
END_DOC
|
||||
read_wf = .True.
|
||||
TOUCH read_wf
|
||||
! !!!
|
||||
call run()
|
||||
! !!!
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
subroutine run
|
||||
implicit none
|
||||
include 'constants.include.F'
|
||||
double precision, allocatable :: A(:,:), U(:,:), V(:,:), D(:)
|
||||
integer :: r, i, j, k, l, m, n, iter, iter_max
|
||||
double precision, allocatable :: Z(:,:), P(:,:), Yt(:,:), UYt(:,:), r1(:,:)
|
||||
! !!!
|
||||
m = n_det_alpha_unique
|
||||
n = n_det_beta_unique
|
||||
r = n
|
||||
print *, 'matrix:', m,'x',n
|
||||
print *, 'N det:', N_det
|
||||
print *, 'rank = ', r
|
||||
iter_max = 20
|
||||
! !!!
|
||||
allocate( Z(m,r) , P(n,r) ) ! Z(m,r) = A(m,n) @ P(n,r)
|
||||
Z(:,:) = 0.d0
|
||||
! !!!
|
||||
! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !
|
||||
! first we apply a RSVD for a pre-fixed rank (r)
|
||||
! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !
|
||||
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l,r1)
|
||||
allocate(r1(N_det,2))
|
||||
!$OMP DO
|
||||
do l=1,r
|
||||
call random_number(r1)
|
||||
r1(:,1) = dsqrt(-2.d0*dlog(r1(:,1)))
|
||||
r1(:,1) = r1(:,1) * dcos(dtwo_pi*r1(:,2))
|
||||
do k=1,N_det
|
||||
i = psi_bilinear_matrix_rows(k)
|
||||
j = psi_bilinear_matrix_columns(k)
|
||||
Z(i,l) = Z(i,l) + psi_bilinear_matrix_values(k,1) * r1(k,1)
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
deallocate(r1)
|
||||
!$OMP END PARALLEL
|
||||
! !!!
|
||||
! Power iterations
|
||||
! !!!
|
||||
do iter=1,iter_max
|
||||
! !!!
|
||||
print *, 'Power iteration ', iter, '/', 20
|
||||
! !!!
|
||||
! P(n,r) = At(n,m) @ Z(m,r)
|
||||
! !!!
|
||||
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l)
|
||||
!$OMP DO
|
||||
do l=1,r
|
||||
P(:,l) = 0.d0
|
||||
do k=1,N_det
|
||||
i = psi_bilinear_matrix_rows(k)
|
||||
j = psi_bilinear_matrix_columns(k)
|
||||
P(j,l) = P(j,l) + psi_bilinear_matrix_values(k,1) * Z(i,l)
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
! !!!
|
||||
! Z(m,r) = A(m,n) @ P(n,r)
|
||||
! !!!
|
||||
!$OMP BARRIER
|
||||
!$OMP DO
|
||||
do l=1,r
|
||||
Z(:,l) = 0.d0
|
||||
do k=1,N_det
|
||||
i = psi_bilinear_matrix_rows(k)
|
||||
j = psi_bilinear_matrix_columns(k)
|
||||
Z(i,l) = Z(i,l) + psi_bilinear_matrix_values(k,1) * P(j,l)
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
! !!!
|
||||
! Compute QR: at return: Q is Z(m,r)
|
||||
! !!!
|
||||
call ortho_qr(Z,size(Z,1),m,r)
|
||||
! !!!
|
||||
enddo
|
||||
! !!!
|
||||
! Y(r,n) = Zt(r,m) @ A(m,n) or Yt(n,r) = At(n,m) @ Z(m,r)
|
||||
! !!!
|
||||
allocate(Yt(n,r))
|
||||
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l)
|
||||
!$OMP DO
|
||||
do l=1,r
|
||||
do k=1,n
|
||||
Yt(k,l) = 0.d0
|
||||
enddo
|
||||
do k=1,N_det
|
||||
i = psi_bilinear_matrix_rows(k)
|
||||
j = psi_bilinear_matrix_columns(k)
|
||||
Yt(j,l) = Yt(j,l) + Z(i,l) * psi_bilinear_matrix_values(k,1)
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
! !!!
|
||||
! Y = UY @ D @ Vt or Yt = V @ Dt @ UYt
|
||||
! !!!
|
||||
allocate(D(r),V(n,r), UYt(r,r))
|
||||
! !!!
|
||||
call svd(Yt,size(Yt,1),V,size(V,1),D,UYt,size(UYt,1),n,r)
|
||||
deallocate(Yt)
|
||||
! !!!
|
||||
! U(m,r) = Z(m,r) @ UY(r,r) or U = Z @ (UYt).T
|
||||
! !!!
|
||||
allocate(U(m,r))
|
||||
call dgemm('N','T',m,r,r,1.d0,Z,size(Z,1),UYt,size(UYt,1),0.d0,U,size(U,1))
|
||||
deallocate(UYt,Z)
|
||||
! !!!
|
||||
!do i=1,r
|
||||
! print *, i, real(D(i)), real(D(i)**2), real(sum(D(1:i)**2))
|
||||
! if (D(i) < 1.d-15) then
|
||||
! k = i
|
||||
! exit
|
||||
! endif
|
||||
!enddo
|
||||
!print *, 'threshold: ', 2.858 * D(k/2)
|
||||
! !!!
|
||||
! Build the new determinant: U @ D @ Vt
|
||||
! !!!
|
||||
!!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l)
|
||||
!!$OMP DO
|
||||
!!
|
||||
!print *, 'ok 1'
|
||||
!N_det = m * n
|
||||
!print *, 'ok 11'
|
||||
!TOUCH N_det
|
||||
!psi_bilinear_matrix_values(:,1) = 0.d0
|
||||
!TOUCH psi_bilinear_matrix_values
|
||||
! print *, size(psi_bilinear_matrix_values,1), size(D), size(U,1), size(U,2), size(V,1), size(V,2)
|
||||
print*, PSI_energy(1) + nuclear_repulsion
|
||||
psi_bilinear_matrix(:,:,:) = 0.d0
|
||||
do r = 1, n
|
||||
call generate_all_alpha_beta_det_products
|
||||
do i = 1, N_det_beta_unique
|
||||
do j = 1, N_det_alpha_unique
|
||||
psi_bilinear_matrix(j,i,1) = 0.d0
|
||||
do l = 1, r
|
||||
psi_bilinear_matrix(j,i,1) = psi_bilinear_matrix(j,i,1) + D(l) * U(j,l) * V(i,l)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
TOUCH psi_bilinear_matrix
|
||||
call update_wf_of_psi_bilinear_matrix(.False.)
|
||||
print*, r, PSI_energy(1) + nuclear_repulsion !CI_energy(1)
|
||||
enddo
|
||||
!!$OMP END DO
|
||||
!!$OMP END PARALLEL
|
||||
deallocate(U,D,V)
|
||||
! !!!
|
||||
end
|
||||
@ -1 +1,2 @@
|
||||
determinants
|
||||
davidson_undressed
|
||||
|
||||
@ -22,13 +22,17 @@ subroutine run
|
||||
|
||||
allocate(Z(m,r))
|
||||
|
||||
! Z(m,r) = A(m,n).P(n,r)
|
||||
Z(:,:) = 0.d0
|
||||
! Z(m,r) = A(m,n).P(n,r)
|
||||
do j=1,r
|
||||
do i=1,m
|
||||
Z(i,j) = 0.d0
|
||||
enddo
|
||||
enddo
|
||||
allocate(P(n,r))
|
||||
|
||||
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l,r1)
|
||||
allocate(r1(N_det,2))
|
||||
!$OMP DO
|
||||
!$OMP DO
|
||||
do l=1,r
|
||||
call random_number(r1)
|
||||
r1(:,1) = dsqrt(-2.d0*dlog(r1(:,1)))
|
||||
@ -36,7 +40,7 @@ subroutine run
|
||||
do k=1,N_det
|
||||
i = psi_bilinear_matrix_rows(k)
|
||||
j = psi_bilinear_matrix_columns(k)
|
||||
Z(i,l) = Z(i,l) + psi_bilinear_matrix_values(k,1) * r1(k,1)
|
||||
Z(i,l) = Z(i,l) + psi_bilinear_matrix_values(k,1) * r1(k,1)
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
@ -98,7 +102,7 @@ subroutine run
|
||||
call svd(Yt,size(Yt,1),V,size(V,1),D,UYt,size(UYt,1),n,r)
|
||||
deallocate(Yt)
|
||||
|
||||
! U(m,r) = Z(m,r).UY(r,r)
|
||||
! U(m,r) = Z(m,r).UY(r,r)
|
||||
allocate(U(m,r))
|
||||
call dgemm('N','T',m,r,r,1.d0,Z,size(Z,1),UYt,size(UYt,1),0.d0,U,size(U,1))
|
||||
deallocate(UYt,Z)
|
||||
|
||||
Loading…
Reference in New Issue
Block a user