function [Xiter, RelErr, RelRes, istop, stopCit, Alpha,  Beta] = fcgls_rp_nn(A, b, x0, m, xex, full, stopC)
%
% setting the stopping criterion options
% allocating memory for the outputs
stopCit = [];
if isstruct(stopC)
    if strcmp(stopC.kind,'stab')
        tau = stopC.tau;
        stopC = 'stab';
        stopCit = [];
        alStop = 0;
    elseif strcmp(stopC.kind,'discr')
        sigma = stopC.noise;
        stopC = 'discr';
        stopCit = [];
        alStop = 0;
    end
end
N = length(x0(:));
Xiter = zeros(N,m);
RelErr = zeros(m,1);
RelRes = zeros(m,1);
istop = m;
Alpha = zeros(m,1);
if full
    Beta = zeros(m);
    Q = zeros(N,m+1);
    AQ = zeros(N,m+1);
else
    Beta = zeros(m,1);
end
if issparse(A) && full
    M = size(A,1);
    AQ = zeros(M,m+1);
else
end
nx = norm(xex(:));
nb = norm(b(:));
%
% initialization
x = x0;
x(x<0)=0;
res = b - A*x; res = res(:);
z = A'*res;
if max(x)~=0
    zb = x(:).*z;
else
    % this means that the initial guess is zero, and we have to perform the
    % first iteration without preconditioning, not to get immediately stucked
    zb = z;
    zb(zb<0) = 0;
end
q = zb;
Azb = A*zb;
Aq = Azb;
if full
    Q(:,1) = q(:);
    AQ(:,1) = Aq(:);
end
% beginning the iterations
for i = 1:m
    w = Aq;
    % FCGLS parameter for the update of the solution
    theta = (res(:)'*w(:))/((w(:))'*w(:));
    % modifying FCGLS to get nonnegativity
    neg_ind = q < 0;
    alpha = min( theta, min( -x(neg_ind)./q(neg_ind) ) );
    if isempty(alpha)
        alpha = theta;
    end
    % detecting stagnation
    if alpha<10^(-15)
        Xiter = Xiter(:,1:i-1); 
        RelErr = RelErr(1:i-1); 
        RelRes = RelRes(1:i-1); 
        Alpha = Alpha(1:i-1);
        istop = i-1;
        Beta = Beta(1:i-1,:);
        return
    end 
    Alpha(i) = alpha;
    x = x + alpha*q;
    res = res - alpha*w;
    % FCGLS parameters for the update of the descent direction
    z = A'*res;
    zb = x(:).*z(:); 
    Azb = A*zb;
    if full
        % full recursion for the update of the descent direction
        beta = zeros(i,1);
        for j=1:i
            beta(j) = -(Azb(:)'*AQ(:,j))/(AQ(:,j)'*AQ(:,j));
        end
        Beta(1:i,i) = beta;
        % truncated recursion (as appropriate) for the update of the descent direction
        if full ~=1
            if full>i
                fulltemp = i;
                q = zb + Q(:,1:fulltemp)*beta(1:fulltemp);
                Aq = Azb + AQ(:,1:fulltemp)*beta(1:fulltemp);
            else
                fulltemp = i-full;
                q = zb + Q(:,fulltemp+1:fulltemp+full)*beta(fulltemp+1:fulltemp+full);
                Aq = Azb + AQ(:,fulltemp+1:fulltemp+full)*beta(fulltemp+1:fulltemp+full);
            end
        else
            q = zb + Q(:,1:i)*beta;
            Aq = Azb + AQ(:,1:i)*beta;
        end
        Q(:,i+1) = q(:);
        AQ(:,i+1) = Aq(:);
    else
        % "strongly" truncated recursion for the update of the descent direction
        % (only the previous direction is considered)
        beta = -(Azb(:)'*w(:))/(w(:)'*w(:));
        Beta(i) = beta;
        q = zb + beta*q;
        Aq = Azb + beta*Aq;
    end
    Xiter(:,i) = real(x(:));
    RelErr(i) = norm(Xiter(:,i)-xex(:))/nx;
    RelRes(i) = norm(real(res(:)))/nb;
    % evaluating if the stopping criterion is satisfied
    if i>1
        if strcmp(stopC,'stab') && alStop == 0
            if (abs(RelRes(i-1)-RelRes(i))/RelRes(i-1))<tau
                stopCit = i;
                alStop = 1;
            end
        elseif strcmp(stopC,'discr') && alStop == 0
            if RelRes(i) < 1.01*sigma
                stopCit = i;
                alStop = 1;
            end
        end
    end
end
