function out = ALM_Mat_Comple_v3(b, Omega, opts)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% b = P_Omega ( M ) is a vector 
% Omega is a vector of index 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% m = length(Omega); 
n1 = opts.n1; n2 = opts.n2; mn1n2 = min(n1,n2);
M = zeros(n1,n2); M(Omega) = b;

mu = opts.mu; sigma = opts.sigma; rho = opts.rho;
Y = zeros(n1,n2); gradgY = zeros(n1,n2);
X = M - Y; Mnorm = norm(M,'fro'); 
Lambda = zeros(n1,n2); sv = opts.sv; unchange_mu = 0; StopCrit = 1;
mux = mu; muy = mu; 

for itr = 1: opts.maxitr
%    if choosvd(mn1n2, sv) == 1
%        [U, gamma, V] = lansvd(mux*gradgY-Y+M,sv,'L');
%    else
        [U, gamma, V] = svd(mux*gradgY-Y+M, 'econ');
%    end
    gamma = diag(gamma);
    gamma_new = gamma-mux*gamma./max(gamma,mux+sigma);
    svp = length(find(gamma > mux));
    if svp < sv
        sv = min(svp + round(0.04*mn1n2),mn1n2);
    elseif svp > sv
        sv = max(svp - round(0.02*mn1n2),10);
    else
        sv = min(svp,mn1n2);
    end
    sv = mn1n2;
    
    Xp = X;
    X = U*diag(gamma_new)*V';
    gradfX = U*diag(min(gamma_new./sigma,1))*V';

    Yp = Y;
    Y = solveSubY(muy,rho,Omega,X,M,gradfX,sigma);
    Yvec = Y(Omega);
    gradgY = zeros(n1,n2); gradgY(Omega) = min(rho,max(-rho,Yvec/sigma));
    
    critp = StopCrit;
    StopCrit = norm(M-X-Y,'fro')/Mnorm;
    critchange = abs(StopCrit - critp)/max([critp,StopCrit,1]);
    relX = norm(X-opts.Xs,'fro')/norm(opts.Xs,'fro');
    diffX = norm(Xp-X,'fro')/max([norm(Xp,'fro'),norm(X,'fro'),1]);
    diffY = norm(Yp-Y,'fro')/max([norm(Yp,'fro'),norm(Y,'fro'),1]);

    RelChg = norm([X-Xp,Y-Yp],'fro')/(norm([Xp,Yp],'fro')+1);
    
    tmpY = Y; tmpY(opts.Omega_c) = 0;
    relY = norm(tmpY-opts.Ys,'fro')/norm(opts.Ys,'fro');
    fprintf('iter: %d, mu: %3.2e, rank(X):%d, RelChg: %3.2e,relX: %3.2e, relY: %3.2e,crit:%3.2e, diffX:%3.2e,diffY:%3.2e\n', ...
        itr, mux, sv, RelChg, relX, relY, norm(M-X-Y,'fro')/Mnorm,diffX,diffY);

    unchange_mu = unchange_mu + 1;
%     if unchange_mu >=3
%     if critchange < 1e-8
%     if relX > 1e-1
%     if itr < 30
%         mux = max(opts.muf, mux*opts.eta_mu);
%         muy = max(opts.muf, muy*opts.eta_mu);
%     if unchange_mu >= 10
       mux = max(opts.muf, mux*opts.eta_mu);
        muy = max(opts.muf, muy*opts.eta_mu);
%         mux = max(mux,10); muy = max(muy,10);
%         mux = max(mux,1); muy = max(muy,1);
%         mu = max(opts.muf, mu*opts.eta_mu);
%         mu = max(mu,500);
%         mu = max(opts.muf, mu*0.1);
%         sigma = max(opts.sigmaf, sigma*opts.eta_sigma);
%         sigma = mu*1e-3;
%         mu = max(1e+1, mu);
        unchange_mu = 0;
%     end
    
    % check stop
    if StopCrit < opts.epsilon
%     if RelChg < 1e-6
        out.X = X; out.iter = itr; out.relX = relX; out.relY = relY; out.StopCrit = StopCrit;
        tmpY = Y; tmpY(opts.Omega_c) = 0; out.Y = tmpY;
        return;
    end
end
out.X = X; out.iter = itr; out.relX = relX; out.relY = relY; out.StopCrit = StopCrit;
tmpY = Y; tmpY(opts.Omega_c) = 0; out.Y = tmpY;

function D = solveSubY(mu,rho,Omega,X,M,Lambda,sigma)
D = M + mu*Lambda - X; 
tmp = D(Omega); 
% tmp = sign(tmp).*max(abs(tmp)-mu*rho,0);
tmp = tmp - mu*min(rho,max(-rho,tmp/(sigma+mu)));
D(Omega) = tmp; 
