%% CGLS/CGQR

% [A,n,nc,entries,rep,field,symm] = mmread('well1033.mtx');
% [b1] = mmread('well1033_rhs1.mtx');
[A,n,nc,entries,rep,field,symm] = mmread('illc1850.mtx');
[b] = mmread('illc1850_rhs1.mtx');

[m,n]  = size(A);
true_sol=0;  %  
if true_sol==1,
   xtrue  = (n : -1: 1)';
   b      = A*xtrue;
end

atol   = 1.0e-6;
btol   = 1.0e-6;
conlim = 1.0e+10;
itnlim = 10*n;
show   = 1;
damp   = 0;
[ x, istop, itn, r1norm, r2norm, Anorm, Acond, Arnorm, xnorm, var ] ...
      = lsqrSOL( m, n, A, b, damp, atol, btol, conlim, itnlim, show );
disp(['lsqr iterations ' int2str(itn)])

  wait
  
%% CGLS

[x1,flag,resNE,iter] = cgls(A,b,damp,atol,n,1);%,x0)                  
disp(['cgls iterations ' int2str(iter)])

figure(1),clf
plot(x),hold
plot(x1,'r')
if true_sol==1,plot(xtrue,'g'),end