function [Z,Ez,Ew,W,Wi] = maLRR(T,S,Dim,Maxiter,alaph,beta)
% maLRR--Multi-Center Adaptation Framework with Low-Rank Representation
% Input:
%       T: Target domain data
%       S: Source domain data
% Output:
%       Z: Low rank representation of the source domain in target domain
%       Ez: error matrix of representation
%       Ew: error matrix of project
%       W: the project matrix of different source data
%       Ws: the shared project matrix
%  The code requires a GPU to run

%The parameters need to be tuned
mu = 1e-5;tol = 1e-8;max_mu = 1e6;rho =1.2;iter = 0;
[d,m] = size(T);M = length(S);Z = cell(1,M);Wi = cell(1,M);
W = eye(Dim,d);Ew = cell(1,M);Ez = cell(1,M);F = cell(1,M);Y1 = cell(1,M);
Y2 = cell(1,M);Y3 = zeros(Dim,d);Y4 = cell(1,M);Y1_tmp = cell(1,M);
Y2_tmp = cell(1,M);Y3_tmp =  zeros(Dim,d);Y4_tmp = cell(1,M);stopC = zeros(1,M);J = zeros(d);
for i=1:M
    Z{1,i} = zeros(m,size(S{1,i},2));
    Wi{1,i} = zeros(Dim,d);Ew{1,i} = zeros(Dim,d);
    Ez{1,i} = zeros(Dim,size(S{1,i},2));
    F{1,i} = zeros(m,size(S{1,i},2));
    Y1{1,i} = zeros(Dim,size(S{1,i},2));
    Y2{1,i} = zeros(Dim,d);
    Y4{1,i} = zeros(m,size(S{1,i},2));
    Y1_tmp{1,i} =zeros(Dim,size(S{1,i},2));
    Y2_tmp{1,i} = zeros(Dim,d);
    Y4_tmp{1,i} = zeros(m,size(S{1,i},2));
end
while iter < Maxiter
    iter
    iter = iter+1;
    for i=1:M
        temp = Z{1,i} + Y4{1,i}/mu;
        gpu_tmp = gpuArray(temp);
        [U_Temp,sigma_Temp,V_Temp] = svd(gpu_tmp,'econ');
        gpu_U_Temp = gather(U_Temp);
        gpu_sigma = gather(sigma_Temp);
        gpu_V_Temp = gather(V_Temp);
        sigma = diag(gpu_sigma);
        svp = length(find(sigma>1/mu));
        if svp>=1
            sigma = sigma(1:svp)-1/mu;
        else
            svp = 1;
            sigma = 0;
        end
        F{1,i} = gpu_U_Temp(:,1:svp)*diag(sigma)*gpu_V_Temp(:,1:svp)';
    end
    for i=1:M
        Wpt1 = S{1,i}*S{1,i}'+eye(d);
        Wpt2 = (W*T*Z{1,i}+Ez{1,i})*S{1,i}'+W+Ew{1,i}-(Y1{1,i}*S{1,i}'+Y2{1,i})/mu;
        temp = Wpt2/ Wpt1;
        Wi{1,i} = temp;
    end
    for i=1:M
        Zpt1 = T'*(W'*W)*T+eye(m);
        Zpt2 = (T'*W'*Y1{1,i}-Y4{1,i})/mu + F{1,i}+T'*W'*(Wi{1,i}*S{1,i}-Ez{1,i});
        temp = Zpt1\Zpt2;
        Z{1,i} = temp;
    end
    temp = W + Y3/mu;
    gpu_tmp = gpuArray(temp);
    [U_Temp,sigma_Temp,V_Temp] = svd(gpu_tmp,'econ');
    gpu_U_Temp = gather(U_Temp);
    gpu_sigma = gather(sigma_Temp);
    gpu_V_Temp = gather(V_Temp);
    sigma = diag(gpu_sigma);
    svp = length(find(sigma>1/mu));
    if svp>=1
        sigma = sigma(1:svp)-1/mu;
    else
        svp = 1;
        sigma = 0;
    end
    J = gpu_U_Temp(:,1:svp)*diag(sigma)*gpu_V_Temp(:,1:svp)';
    for i=1:M
        temp = Wi{1,i}*S{1,i}-W*T*Z{1,i}+ Y1{1,i}/mu;
        Ez{1,i} = max(0,temp - alaph/mu)+min(0,temp + alaph/mu);
    end
    for i=1:M
        temp = Wi{1,i}-W+Y2{1,i}/mu;
        Ew{1,i} = max(0,temp - beta/mu)+min(0,temp + beta/mu);
    end
    temp = zeros(d);
    for i=1:M
        temp = temp + T*Z{1,i}*Z{1,i}'*T'+eye(d);
    end
    Wpt1 = temp + eye(d);
    temp = zeros(Dim,d);
    for i=1:M
        temp = temp + Y1{1,i}*Z{1,i}'*T'+Y2{1,i};
    end
    Wpt2 = (temp+ mu*J-Y3)/mu;
    temp = zeros(Dim,d);
    for i=1:M
        temp = temp + (Wi{1,i}*S{1,i}-Ez{1,i})*Z{1,i}'*T'+Wi{1,i}-Ew{1,i};
    end
    Wpt2 = Wpt2 + temp;
    W = Wpt2/Wpt1;
    W = orth(W);
    Y3_tmp = W-J;
    for i=1:M
        Y1_tmp{1,i} = Wi{1,i}*S{1,i}-W*T*Z{1,i}-Ez{1,i};
        Y2_tmp{1,i} = Wi{1,i}-W-Ew{1,i};
        Y4_tmp{1,i} = Z{1,i}-F{1,i};
        stopC(1,i) = max(max(max(abs(Y1_tmp{1,i}))),max(max(abs(Y2_tmp{1,i}))));
        stopC(1,i) = max(max(max(abs(Y3_tmp))),stopC(1,i));
        stopC(1,i) = max(max(max(abs(Y4_tmp{1,i}))),stopC(1,i));
    end
    if iter==1 || mod(iter,50)==0 || max(stopC)<tol
        disp(['iter ' num2str(iter) ',mu=' num2str(mu,'%2.1e') ...
            ',stopALM=' num2str(max(stopC),'%2.3e')]);
    end
    if max(stopC)<tol
        break;
    else
        Y3 = Y3 + mu*Y3_tmp;
        for i=1:M
            Y1{1,i} = Y1{1,i} + mu*Y1_tmp{1,i};
            Y2{1,i} = Y2{1,i} + mu*Y2_tmp{1,i};
            Y4{1,i} = Y4{1,i} + mu*Y4_tmp{1,i};
        end
        mu = min(max_mu,mu*rho);
    end
end