# !!!
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
# !!!