clear all;
close all;
% Grober Test fuer konjugierte Gradienten
n=1000;  % es gehen ziemlich grosse n 
a=speye(n,n); % wir bauen eine sparse nxn Matrix
for i=1:2*n   % und geben ihr 2n zufaellige Elemente und die Diagonale
    in=1+floor(n*rand(1,1));  % das ist zwischen 1 und n 
    jn=1+floor(n*rand(1,1));  % ebenso
    a(in,jn)=2*rand(1,1)-1;   % und das ist in [-1,1]
end
a=a*a';  % wir brauchen a symmetrisch und positiv definit
spy(a)   % und wir sehen uns die sparsity an.
Dimension=n
VonNullVerschiedeneElemente=nnz(a)
b=rand(n,1); % Eine zufaellige rechte Seite
disp('MATLAB rechnet die Loesung aus')
xdir=a\b;    % und die exakte Loesung, zum Vergleich
disp('Wir rechnen die Loesung aus, vorher ENTER druecken')
pause
% jetzt fangen wir mit CG an, Bezeichnungen wie im Buch 
x=b;
r=b-a*x;
g=r;
iter=1; % wir zaehlen die Iterationen
nr(iter)=norm(r); % und messen die Konvergenz dreifach 
kryval=(xdir-x)'*(a*(xdir-x));
kry(iter)=sqrt(kryval);
infnorm(iter)=norm(x-xdir,Inf);
% wer will, kann die Orthogonalitaet testen
% rr=zeros(n,n);
% gg=zeros(n,n);
while nr(iter)>1.0e-5  % Rechenschlaife
%     rr(:,iter)=r;
    p=a*r; % ab hier wie im Buch, Seite 297
    rho=p'*p;
    sigma=r'*p;
    tau=g'*r;
    t=tau/sigma;
    x=x+t*r;
    g=g-t*p;
%     gg(:,iter)=g;
    gamma=(t*t*rho-tau)/tau;
    r=g+gamma*r;
    iter=iter+1;  % Zaehlen und messen
    kryval=(xdir-x)'*(a*(xdir-x));
    kry(iter)=sqrt(kryval);
    nr(iter)=norm(r);
    infnorm(iter)=norm(x-xdir,Inf);
end
figure
semilogy(1:iter,nr,1:iter,kry,1:iter,infnorm)
legend('Residualnorm','Krylov-Fehlernorm','inf-norm des Fehlers')
% rr=rr(1:n,1:n);
% gg=gg(1:n,1:n);
% nrar=norm(rr'*a*rr-diag(diag(rr'*a*rr)))
% ngg=norm(gg'*gg-diag(diag(gg'*gg)))
letzteResiduumnorm=nr(iter)
Iterationen=iter
