function [Y, errList, snrList, timeList] =FCSA_TC(Y, gamma, mark, maxIter, Y0)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FCSA based Tensor complement algorithm implement 
% Written by Junzhou Huang, Rutgers Universtiy
% Jan 31, 2010
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[m, n, d]=size(Y); dim=[m, n,d];
A=find(mark==0); A=A';
b=Y(A); b=b(:);
errList = zeros(maxIter, 1);
Lf=1; Atb = AtMVMx(A,b,m,n,d); 

x=zeros(m, n, d); tnew=1; xr=x;
tt0=cputime;
for k=1:maxIter,
    if mod(k,20)==0
        fprintf('k=%d \n', k);
    end
     xold = x; told=tnew;
     xg = get_g(xr,m,n,d,A,Atb);    
     xg= xr - (1/Lf)*xg;  
     
     x=zeros(m, n, d);
     for i = 1:ndims(Y)
         site = circshift(dim, [1-i, 1-i]);
         xgm=reshape(shiftdim(xg,i-1), dim(i), []);
         [mid, DD(i)] = Pro2TraceNorm(xgm, ndims(Y)*gamma(i)/Lf);
         xi = shiftdim(reshape(mid, site), -i+1+ndims(Y));
         x=x+xi;
     end
     x=x/ndims(Y);

     tnew=(1+sqrt(1+4*told^2))/2;  
     xr=x+(told-1)/tnew*(x-xold);
     
    Y(mark) = x(mark);
    errList(k) = norm(Y(:)-Y0(:), 2)/norm(Y0(:), 2);
    snrList(k)=snr(Y, Y0);
    timeList(k)=cputime-tt0;
end
%--------------------------------------------------------------------------
% SUBFUNCTION FOR CALCULATING g
%--------------------------------------------------------------------------
    function y = AMVMx(A,x)
        % y = A*x
        y = x(A);
    end

    function y = AtMVMx(A,b,m,n,d)
        % y = A'*b
        y = zeros(m*n*d,1);
        y(A) = b;
        y = reshape(y,m,n,d);
    end


    function g = get_g(x,m,n,d,A,Atb)
        Ax = AMVMx(A,reshape(x,m*n*d,1));
        g = AtMVMx(A,Ax,m,n,d)-Atb;
    end

end