function [f,D,A,E,Errs] = LRSR_RR_A(X,Y,eta,lambda,rho,gamma)
itrN = 10;
errX = zeros(itrN,1);
errD = zeros(itrN,1);
errZ = zeros(itrN,1);
[d_x,n] = size(X);
d_y = size(Y,1);
I = eye(d_x+1);
Z = eye(n);
J = zeros(n,n);
A = X;
D = A*Z;
Dhat = [D;ones(1,n)];
T = Y * Dhat'*inv(Dhat*Dhat'+ gamma*eye(d_x+1));
E = X-D;
Y1 = X/norm(X);
Y2 = Dhat/norm(Dhat);
Y3 = Z/norm(Z);
mu1 = d_x*n/4/norm(X(:),1);
mu2 = d_x*n/4/norm(D(:),1);
mu3 = 100*d_x*n/4/norm(D(:),1);
lambda1 = 1 / sqrt(min(d_x,n));
lambda2 = 0.01 / sqrt(min(d_x,n));
s = svd(X);
nunormX = sum(s);
eta = nunormX / norm(Y-T*Dhat,'fro');
for iter = 1:itrN
    disp([num2str(iter),'...']);
    %LRSR-1: T
    T = (Y * Dhat') / (Dhat*Dhat'+ gamma*eye(d_x+1));
    %LRSR-1: Dhat
    Dhat = (eta*(T'*T)+mu2*I)\(eta*T'*Y-Y2+mu2*[D;ones(1,n)]);
    %LRSR-3: A
    %tau_A1 = vrho(Z*Z') + 0.000001;
    tau_A = norm(Z*Z',2);
    beta_A = (mu1 + mu2) * tau_A / 1;
    FAk = ( -Y1 + mu1 * ( D - (X-E) ) - Y2(1:d_x,:) + mu2 * ( D - Dhat(1:d_x,:) ) ) * Z';
    [U,S,V] = svd(A - FAk / beta_A, 'econ');
    s = diag(S);
    s = max(s-0.5/beta_A,0);
    A = U(:,1:length(s))*diag(s)*V(:,1:length(s))';
    %A = ( Y1 + mu1 * (X-E)  + Y2(1:d_x,:) + mu2 * Dhat(1:d_x,:)  ) * Z' * inv(Z*Z')/(mu1 + mu2);
    D = A*Z;
    %LRSR-4: Z--------------------------------------
    %tau_Z1= vrho(A'*A) + 0.000001;
    %beta_Z1 = (mu1 + mu2 + mu3) * tau_Z1 / 1;
    beta_Z = norm((mu1+mu2)*(A'*A)+mu3*eye(n),2);
    FZk = A'*( -Y1 + mu1 * ( D - (X-E) ) - Y2(1:d_x,:) + mu2 * ( D - Dhat(1:d_x,:) ) ) ...
          + mu3 * (Z - J) + Y3;
    [U,S,V] = svd(Z - FZk / beta_Z, 'econ');
    s = diag(S);
    s = max(s-1/beta_Z,0);
    Z = U(:,1:length(s))*diag(s)*V(:,1:length(s))';
    D = A*Z;
    %LRSR-5: J
    J = max(Z + Y3/mu3 - lambda2/mu3,0);
    %LRSR-6: E
    TempE =  X - D + Y1 / mu1; 
    for cn=1:n
        Ecol = TempE(:,cn);
        if norm(Ecol) < lambda1 / mu1
            Ecol(:) = 0;
        end
        E(:,cn) = Ecol;
    end
    %E1 = So( lambda1 / mu1 , X - D + Y1 / mu1 );
    %LRSR-7: Lag
    Y1 = Y1 + mu1*(X-D-E);
    Y2 = Y2 + mu2*(Dhat-[D;ones(1,n)]);
    Y3 = Y3 + mu3*(Z-J);
    mu1 = min(rho*mu1,1e6);
    mu2 = min(rho*mu2,1e6);
    mu3 = min(rho*mu3,1e6);
    %eta = nunormX / norm(Y-T*Dhat,'fro')^2;
    %LRSR-8: Error Ratios
    Xdif_ratio = norm(X-D-E,'fro')/norm(X,'fro');
    Ddif_ratio = norm(Dhat-[D;ones(1,n)],'fro')/norm(Dhat,'fro');
    Zdif_ratio = norm(Z-J,'fro')/norm(Z,'fro');
    errX(iter) = Xdif_ratio;
    errD(iter) = Ddif_ratio;
    errZ(iter) = Zdif_ratio;
    if  Xdif_ratio <=1e-8 && Ddif_ratio<=1e-8
        break;
    end
end
f = T;
Errs = {errX,errD,errZ};
end
function r = So(tau, X)
    % shrinkage operator
    r = sign(X) .* max(abs(X) - tau, 0);
end

function r = Do(tau, X)
    % shrinkage operator for singular values
    [U, S, V] = svd(X, 'econ');
    r = U*So(tau, S)*V';
end