1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:

function [ x, residuals ] = cgnr( B, c, maxres )
%CGNR Summary of this function goes here
%   Detailed explanation goes here
x = 0;
r = -c;
residuals = [norm(r)];
k = 0;
while norm(r) > maxres
    if k == 0
        d = -r;      
    else
        be = norm(r)^2\norm(z)^2;
        d = -r + be*d;
    end
    residuals = [residuals norm(r)];
    z = r;
    w = B*d;
    a = norm(r)^2\d'*B'*w;  
    x = x + a*d;
    x
    r = r + a*B'*w;
    k = k+1;
end
%disp(length(residuals));
end

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:
77:
78:
79:
80:
81:
82:
83:
84:
85:
86:
87:
88:
89:
90:
91:
92:
93:
94:
95:
96:
97:
98:
99:
100:
101:
102:
103:
104:
105:
106:
107:
108:
109:
110:
111:
112:
113:
114:
115:
116:
117:
118:

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');