function cgX (H:cpx, b:real column, x0:real column=none, f:index=10)
## Conjugate gradient method to solve Hx=b for compressed H.
##
## This is the method of choice for large, sparse matrices. In most
## cases, it will work well, fast, and accurate.
##
## H must be positive definite. Use cpxfit, if it is not.
##
## The accuarcy can be controlled with an additional parameter
## eps. The algorithm stops, when the error gets smaller then eps, or
## after f*n iterations, if the error gets larger. x0 is an optional
## start vector.
##
## H : compressed matrix (nxm)
## b : column vector (mx1)
## x0 : optional start point (mx1)
## f : number of steps, when the method should be restarted
##
## See: cpxfit, cg, cgXnormal
	if isvar("eps") then localepsilon(eps); endif;
	n=cols(H);
	if x0==none then x=zeros(size(b));
	else; x=x0;
	endif;
	loop 1 to 10
		r=b-cpxmult(H,x); p=r; fehler=r'.r;
		loop 1 to f*n
			if sqrt(fehler)~=0 then return x; endif;
			Hp=cpxmult(H,p);
			a=fehler/(p'.Hp);
			x=x+a*p;
			rn=r-a*Hp;
			fehlerneu=rn'.rn;
			p=rn+fehlerneu/fehler*p;
			r=rn; fehler=fehlerneu;
		end;
	end;
	return x;
endfunction