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, s2_values(1) !CI_energy(1)
    call save_wavefunction()
  enddo
  deallocate(U,D,V)
  ! !!!
end