function [x,iter] ...
	= cg(A,x,b,itermax,epsi)
% Input: x-Startvektor
%        b-rechte Seite
%        A-Matrix
%        itermax-max. Iterationszahl
%        epsi-realtive Genauigkeit

% Output: x-Loesungsvektor
%         iter-benoetigte Anzahl Iterationen

% Berechne Startresiduum: 
s=[];
r=[];
w=[];
s=A*x;
r=s-b;
iter=1;
% Vorkonditionierung
w=r;
rw=w'*r;
rwv=rw;
beta=0;
abbr=epsi*epsi*rwv;
disp(['CG Iterationsschritte  <r,w> ']);
while ((abbr < rw) & (iter < itermax))
  iter = iter + 1;
  s    = w    + beta*s;
  w    = A * s;
  sas  = w'*s;
  alpha= -rw/sas;
  x    = x + alpha *s;
  r    = r + alpha *w;          
  w    = r;
  rw   = w'*r;
  beta = rw/rwv;
  rwv  = rw;      
  disp(['             ' num2str([iter,rw])]);
end;  

%disp(['CG Iterationsschritte: ' num2str(iter)]);
