1 /* lrslib.h (vertex enumeration using lexicographic reverse search) */ 2 3 #define TITLE "lrslib " 4 #define VERSION "v.4.2c, 2010.4.26" 5 #define AUTHOR "\n*Copyright (C) 1995,2010, David Avis avis@cs.mcgill.ca " 6 7 /* This program is free software; you can redistribute it and/or modify 8 it under the terms of the GNU General Public License as published by 9 the Free Software Foundation; either version 2 of the License, or 10 (at your option) any later version. 11 12 This program is distributed in the hope that it will be useful, 13 but WITHOUT ANY WARRANTY; without even the implied warranty of 14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15 GNU General Public License for more details. 16 17 You should have received a copy of the GNU General Public License 18 along with this program; if not, write to the Free Software 19 Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 20 */ 21 /*Ver 4.2* library version */ 22 /******************************************************************************/ 23 /* See http://cgm.cs.mcgill.ca/~avis/C/lrs.html for usage instructions */ 24 /******************************************************************************/ 25 /* Selection of arithmetic package */ 26 /*************************************/ 27 #ifdef LONG 28 #define ARITH "lrslong.h" /* lrs long integer arithmetic package */ 29 #else 30 #ifdef GMP 31 #define ARITH "lrsgmp.h" /* lrs wrapper for gmp multiple precsion arithmetic */ 32 #else 33 #define ARITH "lrsmp.h" /* lrs multiple precsion arithmetic */ 34 #define MP 35 #endif 36 #endif 37 38 #include ARITH 39 40 #ifdef SIGNALS 41 #include <signal.h> 42 #include <unistd.h> 43 #define errcheck(s,e) if ((long)(e)==-1L){ perror(s);exit(1);} 44 #endif 45 46 #ifdef TIMES 47 void ptimes (); 48 double get_time(); 49 #endif 50 51 52 #define CALLOC(n,s) xcalloc(n,s,__LINE__,__FILE__) 53 54 /* make this include file includable in a C++ file 55 Note that lrsgmp.h should not be inside an extern "C" {} block 56 */ 57 #ifdef __cplusplus 58 extern "C" { 59 #endif 60 61 62 /*************/ 63 /* typedefs */ 64 /*************/ 65 66 67 /******************************************************************************/ 68 /* Indexing after initialization */ 69 /* Basis Cobasis */ 70 /* --------------------------------------- ----------------------------- */ 71 /* | i |0|1| .... |lastdv|lastdv+1|...|m| | j | 0 | 1 | ... |d-1| d | */ 72 /* |-----|+|+|++++++|++++++|--------|---|-| |----|---|---|-----|---|+++++| */ 73 /* |B[i] |0|1| .... |lastdv|lastdv+1|...|m| |C[j]|m+1|m+2| ... |m+d|m+d+1| */ 74 /* -----|+|+|++++++|++++++|????????|???|?| ----|???|???|-----|???|+++++| */ 75 /* */ 76 /* Row[i] is row location for B[i] Col[j] is column location for C[j] */ 77 /* ----------------------------- ----------------------------- */ 78 /* | i |0|1| ..........|m-1|m| | j | 0 | 1 | ... |d-1| d | */ 79 /* |-------|+|-|-----------|---|-| |------|---|---|--- |---|++++| */ 80 /* |Row[i] |0|1|...........|m-1|m| |Col[j]| 1 | 2 | ... | d | 0 | */ 81 /* --------|+|*|***********|***|*| ------|***|***|*****|***|++++| */ 82 /* */ 83 /* + = remains invariant * = indices may be permuted ? = swapped by pivot */ 84 /* */ 85 /* m = number of input rows n= number of input columns */ 86 /* input dimension inputd = n-1 (H-rep) or n (V-rep) */ 87 /* lastdv = inputd-nredundcol (each redundant column removes a dec. var) */ 88 /* working dimension d=lastdv-nlinearity (an input linearity removes a slack) */ 89 /* obj function in row 0, index 0=B[0] col 0 has index m+d+1=C[d] */ 90 /* H-rep: b-vector in col 0, A matrix in columns 1..n-1 */ 91 /* V-rep: col 0 all zero, b-vector in col 1, A matrix in columns 1..n */ 92 /******************************************************************************/ 93 94 typedef struct lrs_dic_struct /* dynamic dictionary data */ 95 { 96 lrs_mp_matrix A; 97 long m; /* A has m+1 rows, row 0 is cost row */ 98 long m_A; /* =m or m-d if nonnegative flag set */ 99 long d; /* A has d+1 columns, col 0 is b-vector */ 100 long d_orig; /* value of d as A was allocated (E.G.) */ 101 long lexflag; /* true if lexmin basis for this vertex */ 102 long depth; /* depth of basis/vertex in reverse search tree */ 103 long i, j; /* last pivot row and column pivot indices */ 104 lrs_mp det; /* current determinant of basis */ 105 lrs_mp objnum; /* objective numerator value */ 106 lrs_mp objden; /* objective denominator value */ 107 long *B, *Row; /* basis, row location indices */ 108 long *C, *Col; /* cobasis, column location indices */ 109 struct lrs_dic_struct *prev, *next; 110 } 111 lrs_dic; 112 113 typedef struct lrs_dat /* global problem data */ 114 { 115 lrs_mp_vector Gcd; /* Gcd of each row of numerators */ 116 lrs_mp_vector Lcm; /* Lcm for each row of input denominators */ 117 118 lrs_mp sumdet; /* sum of determinants */ 119 lrs_mp Nvolume; /* volume numerator */ 120 lrs_mp Dvolume; /* volume denominator */ 121 lrs_mp boundn; /* objective bound numerator */ 122 lrs_mp boundd; /* objective bound denominator */ 123 long unbounded; /* lp unbounded */ 124 char fname[100]; /* input file name from line 1 of input */ 125 126 long *inequality; /* indices of inequalities corr. to cobasic ind */ 127 /* initially holds order used to find starting */ 128 /* basis, default: m,m-1,...,2,1 */ 129 long *facet; /* cobasic indices for restart in needed */ 130 long *redundcol; /* holds columns which are redundant */ 131 long *linearity; /* holds cobasic indices of input linearities */ 132 long *minratio; /* used for lexicographic ratio test */ 133 long *temparray; /* for sorting indices, dimensioned to d */ 134 long *isave, *jsave; /* arrays for estimator, malloc'ed at start */ 135 long inputd; /* input dimension: n-1 for H-rep, n for V-rep */ 136 137 long m; /* number of rows in input file */ 138 long n; /* number of columns in input file */ 139 long lastdv; /* index of last dec. variable after preproc */ 140 /* given by inputd-nredundcol */ 141 long count[10]; /* count[0]=rays [1]=verts. [2]=base [3]=pivots */ 142 /* count[4]=integer vertices */ 143 long deepest; /* max depth ever reached in search */ 144 long nredundcol; /* number of redundant columns */ 145 long nlinearity; /* number of input linearities */ 146 long totalnodes; /* count total number of tree nodes evaluated */ 147 long runs; /* probes for estimate function */ 148 long seed; /* seed for random number generator */ 149 double cest[10]; /* ests: 0=rays,1=vert,2=bases,3=vol,4=int vert */ 150 /**** flags ********** */ 151 long allbases; /* TRUE if all bases should be printed */ 152 long bound; /* TRUE if upper/lower bound on objective given */ 153 long debug; 154 long dualdeg; /* TRUE if start dictionary is dual degenerate */ 155 long etrace; /* turn off debug at basis # strace */ 156 long frequency; /* frequency to print cobasis indices */ 157 long geometric; /* TRUE if incident vertex prints after each ray */ 158 long getvolume; /* do volume calculation */ 159 long givenstart; /* TRUE if a starting cobasis is given */ 160 long homogeneous; /* TRUE if all entries in column one are zero */ 161 long hull; /* do convex hull computation if TRUE */ 162 long incidence; /* print all tight inequalities (vertices/rays) */ 163 long lponly; /* true if only lp solution wanted */ 164 long maxdepth; /* max depth to search to in treee */ 165 long maximize; /* flag for LP maximization */ 166 long maxoutput; /* if positive, maximum number of output lines */ 167 long minimize; /* flag for LP minimization */ 168 long mindepth; /* do not backtrack above mindepth */ 169 long nash; /* TRUE for computing nash equilibria */ 170 long nonnegative; /* TRUE if last d constraints are nonnegativity */ 171 long polytope; /* TRUE for facet computation of a polytope */ 172 long printcobasis; /* TRUE if all cobasis should be printed */ 173 long printslack; /* TRUE if indices of slack inequal. printed */ 174 long truncate; /* TRUE: truncate tree when moving from opt vert*/ 175 long verbose; /* FALSE for minimalist output */ 176 long restart; /* TRUE if restarting from some cobasis */ 177 long strace; /* turn on debug at basis # strace */ 178 long voronoi; /* compute voronoi vertices by transformation */ 179 180 /* Variables for saving/restoring cobasis, db */ 181 182 long id; /* numbered sequentially */ 183 char *name; /* passed by user */ 184 185 long saved_count[3]; /* How often to print out current cobasis */ 186 long *saved_C; 187 lrs_mp saved_det; 188 long saved_depth; 189 long saved_d; 190 191 long saved_flag; /* There is something in the saved cobasis */ 192 193 /* Variables for cacheing dictionaries, db */ 194 lrs_dic *Qhead, *Qtail; 195 196 } 197 lrs_dat, lrs_dat_p; 198 199 /*********************/ 200 /*global variables */ 201 /*********************/ 202 #define MAX_LRS_GLOBALS 100L 203 #define MAXIMIZE 1L /* maximize the lp */ 204 #define MINIMIZE 0L /* maximize the lp */ 205 #define GE 1L /* constraint is >= */ 206 #define EQ 0L /* constraint is linearity */ 207 208 209 extern FILE *lrs_cfp; /* output file for checkpoint information */ 210 211 extern unsigned long dict_count, dict_limit, cache_tries, cache_misses; 212 extern lrs_dic *PBnew; /* we will save Bob's dictionary in getabasis */ 213 214 215 /*******************************/ 216 /* functions for external use */ 217 /*******************************/ 218 long lrs_main (int argc, char *argv[]); /* lrs driver, argv[1]=input file, [argc-1]=output file */ 219 long redund_main (int argc, char *argv[]); /* redund driver, argv[1]=input file, [2]=output file */ 220 221 lrs_dat *lrs_alloc_dat (const char *name); /* allocate for lrs_dat structure "name" */ 222 lrs_dic *lrs_alloc_dic (lrs_dat * Q); /* allocate for lrs_dic structure corr. to Q */ 223 224 void lrs_estimate (lrs_dic * P, lrs_dat * Q); /* get estimates only */ 225 226 long lrs_read_dat (lrs_dat * Q, int argc, char *argv[]); /* read header and set up lrs_dat */ 227 long lrs_read_dic (lrs_dic * P, lrs_dat * Q); /* read input and set up problem and lrs_dic */ 228 229 long lrs_checkbound (lrs_dic *P, lrs_dat * Q); /* TRUE if current objective value exceeds specified bound */ 230 231 long lrs_getfirstbasis (lrs_dic ** P_p, lrs_dat * Q, lrs_mp_matrix * Lin,long no_output); 232 /* gets first basis, FALSE if none,P may get changed if lin. space Lin found */ 233 /* no_output is TRUE supresses output headers */ 234 235 /* P may get changed if lin. space Lin found */ 236 void lrs_getinput(lrs_dic *P,lrs_dat *Q,long *num,long *den, long m, long d); 237 /* reads input matrix b A in lrs/cdd format */ 238 long lrs_getnextbasis (lrs_dic ** dict_p, lrs_dat * Q, long prune); /* gets next lrs tree basis, FALSE if none */ 239 /* backtrack if prune is TRUE */ 240 241 long lrs_getsolution (lrs_dic * P, lrs_dat * Q, lrs_mp_vector output, long col); 242 long lrs_getray (lrs_dic * P, lrs_dat * Q, long col, long comment, lrs_mp_vector output); 243 long lrs_getvertex (lrs_dic * P, lrs_dat * Q, lrs_mp_vector output); 244 245 void lrs_close (char *name); /* close lrs lib program "name" */ 246 long lrs_init (char *name); /* initialize lrslib and arithmetic package for prog "name" */ 247 long lrs_init_quiet (FILE* in, FILE* out); /* initialize lrslib and arithmetic package for prog "name" */ 248 249 250 void lrs_lpoutput(lrs_dic * P,lrs_dat * Q, lrs_mp_vector output); /* print LP primal and dual solutions */ 251 void lrs_printcobasis (lrs_dic * P, lrs_dat * Q, long col); /* print cobasis for column col(verted or ray) */ 252 void lrs_printoutput (lrs_dat * Q, lrs_mp_vector output); /* print output array */ 253 void lrs_printrow (char name[], lrs_dat * Q, lrs_mp_vector output, long rowd); 254 /*print row of A matrix in output[0..rowd] */ 255 void lrs_printsol (lrs_dic * P, lrs_dat * Q, long col, long comment); /* print out solution from col, comment= */ 256 /* 0=normal,-1=geometric ray,1..inputd=linearity */ 257 void lrs_printtotals (lrs_dic * P, lrs_dat * Q); /* print final totals for lrs */ 258 long lrs_set_digits (long dec_digits ); /* set lrsmp digits to equiv. of decimal dec_digits */ 259 long lrs_solvelp (lrs_dic * P, lrs_dat * Q, long maximize); /* solve primal feas LP:TRUE bounded else FALSE */ 260 261 262 /*******************************/ 263 /* functions for internal use */ 264 /*******************************/ 265 /* basic dictionary functions */ 266 /*******************************/ 267 268 long getabasis (lrs_dic * P, lrs_dat * Q, long order[]); /* Try to find a starting basis */ 269 void getnextoutput (lrs_dic * P, lrs_dat * Q, long i, long col, lrs_mp out); /* get A[B[i][col] and copy to out */ 270 long ismin (lrs_dic * P, lrs_dat * Q, long r, long s); /* test if A[r][s] is a min ratio for col s */ 271 long lexmin (lrs_dic * P, lrs_dat * Q, long col); /* test A to see if current basis is lexmin */ 272 void pivot (lrs_dic * P, lrs_dat * Q, long bas, long cob); /* Qpivot routine for array A */ 273 long primalfeasible (lrs_dic * P, lrs_dat * Q); /* Do dual pivots to get primal feasibility */ 274 long ratio (lrs_dic * P, lrs_dat * Q, long col); /* find lex min. ratio */ 275 long removecobasicindex (lrs_dic * P, lrs_dat * Q, long k); /* remove C[k] from problem */ 276 long restartpivots (lrs_dic * P, lrs_dat * Q); /* restart problem from given cobasis */ 277 long reverse (lrs_dic * P, lrs_dat * Q, long *r, long s); /* TRUE if B[*r] C[s] is a reverse lex-pos pivot */ 278 long selectpivot (lrs_dic * P, lrs_dat * Q, long *r, long *s); /* select pivot indices using lexicographic rule */ 279 long dan_selectpivot (lrs_dic * P, lrs_dat * Q, long *r, long *s); /* select pivot indices using dantzig-lex rule */ 280 void update (lrs_dic * P, lrs_dat * Q, long *i, long *j); /* update the B,C, LOC arrays after a pivot */ 281 void updatevolume (lrs_dic * P, lrs_dat * Q); /* rescale determinant and update the volume */ 282 283 /*******************************/ 284 /* other functions using P,Q */ 285 /*******************************/ 286 287 long lrs_degenerate (lrs_dic * P, lrs_dat * Q); /* TRUE if the dictionary is primal degenerate */ 288 289 void print_basis (FILE * fp, lrs_dat * Q); 290 void printA (lrs_dic * P, lrs_dat * Q); /* raw print of dictionary, bases for debugging */ 291 void pimat (lrs_dic * P, long r, long s, lrs_mp Nt, char name[]); /* print the row r col s of A */ 292 293 long readfacets (lrs_dat * Q, long facet[]); /* read and check facet list */ 294 long readlinearity (lrs_dat * Q); /* read and check linearity list */ 295 void rescaledet (lrs_dic * P, lrs_dat * Q, lrs_mp Vnum, lrs_mp Vden); /* rescale determinant to get its volume */ 296 void rescalevolume (lrs_dic * P, lrs_dat * Q, lrs_mp Vnum, lrs_mp Vden); /* adjust volume for dimension */ 297 298 /***************************************************/ 299 /* Routines for redundancy checking */ 300 /***************************************************/ 301 302 303 long checkredund (lrs_dic * P, lrs_dat * Q); /* solve primal lp to check redund of obj fun. */ 304 /* returns TRUE if redundant, else FALSE */ 305 306 long checkcobasic (lrs_dic * P, lrs_dat * Q, long index); /* TRUE if index is cobasic and nondegenerate */ 307 /* FALSE if basic, or degen. cobasic, where it will get pivoted out */ 308 309 long checkindex (lrs_dic * P, lrs_dat * Q, long index); /* index=0 non-red.,1 red., 2 input linearity */ 310 /* NOTE: row is returned all zero if redundant!! */ 311 312 /***************************************************/ 313 /* Routines for caching and restoring dictionaries */ 314 /***************************************************/ 315 316 void lrs_free_dic ( lrs_dic *P, lrs_dat *Q); 317 void lrs_free_dat ( lrs_dat *Q); 318 void cache_dict (lrs_dic ** D_p, lrs_dat * global, long i, long j); 319 long check_cache (lrs_dic ** D_p, lrs_dat * global, long *i_p, long *j_p); 320 void copy_dict (lrs_dat * global, lrs_dic * dest, lrs_dic * src); 321 void pushQ (lrs_dat * global, long m, long d, long m_A); 322 323 void save_basis (lrs_dic * D, lrs_dat * Q); 324 lrs_dic *alloc_memory (lrs_dat * Q); 325 lrs_dic * lrs_getdic(lrs_dat *Q); 326 lrs_dic *new_lrs_dic (long m, long d, long m_A); 327 lrs_dic *resize (lrs_dic * P, lrs_dat * Q); 328 329 330 /*******************************/ 331 /* signals handling */ 332 /*******************************/ 333 334 void checkpoint (); 335 void die_gracefully (); 336 void setup_signals (); 337 void timecheck (); 338 339 /*******************************/ 340 /* utilities */ 341 /*******************************/ 342 void lprat (const char *name, long Num, long Den); /* Print Num/Den without reducing */ 343 long lreadrat (long *Num, long *Den); /* read a rational string and convert to long integers */ 344 void reorder (long a[], long range); /* reorder array in increasing order with one misplaced element */ 345 void reorder1 (long a[], long b[], long newone, long range); 346 /* reorder array a in increasing order with misplaced element newone */ 347 /* elements of b go along for the ride */ 348 349 /***************************/ 350 /* lp_solve like functions */ 351 /***************************/ 352 long lrs_solve_lp(lrs_dic *P, lrs_dat *Q); /* solve lp only for given dictionary */ 353 354 void lrs_set_row(lrs_dic *P, lrs_dat *Q, long row, long num[], long den[], long ineq); 355 /* load row i of dictionary from num[]/den[] ineq=GE */ 356 void lrs_set_row_mp(lrs_dic *P, lrs_dat *Q, long row, lrs_mp_vector num, lrs_mp_vector den, long ineq); 357 /* same as lrs_set_row except num/den is lrs_mp type */ 358 359 void lrs_set_obj(lrs_dic *P, lrs_dat *Q, long num[], long den[], long max); 360 /* set up objective function with coeffs num[]/den[] max=MAXIMIZE or MINIMIZE */ 361 void lrs_set_obj_mp(lrs_dic *P, lrs_dat *Q, lrs_mp_vector num, lrs_mp_vector den, long max); 362 /* same as lrs_set_obj but num/den has lrs_mp type */ 363 364 365 366 /* end of lrslib.h (vertex enumeration using lexicographic reverse search) */ 367 #ifdef __cplusplus 368 } 369 #endif 370