function [ result ] = PFR_mul_RunRank( mats, Label, varargin )
%PFR_MUL_RUNRANK construct and solve linear system for PFRank

ip = inputParser;
ip.FunctionName = 'PFR_mul_RunRank';
ip.addRequired('mats', @isstruct);
ip.addRequired('Label', @isstruct);

ip.addOptional('lambda1', 1e-2, @isnumeric);
ip.addOptional('lambda2', 1e-2, @isnumeric);
ip.addOptional('lambda3', 1e-2, @isnumeric); %
ip.addOptional('solver', 'backslash', @ischar);
ip.addOptional('lsqr_thresh', 1e-5, @isnumeric);
ip.addOptional('lsqr_maxiter', 300, @isnumeric);

ip.parse(mats, Label, varargin{:});
par = ip.Results;

[p] = size(Label.y,2);
d = mats.d; n = mats.n;

times = mats.times;
ts_linear_system = cputime;

% -- linear system
II = diag(double(Label.mask));
y  = Label.y; y(~Label.mask,:) = 0;

lsA = [par.lambda1 * mats.G + par.lambda2 * mats.B + par.lambda3 * mats.D, -par.lambda1 * mats.beta;  %
       -par.lambda1 * mats.beta'                                         , II + 2*par.lambda1 * mats.L];
lsb = zeros((d+1)*n, p);

% lsA = zeros((d+1)*n);
% lsA(1:d*n,1:d*n)         = par.lambda1 * mats.G + par.lambda2 * mats.B;
% lsA(d*n+1:end,1:d*n)     = -par.lambda1 * mats.beta';
% lsA(1:d*n,d*n+1:end)     = -par.lambda1 * mats.beta;
% lsA(d*n+1:end,d*n+1:end) = II + 2*par.lambda1 * mats.L;

lsb(d*n+1:end,:)         = y;
lsb(1:d*n, :)            = par.lambda3 * mats.H;   %

times.linear_system = cputime - ts_linear_system;

%************************************************
% solving
%************************************************
ts_solve = cputime;

if strcmp(par.solver, 'linsolve')
    lsX = linsolve(lsA, lsb);
elseif strcmp(par.solver, 'lsqr')
    splsA = sparse(lsA);
    lsX = zeros((d+1)*n,p);
    for i = 1:p
        lsX(:,i) = lsqr(splsA, sparse(lsb(:,i)), par.lsqr_thresh, par.lsqr_maxiter);
    end
elseif strcmp(par.solver, 'backslash')
    lsX = sparse(lsA) \ sparse(lsb);
    lsX = full(lsX);
else
    error('Unknown solver: %s', par.solver);
end

times.solve = cputime - ts_solve;

matU = lsX(1:d*n,:);
matF = lsX(d*n+1:end,:);

result = [];
result.f = matF;
result.u = reshape(matU, d, n, p);
result.v = result.u;
result.t = times;
result.T = mats.T;
result.nv = zeros(d, n);
for i = 1:n
    result.nv(:,i)=result.v(:,i)./norm(result.v(:,i));
end

% % only solve f
% Mtmp = -par.lambda1^2 * alpha * ((par.lambda1*G + par.lambda2*B) \  beta) + II + par.lambda1*L;
% result.f = linsolve(Mtmp,y);

end

