Differences From Artifact [f61ecb354a]:

To Artifact [4534575ce1]:


1
2
3
4
5
6
7
8
9
10



11
12
13
14
15
16
17
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
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
  {** eps	 = 0.00000;**}
   iters = 500;
   
VAR
   j,iter,irst: integer;
   rp,gg,gam,eps2,dgg,bsq,anum,aden: real;
   g,h,xi,xj: glnarray;
   
BEGIN
{**   eps2 := n*sqr(eps);
    irst := 0;
1:   irst := irst+1; **}
   mtt_asub(x,xi,n);
   rp := 0.0;
   bsq := 0.0;
   FOR j := 1 TO n DO BEGIN
     {* bsq := bsq+sqr(b[j]);*}
      xi[j] := xi[j]-b[j];
     {* rp := rp+sqr(xi[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;
      IF (aden = 0.0) THEN BEGIN
         writeln('pause in routine SPARSE');
         writeln('very singular matrix'); {***readln ***}
      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);
      {***
       rsq := 0.0;
       **}
      FOR j := 1 TO n DO BEGIN
         xj[j] := xj[j]-b[j];
         {***rsq := rsq+sqr(xj[j])**}
      END;
      {***
      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;
       ***}
      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;
      {**IF (gg = 0.0) THEN GOTO 99;**}
      gam := dgg/gg;
      FOR j := 1 TO n DO BEGIN
         g[j] := -xi[j];
         h[j] := g[j]+gam*h[j]
      END
   END;
   {***writeln('pause in routine SPARSE');
   writeln('too many iterations'); readln; **}

{**99:writeln("---",iter,rsq);**}
END;

MTT: Model Transformation Tools
GitHub | SourceHut | Sourceforge | Fossil RSS ]