1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
|
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
|
