%Es soll das cg-Verfahren ausprobiert werden am Modellproblem Delta u=0
%auf dem Quadrat [0 1]x[0 1]. Die Randwerte sind Eins auf den Quadratseiten
%parallel zur x-Achse und Null auf den Quadratseiten parallel zur y-Achse.
%Das Residuum r wird mit r=residuum(u) und das Produkt  A*h mit
%z=Amalh(h) berechnet, zwei separat formulierten Prozuduren. 
%Wir verwenden die Formeln und Bezeichnungen aus dem Buch, S. 148.
%Die Lösung x hat formal die Form einer Matrix x(j,k), j,k=0,...,n.
%Nur die (n-1)^2 inneren Punke 0<j,k<n werden berechnet. Die Randpunkte x(j,k)
%mit j,k in {0,n} werden vorbesetzt und spielen die Rolle der rechten Seite.
%Da von Matrizen, interpretiert als Vektoren nicht ad hoc Skalarprodukte
%gebildet werden können, muß dafür auch ein separates Programm z=skalarpr(x,y)
%geschrieben werden. Aus matlab-technischen Gründen laufen die Indizes 
%von 1 bis n+1, und nicht von Null bis n.
n=33;%n=100 kommt bei norm(r,inf)<1e-10 mit 213 Schritten zum Ende.
x=zeros(n+1,n+1);r=x;h=x;hilfs=x;        %Startwert Nullvektor
x(1,:)=ones(size(x(1,:)));x(n+1,:)=ones(size(x(1,:))); %Randwerte
h=residuum(x);   r=h;   step=0;          %(5.53)
tic  %Stoppuhr beginnt.
y=norm(r,inf);Y=[y];
while y>=1e-10 & step <=(n-1)^2 
  step=step+1;           hilfs=Amalh(h);
  alpha=skalarpr(r,r)/skalarpr(h,hilfs); %(5.54)
  x=x+alpha*h;    r2hilfs=skalarpr(r,r); %(5.55)
  r=r-alpha*hilfs;                       %(5.56)
  beta=skalarpr(r,r)/r2hilfs;            %(5.57)
  h_alt=h;
  h=r+beta*h;                            %(5.58)
  h0=skalarpr(h_alt,Amalh(h));           %Test auf Orthogonalität
  y=norm(r,inf);Y=[Y,y];
end; %while
toc %Stoppuhr stoppt.
semilogy(Y); print -dps cgfehlera.ps %Y protokolliert Fehlerabfall
    %===========================Ergebnisse=====================================
    %n=33: elapsed_time = 6.613817 Sekunden in 68 Schritten.
    %Fehler des Residuums: 0.00000000008149 = 8.149e-11.
    %Das SOR-Verfahren braucht zum Vergleich bei einem optimalen Parameter
    %omega=1.82 auf demselben Rechner zur Erreichung derselben Genauigkeit
    %161 Schritte und 353.641601 Sekunden, also etwa 53-mal so lange.
    %auf Sun-ultra1: SOR: 161 Schritte bei 22.6717 Sekunden,
    %cg: 68 Schritte in 0.374341 Sekunden, also wieder ähnlicher Faktor 61.
    %n=10: elapsed_time=0.070138 Sekunden (Sun-ultra1, ohne Bildschirmanzeige), 
    %13 Iterationsschritte, Fehler des Residuums: nicht erkennbar.
