diff --git a/devel/svdwf/Evar_TruncSVD.irp.f b/devel/svdwf/Evar_TruncSVD.irp.f
index ce88d64..1a72dc5 100644
--- a/devel/svdwf/Evar_TruncSVD.irp.f
+++ b/devel/svdwf/Evar_TruncSVD.irp.f
@@ -389,7 +389,7 @@ subroutine run
   ! !!!
   !deallocate( Averif, Uverif, Dverif, Vtverif )
   ! !!!
-  low_rank = 12
+  low_rank = 10
   ! !!!
   err_verif = 0.d0
   do j = 1, n
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.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/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=' ')
+    # !!!
+# !!!