Overview
Comment:Commented out the cariabel iteration stuff
Downloads: Tarball | ZIP archive
Timelines: family | ancestors | descendants | both | origin/master | trunk
Files: files | file ages | folders
SHA3-256: f4415c4d81a68b883e92db36dce2cf26b53ee53d9594cf0adefe2e9a64cc6b9a
User & Date: gawthrop@users.sourceforge.net on 1998-08-15 09:30:05.000
Other Links: branch diff | manifest | tags
Context
1998-08-15
09:33:25
Deleted the commented out stuff check-in: 6fe802cf1f user: gawthrop@users.sourceforge.net tags: origin/master, trunk
09:30:05
Commented out the cariabel iteration stuff check-in: f4415c4d81 user: gawthrop@users.sourceforge.net tags: origin/master, trunk
08:26:30
This is variable interations version - now going on to fixed
iterations check-in: 1d3d0b408f user: gawthrop@users.sourceforge.net tags: origin/master, trunk
Changes
Unified Diff Ignore Whitespace Patch
Modified mttroot/mtt/bin/trans/p/mtt_sparse.p from [4258b05f2c] to [f61ecb354a]. 
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
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

PROCEDURE mtt_sparse(	 b   : glnarray;
			 n   : integer;
		     VAR x   : glnarray;
		     VAR rsq : real);
{*
###############################################################
## Version control history
###############################################################
## $Id$
## $Log$




## Revision 1.1  1998/08/13 14:58:35  peterg
## Initial revision
##
## Revision 1.1  1998/08/13 08:40:40  peterg
## Initial revision
##
###############################################################
*}

(* Programs using routine SPARSE must define the type
TYPE
   glnarray = ARRAY [1..n] OF real;
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.000001;

   
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 10*n 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;











>
>
>
>



















|
>





>

|
|
|




|

|






|

















>
|
>


|

>


>

|
|
>




>







|









|


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
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

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
##
## Revision 1.1  1998/08/13 08:40:40  peterg
## Initial revision
##
###############################################################
*}

(* Programs using routine SPARSE must define the type
TYPE
   glnarray = ARRAY [1..n] OF real;
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 ]