40 lines
1.1 KiB
Plaintext
40 lines
1.1 KiB
Plaintext
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
|