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'
  ! !!!
  integer                       :: m, n, i_state
  double precision              :: error_thr, error_RRRSVD, norm_psi, norm_SVD, err_verif, err_tmp
  integer                       :: i, j, k, l, It, PowerIt
  integer                       :: It_max, PowerIt_max, nb_oversamp
  integer                       :: r_init, delta_r, low_rank
  double precision, allocatable :: B_old(:,:), Q_old(:,:)
  double precision, allocatable :: UB(:,:), D(:), Vt(:,:), U(:,:)
  double precision, allocatable :: P(:,:), r1(:,:)
  double precision, allocatable :: Q_new(:,:), Q_tmp(:,:), Q_mult(:,:)
  double precision, allocatable :: B_new(:,:), B_tmp(:,:)
  double precision, allocatable :: URSVD(:,:), DRSVD(:), VtRSVD(:,:)
  double precision, allocatable :: Uverif(:,:), Dverif(:), Vtverif(:,:), Averif(:,:)
  double precision, allocatable :: Uezfio(:,:,:), Dezfio(:,:), Vezfio(:,:,:)
  double precision :: tmp
  ! !!!
  i_state = 1
  m       = n_det_alpha_unique
  n       = n_det_beta_unique
  ! !!!
  !open(111, file = 'data_to_python.txt', action = 'WRITE' )  
  !  do k = 1, N_det
  !    write(111, '(I8, 5X, I8, 5X, E15.7)' ) psi_bilinear_matrix_rows(k), psi_bilinear_matrix_columns(k), psi_bilinear_matrix_values(k,i_state)
  !  end do  
  !close(111)
  ! !!!
  ! !!!
  print *, 'matrix:', m,'x',n
  print *, 'N det:', N_det
  ! !!!
  r_init  = 78
  delta_r = 5
  ! !!!
  PowerIt_max = 10
  nb_oversamp = 10
  ! !!!
  It_max    = 10
  error_thr = 1.d-3 !! don't touche it but rather increase r_init
  ! !!!
  ! !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! !
  ! !!! ~      Rank Revealing Randomized SVD      ~ !!! !
  ! !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! !
  ! !!!
  norm_psi = 0.d0
  do k = 1, N_det
    norm_psi = norm_psi + psi_bilinear_matrix_values(k,i_state) * psi_bilinear_matrix_values(k,i_state)
  enddo
  ! !!!
  ! !!! !!! !!! !!! !!! !!! !!! !!! !!! !!! !
  ! !!!  build initial QB decomposition !!! !
  ! !!! !!! !!! !!! !!! !!! !!! !!! !!! !!! !
  ! !!!
  allocate( Q_old(m, r_init) )
  !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,k,l,r1)
  allocate( r1(N_det,2) )
  !$OMP DO
  do l = 1, r_init
    Q_old(:,l) = 0.d0
    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)
      Q_old(i,l) = Q_old(i,l) + psi_bilinear_matrix_values(k,i_state) * r1(k,1)
    enddo
  enddo
  !$OMP END DO
  deallocate(r1)
  !$OMP END PARALLEL
  ! !!!
  ! power scheme
  ! !!!
  allocate( P(n, r_init) ) 
  do PowerIt = 1, PowerIt_max
    ! !!!
    call my_ortho_qr(Q_old, size(Q_old,1), m, r_init)
    ! !!!
    !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l)
    !$OMP DO
    do l = 1, r_init
      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,i_state) * Q_old(i,l)
      enddo
    enddo
    !$OMP END DO
    !$OMP END PARALLEL
    ! !!!
    call my_ortho_qr(P, size(P,1), n, r_init)
    ! !!!
    !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l)
    !$OMP DO
    do l = 1, r_init
      Q_old(:,l) = 0.d0
      do k = 1 , N_det
        i          = psi_bilinear_matrix_rows(k)
        j          = psi_bilinear_matrix_columns(k)
        Q_old(i,l) = Q_old(i,l) + psi_bilinear_matrix_values(k,i_state) * P(j,l)
      enddo
    enddo
    !$OMP END DO
    !$OMP END PARALLEL
    ! !!!
  enddo
  deallocate( P )
  ! !!!
  call my_ortho_qr(Q_old, size(Q_old,1), m, r_init)
  ! !!!
  allocate( B_old(r_init,n) )
  !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l)
  !$OMP DO
  do l = 1, r_init
    B_old(l,:) = 0.d0
    do k = 1, N_det
      i          = psi_bilinear_matrix_rows(k)
      j          = psi_bilinear_matrix_columns(k)
      B_old(l,j) = B_old(l,j) + Q_old(i,l) * psi_bilinear_matrix_values(k,i_state)
    enddo
  enddo
  !$OMP END DO
  !$OMP END PARALLEL
  norm_SVD = 0.d0
  do j = 1, n
    do l = 1, r_init
      norm_SVD = norm_SVD + B_old(l,j) * B_old(l,j)
    enddo
  enddo
  error_RRRSVD  = dabs( norm_psi - norm_SVD ) / norm_psi
  It       = 1
  low_rank = r_init
  print *, It, low_rank, error_RRRSVD
  ! !!!
  ! !!! !!! !!! !!! !!! !!! !!! !!! !!! !!! !!! !!! !
  ! !!! incrementally build up QB decomposition !!! !
  ! !!! !!! !!! !!! !!! !!! !!! !!! !!! !!! !!! !!! !
  ! !!!
  do while( ( error_RRRSVD.gt.error_thr ).and.( It.lt.It_max ).and.( low_rank.lt.(min(m,n)-delta_r) ) )
    ! !!!
    allocate( Q_new(m,delta_r) ) 
    allocate( r1(N_det, 2) )
    do l = 1, delta_r
      Q_new(:,l) = 0.d0
      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)
        Q_new(i,l) = Q_new(i,l) + psi_bilinear_matrix_values(k,i_state) * r1(k,1)
      enddo
    enddo
    deallocate(r1)
    ! !!!
    ! orthogonalization with Q_old:  Q_new = Q_new - Q_old @ Q_old.T @ Q_new
    ! !!!
    !allocate( Q_mult(m,m) )
    !call dgemm( 'N', 'T', m, m, low_rank, +1.d0, Q_old,  size(Q_old,1),  Q_old, size(Q_old,1), 0.d0, Q_mult, size(Q_mult,1) )
    !!do i = 1, m
    !!  do j = 1, m
    !!    Q_mult(j,i) = 0.d0
    !!    do l = 1, low_rank
    !!      Q_mult(j,i) = Q_mult(j,i) +  Q_old(i,l) * Q_old(j,l)
    !!    enddo
    !!  enddo
    !!enddo
    !!call dgemm( 'N', 'N', m, delta_r, m,  -1.d0, Q_mult, size(Q_mult,1), Q_new, size(Q_new,1), 1.d0, Q_new,  size(Q_new,1)  )
    !do l = 1, delta_r
    !  do i = 1, m
    !    tmp = 0.d0
    !    do j = 1, m
    !      tmp = tmp + Q_mult(i,j) * Q_new(j,l)
    !    enddo
    !    Q_new(i,l) = Q_new(i,l) - tmp
    !  enddo
    !enddo
    !deallocate( Q_mult )
    ! !!!
    ! power scheme
    ! !!!
    allocate( P(n, delta_r) )
    do PowerIt = 1, PowerIt_max
      ! !!!
      call my_ortho_qr(Q_new, size(Q_new,1), m, delta_r)
      ! !!!
      do l = 1, delta_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,i_state) * Q_new(i,l)
        enddo
      enddo
      ! !!!
      call my_ortho_qr(P, size(P,1), n, delta_r)
      ! !!!
      do l = 1, delta_r
        Q_new(:,l) = 0.d0
        do k = 1, N_det
          i          = psi_bilinear_matrix_rows(k)
          j          = psi_bilinear_matrix_columns(k)
          Q_new(i,l) = Q_new(i,l) + psi_bilinear_matrix_values(k,i_state) * P(j,l)
        enddo
      enddo
      ! !!!
    enddo
    deallocate( P )
    ! !!!
    ! orthogonalization with Q_old:  Q_new = Q_new - Q_old @ Q_old.T @ Q_new
    ! !!!
    allocate( Q_mult(m,m) )
    call dgemm( 'N', 'T', m, m, low_rank, +1.d0, Q_old,  size(Q_old,1),  Q_old, size(Q_old,1), 0.d0, Q_mult, size(Q_mult,1) )
    !do i = 1, m
    !  do j = 1, m
    !    Q_mult(j,i) = 0.d0
    !    do l = 1, low_rank
    !      Q_mult(j,i) = Q_mult(j,i) +  Q_old(i,l) * Q_old(j,l)
    !    enddo
    !  enddo
    !enddo
    !call dgemm( 'N', 'N', m, delta_r, m,  -1.d0, Q_mult, size(Q_mult,1), Q_new, size(Q_new,1), 1.d0, Q_new,  size(Q_new,1)  )
    do l = 1, delta_r
      do i = 1, m
        tmp = 0.d0
        do j = 1, m
          tmp = tmp + Q_mult(i,j) * Q_new(j,l)
        enddo
        Q_new(i,l) = Q_new(i,l) - tmp
      enddo
    enddo
    deallocate( Q_mult )
    ! !!!
    call my_ortho_qr(Q_new, size(Q_new,1), m, delta_r)
    ! !!!
    allocate( B_new(delta_r,n) )
    do l = 1 , delta_r
      B_new(l,:) = 0.d0
      do k = 1 , N_det
        i          = psi_bilinear_matrix_rows(k)
        j          = psi_bilinear_matrix_columns(k)
        B_new(l,j) = B_new(l,j) + Q_new(i,l) * psi_bilinear_matrix_values(k,i_state)
      enddo
    enddo
    ! !!!
    do j = 1, n
      do l = 1, delta_r
        norm_SVD = norm_SVD + B_new(l,j) * B_new(l,j)
      enddo
    enddo
    ! !!!
    error_RRRSVD = dabs( norm_psi - norm_SVD ) / norm_psi
    It           = It + 1
    low_rank     = low_rank + delta_r
    ! !!!
    ! build up approximate basis:
    !         Q_old  = np.append(Q_new, Q_old, axis=1)
    !         B      = np.append(B_new, B, axis=0)
    ! !!!
    allocate( Q_tmp(m,low_rank) , B_tmp(low_rank,n) )
    ! !!!
    !do l = 1, delta_r
    !  do i = 1, m
    !    Q_tmp(i,l) = Q_new(i,l)
    !  enddo
    !enddo
    !do l = 1 , low_rank-delta_r
    !  do i = 1, m
    !    Q_tmp(i,l+delta_r) = Q_old(i,l)
    !  enddo
    !enddo
    !do i = 1, n
    !  do l = 1 , delta_r
    !    B_tmp(l,i) = B_new(l,i)
    !  enddo
    !enddo
    !do i = 1, n
    !  do l = 1 , low_rank-delta_r
    !    B_tmp(l+delta_r,i) = B_old(l,i)
    !  enddo
    !enddo
    ! !!!
    do i = 1, m
      do l = 1, low_rank-delta_r
        Q_tmp(i,l) = Q_old(i,l)
      enddo
      do l = 1, delta_r
        Q_tmp(i,l+low_rank-delta_r) = Q_new(i,l)
      enddo
    enddo
    do i = 1, n
      do l = 1 , low_rank-delta_r
        B_tmp(l,i) = B_old(l,i)
      enddo
      do l = 1 , delta_r
        B_tmp(l+low_rank-delta_r,i) = B_new(l,i)
      enddo
    enddo
    ! !!!
    deallocate( Q_old, B_old, Q_new, B_new )
    allocate( Q_old(m,low_rank) , B_old(low_rank,n) )
    ! !!!
    do l = 1, low_rank
      do i = 1, m
        Q_old(i,l) = Q_tmp(i,l)
      enddo
    enddo
    do l = 1, n
      do i = 1, low_rank
        B_old(i,l) = B_tmp(i,l)
      enddo
    enddo
    deallocate( Q_tmp , B_tmp )
    ! !!!
  print *, It, low_rank, error_RRRSVD
  ! !!!
  enddo
  ! !!!
  ! !!! !!! !!! !!! !!! !!! !!! !!! !!! !!! !!! !!! !
  ! !!! !!! !!! !!! !!! !!! !!! !!! !!! !!! !!! !!! !
  ! !!!
  allocate( UB(low_rank,low_rank), D(low_rank), Vt(low_rank,n) )
  call svd_s(B_old, size(B_old,1), UB, size(UB,1), D, Vt, size(Vt,1), low_rank, n)
  deallocate(B_old)
  print*, 'ok 1'
  ! !!!
  allocate( U(m,low_rank) )
  call dgemm('N', 'N', m, low_rank, low_rank, 1.d0, Q_old, size(Q_old,1), UB, size(UB,1), 0.d0, U, size(U,1))
  deallocate( Q_old,UB )
  print*, 'ok 2'
  ! !!!
  ! !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! !
  ! !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! !
  ! !!!
  allocate( URSVD(m,low_rank), DRSVD(low_rank), VtRSVD(low_rank,n) )
  call RSVD( i_state, low_rank, PowerIt_max, nb_oversamp, URSVD, DRSVD, VtRSVD )
  print*, 'ok 3'
  ! !!!
  ! !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! !
  ! !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! !
  ! !!!
  allocate( Averif(m,n), Uverif(m,n), Dverif(n), Vtverif(n,n) )
  do i = 1, n
    Averif(:,i) = 1d-16
  enddo
  do k = 1, N_det
    i           = psi_bilinear_matrix_rows(k)
    j           = psi_bilinear_matrix_columns(k)
    Averif(i,j) = psi_bilinear_matrix_values(k,i_state)
  enddo
  call svd_s( Averif, size(Averif,1), Uverif, size(Uverif,1), Dverif, Vtverif, size(Vtverif,1), m, n)
  print*, 'ok 4'
  ! !!!
  err_verif = 0.d0
  do j = 1, n
    do i = 1, m
    err_tmp = 0.d0
    do l = 1, low_rank
      err_tmp = err_tmp + Dverif(l) * Uverif(i,l) * Vtverif(l,j)
    enddo
    err_verif = err_verif + ( Averif(i,j) - err_tmp )**2.d0
    enddo
  enddo
  print*, 'err verif (%) = ', 100.d0 * dsqrt(err_verif/norm_psi)
  ! !!!
  ! !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! !
  ! !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! !
  ! !!!
  open(111, file = 'singular_values.txt', action = 'WRITE' )  
    do i = 1, low_rank  
      write(111, '(I8, 5X, E15.7, 2(5X, E15.7, E15.7) )' ) i, Dverif(i), D(i), 100.d0*dabs(D(i)-Dverif(i))/Dverif(i), DRSVD(i), 100.d0*dabs(DRSVD(i)-Dverif(i))/Dverif(i)
    end do  
  close(111)
  ! !!!
  !deallocate( Averif, Uverif, Dverif, Vtverif )
  ! !!!
  low_rank = 10
  ! !!!
  err_verif = 0.d0
  do j = 1, n
    do i = 1, m
    err_tmp = 0.d0
    do l = 1, low_rank
      err_tmp = err_tmp + Dverif(l) * Uverif(i,l) * Vtverif(l,j)
    enddo
    err_verif = err_verif + ( Averif(i,j) - err_tmp )**2.d0
    enddo
  enddo
  print*, 'err verif (%) = ', 100.d0 * dsqrt(err_verif/norm_psi)
  ! !!!
  print*, 'low_rank =', low_rank
  allocate(Uezfio(m,low_rank,1), Dezfio(low_rank,1), Vezfio(n,low_rank,1))
  do l = 1, low_rank
    Dezfio(l,1) = Dverif(l)
    do j = 1, m
      Uezfio(j,l,1) = Uverif(j,l)
    enddo
    do j = 1, n
      Vezfio(j,l,1) = Vtverif(l,j)
    enddo
  enddo
  deallocate( U, D, Vt )
  ! !!!
  call ezfio_set_spindeterminants_n_det(N_det)
  call ezfio_set_spindeterminants_n_states(N_states)
  call ezfio_set_spindeterminants_n_det_alpha(n_det_alpha_unique)
  call ezfio_set_spindeterminants_n_det_beta(n_det_beta_unique)
  call ezfio_set_spindeterminants_psi_coef_matrix_rows(psi_bilinear_matrix_rows)
  call ezfio_set_spindeterminants_psi_coef_matrix_columns(psi_bilinear_matrix_columns)
  call ezfio_set_spindeterminants_psi_coef_matrix_values(psi_bilinear_matrix_values)
  ! !!!
  call ezfio_set_spindeterminants_n_svd_coefs(low_rank)
  call ezfio_set_spindeterminants_psi_svd_alpha(Uezfio)
  call ezfio_set_spindeterminants_psi_svd_beta(Vezfio )
  call ezfio_set_spindeterminants_psi_svd_coefs(Dezfio)
  deallocate(Uezfio, Dezfio, Vezfio)
  ! !!!
  !print*, PSI_energy(1) + nuclear_repulsion
  !psi_bilinear_matrix(:,:,:) = 0.d0
  !do low_rank = n, 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*, low_rank, PSI_energy(1) + nuclear_repulsion !CI_energy(1)
  !enddo
  !deallocate(U,D,V)
  ! !!!
end