archive/share/trash/kach.mac

   1 /*      Hacijan Linear Programming Algorithm            */
   2
   3 /*      Implementation by Jim Purtilo ( JIMP@MIT-MC )
   4                           Kent State University
   5                           November 1979
   6
   7         First added to share directory January 5, 1980
   8                                                         */
   9
  10
  11 /*      This version has more hooks for adding users definable
  12         criterion for stopping and restarting in mid alg;
  13         trying to solve the difficulty in loosing positive
  14         definitness in time.                            */
  15
  16 /*                                                      */
  17
  18 /*      Primitives defined are:
  19
  20         HACH( a , b , n , m , l )
  21
  22                 where a is the matrix and b is the
  23                 column vector such that a.x<b is the
  24                 problem to be solved; n is number variables;
  25                 l is as described in the theory, roughly
  26                 a length, and also what, for instance, might
  27                 be output from GETL.  This is the main alg-
  28                 orithm, and return either a vector x such
  29                 that the above inequality is satisfied, or,
  30                 hopefully, returns why not. m is the number
  31                 of inequalities.
  32
  33         CHECK( a , b , x , m )
  34
  35                 is a primitive which CHECKs the inequality
  36                 a.x<b, returning either 0 ( when no rows
  37                 in a.x violate the equality ) or j ( where
  38                 j is in fact a row which does violate the
  39                 inequality ).  This is the new check, which
  40                 will always check first that the last row iterated
  41                 on is satisfied in the new inequality, before
  42                 giving the rest of the rows a chance to be
  43                 iterated on.  New strategies for row choice
  44                 will be implemented as the theory developes
  45                 at Kent CSR Group.   m is the number rows
  46                 ( inequalities ).
  47
  48
  49
  50         GETL( a , b , n , m )
  51
  52                 returns a value l ( for length ) as required
  53                 by and defined in the theory. a and b are as
  54                 above, also n is number variables, and m is
  55                 the number inequalities in a.
  56
  57 Several utilities are here-in defined:
  58
  59                 FX              FQ              FH
  60                 ITERATE         Z               LOG2
  61                 ACCELERATE      CRITERION
  62
  63                                                         */
  64 /*      Note that the desired precision in bigfloat desired
  65         for the calculations of the elipsoids should be set
  66         prior to any calculations by fpprec:<num digits>        */
  67 /*      Note that there is a global variable show%all[false]
  68         which when true will print the vector x at each
  69         iteration.                      */
  70
  71 declare(z,special)$
  72
  73 define_variable(show%all,false,boolean)$
  74
  75 define_variable(last%row,0,any)$
  76
  77 hach(a,b,n,m,l):=
  78         block([x,i,k,h,q,last%row],
  79
  80                 last%row:0,
  81                 h:bfloat(2^(2*l)*ident(n)),
  82                 x:genmatrix(z,n,1),
  83
  84                 /*  Note this is a "do forever...", not
  85                 automatically checking the counter k against
  86                 in the stopping criterion of the theory.
  87                 The reason for this is that oft' times we
  88                 like to play with whacky L's, and don't want
  89                 to mess with the resulting random stoppage.
  90                 Next update will have a settable switch...
  91                 and possibly also an improved polynomial
  92                 stopping time...                */
  93
  94                 for k:0 next k+1
  95                 do  if (i:check(a,b,x,m))=0
  96                     then return(["win",k,"iterations",x])
  97                     else (
  98                         if show%all then ( print(k,x),print(" ")),
  99                         if not criterion(a,b,x,n,m,l,k)
 100                                  then  iterate(x,a,h,q,i,n)
 101                                 else
 102
 103                 /*      In here would be code for
 104                         resetting, say, h, x , lastrow, etc.
 105                         Then we start the iteration process
 106                 again, hopefully accellerated.
 107                 Note that for debugging purposes, code here
 108                 could be a break...
 109                 */
 110                         (       accelerate(a,b,x,n,m,l,k),
 111                                 iterate(x,a,h,q,i,n) ) )
 112
 113
 114         )$
 115 accelerate(a,b,x,n,m,l,k):=false$
 116
 117 criterion(a,b,x,n,m,l,k):=true$
 118 /* for now the above will keep things runningtill the users defined
 119         what he wants... */
 120
 121 check(a,b,x,m):=
 122         block([y,j,j1],
 123                 j1:0,
 124                 y:a.x,
 125                 if numberp( last%row ) then
 126                         if last%row#0 then
 127                                 if ">="(y[last%row,1],b[last%row,1])
 128                                         then j1:last%row,
 129                 if j1#0 then return( j1 ),
 130                 for j:1 thru m
 131                         do if ">="(y[j,1],b[j,1]) then j1:j,
 132                 j1
 133 /*      Note that this strategy is that if last%row is ok, or
 134         is undefined due to a user CHECKing a system, the
 135         the first row violating a.x<b is used           */
 136         )$
 137
 138 getl(a,b,n,m):=
 139         block([i,j,l],
 140                 l:product(
 141                         (abs(b[i,1])+1)*product((abs(a[j,i])+1),j,1,m),
 142                           i,1,n ),
 143                 l:1+log2(n*m*l)
 144         )$
 145
 146 log2(x):=bfloat(log(x)/log(2))$
 147
 148 z[i,j]:=0$
 149
 150 iterate(x,a,h,q,i,n):=( q:fq(a,h,i,n),
 151                 x:fx(x,a,h,q,i,n),
 152                 h:fh(a,h,q,i,n)
 153
 154            )$
 155
 156 fx(x,a,h,q,i,n):=x-(n+1)^(-1)*q^(-1/2)*h.transpose(row(a,i))$
 157
 158 fh(a,h,q,i,n):=n^2/(n^2-1)*(h-2/(n+1)*q^-1*
 159                 h.transpose(row(a,i)).row(a,i).transpose(h))$
 160
 161 fq(a,h,i,n):=row(a,i).h.transpose(row(a,i))$
 162