program RSVD_ijHkl_det

  BEGIN_DOC
  !
  ! decompose < di dj | H | dk dl > by a RSVD:
  !
  ! H_{ik jl} = < di dj | H | dk dl >
  !           = \sum_{s=1}^{n_svd} U_{ik s} V_{jl s}
  !           = \sum_{s=1}^{n_svd} U_{ik s} U_{jl s}
  !
  ! where n_svd << min( NaxNa , NbxNb )
  !
  END_DOC
  
  read_wf = .true.
  TOUCH read_wf

  call run()

end 

  


subroutine run()

  implicit none

  integer(bit_kind)             :: det1(N_int,2), det2(N_int,2)
  integer                       :: degree, i_state
  double precision              :: h12

  integer*8                     :: dim_Hdet
  integer*8,        allocatable :: Hdet_ik(:), Hdet_jl(:)
  double precision, allocatable :: Hdet_v(:)

  integer                       :: dim_U , dim_V, dim_RSVD

  double precision              :: t_beg, t_end

  call wall_time(t_beg)

  i_state = 1
  det1(:,1) = psi_det_alpha_unique(:,1)
  det2(:,1) = psi_det_alpha_unique(:,1)
  det1(:,2) = psi_det_beta_unique(:,1)
  det2(:,2) = psi_det_beta_unique(:,1)
  call get_excitation_degree_spin(det1(1,1),det2(1,1),degree,N_int)
  call get_excitation_degree(det1,det2,degree,N_int)
  call i_H_j(det1, det2, N_int, h12)


  ! sparse representation of < di dj | H | dk dl >
  !dim_Hdet = 7486379362 
  dim_Hdet = 1000000
  allocate( Hdet_ik(dim_Hdet), Hdet_jl(dim_Hdet), Hdet_v(dim_Hdet) )
  call const_ij_hdet_kl_sparse(dim_hdet, hdet_ik, hdet_jl, hdet_v)
 

  ! decompose Hdet by a Randomized SVD
  dim_U    = n_det_alpha_unique * n_det_alpha_unique
  dim_V    = n_det_beta_unique  * n_det_beta_unique
  dim_RSVD = min(dim_U, dim_V)
  dim_RSVD = min(5000, dim_RSVD)

  call perform_RSVD(dim_U, dim_V, dim_RSVD, dim_Hdet, Hdet_ik, Hdet_jl, Hdet_v)
  deallocate( Hdet_ik, Hdet_jl, Hdet_v )

  call wall_time(t_end)
  print *, ' end after (min)', (t_end-t_beg)/60.d0
 
end 
! _________________________________________________________________________________________________
! _________________________________________________________________________________________________
! _________________________________________________________________________________________________






! _________________________________________________________________________________________________
!
subroutine const_ij_Hdet_kl_sparse(dim_Hdet, Hdet_ik, Hdet_jl, Hdet_v)

  implicit none

  integer*8,         intent(in) :: dim_Hdet

  integer*8,        intent(out) :: Hdet_ik(dim_Hdet), Hdet_jl(dim_Hdet)
  double precision, intent(out) :: Hdet_v(dim_Hdet)

  integer(bit_kind)             :: det1(N_int,2), det2(N_int,2)
  integer                       :: degree

  integer                       :: na, nb, i, j, k, l, ii
  double precision              :: h12

  double precision              :: t1, t2

  print *, ""
  print *, " start const_ij_Hdet_kl_sparse"
  call wall_time(t1)

  na = n_det_alpha_unique
  nb = n_det_beta_unique

  ii = 0
 !$OMP PARALLEL DEFAULT(NONE)                                          &
 !$OMP          PRIVATE(i,j,k,l,ii,h12,det1,det2,degree)               &
 !$OMP          SHARED(na,nb,psi_det_alpha_unique,psi_det_beta_unique, &
 !$OMP                 N_int,Hdet_ik,Hdet_jl,Hdet_v)
 !$OMP DO
  do l = 1, nb
    det2(:,2) = psi_det_beta_unique(:,l)
    do j = 1, nb
      det1(:,2) = psi_det_beta_unique(:,j)

      call get_excitation_degree_spin(det1(1,2),det2(1,2),degree,N_int)
      if(degree .gt. 2) cycle

      do k = 1, na
        det2(:,1) = psi_det_alpha_unique(:,k)
        do i = 1, na
          det1(:,1) = psi_det_alpha_unique(:,i)

          call get_excitation_degree(det1,det2,degree,N_int)
          if(degree .gt. 2) cycle

          call i_H_j(det1, det2, N_int, h12)
          if(dabs(h12) .le. (1d-15)) cycle

          ii          = ii + 1
          Hdet_ik(ii) = (i-1)*na + k
          Hdet_jl(ii) = (j-1)*na + l
          Hdet_v (ii) = h12
        enddo
      enddo
    enddo
  enddo
 !$OMP END DO
 !$OMP END PARALLEL

  call wall_time(t2)
  print *, " end const_ij_Hdet_kl_sparse after (min) ", (t2-t1)/60.
  print *, ""

  if( ii .ne. dim_Hdet) then
    print*, ' error in const_ij_Hdet_kl_sparse'
    print*, ' dim supposed = ', dim_Hdet
    print*, ' dim foundedi = ', ii
    stop
  end if

  return
end subroutine const_ij_Hdet_kl_sparse
! _________________________________________________________________________________________________
! _________________________________________________________________________________________________





