function [Xiter, RelErr, RelRes, Alpha, istop, Beta, stopCit] = fcglsPoiss_rp_nn(A, b, x0, m, sigmaSq, beta, xex, full, c, stopC)
%
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);
Alpha = zeros(m,1);
nx = norm(xex(:));
nb = norm(b(:));
istop = m;
%
x = x0;
x(x<0)=0;
if nargin == 7
    full = 0;
    c = b + sigmaSq;
elseif nargin == 8
    c = b+sigmaSq;
end
if full
    Beta = zeros(m);
    Q = zeros(N,m+1);
    AQ = zeros(N,m+1);
else
    Beta = zeros(m,1);
end
wt = sqrt(c);
b = b - beta;
res = b - A*x; res = res(:); resp = res./wt;
z = A'*(resp./wt);
if max(x)~=0
    zb = x(:).*z;
else
    zb = z;
end
q = zb;
Azb = (A*zb)./wt;
Aq = Azb;
if full
    Q(:,1) = q(:);
    AQ(:,1) = Aq(:);
end
for i = 1:m
    w = Aq;
    theta = (resp(:)'*w(:))/((w(:))'*w(:));
    neg_ind = q < 0;
    alpha = min( theta, min( -x(neg_ind)./q(neg_ind) ) );
    if isempty(alpha)
        alpha = theta;
    end
    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;
        return
    end 
    Alpha(i) = alpha;
    x = x + alpha*q;
    res = res - alpha*(w.*wt);
    resp = res./wt;
    z = A'*(resp./wt);
    zb = x(:).*z(:); 
    Azb = (A*zb)./wt;
    if full
        beta = zeros(i,1);
        for j=1:i
            beta(j) = -(Azb(:)'*AQ(:,j))/(AQ(:,j)'*AQ(:,j));
        end
        Beta(1:i,i) = beta;
        if full ~=1
            q = zb + Q(:,1:full)*beta(1:full);
            Aq = Azb + AQ(:,1:full)*beta(1:full);
        else
            q = zb + Q(:,1:i)*beta;
            Aq = Azb + AQ(:,1:i)*beta;
        end
        Q(:,i+1) = q(:);
        AQ(:,i+1) = Aq(:);
    else
        beta = -(Azb(:)'*w(:))/(w(:)'*w(:));
        Beta(i) = beta;
        q = zb + beta*q;
        Aq = Azb + beta*Aq;
    end
    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
    Xiter(:,i) = real(x(:));
    RelErr(i) = norm(Xiter(:,i)-xex(:))/nx;
    RelRes(i) = norm(real(res(:)))/nb;
end
