1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
|
PROCEDURE mtt_sparse( b : glnarray;
n : integer;
VAR x : glnarray;
VAR rsq : real);
{*
###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.2 1998/08/15 08:26:30 peterg
## This is variable interations version - now going on to fixed
## iterations
##
## Revision 1.1 1998/08/13 14:58:35 peterg
## Initial revision
##
|
>
>
>
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
|
PROCEDURE mtt_sparse( b : glnarray;
n : integer;
VAR x : glnarray;
VAR rsq : real);
{*
###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.3 1998/08/15 09:30:05 peterg
## Commented out the cariabel iteration stuff
##
## Revision 1.2 1998/08/15 08:26:30 peterg
## This is variable interations version - now going on to fixed
## iterations
##
## Revision 1.1 1998/08/13 14:58:35 peterg
## Initial revision
##
|
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
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
|
in the main routine. They must also provide two routines,
PROCEDURE asub(x: glnarray; VAR y: glnarray; n: integer);
and
PROCEDURE atsub(x: glnarray; VAR z: glnarray; n: integer);
which calculate A*x and (A transpose)*x *)
LABEL 1,99;
CONST
{** eps = 0.00000;**}
{** eps2 := n*sqr(eps);
irst := 0;
1: irst := irst+1; **}
{* bsq := bsq+sqr(b[j]);*}
{* rp := rp+sqr(xi[j]) *}
IF (aden = 0.0) THEN BEGIN
writeln('pause in routine SPARSE');
writeln('very singular matrix'); {***readln ***}
END;
{***
rsq := 0.0;
**}
{***rsq := rsq+sqr(xj[j])**}
{***
IF ((rsq = rp) OR (rsq <= bsq*eps2)) THEN GOTO 99;
IF (rsq > rp) THEN BEGIN
***}{**
FOR j := 1 TO n DO BEGIN
x[j] := xi[j]
END;**}
{***
IF (irst >= 3) THEN GOTO 99;
GOTO 1
END;
rp := rsq;
***}
{**IF (gg = 0.0) THEN GOTO 99;**}
{***writeln('pause in routine SPARSE');
writeln('too many iterations'); readln; **}
{**99:writeln("---",iter,rsq);**}
|
<
<
<
<
<
<
<
<
<
<
<
<
<
<
<
<
<
<
<
<
<
<
<
<
<
<
<
<
<
<
<
<
|
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
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
|
in the main routine. They must also provide two routines,
PROCEDURE asub(x: glnarray; VAR y: glnarray; n: integer);
and
PROCEDURE atsub(x: glnarray; VAR z: glnarray; n: integer);
which calculate A*x and (A transpose)*x *)
LABEL 1,99;
CONST
iters = 500;
VAR
j,iter,irst: integer;
rp,gg,gam,eps2,dgg,bsq,anum,aden: real;
g,h,xi,xj: glnarray;
BEGIN
mtt_asub(x,xi,n);
rp := 0.0;
bsq := 0.0;
FOR j := 1 TO n DO BEGIN
xi[j] := xi[j]-b[j];
END;
mtt_atsub(xi,g,n);
FOR j := 1 TO n DO BEGIN
g[j] := -g[j];
h[j] := g[j]
END;
FOR iter := 1 TO iters DO BEGIN
mtt_asub(h,xi,n);
anum := 0.0;
aden := 0.0;
FOR j := 1 TO n DO BEGIN
anum := anum+g[j]*h[j];
aden := aden+sqr(xi[j])
END;
anum := anum/aden;
FOR j := 1 TO n DO BEGIN
xi[j] := x[j];
x[j] := x[j]+anum*h[j]
END;
mtt_asub(x,xj,n);
FOR j := 1 TO n DO BEGIN
xj[j] := xj[j]-b[j];
END;
mtt_atsub(xj,xi,n);
gg := 0.0;
dgg := 0.0;
FOR j := 1 TO n DO BEGIN
gg := gg+sqr(g[j]);
dgg := dgg+(xi[j]+g[j])*xi[j]
END;
gam := dgg/gg;
FOR j := 1 TO n DO BEGIN
g[j] := -xi[j];
h[j] := g[j]+gam*h[j]
END
END;
END;
|