
repeat = 1;
ns = unique(floor(exp(log(1):0.2:log(100000))));

max_time = 0.3;

maxrelres = 1e-5;

time_abgabe = zeros(1, length(ns));
time_lapack = zeros(1, length(ns));
res         = zeros(1, length(ns));
max_res      = zeros(1, length(ns));
iterations   = zeros(1, length(ns));

act_pos = 0;

for i=1:length(ns)
   act_pos = act_pos+1;

   n = ns(i);

   sumtime_abgabe = 0;
   sumtime_lapack = 0;
   sumres         = 0;
   summax_res      = 0;
   sumiter = 0;

   for j=1:repeat
      a = randn(n,1)+5;
      b = randn(n,1);
      b(end) = inf;
      c = randn(n,1);
      c(1) = inf;
      y = randn(n,1);

      B = spdiags([[c(2:end); inf] a [inf; b(1:end-1)]], -1:1, n, n);

      c = randn(n,1);
      nBc = norm(B'*c);
      funfact = (1+randn(1)*0.5);
      if( funfact < 0.5 )
         funfact = 0.5;
      end
      maxres = maxrelres*nBc*funfact;

      % run Abgabe
      tic;
      [x, residuals] = cgnr(B, c, maxres);
      sumtime_abgabe = sumtime_abgabe + toc;

      sumiter = sumiter + length(residuals);

      % run LAPACK
      tic;
      x_lapack = B\c;
      sumtime_lapack = sumtime_lapack + toc;

      realres = norm(B'*(B*x - c));
      if( realres > maxres )
         realres
         maxres
         save('fehler_beispiel.mat', 'B', 'c', 'maxres');
         disp('Die von "cgnr" ausgerechnete Approximation hat zu grosses Residuum !');
         disp('Es ist realres > maxres!');
         disp(' Geben sie "load(''fehler_beispiel.mat'')" ein um B, c, maxres zu laden, welche auf das Problem gefuehrt haben.');
         error('Fehler!');
      end

      sumres = sumres + realres;
      summax_res = summax_res + maxres;

   end

   time_abgabe(i) = sumtime_abgabe/repeat;
   time_lapack(i) = sumtime_lapack/repeat;
   res(i)         = sumres/repeat;
   max_res(i)     = summax_res/repeat;

   iterations(i) = sumiter/repeat;

   if( (sumtime_abgabe+sumtime_lapack)>max_time )
      break;
   end
end

time_abgabe = time_abgabe(1:act_pos);
time_lapack = time_lapack(1:act_pos);
res    = res(1:act_pos);
max_res     = max_res(1:act_pos);
iterations = iterations(1:act_pos);
ns = ns(1:act_pos);

subplot(2,2,1);
loglog(ns, time_abgabe, 'b-', ns, time_lapack, 'r--');
title(['Laufzeiten']);
legend('Abgabe', 'Direkter Loeser (LAPACK)', ...
         'Location', 'NorthWest');

subplot(2,2,2);
big_ns = (ns>20);
loglog(ns(big_ns), res(big_ns), 'b-', 'LineWidth', 2);
hold on;
loglog( ns(big_ns), max_res(big_ns), 'r-' )
legend('Wirkliches Residuum der Abgabe', 'maxres Parameter', 'Location', 'NorthWest');

subplot(2,2,3);
plot( ns, iterations, 'b-' );
title('Mittlere Anzahl der Iterationen von CGNR');


subplot(2,2,4);
iters = length(residuals);
lastn = ns(end);
loglog(1:iters, residuals, 'b-', [1, iters], [maxres, maxres], 'k--');
title(['Konvergenzhistorie des letzten Beispiels; n = ', int2str(lastn)]);
legend('residuals', 'maxres', 'Location', 'NorthEast');

