program svdwf
  implicit none
  BEGIN_DOC
! TODO : Make the SVD of the alpha-beta wave function and print singular values.
  END_DOC
  read_wf = .True.
  TOUCH read_wf
  call run()
end

subroutine run
  implicit none
  double precision, allocatable :: U(:,:), Vt(:,:), D(:), A(:,:)
  integer :: i, j, k, p, q
  double precision :: entropy

  allocate( A (n_det_alpha_unique, n_det_beta_unique), &
            U (n_det_alpha_unique, n_det_alpha_unique), &
            Vt(n_det_beta_unique, n_det_beta_unique),  &
            D(max(n_det_beta_unique,n_det_alpha_unique)) )

  A = 0.D0
  do k=1,N_det
    i = psi_bilinear_matrix_rows(k)
    j = psi_bilinear_matrix_columns(k)
    A(i,j) = psi_bilinear_matrix_values(k,1)
  enddo

  call randomized_svd(A, size(A,1), &
           U, size(U,1), D, Vt, size(Vt,1), n_det_alpha_unique, n_det_beta_unique, &
           6,1000)

  entropy = 0.d0
  do i=1,n_det_beta_unique
    print *, i, real(D(i)), real(D(i)**2), real(sum(D(1:i)**2))
    entropy -= D(i) * dlog(D(i))
    if (D(i) < 1.d-15) then
      k = i
      exit
    endif
  enddo
  print *, 'threshold: ', 2.858 * D(k/2)
  print *, 'Entropy  : ', entropy
!  do i=1,n_det_alpha_unique
!    print '(I6,4(X,F12.8))', i, U(i,1:4)
!  enddo
!  print *, ''
!  do i=1,n_det_beta_unique
!    print '(I6,4(X,F12.8))', i, Vt(1:4,i)
!  enddo
end