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@MITMC )
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 herein 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^21)*(h2/(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))$