! _________________________________________________________________________________________________
!
subroutine perform_RSVD(dim_U, dim_V, dim_RSVD, dim_Hdet, Hdet_ik, Hdet_jl, Hdet_v)

  include 'constants.include.F'

  implicit none

  integer,          intent(in)  :: dim_U, dim_V, dim_RSVD
  integer*8,        intent(in)  :: dim_Hdet
  integer*8,        intent(in)  :: Hdet_ik(dim_Hdet), Hdet_jl(dim_Hdet)
  double precision, intent(in)  :: Hdet_v(dim_Hdet)

  integer                       :: na, nb
  integer                       :: iter, iter_max
  integer                       :: ii, rr, ik, jl
  double precision, allocatable :: GR(:,:)
  double precision, allocatable :: Z(:,:), P(:,:), Yt(:,:), UYt(:,:)
  double precision, allocatable :: U(:,:), V(:,:), D(:)

  na = n_det_alpha_unique
  nb = n_det_beta_unique

  ! Z = Hdet x G
  ! G: random gauss dist
  allocate( Z(dim_U,dim_RSVD) )
  Z(:,:) = 0
 !$OMP PARALLEL DEFAULT(NONE) PRIVATE(ii, rr, ik, GR)        &
 !$OMP          SHARED(dim_Hdet, dim_RSVD, Z, Hdet_ik, Hdet_v)
  allocate( GR(dim_Hdet,2) )
 !$OMP DO
  do rr = 1, dim_RSVD
   call random_number(GR)
   GR(:,1) = dsqrt(-2.d0*dlog(GR(:,1)))
   GR(:,1) = GR(:,1) * dcos(dtwo_pi*GR(:,2))
   do ii = 1, dim_Hdet
     ik       = Hdet_ik(ii)
     !jl       = Hdet_jl(ii)
     Z(ik,rr) = Z(ik,rr) + Hdet_v(ii) * GR(ii,1)
   enddo
  enddo
 !$OMP END DO
  deallocate(GR)
 !$OMP END PARALLEL


  ! parameter
  iter_max = 10
  allocate( P(dim_V,dim_RSVD) )
  ! ------------------------------------------------------------------------------------------
  ! Power iterations
  !
  do iter = 1, iter_max

    print *, ' power iteration ', iter, '/', iter_max

 !$OMP PARALLEL DEFAULT(NONE) PRIVATE(ii, rr, ik, jl)                    &
 !$OMP          SHARED(dim_Hdet, dim_RSVD, Z, P, Hdet_ik, Hdet_jl, Hdet_v)

 !$OMP DO
    ! P = Hdet.T x Z
    do rr = 1, dim_RSVD
      P(:,rr) = 0.d0
      do ii = 1, dim_Hdet
        ik       = Hdet_ik(ii)
        jl       = Hdet_jl(ii)
        P(jl,rr) = P(jl,rr) + Hdet_v(ii) * Z(ik,rr)
      enddo
    enddo
 !$OMP END DO

 !$OMP BARRIER

 !$OMP DO
    do rr = 1, dim_RSVD
      Z(:,rr) = 0.d0
      do ii = 1, dim_Hdet
        ik       = Hdet_ik(ii)
        jl       = Hdet_jl(ii)
        Z(ik,rr) = Z(ik,rr) + Hdet_v(ii) * P(jl,rr)
      enddo
    enddo
 !$OMP END DO

 !$OMP END PARALLEL

    call ortho_qr(Z, size(Z,1), dim_U, dim_RSVD)

  enddo
  !
  ! ------------------------------------------------------------------------------------------

  ! Y  = Z.T x Hdet
  ! Yt = Hdet.T x Z
  allocate( Yt(dim_V,dim_RSVD) )
 !$OMP PARALLEL DEFAULT(NONE) PRIVATE(ii, rr, ik, jl)                     &
 !$OMP          SHARED(dim_Hdet, dim_RSVD, Z, Yt, Hdet_ik, Hdet_jl, Hdet_v)
 !$OMP DO
  do rr = 1, dim_RSVD
    Yt(:,rr) = 0.d0
    do ii = 1, dim_Hdet
      ik       = Hdet_ik(ii)
      jl       = Hdet_jl(ii)
      Yt(jl,rr) = Yt(jl,rr) + Z(ik,rr) * Hdet_v(ii)
    enddo
  enddo
 !$OMP END DO
 !$OMP END PARALLEL

  ! Y  = U x D x Vt 
  ! Yt = V x D x Ut 
  allocate( D(dim_RSVD), V(dim_V,dim_RSVD), UYt(dim_RSVD,dim_RSVD) )
  call svd(Yt, size(Yt,1), V, size(V,1), D, UYt, size(UYt,1), dim_V, dim_RSVD)
  deallocate( Yt )

  ! U = Z x UY
  allocate( U(dim_U,dim_RSVD) )
  call dgemm('N', 'T', dim_U, dim_RSVD, dim_RSVD, 1.d0, Z, size(Z,1), UYt, size(UYt,1), 0.d0, U, size(U,1))
  deallocate( UYt, Z )

  open(unit=41, file='u_det.txt', action='write')
    do rr = 1, dim_RSVD
      do ii = 1, dim_U
        write(41,*) U(ii,rr)
      enddo 
    enddo 
  close(41)

  ! v = u because H is symmetric
  !open(unit=41, file='v_det.txt', action='write')
  !  do rr = 1, dim_RSVD
  !    do ii = 1, dim_V
  !      write(41,*) V(ii,rr)
  !    enddo 
  !  enddo 
  !close(41)

  open(unit=41, file='d_det.txt', action='write')
    do rr = 1, dim_RSVD
      write(41,*) D(rr), sum( D(1:rr)**2 )
    enddo 
  close(41)

  deallocate( U, D, V )

end subroutine perform_RSVD
! _________________________________________________________________________________________________
! _________________________________________________________________________________________________