Merge branch 'master' of gitlab.com:scemama/qp_plugins_scemama

devel/svdwf/Evar_TruncSVD.irp.f

@@ -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

devel/svdwf/NEED

@@ -1 +1,2 @@
 determinants
+davidson_undressed

devel/svdwf/random_svd.irp.f

@@ -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
@@ -100,7 +104,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)