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
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
|
/* Hacijan Linear Programming Algorithm */
/* Implementation by Jim Purtilo ( JIMP@MIT-MC )
Kent State University
November 1979
First added to share directory January 5, 1980
*/
/* This version has more hooks for adding users definable
criterion for stopping and restarting in mid alg;
trying to solve the difficulty in loosing positive
definitness in time. */
/* */
/* Primitives defined are:
HACH( a , b , n , m , l )
where a is the matrix and b is the
column vector such that a.x<b is the
problem to be solved; n is number variables;
l is as described in the theory, roughly
a length, and also what, for instance, might
be output from GETL. This is the main alg-
orithm, and return either a vector x such
that the above inequality is satisfied, or,
hopefully, returns why not. m is the number
of inequalities.
CHECK( a , b , x , m )
is a primitive which CHECKs the inequality
a.x<b, returning either 0 ( when no rows
in a.x violate the equality ) or j ( where
j is in fact a row which does violate the
inequality ). This is the new check, which
will always check first that the last row iterated
on is satisfied in the new inequality, before
giving the rest of the rows a chance to be
iterated on. New strategies for row choice
will be implemented as the theory developes
at Kent CSR Group. m is the number rows
( inequalities ).
GETL( a , b , n , m )
returns a value l ( for length ) as required
by and defined in the theory. a and b are as
above, also n is number variables, and m is
the number inequalities in a.
Several utilities are here-in defined:
FX FQ FH
ITERATE Z LOG2
ACCELERATE CRITERION
*/
/* Note that the desired precision in bigfloat desired
for the calculations of the elipsoids should be set
prior to any calculations by fpprec:<num digits> */
/* Note that there is a global variable show%all[false]
which when true will print the vector x at each
iteration. */
declare(z,special)$
define_variable(show%all,false,boolean)$
define_variable(last%row,0,any)$
hach(a,b,n,m,l):=
block([x,i,k,h,q,last%row],
last%row:0,
h:bfloat(2^(2*l)*ident(n)),
x:genmatrix(z,n,1),
/* Note this is a "do forever...", not
automatically checking the counter k against
in the stopping criterion of the theory.
The reason for this is that oft' times we
like to play with whacky L's, and don't want
to mess with the resulting random stoppage.
Next update will have a settable switch...
and possibly also an improved polynomial
stopping time... */
for k:0 next k+1
do if (i:check(a,b,x,m))=0
then return(["win",k,"iterations",x])
else (
if show%all then ( print(k,x),print(" ")),
if not criterion(a,b,x,n,m,l,k)
then iterate(x,a,h,q,i,n)
else
/* In here would be code for
resetting, say, h, x , lastrow, etc.
Then we start the iteration process
again, hopefully accellerated.
Note that for debugging purposes, code here
could be a break...
*/
( accelerate(a,b,x,n,m,l,k),
iterate(x,a,h,q,i,n) ) )
)$
accelerate(a,b,x,n,m,l,k):=false$
criterion(a,b,x,n,m,l,k):=true$
/* for now the above will keep things runningtill the users defined
what he wants... */
check(a,b,x,m):=
block([y,j,j1],
j1:0,
y:a.x,
if numberp( last%row ) then
if last%row#0 then
if ">="(y[last%row,1],b[last%row,1])
then j1:last%row,
if j1#0 then return( j1 ),
for j:1 thru m
do if ">="(y[j,1],b[j,1]) then j1:j,
j1
/* Note that this strategy is that if last%row is ok, or
is undefined due to a user CHECKing a system, the
the first row violating a.x<b is used */
)$
getl(a,b,n,m):=
block([i,j,l],
l:product(
(abs(b[i,1])+1)*product((abs(a[j,i])+1),j,1,m),
i,1,n ),
l:1+log2(n*m*l)
)$
log2(x):=bfloat(log(x)/log(2))$
z[i,j]:=0$
iterate(x,a,h,q,i,n):=( q:fq(a,h,i,n),
x:fx(x,a,h,q,i,n),
h:fh(a,h,q,i,n)
)$
fx(x,a,h,q,i,n):=x-(n+1)^(-1)*q^(-1/2)*h.transpose(row(a,i))$
fh(a,h,q,i,n):=n^2/(n^2-1)*(h-2/(n+1)*q^-1*
h.transpose(row(a,i)).row(a,i).transpose(h))$
fq(a,h,i,n):=row(a,i).h.transpose(row(a,i))$
|
