comment
The CG method
endcomment

function cg (A,b,x0=0,f=4)
## Solves the linear system Ab=x, iterating with starting point x0.
## If x0 is missing the algorithm starts with 0.
## A must be positive definite.
	n=length(A);
	if argn==2; x=zeros(size(b)); else; x=x0; endif;
	h=b-A.x; r=h;
	alt=1e300;
	count=0;
	repeat
		fehler=sqrt(r'.r);
		if fehler~=0; return x; endif;
		if count>f*n;
			if fehler>alt; return x; endif;
		endif;
		count=count+1;
		alt=fehler;
		a=(r'.r)/(h'.A.h);
		x=x+a*h;
		rn=r-a*A.h;
		h=rn+(rn'.rn)/(r'.r)*h;
		r=rn;
	end;
endfunction

function cgn (A,b,x0=0)
## Solves the linear system Ab=x, iterating with starting point x0.
## If x0 is missing the algorithm starts with 0.
## It iterates only dim(A) times.
	n=length(A);
	if argn==2; x=zeros(size(b)); else; x=x0; endif;
	p=b-A.x; r=p;
	loop 1 to n;
		fehler=r'.r;
		if fehler~=0; return {x,#}; endif;
		a=(r'.r)/(p'.A.p);
		x=x+a*p;
		rn=r-a*A.p;
		p=rn+(rn'.rn)/(r'.r)*p;
		r=rn;
	end;
	return {x,n};
endfunction

function cginv(A)
	B=A; n=cols(A); Id=id(n); null=zeros(n,1);
	loop 1 to n;
		B[:,#]=cg(A,Id[:,#],null);
	end;
	return B;
endfunction