function [x,X,xnorm,resnorm]=cgls(A,b,N)
%                                K,g,N
    x=zeros(length(A),1);
    r=b-A*x;
    d=A'*r;
    X=zeros(length(A),length(A))
    for k=1:N
         a=(norm(A'*r,2)/norm(A*d),2))^2;
         x2=x+a*d;
         r2=r-a*A*d;
         d2=A'*r2+((norm(A'*r2,2)/norm(A'*r,2))^2)*d;

         x=x2;
         r=r2;
         d=d2;
         X(:,k)=x;
    end

    xnorm=zeros(N,1);
    for j=1:N
        xnorm(j,:)=(norm(x(j,:),2));
    end
    resnorm=zeros(N,1);
    res=zeros(length(K),N);
    for j=1:N
        res(:,j)=g-K*x(:,j);
        resnorm(j,:)=(norm(res(j,:),2));
    end
end