function T=RS_img_Fusion_Fast(HSMS,LSHS,D,alphia,mu,beta,gamma,SRF,maxIter,ratio)

[H0,W0,diml]=size(HSMS);
[Hl,Wl,Dim]=size(LSHS);
%[TH,TW]=size(PSF);
%%%%%padding or cut HSMS
% if H0<Hl*TH
%     paddingH=Hl*TH-H0;
%     HSMS=padarray(HSMS,[paddingH,0],'replicate','post');
% else if H0>Hl*TH
%         HSMS=HSMS(1:Hl*TH,:);
%     end
% end
% 
% if W0<Wl*TW
%     paddingW=Wl*TW-W0;
%     HSMS=padarray(HSMS,[0,paddingW],'replicate','post');
% else if W0>Wl*TW
%         HSMS=HSMS(:,1:Wl*TW);
%     end
% end
[H,W,~]=size(HSMS);
%%%%%initialzation
T=zeros(H,W,Dim);

%D=spatial_samples(PSF,H,W)';%%%%%computing sample matrix D by PSF

%%%%% initalzation T by LSHS
T=imresize(LSHS,ratio);

norm_init_T=norm(T(:),'fro');

M1=Unfold(T,[H,W,Dim],1);
M2=Unfold(T,[H,W,Dim],2);
M3=Unfold(T,[H,W,Dim],3);

HSMS_3=Unfold(HSMS,[H,W,diml],3);%[H0,W0,diml]=size(HSMS);%Tensor_unfold(HSMS,3);
LSHS_3=Unfold(LSHS,[Hl,Wl,Dim],3);
norm_H=norm(HSMS_3,'fro');
norm_L=norm(LSHS_3,'fro');
%S=zeros(Dim,diml);
%%%%% initalzation S
% for i=1:diml
%     S(:,i)=lsqnonneg(M3,HSMS_3(:,i));
% end%M3\HSMS_3;    
%%%%% initalzation laplacian variable
Y1=zeros(H,W*Dim);
Y2=zeros(W,H*Dim);
Y3=zeros(Dim,H*W);
Y_1=zeros(Dim,Hl*Wl);
Y_2=zeros(diml,H*W);

DTD=D*D';
DTL=LSHS_3*D';
STS=SRF'*SRF;
HST=SRF'*HSMS_3;

% normfX = norm(M3,'fro');
tol1 = 1e-4;%threshold for the error in constraint
tol2 = 1e-5;%threshold for the change in the solutions
% opt.tol = tol2;%precision for computing the partial SVD
%maxIter=100;
s1=svds(D,1);
s2=svds(SRF,1);
tau=s1^2*beta+gamma*s2^2+mu;

%% Start main loop

% mu_max=5e-2;
% rho=1.1;
convergenced = 0;
iter = 0;
while iter<maxIter
    iter=iter+1;
    Tk=T;
    %%update M1
    T1=Unfold(T,[H,W,Dim],1);
    
    %[M1,sv1]=svd_threshold(T1+Y1/mu,sv1,alphia(1)/mu,opt);
    [M1, ~, ~] = Pro2TraceNorm(T1+Y1/mu,alphia(1)/mu);
    %%update M2
    T2=Unfold(T,[H,W,Dim],2);
    [M2, ~, ~]=Pro2TraceNorm(T2+Y2/mu,alphia(2)/mu);
    %%update M3
    T3=Unfold(T,[H,W,Dim],3);
    delta_f_M3=-mu*(T3-M3+Y3/mu)+beta*(M3*DTD-DTL+Y_1*D'/beta)+gamma*(STS*M3-HST+SRF'*Y_2/gamma);
    [M3, ~, ~]=Pro2TraceNorm(M3-delta_f_M3/tau,alphia(3)/tau);
    %%update T
    T=(Fold(M1-Y1/mu,[H,W,Dim],1)+Fold(M2-Y2/mu,[H,W,Dim],2)+Fold(M3-Y3/mu,[H,W,Dim],3))/3;    
    %%update S
%     for i=1:diml
%         S(:,i)=lsqnonneg(M3,HSMS_3(:,i)-Y_2(:,i)/gamma);
%     end
    %S=M3\(HSMS_3-Y_2/gamma);
    %%update laplacian variable
    Y1=Y1+mu*(Unfold(T,[H,W,Dim],1)-M1);
    Y2=Y2+mu*(Unfold(T,[H,W,Dim],2)-M2);
    Y3=Y3+mu*(Unfold(T,[H,W,Dim],3)-M3);
    Y_1=Y_1+beta*(M3*D-LSHS_3);
    Y_2=Y_2+gamma*(SRF*M3-HSMS_3);
    
    %mu update
   %mu=min(mu_max,rho*mu);
    
    %%%% convergence condition judgement
    diffT=T-Tk;
    relChgT=sqrt(diffT(:)'*diffT(:))/norm_init_T;
    T3=Unfold(T,[H,W,Dim],3);
    recErr=max(norm(T3*D-LSHS_3,'fro')/norm_L,norm(SRF*T3-HSMS_3,'fro')/norm_H);
    convergenced = recErr <tol1 && relChgT < tol2;
    
%     if iter==1 || mod(iter,10)==0 || convergenced
%         disp(['iter ' num2str(iter) ',relChT=' num2str(relChgT)...
%             ',recErr=' num2str(recErr)]);
%         RMS=norm(T(:)-T0(:),'fro')/norm(T0(:),'fro')
%       % RMS=0;
%     end
    
    if convergenced
        break;
    end
end



end
