diff --git a/devel/svdwf/Evar_TruncSVD.irp.f b/devel/svdwf/Evar_TruncSVD.irp.f
new file mode 100644
index 0000000..aa480e1
--- /dev/null
+++ b/devel/svdwf/Evar_TruncSVD.irp.f
@@ -0,0 +1,451 @@
+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
diff --git a/devel/svdwf/NEED b/devel/svdwf/NEED
index d3d4d2c..49488c9 100644
--- a/devel/svdwf/NEED
+++ b/devel/svdwf/NEED
@@ -1 +1,2 @@
 determinants
+davidson_undressed 
diff --git a/devel/svdwf/QR.py b/devel/svdwf/QR.py
new file mode 100644
index 0000000..262ee4e
--- /dev/null
+++ b/devel/svdwf/QR.py
@@ -0,0 +1,10 @@
+# !!!
+import numpy as np
+# !!!
+def QR_fact(X):
+    Q, R     = np.linalg.qr(X, mode="reduced")
+    D        = np.diag( np.sign( np.diag(R) ) )
+    Qunique  = np.dot(Q,D)
+    #Runique = np.dot(D,R)
+    return(Qunique)
+# !!!
\ No newline at end of file
diff --git a/devel/svdwf/R3SVD_AMMAR.py b/devel/svdwf/R3SVD_AMMAR.py
new file mode 100644
index 0000000..6a719b7
--- /dev/null
+++ b/devel/svdwf/R3SVD_AMMAR.py
@@ -0,0 +1,68 @@
+# !!!
+import numpy as np
+from QR   import QR_fact
+from RSVD import powit_RSVD
+# !!!
+def R3SVD_AMMAR(A, t, delta_t, npow, nb_oversamp, err_thr, maxit, tol):
+    # !!!
+    # build initial QB decomposition
+    # !!!
+    n = A.shape[1]
+    G = np.random.randn(n, t)
+    normA  = np.linalg.norm(A, ord='fro')**2
+    i_it = 0
+    rank = 0
+    Y = np.dot(A,G)
+    # The power scheme
+    for j in range(npow):
+        Q = QR_fact(Y)
+        Q = QR_fact( np.dot(A.T,Q) )
+        Y = np.dot(A,Q)
+    # orthogonalization process
+    Q_old = QR_fact(Y)
+    B     = np.dot(Q_old.T,A)
+    normB = np.linalg.norm(B, ord='fro')**2
+    # error percentage
+    errpr = abs( normA - normB ) / normA
+    rank += t
+    i_it += 1
+    print("iteration = {}, rank = {}, error = {}".format(i_it, rank, errpr))
+    # !!!
+    # incrementally build up QB decomposition
+    # !!!
+    while ( (errpr>err_thr) and (i_it<maxit) and (rank<=min(A.shape)-delta_t) ): #
+        G = np.random.randn(n, delta_t)
+        Y = np.dot(A,G)
+        #Y = Y - np.dot(Q_old, np.dot(Q_old.T,Y) ) # orthogonalization with Q
+        # power scheme
+        for j in range(npow):
+            Q = QR_fact(Y)
+            Q = QR_fact( np.dot(A.T,Q) )
+            Y = np.dot(A,Q)
+        Y  = Y - np.dot(Q_old, np.dot(Q_old.T,Y) ) # orthogonalization with Q
+        Q_new = QR_fact(Y)
+        B_new = np.dot(Q_new.T,A)
+        # build up approximate basis
+        Q_old  = np.append(Q_new, Q_old, axis=1)
+        #B      = np.append(B_new, B, axis=0)
+        normB += np.linalg.norm(B_new, ord='fro')**2
+        errpr  = abs( normA - normB ) / normA
+        rank  += delta_t
+        i_it  += 1
+        print("iteration = {}, rank = {}, error = {}".format(i_it, rank, errpr))
+    # !!!
+    #UL, SL, VLT = np.linalg.svd(B, full_matrices=0)
+    #UL = np.dot(Q_old,UL)
+    # !!!
+    print("iteration = {}, rank = {}, error = {}".format(i_it, rank, errpr))
+    UL, SL, VLT = powit_RSVD(A, rank, npow, nb_oversamp)
+    #return UL, SL, VLT
+    # !!!
+    rank  = SL.shape[0]
+    new_r = rank
+    for i in range(rank):
+        if( SL[i] <= tol ):
+            new_r = i
+            break
+    return UL[:,:(new_r)], SL[:(new_r)], VLT[:(new_r),:]
+# !!!
diff --git a/devel/svdwf/R3SVD_LiYu.py b/devel/svdwf/R3SVD_LiYu.py
new file mode 100644
index 0000000..88221e4
--- /dev/null
+++ b/devel/svdwf/R3SVD_LiYu.py
@@ -0,0 +1,58 @@
+# !!!
+import numpy as np
+from QR import QR_fact
+# !!!
+def R3SVD_LiYu(A, t, delta_t, npow, err_thr, maxit):
+    # !!!
+    # build initial QB decomposition
+    # !!!
+    n = A.shape[1]
+    G = np.random.randn(n, t) # n x t Gaussian random matrix
+    normA  = np.linalg.norm(A, ord='fro')**2
+    i_it = 0
+    rank = 0
+    Y = np.dot(A,G)
+    # The power scheme
+    for j in range(npow):
+        Q = QR_fact(Y)
+        Q = QR_fact( np.dot(A.T,Q) )
+        Y = np.dot(A,Q)
+    # orthogonalization process
+    Q_old = QR_fact(Y)
+    B     = np.dot(Q_old.T,A)
+    normB = np.linalg.norm(B, ord='fro')**2
+    # error percentage
+    errpr = abs( normA - normB ) / normA
+    rank += t
+    i_it += 1
+    print("iteration = {}, rank = {}, error = {}".format(i_it, rank, errpr))
+    # !!!
+    # incrementally build up QB decomposition
+    # !!!
+    while ( (errpr>err_thr) and (i_it<maxit) and (rank<=min(A.shape)-delta_t) ): #
+        G = np.random.randn(n, delta_t) # n x delta_t Gaussian random matrix
+        Y = np.dot(A,G)
+        Y = Y - np.dot(Q_old, np.dot(Q_old.T,Y) ) # orthogonalization with Q
+        # power scheme
+        for j in range(npow):
+            Q = QR_fact(Y)
+            Q = QR_fact( np.dot(A.T,Q) )
+            Y = np.dot(A,Q)
+        Y  = Y - np.dot(Q_old, np.dot(Q_old.T,Y) ) # orthogonalization with Q
+        Q_new = QR_fact(Y)
+        B_new = np.dot(Q_new.T,A)
+        # build up approximate basis
+        Q_old  = np.append(Q_new, Q_old, axis=1)
+        B      = np.append(B_new, B, axis=0)
+        rank  += delta_t
+        i_it  += 1
+        normB += np.linalg.norm(B_new, ord='fro')**2
+        errpr  = abs( normA - normB ) / normA
+        print("iteration = {}, rank = {}, error = {}".format(i_it, rank, errpr))
+    # !!!
+    print("iteration = {}, rank = {}, error = {}".format(i_it, rank, errpr))
+    UL, SL, VLT = np.linalg.svd(B, full_matrices=0)
+    UL = np.dot(Q_old,UL)
+    # !!!
+    return UL, SL, VLT
+# !!!
diff --git a/devel/svdwf/RSVD.irp.f b/devel/svdwf/RSVD.irp.f
new file mode 100644
index 0000000..f25ad1f
--- /dev/null
+++ b/devel/svdwf/RSVD.irp.f
@@ -0,0 +1,114 @@
+! !!!
+subroutine RSVD( i_state, low_rank, PowerIt_max, nb_oversamp, URSVD, DRSVD, VtRSVD )
+  ! !!!
+  BEGIN_DOC
+  ! standard RSVD for a prefixed rank
+  END_DOC
+  ! !!!
+  implicit none
+  include 'constants.include.F'
+  ! !!!
+  integer, intent(in)           :: i_state, low_rank, PowerIt_max, nb_oversamp
+  double precision, intent(out) :: URSVD(n_det_alpha_unique,low_rank), DRSVD(low_rank), VtRSVD(low_rank,n_det_beta_unique)
+  ! !!!
+  integer                       :: i, j, k, l, PowerIt, m, n 
+  double precision, allocatable :: r1(:,:), Q(:,:), P(:,:), B(:,:)
+  double precision, allocatable :: UB(:,:), D(:), Vt(:,:), U(:,:)
+  ! !!!
+  m = n_det_alpha_unique
+  n = n_det_beta_unique
+  ! !!!
+  ! !!! !!! !!! !!! !!! !!! !!! !!! !!! !!! !
+  ! !!!
+  allocate( Q(m, low_rank+nb_oversamp) )
+  !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,k,l,r1)
+  allocate( r1(N_det,2) )
+  !$OMP DO
+  do l = 1, low_rank+nb_oversamp
+    Q(:,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(i,l) = Q(i,l) + psi_bilinear_matrix_values(k,i_state) * r1(k,1)
+    enddo
+  enddo
+  !$OMP END DO
+  deallocate(r1)
+  !$OMP END PARALLEL
+  ! !!!
+  call ortho_qr(Q, size(Q,1), m, low_rank+nb_oversamp)
+  ! !!!
+  ! power scheme
+  ! !!!
+  allocate( P(n, low_rank+nb_oversamp) ) 
+  do PowerIt = 1, PowerIt_max
+    ! !!!
+    !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l)
+    !$OMP DO
+    do l = 1, low_rank+nb_oversamp
+      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(i,l)
+      enddo
+    enddo
+    !$OMP END DO
+    !$OMP END PARALLEL
+    ! !!!
+    call ortho_qr(P, size(P,1), n, low_rank+nb_oversamp)
+    ! !!!
+    !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l)
+    !$OMP DO
+    do l = 1, low_rank+nb_oversamp
+      Q(:,l) = 0.d0
+      do k = 1 , N_det
+        i      = psi_bilinear_matrix_rows(k)
+        j      = psi_bilinear_matrix_columns(k)
+        Q(i,l) = Q(i,l) + psi_bilinear_matrix_values(k,i_state) * P(j,l)
+      enddo
+    enddo
+    !$OMP END DO
+    !$OMP END PARALLEL
+    ! !!!
+    call ortho_qr(Q, size(Q,1), m, low_rank+nb_oversamp)
+    ! !!!
+  enddo
+  deallocate(P)
+  ! !!!
+  allocate( B(low_rank+nb_oversamp,n) )
+  !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l)
+  !$OMP DO
+  do l = 1, low_rank+nb_oversamp
+    B(l,:) = 0.d0
+    do k = 1, N_det
+      i      = psi_bilinear_matrix_rows(k)
+      j      = psi_bilinear_matrix_columns(k)
+      B(l,j) = B(l,j) + psi_bilinear_matrix_values(k,i_state) * Q(i,l)
+    enddo
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+  ! !!!
+  allocate( UB(low_rank+nb_oversamp,low_rank+nb_oversamp), D(low_rank+nb_oversamp), Vt(low_rank+nb_oversamp,n) )
+  call svd_s( B, size(B,1), UB, size(UB,1), D, Vt, size(Vt,1), low_rank+nb_oversamp, n)
+  deallocate(B)
+  allocate( U(m,low_rank+nb_oversamp) )
+  call dgemm('N', 'N', m, low_rank+nb_oversamp, low_rank+nb_oversamp, 1.d0, Q, size(Q,1), UB, size(UB,1), 0.d0, U, size(U,1))
+  deallocate( Q,UB )
+  ! !!!
+  do l = 1, low_rank
+    DRSVD(l) = D(l)
+    do i = 1, m
+      URSVD(i,l) = U(i,l)
+    enddo
+    do i = 1, n
+      VtRSVD(l,i) = Vt(l,i)
+    enddo
+  enddo
+  ! !!!
+  return
+  ! !!!
+end
diff --git a/devel/svdwf/RSVD.py b/devel/svdwf/RSVD.py
new file mode 100644
index 0000000..4f59f39
--- /dev/null
+++ b/devel/svdwf/RSVD.py
@@ -0,0 +1,20 @@
+# !!!
+import numpy as np
+from QR import QR_fact
+# !!!
+def powit_RSVD(X, new_r, nb_powit, nb_oversamp):
+    # !!!
+    G  = np.random.randn(X.shape[1], new_r+nb_oversamp)
+    Q  = QR_fact( np.dot(X,G) )
+    # !!!
+    for _ in range(nb_powit):
+        Q = QR_fact( np.dot(X.T,Q) )
+        Q = QR_fact( np.dot(X,Q) )
+        # !!!
+    Y = np.dot(Q.T,X)
+    # !!!
+    U, S, VT = np.linalg.svd(Y, full_matrices=0)
+    U = np.dot(Q,U)
+    return U[:,:(new_r)], S[:(new_r)], VT[:(new_r),:]
+    # !!!
+# !!!
\ No newline at end of file
diff --git a/devel/svdwf/linear_algebra.irp.f b/devel/svdwf/linear_algebra.irp.f
new file mode 100644
index 0000000..ed675f9
--- /dev/null
+++ b/devel/svdwf/linear_algebra.irp.f
@@ -0,0 +1,127 @@
+subroutine svd_s(A, LDA, U, LDU, D, Vt, LDVt, m, n)
+  implicit none
+  BEGIN_DOC
+  ! !!!
+  ! DGESVD computes the singular value decomposition (SVD) of a real
+  ! M-by-N matrix A, optionally computing the left and/or right singular
+  ! vectors. The SVD is written:
+  !      A = U * SIGMA * transpose(V)
+  ! where SIGMA is an M-by-N matrix which is zero except for its
+  ! min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
+  ! V is an N-by-N orthogonal matrix.  The diagonal elements of SIGMA
+  ! are the singular values of A; they are real and non-negative, and
+  ! are returned in descending order.  The first min(m,n) columns of
+  ! U and V are the left and right singular vectors of A.
+  !
+  ! Note that the routine returns V**T, not V.
+  ! !!!
+  END_DOC
+
+  integer, intent(in)             :: LDA, LDU, LDVt, m, n
+  double precision, intent(in)    :: A(LDA,n)
+  double precision, intent(out)   :: U(LDU,min(m,n)), Vt(LDVt,n), D(min(m,n))
+  double precision,allocatable    :: work(:), A_tmp(:,:)
+  integer                         :: info, lwork, i, j, k
+
+
+   allocate (A_tmp(LDA,n))
+   do k=1,n
+     do i=1,m
+       A_tmp(i,k) = A(i,k) + 1d-16
+     enddo
+   enddo
+ 
+   ! Find optimal size for temp arrays
+   allocate(work(1))
+   lwork = -1
+   call dgesvd('A', 'S', m, n, A_tmp, LDA,                             &
+       D, U, LDU, Vt, LDVt, work, lwork, info)
+   ! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648
+   lwork = max(int(work(1)), 5*MIN(M,N))
+   deallocate(work)
+
+   allocate(work(lwork))
+   call dgesvd('A','S', m, n, A_tmp, LDA,                             &
+       D, U, LDU, Vt, LDVt, work, lwork, info)
+   deallocate(A_tmp,work)
+
+   if (info /= 0) then
+     print *,  info, ': SVD in svds failed'
+     stop
+   endif
+ 
+   do j=1,min(m,n)
+     do i=1,m
+       if (dabs(U(i,j)) < 1.d-14)  U(i,j) = 0.d0
+     enddo
+   enddo
+ 
+   do j = 1, n
+     do i = 1, LDVt
+       if (dabs(Vt(i,j)) < 1.d-14) Vt(i,j) = 0.d0
+     enddo
+   enddo
+ 
+end
+
+
+
+
+
+
+subroutine my_ortho_qr(A,LDA,m,n)
+  implicit none
+  BEGIN_DOC
+  ! Orthogonalization using Q.R factorization
+  !
+  ! A : matrix to orthogonalize
+  !
+  ! LDA : leftmost dimension of A
+  !
+  ! m : Number of rows of A
+  !
+  ! n : Number of columns of A
+  !
+  END_DOC
+  integer, intent(in)             :: m, n, LDA
+  double precision, intent(inout) :: A(LDA,n)
+  integer                         :: LWORK, INFO, nTAU, ii, jj
+  double precision, allocatable   :: TAU(:), WORK(:)
+  double precision                :: Adorgqr(LDA,n)
+
+  allocate (TAU(min(m,n)), WORK(1))
+
+  LWORK=-1
+  call dgeqrf( m, n, A, LDA, TAU, WORK, LWORK, INFO )
+  ! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648
+  LWORK=max(n,int(WORK(1)))
+
+  deallocate(WORK)
+  allocate(WORK(LWORK))
+  call dgeqrf(m, n, A, LDA, TAU, WORK, LWORK, INFO )
+  if(INFO.ne.0 ) then
+   print*, 'problem in DGEQRF'
+  endif
+
+  nTAU = size(TAU)
+  do jj = 1, n
+    do ii = 1, LDA
+      Adorgqr(ii,jj) = A(ii,jj)
+    enddo
+  enddo
+
+  LWORK=-1
+  call dorgqr(m, n, nTAU, Adorgqr, LDA, TAU, WORK, LWORK, INFO)
+  ! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648
+  LWORK=max(n,int(WORK(1)))
+
+  deallocate(WORK)
+  allocate(WORK(LWORK))
+  call dorgqr(m, n, nTAU, A, LDA, TAU, WORK, LWORK, INFO)
+  if(INFO.ne.0 ) then
+   print*, 'problem in DORGQR'
+  endif
+
+
+  deallocate(WORK,TAU)
+end
diff --git a/devel/svdwf/pyth_RSVD.py b/devel/svdwf/pyth_RSVD.py
new file mode 100644
index 0000000..3729a33
--- /dev/null
+++ b/devel/svdwf/pyth_RSVD.py
@@ -0,0 +1,123 @@
+#!/usr/bin/env python3
+# !!!
+import os, sys
+# !!!
+#QP_PATH=os.environ["QMCCHEM_PATH"]
+#sys.path.insert(0,QMCCHEM_PATH+"/EZFIO/Python/")
+# !!!
+from ezfio import ezfio
+from datetime import datetime
+import numpy as np
+from scipy.sparse.linalg import svds
+from R3SVD_LiYu  	       import R3SVD_LiYu
+from RSVD                import powit_RSVD  
+from R3SVD_AMMAR         import R3SVD_AMMAR
+import time 
+# !!!
+fmt = '%5d' + 2 * ' %e'
+# !!!
+if __name__ == "__main__":
+    # !!!
+    if len(sys.argv) != 2:
+        print("Usage: %s <EZFIO_DIRECTORY>"%sys.argv[0])
+        sys.exit(1)
+    filename = sys.argv[1]
+    ezfio.set_file(filename)
+    # !!!
+    N_det  = ezfio.get_spindeterminants_n_det()
+    A_rows = np.array(ezfio.get_spindeterminants_psi_coef_matrix_rows())
+    A_cols = np.array(ezfio.get_spindeterminants_psi_coef_matrix_columns())
+    A_vals = np.array(ezfio.get_spindeterminants_psi_coef_matrix_values())
+    nrows, ncols  = ezfio.get_spindeterminants_n_det_alpha(), ezfio.get_spindeterminants_n_det_beta()
+    Y = np.zeros( (nrows, ncols) )
+    for k in range(N_det):
+        i = A_rows[k] - 1
+        j = A_cols[k] - 1
+        Y[i,j] = A_vals[0][k]
+    print("# # # # # # # # # # # # # # # # # # # # # #")
+    print('matrix dimensions = {} x {}'.format(nrows, ncols))
+    print("# # # # # # # # # # # # # # # # # # # # # # \n")
+    normY = np.linalg.norm(Y, ord='fro')
+    print( normY )
+    # !!!
+    print('Full SVD:')
+    t_beg = time.time()
+    U, S_FSVD, VT = np.linalg.svd(Y, full_matrices=0)
+    t_end         = time.time()
+    rank          = S_FSVD.shape[0]
+    energy        = np.sum(np.square(S_FSVD)) / normY**2
+    err_SVD       = 100. * np.linalg.norm(Y - np.dot(U,np.dot(np.diag(S_FSVD),VT)), ord='fro') / normY
+    print('rank = {}, energy = {}, error = {}%, CPU time = {} \n'.format(rank, energy, err_SVD, t_end-t_beg))
+    # !!!
+    np.savetxt('results_python/h2o_pseudo/SingValues_FullSVD.txt', np.transpose([ np.array(range(rank))+1, S_FSVD  ]), fmt='%5d' + ' %e', delimiter=' ')
+    # !!!
+    t       = 50
+    delta_t = 10
+    npow    = 15
+    err_thr = 1e-3
+    maxit   = 10
+    # !!!
+    print('RRR SVD Li & Yu:')
+    t_beg          = time.time()
+    U, S_R3SVD, VT = R3SVD_LiYu(Y, t, delta_t, npow, err_thr, maxit)
+    t_end          = time.time()
+    rank           = S_R3SVD.shape[0]
+    energy         = np.sum( np.square(S_R3SVD) ) / normY**2
+    err_SVD        = 100. * np.linalg.norm(Y - np.dot(U,np.dot(np.diag(S_R3SVD),VT)), ord='fro') / normY
+    print('nb Pow It = {}, rank = {}, energy = {}, error = {} %, CPU time = {}\n'.format(npow, rank, energy, err_SVD, t_end-t_beg))
+    # !!!
+    err_R3SVD = np.zeros( (rank) )
+    for i in range(rank):
+        err_R3SVD[i] = 100.0 * abs( S_FSVD[i] - S_R3SVD[i] ) / S_FSVD[i]
+    np.savetxt('results_python/h2o_pseudo/SingValues_R3SVD.txt', np.transpose([ np.array(range(rank))+1, S_R3SVD, err_R3SVD ]), fmt=fmt, delimiter=' ')
+    # !!!
+    nb_oversamp = 10
+    tol_SVD     = 1e-10
+    print('RRR SVD my version:')
+    t_beg          = time.time()
+    U, S_MRSVD, VT = R3SVD_AMMAR(Y, t, delta_t, npow, nb_oversamp, err_thr, maxit, tol_SVD)
+    t_end          = time.time()
+    rank           = S_MRSVD.shape[0]
+    energy         = np.sum( np.square(S_MRSVD) ) / normY**2
+    err_SVD        = 100. * np.linalg.norm(Y - np.dot(U,np.dot(np.diag(S_MRSVD),VT)), ord='fro') / normY
+    print('nb Pow It = {}, rank = {}, energy = {}, error = {} %, CPU time = {}\n'.format(npow, rank, energy, err_SVD, t_end-t_beg))
+    # !!!
+    err_MRSVD = np.zeros( (rank) )
+    for i in range(rank):
+        err_MRSVD[i] = 100.0 * abs( S_FSVD[i] - S_MRSVD[i] ) / S_FSVD[i]
+    np.savetxt('results_python/h2o_pseudo/SingValues_MRSVD.txt', np.transpose([ np.array(range(rank))+1, S_MRSVD, err_MRSVD ]), fmt=fmt, delimiter=' ')
+    # !!!
+    trank = rank
+    print("Truncated RSVD (pre-fixed rank = {} & oversampling parameter = {}):".format(trank,nb_oversamp))
+    t_beg         = time.time()
+    U, S_RSVD, VT = powit_RSVD(Y, trank, npow, nb_oversamp)
+    t_end         = time.time()
+    rank          = S_RSVD.shape[0]
+    energy        = np.sum( np.square(S_RSVD) ) / normY**2
+    err_SVD       = 100. * np.linalg.norm( Y - np.dot(U,np.dot(np.diag(S_RSVD),VT)), ord="fro") / normY
+    print('nb Pow It = {}, rank = {}, energy = {}, error = {} %, CPU time = {}\n'.format(npow, rank, energy, err_SVD, t_end-t_beg))
+    # !!!
+    err_RSVD = np.zeros( (rank) )
+    for i in range(rank):
+        err_RSVD[i] = 100.0 * abs( S_FSVD[i] - S_RSVD[i] ) / S_FSVD[i]
+    np.savetxt('results_python/h2o_pseudo/SingValues_RSVD.txt', np.transpose([ np.array(range(rank))+1, S_RSVD, err_RSVD ]), fmt=fmt, delimiter=' ')
+    # !!!
+    print("Truncated SVD (scipy):")
+    t_beg         = time.time()
+    U, S_TSVD, VT = svds(Y, min(trank, min(Y.shape[0],Y.shape[1])-1 ), which='LM')
+    t_end         = time.time()
+    rank          = S_TSVD.shape[0]
+    energy        = np.sum( np.square(S_TSVD) ) / normY**2
+    err_SVD       = 100. * np.linalg.norm( Y - np.dot(U, np.dot(np.diag(S_TSVD),VT) ), ord="fro") / normY
+    print('rank = {}, energy = {}, error = {} %, CPU time = {}\n'.format(rank, energy, err_SVD, t_end-t_beg))
+    # !!!
+    err_TSVD = np.zeros( (rank) )
+    for i in range(rank-1):
+        for j in range(i+1,rank):
+            if( S_TSVD[j] > S_TSVD[i]):
+                S_TSVD[i], S_TSVD[j] = S_TSVD[j], S_TSVD[i]
+    for i in range(rank):
+        err_TSVD[i] = 100.0 * abs( S_FSVD[i] - S_TSVD[i] ) / S_FSVD[i]
+    np.savetxt('results_python/h2o_pseudo/SingValues_TSVD.txt', np.transpose([ np.array(range(rank))+1, S_TSVD, err_TSVD ]), fmt=fmt, delimiter=' ')
+    # !!!
+# !!!