clear;
tic;
n=15;N=n^2;
A=gallery('poisson',n);
%e=ones(N,1);
%A=spdiags([0.16*e e -1.16*e],-1:1,N,N);
b=zeros(N,1);
rtol=1e-6;
maxit=250000;
xex = zeros(N,1);
x0=rand(N,1);
disp(sprintf('setup Ax=b:    %f',toc));

%runlist={'MR'};
runlist={'SD','MR','RNSD'};

% prediction
lambdamax=eigs(A,1,'LM');
lambdamin=eigs(A,1,'SM');
efactor = (lambdamax-lambdamin)/(lambdamax+lambdamin);
mu = eigs((A+A')/2,1,'SM');
sig = normest(A);
rfactor=sqrt(1-((mu^2)/(sig^2)));


for ri=1:length(runlist)
  if strcmp(runlist{ri},'SD')
    tic;
    [x{ri},rhist{ri},ehist{ri}]=algSD(A,b,rtol,maxit,x0,xex);
    iter{ri}=length(rhist{ri});
    disp(sprintf('ran %d iterations of algSD:    %f',iter{ri},toc));
    epredict{ri}=efactor.^(1:iter{ri});
    rpredict{ri}=ones(1,iter{ri});
  end

  if strcmp(runlist{ri},'MR')
    tic;
    [x{ri},rhist{ri},ehist{ri}]=algMR(A,b,rtol,maxit,x0,xex);
    iter{ri}=length(rhist{ri});
    disp(sprintf('ran %d iterations of algMR:    %f',iter{ri},toc));
    epredict{ri}=ones(1,iter{ri});
    rpredict{ri}=rfactor.^(1:iter{ri});
  end

  if strcmp(runlist{ri},'RNSD')
    tic;
    [x{ri},rhist{ri},ehist{ri}]=algRNSD(A,b,rtol,maxit,x0,xex);
    iter{ri}=length(rhist{ri});
    disp(sprintf('ran %d iterations of algRNSD:    %f',iter{ri},toc));
    epredict{ri}=efactor.^(1:iter{ri});
    rpredict{ri}=ones(1,iter{ri});
  end
end


%%%%% plot
figure(1);clf;title('residuals');xlabel('iteration');ylabel('||r||');hold on;
figure(2);clf;title('errors');xlabel('iteration');ylabel('||e||');hold on;
colors={'b','g','r','k','c','m'};
for ri=1:length(runlist)
  figure(1);
  loglog(1:iter{ri},rhist{ri}/rhist{ri}(1),colors{ri});
  figure(2);
  loglog(1:iter{ri},ehist{ri}/ehist{ri}(1),colors{ri});
end
figure(1);
legend(runlist);
figure(2);
legend(runlist);
if(0)
for ri=1:length(runlist)
  figure(1);
  loglog(1:iter{ri},rpredict{ri},strcat(colors{ri},'--'));
  figure(2);
  loglog(1:iter{ri},epredict{ri},strcat(colors{ri},'--'));
end
end
figure(1);
set(gca,'XScale','log','YScale','log');
figure(2);
set(gca,'XScale','log','YScale','log');