program svdwf
  implicit none
  BEGIN_DOC
! 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
  include 'constants.include.F'
  double precision, allocatable :: U(:,:), V(:,:), D(:), A(:,:)
  integer :: i, j, k, l, q, r, m, n, iter
  double precision,allocatable    :: Z(:,:), P(:,:), Yt(:,:), UYt(:,:)
  double precision :: r1,r2

  m = n_det_alpha_unique
  n = n_det_beta_unique

  r = min(1000,n)

  allocate(Z(m,r))

  ! Z(m,r) = A(m,n).P(n,r)                                                                    
  Z(:,:) = 0.d0
  do l=1,r
   do k=1,N_det
     i = psi_bilinear_matrix_rows(k)
     j = psi_bilinear_matrix_columns(k)
     call random_number(r1)
     call random_number(r2)
     r1 = dsqrt(-2.d0*dlog(r1))
     r2 = dtwo_pi*r2
     Z(i,l) = Z(i,l) + psi_bilinear_matrix_values(k,1) * r1*dcos(r2)
   enddo
  enddo

  ! Power iterations
  allocate(P(n,r))
  do iter=1,20
   ! P(n,r) = At(n,m).Z(m,r)
   P(:,:) = 0.d0
   do l=1,r
    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
   Z(:,:) = 0.d0
   do l=1,r
    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
   ! Compute QR
   call ortho_qr(Z,size(Z,1),m,r)
  enddo

  ! Y(r,n) = Zt(r,m).A(m,n)
  allocate(Yt(n,r))
  Yt(:,:) = 0.d0
  do l=1,r
    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

  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)                                                                   
  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)
  do i=1,m
    print '(I6,4(X,F12.8))', i, U(i,1:4)
  enddo
  print *, ''
  do i=1,n
    print '(I6,4(X,F12.8))', i, V(i,1:4)
  enddo

  deallocate(U,D,V)
end