1 /* lrslib.c library code for lrs */
2
3
4 /* modified by Gary Roumanis for multithread plrs compatability */
5 /* truncate needs mod to supress last pivot */
6 /* need to add a test for non-degenerate pivot step in reverse I guess */
7 /* Copyright: David Avis 2005,2011 avis@cs.mcgill.ca */
8
9 /* This program is free software; you can redistribute it and/or modify
10 it under the terms of the GNU General Public License as published by
11 the Free Software Foundation; either version 2 of the License, or
12 (at your option) any later version.
13
14 This program is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 GNU General Public License for more details.
18
19 You should have received a copy of the GNU General Public License
20 along with this program; if not, write to the Free Software
21 Foundation, Inc., 51 Franklin Street, Suite 500, Boston, MA 02110-1335, USA.
22 */
23
24
25
26 #include <stdio.h>
27 #include <string.h>
28 #include <setjmp.h>
29 #include <limits.h>
30 #include "lrsrestart.h"
31 #include "lrslib.h"
32
33 static unsigned long dict_count, dict_limit, cache_tries, cache_misses;
34
35 /* Variables and functions global to this file only */
36
37 static long lrs_checkpoint_seconds = 0;
38
39 static long lrs_global_count = 0; /* Track how many lrs_dat records are
40 allocated */
41 static size_t infileLen; /* length of cache of input file */
42 static char *infile = NULL; /* cache of input for restart */
43 static char infilename[PATH_MAX];
44 static char outfilename[PATH_MAX];
45 static char tmpfilename[PATH_MAX];
46 static int tmpfd;
47 static long overflow=0; /* =0 no overflow =1 restart overwrite =2 restart append */
48 static long pivoting=FALSE; /* =0 no overflow =1 restart overwrite =2 restart append */
49
50 static jmp_buf buf1;
51
52 static lrs_dat_p *lrs_global_list[MAX_LRS_GLOBALS + 1];
53
54 static lrs_dic *new_lrs_dic (long m, long d, long m_A);
55
56
57 static void cache_dict (lrs_dic ** D_p, lrs_dat * global, long i, long j);
58 static long check_cache (lrs_dic ** D_p, lrs_dat * global, long *i_p, long *j_p);
59 static void save_basis (lrs_dic * D, lrs_dat * Q);
60
61 static void lrs_dump_state ();
62
63 static void pushQ (lrs_dat * global, long m, long d, long m_A);
64
65 #ifndef TIMES
66 static void ptimes (void);
67 static double get_time(void);
68 #endif
69
70
71 /*******************************/
72 /* signals handling */
73 /*******************************/
74 #ifndef SIGNALS
75 static void checkpoint ();
76 static void die_gracefully ();
77 static void setup_signals (void);
78 static void timecheck ();
79 #endif
80
81 /*******************************/
82 /* functions for external use */
83 /*******************************/
84
85 /******************************************************************/
86 /* lrs_run is the main reverse search part of lrs */
87 /* should be called by lrsv2_main which does setup and close also */
88 /******************************************************************/
89 long
lrs_run(lrs_dic * P,lrs_dat * Q)90 lrs_run ( lrs_dic *P, lrs_dat * Q)
91
92 {
93
94
95 lrs_mp_matrix Lin; /* holds input linearities if any are found */
96 long col; /* output column index for dictionary */
97 long startcol = 0;
98 long prune = FALSE; /* if TRUE, getnextbasis will prune tree and backtrack */
99
100
101 /*********************************************************************************/
102 /* Step 1: Find a starting cobasis from default of specified order */
103 /* P is created to hold active dictionary data and may be cached */
104 /* Lin is created if necessary to hold linearity space */
105 /* Print linearity space if any, and retrieve output from first dict. */
106 /*********************************************************************************/
107
108 if (!lrs_getfirstbasis (&P, Q, &Lin, FALSE))
109 return 1;
110
111 /* Pivot to a starting dictionary */
112 /* There may have been column redundancy */
113 /* If so the linearity space is obtained and redundant */
114 /* columns are removed. User can access linearity space */
115 /* from lrs_mp_matrix Lin dimensions nredundcol x d+1 */
116
117
118
119 if (Q->homogeneous && Q->hull)
120 startcol++; /* col zero not treated as redundant */
121
122 if(!Q->restart)
123 for (col = startcol; col < Q->nredundcol; col++) /* print linearity space */
124 lrs_printoutput (Q, Lin[col]); /* Array Lin[][] holds the coeffs. */
125
126
127 if(Q->nredundcol > 0)
128 lrs_clear_mp_matrix(Lin,Q->nredundcol,Q->n);
129
130
131 /*********************************************************************************/
132 /* Step 3: Terminate if lponly option set, otherwise initiate a reverse */
133 /* search from the starting dictionary. Get output for each new dict. */
134 /*********************************************************************************/
135
136
137
138 /* We initiate reverse search from this dictionary */
139 /* getting new dictionaries until the search is complete */
140 /* User can access each output line from output which is */
141 /* vertex/ray/facet from the lrs_mp_vector output */
142 /* prune is TRUE if tree should be pruned at current node */
143 do
144 {
145
146 //2015.6.5 after maxcobases reached, generate subtrees that have not been enumerated
147 //2018.1.19 fix printcobasis bug when maxcobases set
148 //2019.5.8 new givoutput flag to avoid printing restart cobases
149
150 prune=lrs_checkbound(P,Q);
151
152 if (!prune && Q->giveoutput)
153 {
154 lrs_open_outputblock(); /* keeps output together when using mplrs */
155
156 for (col = 0; col <= P->d; col++) /* print output if any */
157 if (lrs_getsolution (P, Q, Q->output, col))
158 lrs_printoutput (Q, Q->output);
159
160 lrs_close_outputblock();
161 }
162 else
163 Q->giveoutput=TRUE; /* first output supressed for restart */
164
165 /*2020.3.9 bounds on objective function check corrected */
166
167
168 if ((Q->maxcobases > 0) && (Q->count[2] >=Q->maxcobases))
169 {
170 prune=TRUE;
171 if( !lrs_leaf(P,Q)) /* do not return cobases of a leaf */
172 lrs_return_unexplored(P,Q);
173 }
174
175 save_basis(P,Q);
176
177 }while (!Q->lponly && lrs_getnextbasis (&P, Q, prune)); // do ...
178
179 if (Q->lponly)
180 lrs_lpoutput(P,Q,Q->output);
181 else
182 lrs_printtotals (P, Q); /* print final totals, including estimates */
183
184 Q->m=P->m;
185 lrs_free_dic(P,Q); /* note Q is not free here and can be reused */
186
187 return 0;
188 }
189 /*********************************************/
190 /* end of model test program for lrs library */
191 /*********************************************/
192
193
194 /*******************************************************/
195 /* redund_run is main loop for redundancy removal */
196 /*******************************************************/
197 long
redund_run(lrs_dic * P,lrs_dat * Q)198 redund_run ( lrs_dic *P, lrs_dat * Q)
199
200 {
201 lrs_mp_matrix Ain; /* holds a copy of the input matrix to output at the end */
202
203 long ineq; /* input inequality number of current index */
204 long *redineq;
205
206 lrs_mp_matrix Lin; /* holds input linearities if any are found */
207
208 long i, j, d, m;
209 long nlinearity; /* number of linearities in input file */
210 long lastdv;
211 long debug;
212 long index; /* basic index for redundancy test */
213
214 /*********************************************************************************/
215
216
217 /* if non-negative flag is set, non-negative constraints are not input */
218 /* explicitly, and are not checked for redundancy */
219
220
221 m = P->m_A; /* number of rows of A matrix */
222 d = P->d;
223 redineq = Q->redineq;
224 debug = Q->debug;
225
226 Q->Ain = lrs_alloc_mp_matrix (m, d); /* make a copy of A matrix for output later */
227 Ain=Q->Ain;
228
229 for (i = 1; i <= m; i++)
230 {
231 for (j = 0; j <= d; j++)
232 copy (Ain[i][j], P->A[i][j]);
233
234 if (debug)
235 lrs_printrow ("*", Q, Ain[i], d);
236 }
237
238
239 /*********************************************************************************/
240 /* Step 1: Find a starting cobasis from default of specified order */
241 /* Lin is created if necessary to hold linearity space */
242 /*********************************************************************************/
243
244 if (!lrs_getfirstbasis (&P, Q, &Lin, TRUE))
245 return 1;
246
247 /* Pivot to a starting dictionary */
248 /* There may have been column redundancy */
249 /* If so the linearity space is obtained and redundant */
250 /* columns are removed. User can access linearity space */
251 /* from lrs_mp_matrix Lin dimensions nredundcol x d+1 */
252
253
254 /*********************************************************************************/
255 /* Step 2: Test rows i where redineq[i]=1 for redundancy */
256 /*********************************************************************************/
257
258 /* note some of these may have been changed in getting initial dictionary */
259 m = P->m_A;
260 d = P->d;
261 nlinearity = Q->nlinearity;
262 lastdv = Q->lastdv;
263
264 /* linearities are not considered for redundancy */
265
266 for (i = 0; i < nlinearity; i++)
267 redineq[Q->linearity[i]] = 2L;
268
269 if(Q->debug)
270 fprintf (lrs_ofp, "\nredundcheck=%ld verifyredund=%ld",Q->noredundcheck, Q->verifyredund);
271
272
273 /* Q->verifyredund always false in lrs, set by mplrs to check duplicated redundancy removal */
274 /* Q->noredundcheck overides this to skip verification */
275
276 if(Q->noredundcheck && Q->verifyredund)
277 goto done;
278
279 /* mplrs sets redineq[i]==-1 for guaranteed redundant inequalities */
280 /* these rows must be zeroed out before testing the others */
281
282
283
284 if (Q->verifyredund) /* this is never run by lrs, final step of mplrs redund */
285 for (index = lastdv + Q->redineq[0]; index <= m + d; index++)
286 {
287 ineq = Q->inequality[index - lastdv]; /* the input inequality number corr. to this index */
288 if( redineq[ineq]== -1 )
289 checkindex (P, Q, -index); /* used to zero correct row of A no LP solved */
290 }
291
292
293 /* rows 0..lastdv are cost, decision variables, or linearities */
294 /* other rows need to be tested */
295
296 for (index = lastdv + Q->redineq[0]; index <= m + d; index++)
297 {
298 ineq = Q->inequality[index - lastdv]; /* the input inequality number corr. to this index */
299 Q->redineq[0] = ineq; /* used for restarting after arithmetic overflow */
300
301 if( redineq[ineq]==1 )
302 {
303 redineq[ineq] = checkindex (P, Q, index);
304
305 if (debug)
306 fprintf (lrs_ofp, "\ncheck index=%ld, inequality=%ld, redineq=%ld", index, ineq, redineq[ineq]);
307 if(!Q->mplrs && Q->verbose)
308 {
309 if( redineq[ineq]==1 )
310 lrs_printrow ("*re ", Q, Ain[ineq], Q->inputd);
311 else
312 lrs_printrow ("*nr ", Q, Ain[ineq], Q->inputd);
313 }
314 }
315
316 } /* end for index ..... */
317
318 done:
319
320
321
322 if(Q->mplrs && !Q->verifyredund)
323 { /* return array redineq to consumer */
324 char *ss;
325 int len=0;
326 ss=(char *)malloc(20*m*sizeof(char));
327
328 for (i=1; i<=m; i++)
329 if(redineq[i]==1)
330 len=len+sprintf(ss+len," %ld",i);
331 if(len>0)
332 lrs_post_output("redund", ss);
333
334 free(ss);
335 lrs_clear_mp_matrix(Ain,P->m_A,P->d);
336 Q->m=P->m;
337 lrs_free_dic(P,Q); /* note Q is not free here and can be reused */
338
339 return 0;
340 }
341
342
343 if (Q->verbose || Q->debug)
344 {
345 fprintf (lrs_ofp, "\n*redineq:");
346 for (i = 1; i <= m; i++)
347 fprintf (lrs_ofp, " %ld", redineq[i]);
348 }
349
350 redund_print(Ain,P,Q);
351
352 lrs_clear_mp_matrix(Ain,P->m_A,P->d);
353 Q->m=P->m;
354 lrs_free_dic(P,Q); /* note Q is not free here and can be reused */
355
356
357 return 0;
358 }
359 /*********************************************/
360
redund_print(lrs_mp_matrix Ain,lrs_dic * P,lrs_dat * Q)361 void redund_print(lrs_mp_matrix Ain,lrs_dic *P,lrs_dat *Q)
362 {
363 long i, m;
364 long nlinearity; /* number of linearities in input file */
365 long nredund; /* number of redundant rows in input file */
366 long *redineq=Q->redineq;
367
368 m = P->m_A; /* number of rows of A matrix */
369 nlinearity = Q->nlinearity;
370
371 /* restore as mplrs loses this */
372 for (i = 0; i < nlinearity; i++)
373 redineq[Q->linearity[i]]=2;
374
375 /*
376 fprintf(lrs_ofp,"\nQ->red");
377 for (i = 1; i <= m; i++)
378 fprintf(lrs_ofp," %ld",Q->redineq[i]);
379 */
380
381
382 if (!Q->hull)
383 fprintf (lrs_ofp, "\nH-representation");
384 else
385 fprintf (lrs_ofp, "\nV-representation");
386
387 /* linearities will be printed first in output */
388
389 if (nlinearity > 0)
390 {
391 fprintf (lrs_ofp, "\nlinearity %ld", nlinearity);
392 for (i = 1; i <= nlinearity; i++)
393 fprintf (lrs_ofp, " %ld", i);
394
395 }
396
397 nredund = 0; /* count number of non-redundant inequalities */
398
399 for (i = 1; i <= m; i++)
400 if (redineq[i] == 0)
401 nredund++;
402
403 fprintf (lrs_ofp, "\nbegin");
404 fprintf (lrs_ofp, "\n%ld %ld rational", nlinearity+nredund, Q->n);
405
406 pivoting=TRUE;
407 /* print the linearities first */
408
409 for (i = 0; i < nlinearity; i++)
410 lrs_printrow ("", Q, Ain[Q->linearity[i]], Q->inputd);
411
412 for (i = 1; i <= m; i++)
413 if (redineq[i] == 0)
414 lrs_printrow ("", Q, Ain[i], Q->inputd);
415
416 fprintf (lrs_ofp, "\nend");
417 fprintf (lrs_ofp, "\n*Input had %ld rows and %ld columns", m, Q->n);
418
419 if( m==nredund)
420 fprintf (lrs_ofp, "\n*No redundant rows found");
421 else
422 {
423 fprintf (lrs_ofp, "\n* %ld redundant row(s) found:\n", m - nredund-nlinearity);
424 for (i=1; i<=m; i++)
425 if(redineq[i]==1 || redineq[i]==-1)
426 fprintf (lrs_ofp, " %ld",i);
427 if (Q->noredundcheck)
428 fprintf (lrs_ofp, "\n*Warning: not verified - input should be full dimensional and duplicate free");
429 }
430 fprintf (lrs_ofp, "\n");
431 return;
432 } /* end of redund_print */
433
434
435 /*******************/
436 /* lrs_printoutput */
437 /* one line only */
438 /*******************/
439 void
lrs_printoutput(lrs_dat * Q,lrs_mp_vector output)440 lrs_printoutput (lrs_dat * Q, lrs_mp_vector output)
441 {
442 char *sss;
443 char **ss;
444
445 long i;
446 long len=0;
447
448 if (Q->countonly)
449 return;
450
451 ss = (char **)malloc((1+Q->n) * sizeof(char*));
452
453 if (Q->hull || zero (output[0])) /*non vertex */
454 for (i = 0; i < Q->n; i++)
455 {
456 ss[i]=cpmp ("", output[i]);
457 len=len+snprintf(NULL, 0, "%s ", ss[i] );
458 }
459 else
460 for (i = 1; i < Q->n; i++)
461 {
462 ss[i]=cprat("", output[i], output[0]);
463 len=len+snprintf(NULL, 0, "%s ", ss[i] );
464 }
465
466 sss=(char*)malloc((len+5)* sizeof(char*));
467 len=0;
468
469 if (Q->hull || zero (output[0])) /*non vertex */
470 for (i = 0; i < Q->n; i++)
471 {
472 len=len+sprintf(sss+len,"%s ",ss[i]);
473 free(ss[i]);
474 }
475 else
476 { /* vertex */
477 len=sprintf (sss, " 1 ");
478 for (i = 1; i < Q->n; i++)
479 {
480 len=len+sprintf(sss+len, "%s ", ss[i] );
481 free(ss[i]);
482 }
483 }
484
485 if(Q->mplrs)
486 lrs_post_output("vertex",sss);
487 else
488 fprintf (lrs_ofp, "\n%s",sss);
489
490 free(ss);
491 free(sss);
492
493 }
494 /**************************/
495 /* end of lrs_printoutput */
496 /**************************/
497
498 /****************/
499 /* lrs_lpoutput */
500 /****************/
lrs_lpoutput(lrs_dic * P,lrs_dat * Q,lrs_mp_vector output)501 void lrs_lpoutput(lrs_dic * P,lrs_dat * Q, lrs_mp_vector output)
502 {
503
504 if(Q->unbounded || !Q->messages)
505 return;
506
507 lrs_mp Temp1, Temp2;
508 long i;
509
510 lrs_alloc_mp (Temp1);
511 lrs_alloc_mp (Temp2);
512
513 prat ("\n*Obj=",P->objnum, P->objden);
514 fprintf (lrs_ofp, " pivots=%ld ",Q->count[3]);
515 if(Q->verbose)
516 {
517 fprintf (lrs_ofp, "\n\n*Primal: ");
518 for (i = 1; i < Q->n; i++)
519 {
520 fprintf(lrs_ofp,"x_%ld=",i);
521 prat ("", output[i], output[0]);
522 }
523 if(Q->nlinearity > 0)
524 fprintf (lrs_ofp, "\n\n*Linearities in input file - partial dual solution only");
525 fprintf (lrs_ofp, "\n\n*Dual: ");
526
527 for (i = 0; i < P->d; i++)
528 {
529 fprintf(lrs_ofp,"y_%ld=",Q->inequality[P->C[i]-Q->lastdv]);
530 changesign(P->A[0][P->Col[i]]);
531 mulint(Q->Lcm[P->Col[i]],P->A[0][P->Col[i]],Temp1);
532 mulint(Q->Gcd[P->Col[i]],P->det,Temp2);
533 prat("",Temp1,Temp2);
534 changesign(P->A[0][P->Col[i]]);
535 }
536 }
537 fprintf (lrs_ofp, "\n");
538 lrs_clear_mp (Temp1);
539 lrs_clear_mp (Temp2);
540 }
541 /***********************/
542 /* end of lrs_lpoutput */
543 /***********************/
544 void
lrs_printrow(const char * name,lrs_dat * Q,lrs_mp_vector output,long rowd)545 lrs_printrow (const char *name, lrs_dat * Q, lrs_mp_vector output, long rowd)
546 /* print a row of A matrix in output in "original" form */
547 /* rowd+1 is the dimension of output vector */
548 /* if input is H-rep. output[0] contains the RHS */
549 /* if input is V-rep. vertices are scaled by 1/output[1] */
550 {
551 long i;
552 fprintf (lrs_ofp, "\n%s", name);
553 if (!Q->hull) /* input was inequalities, print directly */
554
555 {
556 for (i = 0; i <= rowd; i++)
557 pmp ("", output[i]);
558 return;
559 }
560
561 /* input was vertex/ray */
562
563 if (zero (output[1])) /*non-vertex */
564 {
565 for (i = 1; i <= rowd; i++)
566 pmp ("", output[i]);
567 }
568 else
569 { /* vertex */
570 fprintf (lrs_ofp, " 1 ");
571 for (i = 2; i <= rowd; i++)
572 prat ("", output[i], output[1]);
573 }
574
575 return;
576
577 } /* end of lrs_printrow */
578
579 long
lrs_getsolution(lrs_dic * P,lrs_dat * Q,lrs_mp_vector output,long col)580 lrs_getsolution (lrs_dic * P, lrs_dat * Q, lrs_mp_vector output, long col)
581 /* check if column indexed by col in this dictionary */
582 /* contains output */
583 /* col=0 for vertex 1....d for ray/facet */
584 {
585
586
587 long j; /* cobasic index */
588
589 lrs_mp_matrix A = P->A;
590 long *Row = P->Row;
591
592 if (col == ZERO) /* check for lexmin vertex */
593 return lrs_getvertex (P, Q, output);
594
595 /* check for rays: negative in row 0 , positive if lponly */
596
597 if (Q->lponly)
598 {
599 if (!positive (A[0][col]))
600 return FALSE;
601 }
602
603 else if (!negative (A[0][col]))
604 return FALSE;
605
606 /* and non-negative for all basic non decision variables */
607
608 j = Q->lastdv + 1;
609 while (j <= P->m && !negative (A[Row[j]][col]))
610 j++;
611
612 if (j <= P->m)
613 return FALSE;
614
615 if (Q->geometric || Q->allbases || lexmin (P, Q, col) || Q->lponly)
616 return lrs_getray (P, Q, col, Q->n, output);
617
618 return FALSE; /* no more output in this dictionary */
619
620 } /* end of lrs_getsolution */
621
622 void
lrs_print_header(const char * name)623 lrs_print_header(const char *name)
624 {
625 if(lrs_ofp == NULL)
626 lrs_ofp=stdout;
627 #ifdef LRS_QUIET
628 return;
629 #endif
630 fprintf (lrs_ofp,"%s", name);
631 fprintf (lrs_ofp,TITLE);
632 fprintf (lrs_ofp,VERSION);
633 fprintf (lrs_ofp,"(");
634 fprintf (lrs_ofp,BIT);
635 fprintf (lrs_ofp,",");
636 fprintf (lrs_ofp,ARITH);
637 #ifdef MA
638 fprintf (lrs_ofp,",hybrid arithmetic");
639 #endif
640 #ifdef LRSLONG
641 #ifndef SAFE
642 fprintf (lrs_ofp,",no overflow checking");
643 #endif
644 #endif
645 fprintf (lrs_ofp,")");
646 if(overflow != 2)
647 {
648 #ifdef GMP
649 fprintf(lrs_ofp," gmp v.%d.%d",__GNU_MP_VERSION,__GNU_MP_VERSION_MINOR);
650 #elif defined(FLINT)
651 fprintf(lrs_ofp," %dbit flint v.%s", FLINT_BITS, FLINT_VERSION);
652 #endif
653 }
654 }
655
656 long
lrs_init(const char * name)657 lrs_init (const char *name) /* returns TRUE if successful, else FALSE */
658 {
659
660 #ifndef PLRS
661 #ifndef LRS_QUIET
662 if(overflow!=2)
663 lrs_print_header(name);
664 #endif
665 #endif
666
667 if (!lrs_mp_init (ZERO, stdin, stdout)) /* initialize arithmetic */
668 return FALSE;
669
670 lrs_global_count = 0;
671 lrs_checkpoint_seconds = 0;
672 #ifndef SIGNALS
673 setup_signals ();
674 #endif
675 return TRUE;
676 }
677
678 void
lrs_close(const char * name)679 lrs_close (const char *name)
680 {
681
682 #ifdef PLRS
683 return;
684 #endif
685
686 #ifdef LRS_QUIET
687 fprintf (lrs_ofp, "\n");
688 if (lrs_ofp != stdout)
689 {
690 fclose (lrs_ofp);
691 lrs_ofp=NULL;
692 }
693 return;
694 #endif
695
696 #ifdef LRSLONG
697 #ifdef SAFE
698 fprintf (lrs_ofp, "\n*Overflow checking on lrslong arithmetic");
699 #else
700 fprintf (lrs_ofp, "\n*Caution: no overflow checking on long integer arithemtic");
701 #endif
702 #endif
703
704 fprintf (lrs_ofp, "\n*%s", name);
705 fprintf (lrs_ofp, TITLE);
706 fprintf (lrs_ofp, VERSION);
707 fprintf (lrs_ofp, "(");
708 fprintf (lrs_ofp, BIT);
709 fprintf (lrs_ofp, ",");
710 fprintf (lrs_ofp, ARITH);
711 #ifdef MA
712 fprintf (lrs_ofp, ",hybrid arithmetic");
713 #endif
714 fprintf (lrs_ofp, ")");
715
716 #ifdef MP
717 fprintf (lrs_ofp, " max digits=%ld/%ld", DIG2DEC (lrs_record_digits), DIG2DEC (lrs_digits));
718 #endif
719
720 #ifndef TIMES
721 ptimes ();
722 #endif
723
724 if (lrs_ofp != stdout)
725 {
726 fclose (lrs_ofp);
727 lrs_ofp=NULL;
728 }
729 }
730
731 /***********************************/
732 /* allocate and initialize lrs_dat */
733 /***********************************/
734 lrs_dat *
lrs_alloc_dat(const char * name)735 lrs_alloc_dat (const char *name)
736 {
737 lrs_dat *Q;
738 long i;
739
740
741 if (lrs_global_count >= MAX_LRS_GLOBALS)
742 {
743 fprintf (stderr,
744 "Fatal: Attempt to allocate more than %ld global data blocks\n", MAX_LRS_GLOBALS);
745 return NULL;
746
747 }
748
749 Q = (lrs_dat *) malloc (sizeof (lrs_dat));
750 if (Q == NULL)
751 return Q; /* failure to allocate */
752
753 lrs_global_list[lrs_global_count] = Q;
754 Q->id = lrs_global_count;
755 lrs_global_count++;
756 Q->name=(char *) CALLOC ((unsigned) strlen(name)+1, sizeof (char));
757 strcpy(Q->name,name);
758
759 /* initialize variables */
760 Q->mplrs=FALSE;
761 Q->messages=TRUE;
762 #ifdef PLRS
763 Q->mplrs=TRUE;
764 #endif
765 #ifdef LRS_QUIET
766 Q->messages=FALSE;
767 #endif
768 strcpy(Q->fname,""); /* name of program, filled in later */
769 Q->m = 0L;
770 Q->n = 0L;
771 Q->inputd = 0L;
772 Q->deepest = 0L;
773 Q->nlinearity = 0L;
774 Q->nredundcol = 0L;
775 Q->runs = 0L;
776 Q->subtreesize=MAXD;
777 Q->seed = 1234L;
778 Q->totalnodes = 0L;
779 for (i = 0; i < 10; i++)
780 {
781 Q->count[i] = 0L;
782 Q->cest[i] = 0.0;
783 if(i < 5)
784 Q->startcount[i] = 0L;
785 }
786 Q->count[2] = 1L; /* basis counter */
787 Q->startcount[2] = 0L; /* starting basis counter */
788 /* initialize flags */
789 Q->allbases = FALSE;
790 Q->bound = FALSE; /* upper/lower bound on objective function given */
791 Q->countonly = FALSE; /* produce the usual output */
792 Q->debug = FALSE;
793 Q->frequency = 0L;
794 Q->dualdeg = FALSE; /* TRUE if dual degenerate starting dictionary */
795 Q->geometric = FALSE;
796 Q->getvolume = FALSE;
797 Q->homogeneous = TRUE;
798 Q->polytope = FALSE;
799 Q->triangulation = FALSE;
800 Q->hull = FALSE;
801 Q->incidence = FALSE;
802 Q->lponly = FALSE;
803 Q->maxdepth = MAXD;
804 Q->mindepth = -MAXD;
805 Q->maxoutput = 0L;
806 Q->maxcobases = 0L; /* after maxcobases have been found unexplored subtrees reported */
807 Q->nash = FALSE;
808 Q->nonnegative = FALSE;
809 Q->printcobasis = FALSE;
810 Q->printslack = FALSE;
811 Q->truncate = FALSE; /* truncate tree when moving from opt vertex */
812 Q->extract=FALSE;
813 Q->verbose=FALSE;
814 Q->voronoi = FALSE;
815 Q->maximize = FALSE; /*flag for LP maximization */
816 Q->minimize = FALSE; /*flag for LP minimization */
817 Q->restart = FALSE; /* TRUE if restarting from some cobasis */
818 Q->givenstart = FALSE; /* TRUE if a starting cobasis is given */
819 Q->strace = -1L; /* turn on debug at basis # strace */
820 Q->etrace = -1L; /* turn off debug at basis # etrace */
821 Q->newstart=FALSE;
822 Q->giveoutput=TRUE; /* set to false for first output after restart */
823 Q->verifyredund=FALSE; /* set to true when mplrs verifies redund output */
824 Q->noredundcheck=FALSE; /* set to true when mplrs skips verifying output */
825 Q->nextineq=15; /* start redundancy testing from this row */
826 Q->startcob=NULL;
827
828 Q->saved_flag = 0; /* no cobasis saved initially, db */
829 lrs_alloc_mp (Q->Nvolume);
830 lrs_alloc_mp (Q->Dvolume);
831 lrs_alloc_mp (Q->sumdet);
832 lrs_alloc_mp (Q->saved_det);
833 lrs_alloc_mp (Q->boundn);
834 lrs_alloc_mp (Q->boundd);
835 itomp (ZERO, Q->Nvolume);
836 itomp (ONE, Q->Dvolume);
837 itomp (ZERO, Q->sumdet);
838 Q->unbounded = FALSE;
839
840 return Q;
841 } /* end of allocate and initialize lrs_dat */
842
843
844 /*******************************/
845 /* lrs_read_dat */
846 /*******************************/
847 long
lrs_read_dat(lrs_dat * Q,int argc,char * argv[])848 lrs_read_dat (lrs_dat * Q, int argc, char *argv[])
849 {
850 char name[1000];
851 char writemode[2]="w"; /* will be set to "a" (append) for overflow or newstart */
852 long dec_digits = DEFAULT_DIGITS;
853 long infilenum=0; /*input file number to open if any */
854 long firstline = TRUE; /*flag for picking off name at line 1 */
855 long i;
856 int c; /* for fgetc */
857 int messages = Q->messages; /* print output for each option */
858
859 *tmpfilename='\0';
860 if(overflow==2) /* otherwise overwrite output */
861 strcpy(writemode,"a");
862
863 strcpy(outfilename, "\0");
864
865 if(argc > 1 )
866 {
867 infilenum=1;
868 if(Q->nash && argc ==2) /* legacy code to open second nash input file */
869 infilenum=2;
870 if(Q->nash && argc ==4) /* legacy code for nash output file */
871 strcpy(outfilename,argv[3]);
872 }
873
874 if (infilenum > 0 && (lrs_ifp = fopen (argv[infilenum], "r")) == NULL) /* command line overides stdin */
875 {
876 fprintf (stderr,"\n*bad input file name\n");
877 return (FALSE);
878 }
879
880 if (infilenum==1)
881 {
882 strcpy(infilename,argv[1]);
883 if(!Q->mplrs && messages && overflow == 0 )
884 printf ("\n*Input taken from %s", infilename);
885 fflush(stdout);
886 }
887
888 #ifdef LRSLONG
889 if(infilenum==0) /* stdin gets written to a temporary file */
890 {
891 strcpy(tmpfilename,"/tmp/lrs_stdinXXXXXX");
892 mkstemp(tmpfilename);
893 lrs_stdin_to_file(tmpfilename);
894 lrs_ifp=fopen (tmpfilename, "r");
895 strcpy(infilename,tmpfilename);
896 }
897
898 lrs_file_to_cache(lrs_ifp);
899
900 #endif
901
902 if(argc > 2) /* lrs has commandline arguments for newstart */
903 {
904 if (!Q->nash )
905 {
906 i=2;
907 while (i < argc) /* add command line arguments here */
908 {
909 if(strcmp(argv[i],"-newstart")==0) /* newstart not currently used ... */
910 {
911 strcpy(writemode,"a");
912 Q->newstart=TRUE;
913 }
914 else /* command line argument overides stdout */
915 strcpy(outfilename,argv[i++]);
916 }
917 }
918 if(strcmp(outfilename,"\0") != 0 )
919 {
920 if ((lrs_ofp = fopen (outfilename, writemode)) == NULL)
921 {
922 fprintf (stderr,"\n*bad output file name %s\n",outfilename);
923 return (FALSE);
924 }
925 else
926 if(overflow == 0)
927 printf ("\n*Output sent to file %s\n", outfilename);
928 }
929 }
930
931 /*2020.5.19 new redund handling, thanks to DB */
932 /* symbolic link from redund to lrs needed */
933 /* similar links if lrs1, lrs2 or lrsgmp used */
934 /* any redund option in input will overide */
935
936 if(!Q->mplrs && lrs_ofp != stdout && overflow != 2 ) /* headers for the output file also */
937 {
938 char *name;
939 name=(char *) malloc(strlen(Q->fname)+5);
940 strcpy(name,"*");
941 strcat(name,Q->fname);
942 strcat(name,":");
943 lrs_print_header(name);
944 free(name);
945 }
946
947 /* process input file */
948 if( fscanf (lrs_ifp, "%s", name) == EOF)
949 {
950 fprintf (stderr, "\n*no begin line");
951 return (FALSE);
952 }
953
954 while (strcmp (name, "begin") != 0) /*skip until "begin" found processing options */
955 {
956 if (strncmp (name, "*", 1) == 0) /* skip any line beginning with * */
957 {
958 c = name[0];
959 while (c != EOF && c != '\n')
960 c = fgetc (lrs_ifp);
961 }
962
963 else if (strcmp (name, "H-representation") == 0)
964 Q->hull = FALSE;
965 else if ((strcmp (name, "hull") == 0) || (strcmp (name, "V-representation") == 0))
966 {
967 Q->hull = TRUE;
968 Q->polytope = TRUE; /* will be updated as input read */
969 }
970 else if (strcmp (name, "digits") == 0)
971 {
972 if (fscanf (lrs_ifp, "%ld", &dec_digits) == EOF)
973 {
974 fprintf (stderr, "\n*no begin line");
975 return (FALSE);
976 }
977 if (!lrs_set_digits(dec_digits))
978 return (FALSE);
979 }
980 else if (strcmp (name, "linearity") == 0)
981 {
982 if (!readlinearity (Q))
983 return FALSE;
984 }
985 else if (strcmp (name, "nonnegative") == 0)
986 {
987 if(Q->nash)
988 fprintf (stderr, "\nNash incompatibile with nonnegative option - skipped");
989 else
990 Q->nonnegative = TRUE;
991 }
992 else if (firstline)
993 {
994 // printf("\nov=%ld mess=%ld",overflow,Q->messages);
995 if(overflow != 2)
996 lrs_warning(Q,"warning",name);
997 firstline = FALSE;
998 }
999
1000 if (fscanf (lrs_ifp, "%s", name) == EOF)
1001 {
1002 fprintf (stderr, "\n*no begin line\n");
1003 return (FALSE);
1004 }
1005
1006 } /* end of while */
1007
1008
1009 if (fscanf (lrs_ifp, "%ld %ld %s", &Q->m, &Q->n, name) == EOF)
1010 {
1011 fprintf (stderr, "\n*no data in file\n");
1012 return (FALSE);
1013 }
1014 if (strcmp (name, "integer") != 0 && strcmp (name, "rational") != 0)
1015 {
1016 fprintf (stderr,"\n*data type must be integer of rational\n");
1017 return (FALSE);
1018 }
1019
1020
1021 if (Q->m == 0)
1022 {
1023 fprintf (stderr, "\n*no input given\n"); /* program dies ungracefully */
1024 return (FALSE);
1025 }
1026
1027
1028 /* inputd may be reduced in preprocessing of linearities and redund cols */
1029
1030 return TRUE;
1031 } /* end of lrs_read_dat */
1032
1033 /****************************/
1034 /* set up lrs_dic structure */
1035 /****************************/
1036 long
lrs_read_dic(lrs_dic * P,lrs_dat * Q)1037 lrs_read_dic (lrs_dic * P, lrs_dat * Q)
1038 /* read constraint matrix and set up problem and dictionary */
1039
1040 {
1041 lrs_mp Temp,Tempn,Tempd, mpone, mptwo;
1042
1043 long i, j, m, d;;
1044 char name[100];
1045 char mess[100];
1046 char *ss;
1047 int c; /* fgetc actually returns an int. db */
1048 long redundstart, redundend;
1049 long dualperturb=FALSE; /* dualperturb=TRUE: objective function perturbed */
1050 /* assign local variables to structures */
1051
1052 lrs_mp_matrix A;
1053 lrs_mp_vector Gcd, Lcm;
1054
1055 lrs_mp_vector oD=Q->output; /* Denom for obj fun, temp use of output */
1056 long hull = Q->hull;
1057 long messages=Q->messages;
1058
1059 lrs_alloc_mp(Temp); lrs_alloc_mp(mpone);
1060 lrs_alloc_mp(Tempn); lrs_alloc_mp(Tempd); lrs_alloc_mp(mptwo);
1061 A = P->A;
1062 m = Q->m;
1063 d = Q->inputd;
1064 Gcd = Q->Gcd;
1065 Lcm = Q->Lcm;
1066
1067 itomp (ONE, mpone);
1068 itomp(2L,mptwo);
1069 itomp (ONE, A[0][0]);
1070 itomp (ONE, Lcm[0]);
1071 itomp (ONE, Gcd[0]);
1072
1073 for (i = 1; i <= m; i++) /* read in input matrix row by row */
1074 {
1075 itomp (ONE, Lcm[i]); /* Lcm of denominators */
1076 itomp (ZERO, Gcd[i]); /* Gcd of numerators */
1077 for (j = hull; j <= d; j++) /* hull data copied to cols 1..d */
1078 {
1079 if (readrat (A[i][j], A[0][j]))
1080 lcm (Lcm[i], A[0][j]); /* update lcm of denominators */
1081 copy (Temp, A[i][j]);
1082 gcd (Gcd[i], Temp); /* update gcd of numerators */
1083 }
1084
1085 if (hull)
1086 {
1087 itomp (ZERO, A[i][0]); /*for hull, we have to append an extra column of zeroes */
1088 if (!one (A[i][1]) || !one (A[0][1])) /* all rows must have a one in column one */
1089 Q->polytope = FALSE;
1090 }
1091 if (!zero (A[i][hull])) /* for H-rep, are zero in column 0 */
1092 Q->homogeneous = FALSE; /* for V-rep, all zero in column 1 */
1093
1094 storesign (Gcd[i], POS);
1095 storesign (Lcm[i], POS);
1096 if (mp_greater (Gcd[i], mpone) || mp_greater (Lcm[i], mpone))
1097 for (j = 0; j <= d; j++)
1098 {
1099 exactdivint (A[i][j], Gcd[i], Temp); /*reduce numerators by Gcd */
1100 mulint (Lcm[i], Temp, Temp); /*remove denominators */
1101 exactdivint (Temp, A[0][j], A[i][j]); /*reduce by former denominator */
1102 }
1103
1104 } /* end of for i= */
1105
1106 /* 2010.4.26 patch */
1107 if(Q->nonnegative) /* set up Gcd and Lcm for nonexistent nongative inequalities */
1108 for (i=m+1;i<=m+d;i++)
1109 { itomp (ONE, Lcm[i]);
1110 itomp (ONE, Gcd[i]);
1111 }
1112
1113 if (Q->homogeneous && Q->verbose && overflow != 2)
1114 {
1115 lrs_warning(Q,"warning","*Input is homogeneous, column 1 not treated as redundant");
1116 }
1117
1118
1119 /* read in flags */
1120 while (fscanf (lrs_ifp, "%s", name) != EOF)
1121 {
1122 if (strncmp (name, "*", 1) == 0) /* skip any line beginning with * */
1123 {
1124 c = name[0];
1125 while (c != EOF && c != '\n')
1126 c = fgetc (lrs_ifp);
1127 }
1128
1129
1130 if (strcmp (name, "checkpoint") == 0)
1131 {
1132 long seconds;
1133
1134 if(fscanf (lrs_ifp, "%ld", &seconds) == EOF)
1135 {
1136 fprintf (lrs_ofp, "\nInvalid checkpoint option");
1137 return (FALSE);
1138 }
1139
1140 #ifndef SIGNALS
1141 if (seconds > 0)
1142 {
1143 lrs_checkpoint_seconds = seconds;
1144 errcheck ("signal", signal (SIGALRM, timecheck));
1145 alarm (lrs_checkpoint_seconds);
1146 }
1147 #endif
1148 }
1149
1150 if (strcmp (name, "debug") == 0)
1151 {
1152 Q->etrace =0;
1153 if(fscanf (lrs_ifp, "%ld %ld", &Q->strace, &Q->etrace)==EOF)
1154 Q->strace =0;
1155 if(!Q->mplrs)
1156 fprintf (lrs_ofp, "\n*%s from B#%ld to B#%ld", name, Q->strace, Q->etrace);
1157 Q->debug=TRUE;
1158 if (Q->strace <= 1)
1159 Q->debug = TRUE;
1160 }
1161 if (strcmp (name, "startingcobasis") == 0)
1162 {
1163 if(Q->nonnegative)
1164 lrs_warning(Q,"warning", "*startingcobasis incompatible with nonnegative option:skipped");
1165 else
1166 {
1167 if(!Q->mplrs && messages && overflow != 2)
1168 fprintf (lrs_ofp, "\n*startingcobasis");
1169 Q->givenstart = TRUE;
1170 if (!readfacets (Q, Q->inequality))
1171 return FALSE;
1172 }
1173 }
1174
1175
1176 if (strcmp (name, "restart")==0 )
1177 {
1178 if(Q->mplrs)
1179 {
1180 fprintf (lrs_ofp, "\n\n*** %s is an lrs option only\n", name);
1181 return(FALSE);
1182 }
1183 Q->restart = TRUE;
1184 Q->giveoutput = FALSE; /* first time only */
1185 if(Q->voronoi)
1186 {
1187 if(fscanf (lrs_ifp, "%ld %ld %ld %ld", &Q->count[1], &Q->count[0], &Q->count[2], &P->depth)==EOF)
1188 return FALSE;
1189 if(!Q->mplrs && messages)
1190 fprintf (lrs_ofp, "\n*%s V#%ld R#%ld B#%ld h=%ld data points", name, Q->count[1], Q->count[0], Q->count[2], P->depth);
1191 }
1192 else if(hull)
1193 {
1194 if( fscanf (lrs_ifp, "%ld %ld %ld", &Q->count[0], &Q->count[2], &P->depth)==EOF)
1195 return(FALSE);
1196 if(!Q->mplrs && messages && overflow != 2)
1197 fprintf (lrs_ofp, "\n*%s F#%ld B#%ld h=%ld vertices/rays", name, Q->count[0], Q->count[2], P->depth);
1198 }
1199 else
1200 {
1201 if(fscanf (lrs_ifp, "%ld %ld %ld %ld", &Q->count[1], &Q->count[0], &Q->count[2], &P->depth)==EOF)
1202 return FALSE;
1203 if(!Q->mplrs && messages && overflow != 2)
1204 fprintf (lrs_ofp, "\n*%s V#%ld R#%ld B#%ld h=%ld facets", name, Q->count[1], Q->count[0], Q->count[2], P->depth);
1205 }
1206 /* store starting counts to calculate totals of plrs/mplrs subjob */
1207 for (i = 0; i<5; i++)
1208 Q->startcount[i] = Q->count[i];
1209 if (!readfacets (Q, Q->facet))
1210 return FALSE;
1211 } /* end of restart */
1212
1213 /* The next flags request a LP solution only */
1214 if(Q->mplrs)
1215 {
1216 if (strncmp (name,"lponly",6)== 0)
1217 {
1218 fprintf (lrs_ofp, "\n\n*** %s is an lrs option only\n", name);
1219 return(FALSE);
1220 }
1221 }
1222 else
1223 {
1224 if (strcmp (name, "lponly") == 0 || strcmp (name, "lponly_d") == 0 )
1225
1226 {
1227 if (Q->hull)
1228 fprintf (lrs_ofp, "\n*lponly option not valid for V-representation-skipped");
1229 else
1230 {
1231 Q->lponly = 1; /*Dantzig's rule is default */
1232 fprintf (lrs_ofp, "\n*Dantzig's rule");
1233 }
1234 }
1235
1236 if (strcmp (name, "lponly_r") == 0)
1237
1238 {
1239 if (Q->hull)
1240 fprintf (lrs_ofp, "\n*lponly option not valid for V-representation-skipped");
1241 else
1242 {
1243 Q->lponly = 2; /*random edge rule */
1244 fprintf (lrs_ofp, "\n*random edge rule");
1245 }
1246 }
1247
1248 if (strcmp (name, "lponly_rd") == 0)
1249
1250 {
1251 if (Q->hull)
1252 fprintf (lrs_ofp, "\n*lponly option not valid for V-representation-skipped");
1253 else
1254 {
1255 Q->lponly = 3; /*random edge/Dantzig hybrid */
1256 fprintf (lrs_ofp, "\n*random edge/Dantzig hybrid");
1257 }
1258 }
1259
1260 if (strcmp (name, "lponly_b") == 0)
1261
1262 {
1263 if (Q->hull)
1264 fprintf (lrs_ofp, "\n*lponly option not valid for V-representation-skipped");
1265 else
1266 {
1267 Q->lponly = 4; /*Bland' rule */
1268 fprintf (lrs_ofp, "\n*Bland's rule");
1269 }
1270 }
1271 } /* else !Q->mplrs */
1272 /* The LP will be solved after initialization to get starting vertex */
1273
1274
1275
1276 if (strcmp (name, "maximize") == 0 || strcmp (name, "minimize") == 0)
1277 {
1278 if (Q->hull )
1279 lrs_warning(Q,"warning","*minimize/maximize options not valid for V-representation-skipped");
1280 else
1281 {
1282 if (strcmp (name, "maximize") == 0)
1283 Q->maximize = TRUE;
1284 else
1285 Q->minimize = TRUE;
1286 if(overflow != 2)
1287 lrs_warning(Q,"warning",name);
1288
1289 if(dualperturb && overflow != 2) /* apply a perturbation to objective function */
1290 {
1291 lrs_warning(Q,"warning","*Objective function perturbed");
1292 copy(Temp,mptwo);
1293 }
1294
1295 for (j = 0; j <= d; j++)
1296 {
1297 if (readrat (A[0][j], oD[j]) || dualperturb )
1298 {
1299 if(dualperturb && j > 0 && j < d )
1300 {
1301 linrat(A[0][j], oD[j],ONE,mpone,Temp,ONE,Tempn,Tempd);
1302 copy(A[0][j],Tempn);
1303 copy(oD[j],Tempd);
1304 mulint(mptwo,Temp,Temp);
1305 }
1306 reduce (A[0][j], oD[j]);
1307 lcm (Q->Lcm[0], oD[j]); /* update lcm of denominators */
1308 }
1309 }
1310
1311
1312 for (j = 0; j <= d; j++)
1313 if (!Q->maximize)
1314 changesign (A[0][j]);
1315 storesign (Q->Lcm[0], POS);
1316 if (mp_greater (Q->Lcm[0], mpone))
1317 for (j = 0; j <= d; j++)
1318 {
1319 mulint (Q->Lcm[0], A[0][j], A[0][j]); /*remove denominators */
1320 copy (Temp, A[0][j]);
1321 exactdivint (Temp, oD[j], A[0][j]);
1322 }
1323 if(messages && overflow !=2 )
1324 lrs_printoutput(Q,A[0]);
1325
1326
1327 if (Q->debug)
1328 printA (P, Q);
1329 } /* else */
1330 } /* end of LP setup */
1331 if (strcmp (name, "volume") == 0)
1332 {
1333 if(overflow != 2)
1334 lrs_warning(Q,"warning", "*volume");
1335 Q->getvolume = TRUE;
1336 }
1337 if (strcmp (name, "geometric") == 0)
1338 {
1339 if (hull && !Q->voronoi)
1340 lrs_warning(Q,"warning", "*geometric - option for H-representation or voronoi only, skipped");
1341 else
1342 {
1343 lrs_warning(Q,"warning", "*geometric");
1344 Q->geometric = TRUE;
1345 }
1346 }
1347 if (strcmp (name, "allbases") == 0)
1348 {
1349 if(overflow != 2)
1350 lrs_warning (Q,"warning", "*allbases");
1351 Q->allbases = TRUE;
1352 }
1353
1354 if (strcmp (name, "countonly") == 0)
1355 {
1356 if(overflow != 2)
1357 lrs_warning (Q,"warning", "*countonly");
1358 Q->countonly = TRUE;
1359 }
1360
1361 if (strcmp (name, "triangulation") == 0)
1362 {
1363 if (hull)
1364 {
1365 if(overflow != 2)
1366 lrs_warning (Q,"warning","*triangulation");
1367 Q->triangulation = TRUE;
1368 Q->getvolume = TRUE;
1369 }
1370 else
1371 printf ("\n*triangulation only valid for V-representations: skipped");
1372 }
1373 if (strcmp (name, "dualperturb") == 0)
1374 {
1375 dualperturb = TRUE;
1376 }
1377
1378 if (strcmp (name, "incidence") == 0)
1379 {
1380 if(overflow != 2)
1381 lrs_warning (Q,"warning", "*incidence");
1382 Q->incidence = TRUE;
1383 }
1384
1385 if (strcmp (name, "#incidence") == 0) /* number of incident inequalities only */
1386 {
1387 Q->printcobasis = TRUE;
1388 }
1389
1390 if (strcmp (name, "printcobasis") == 0)
1391 {
1392 if(fscanf (lrs_ifp, "%ld", &Q->frequency)==EOF)
1393 /*2010.7.7 set default to zero = print only when outputting vertex/ray/facet */
1394 Q->frequency=0;
1395 if(overflow != 2)
1396 {
1397 if (Q->frequency > 0)
1398 sprintf(mess,"*%s %ld",name, Q->frequency);
1399 else
1400 sprintf(mess,"%s",name);
1401 lrs_warning(Q,"warning",mess);
1402 }
1403 Q->printcobasis = TRUE;
1404 }
1405
1406 if (strcmp (name, "integervertices") == 0) /* when restarting reinitialize */
1407 {
1408 if(fscanf (lrs_ifp, "%ld", &Q->count[4])==EOF)
1409 Q->count[4]=0;
1410 }
1411
1412 if (strcmp (name, "printslack") == 0)
1413 {
1414 Q->printslack = TRUE;
1415 }
1416
1417 if (strcmp (name, "cache") == 0)
1418 {
1419 if(fscanf (lrs_ifp, "%ld", &dict_limit)==EOF)
1420 dict_limit=1;
1421 if (dict_limit < 1)
1422 dict_limit = 1;
1423 if(overflow != 2)
1424 {
1425 sprintf(mess,"*%s %ld",name, dict_limit);
1426 lrs_warning(Q,"warning",mess);
1427 }
1428
1429 }
1430 if (strcmp (name, "linearity") == 0)
1431 {
1432 if (!readlinearity (Q))
1433 return FALSE;
1434 }
1435
1436 if (strcmp (name, "maxdepth") == 0)
1437 {
1438 Q->maxdepth=MAXD;
1439 if(fscanf (lrs_ifp, "%lld", &Q->maxdepth)==EOF)
1440 Q->maxdepth=MAXD;
1441 if(overflow != 2)
1442 {
1443 if(Q->mplrs) /* taken from control line */
1444 lrs_warning(Q,"warning","*maxdepth option skipped - supplied on control line in mplrs");
1445 else
1446 fprintf (lrs_ofp, "\n*%s %lld", name, Q->maxdepth);
1447 }
1448 }
1449
1450 if (strcmp (name, "mindepth") == 0)
1451 {
1452 Q->mindepth=0;
1453 if(fscanf (lrs_ifp, "%lld", &Q->mindepth)==EOF)
1454 Q->mindepth=0;
1455 if(overflow != 2)
1456 {
1457 if(Q->mplrs) /* taken from control line */
1458 lrs_warning(Q,"warning","*mindepth option skipped in mplrs");
1459 else
1460 fprintf (lrs_ofp, "\n*%s %lld", name, Q->mindepth);
1461 }
1462 }
1463
1464 if (strcmp (name, "maxoutput") == 0)
1465 {
1466 if(fscanf (lrs_ifp, "%ld", &Q->maxoutput)==EOF)
1467 Q->maxoutput = 100;
1468 if(overflow != 2)
1469 {
1470 if(Q->mplrs) /* taken from control line */
1471 lrs_warning(Q,"warning","*maxoutput option skipped in mplrs");
1472 else
1473 fprintf (lrs_ofp, "\n*%s %ld", name, Q->maxoutput);
1474 }
1475 }
1476
1477 if (strcmp (name, "maxcobases") == 0)
1478 {
1479 if(fscanf (lrs_ifp, "%ld", &Q->maxcobases)==EOF)
1480 Q->maxcobases = 1000;
1481 if(overflow != 2)
1482 {
1483 if(Q->mplrs) /* taken from control line */
1484 lrs_warning(Q,"warning","*maxcobases option skipped - supplied on control line in mplrs");
1485 else
1486 fprintf (lrs_ofp, "\n*%s %ld", name, Q->maxcobases);
1487 }
1488
1489 }
1490
1491
1492 /*2019.8.24 bounds for redund */
1493 if (strcmp (name, "redund") == 0)
1494 {
1495 strcpy(Q->fname,"redund");
1496 for (i = 1; i <= Q->m; i++) /*reset any previous redund option except =2 values */
1497 if (Q->redineq[i] != 2)
1498 Q->redineq[i]=0;
1499 Q->redineq[0]=1;
1500
1501 if( fscanf (lrs_ifp, "%ld %ld", &redundstart, &redundend)==EOF)
1502 {
1503 redundstart=1;
1504 redundend=m;
1505 }
1506 if (redundstart <1 || redundstart > redundend )
1507 redundstart=1;
1508 if (redundend < 1 || redundend > Q->m )
1509 redundend=Q->m;
1510 for (i=redundstart;i<=redundend;i++)
1511 Q->redineq[i]=1;
1512 if(overflow != 2 )
1513 {
1514 sprintf (mess, "%s %ld %ld", name, redundstart, redundend);
1515 lrs_warning(Q,"warning",mess);
1516 }
1517
1518 }
1519
1520 if (strcmp (name, "redund_list") == 0)
1521 {
1522 strcpy(Q->fname,"redund");
1523 readredund(Q);
1524 }
1525
1526 if (strcmp (name, "noredundcheck") == 0)
1527 {
1528 if(Q->mplrs)
1529 {
1530 if(messages && overflow != 2)
1531 lrs_warning(Q,"warning","*noredundcheck");
1532 Q->noredundcheck = TRUE;
1533 }
1534 }
1535
1536 if (strcmp (name, "truncate") == 0)
1537 {
1538 if (hull)
1539 {
1540 if(overflow != 2)
1541 lrs_warning(Q,"warning","*truncate - option for H-representation only, skipped");
1542 }
1543 else
1544 {
1545 Q->truncate = TRUE;
1546 lrs_warning(Q,"warning","*truncate");
1547 }
1548 }
1549 if (strcmp (name, "extract") == 0)
1550 {
1551 if(Q->mplrs)
1552 {
1553 fprintf (lrs_ofp, "\n\n*** %s is an lrs option only\n", name);
1554 return(FALSE);
1555 }
1556 Q->extract = TRUE;
1557 readremain(Q);
1558 }
1559
1560 if (strcmp (name, "verbose") == 0)
1561 Q->verbose = TRUE;
1562
1563 if (strcmp (name, "bound") == 0)
1564 {
1565 readrat(Q->boundn,Q->boundd);
1566 Q->bound = TRUE;
1567 }
1568
1569 if (strcmp (name, "nonnegative") == 0)
1570 if(overflow != 2)
1571 lrs_warning(Q,"warning","*nonnegative - option must come before begin line - skipped");
1572
1573 if (strcmp (name, "seed") == 0)
1574 {
1575 if(fscanf (lrs_ifp, "%ld", &Q->seed)==EOF)
1576 Q->seed = 3142;
1577 sprintf(mess,"*seed=%ld",Q->seed);
1578 if(overflow != 2)
1579 lrs_warning(Q,"warning",mess);
1580 srandom(Q->seed);
1581 }
1582
1583 if (strcmp (name, "estimates") == 0)
1584 {
1585 if(Q->mplrs)
1586 {
1587 fprintf (lrs_ofp, "\n\n*** %s is an lrs option only\n", name);
1588 return(FALSE);
1589 }
1590 if(fscanf (lrs_ifp, "%ld", &Q->runs)==EOF)
1591 Q->runs=1;
1592 if(messages && overflow != 2)
1593 fprintf (lrs_ofp, "\n*%ld %s", Q->runs, name);
1594 }
1595
1596 // 2015.2.9 Estimates will continue until estimate is less than subtree size
1597 if (strcmp (name, "subtreesize") == 0)
1598 {
1599 if(fscanf (lrs_ifp, "%lld", &Q->subtreesize)==EOF)
1600 Q->subtreesize=MAXD;
1601 if(messages && overflow != 2)
1602 fprintf (lrs_ofp, "\n*%s %lld", name, Q->subtreesize);
1603 if(Q->mplrs)
1604 {
1605 fprintf (lrs_ofp, "\n\n*** %s is an lrs option only\n", name);
1606 return(FALSE);
1607 }
1608 }
1609
1610
1611
1612 if ((strcmp (name, "voronoi") == 0) || (strcmp (name, "Voronoi") == 0))
1613 {
1614 if (!hull)
1615 {
1616 if(overflow != 2)
1617 lrs_warning(Q,"warning","*voronoi requires V-representation - option skipped");
1618 }
1619 else
1620 {
1621 lrs_warning(Q,"warning","*voronoi");
1622 Q->voronoi = TRUE;
1623 Q->polytope = FALSE;
1624 }
1625 }
1626
1627 } /* end of while for reading flags */
1628
1629 if(Q->bound)
1630 {
1631 if(Q->maximize)
1632 ss=cprat("*Lower bound on objective function:",Q->boundn,Q->boundd);
1633 else
1634 ss=cprat("*Upper bound on objective function:",Q->boundn,Q->boundd);
1635 lrs_warning(Q,"warning",ss);
1636 free(ss);
1637 }
1638
1639 /* Certain options are incompatible, this is fixed here */
1640
1641 if (Q->restart && Q->maxcobases > 0) //2015.4.3 adjust for restart
1642 Q->maxcobases = Q->maxcobases + Q->count[2];
1643
1644 if (Q->incidence)
1645 Q->printcobasis = TRUE;
1646
1647 if (Q->debug)
1648 {
1649 printA (P, Q);
1650 fprintf (lrs_ofp, "\nexiting lrs_read_dic");
1651 }
1652 fflush(lrs_ofp); fflush(stdout);
1653
1654 /*removing tmpfiles */
1655
1656 fclose(lrs_ifp);
1657 lrs_ifp=NULL;
1658
1659 if ( overflow > 0 ) /* we made a temporary file for overflow or stdin */
1660 if(remove(infilename) != 0)
1661 lrs_warning(Q,"warning","*Could not delete temporary file");
1662
1663 if (*tmpfilename != '\0' ) /* we made a temporary file for stdin */
1664 if(remove(tmpfilename) != 0)
1665 lrs_warning(Q,"warning", "*Could not delete temporary file");
1666 *tmpfilename = '\0';
1667
1668 lrs_clear_mp(Temp); lrs_clear_mp(mpone);
1669 lrs_clear_mp(Tempn); lrs_clear_mp(Tempd); lrs_clear_mp(mptwo);
1670 return TRUE;
1671
1672 } /* lrs_read_dic */
1673
1674
1675 /* In lrs_getfirstbasis and lrs_getnextbasis we use D instead of P */
1676 /* since the dictionary P may change, ie. &P in calling routine */
1677
1678 #define D (*D_p)
1679
1680 long
lrs_getfirstbasis(lrs_dic ** D_p,lrs_dat * Q,lrs_mp_matrix * Lin,long no_output)1681 lrs_getfirstbasis (lrs_dic ** D_p, lrs_dat * Q, lrs_mp_matrix * Lin, long no_output)
1682 /* gets first basis, FALSE if none */
1683 /* P may get changed if lin. space Lin found */
1684 /* no_output is TRUE supresses output headers */
1685 /* 2017.12.22 could use no_output=2 to get early exit for criss-cross method */
1686 {
1687 lrs_mp scale, Temp;
1688
1689 long i, j, k;
1690
1691 /* assign local variables to structures */
1692
1693 lrs_mp_matrix A;
1694 long *B, *C, *Col;
1695 long *inequality;
1696 long *linearity;
1697 long hull = Q->hull;
1698 long m, d, lastdv, nlinearity, nredundcol;
1699
1700 lrs_alloc_mp(Temp); lrs_alloc_mp(scale);
1701
1702 if (Q->lponly)
1703 no_output = TRUE;
1704 m = D->m;
1705 d = D->d;
1706
1707
1708
1709 lastdv = Q->lastdv;
1710
1711 nredundcol = 0L; /* will be set after getabasis */
1712 nlinearity = Q->nlinearity; /* may be reset if new linearity read or in getabasis*/
1713 linearity = Q->linearity;
1714
1715
1716
1717 A = D->A;
1718 B = D->B;
1719 C = D->C;
1720 Col = D->Col;
1721 inequality = Q->inequality;
1722
1723
1724 /**2020.2.5 no linearities to remove so just select columns and quit */
1725
1726 if (Q->extract && Q->nlinearity == 0)
1727 return extractcols(D,Q);
1728
1729 if(Q->verbose && overflow != 2)
1730 {
1731 if (Q->nlinearity > 0 && Q->nonnegative)
1732 {
1733 fprintf (lrs_ofp, "\n*linearity and nonnegative options incompatible");
1734 fprintf (lrs_ofp, " - all linearities are skipped");
1735 fprintf (lrs_ofp, "\n*add nonnegative constraints explicitly and ");
1736 fprintf (lrs_ofp, " remove nonnegative option");
1737 }
1738
1739 if (Q->nlinearity && Q->voronoi)
1740 fprintf (lrs_ofp, "\n*linearity and Voronoi options set - results unpredictable");
1741
1742 if (Q->lponly && !Q->maximize && !Q->minimize)
1743 fprintf (lrs_ofp, "\n*LP has no objective function given - assuming all zero");
1744
1745 }
1746
1747 if (Q->runs > 0) /* arrays for estimator */
1748 {
1749 Q->isave = (long *) CALLOC ((unsigned) (m * d), sizeof (long));
1750 Q->jsave = (long *) CALLOC ((unsigned) (m * d), sizeof (long));
1751 }
1752 /* default is to look for starting cobasis using linearies first, then */
1753 /* filling in from last rows of input as necessary */
1754 /* linearity array is assumed sorted here */
1755 /* note if restart/given start inequality indices already in place */
1756 /* from nlinearity..d-1 */
1757 for (i = 0; i < nlinearity; i++) /* put linearities first in the order */
1758 inequality[i] = linearity[i];
1759
1760 k = 0; /* index for linearity array */
1761
1762 if (Q->givenstart)
1763 k = d;
1764 else
1765 k = nlinearity;
1766 for (i = m; i >= 1; i--)
1767 {
1768 j = 0;
1769 while (j < k && inequality[j] != i)
1770 j++; /* see if i is in inequality */
1771 if (j == k)
1772 inequality[k++] = i;
1773 }
1774 if (Q->debug)
1775 {
1776 fprintf (lrs_ofp, "\n*Starting cobasis uses input row order");
1777 for (i = 0; i < m; i++)
1778 fprintf (lrs_ofp, " %ld", inequality[i]);
1779
1780 }
1781 /* for voronoi convert to h-description using the transform */
1782 /* a_0 .. a_d-1 -> (a_0^2 + ... a_d-1 ^2)-2a_0x_0-...-2a_d-1x_d-1 + x_d >= 0 */
1783 /* note constant term is stored in column d, and column d-1 is all ones */
1784 /* the other coefficients are multiplied by -2 and shifted one to the right */
1785 if (Q->debug)
1786 printA (D, Q);
1787 if (Q->voronoi)
1788 {
1789 Q->hull = FALSE;
1790 hull = FALSE;
1791 for (i = 1; i <= m; i++)
1792 {
1793 if (zero (A[i][1]))
1794 {
1795 printf("\nWith voronoi option column one must be all one\n");
1796 return (FALSE);
1797 }
1798 copy (scale, A[i][1]); /*adjust for scaling to integers of rationals */
1799 itomp (ZERO, A[i][0]);
1800 for (j = 2; j <= d; j++) /* transform each input row */
1801 {
1802 copy (Temp, A[i][j]);
1803 mulint (A[i][j], Temp, Temp);
1804 linint (A[i][0], ONE, Temp, ONE);
1805 linint (A[i][j - 1], ZERO, A[i][j], -TWO);
1806 mulint (scale, A[i][j - 1], A[i][j - 1]);
1807 } /* end of for (j=1;..) */
1808 copy (A[i][d], scale);
1809 mulint (scale, A[i][d], A[i][d]);
1810 }/* end of for (i=1;..) */
1811 if (Q->debug)
1812 printA (D, Q);
1813 } /* end of if(voronoi) */
1814 if (!Q->maximize && !Q->minimize)
1815 for (j = 0; j <= d; j++)
1816 itomp (ZERO, A[0][j]);
1817
1818 /* Now we pivot to standard form, and then find a primal feasible basis */
1819 /* Note these steps MUST be done, even if restarting, in order to get */
1820 /* the same index/inequality correspondance we had for the original prob. */
1821 /* The inequality array is used to give the insertion order */
1822 /* and is defaulted to the last d rows when givenstart=FALSE */
1823
1824 if(Q->nonnegative)
1825 {
1826 /* no need for initial pivots here, labelling already done */
1827 Q->lastdv = d;
1828 Q->nredundcol = 0;
1829 }
1830 else
1831 {
1832 if (!getabasis (D, Q, inequality))
1833 return FALSE;
1834 /* bug fix 2009.12.2 */
1835 nlinearity=Q->nlinearity; /*may have been reset if some lins are redundant*/
1836 }
1837
1838 /* 2020.2.2 */
1839 /* extract option asked to remove all linearities and output the reduced A matrix */
1840 /* should be followed by redund to get minimum representation */
1841
1842 if (Q->extract)
1843 return linextractcols(D,Q);
1844
1845 if(Q->debug)
1846 {
1847 fprintf(lrs_ofp,"\nafter getabasis");
1848 printA(D, Q);
1849 }
1850 nredundcol = Q->nredundcol;
1851 lastdv = Q->lastdv;
1852 d = D->d;
1853
1854
1855
1856 /********************************************************************/
1857 /* now we start printing the output file unless no output requested */
1858 /********************************************************************/
1859
1860
1861 if (Q->count[2]==1 && (no_output==0 || Q->debug)) /* don't reprint after newstart */
1862 {
1863 int len=0;
1864 char *header;
1865
1866 header=(char *)malloc((100+20*Q->n)*sizeof(char));
1867
1868 if (Q->voronoi)
1869 len=sprintf (header, "*Voronoi Diagram: Voronoi vertices and rays are output");
1870 else
1871 {
1872 if (hull)
1873 len=sprintf (header, "H-representation");
1874 else
1875 len=sprintf (header, "V-representation");
1876 }
1877
1878 /* Print linearity space */
1879 /* Don't print linearity if first column zero in hull computation */
1880
1881 if (hull && Q->homogeneous)
1882 k = 1; /* 0 normally, 1 for homogeneous case */
1883 else
1884 k = 0;
1885
1886 if (nredundcol > k)
1887 {
1888 len=len+sprintf (header+len, "\nlinearity %ld ", nredundcol - k); /*adjust nredundcol for homog. */
1889 for (i = 1; i <= nredundcol - k; i++)
1890 len=len+sprintf (header+len, " %ld", i);
1891 } /* end print of linearity space */
1892
1893 len=len+sprintf (header+len, "\nbegin");
1894 len=len+sprintf (header+len, "\n***** %ld rational", Q->n);
1895
1896 if(Q->mplrs)
1897 lrs_post_output("header",header);
1898 else
1899 if(Q->messages)
1900 fprintf(lrs_ofp,"\n%s",header);
1901
1902 free(header);
1903 }
1904
1905 /* end of if !no_output ....... */
1906
1907
1908 /* Reset up the inequality array to remember which index is which input inequality */
1909 /* inequality[B[i]-lastdv] is row number of the inequality with index B[i] */
1910 /* inequality[C[i]-lastdv] is row number of the inequality with index C[i] */
1911
1912 for (i = 1; i <= m; i++)
1913 inequality[i] = i;
1914 if (nlinearity > 0) /* some cobasic indices will be removed */
1915 {
1916 for (i = 0; i < nlinearity; i++) /* remove input linearity indices */
1917 inequality[linearity[i]] = 0;
1918 k = 1; /* counter for linearities */
1919 for (i = 1; i <= m - nlinearity; i++)
1920 {
1921 while (k <= m && inequality[k] == 0)
1922 k++; /* skip zeroes in corr. to linearity */
1923 inequality[i] = inequality[k++];
1924 }
1925 }
1926 /* end if linearity */
1927
1928 if (Q->debug)
1929 {
1930 fprintf (lrs_ofp, "\ninequality array initialization:");
1931 for (i = 1; i <= m - nlinearity; i++)
1932 fprintf (lrs_ofp, " %ld", inequality[i]);
1933
1934 }
1935
1936
1937
1938
1939 if (nredundcol > 0)
1940 {
1941 const unsigned int Qn = Q->n;
1942 *Lin = lrs_alloc_mp_matrix (nredundcol, Qn);
1943
1944 for (i = 0; i < nredundcol; i++)
1945 {
1946
1947
1948 if (!(Q->homogeneous && Q->hull && i == 0)) /* skip redund col 1 for homog. hull */
1949 {
1950
1951 lrs_getray (D, Q, Col[0], D->C[0] + i - hull, (*Lin)[i]); /* adjust index for deletions */
1952 }
1953
1954 if (!removecobasicindex (D, Q, 0L))
1955 {
1956 lrs_clear_mp_matrix (*Lin, nredundcol, Qn);
1957 return FALSE;
1958 }
1959 }
1960
1961
1962 } /* end if nredundcol > 0 */
1963
1964
1965
1966 if (Q->lponly || Q->nash ){
1967 if (Q->verbose )
1968 {
1969 fprintf (lrs_ofp, "\nNumber of pivots for starting dictionary: %ld",Q->count[3]);
1970 if(Q->lponly && Q->debug)
1971 printA (D, Q);
1972 }
1973 }
1974
1975 /*2017.12.22 If you want to do criss-cross now is the time ! */
1976
1977
1978 /* Do dual pivots to get primal feasibility */
1979 if (!primalfeasible (D, Q))
1980 {
1981 if(!Q->mplrs)
1982 fprintf (lrs_ofp, "\nend");
1983 lrs_warning(Q,"finalwarn", "\nNo feasible solution\n");
1984 if (Q->nash && Q->verbose )
1985 {
1986 fprintf (lrs_ofp, "\nNumber of pivots for feasible solution: %ld",Q->count[3]);
1987 fprintf (lrs_ofp, " - No feasible solution");
1988 }
1989 return FALSE;
1990 }
1991
1992 if (Q->lponly || Q->nash )
1993 if (Q->verbose)
1994 {
1995 fprintf (lrs_ofp, "\nNumber of pivots for feasible solution: %ld",Q->count[3]);
1996 if(Q->lponly && Q->debug)
1997 printA (D, Q);
1998 }
1999
2000
2001
2002 /* Now solve LP if objective function was given */
2003 if (Q->maximize || Q->minimize)
2004 {
2005 Q->unbounded = !lrs_solvelp (D, Q, Q->maximize);
2006 if (Q->lponly)
2007 {
2008
2009 if (Q->verbose)
2010 {
2011 fprintf (lrs_ofp, "\nNumber of pivots for optimum solution: %ld",Q->count[3]);
2012 if(Q->debug)
2013 printA (D, Q);
2014 }
2015 lrs_clear_mp(Temp); lrs_clear_mp(scale);
2016 return TRUE;
2017 }
2018
2019 else /* check to see if objective is dual degenerate */
2020 {
2021 j = 1;
2022 while (j <= d && !zero (A[0][j]))
2023 j++;
2024 if (j <= d)
2025 Q->dualdeg = TRUE;
2026 }
2027 }
2028 else
2029 /* re-initialize cost row to -det */
2030 {
2031 for (j = 1; j <= d; j++)
2032 {
2033 copy (A[0][j], D->det);
2034 storesign (A[0][j], NEG);
2035 }
2036
2037 itomp (ZERO, A[0][0]); /* zero optimum objective value */
2038 }
2039
2040
2041 /* reindex basis to 0..m if necessary */
2042 /* we use the fact that cobases are sorted by index value */
2043 if (Q->debug)
2044 printA (D, Q);
2045
2046
2047
2048 while (C[0] <= m)
2049 {
2050 i = C[0];
2051 j = inequality[B[i] - lastdv];
2052 inequality[B[i] - lastdv] = inequality[C[0] - lastdv];
2053 inequality[C[0] - lastdv] = j;
2054 C[0] = B[i];
2055 B[i] = i;
2056 reorder1 (C, Col, ZERO, d);
2057 }
2058
2059
2060
2061 if (Q->debug)
2062 {
2063 fprintf (lrs_ofp, "\n*Inequality numbers for indices %ld .. %ld : ", lastdv + 1, m + d);
2064 for (i = 1; i <= m - nlinearity; i++)
2065 fprintf (lrs_ofp, " %ld ", inequality[i]);
2066 printA (D, Q);
2067 }
2068
2069
2070
2071 if (Q->restart)
2072 {
2073 if (Q->debug)
2074 fprintf (lrs_ofp, "\nPivoting to restart co-basis");
2075 if (!restartpivots (D, Q))
2076 return FALSE;
2077 D->lexflag = lexmin (D, Q, ZERO); /* see if lexmin basis */
2078 if (Q->debug)
2079 printA (D, Q);
2080 }
2081
2082 /* Check to see if necessary to resize */
2083 /* bug fix 2018.6.7 new value of d required below */
2084 if (Q->inputd > D->d)
2085 *D_p = resize (D, Q);
2086
2087
2088 lrs_clear_mp(Temp); lrs_clear_mp(scale);
2089 return TRUE;
2090 }
2091 /********* end of lrs_getfirstbasis ***************/
2092
2093 /*****************************************/
2094 /* getnextbasis in reverse search order */
2095 /*****************************************/
2096
2097
2098 long
lrs_getnextbasis(lrs_dic ** D_p,lrs_dat * Q,long backtrack)2099 lrs_getnextbasis (lrs_dic ** D_p, lrs_dat * Q, long backtrack)
2100 /* gets next reverse search tree basis, FALSE if none */
2101 /* switches to estimator if maxdepth set */
2102 /* backtrack TRUE means backtrack from here */
2103
2104 {
2105 /* assign local variables to structures */
2106 long i = 0L, j = 0L;
2107 long m = D->m;
2108 long d = D->d;
2109 long saveflag;
2110 long cob_est=0; /* estimated number of cobases in subtree from current node */
2111
2112
2113 if (backtrack && D->depth == 0)
2114 return FALSE; /* cannot backtrack from root */
2115
2116 if (Q->maxoutput > 0 && Q->count[0]+Q->count[1]-Q->hull >= Q->maxoutput)
2117 return FALSE; /* output limit reached */
2118
2119 while ((j < d) || (D->B[m] != m)) /*main while loop for getnextbasis */
2120 {
2121 if (D->depth >= Q->maxdepth)
2122 {
2123 if (Q->runs > 0 && !backtrack ) /*get an estimate of remaining tree */
2124 {
2125
2126 //2015.2.9 do iterative estimation backtracking when estimate is small
2127
2128 saveflag=Q->printcobasis;
2129 Q->printcobasis=FALSE;
2130 cob_est=lrs_estimate (D, Q);
2131 Q->printcobasis=saveflag;
2132 if(cob_est <= Q->subtreesize) /* stop iterative estimation */
2133 {
2134 if(Q->verbose && cob_est > 0) /* when zero we are at a leaf */
2135 { lrs_printcobasis(D,Q,ZERO);
2136 fprintf(lrs_ofp," cob_est= %ld *subtree",cob_est);
2137 }
2138 backtrack=TRUE;
2139 }
2140
2141 }
2142 else // either not estimating or we are backtracking
2143
2144 //2018.1.19 if (!backtrack && !Q->printcobasis)
2145 if (!backtrack )
2146 if(!lrs_leaf(D,Q)) /* 2015.6.5 cobasis returned if not a leaf */
2147 lrs_return_unexplored(D,Q);
2148
2149 backtrack = TRUE;
2150
2151
2152 if (Q->maxdepth == 0 && cob_est <= Q->subtreesize) /* root estimate only */
2153 return FALSE; /* no nextbasis */
2154 } // if (D->depth >= Q->maxdepth)
2155
2156
2157 /* if ( Q->truncate && negative(D->A[0][0]))*/ /* truncate when moving from opt. vertex */
2158 /* backtrack = TRUE; 2011.7.14 */
2159
2160 if (backtrack) /* go back to prev. dictionary, restore i,j */
2161 {
2162 backtrack = FALSE;
2163
2164 if (check_cache (D_p, Q, &i, &j))
2165 {
2166 if (Q->debug)
2167 fprintf (lrs_ofp, "\n Cached Dict. restored to depth %ld\n", D->depth);
2168 }
2169 else
2170 {
2171 D->depth--;
2172 selectpivot (D, Q, &i, &j);
2173 pivot (D, Q, i, j);
2174 update (D, Q, &i, &j); /*Update B,C,i,j */
2175 }
2176
2177 if (Q->debug)
2178 {
2179 fprintf (lrs_ofp, "\n Backtrack Pivot: indices i=%ld j=%ld depth=%ld", i, j, D->depth);
2180 printA (D, Q);
2181 };
2182
2183 j++; /* go to next column */
2184 } /* end of if backtrack */
2185
2186 if (D->depth < Q->mindepth)
2187 break;
2188
2189 /* try to go down tree */
2190
2191 /* 2011.7.14 patch */
2192 while ((j < d) && (!reverse (D, Q, &i, j) || (Q->truncate && Q->minratio[D->m]==1)))
2193 j++;
2194 if (j == d )
2195 backtrack = TRUE;
2196 else
2197 /*reverse pivot found */
2198 {
2199 cache_dict (D_p, Q, i, j);
2200 /* Note that the next two lines must come _after_ the
2201 call to cache_dict */
2202
2203 D->depth++;
2204 if (D->depth > Q->deepest)
2205 Q->deepest++;
2206
2207 pivot (D, Q, i, j);
2208 update (D, Q, &i, &j); /*Update B,C,i,j */
2209
2210 D->lexflag = lexmin (D, Q, ZERO); /* see if lexmin basis */
2211 Q->count[2]++;
2212 Q->totalnodes++;
2213
2214 if (Q->strace == Q->count[2])
2215 Q->debug = TRUE;
2216 if (Q->etrace == Q->count[2])
2217 Q->debug = FALSE;
2218 return TRUE; /*return new dictionary */
2219 }
2220
2221
2222 } /* end of main while loop for getnextbasis */
2223 return FALSE; /* done, no more bases */
2224 } /*end of lrs_getnextbasis */
2225
2226 /*************************************/
2227 /* print out one line of output file */
2228 /*************************************/
2229 long
lrs_getvertex(lrs_dic * P,lrs_dat * Q,lrs_mp_vector output)2230 lrs_getvertex (lrs_dic * P, lrs_dat * Q, lrs_mp_vector output)
2231 /*Print out current vertex if it is lexmin and return it in output */
2232 /* return FALSE if no output generated */
2233 {
2234 lrs_mp_matrix A = P->A;
2235 long i;
2236 long ind; /* output index */
2237 long ired; /* counts number of redundant columns */
2238 /* assign local variables to structures */
2239 long *redundcol = Q->redundcol;
2240 long *count = Q->count;
2241 long *B = P->B;
2242 long *Row = P->Row;
2243
2244 long lastdv = Q->lastdv;
2245
2246 long hull;
2247 long lexflag;
2248
2249
2250 hull = Q->hull;
2251 lexflag = P->lexflag;
2252 if (lexflag || Q->allbases)
2253 ++(Q->count[1]);
2254
2255 if (Q->debug)
2256 printA (P, Q);
2257
2258 if (Q->getvolume)
2259 {
2260 linint (Q->sumdet, 1, P->det, 1);
2261 updatevolume (P, Q);
2262 }
2263 if(Q->triangulation) /* this will print out a triangulation */
2264 lrs_printcobasis(P,Q,ZERO);
2265
2266
2267 /*print cobasis if printcobasis=TRUE and count[2] a multiple of frequency */
2268 /* or for lexmin basis, except origin for hull computation - ugly! */
2269
2270 if (Q->printcobasis)
2271 if ((lexflag && !hull) || ((Q->frequency > 0) && (count[2] == (count[2] / Q->frequency) * Q->frequency)))
2272 lrs_printcobasis(P,Q,ZERO);
2273
2274 if (hull)
2275 return FALSE; /* skip printing the origin */
2276
2277 if (!lexflag && !Q->allbases && !Q->lponly) /* not lexmin, and not printing forced */
2278 return FALSE;
2279
2280
2281 /* copy column 0 to output */
2282
2283 i = 1;
2284 ired = 0;
2285 copy (output[0], P->det);
2286
2287 for (ind = 1; ind < Q->n; ind++) /* extract solution */
2288
2289 if ((ired < Q->nredundcol) && (redundcol[ired] == ind))
2290 /* column was deleted as redundant */
2291 {
2292 itomp (ZERO, output[ind]);
2293 ired++;
2294 }
2295 else
2296 /* column not deleted as redundant */
2297 {
2298 getnextoutput (P, Q, i, ZERO, output[ind]);
2299 i++;
2300 }
2301
2302 reducearray (output, Q->n);
2303 if (lexflag && one(output[0]))
2304 ++Q->count[4]; /* integer vertex */
2305
2306
2307 /* uncomment to print nonzero basic variables
2308
2309 printf("\n nonzero basis: vars");
2310 for(i=1;i<=lastdv; i++)
2311 {
2312 if ( !zero(A[Row[i]][0]) )
2313 printf(" %ld ",B[i]);
2314 }
2315 */
2316
2317 /* printslack inequality indices */
2318
2319 if (Q->printslack)
2320 {
2321 fprintf(lrs_ofp,"\nslack ineq:");
2322 for(i=lastdv+1;i<=P->m; i++)
2323 {
2324 if (!zero(A[Row[i]][0]))
2325 fprintf(lrs_ofp," %ld ", Q->inequality[B[i]-lastdv]);
2326 }
2327 }
2328
2329
2330
2331 return TRUE;
2332 } /* end of lrs_getvertex */
2333
2334 long
lrs_getray(lrs_dic * P,lrs_dat * Q,long col,long redcol,lrs_mp_vector output)2335 lrs_getray (lrs_dic * P, lrs_dat * Q, long col, long redcol, lrs_mp_vector output)
2336 /*Print out solution in col and return it in output */
2337 /*redcol =n for ray/facet 0..n-1 for linearity column */
2338 /*hull=1 implies facets will be recovered */
2339 /* return FALSE if no output generated in column col */
2340 {
2341 long i;
2342 long ind; /* output index */
2343 long ired; /* counts number of redundant columns */
2344 /* assign local variables to structures */
2345 long *redundcol = Q->redundcol;
2346 long *count = Q->count;
2347 long hull = Q->hull;
2348 long n = Q->n;
2349
2350 long *B = P->B;
2351 long *Row = P->Row;
2352 long lastdv = Q->lastdv;
2353
2354 if (Q->debug)
2355 {
2356 printA (P, Q);
2357 for (i = 0; i < Q->nredundcol; i++)
2358 fprintf (lrs_ofp, " %ld", redundcol[i]);
2359 fflush(lrs_ofp);
2360 }
2361
2362 if (redcol == n)
2363 {
2364 ++count[0];
2365 if (Q->printcobasis)
2366 lrs_printcobasis(P,Q,col);
2367
2368 }
2369
2370 i = 1;
2371 ired = 0;
2372
2373 for (ind = 0; ind < n; ind++) /* print solution */
2374 {
2375 if (ind == 0 && !hull) /* must have a ray, set first column to zero */
2376 itomp (ZERO, output[0]);
2377
2378 else if ((ired < Q->nredundcol) && (redundcol[ired] == ind))
2379 /* column was deleted as redundant */
2380 {
2381 if (redcol == ind) /* true for linearity on this cobasic index */
2382 /* we print reduced determinant instead of zero */
2383 copy (output[ind], P->det);
2384 else
2385 itomp (ZERO, output[ind]);
2386 ired++;
2387 }
2388 else
2389 /* column not deleted as redundant */
2390 {
2391 getnextoutput (P, Q, i, col, output[ind]);
2392 i++;
2393 }
2394
2395 }
2396 reducearray (output, n);
2397 /* printslack for rays: 2006.10.10 */
2398 /* printslack inequality indices */
2399
2400 if (Q->printslack)
2401 {
2402 fprintf(lrs_ofp,"\nslack ineq:");
2403 for(i=lastdv+1;i<=P->m; i++)
2404 {
2405 if (!zero(P->A[Row[i]][col]))
2406 fprintf(lrs_ofp," %ld ", Q->inequality[B[i]-lastdv]);
2407 }
2408 }
2409
2410
2411
2412 return TRUE;
2413 } /* end of lrs_getray */
2414
2415
2416 void
getnextoutput(lrs_dic * P,lrs_dat * Q,long i,long col,lrs_mp out)2417 getnextoutput (lrs_dic * P, lrs_dat * Q, long i, long col, lrs_mp out)
2418 /* get A[B[i][col] and copy to out */
2419 {
2420 long row;
2421 long m = P->m;
2422 long d = P->d;
2423 long lastdv = Q->lastdv;
2424 lrs_mp_matrix A = P->A;
2425 long *B = P->B;
2426 long *Row = P->Row;
2427 long j;
2428
2429 if (i == d && Q->voronoi)
2430 return; /* skip last column if voronoi set */
2431
2432 row = Row[i];
2433
2434 if (Q->nonnegative) /* if m+i basic get correct value from dictionary */
2435 /* the slack for the inequality m-d+i contains decision */
2436 /* variable x_i. We first see if this is in the basis */
2437 /* otherwise the value of x_i is zero, except for a ray */
2438 /* when it is one (det/det) for the actual column it is in */
2439 {
2440 for (j = lastdv+ 1; j <= m; j++)
2441 {
2442 if ( Q->inequality[B[j]-lastdv] == m-d+i )
2443 {
2444 copy (out, A[Row[j]][col]);
2445 return;
2446 }
2447 }
2448 /* did not find inequality m-d+i in basis */
2449 if ( i == col )
2450 copy(out,P->det);
2451 else
2452 itomp(ZERO,out);
2453
2454 }
2455 else
2456 copy (out, A[row][col]);
2457
2458 } /* end of getnextoutput */
2459
2460
2461 void
lrs_printcobasis(lrs_dic * P,lrs_dat * Q,long col)2462 lrs_printcobasis (lrs_dic * P, lrs_dat * Q, long col)
2463 /* col is output column being printed */
2464 {
2465 char *ss, *sdet, *sin_det, *sz;
2466 long i;
2467 long rflag; /* used to find inequality number for ray column */
2468 /* assign local variables to structures */
2469 lrs_mp_matrix A = P->A;
2470 lrs_mp Nvol, Dvol; /* hold rescaled det of current basis */
2471 long *B = P->B;
2472 long *C = P->C;
2473 long *Col = P->Col;
2474 long *Row = P->Row;
2475 long *inequality = Q->inequality;
2476 long *temparray = Q->temparray;
2477 long *count = Q->count;
2478 long hull = Q->hull;
2479 long d = P->d;
2480 long lastdv = Q->lastdv;
2481 long m=P->m;
2482 long firstime=TRUE;
2483 long nincidence; /* count number of tight inequalities */
2484 long len=0;
2485
2486 lrs_alloc_mp(Nvol); lrs_alloc_mp(Dvol);
2487
2488 /* convert lrs_mp to char, compute length of string ss and malloc*/
2489
2490 sdet=cpmp(" det=", P->det);
2491
2492 rescaledet (P, Q, Nvol, Dvol); /* scales determinant in case input rational */
2493 sin_det=cprat("in_det=", Nvol,Dvol);
2494
2495 sz=cprat("z=", P->objnum, P->objden);
2496
2497 len = snprintf(NULL, 0, "%s%s%s", sdet,sin_det,sz);
2498
2499 ss=(char*)malloc(len+(d+m)*20);
2500
2501 /* build the printcobasis string */
2502
2503 len=0;
2504 if (hull)
2505 len=len+sprintf (ss+len, "F#%ld B#%ld h=%ld vertices/rays ", count[0], count[2], P->depth);
2506 else if (Q->voronoi)
2507 len=len+sprintf (ss+len, "V#%ld R#%ld B#%ld h=%ld data points ", count[1], count[0], count[2], P->depth);
2508 else
2509 len=len+sprintf (ss+len, "V#%ld R#%ld B#%ld h=%ld facets ", count[1], count[0], count[2], P->depth);
2510
2511 rflag = (-1);
2512 for (i = 0; i < d; i++)
2513 {
2514 temparray[i] = inequality[C[i] - lastdv];
2515 if (Col[i] == col)
2516 rflag = temparray[i]; /* look for ray index */
2517 }
2518 for (i = 0; i < d; i++)
2519 reorder (temparray, d);
2520 for (i = 0; i < d; i++)
2521 {
2522 len=len+sprintf (ss+len, " %ld", temparray[i]);
2523
2524 if (!(col == ZERO) && (rflag == temparray[i])) /* missing cobasis element for ray */
2525 len=len+sprintf (ss+len, "*");
2526
2527 }
2528
2529 /* get and print incidence information */
2530 if ( col == 0 )
2531 nincidence = d;
2532 else
2533 nincidence = d-1;
2534
2535 for(i=lastdv+1;i<=m;i++)
2536 if ( zero (A[Row[i]][0] ))
2537 if( ( col == ZERO ) || zero (A[Row[i]] [col]) )
2538 {
2539 nincidence++;
2540 if( Q->incidence )
2541 {
2542 if (firstime)
2543 {
2544 len=len+sprintf (ss+len," :");
2545 firstime = FALSE;
2546 }
2547 len=len+sprintf(ss+len," %ld",inequality[B[i] - lastdv ] );
2548 }
2549 }
2550
2551 len=len+sprintf(ss+len," I#%ld",nincidence);
2552
2553 sprintf (ss+len,"%s %s %s ",sdet,sin_det,sz);
2554
2555 if(Q->mplrs)
2556 lrs_post_output("cobasis", ss);
2557 else
2558 fprintf(lrs_ofp,"\n%s",ss);
2559
2560 free(ss); free(sdet); free(sin_det); free(sz);
2561 lrs_clear_mp(Nvol); lrs_clear_mp(Dvol);
2562
2563
2564 } /* end of lrs_printcobasis */
2565
2566
2567 /*********************/
2568 /* print final totals */
2569 /*********************/
2570 void
lrs_printtotals(lrs_dic * P,lrs_dat * Q)2571 lrs_printtotals (lrs_dic * P, lrs_dat * Q)
2572 {
2573 static int first_time=1;
2574 /* print warnings */
2575
2576 if(first_time)
2577 {
2578 first_time=0;
2579
2580 if (!Q->mplrs)
2581 fprintf(lrs_ofp,"\nend");
2582
2583 if (Q->dualdeg)
2584 {
2585 lrs_warning(Q,"finalwarn","*Warning: Starting dictionary is dual degenerate");
2586 lrs_warning(Q,"finalwarn","*Complete enumeration may not have been produced");
2587 if (Q->maximize)
2588 lrs_warning(Q,"finalwarn","*Recommendation: Add dualperturb option before maximize in input file\n");
2589 else
2590 lrs_warning(Q,"finalwarn","*Recommendation: Add dualperturb option before minimize in input file\n");
2591 }
2592
2593 if (Q->unbounded)
2594 {
2595 lrs_warning(Q,"finalwarn","*Warning: Starting dictionary contains rays");
2596 lrs_warning(Q,"finalwarn","*Complete enumeration may not have been produced");
2597 if (Q->maximize)
2598 lrs_warning(Q,"finalwarn","*Recommendation: Change or remove maximize option or add bounds\n");
2599 else
2600 lrs_warning(Q,"finalwarn","*Recommendation: Change or remove minimize option or add bounds\n");
2601 }
2602
2603 if (Q->truncate)
2604 lrs_warning(Q,"finalwarn","*Tree truncated at each new vertex");
2605 }
2606
2607 if(!Q->hull)
2608 {
2609 if ( Q->allbases )
2610 lrs_warning(Q,"finalwarn","*Note! Duplicate vertices/rays may be present");
2611 else if (Q->count[0] > 1 && !Q->homogeneous)
2612 lrs_warning(Q,"finalwarn","*Note! Duplicate rays may be present");
2613 }
2614
2615 if(Q->mplrs)
2616 {
2617 char *vol;
2618 if (Q->hull && Q->getvolume)
2619 {
2620 rescalevolume (P, Q, Q->Nvolume, Q->Dvolume);
2621 vol=cprat("",Q->Nvolume, Q->Dvolume);
2622 lrs_post_output("volume", vol);
2623 free(vol);
2624 }
2625 return;
2626 }
2627
2628 if(!Q->messages)
2629 return;
2630
2631 long i;
2632 double x;
2633 /* local assignments */
2634 double *cest = Q->cest;
2635 long *count = Q->count;
2636 long *inequality = Q->inequality;
2637 long *linearity = Q->linearity;
2638 long *temparray = Q->temparray;
2639
2640 long *C = P->C;
2641
2642 long hull = Q->hull;
2643 long homogeneous = Q->homogeneous;
2644 long nlinearity = Q->nlinearity;
2645 long nredundcol = Q->nredundcol;
2646 long d, lastdv;
2647 d = P->d;
2648 lastdv = Q->lastdv;
2649 if(Q->hull) /* count[1] stores the number of linearities */
2650 Q->count[1]=Q->nredundcol - Q->homogeneous;
2651
2652 /* warnings for lrs only */
2653 if (Q->maxdepth < MAXD)
2654 fprintf (lrs_ofp, "\n*Tree truncated at depth %lld", Q->maxdepth);
2655 if (Q->maxcobases > 0L)
2656 fprintf (lrs_ofp, "\n*maxcobases = %ld", Q->maxcobases-Q->startcount[2]);
2657 if (Q->maxoutput > 0L)
2658 fprintf (lrs_ofp, "\n*Maximum number of output lines = %ld", Q->maxoutput);
2659
2660
2661 /* next block with volume rescaling must come before estimates are printed */
2662
2663 if (Q->getvolume && Q->runs == 0)
2664 {
2665
2666 if( Q->debug)
2667 {
2668 fprintf (lrs_ofp, "\n*Sum of det(B)=");
2669 pmp ("", Q->sumdet);
2670 }
2671
2672 rescalevolume (P, Q, Q->Nvolume, Q->Dvolume);
2673
2674 if (Q->polytope)
2675 prat ("\n*Volume=", Q->Nvolume, Q->Dvolume);
2676 else
2677 prat ("\n*Pseudovolume=", Q->Nvolume, Q->Dvolume);
2678 }
2679
2680 if (hull) /* output things that are specific to hull computation */
2681 {
2682 fprintf (lrs_ofp, "\n*Totals: facets=%ld bases=%ld", count[0], count[2]);
2683
2684 if (count[1] > 0 )
2685 {
2686 fprintf (lrs_ofp, " linearities=%ld", count[1]);
2687 fprintf (lrs_ofp, " facets+linearities=%ld",count[1]+count[0]);
2688 }
2689 if(lrs_ofp != stdout)
2690 {
2691 printf ("\n*Totals: facets=%ld bases=%ld", count[0], count[2]);
2692
2693 if (count[1] > 0)
2694 {
2695 printf (" linearities=%ld", count[1]);
2696 printf (" facets+linearities=%ld",count[1]+count[0]);
2697 }
2698 }
2699
2700
2701 if(Q->runs > 0)
2702 {
2703 fprintf (lrs_ofp, "\n*Estimates: facets=%.0f bases=%.0f", count[0] + cest[0], count[2] + cest[2]);
2704 if (Q->getvolume)
2705 {
2706 rescalevolume (P, Q, Q->Nvolume, Q->Dvolume);
2707 rattodouble (Q->Nvolume, Q->Dvolume, &x);
2708 for (i = 2; i < d; i++)
2709 cest[3] = cest[3] / i; /*adjust for dimension */
2710 if(cest[3]==0)
2711 prat (" volume=", Q->Nvolume, Q->Dvolume);
2712 else
2713 fprintf (lrs_ofp, " volume=%g", cest[3] + x);
2714
2715 }
2716
2717 fprintf (lrs_ofp, "\n*Total number of tree nodes evaluated: %ld", Q->totalnodes);
2718 #ifndef TIMES
2719 fprintf (lrs_ofp, "\n*Estimated total running time=%.1f secs ",(count[2]+cest[2])/Q->totalnodes*get_time () );
2720 #endif
2721 }
2722
2723
2724 }
2725 else /* output things specific to vertex/ray computation */
2726 {
2727 fprintf (lrs_ofp, "\n*Totals: vertices=%ld rays=%ld bases=%ld", count[1], count[0], count[2]);
2728
2729 fprintf (lrs_ofp, " integer_vertices=%ld ",count[4]);
2730
2731
2732 if (nredundcol > 0)
2733 fprintf (lrs_ofp, " linearities=%ld", nredundcol);
2734 if ( count[0] + nredundcol > 0 )
2735 {
2736 fprintf (lrs_ofp, " vertices+rays");
2737 if ( nredundcol > 0 )
2738 fprintf (lrs_ofp, "+linearities");
2739 fprintf (lrs_ofp, "=%ld",nredundcol+count[0]+count[1]);
2740 }
2741
2742 if(lrs_ofp != stdout)
2743 {
2744 printf ("\n*Totals: vertices=%ld rays=%ld bases=%ld", count[1], count[0], count[2]);
2745
2746 printf (" integer_vertices=%ld ",count[4]);
2747
2748 if (nredundcol > 0)
2749 printf (" linearities=%ld", nredundcol);
2750 if ( count[0] + nredundcol > 0 )
2751 {
2752 printf (" vertices+rays");
2753 if ( nredundcol > 0 )
2754 printf ("+linearities");
2755 printf ("=%ld",nredundcol+count[0]+count[1]);
2756 }
2757 } /* end lrs_ofp != stdout */
2758
2759
2760 if (Q->runs > 0)
2761 {
2762 fprintf (lrs_ofp, "\n*Estimates: vertices=%.0f rays=%.0f", count[1]+cest[1], count[0]+cest[0]);
2763 fprintf (lrs_ofp, " bases=%.0f integer_vertices=%.0f ",count[2]+cest[2], count[4]+cest[4]);
2764
2765 if (Q->getvolume)
2766 {
2767 rattodouble (Q->Nvolume, Q->Dvolume, &x);
2768 for (i = 2; i <= d-homogeneous; i++)
2769 cest[3] = cest[3] / i; /*adjust for dimension */
2770 fprintf (lrs_ofp, " pseudovolume=%g", cest[3] + x);
2771 }
2772 fprintf (lrs_ofp, "\n*Total number of tree nodes evaluated: %ld", Q->totalnodes);
2773 #ifndef TIMES
2774 fprintf (lrs_ofp, "\n*Estimated total running time=%.1f secs ",(count[2]+cest[2])/Q->totalnodes*get_time () );
2775 #endif
2776 }
2777
2778 } /* end of output for vertices/rays */
2779
2780 fprintf (lrs_ofp, "\n*Dictionary Cache: max size= %ld misses= %ld/%ld Tree Depth= %ld", dict_count, cache_misses, cache_tries, Q->deepest);
2781 if(lrs_ofp != stdout)
2782 printf ("\n*Dictionary Cache: max size= %ld misses= %ld/%ld Tree Depth= %ld", dict_count, cache_misses, cache_tries, Q->deepest);
2783
2784 if(Q->debug)
2785 {
2786 fprintf (lrs_ofp, "\n*Input size m=%ld rows n=%ld columns", P->m, Q->n);
2787
2788 if (hull)
2789 fprintf (lrs_ofp, " working dimension=%ld", d - 1 + homogeneous);
2790 else
2791 fprintf (lrs_ofp, " working dimension=%ld", d);
2792
2793 fprintf (lrs_ofp, "\n*Starting cobasis defined by input rows");
2794 for (i = 0; i < nlinearity; i++)
2795 temparray[i] = linearity[i];
2796 for (i = nlinearity; i < lastdv; i++)
2797 temparray[i] = inequality[C[i - nlinearity] - lastdv];
2798 for (i = 0; i < lastdv; i++)
2799 reorder (temparray, lastdv);
2800 for (i = 0; i < lastdv; i++)
2801 fprintf (lrs_ofp, " %ld", temparray[i]);
2802 }
2803 return;
2804
2805
2806 } /* end of lrs_printtotals */
2807 /************************/
2808 /* Estimation function */
2809 /************************/
2810 long
lrs_estimate(lrs_dic * P,lrs_dat * Q)2811 lrs_estimate (lrs_dic * P, lrs_dat * Q)
2812 /*returns estimate of subtree size (no. cobases) from current node */
2813 /*current node is not counted. */
2814 /*cest[0]rays [1]vertices [2]bases [3]volume */
2815 /* [4] integer vertices */
2816 {
2817
2818 lrs_mp_vector output; /* holds one line of output; ray,vertex,facet,linearity */
2819 lrs_mp Nvol, Dvol; /* hold volume of current basis */
2820 long estdepth = 0; /* depth of basis/vertex in subtree for estimate */
2821 long i = 0, j = 0, k, nchild, runcount, col;
2822 double prod = 0.0;
2823 double cave[] =
2824 {0.0, 0.0, 0.0, 0.0, 0.0};
2825 double nvertices, nbases, nrays, nvol, nivertices;
2826 long rays = 0;
2827 double newvol = 0.0;
2828 /* assign local variables to structures */
2829 lrs_mp_matrix A = P->A;
2830 long *isave = Q->isave;
2831 long *jsave = Q->jsave;
2832 double *cest = Q->cest;
2833 long d = P->d;
2834 lrs_alloc_mp(Nvol); lrs_alloc_mp(Dvol);
2835 /* Main Loop of Estimator */
2836
2837 output = lrs_alloc_mp_vector (Q->n); /* output holds one line of output from dictionary */
2838
2839 for (runcount = 1; runcount <= Q->runs; runcount++)
2840 { /* runcount counts number of random probes */
2841 j = 0;
2842 nchild = 1;
2843 prod = 1;
2844 nvertices = 0.0;
2845 nbases = 0.0;
2846 nrays = 0.0;
2847 nvol = 0.0;
2848 nivertices =0.0;
2849
2850 while (nchild != 0) /* while not finished yet */
2851 {
2852
2853 nchild = 0;
2854 while (j < d)
2855 {
2856 if (reverse (P, Q, &i, j))
2857 {
2858 isave[nchild] = i;
2859 jsave[nchild] = j;
2860 nchild++;
2861 }
2862 j++;
2863 }
2864
2865 if (estdepth == 0 && nchild == 0)
2866 {
2867 cest[0] = cest[0] + rays; /* may be some rays here */
2868 lrs_clear_mp(Nvol); lrs_clear_mp(Dvol);
2869 lrs_clear_mp_vector(output, Q->n);
2870 return(0L); /*subtree is a leaf */
2871 }
2872
2873 prod = prod * nchild;
2874 nbases = nbases + prod;
2875 if (Q->debug)
2876 {
2877 fprintf (lrs_ofp, " degree= %ld ", nchild);
2878 fprintf (lrs_ofp, "\nPossible reverse pivots: i,j=");
2879 for (k = 0; k < nchild; k++)
2880 fprintf (lrs_ofp, "%ld,%ld ", isave[k], jsave[k]);
2881 }
2882
2883 if (nchild > 0) /*reverse pivot found choose random child */
2884 {
2885 k = myrandom (Q->seed, nchild);
2886 Q->seed = myrandom (Q->seed, 977L);
2887 i = isave[k];
2888 j = jsave[k];
2889 if (Q->debug)
2890 fprintf (lrs_ofp, " selected pivot k=%ld seed=%ld ", k, Q->seed);
2891 estdepth++;
2892 Q->totalnodes++; /* calculate total number of nodes evaluated */
2893 pivot (P, Q, i, j);
2894 update (P, Q, &i, &j); /*Update B,C,i,j */
2895 if (lexmin (P, Q, ZERO)) /* see if lexmin basis for vertex */
2896 {
2897 nvertices = nvertices + prod;
2898 /* integer vertex estimate */
2899 if( lrs_getvertex(P,Q,output))
2900 { --Q->count[1];
2901 if (one(output[0] ))
2902 { --Q->count[4];
2903 nivertices = nivertices + prod;
2904 }
2905 }
2906 }
2907
2908 rays = 0;
2909 for (col = 1; col <= d; col++)
2910 if (negative (A[0][col]) && (lrs_ratio (P, Q, col) == 0) && lexmin (P, Q, col))
2911 rays++;
2912 nrays = nrays + prod * rays; /* update ray info */
2913
2914 if (Q->getvolume)
2915 {
2916 rescaledet (P, Q, Nvol, Dvol); /* scales determinant in case input rational */
2917 rattodouble (Nvol, Dvol, &newvol);
2918 nvol = nvol + newvol * prod; /* adjusts volume for degree */
2919 }
2920 j = 0;
2921 }
2922 }
2923 cave[0] = cave[0] + nrays;
2924 cave[1] = cave[1] + nvertices;
2925 cave[2] = cave[2] + nbases;
2926 cave[3] = cave[3] + nvol;
2927 cave[4] = cave[4] + nivertices;
2928
2929 /* backtrack to root and do it again */
2930
2931 while (estdepth > 0)
2932 {
2933 estdepth = estdepth - 1;
2934 selectpivot (P, Q, &i, &j);
2935 pivot (P, Q, i, j);
2936 update (P, Q, &i, &j); /*Update B,C,i,j */
2937 /*fprintf(lrs_ofp,"\n0 +++"); */
2938 if (Q->debug)
2939 {
2940 fprintf (lrs_ofp, "\n Backtrack Pivot: indices i,j %ld %ld ", i, j);
2941 printA (P, Q);
2942 }
2943 j++;
2944 }
2945
2946 } /* end of for loop on runcount */
2947
2948 j=(long) cave[2]/Q->runs;
2949
2950 //2015.2.9 Do not update totals if we do iterative estimating and subtree is too big
2951 if(Q->subtreesize == 0 || j <= Q->subtreesize )
2952 for (i = 0; i < 5; i++)
2953 cest[i] = cave[i] / Q->runs + cest[i];
2954
2955
2956 lrs_clear_mp(Nvol); lrs_clear_mp(Dvol);
2957 lrs_clear_mp_vector(output, Q->n);
2958 return(j);
2959 } /* end of lrs_estimate */
2960
2961
2962 /*********************************/
2963 /* Internal functions */
2964 /*********************************/
2965 /* Basic Dictionary functions */
2966 /******************************* */
2967
2968 long
reverse(lrs_dic * P,lrs_dat * Q,long * r,long s)2969 reverse (lrs_dic * P, lrs_dat * Q, long *r, long s)
2970 /* find reverse indices */
2971 /* TRUE if B[*r] C[s] is a reverse lexicographic pivot */
2972 {
2973 long i, j, enter, row, col;
2974
2975 /* assign local variables to structures */
2976 lrs_mp_matrix A = P->A;
2977 long *B = P->B;
2978 long *C = P->C;
2979 long *Row = P->Row;
2980 long *Col = P->Col;
2981 long d = P->d;
2982
2983 enter = C[s];
2984 col = Col[s];
2985 if (Q->debug)
2986 {
2987 fprintf (lrs_ofp, "\n+reverse: col index %ld C %ld Col %ld ", s, enter, col);
2988 fflush (lrs_ofp);
2989 }
2990 if (!negative (A[0][col]))
2991 {
2992 if (Q->debug)
2993 fprintf (lrs_ofp, " Pos/Zero Cost Coeff");
2994 Q->minratio[P->m]=0; /* 2011.7.14 */
2995 return (FALSE);
2996 }
2997
2998 *r = lrs_ratio (P, Q, col);
2999 if (*r == 0) /* we have a ray */
3000 {
3001 if (Q->debug)
3002 fprintf (lrs_ofp, " Pivot col non-negative: ray found");
3003 Q->minratio[P->m]=0; /* 2011.7.14 */
3004 return (FALSE);
3005 }
3006
3007 row = Row[*r];
3008
3009 /* check cost row after "pivot" for smaller leaving index */
3010 /* ie. j s.t. A[0][j]*A[row][col] < A[0][col]*A[row][j] */
3011 /* note both A[row][col] and A[0][col] are negative */
3012
3013 for (i = 0; i < d && C[i] < B[*r]; i++)
3014 if (i != s)
3015 {
3016 j = Col[i];
3017 if (positive (A[0][j]) || negative (A[row][j])) /*or else sign test fails trivially */
3018 if ((!negative (A[0][j]) && !positive (A[row][j])) ||
3019 comprod (A[0][j], A[row][col], A[0][col], A[row][j]) == -1)
3020 { /*+ve cost found */
3021 if (Q->debug)
3022 {
3023 fprintf (lrs_ofp, "\nPositive cost found: index %ld C %ld Col %ld", i, C[i], j);
3024 fflush(lrs_ofp);
3025 }
3026 Q->minratio[P->m]=0; /* 2011.7.14 */
3027
3028 return (FALSE);
3029 }
3030 }
3031 if (Q->debug)
3032 {
3033 fprintf (lrs_ofp, "\n+end of reverse : indices r %ld s %ld \n", *r, s);
3034 fflush (stdout);
3035 }
3036 return (TRUE);
3037 } /* end of reverse */
3038
3039 long
selectpivot(lrs_dic * P,lrs_dat * Q,long * r,long * s)3040 selectpivot (lrs_dic * P, lrs_dat * Q, long *r, long *s)
3041 /* select pivot indices using lexicographic rule */
3042 /* returns TRUE if pivot found else FALSE */
3043 /* pivot variables are B[*r] C[*s] in locations Row[*r] Col[*s] */
3044 {
3045 long j, col;
3046 /* assign local variables to structures */
3047 lrs_mp_matrix A = P->A;
3048 long *Col = P->Col;
3049 long d = P->d;
3050
3051 *r = 0;
3052 *s = d;
3053 j = 0;
3054
3055 /*find positive cost coef */
3056 while ((j < d) && (!positive (A[0][Col[j]])))
3057 j++;
3058
3059 if (j < d) /* pivot column found! */
3060 {
3061 *s = j;
3062 col = Col[j];
3063
3064 /*find min index ratio */
3065 *r = lrs_ratio (P, Q, col);
3066 if (*r != 0)
3067 return (TRUE); /* unbounded */
3068 }
3069 return (FALSE);
3070 } /* end of selectpivot */
3071 /******************************************************* */
3072
3073 void
pivot(lrs_dic * P,lrs_dat * Q,long bas,long cob)3074 pivot (lrs_dic * P, lrs_dat * Q, long bas, long cob)
3075 /* Qpivot routine for array A */
3076 /* indices bas, cob are for Basis B and CoBasis C */
3077 /* corresponding to row Row[bas] and column */
3078 /* Col[cob] respectively */
3079 {
3080 long r, s;
3081 long i, j;
3082 lrs_mp Ns, Nt, Ars;
3083 /* assign local variables to structures */
3084
3085 lrs_mp_matrix A = P->A;
3086 long *B = P->B;
3087 long *C = P->C;
3088 long *Row = P->Row;
3089 long *Col = P->Col;
3090 long d, m_A;
3091
3092 lrs_alloc_mp(Ns); lrs_alloc_mp(Nt); lrs_alloc_mp(Ars);
3093
3094 d = P->d;
3095 m_A = P->m_A;
3096 Q->count[3]++; /* count the pivot */
3097
3098 r = Row[bas];
3099 s = Col[cob];
3100
3101 /* Ars=A[r][s] */
3102 if (Q->debug)
3103 {
3104 fprintf (lrs_ofp, "\n pivot B[%ld]=%ld C[%ld]=%ld ", bas, B[bas], cob, C[cob]);
3105 fflush (stdout);
3106 }
3107 copy (Ars, A[r][s]);
3108 storesign (P->det, sign (Ars)); /*adjust determinant to new sign */
3109
3110
3111 for (i = 0; i <= m_A; i++)
3112 if (i != r)
3113 for (j = 0; j <= d; j++)
3114 if (j != s)
3115 {
3116 /* A[i][j]=(A[i][j]*Ars-A[i][s]*A[r][j])/P->det; */
3117
3118 mulint (A[i][j], Ars, Nt);
3119 mulint (A[i][s], A[r][j], Ns);
3120 #ifdef LRSLONG
3121 unchecked_decint (Nt, Ns); /* overflow cannot happen */
3122 #else
3123 decint (Nt, Ns);
3124 #endif
3125 exactdivint (Nt, P->det, A[i][j]);
3126 } /* end if j .... */
3127
3128 if (sign (Ars) == POS)
3129 {
3130 for (j = 0; j <= d; j++) /* no need to change sign if Ars neg */
3131 /* A[r][j]=-A[r][j]; */
3132 if (!zero (A[r][j]))
3133 changesign (A[r][j]);
3134 } /* watch out for above "if" when removing this "}" ! */
3135 else
3136 for (i = 0; i <= m_A; i++)
3137 if (!zero (A[i][s]))
3138 changesign (A[i][s]);
3139
3140 /* A[r][s]=P->det; */
3141
3142 copy (A[r][s], P->det); /* restore old determinant */
3143 copy (P->det, Ars);
3144 storesign (P->det, POS); /* always keep positive determinant */
3145
3146
3147 if (Q->debug)
3148 {
3149 fprintf (lrs_ofp, " depth=%ld ", P->depth);
3150 pmp ("det=", P->det);
3151 fflush(stdout);
3152 }
3153 /* set the new rescaled objective function value */
3154
3155 mulint (P->det, Q->Lcm[0], P->objden);
3156 mulint (Q->Gcd[0], A[0][0], P->objnum);
3157
3158 if (!Q->maximize)
3159 changesign (P->objnum);
3160 if (zero (P->objnum))
3161 storesign (P->objnum, POS);
3162 else
3163 reduce (P->objnum,P->objden);
3164
3165
3166 lrs_clear_mp(Ns); lrs_clear_mp(Nt); lrs_clear_mp(Ars);
3167 } /* end of pivot */
3168
3169 long
primalfeasible(lrs_dic * P,lrs_dat * Q)3170 primalfeasible (lrs_dic * P, lrs_dat * Q)
3171 /* Do dual pivots to get primal feasibility */
3172 /* Note that cost row is all zero, so no ratio test needed for Dual Bland's rule */
3173 {
3174 long primalinfeasible = TRUE;
3175 long i, j;
3176 /* assign local variables to structures */
3177 lrs_mp_matrix A = P->A;
3178 long *Row = P->Row;
3179 long *Col = P->Col;
3180 long m, d, lastdv;
3181 m = P->m;
3182 d = P->d;
3183 lastdv = Q->lastdv;
3184
3185 /*temporary: try to get new start after linearity */
3186
3187 while (primalinfeasible)
3188 {
3189 i=lastdv+1;
3190 while (i <= m && !negative (A[Row[i]][0]) )
3191 i++;
3192 if (i <= m )
3193 {
3194 j = 0; /*find a positive entry for in row */
3195 while (j < d && !positive (A[Row[i]][Col[j]]))
3196 j++;
3197 if (j >= d)
3198 return (FALSE); /* no positive entry */
3199 pivot (P, Q, i, j);
3200 update (P, Q, &i, &j);
3201 }
3202 else
3203 primalinfeasible = FALSE;
3204 } /* end of while primalinfeasibile */
3205 return (TRUE);
3206 } /* end of primalfeasible */
3207
3208
3209 long
lrs_solvelp(lrs_dic * P,lrs_dat * Q,long maximize)3210 lrs_solvelp (lrs_dic * P, lrs_dat * Q, long maximize)
3211 /* Solve primal feasible lp by Dantzig`s rule and lexicographic ratio test */
3212 /* return TRUE if bounded, FALSE if unbounded */
3213 {
3214 long i, j, k=0L;
3215 long notdone=TRUE;
3216 /* assign local variables to structures */
3217 long d = P->d;
3218
3219 /* lponly=1 Dantzig, =2 random, =3 hybrid, =4 Bland */
3220
3221 if(Q->lponly <=1) /* Dantzig's rule */
3222 while (dan_selectpivot (P, Q, &i, &j))
3223 {
3224 pivot (P, Q, i, j);
3225 update (P, Q, &i, &j); /*Update B,C,i,j */
3226 }
3227
3228 if(Q->lponly ==2) /* random edge rule */
3229 while (ran_selectpivot (P, Q, &i, &j))
3230 {
3231 pivot (P, Q, i, j);
3232 update (P, Q, &i, &j); /*Update B,C,i,j */
3233 }
3234
3235 if(Q->lponly ==3) /* alternate Dantzig/randome rules */
3236 while (notdone)
3237 {
3238 if(k % 2) /* odd for dantzig even for random */
3239 notdone=dan_selectpivot (P, Q, &i, &j);
3240 else
3241 notdone=ran_selectpivot (P, Q, &i, &j);
3242
3243 if(notdone)
3244 {
3245 pivot (P, Q, i, j);
3246 update (P, Q, &i, &j); /*Update B,C,i,j */
3247 }
3248 k++;
3249 }
3250
3251 if(Q->lponly ==4) /* Bland's rule - used for vertex enumeration */
3252 while (selectpivot (P, Q, &i, &j))
3253 {
3254 pivot (P, Q, i, j);
3255 update (P, Q, &i, &j); /*Update B,C,i,j */
3256 }
3257
3258
3259 if (Q->debug)
3260 printA (P, Q);
3261
3262 if (j < d && i == 0) /* selectpivot gives information on unbounded solution */
3263 {
3264 if (Q->lponly && Q->messages)
3265 fprintf (lrs_ofp, "\n*Unbounded solution");
3266 return FALSE;
3267 }
3268 return TRUE;
3269 } /* end of lrs_solvelp */
3270
3271 long
getabasis(lrs_dic * P,lrs_dat * Q,long order[])3272 getabasis (lrs_dic * P, lrs_dat * Q, long order[])
3273
3274 /* Pivot Ax<=b to standard form */
3275 /*Try to find a starting basis by pivoting in the variables x[1]..x[d] */
3276 /*If there are any input linearities, these appear first in order[] */
3277 /* Steps: (a) Try to pivot out basic variables using order */
3278 /* Stop if some linearity cannot be made to leave basis */
3279 /* (b) Permanently remove the cobasic indices of linearities */
3280 /* (c) If some decision variable cobasic, it is a linearity, */
3281 /* and will be removed. */
3282
3283 {
3284 long i, j, k;
3285 /* assign local variables to structures */
3286 lrs_mp_matrix A = P->A;
3287 long *B = P->B;
3288 long *C = P->C;
3289 long *Row = P->Row;
3290 long *Col = P->Col;
3291 long *linearity = Q->linearity;
3292 long *redundcol = Q->redundcol;
3293 long m, d, nlinearity;
3294 long nredundcol = 0L; /* will be calculated here */
3295 char mess[100];
3296 m = P->m;
3297 d = P->d;
3298 nlinearity = Q->nlinearity;
3299
3300 if (Q->debug)
3301 {
3302 fprintf (lrs_ofp, "\ngetabasis from inequalities given in order");
3303 for (i = 0l; i < m; i++)
3304 fprintf (lrs_ofp, " %ld", order[i]);
3305 }
3306 for (j = 0l; j < m; j++)
3307 {
3308 i = 0l;
3309 while (i <= m && B[i] != d + order[j])
3310 i++; /* find leaving basis index i */
3311 if (j < nlinearity && i > m) /* cannot pivot linearity to cobasis */
3312 {
3313 if (Q->debug)
3314 printA (P, Q);
3315 if(Q->messages)
3316 fprintf (lrs_ofp, "\nCannot find linearity in the basis");
3317 return FALSE;
3318 }
3319 if (i <= m)
3320 { /* try to do a pivot */
3321 k = 0l;
3322 while (C[k] <= d && zero (A[Row[i]][Col[k]])){
3323 k++;
3324 }
3325 if (C[k] <= d)
3326 {
3327
3328 pivot (P, Q, i, k);
3329 update (P, Q, &i, &k);
3330 }
3331 else if (j < nlinearity)
3332 { /* cannot pivot linearity to cobasis */
3333 if (zero (A[Row[i]][0]))
3334 {
3335 if(Q->messages && overflow != 2)
3336 {
3337 sprintf (mess,"*Input linearity in row %ld is redundant--converted to inequality", order[j]);
3338 lrs_warning(Q,"warning",mess);
3339 }
3340 linearity[j]=0l;
3341 Q->redineq[j]=1; /* check for redundancy if running redund */
3342 }
3343 else
3344 {
3345 if (Q->debug)
3346 printA (P, Q);
3347 lrs_warning(Q,"warning","*No feasible solution");
3348 return FALSE;
3349 }
3350 }
3351
3352
3353 }
3354 }
3355
3356
3357 /* update linearity array to get rid of redundancies */
3358 i = 0;
3359 k = 0; /* counters for linearities */
3360 while (k < nlinearity)
3361 {
3362 while (k < nlinearity && linearity[k] == 0)
3363 k++;
3364 if (k < nlinearity)
3365 linearity[i++] = linearity[k++];
3366 }
3367
3368 nlinearity = i;
3369 /* bug fix, 2009.6.27 */ Q->nlinearity = i;
3370
3371 /* column dependencies now can be recorded */
3372 /* redundcol contains input column number 0..n-1 where redundancy is */
3373 k = 0;
3374 while (k < d && C[k] <= d)
3375 {
3376 if (C[k] <= d){ /* decision variable still in cobasis */
3377 redundcol[nredundcol++] = C[k] - Q->hull; /* adjust for hull indices */
3378
3379 }
3380 k++;
3381 }
3382
3383 /* now we know how many decision variables remain in problem */
3384 Q->nredundcol = nredundcol;
3385 Q->lastdv = d - nredundcol;
3386 if (Q->debug)
3387 {
3388 fprintf (lrs_ofp, "\nend of first phase of getabasis: ");
3389 fprintf (lrs_ofp, "lastdv=%ld nredundcol=%ld", Q->lastdv, Q->nredundcol);
3390 fprintf (lrs_ofp, "\nredundant cobases:");
3391 for (i = 0; i < nredundcol; i++)
3392 fprintf (lrs_ofp, " %ld", redundcol[i]);
3393 printA (P, Q);
3394 }
3395
3396 /* Remove linearities from cobasis for rest of computation */
3397 /* This is done in order so indexing is not screwed up */
3398
3399 for (i = 0; i < nlinearity; i++)
3400 { /* find cobasic index */
3401 k = 0;
3402 while (k < d && C[k] != linearity[i] + d)
3403 k++;
3404 if (k >= d)
3405 {
3406 lrs_warning(Q,"warning","\nError removing linearity");
3407 return FALSE;
3408 }
3409 if (!removecobasicindex (P, Q, k))
3410 return FALSE;
3411 d = P->d;
3412 }
3413 if (Q->debug && nlinearity > 0)
3414 printA (P, Q);
3415
3416 /* Check feasability */
3417 if (Q->givenstart)
3418 {
3419 i = Q->lastdv + 1;
3420 while (i <= m && !negative (A[Row[i]][0]))
3421 i++;
3422 if (i <= m)
3423 fprintf (lrs_ofp, "\n*Infeasible startingcobasis - will be modified");
3424 }
3425 return TRUE;
3426 } /* end of getabasis */
3427
3428 long
removecobasicindex(lrs_dic * P,lrs_dat * Q,long k)3429 removecobasicindex (lrs_dic * P, lrs_dat * Q, long k)
3430 /* remove the variable C[k] from the problem */
3431 /* used after detecting column dependency */
3432 {
3433 long i, j, cindex, deloc;
3434 /* assign local variables to structures */
3435 lrs_mp_matrix A = P->A;
3436 long *B = P->B;
3437 long *C = P->C;
3438 long *Col = P->Col;
3439 long m, d;
3440 m = P->m;
3441 d = P->d;
3442
3443 if (Q->debug)
3444 fprintf (lrs_ofp, "\nremoving cobasic index k=%ld C[k]=%ld Col[k]=%ld", k, C[k],Col[k]);
3445 cindex = C[k]; /* cobasic index to remove */
3446 deloc = Col[k]; /* matrix column location to remove */
3447
3448 for (i = 1; i <= m; i++) /* reduce basic indices by 1 after index */
3449 if (B[i] > cindex)
3450 B[i]--;
3451
3452 for (j = k; j < d; j++) /* move down other cobasic variables */
3453 {
3454 C[j] = C[j + 1] - 1; /* cobasic index reduced by 1 */
3455 Col[j] = Col[j + 1];
3456 }
3457
3458 if (deloc != d)
3459 {
3460 /* copy col d to deloc */
3461 for (i = 0; i <= m; i++)
3462 copy (A[i][deloc], A[i][d]);
3463
3464 /* reassign location for moved column */
3465 j = 0;
3466 while (Col[j] != d)
3467 j++;
3468
3469 Col[j] = deloc;
3470 }
3471
3472 P->d--;
3473 if (Q->debug)
3474 printA (P, Q);
3475 return TRUE;
3476 } /* end of removecobasicindex */
3477
3478 lrs_dic *
resize(lrs_dic * P,lrs_dat * Q)3479 resize (lrs_dic * P, lrs_dat * Q)
3480 /* resize the dictionary after some columns are deleted, ie. inputd>d */
3481 /* a new lrs_dic record is created with reduced size, and items copied over */
3482 {
3483 lrs_dic *P1; /* to hold new dictionary in case of resizing */
3484
3485 long i, j;
3486 long m, d, m_A;
3487
3488
3489 m = P->m;
3490 d = P->d;
3491 m_A = P->m_A;
3492
3493 /* get new dictionary record */
3494
3495 P1 = new_lrs_dic (m, d, m_A);
3496
3497 /* copy data from P to P1 */
3498 P1->i = P->i;
3499 P1->j = P->j;
3500 P1->depth = P->depth;
3501 P1->m = P->m;
3502 P1->d = P1->d_orig = d;
3503 P1->lexflag = P->lexflag;
3504 P1->m_A = P->m_A;
3505 copy (P1->det, P->det);
3506 copy (P1->objnum, P->objnum);
3507 copy (P1->objden, P->objden);
3508
3509 for (i = 0; i <= m; i++)
3510 {
3511 P1->B[i] = P->B[i];
3512 P1->Row[i] = P->Row[i];
3513 }
3514 for (i = 0; i <= m_A; i++)
3515 {
3516 for (j = 0; j <= d; j++)
3517 copy (P1->A[i][j], P->A[i][j]);
3518 }
3519
3520
3521 for (j = 0; j <= d; j++)
3522 {
3523 P1->Col[j] = P->Col[j];
3524 P1->C[j] = P->C[j];
3525 }
3526
3527 if (Q->debug)
3528 {
3529 fprintf (lrs_ofp, "\nDictionary resized from d=%ld to d=%ld", Q->inputd, P->d);
3530 printA (P1, Q);
3531 }
3532
3533 lrs_free_dic (P,Q);
3534
3535 /* Reassign cache pointers */
3536
3537 Q->Qhead = P1;
3538 Q->Qtail = P1;
3539 P1->next = P1;
3540 P1->prev = P1;
3541
3542 return P1;
3543
3544 }
3545 /********* resize ***************/
3546
3547
3548 long
restartpivots(lrs_dic * P,lrs_dat * Q)3549 restartpivots (lrs_dic * P, lrs_dat * Q)
3550 /* facet contains a list of the inequalities in the cobasis for the restart */
3551 /* inequality contains the relabelled inequalities after initialization */
3552 {
3553 long i, j, k;
3554 long *Cobasic; /* when restarting, Cobasic[j]=1 if j is in cobasis */
3555 /* assign local variables to structures */
3556 lrs_mp_matrix A = P->A;
3557 long *B = P->B;
3558 long *C = P->C;
3559 long *Row = P->Row;
3560 long *Col = P->Col;
3561 long *inequality = Q->inequality;
3562 long *facet = Q->facet;
3563 long nlinearity = Q->nlinearity;
3564 long m, d, lastdv;
3565 m = P->m;
3566 d = P->d;
3567 lastdv = Q->lastdv;
3568
3569 Cobasic = (long *) CALLOC ((unsigned) m + d + 2, sizeof (long));
3570
3571 if (Q->debug)
3572 fprintf(lrs_ofp,"\nCobasic flags in restartpivots");
3573 /* set Cobasic flags */
3574 for (i = 0; i < m + d + 1; i++)
3575 Cobasic[i] = 0;
3576 for (i = 0; i < d; i++) /* find index corresponding to facet[i] */
3577 {
3578 j = 1;
3579 while (facet[i +nlinearity] != inequality[j])
3580 j++;
3581 Cobasic[j + lastdv] = 1;
3582 if (Q->debug)
3583 fprintf(lrs_ofp," %ld %ld;",facet[i+nlinearity],j+lastdv);
3584 }
3585
3586 /* Note that the order of doing the pivots is important, as */
3587 /* the B and C vectors are reordered after each pivot */
3588
3589 /* Suggested new code from db starts */
3590 i=m;
3591 while (i>d){
3592 while(Cobasic[B[i]]){
3593 k = d - 1;
3594 while ((k >= 0) && (zero (A[Row[i]][Col[k]]) || Cobasic[C[k]])){
3595 k--;
3596 }
3597 if (k >= 0) {
3598 /*db asks: should i really be modified here? (see old code) */
3599 /*da replies: modifying i only makes is larger, and so */
3600 /*the second while loop will put it back where it was */
3601 /*faster (and safer) as done below */
3602 long ii=i;
3603 pivot (P, Q, ii, k);
3604 update (P, Q, &ii, &k);
3605 } else {
3606 lrs_warning(Q,"warning","\nInvalid Co-basis - does not have correct rank");
3607 free(Cobasic);
3608 return FALSE;
3609 }
3610 }
3611 i--;
3612 }
3613 /* Suggested new code from db ends */
3614
3615 /* check restarting from a primal feasible dictionary */
3616 for (i = lastdv + 1; i <= m; i++)
3617 if (negative (A[Row[i]][0]))
3618 {
3619 lrs_warning(Q,"warning","\nTrying to restart from infeasible dictionary");
3620 free(Cobasic);
3621 return FALSE;
3622 }
3623 free(Cobasic);
3624 return TRUE;
3625
3626 } /* end of restartpivots */
3627
3628
3629 long
lrs_ratio(lrs_dic * P,lrs_dat * Q,long col)3630 lrs_ratio (lrs_dic * P, lrs_dat * Q, long col) /*find lex min. ratio */
3631 /* find min index ratio -aig/ais, ais<0 */
3632 /* if multiple, checks successive basis columns */
3633 /* recoded Dec 1997 */
3634 {
3635 long i, j, comp, ratiocol, basicindex, start, nstart, cindex, bindex;
3636 long firstime; /*For ratio test, true on first pass,else false */
3637 lrs_mp Nmin, Dmin;
3638 long degencount, ndegencount;
3639 /* assign local variables to structures */
3640 lrs_mp_matrix A = P->A;
3641 long *B = P->B;
3642 long *Row = P->Row;
3643 long *Col = P->Col;
3644 long *minratio = Q->minratio;
3645 long m, d, lastdv;
3646
3647 m = P->m;
3648 d = P->d;
3649 lastdv = Q->lastdv;
3650
3651
3652 nstart=0;
3653 ndegencount=0;
3654 degencount = 0;
3655 minratio[P->m]=1; /*2011.7.14 non-degenerate pivot flag */
3656
3657 for (j = lastdv + 1; j <= m; j++)
3658 {
3659 /* search rows with negative coefficient in dictionary */
3660 /* minratio contains indices of min ratio cols */
3661 if (negative (A[Row[j]][col]))
3662 {
3663 minratio[degencount++] = j;
3664 if(zero (A[Row[j]][0]))
3665 minratio[P->m]=0; /*2011.7.14 degenerate pivot flag */
3666 }
3667 } /* end of for loop */
3668 if (Q->debug)
3669 {
3670 fprintf (lrs_ofp, " Min ratios: ");
3671 for (i = 0; i < degencount; i++)
3672 fprintf (lrs_ofp, " %ld ", B[minratio[i]]);
3673 }
3674 if (degencount == 0)
3675 return (degencount); /* non-negative pivot column */
3676
3677 lrs_alloc_mp(Nmin); lrs_alloc_mp(Dmin);
3678 ratiocol = 0; /* column being checked, initially rhs */
3679 start = 0; /* starting location in minratio array */
3680 bindex = d + 1; /* index of next basic variable to consider */
3681 cindex = 0; /* index of next cobasic variable to consider */
3682 basicindex = d; /* index of basis inverse for current ratio test, except d=rhs test */
3683 while (degencount > 1) /*keep going until unique min ratio found */
3684 {
3685 if (B[bindex] == basicindex) /* identity col in basis inverse */
3686 {
3687 if (minratio[start] == bindex)
3688 /* remove this index, all others stay */
3689 {
3690 start++;
3691 degencount--;
3692 }
3693 bindex++;
3694 }
3695 else
3696 /* perform ratio test on rhs or column of basis inverse */
3697 {
3698 firstime = TRUE;
3699 /*get next ratio column and increment cindex */
3700 if (basicindex != d)
3701 ratiocol = Col[cindex++];
3702 for (j = start; j < start + degencount; j++)
3703 {
3704 i = Row[minratio[j]]; /* i is the row location of the next basic variable */
3705 comp = 1; /* 1: lhs>rhs; 0:lhs=rhs; -1: lhs<rhs */
3706 if (firstime)
3707 firstime = FALSE; /*force new min ratio on first time */
3708 else
3709 {
3710 if (positive (Nmin) || negative (A[i][ratiocol]))
3711 {
3712 if (negative (Nmin) || positive (A[i][ratiocol]))
3713 comp = comprod (Nmin, A[i][col], A[i][ratiocol], Dmin);
3714 else
3715 comp = -1;
3716 }
3717
3718 else if (zero (Nmin) && zero (A[i][ratiocol]))
3719 comp = 0;
3720
3721 if (ratiocol == ZERO)
3722 comp = -comp; /* all signs reversed for rhs */
3723 }
3724 if (comp == 1)
3725 { /*new minimum ratio */
3726 nstart = j;
3727 copy (Nmin, A[i][ratiocol]);
3728 copy (Dmin, A[i][col]);
3729 ndegencount = 1;
3730 }
3731 else if (comp == 0) /* repeated minimum */
3732 minratio[nstart + ndegencount++] = minratio[j];
3733
3734 } /* end of for (j=start.... */
3735 degencount = ndegencount;
3736 start = nstart;
3737 } /* end of else perform ratio test statement */
3738 basicindex++; /* increment column of basis inverse to check next */
3739 if (Q->debug)
3740 {
3741 fprintf (lrs_ofp, " ratiocol=%ld degencount=%ld ", ratiocol, degencount);
3742 fprintf (lrs_ofp, " Min ratios: ");
3743 for (i = start; i < start + degencount; i++)
3744 fprintf (lrs_ofp, " %ld ", B[minratio[i]]);
3745 }
3746 } /*end of while loop */
3747 lrs_clear_mp(Nmin); lrs_clear_mp(Dmin);
3748 return (minratio[start]);
3749 } /* end of ratio */
3750
3751
3752
3753 long
lexmin(lrs_dic * P,lrs_dat * Q,long col)3754 lexmin (lrs_dic * P, lrs_dat * Q, long col)
3755 /*test if basis is lex-min for vertex or ray, if so TRUE */
3756 /* FALSE if a_r,g=0, a_rs !=0, r > s */
3757 {
3758 /*do lexmin test for vertex if col=0, otherwise for ray */
3759 long r, s, i, j;
3760 /* assign local variables to structures */
3761 lrs_mp_matrix A = P->A;
3762 long *B = P->B;
3763 long *C = P->C;
3764 long *Row = P->Row;
3765 long *Col = P->Col;
3766 long m = P->m;
3767 long d = P->d;
3768
3769 for (i = Q->lastdv + 1; i <= m; i++)
3770 {
3771 r = Row[i];
3772 if (zero (A[r][col])) /* necessary for lexmin to fail */
3773 for (j = 0; j < d; j++)
3774 {
3775 s = Col[j];
3776 if (B[i] > C[j]) /* possible pivot to reduce basis */
3777 {
3778 if (zero (A[r][0])) /* no need for ratio test, any pivot feasible */
3779 {
3780 if (!zero (A[r][s]))
3781 return (FALSE);
3782 }
3783 else if (negative (A[r][s]) && ismin (P, Q, r, s))
3784 {
3785 return (FALSE);
3786 }
3787 } /* end of if B[i] ... */
3788 }
3789 }
3790 if ((col != ZERO) && Q->debug)
3791 {
3792 fprintf (lrs_ofp, "\n lexmin ray in col=%ld ", col);
3793 printA (P, Q);
3794 }
3795 return (TRUE);
3796 } /* end of lexmin */
3797
3798 long
ismin(lrs_dic * P,lrs_dat * Q,long r,long s)3799 ismin (lrs_dic * P, lrs_dat * Q, long r, long s)
3800 /*test if A[r][s] is a min ratio for col s */
3801 {
3802 long i;
3803 /* assign local variables to structures */
3804 lrs_mp_matrix A = P->A;
3805 long m_A = P->m_A;
3806
3807 for (i = 1; i <= m_A; i++)
3808 if ((i != r) &&
3809 negative (A[i][s]) && comprod (A[i][0], A[r][s], A[i][s], A[r][0]))
3810 {
3811 return (FALSE);
3812 }
3813
3814 return (TRUE);
3815 }
3816
3817 void
update(lrs_dic * P,lrs_dat * Q,long * i,long * j)3818 update (lrs_dic * P, lrs_dat * Q, long *i, long *j)
3819 /*update the B,C arrays after a pivot */
3820 /* involving B[bas] and C[cob] */
3821 {
3822
3823 long leave, enter;
3824 /* assign local variables to structures */
3825 long *B = P->B;
3826 long *C = P->C;
3827 long *Row = P->Row;
3828 long *Col = P->Col;
3829 long m = P->m;
3830 long d = P->d;
3831
3832 leave = B[*i];
3833 enter = C[*j];
3834 B[*i] = enter;
3835 reorder1 (B, Row, *i, m + 1);
3836 C[*j] = leave;
3837 reorder1 (C, Col, *j, d);
3838 /* restore i and j to new positions in basis */
3839 for (*i = 1; B[*i] != enter; (*i)++); /*Find basis index */
3840 for (*j = 0; C[*j] != leave; (*j)++); /*Find co-basis index */
3841 if (Q->debug)
3842 printA(P,Q);
3843 } /* end of update */
3844
3845 long
lrs_degenerate(lrs_dic * P,lrs_dat * Q)3846 lrs_degenerate (lrs_dic * P, lrs_dat * Q)
3847 /* TRUE if the current dictionary is primal degenerate */
3848 /* not thoroughly tested 2000/02/15 */
3849 {
3850 long i;
3851 long *Row;
3852
3853 lrs_mp_matrix A = P->A;
3854 long d = P->d;
3855 long m = P->m;
3856
3857 Row = P->Row;
3858
3859 for (i = d + 1; i <= m; i++)
3860 if (zero (A[Row[i]][0]))
3861 return TRUE;
3862
3863 return FALSE;
3864 }
3865
3866
3867 /*********************************************************/
3868 /* Miscellaneous */
3869 /******************************************************* */
3870
3871 void
reorder(long a[],long range)3872 reorder (long a[], long range)
3873 /*reorder array in increasing order with one misplaced element */
3874 {
3875 long i, temp;
3876 for (i = 0; i < range - 1; i++)
3877 if (a[i] > a[i + 1])
3878 {
3879 temp = a[i];
3880 a[i] = a[i + 1];
3881 a[i + 1] = temp;
3882 }
3883 for (i = range - 2; i >= 0; i--)
3884 if (a[i] > a[i + 1])
3885 {
3886 temp = a[i];
3887 a[i] = a[i + 1];
3888 a[i + 1] = temp;
3889 }
3890
3891 } /* end of reorder */
3892
3893 void
reorder1(long a[],long b[],long newone,long range)3894 reorder1 (long a[], long b[], long newone, long range)
3895 /*reorder array a in increasing order with one misplaced element at index newone */
3896 /*elements of array b are updated to stay aligned with a */
3897 {
3898 long temp;
3899 while (newone > 0 && a[newone] < a[newone - 1])
3900 {
3901 temp = a[newone];
3902 a[newone] = a[newone - 1];
3903 a[newone - 1] = temp;
3904 temp = b[newone];
3905 b[newone] = b[newone - 1];
3906 b[--newone] = temp;
3907 }
3908 while (newone < range - 1 && a[newone] > a[newone + 1])
3909 {
3910 temp = a[newone];
3911 a[newone] = a[newone + 1];
3912 a[newone + 1] = temp;
3913 temp = b[newone];
3914 b[newone] = b[newone + 1];
3915 b[++newone] = temp;
3916 }
3917 } /* end of reorder1 */
3918
3919
3920 void
rescaledet(lrs_dic * P,lrs_dat * Q,lrs_mp Vnum,lrs_mp Vden)3921 rescaledet (lrs_dic * P, lrs_dat * Q, lrs_mp Vnum, lrs_mp Vden)
3922 /* rescale determinant to get its volume */
3923 /* Vnum/Vden is volume of current basis */
3924 {
3925 lrs_mp gcdprod; /* to hold scale factors */
3926 long i;
3927 /* assign local variables to structures */
3928 long *C = P->C;
3929 long *B = P->B;
3930 long m, d, lastdv;
3931
3932 lrs_alloc_mp(gcdprod);
3933 m = P->m;
3934 d = P->d;
3935 lastdv = Q->lastdv;
3936
3937 itomp (ONE, gcdprod);
3938 itomp (ONE, Vden);
3939
3940 for (i = 0; i < d; i++)
3941 if (B[i] <= m)
3942 {
3943 mulint (Q->Gcd[Q->inequality[C[i] - lastdv]], gcdprod, gcdprod);
3944 mulint (Q->Lcm[Q->inequality[C[i] - lastdv]], Vden, Vden);
3945 }
3946 mulint (P->det, gcdprod, Vnum);
3947 // reduce (Vnum, Vden);
3948 lrs_clear_mp(gcdprod);
3949 } /* end rescaledet */
3950
3951 void
rescalevolume(lrs_dic * P,lrs_dat * Q,lrs_mp Vnum,lrs_mp Vden)3952 rescalevolume (lrs_dic * P, lrs_dat * Q, lrs_mp Vnum, lrs_mp Vden)
3953 /* adjust volume for dimension */
3954 {
3955 lrs_mp temp, dfactorial;
3956 /* assign local variables to structures */
3957 long lastdv = Q->lastdv;
3958
3959 lrs_alloc_mp(temp); lrs_alloc_mp(dfactorial);
3960
3961 /*reduce Vnum by d factorial */
3962 getfactorial (dfactorial, lastdv);
3963 mulint (dfactorial, Vden, Vden);
3964 if (Q->hull && !Q->homogeneous)
3965 { /* For hull option multiply by d to correct for lifting */
3966 itomp (lastdv, temp);
3967 mulint (temp, Vnum, Vnum);
3968 }
3969
3970 reduce (Vnum, Vden);
3971 lrs_clear_mp(temp); lrs_clear_mp(dfactorial);
3972 }
3973
3974
3975 void
updatevolume(lrs_dic * P,lrs_dat * Q)3976 updatevolume (lrs_dic * P, lrs_dat * Q) /* rescale determinant and update the volume */
3977 {
3978 lrs_mp tN, tD, Vnum, Vden;
3979 lrs_alloc_mp(tN); lrs_alloc_mp(tD); lrs_alloc_mp(Vnum); lrs_alloc_mp(Vden);
3980 rescaledet (P, Q, Vnum, Vden);
3981 copy (tN, Q->Nvolume);
3982 copy (tD, Q->Dvolume);
3983 linrat (tN, tD, ONE, Vnum, Vden, ONE, Q->Nvolume, Q->Dvolume);
3984 if (Q->debug)
3985 {
3986 prat ("\n*Volume=", Q->Nvolume, Q->Dvolume);
3987 pmp (" Vnum=", Vnum);
3988 pmp (" Vden=", Vden);
3989 }
3990 lrs_clear_mp(tN); lrs_clear_mp(tD); lrs_clear_mp(Vnum); lrs_clear_mp(Vden);
3991
3992 } /* end of updatevolume */
3993
3994
3995
3996 /***************************************************/
3997 /* Routines for redundancy checking */
3998 /***************************************************/
3999
4000 long
checkredund(lrs_dic * P,lrs_dat * Q)4001 checkredund (lrs_dic * P, lrs_dat * Q)
4002 /* Solve primal feasible lp by least subscript and lex min basis method */
4003 /* to check redundancy of a row in objective function */
4004 /* returns TRUE if redundant, else FALSE */
4005 {
4006 lrs_mp Ns, Nt;
4007 long i, j;
4008 long r, s;
4009
4010 /* assign local variables to structures */
4011 lrs_mp_matrix A = P->A;
4012 long *Row, *Col;
4013 long d = P->d;
4014
4015 lrs_alloc_mp(Ns); lrs_alloc_mp(Nt);
4016 Row = P->Row;
4017 Col = P->Col;
4018
4019 while (selectpivot (P, Q, &i, &j))
4020 {
4021 Q->count[2]++;
4022
4023 /* sign of new value of A[0][0] */
4024 /* is A[0][s]*A[r][0]-A[0][0]*A[r][s] */
4025
4026 r = Row[i];
4027 s = Col[j];
4028
4029 mulint (A[0][s], A[r][0], Ns);
4030 mulint (A[0][0], A[r][s], Nt);
4031
4032 if (mp_greater (Ns, Nt))
4033 {
4034 lrs_clear_mp(Ns); lrs_clear_mp(Nt);
4035 return FALSE; /* non-redundant */
4036 }
4037
4038 pivot (P, Q, i, j);
4039 update (P, Q, &i, &j); /*Update B,C,i,j */
4040
4041 }
4042
4043 lrs_clear_mp(Ns); lrs_clear_mp(Nt);
4044
4045 return !(j < d && i == 0); /* unbounded is also non-redundant */
4046
4047
4048 } /* end of checkredund */
4049
4050 long
checkcobasic(lrs_dic * P,lrs_dat * Q,long index)4051 checkcobasic (lrs_dic * P, lrs_dat * Q, long index)
4052 /* TRUE if index is cobasic and nonredundant */
4053 /* FALSE if basic, or degen. cobasic, where it will get pivoted out */
4054
4055 {
4056
4057 /* assign local variables to structures */
4058
4059 lrs_mp_matrix A = P->A;
4060 long *B, *C, *Row, *Col;
4061 long d = P->d;
4062 long m = P->m;
4063 long debug = Q->debug;
4064 long i = 0;
4065 long j = 0;
4066 long s;
4067
4068 B = P->B;
4069 C = P->C;
4070 Row = P->Row;
4071 Col = P->Col;
4072
4073
4074 while ((j < d) && C[j] != index)
4075 j++;
4076
4077 if (j == d)
4078 return FALSE; /* not cobasic index */
4079
4080
4081 /* index is cobasic */
4082
4083 if (debug)
4084 fprintf (lrs_ofp, "\nindex=%ld cobasic", index);
4085
4086 s = Col[j];
4087 i = Q->lastdv + 1;
4088
4089 while ((i <= m) &&
4090 (zero (A[Row[i]][s]) || !zero (A[Row[i]][0])))
4091 i++;
4092
4093 if (i > m)
4094 {
4095 if (debug)
4096 fprintf (lrs_ofp, " is non-redundant");
4097 return TRUE;
4098 }
4099 if (debug)
4100 fprintf (lrs_ofp, " is degenerate B[i]=%ld", B[i]);
4101
4102 pivot (P, Q, i, j);
4103 update (P, Q, &i, &j); /*Update B,C,i,j */
4104
4105 return FALSE; /*index is no longer cobasic */
4106
4107 } /* end of checkcobasic */
4108
4109 long
checkindex(lrs_dic * P,lrs_dat * Q,long index)4110 checkindex (lrs_dic * P, lrs_dat * Q, long index)
4111 /* 0 if index is non-redundant inequality */
4112 /* 1 if index is redundant inequality */
4113 /* 2 if index is input linearity */
4114 /*NOTE: row is returned all zero if redundant!! */
4115 {
4116 long i, j;
4117
4118 lrs_mp_matrix A = P->A;
4119 long *Row = P->Row;
4120 long *B = P->B;
4121 long d = P->d;
4122 long m = P->m;
4123 long zeroonly=0;
4124
4125 if(index < 0) /* used to zero out known redundant rows in mplrs verifyredund */
4126 {
4127 zeroonly=1;
4128 index=-index;
4129 }
4130
4131 if (Q->debug)
4132 printA (P, Q);
4133
4134 /* each slack index must be checked for redundancy */
4135 /* if in cobasis, it is pivoted out if degenerate */
4136 /* else it is non-redundant */
4137
4138 if (checkcobasic (P, Q, index))
4139 {
4140 return ZERO;
4141 }
4142 /* index is basic */
4143 j = 1;
4144 while ((j <= m) && (B[j] != index))
4145 j++;
4146
4147 i = Row[j];
4148
4149 /* copy row i to cost row, and set it to zero */
4150
4151 for (j = 0; j <= d; j++)
4152 {
4153 copy (A[0][j], A[i][j]);
4154 changesign (A[0][j]);
4155 itomp (ZERO, A[i][j]);
4156 }
4157
4158
4159 if (zeroonly || checkredund (P, Q))
4160 return ONE;
4161
4162 /* non-redundant, copy back and change sign */
4163
4164 for (j = 0; j <= d; j++)
4165 {
4166 copy (A[i][j], A[0][j]);
4167 changesign (A[i][j]);
4168 }
4169
4170 return ZERO;
4171
4172 } /* end of checkindex */
4173
4174
4175 /***************************************************************/
4176 /* */
4177 /* f I/O routines */
4178 /* */
4179 /***************************************************************/
4180
4181 void
lprat(const char * name,long Nt,long Dt)4182 lprat (const char *name, long Nt, long Dt)
4183 /*print the long precision rational Nt/Dt without reducing */
4184 {
4185 if ( Nt > 0 )
4186 fprintf (lrs_ofp, " ");
4187 fprintf (lrs_ofp, "%s%ld", name, Nt);
4188 if (Dt != 1)
4189 fprintf (lrs_ofp, "/%ld", Dt);
4190 fprintf (lrs_ofp, " ");
4191 } /* lprat */
4192
4193 long
lreadrat(long * Num,long * Den)4194 lreadrat (long *Num, long *Den)
4195 /* read a rational string and convert to long */
4196 /* returns true if denominator is not one */
4197 {
4198 char in[MAXINPUT], num[MAXINPUT], den[MAXINPUT];
4199 if(fscanf (lrs_ifp, "%s", in) == EOF)
4200 return(FALSE);
4201 atoaa (in, num, den); /*convert rational to num/dem strings */
4202 *Num = atol (num);
4203 if (den[0] == '\0')
4204 {
4205 *Den = 1L;
4206 return (FALSE);
4207 }
4208 *Den = atol (den);
4209 return (TRUE);
4210 }
4211
4212 void
lrs_getinput(lrs_dic * P,lrs_dat * Q,long * num,long * den,long m,long d)4213 lrs_getinput(lrs_dic *P,lrs_dat *Q,long *num,long *den, long m, long d)
4214 /* code for reading data matrix in lrs/cdd format */
4215 {
4216 long j,row;
4217
4218 printf("\nEnter each row: b_i a_ij j=1..%ld",d);
4219 for (row=1;row<=m;row++)
4220 {
4221 printf("\nEnter row %ld: ",row );
4222 for(j=0;j<=d;j++)
4223 {
4224 lreadrat(&num[j],&den[j]);
4225 lprat(" ",num[j],den[j]);
4226 }
4227
4228 lrs_set_row(P,Q,row,num,den,GE);
4229 }
4230
4231 printf("\nEnter objective row c_j j=1..%ld: ",d);
4232 num[0]=0; den[0]=1;
4233 for(j=1;j<=d;j++)
4234 {
4235 lreadrat(&num[j],&den[j]);
4236 lprat(" ",num[j],den[j]);
4237 }
4238
4239 lrs_set_obj(P,Q,num,den,MAXIMIZE);
4240 }
4241
4242
4243 long
readlinearity(lrs_dat * Q)4244 readlinearity (lrs_dat * Q) /* read in and check linearity list */
4245 {
4246 long i, j;
4247 long nlinearity;
4248 if(fscanf (lrs_ifp, "%ld", &nlinearity)==EOF )
4249 {
4250 lrs_warning(Q,"warning","\nLinearity option invalid, no indices ");
4251 return (FALSE);
4252 }
4253 if (nlinearity < 1)
4254 {
4255 lrs_warning(Q,"warning","\nLinearity option invalid, indices must be positive");
4256 return (FALSE);
4257 }
4258
4259 Q->linearity = (long int*) CALLOC ((nlinearity + 1), sizeof (long));
4260
4261 for (i = 0; i < nlinearity; i++)
4262 {
4263 if(fscanf (lrs_ifp, "%ld", &j)==EOF)
4264 {
4265 lrs_warning(Q,"warning","\nLinearity option invalid, missing indices");
4266 return (FALSE);
4267 }
4268 Q->linearity[i] = j;
4269
4270 }
4271 for (i = 1; i < nlinearity; i++) /*sort in order */
4272 reorder (Q->linearity, nlinearity);
4273
4274 Q->nlinearity = nlinearity;
4275 Q->polytope = FALSE;
4276 return TRUE;
4277 } /* end readlinearity */
4278
4279 long
readredund(lrs_dat * Q)4280 readredund (lrs_dat * Q) /* read in and check linearity list */
4281 {
4282 long i,j,k;
4283 char *mess;
4284 int len=0;
4285
4286 if(fscanf (lrs_ifp, "%ld", &k)==EOF )
4287 {
4288 lrs_warning(Q,"warning","\nredund_list option invalid: no indices ");
4289 return (FALSE);
4290 }
4291 if ( k < 0 )
4292 {
4293 lrs_warning(Q,"warning","\nredund_list option invalid, first index must be >= 0");
4294 return (FALSE);
4295 }
4296
4297 for (i = 1; i <= Q->m; i++) /*reset any previous redund option except =2 values */
4298 if (Q->redineq[i] != 2)
4299 Q->redineq[i]=0;
4300 Q->redineq[0]=1;
4301
4302 for (i = 0; i < k; i++)
4303 {
4304 if(fscanf (lrs_ifp, "%ld", &j)==EOF)
4305 {
4306 lrs_warning(Q,"warning","\nredund_list option invalid: missing indices");
4307 fflush(lrs_ofp);
4308 return (FALSE);
4309 }
4310
4311 if( j< 0 || j > Q->m)
4312 {
4313 fprintf (lrs_ofp,"\nredund_list option invalid: indices not between 1 and %ld",Q->m);
4314 return (FALSE);
4315 }
4316 Q->redineq[j] = 1;
4317
4318 }
4319
4320 //if(!Q->mplrs && overflow != 2 )
4321 if( overflow != 2 )
4322 {
4323 mess=(char *)malloc(20*Q->m*sizeof(char));
4324 len=sprintf(mess,"redund_list %ld ",k);
4325 for (i=1;i<=Q->m;i++)
4326 if(Q->redineq[i] == 1)
4327 len=len+sprintf(mess+len," %ld",i);
4328 lrs_warning(Q,"warning",mess);
4329 free(mess);
4330 }
4331 return TRUE;
4332 } /* end readredund */
4333
4334
4335
4336 long
readfacets(lrs_dat * Q,long facet[])4337 readfacets (lrs_dat * Q, long facet[])
4338 /* read and check facet list for obvious errors during start/restart */
4339 /* this must be done after linearity option is processed!! */
4340 {
4341 long i, j;
4342 char str[1000000],*p,*e;
4343
4344 /* assign local variables to structures */
4345 long m, d;
4346 long *linearity = Q->linearity;
4347 m = Q->m;
4348 d = Q->inputd;
4349
4350 /* modified 2018.6.7 to fix bug restarting with less than full dimension input */
4351 /* number of restart indices is not known at this point */
4352
4353 j=Q->nlinearity; /* note we place these after the linearity indices */
4354
4355 fgets(str,1000000,lrs_ifp); /* pick up indices from the input line */
4356 for (p = str; ; p = e) {
4357 facet[j] = strtol(p, &e, 10);
4358 if (p == e)
4359 break;
4360 if(!Q->mplrs && Q->verbose && overflow != 2)
4361 fprintf(lrs_ofp," %ld",facet[j] );
4362
4363 //fprintf(lrs_ofp,"\n j %ld d %ld facet %ld",j,d,facet[j]);
4364
4365 /* 2010.4.26 nonnegative option needs larger range of indices */
4366 if(Q->nonnegative)
4367 if (facet[j] < 1 || facet[j] > m+d)
4368 {
4369 fprintf (lrs_ofp,"\n Start/Restart cobasic indices must be in range 1 .. %ld ", m+d);
4370 return FALSE;
4371 }
4372
4373 if(!Q->nonnegative)
4374 if (facet[j] < 1 || facet[j] > m)
4375 {
4376 fprintf (lrs_ofp, "\n Start/Restart cobasic indices must be in range 1 .. %ld ", m);
4377 return FALSE;
4378 }
4379
4380 for (i = 0; i < Q->nlinearity; i++)
4381 if (linearity[i] == facet[j])
4382 {
4383 fprintf (lrs_ofp, "\n Start/Restart cobasic indices should not include linearities");
4384 return FALSE;
4385 }
4386 /* bug fix 2011.8.1 reported by Steven Wu*/
4387 for (i = Q->nlinearity; i < j; i++)
4388 /* end bug fix 2011.8.1 */
4389
4390 if (facet[i] == facet[j])
4391 {
4392 fprintf (lrs_ofp, "\n Start/Restart cobasic indices must be distinct");
4393 return FALSE;
4394 }
4395 j++;
4396 }
4397 return TRUE;
4398 } /* end of readfacets */
4399
4400 long
readremain(lrs_dat * Q)4401 readremain (lrs_dat * Q) /* read in and check ordered list of vars for extract */
4402 {
4403 long i, j, kk;
4404 long nremain;
4405 long k=0;
4406 long n=Q->n;
4407 long *remain;
4408
4409 Q->remain = (long int*) CALLOC ((n + 2), sizeof (long));
4410
4411 remain=Q->remain;
4412 for (i=0;i<n-1;i++)
4413 remain[i]=0;
4414
4415 if(fscanf (lrs_ifp, "%ld", &nremain)==EOF )
4416 {
4417 for (i=1;i<n;i++)
4418 remain[i-1]=i;
4419 return 0;
4420 }
4421
4422 if( nremain > n-1)
4423 {
4424 nremain = n-1;
4425 fprintf (lrs_ofp, "\n*extract: too many indices, first %ld taken",n-1);
4426 }
4427
4428 for (i = 0; i < nremain; i++)
4429 {
4430 if(fscanf (lrs_ifp, "%ld", &j)==EOF)
4431 {
4432 fprintf (lrs_ofp, "\n*extract: missing indices");
4433 break;
4434 }
4435 if(j>0 && j<n)
4436 {
4437 kk=0;
4438 while (kk<k && remain[kk] != j)
4439 kk++;
4440 if(kk==k)
4441 remain[k++] = j;
4442 else
4443 fprintf (lrs_ofp, "\n*extract: duplicate index %ld skipped",j);
4444 }
4445 else
4446 fprintf (lrs_ofp, "\n*extract: index %ld out of range 1 to %ld",j,n-1);
4447 }
4448 nremain=0;
4449 while (nremain < n && remain[nremain] != 0 )
4450 nremain++;
4451
4452 /* if nlinearity>0 fill up list with remaining decision variables */
4453
4454 if(Q->nlinearity > 0)
4455 for (i=1;i<n;i++)
4456 {
4457 j=0;
4458 while (j<nremain && remain[j] != i)
4459 j++;
4460 if (j == nremain)
4461 remain[nremain++]=i;
4462 }
4463 return 0;
4464 } /* readremain */
4465
4466
extractcols(lrs_dic * P,lrs_dat * Q)4467 long extractcols (lrs_dic * P, lrs_dat * Q)
4468 {
4469 /* 2020.2.2 */
4470 /* when there are no linearities extract just prints the chosen columns verbatim */
4471 long i,j,m,n;
4472 long ncols=0;
4473 lrs_mp_matrix A;
4474 long *Col, *Row, *remain, *output;
4475
4476 A = P->A;
4477 Col = P->Col;
4478 Row = P->Row;
4479 remain=Q->remain;
4480 output=Q->temparray;
4481 m=P->m;
4482 n=Q->n;
4483
4484 for(j=0;j<n;j++)
4485 {
4486 output[j]=0;
4487 if (remain[j])
4488 ncols++;
4489 }
4490
4491 for(j=0;j<n;j++)
4492 output[remain[j]]=1;
4493
4494 fprintf(lrs_ofp,"\n*output");
4495 for (j=0;j<n;j++)
4496 fprintf(lrs_ofp," %ld",output[j]);
4497 fprintf(lrs_ofp,"\n*columns retained:");
4498 for (j=0;j<n;j++)
4499 if(output[j])
4500 fprintf(lrs_ofp," %ld",j);
4501
4502 if(Q->hull)
4503 fprintf(lrs_ofp,"\nV-representation\nbegin");
4504 else
4505 fprintf(lrs_ofp,"\nH-representation\nbegin");
4506 fprintf(lrs_ofp,"\n%ld %ld rational",m,ncols+1);
4507
4508 for(i=1;i<=m;i++ )
4509 {
4510 reducearray(A[Row[i]],n);
4511 fprintf(lrs_ofp,"\n");
4512 if(Q->hull)
4513 {
4514 for(j=0;j<n;j++)
4515 if(output[j])
4516 pmp("",A[Row[i]][Col[j]]);
4517 }
4518 else
4519 {
4520 pmp("",A[Row[i]][0]);
4521 for(j=1;j<n;j++)
4522 if(output[j])
4523 pmp("",A[Row[i]][Col[j-1]]);
4524 }
4525 }
4526 fprintf(lrs_ofp,"\nend");
4527
4528 fprintf(lrs_ofp,"\n*columns retained:");
4529 for (j=0;j<n;j++)
4530 if(output[j])
4531 fprintf(lrs_ofp," %ld",j);
4532
4533 fprintf(lrs_ofp,"\n");
4534
4535 if(Q->debug)
4536 printA(P,Q);
4537
4538 return 0;
4539 } /* extractcols */
4540
linextractcols(lrs_dic * P,lrs_dat * Q)4541 long linextractcols (lrs_dic * P, lrs_dat * Q)
4542 /* 2020.2.2 */
4543 /* extract option to output the reduced A matrix after linearities are removed */
4544 /* should be followed by redund to get minimum representation */
4545 {
4546 long d,i,j,k,m,n;
4547 long ii,jj;
4548 long nlinearity=Q->nlinearity;
4549 lrs_mp_matrix A;
4550 long *B, *C, *Col, *Row, *remain;
4551
4552 A = P->A;
4553 B = P->B;
4554 C = P->C;
4555 Col = P->Col;
4556 Row = P->Row;
4557 remain=Q->remain;
4558
4559 m=P->m;
4560 n=Q->n;
4561 d=Q->inputd;
4562
4563 fprintf(lrs_ofp,"\n*extract col order: ");
4564
4565 for(j=0;j<n-1;j++)
4566 fprintf (lrs_ofp, "%ld ",remain[j]);
4567
4568 for(k=0;k<n-1;k++) /* go through input order for vars to remain */
4569 {
4570 i=1;
4571 while( i<= m)
4572 {
4573 if(B[i]== remain[k])
4574 {
4575 j=0;
4576 while(j+nlinearity<d && (C[j]<= d || zero(A[Row[i]][Col[j]])))
4577 j++;
4578 if (j+nlinearity<d)
4579 {
4580 ii=i;
4581 jj=j;
4582 pivot(P,Q,ii,jj);
4583 update(P,Q,&ii,&jj);
4584 i=0;
4585 }
4586 } /* if B[i]... */
4587 i++;
4588 } /* while i <= m */
4589 } /* for k=0 */
4590
4591 if(Q->hull)
4592 fprintf(lrs_ofp,"\n*columns retained:");
4593 else
4594 fprintf(lrs_ofp,"\n*columns retained: 0");
4595 for (j=0;j<d-nlinearity;j++)
4596 fprintf(lrs_ofp," %ld",C[j]-Q->hull);
4597
4598 if(Q->hull)
4599 fprintf(lrs_ofp,"\nV-representation\nbegin");
4600 else
4601 fprintf(lrs_ofp,"\nH-representation\nbegin");
4602 fprintf(lrs_ofp,"\n%ld %ld rational",m-nlinearity,1+P->d-Q->hull);
4603
4604 for(i=nlinearity+1;i<=m;i++ )
4605 {
4606 reducearray(A[Row[i]],n-nlinearity);
4607 fprintf(lrs_ofp,"\n");
4608 if(!Q->hull)
4609 pmp("",A[Row[i]][0]);
4610 for(j=0;j<d-nlinearity;j++)
4611 pmp("",A[Row[i]][Col[j]]);
4612 }
4613 fprintf(lrs_ofp,"\nend");
4614 if(Q->hull)
4615 fprintf(lrs_ofp,"\n*columns retained:");
4616 else
4617 fprintf(lrs_ofp,"\n*columns retained: 0");
4618 for (j=0;j<d-nlinearity;j++)
4619 fprintf(lrs_ofp," %ld",C[j]-Q->hull);
4620 fprintf(lrs_ofp,"\n");
4621
4622 if(Q->debug)
4623 printA(P,Q);
4624
4625 return 0;
4626 } /* linextractcols */
4627
4628
4629 void
printA(lrs_dic * P,lrs_dat * Q)4630 printA (lrs_dic * P, lrs_dat * Q) /* print the integer m by n array A
4631 with B,C,Row,Col vectors */
4632 {
4633 long i, j;
4634 /* assign local variables to structures */
4635 lrs_mp_matrix A = P->A;
4636 long *B = P->B;
4637 long *C = P->C;
4638 long *Row = P->Row;
4639 long *Col = P->Col;
4640 long m, d;
4641 m = P->m;
4642 d = P->d;
4643
4644 fprintf (lrs_ofp, "\n Basis ");
4645 for (i = 0; i <= m; i++)
4646 fprintf (lrs_ofp, "%ld ", B[i]);
4647 fprintf (lrs_ofp, " Row ");
4648 for (i = 0; i <= m; i++)
4649 fprintf (lrs_ofp, "%ld ", Row[i]);
4650 fprintf (lrs_ofp, "\n Co-Basis ");
4651 for (i = 0; i <= d; i++)
4652 fprintf (lrs_ofp, "%ld ", C[i]);
4653 fprintf (lrs_ofp, " Column ");
4654 for (i = 0; i <= d; i++)
4655 fprintf (lrs_ofp, "%ld ", Col[i]);
4656 pmp (" det=", P->det);
4657 fprintf (lrs_ofp, "\n");
4658 i=0;
4659 while ( i<= m )
4660 {
4661 for (j = 0; j <= d; j++)
4662 pimat (P, i, j, A[Row[i]][Col[j]], "A");
4663 fprintf (lrs_ofp, "\n");
4664 if (i==0 && Q->nonnegative) /* skip basic rows - don't exist! */
4665 i=d;
4666 i++;
4667 fflush (stdout);
4668 }
4669 fflush (stdout);
4670 }
4671
4672
4673 void
pimat(lrs_dic * P,long r,long s,lrs_mp Nt,const char * name)4674 pimat (lrs_dic * P, long r, long s, lrs_mp Nt, const char *name)
4675 /*print the long precision integer in row r col s of matrix A */
4676 {
4677 long *B = P->B;
4678 long *C = P->C;
4679 if (s == 0)
4680 fprintf (lrs_ofp, "%s[%ld][%ld]=", name, B[r], C[s]);
4681 else
4682 fprintf (lrs_ofp, "[%ld]=", C[s]);
4683 pmp ("", Nt);
4684
4685 }
4686
4687 /***************************************************************/
4688 /* */
4689 /* Routines for caching, allocating etc. */
4690 /* */
4691 /***************************************************************/
4692
4693 /* From here mostly Bremner's handiwork */
4694
4695 static void
cache_dict(lrs_dic ** D_p,lrs_dat * global,long i,long j)4696 cache_dict (lrs_dic ** D_p, lrs_dat * global, long i, long j)
4697 {
4698
4699 if (dict_limit > 1)
4700 {
4701 /* save row, column indicies */
4702 (*D_p)->i = i;
4703 (*D_p)->j = j;
4704
4705 /* Make a new, blank spot at the end of the queue to copy into */
4706
4707
4708 pushQ (global, (*D_p)->m, (*D_p)->d, (*D_p)->m_A);
4709
4710 /*2019.6.7 This ought not to happen but it does */
4711
4712 if( global->Qtail== *D_p)
4713 return;
4714
4715 copy_dict (global, global->Qtail, *D_p); /* Copy current dictionary */
4716 }
4717 *D_p = global->Qtail;
4718
4719 }
4720
4721 void
copy_dict(lrs_dat * global,lrs_dic * dest,lrs_dic * src)4722 copy_dict (lrs_dat * global, lrs_dic * dest, lrs_dic * src)
4723 {
4724 long m = src->m;
4725 long m_A = src->m_A; /* number of rows in A */
4726 long d = src->d;
4727 long r,s;
4728
4729 if( dest == src)
4730 {
4731 if(global->mplrs)
4732 lrs_post_output("warning", "*copy_dict has dest=src -ignoring copy");
4733 else
4734 fprintf(stderr,"*copy_dict has dest=src -ignoring copy");
4735 return;
4736 }
4737
4738 #if defined(GMP) || defined(FLINT)
4739
4740 for ( r=0;r<=m_A;r++)
4741 for( s=0;s<=d;s++)
4742 copy(dest->A[r][s],src->A[r][s]);
4743
4744 #else /* fast copy for MP and LRSLONG arithmetic */
4745 /* Note that the "A" pointer trees need not be copied, since they
4746 always point to the same places within the corresponding space
4747 */
4748 /* I wish I understood the above remark. For the time being, do it the easy way for Nash */
4749 /* Looking at lrs_alloc_mp_matrix for MP and LRSLONG, A[0][0] is the
4750 * start of araw, which holds the actual values and so the memcpy below
4751 * copies the values. The pointer trees (A, A[i], A[i][j]) already point
4752 * to the appropriate places: we don't want to change the pointers, only the
4753 * values. lrs_alloc_mp_matrix is different for GMP.
4754 */
4755 if(global->nash)
4756 {
4757 for ( r=0;r<=m_A;r++)
4758 for( s=0;s<=d;s++)
4759 copy(dest->A[r][s],src->A[r][s]);
4760 }
4761 else
4762 #ifdef B128
4763 memcpy (dest->A[0][0], (global->Qtail->prev)->A[0][0],
4764 (d + 1) * (lrs_digits + 1) * (m_A + 1) * sizeof (__int128));
4765 #else
4766 memcpy (dest->A[0][0], (global->Qtail->prev)->A[0][0],
4767 (d + 1) * (lrs_digits + 1) * (m_A + 1) * sizeof (long long));
4768 #endif
4769
4770 #endif
4771
4772 dest->i = src->i;
4773 dest->j = src->j;
4774 dest->m = m;
4775 dest->d = d;
4776 dest->d_orig = src->d_orig;
4777 dest->m_A = src->m_A;
4778
4779 dest->depth = src->depth;
4780 dest->lexflag = src->lexflag;
4781
4782 copy (dest->det, src->det);
4783 copy (dest->objnum, src->objnum);
4784 copy (dest->objden, src->objden);
4785
4786 if (global->debug)
4787 fprintf (lrs_ofp, "\nSaving dict at depth %ld\n", src->depth);
4788
4789 memcpy (dest->B, src->B, (m + 1) * sizeof (long));
4790 memcpy (dest->C, src->C, (d + 1) * sizeof (long));
4791 memcpy (dest->Row, src->Row, (m + 1) * sizeof (long));
4792 memcpy (dest->Col, src->Col, (d + 1) * sizeof (long));
4793 }
4794
4795 /*
4796 * pushQ(lrs_dat *globals,m,d):
4797 * this routine ensures that Qtail points to a record that
4798 * may be copied into.
4799 *
4800 * It may create a new record, or it may just move the head pointer
4801 * forward so that know that the old record has been overwritten.
4802 */
4803 #if 0
4804 #define TRACE(s) fprintf(stderr,"\n%s %p %p\n",s,global->Qhead,global->Qtail);
4805 #else
4806 #define TRACE(s)
4807 #endif
4808
4809 static void
pushQ(lrs_dat * global,long m,long d,long m_A)4810 pushQ (lrs_dat * global, long m, long d ,long m_A)
4811 {
4812
4813 if ((global->Qtail->next) == global->Qhead)
4814 {
4815 /* the Queue is full */
4816 if (dict_count < dict_limit)
4817 {
4818 /* but we are allowed to create more */
4819 lrs_dic *p;
4820
4821 p = new_lrs_dic (m, d, m_A);
4822
4823 if (p)
4824 {
4825
4826 /* we successfully created another record */
4827
4828 p->next = global->Qtail->next;
4829 (global->Qtail->next)->prev = p;
4830 (global->Qtail->next) = p;
4831 p->prev = global->Qtail;
4832
4833 dict_count++;
4834 global->Qtail = p;
4835
4836 TRACE ("Added new record to Q");
4837
4838 }
4839 else
4840 {
4841 /* virtual memory exhausted. bummer */
4842 global->Qhead = global->Qhead->next;
4843 global->Qtail = global->Qtail->next;
4844
4845 TRACE ("VM exhausted");
4846 }
4847 }
4848 else
4849 {
4850 /*
4851 * user defined limit reached. start overwriting the
4852 * beginning of Q
4853 */
4854 global->Qhead = global->Qhead->next;
4855 global->Qtail = global->Qtail->next;
4856 TRACE ("User limit");
4857
4858 }
4859 }
4860
4861 else
4862 {
4863 global->Qtail = global->Qtail->next;
4864 TRACE ("Reusing");
4865 }
4866 }
4867
4868 lrs_dic *
lrs_getdic(lrs_dat * Q)4869 lrs_getdic(lrs_dat *Q)
4870 /* create another dictionary for Q without copying any values */
4871 /* derived from lrs_alloc_dic, used by nash.c */
4872 {
4873 lrs_dic *p;
4874
4875 long m;
4876
4877 m = Q->m;
4878
4879 /* nonnegative flag set means that problem is d rows "bigger" */
4880 /* since nonnegative constraints are not kept explicitly */
4881
4882
4883 if(Q->nonnegative)
4884 m = m+Q->inputd;
4885
4886 p = new_lrs_dic (m, Q->inputd, Q->m);
4887 if (!p)
4888 return NULL;
4889
4890 p->next = p;
4891 p->prev = p;
4892 Q->Qhead = p;
4893 Q->Qtail = p;
4894
4895 return p;
4896 }
4897
4898
4899 #define NULLRETURN(e) if (!(e)) return NULL;
4900
4901 static lrs_dic *
new_lrs_dic(long m,long d,long m_A)4902 new_lrs_dic (long m, long d, long m_A)
4903 {
4904 lrs_dic *p;
4905
4906 NULLRETURN (p = (lrs_dic *) malloc (sizeof (lrs_dic)));
4907
4908
4909 NULLRETURN (p->B = (long int*) calloc ((m + 1), sizeof (long)));
4910 NULLRETURN (p->Row = (long int*) calloc ((m + 1), sizeof (long)));
4911
4912 NULLRETURN (p->C = (long int*) calloc ((d + 1), sizeof (long)));
4913 NULLRETURN (p->Col = (long int*) calloc ((d + 1), sizeof (long)));
4914
4915 #if defined(GMP) || defined(FLINT)
4916 lrs_alloc_mp(p->det);
4917 lrs_alloc_mp(p->objnum);
4918 lrs_alloc_mp(p->objden);
4919 #endif
4920
4921 p->d_orig=d;
4922 p->A=lrs_alloc_mp_matrix(m_A,d);
4923
4924
4925 return p;
4926 }
4927
4928 void
lrs_free_dic(lrs_dic * P,lrs_dat * Q)4929 lrs_free_dic (lrs_dic * P, lrs_dat *Q)
4930 {
4931 /* do the same steps as for allocation, but backwards */
4932 /* gmp variables cannot be cleared using free: use lrs_clear_mp* */
4933 lrs_dic *P1;
4934
4935 if (Q == NULL )
4936 {
4937 if(Q->mplrs)
4938 lrs_post_output("warning","*lrs_free_dic trying to free null Q : skipped");
4939 else
4940 fprintf(stderr,"*lrs_free_dic trying to free null Q : skipped");
4941 return;
4942 }
4943
4944 if (P == NULL )
4945 {
4946 if(Q->mplrs)
4947 lrs_post_output("warning","*lrs_free_dic trying to free null P : skipped");
4948 else
4949 fprintf(stderr,"*lrs_free_dic trying to free null P : skipped");
4950 return;
4951 }
4952 /* repeat until cache is empty */
4953
4954 do
4955 {
4956 /* I moved these here because I'm not certain the cached dictionaries
4957 need to be the same size. Well, it doesn't cost anything to be safe. db */
4958
4959 long d = P->d_orig;
4960 long m_A = P->m_A;
4961
4962 lrs_clear_mp_matrix (P->A,m_A,d);
4963
4964 /* "it is a ghastly error to free something not assigned my malloc" KR167 */
4965 /* so don't try: free (P->det); */
4966
4967 lrs_clear_mp (P->det);
4968 lrs_clear_mp (P->objnum);
4969 lrs_clear_mp (P->objden);
4970
4971 free (P->Row);
4972 free (P->Col);
4973 free (P->C);
4974 free (P->B);
4975
4976 /* go to next record in cache if any */
4977 P1 =P->next;
4978 free (P);
4979 P=P1;
4980
4981 } while (Q->Qhead != P );
4982
4983
4984 }
4985
4986 void
lrs_free_dic2(lrs_dic * P,lrs_dat * Q)4987 lrs_free_dic2 (lrs_dic * P, lrs_dat *Q)
4988 {
4989 /* do the same steps as for allocation, but backwards */
4990 /* same as lrs_free_dic except no cache for P */
4991 /* I moved these here because I'm not certain the cached dictionaries
4992 need to be the same size. Well, it doesn't cost anything to be safe. db */
4993
4994 long d = P->d_orig;
4995 long m_A = P->m_A;
4996
4997
4998 lrs_clear_mp_matrix (P->A,m_A,d);
4999
5000 /* "it is a ghastly error to free something not assigned my malloc" KR167 */
5001 /* so don't try: free (P->det); */
5002
5003 lrs_clear_mp (P->det);
5004 lrs_clear_mp (P->objnum);
5005 lrs_clear_mp (P->objden);
5006
5007 free (P->Row);
5008 free (P->Col);
5009 free (P->C);
5010 free (P->B);
5011
5012 free (P);
5013
5014 }
5015
5016 void
lrs_free_dat(lrs_dat * Q)5017 lrs_free_dat ( lrs_dat *Q )
5018 {
5019
5020 if (Q == NULL)
5021 {
5022 if(Q->mplrs)
5023 lrs_post_output("warning","*lrs_free_dat trying tor free null Q : skipped");
5024 else
5025 fprintf(stderr,"*lrs_free_dat trying tor free null Q : skipped");
5026 return;
5027 }
5028
5029 if(Q->extract)
5030 free(Q->remain);
5031
5032 /* most of these items were allocated in lrs_alloc_dic */
5033
5034 lrs_clear_mp_vector (Q->Gcd,Q->m);
5035 lrs_clear_mp_vector (Q->Lcm,Q->m);
5036 lrs_clear_mp_vector (Q->output,Q->n);
5037
5038 lrs_clear_mp (Q->sumdet);
5039 lrs_clear_mp (Q->Nvolume);
5040 lrs_clear_mp (Q->Dvolume);
5041 lrs_clear_mp (Q->saved_det);
5042 lrs_clear_mp (Q->boundd);
5043 lrs_clear_mp (Q->boundn);
5044
5045 free (Q->facet);
5046 free (Q->redundcol);
5047 free (Q->inequality);
5048 free (Q->linearity);
5049 free (Q->redineq);
5050 free (Q->minratio);
5051 free (Q->temparray);
5052 free (Q->startcob);
5053
5054 free (Q->name);
5055 free (Q->saved_C);
5056
5057 lrs_global_count--;
5058
5059 free(Q);
5060 }
5061
5062
5063 static long
check_cache(lrs_dic ** D_p,lrs_dat * global,long * i_p,long * j_p)5064 check_cache (lrs_dic ** D_p, lrs_dat * global, long *i_p, long *j_p)
5065 {
5066 /* assign local variables to structures */
5067
5068
5069
5070 cache_tries++;
5071
5072 if (global->Qtail == global->Qhead)
5073 {
5074 TRACE ("cache miss");
5075 /* Q has only one element */
5076 cache_misses++;
5077 return 0;
5078
5079 }
5080 else
5081 {
5082 global->Qtail = global->Qtail->prev;
5083
5084 *D_p = global->Qtail;
5085
5086 *i_p = global->Qtail->i;
5087 *j_p = global->Qtail->j;
5088
5089 TRACE ("restoring dict");
5090 return 1;
5091 }
5092 }
5093
5094
5095 lrs_dic *
lrs_alloc_dic(lrs_dat * Q)5096 lrs_alloc_dic (lrs_dat * Q)
5097 /* allocate and initialize lrs_dic */
5098 {
5099
5100 lrs_dic *p;
5101 long i, j;
5102 long m, d, m_A;
5103
5104 if (Q->hull) /* d=col dimension of A */
5105 Q->inputd = Q->n; /* extra column for hull */
5106 else
5107 Q->inputd = Q->n - 1;
5108
5109 m = Q->m;
5110 d = Q->inputd;
5111 m_A = m; /* number of rows in A */
5112
5113 /* nonnegative flag set means that problem is d rows "bigger" */
5114 /* since nonnegative constraints are not kept explicitly */
5115
5116 if(Q->nonnegative)
5117 m = m+d;
5118
5119 p = new_lrs_dic (m, d, m_A);
5120 if (!p)
5121 return NULL;
5122
5123 p->next = p;
5124 p->prev = p;
5125 Q->Qhead = p;
5126 Q->Qtail = p;
5127
5128
5129 dict_count = 1;
5130 dict_limit = 50;
5131 cache_tries = 0;
5132 cache_misses = 0;
5133
5134 /* Initializations */
5135
5136 p->d = p->d_orig = d;
5137 p->m = m;
5138 p->m_A = m_A;
5139 p->depth = 0L;
5140 p->lexflag = TRUE;
5141 itomp (ONE, p->det);
5142 itomp (ZERO, p->objnum);
5143 itomp (ONE, p->objden);
5144
5145 /*m+d+1 is the number of variables, labelled 0,1,2,...,m+d */
5146 /* initialize array to zero */
5147 for (i = 0; i <= m_A; i++)
5148 for (j = 0; j <= d; j++)
5149 itomp (ZERO, p->A[i][j]);
5150
5151 if (Q->nlinearity == ZERO) /* linearity may already be allocated */
5152 Q->linearity = (long int*) CALLOC ((m + d + 1), sizeof (long));
5153
5154 Q->inequality = (long int*) CALLOC ((m + d + 1), sizeof (long));
5155 Q->facet = (long int*) CALLOC ((unsigned) m + d + 1, sizeof (long));
5156 Q->redundcol = (long int*) CALLOC ((m + d + 1), sizeof (long));
5157 Q->minratio = (long int*) CALLOC ((m+d + 1), sizeof (long));
5158 /* 2011.7.14 minratio[m]=0 for degen =1 for nondegen pivot*/
5159 Q->redineq = (long int*) CALLOC ((m + d + 1), sizeof (long));
5160 Q->temparray = (long int*) CALLOC ((unsigned) m + d + 1, sizeof (long));
5161
5162 Q->inequality[0] = 2L;
5163 Q->Gcd = lrs_alloc_mp_vector(m);
5164 Q->Lcm = lrs_alloc_mp_vector(m);
5165 Q->output = lrs_alloc_mp_vector(Q->n);
5166 Q->saved_C = (long int*) CALLOC (d + 1, sizeof (long));
5167
5168 Q->lastdv = d; /* last decision variable may be decreased */
5169 /* if there are redundant columns */
5170
5171 for (i = 0; i < m+d+1; i++)
5172 {
5173 Q->redineq[i]=1;
5174 Q->inequality[i]=0;
5175 }
5176
5177
5178 /*initialize basis and co-basis indices, and row col locations */
5179 /*if nonnegative, we label differently to avoid initial pivots */
5180 /* set basic indices and rows */
5181 if(Q->nonnegative)
5182 for (i = 0; i <= m; i++)
5183 {
5184 p->B[i] = i;
5185 if (i <= d )
5186 p->Row[i]=0; /* no row for decision variables */
5187 else
5188 p->Row[i]=i-d;
5189 }
5190 else
5191 for (i = 0; i <= m; i++)
5192 {
5193 if (i == 0 )
5194 p->B[0]=0;
5195 else
5196 p->B[i] = d + i;
5197 p->Row[i] = i;
5198 }
5199
5200 for (j = 0; j < d; j++)
5201 {
5202 if(Q->nonnegative)
5203 p->C[j] = m+j+1;
5204 else
5205 p->C[j] = j + 1;
5206 p->Col[j] = j + 1;
5207 }
5208 p->C[d] = m + d + 1;
5209 p->Col[d] = 0;
5210 return p;
5211 } /* end of lrs_alloc_dic */
5212
5213
5214 /*
5215 this routine makes a copy of the information needed to restart,
5216 so that we can guarantee that if a signal is received, we
5217 can guarantee that nobody is messing with it.
5218 This as opposed to adding all kinds of critical regions in
5219 the main line code.
5220
5221 It is also used to make sure that in case of overflow, we
5222 have a valid cobasis to restart from.
5223 */
5224 static void
save_basis(lrs_dic * P,lrs_dat * Q)5225 save_basis (lrs_dic * P, lrs_dat * Q)
5226 {
5227 int i;
5228 /* assign local variables to structures */
5229 long *C = P->C;
5230 long d;
5231
5232 #ifndef SIGNALS
5233 sigset_t oset, blockset;
5234 sigemptyset (&blockset);
5235 sigaddset (&blockset, SIGTERM);
5236 sigaddset (&blockset, SIGHUP);
5237 sigaddset (&blockset, SIGUSR1);
5238
5239 errcheck ("sigprocmask", sigprocmask (SIG_BLOCK, &blockset, &oset));
5240 #endif
5241 d = P->d;
5242
5243 Q->saved_flag = 1;
5244
5245 for (i = 0; i < 5; i++)
5246 Q->saved_count[i] = Q->count[i];
5247
5248 for (i = 0; i < d + 1; i++)
5249 Q->saved_C[i] = C[i];
5250
5251 copy (Q->saved_det, P->det);
5252
5253 Q->saved_d = P->d;
5254 Q->saved_depth = P->depth;
5255
5256 #ifndef SIGNALS
5257 errcheck ("sigprocmask", sigprocmask (SIG_SETMASK, &oset, 0));
5258 #endif
5259 }
5260
5261 /* digits overflow is a call from lrs_mp package */
5262
5263 void
digits_overflow()5264 digits_overflow ()
5265 {
5266 fprintf (lrs_ofp, "\nOverflow at digits=%ld", DIG2DEC (lrs_digits));
5267 fprintf (lrs_ofp, "\nRerun with option: digits n, where n > %ld\n", DIG2DEC (lrs_digits));
5268 lrs_dump_state ();
5269
5270 notimpl("");
5271 }
5272
5273 static void
lrs_dump_state()5274 lrs_dump_state ()
5275 {
5276 long i;
5277
5278 fprintf (lrs_ofp, "\n\nlrs_lib: checkpointing:\n");
5279
5280 #ifdef MP
5281 fprintf (stderr, "lrs_lib: Current digits at %ld out of %ld\n",
5282 DIG2DEC (lrs_record_digits),
5283 DIG2DEC (lrs_digits));
5284 #endif
5285
5286 for (i = 0; i < lrs_global_count; i++)
5287 {
5288 print_basis (lrs_ofp, lrs_global_list[i]);
5289 }
5290 fprintf (lrs_ofp, "lrs_lib: checkpoint finished\n");
5291 }
5292
5293
5294 /* print out the saved copy of the basis */
5295 void
print_basis(FILE * fp,lrs_dat * global)5296 print_basis (FILE * fp, lrs_dat * global)
5297 {
5298 int i;
5299 /* assign local variables to structures */
5300 fprintf (fp, "lrs_lib: State #%ld: (%s)\t", global->id, global->name);
5301
5302 if (global->saved_flag)
5303 {
5304
5305 /* legacy output which is not actually correct for V-representations as V# is not used */
5306 /*
5307 fprintf (fp, "V#%ld R#%ld B#%ld h=%ld facets ",
5308 global->saved_count[1],
5309 global->saved_count[0],
5310 global->saved_count[2],
5311 global->saved_depth);
5312 for (i = 0; i < global->saved_d; i++)
5313 fprintf (fp, "%ld ",
5314 global->inequality[global->saved_C[i] - global->lastdv]);
5315 pmp (" det=", global->saved_det);
5316 fprintf (fp, "\n");
5317 */
5318
5319 if( global->hull)
5320 fprintf (fp, "\nrestart %ld %ld %ld ",
5321 global->saved_count[0],
5322 global->saved_count[2],
5323 global->saved_depth);
5324 else
5325 fprintf (fp, "\nrestart %ld %ld %ld %ld ",
5326 global->saved_count[1],
5327 global->saved_count[0],
5328 global->saved_count[2],
5329 global->saved_depth);
5330
5331 for (i = 0; i < global->saved_d; i++)
5332 fprintf (fp, "%ld ",
5333 global->inequality[global->saved_C[i] - global->lastdv]);
5334 if(global->saved_count[4] >0)
5335 fprintf (fp, "\nintegervertices %ld", global->saved_count[4]);
5336 fprintf (fp, "\n");
5337
5338
5339 }
5340 else
5341 {
5342 fprintf (fp, "lrs_lib: Computing initial basis\n");
5343 }
5344
5345
5346 fflush (fp);
5347 }
5348
5349 #ifndef SIGNALS
5350
5351 /*
5352 If given a signal
5353 USR1 print current cobasis and continue
5354 TERM print current cobasis and terminate
5355 INT (ctrl-C) ditto
5356 HUP ditto
5357 */
5358 static void
setup_signals()5359 setup_signals ()
5360 {
5361 errcheck ("signal", signal (SIGTERM, die_gracefully));
5362 errcheck ("signal", signal (SIGALRM, timecheck));
5363 errcheck ("signal", signal (SIGHUP, die_gracefully));
5364 errcheck ("signal", signal (SIGINT, die_gracefully));
5365 errcheck ("signal", signal (SIGUSR1, checkpoint));
5366 }
5367
5368 static void
timecheck()5369 timecheck ()
5370 {
5371 lrs_dump_state ();
5372 errcheck ("signal", signal (SIGALRM, timecheck));
5373 alarm (lrs_checkpoint_seconds);
5374 }
5375
5376 static void
checkpoint()5377 checkpoint ()
5378 {
5379 lrs_dump_state ();
5380 errcheck ("signal", signal (SIGUSR1, checkpoint));
5381 }
5382
5383 static void
die_gracefully()5384 die_gracefully ()
5385 {
5386 lrs_dump_state ();
5387
5388 exit (1);
5389 }
5390
5391 #endif
5392
5393 #ifndef TIMES
5394 /*
5395 * Not sure about the portability of this yet,
5396 * - db
5397 */
5398 #include <sys/resource.h>
5399 #define double_time(t) ((double)(t.tv_sec)+(double)(t.tv_usec)/1000000)
5400
5401 static void
ptimes()5402 ptimes ()
5403 {
5404 struct rusage rusage;
5405 getrusage (RUSAGE_SELF, &rusage);
5406 fprintf (lrs_ofp, "\n*%0.3fu %0.3fs %ldKb %ld flts %ld swaps %ld blks-in %ld blks-out \n",
5407 double_time (rusage.ru_utime),
5408 double_time (rusage.ru_stime),
5409 rusage.ru_maxrss, rusage.ru_majflt, rusage.ru_nswap,
5410 rusage.ru_inblock, rusage.ru_oublock);
5411 if(lrs_ofp != stdout)
5412 printf ("\n*%0.3fu %0.3fs %ldKb %ld flts %ld swaps %ld blks-in %ld blks-out \n",
5413 double_time (rusage.ru_utime),
5414 double_time (rusage.ru_stime),
5415 rusage.ru_maxrss, rusage.ru_majflt, rusage.ru_nswap,
5416 rusage.ru_inblock, rusage.ru_oublock);
5417 }
5418
get_time()5419 static double get_time()
5420 {
5421 struct rusage rusage;
5422 getrusage (RUSAGE_SELF, &rusage);
5423 return ( double_time (rusage.ru_utime));
5424
5425 }
5426
5427 #endif
5428
5429 /* Routines based on lp_solve */
5430
5431 void
lrs_set_row(lrs_dic * P,lrs_dat * Q,long row,long num[],long den[],long ineq)5432 lrs_set_row(lrs_dic *P, lrs_dat *Q, long row, long num[], long den[], long ineq)
5433 /* convert to lrs_mp then call lrs_set_row */
5434 {
5435 lrs_mp_vector Num, Den;
5436 long d;
5437 long j;
5438
5439 d = P->d;
5440
5441 Num=lrs_alloc_mp_vector(d+1);
5442 Den=lrs_alloc_mp_vector(d+1);
5443
5444 for (j=0;j<=d;j++)
5445 {
5446 itomp(num[j],Num[j]);
5447 itomp(den[j],Den[j]);
5448 }
5449
5450 lrs_set_row_mp(P,Q,row,Num,Den,ineq);
5451
5452 lrs_clear_mp_vector(Num,d+1);
5453 lrs_clear_mp_vector(Den,d+1);
5454
5455 }
5456
5457 void
lrs_set_row_mp(lrs_dic * P,lrs_dat * Q,long row,lrs_mp_vector num,lrs_mp_vector den,long ineq)5458 lrs_set_row_mp(lrs_dic *P, lrs_dat *Q, long row, lrs_mp_vector num, lrs_mp_vector den, long ineq)
5459 /* set row of dictionary using num and den arrays for rational input */
5460 /* ineq = 1 (GE) - ordinary row */
5461 /* = 0 (EQ) - linearity */
5462 {
5463 lrs_mp Temp, mpone;
5464 lrs_mp_vector oD; /* denominator for row */
5465
5466 long i, j;
5467
5468 /* assign local variables to structures */
5469
5470 lrs_mp_matrix A;
5471 lrs_mp_vector Gcd, Lcm;
5472 long hull;
5473 long m, d;
5474 lrs_alloc_mp(Temp); lrs_alloc_mp(mpone);
5475 hull = Q->hull;
5476 A = P->A;
5477 m = P->m;
5478 d = P->d;
5479 Gcd = Q->Gcd;
5480 Lcm = Q->Lcm;
5481
5482 oD = lrs_alloc_mp_vector (d);
5483 itomp (ONE, mpone);
5484 itomp (ONE, oD[0]);
5485
5486 i=row;
5487 itomp (ONE, Lcm[i]); /* Lcm of denominators */
5488 itomp (ZERO, Gcd[i]); /* Gcd of numerators */
5489 for (j = hull; j <= d; j++) /* hull data copied to cols 1..d */
5490 {
5491 copy( A[i][j],num[j-hull]);
5492 copy(oD[j],den[j-hull]);
5493 if (!one(oD[j]))
5494 lcm (Lcm[i], oD[j]); /* update lcm of denominators */
5495 copy (Temp, A[i][j]);
5496 gcd (Gcd[i], Temp); /* update gcd of numerators */
5497 }
5498
5499 if (hull)
5500 {
5501 itomp (ZERO, A[i][0]); /*for hull, we have to append an extra column of zeroes */
5502 if (!one (A[i][1]) || !one (oD[1])) /* all rows must have a one in column one */
5503 Q->polytope = FALSE;
5504 }
5505 if (!zero (A[i][hull])) /* for H-rep, are zero in column 0 */
5506 Q->homogeneous = FALSE; /* for V-rep, all zero in column 1 */
5507
5508 storesign (Gcd[i], POS);
5509 storesign (Lcm[i], POS);
5510 if (mp_greater (Gcd[i], mpone) || mp_greater (Lcm[i], mpone))
5511 for (j = 0; j <= d; j++)
5512 {
5513 exactdivint (A[i][j], Gcd[i], Temp); /*reduce numerators by Gcd */
5514 mulint (Lcm[i], Temp, Temp); /*remove denominators */
5515 exactdivint (Temp, oD[j], A[i][j]); /*reduce by former denominator */
5516 }
5517
5518 if ( ineq == EQ ) /* input is linearity */
5519 {
5520 Q->linearity[Q->nlinearity]=row;
5521 Q->nlinearity++;
5522 }
5523
5524 /* 2010.4.26 Set Gcd and Lcm for the non-existant rows when nonnegative set */
5525
5526
5527 if(Q->nonnegative && row==m)
5528 for(j=1;j<=d;j++)
5529 { itomp (ONE, Lcm[m+j]);
5530 itomp (ONE, Gcd[m+j]);
5531 }
5532
5533
5534 lrs_clear_mp_vector (oD,d);
5535 lrs_clear_mp(Temp); lrs_clear_mp(mpone);
5536 } /* end of lrs_set_row_mp */
5537
5538 void
lrs_set_obj(lrs_dic * P,lrs_dat * Q,long num[],long den[],long max)5539 lrs_set_obj(lrs_dic *P, lrs_dat *Q, long num[], long den[], long max)
5540 {
5541 long i;
5542
5543 if (max == MAXIMIZE)
5544 Q->maximize=TRUE;
5545 else
5546 {
5547 Q->minimize=TRUE;
5548 for(i=0;i<=P->d;i++)
5549 num[i]=-num[i];
5550 }
5551
5552 lrs_set_row(P,Q,0L,num,den,GE);
5553 }
5554
5555 void
lrs_set_obj_mp(lrs_dic * P,lrs_dat * Q,lrs_mp_vector num,lrs_mp_vector den,long max)5556 lrs_set_obj_mp(lrs_dic *P, lrs_dat *Q, lrs_mp_vector num, lrs_mp_vector den, long max)
5557 {
5558 long i;
5559
5560 if (max == MAXIMIZE)
5561 Q->maximize=TRUE;
5562 else
5563 {
5564 Q->minimize=TRUE;
5565 for(i=0;i<=P->d;i++)
5566 changesign(num[i]);
5567 }
5568
5569 lrs_set_row_mp(P,Q,0L,num,den,GE);
5570 }
5571
5572
5573 long
lrs_solve_lp(lrs_dic * P,lrs_dat * Q)5574 lrs_solve_lp(lrs_dic *P, lrs_dat *Q)
5575 /* user callable function to solve lp only */
5576 {
5577 lrs_mp_matrix Lin; /* holds input linearities if any are found */
5578 long col;
5579
5580 Q->lponly = TRUE;
5581
5582 if (!lrs_getfirstbasis (&P, Q, &Lin, FALSE))
5583 return FALSE;
5584
5585 /* There may have been column redundancy */
5586 /* If so the linearity space is obtained and redundant */
5587 /* columns are removed. User can access linearity space */
5588 /* from lrs_mp_matrix Lin dimensions nredundcol x d+1 */
5589
5590 for (col = 0; col < Q->nredundcol; col++) /* print linearity space */
5591 lrs_printoutput (Q, Lin[col]); /* Array Lin[][] holds the coeffs. */
5592
5593 return TRUE;
5594 } /* end of lrs_solve_lp */
5595
5596 long
dan_selectpivot(lrs_dic * P,lrs_dat * Q,long * r,long * s)5597 dan_selectpivot (lrs_dic * P, lrs_dat * Q, long *r, long *s)
5598 /* select pivot indices using dantzig simplex method */
5599 /* largest coefficient with lexicographic rule to avoid cycling */
5600 /* Bohdan Kaluzny's handiwork */
5601 /* returns TRUE if pivot found else FALSE */
5602 /* pivot variables are B[*r] C[*s] in locations Row[*r] Col[*s] */
5603 {
5604 long j,k,col;
5605 lrs_mp coeff;
5606 /* assign local variables to structures */
5607 lrs_mp_matrix A = P->A;
5608 long *Col = P->Col;
5609 long d = P->d;
5610
5611 /* printf("\n*dantzig"); */
5612 lrs_alloc_mp (coeff);
5613 *r = 0;
5614 *s = d;
5615 j = 0;
5616 k = 0;
5617
5618 itomp(0,coeff);
5619 /*find positive cost coef */
5620 while (k < d)
5621 {
5622 if(mp_greater(A[0][Col[k]],coeff))
5623 {
5624 j = k;
5625 copy(coeff,A[0][Col[j]]);
5626 }
5627 k++;
5628 }
5629
5630 if (positive(coeff)) /* pivot column found! */
5631 {
5632 *s = j;
5633 col = Col[j];
5634
5635 /*find min index ratio */
5636 *r = lrs_ratio (P, Q, col);
5637 if (*r != 0)
5638 {
5639 lrs_clear_mp(coeff);
5640 return (TRUE); /* pivot found */
5641 }
5642 }
5643 lrs_clear_mp(coeff);
5644 return (FALSE);
5645 } /* end of dan_selectpivot */
5646
5647
5648
5649 long
ran_selectpivot(lrs_dic * P,lrs_dat * Q,long * r,long * s)5650 ran_selectpivot (lrs_dic * P, lrs_dat * Q, long *r, long *s)
5651 /* select pivot indices using random edge rule */
5652 /* largest coefficient with lexicographic rule to avoid cycling */
5653 /* pivot variables are B[*r] C[*s] in locations Row[*r] Col[*s] */
5654 {
5655 long i,j,k,col,t;
5656 /* assign local variables to structures */
5657 lrs_mp_matrix A = P->A;
5658 long *Col = P->Col;
5659 long d = P->d;
5660 long *perm;
5661
5662 perm = (long *) calloc ((d + 1), sizeof (long));
5663 *r = 0;
5664 *s = d;
5665 k = 0;
5666 /* printf("\n*random edge"); */
5667
5668
5669 /* generate random permutation of 0..d-1 */
5670 for (i = 0; i < d; i++) perm[i] = i;
5671
5672 for ( i = 0; i < d; i++)
5673 {
5674 j = random() % (d-i) + i;
5675 t = perm[j]; perm[j] = perm[i]; perm[i] = t; // Swap i and j
5676 }
5677 if(Q->debug)
5678 {
5679 printf("\n perm: ");
5680 for (i = 0; i < d; i++) printf(" %ld",perm[i]);
5681 }
5682
5683 /*find first positive cost coef according to perm */
5684 while (k < d && !positive(A[0][Col[perm[k]]]))
5685 k++;
5686
5687 if ( k<d ) /* pivot column found! */
5688 {
5689 j=perm[k];
5690 *s = j;
5691 col = Col[j];
5692
5693 /*find min index ratio */
5694 *r = lrs_ratio (P, Q, col);
5695 if (*r != 0)
5696 {
5697 free(perm);
5698 return (TRUE); /* pivot found */
5699 }
5700 }
5701 free(perm);
5702 return (FALSE);
5703 } /* end of ran_selectpivot */
5704
5705
5706 long
phaseone(lrs_dic * P,lrs_dat * Q)5707 phaseone (lrs_dic * P, lrs_dat * Q)
5708 /* Do a dual pivot to get primal feasibility (pivot in X_0)*/
5709 /* Bohdan Kaluzny's handiwork */
5710 {
5711 long i, j, k;
5712 /* assign local variables to structures */
5713 lrs_mp_matrix A = P->A;
5714 long *Row = P->Row;
5715 long *Col = P->Col;
5716 long m, d;
5717 lrs_mp b_vector;
5718 lrs_alloc_mp (b_vector);
5719 m = P->m;
5720 d = P->d;
5721 i = 0;
5722 k = d+1;
5723
5724 itomp(0,b_vector);
5725
5726 fprintf (lrs_ofp, "\nLP: Phase One: Dual pivot on artificial variable");
5727
5728 /*find most negative b vector */
5729 while (k <= m)
5730 {
5731 if(mp_greater(b_vector,A[Row[k]][0]))
5732 {
5733 i = k;
5734 copy(b_vector,A[Row[i]][0]);
5735 }
5736 k++;
5737 }
5738
5739 if (negative(b_vector)) /* pivot row found! */
5740 {
5741 j = 0; /*find a positive entry for in row */
5742 while (j < d && !positive (A[Row[i]][Col[j]]))
5743 j++;
5744 if (j >= d)
5745 {
5746 lrs_clear_mp (b_vector);
5747 return (FALSE); /* no positive entry */
5748 }
5749 pivot (P, Q, i, j);
5750 update (P, Q, &i, &j);
5751 }
5752 lrs_clear_mp (b_vector);
5753 return (TRUE);
5754 }
5755
5756
5757 long
lrs_set_digits(long dec_digits)5758 lrs_set_digits(long dec_digits)
5759 {
5760 /* convert user specified decimal digits to mp digits */
5761
5762 if (dec_digits > 0)
5763 lrs_digits = DEC2DIG (dec_digits);
5764 if (lrs_digits > MAX_DIGITS)
5765 {
5766 fprintf (lrs_ofp, "\nDigits must be at most %ld\nChange MAX_DIGITS and recompile",
5767 DIG2DEC (MAX_DIGITS));
5768 fflush(stdout);
5769 return (FALSE);
5770 }
5771 return (TRUE);
5772 }
5773
5774 long
lrs_checkbound(lrs_dic * P,lrs_dat * Q)5775 lrs_checkbound(lrs_dic *P, lrs_dat *Q)
5776 {
5777 /* check bound on objective and return TRUE if exceeded */
5778
5779 if(!Q->bound)
5780 return FALSE;
5781
5782 if( Q->maximize && comprod(Q->boundn,P->objden,P->objnum,Q->boundd) == 1 )
5783 {
5784 if(Q->verbose)
5785 {
5786 prat(" \nObj value: ",P->objnum,P->objden);
5787 fprintf(lrs_ofp," Pruning ");
5788 }
5789 return TRUE;
5790 }
5791 if( Q->minimize && comprod(Q->boundn,P->objden,P->objnum,Q->boundd) == -1 )
5792 {
5793 if(Q->verbose)
5794 {
5795 prat(" \nObj value: ",P->objnum,P->objden);
5796 fprintf(lrs_ofp," Pruning ");
5797 }
5798 return TRUE;
5799 }
5800 return FALSE;
5801 }
5802
5803
5804 long
lrs_leaf(lrs_dic * P,lrs_dat * Q)5805 lrs_leaf(lrs_dic *P, lrs_dat *Q)
5806 {
5807 /* check if current dictionary is a leaf of reverse search tree */
5808 long col=0;
5809 long tmp=0;
5810
5811 while (col < P->d && !reverse(P,Q,&tmp,col))
5812 col++;
5813 if(col < P->d)
5814 return 0; /* dictionary is not a leaf */
5815 else
5816 return 1;
5817 }
5818
5819 /* prevent output flushes in mplrs */
lrs_open_outputblock(void)5820 void lrs_open_outputblock(void)
5821 {
5822 #ifdef PLRS
5823 open_outputblock();
5824 #endif
5825 }
5826
5827 /* re-enable output flushes in mplrs */
lrs_close_outputblock(void)5828 void lrs_close_outputblock(void)
5829 {
5830 #ifdef PLRS
5831 close_outputblock();
5832 #endif
5833 }
5834
lrs_post_output(const char * type,const char * data)5835 void lrs_post_output(const char *type, const char *data)
5836 {
5837 #ifdef PLRS
5838 post_output(type,data);
5839 #endif
5840 }
5841
lrs_return_unexplored(lrs_dic * P,lrs_dat * Q)5842 void lrs_return_unexplored(lrs_dic *P,lrs_dat *Q) /* send cobasis data for unexplored nodes */
5843
5844 {
5845
5846 #ifdef PLRS
5847 lrs_restart_dat R;
5848 int i;
5849 for (i = 0; i < P->d; i++)
5850 Q->temparray[i] = Q->inequality[P->C[i] - Q->lastdv];
5851 R.facet=Q->temparray;
5852 R.d=P->d;
5853 R.depth=P->depth;
5854 update_R(P,Q,&R);
5855 post_R(&R);
5856 #else
5857 if(Q->verbose)
5858 {
5859 lrs_printcobasis(P,Q,ZERO);
5860 fprintf(lrs_ofp," *unexplored");
5861 }
5862 #endif
5863 }
5864
lrs_overflow(int parm)5865 void lrs_overflow(int parm)
5866 {
5867 lrs_dic *P;
5868 lrs_dat *Q;
5869 char *restart;
5870 char *part;
5871
5872 int i;
5873 int try_restart=FALSE;
5874
5875 if (lrs_global_list[0] == NULL)
5876 {
5877 #ifdef PLRS
5878 post_output("warning","*lrs_overflow has null Q ");
5879 #else
5880 fprintf(stderr,"*lrs_overflow has null Q ");
5881 #endif
5882 lrs_exit(parm);
5883 }
5884
5885 Q = lrs_global_list[0]; /* db's cunningly hidden locations */
5886 P = Q->Qhead;
5887
5888 /* mplrs overflow handling */
5889
5890
5891 #if defined(PLRS) && !defined(GMP)
5892
5893 lrs_free_dic(P,Q); /* note Q is not freed here and is needed again */
5894 overflow=1;
5895 longjmp(buf1,1); /* return to lrsv2_main */
5896
5897 #elif defined(GMP) || defined(FLINT) /* should not be here, but just in case ... */
5898
5899 printf("\n*integer overflow for gmp or flint !?");
5900 lrs_free_all_memory (P,Q);
5901 lrs_exit(parm);
5902
5903 #endif
5904
5905 /* non mplrs overflow handling */
5906 /* lrs and redund restarted at the moment */
5907
5908 #ifdef MA
5909 if (strcmp(Q->fname,"lrs") == 0 || strcmp(Q->fname,"redund")==0)
5910 try_restart=TRUE;
5911 #endif
5912
5913 if(lrs_ifp != NULL)
5914 fclose(lrs_ifp);
5915
5916 if (!try_restart ) /* hard exit */
5917 {
5918 if (strcmp(BIT,"64bit")==0 )
5919 {
5920 fprintf(stderr,"\n*64bit integer overflow: try running 128bit or gmp versions\n");
5921 if (lrs_ofp != stdout)
5922 fprintf(lrs_ofp,"\n*64bit integer overflow: try running 128bit or gmp versions\n");
5923 }
5924 else
5925 {
5926 fprintf(stderr,"\n*128bit integer overflow: try running gmp version\n");
5927 if (lrs_ofp != stdout)
5928 fprintf(lrs_ofp,"\n*128bit integer overflow: try running gmp version\n");
5929 }
5930 lrs_exit(parm);
5931 }
5932
5933 /* try to restart */
5934
5935 if(overflow == 0) /* first overflow */
5936 {
5937 if (*tmpfilename != '\0' ) /* we made a temporary file for stdin */
5938 if(remove(tmpfilename) != 0)
5939 fprintf (lrs_ofp, "\nCould not delete temporary file");
5940 strncpy(tmpfilename,"/tmp/lrs_restartXXXXXX",PATH_MAX);
5941 /* XXX in principle this file descriptor should be used instead of the name */
5942 tmpfd = mkstemp(tmpfilename);
5943 }
5944 else
5945 strcpy(tmpfilename,infilename);
5946
5947 if( !pivoting || strcmp(Q->fname,"redund") == 0 || Q->getvolume) /* we make restart from original input */
5948 {
5949 overflow = 1L;
5950 lrs_cache_to_file(tmpfilename," ");
5951 }
5952 else
5953 {
5954 restart = (char *) malloc( Q->saved_d * 20+100);
5955 part = (char *) malloc( Q->saved_d * 20+100);
5956 overflow=2L;
5957 if(Q->hull)
5958 sprintf (restart," %ld %ld %ld ",
5959 Q->saved_count[2],Q->saved_count[0], Q->saved_depth);
5960 else
5961 sprintf (restart," %ld %ld %ld %ld ",
5962 Q->saved_count[1],Q->saved_count[0],Q->saved_count[2], Q->saved_depth);
5963
5964 for (i = 0; i < Q->saved_d; i++)
5965 {
5966 sprintf (part,"%ld ", Q->inequality[Q->saved_C[i] - Q->lastdv]);
5967 strcat(restart,part);
5968 }
5969 sprintf(part,"\nintegervertices %ld",Q->saved_count[4]);
5970 strcat(restart,part);
5971
5972 lrs_cache_to_file(tmpfilename,restart);
5973 free(restart); free(part);
5974 }
5975
5976 Q->m=P->m;
5977
5978 lrs_free_dic(P,Q); /* note Q is not freed here and is needed again */
5979
5980 if (lrs_ofp != NULL && lrs_ofp != stdout )
5981 {
5982 fclose (lrs_ofp);
5983 lrs_ofp=NULL;
5984 }
5985 close(tmpfd);
5986
5987 longjmp(buf1,1); /* return to lrsv2_main */
5988
5989 lrs_exit(parm); /* should not happen */
5990
5991 }
5992
5993
lrs_exit(int i)5994 void lrs_exit(int i)
5995 {
5996 fflush(stdout);
5997 exit(i);
5998 }
5999
lrs_free_all_memory(lrs_dic * P,lrs_dat * Q)6000 void lrs_free_all_memory(lrs_dic * P, lrs_dat * Q)
6001 {
6002
6003 if(Q->runs > 0)
6004 {
6005 free(Q->isave);
6006 free(Q->jsave);
6007 }
6008 if(P != NULL) /* may not have allocated P yet */
6009 {
6010 long savem=P->m; /* need this to clear Q*/
6011 lrs_free_dic (P,Q); /* deallocate lrs_dic */
6012 Q->m=savem;
6013 }
6014
6015 lrs_free_dat (Q); /* deallocate lrs_dat */
6016 #ifdef LRSLONG
6017 free(infile); /* we cached input file for possible restart */
6018 #endif
6019 return;
6020 }
6021
lrs_stdin_to_file(char * filename)6022 long lrs_stdin_to_file(char *filename)
6023 {
6024 FILE *fptr1, *fptr2;
6025 char c;
6026
6027 fptr1 = stdin;
6028 fptr2 = fopen(filename, "w");
6029 if (fptr2 == NULL)
6030 {
6031 printf("Cannot open file %s \n", filename);
6032 exit(0);
6033 }
6034
6035 c = fgetc(fptr1);
6036 while (c != EOF)
6037 {
6038 fputc(c, fptr2);
6039 c = fgetc(fptr1);
6040 }
6041
6042 fclose(fptr2);
6043 fptr2=NULL;
6044
6045 return 0;
6046 }
6047
lrs_file_to_cache(FILE * ifp)6048 long lrs_file_to_cache(FILE *ifp)
6049 {
6050 long ret;
6051
6052 if (ifp != NULL)
6053 if (fseek(ifp, 0L, SEEK_END) == 0)
6054 {
6055 ret = ftell(ifp);
6056 if (ret == -1)
6057 {
6058 fputs("*Error reading file", stderr);
6059 return 1;
6060 }
6061 infileLen = ret;
6062 infile = (char *) malloc(sizeof(char) * (infileLen + 1));
6063
6064 if (fseek(ifp, 0L, SEEK_SET) != 0)
6065 {
6066 fputs("*Error resetting input file", stderr);
6067 return 1;
6068 }
6069 infileLen = fread(infile, sizeof(char), infileLen, ifp);
6070 if ( ferror( ifp ) != 0 )
6071 {
6072 fputs("*Error reading input file", stderr);
6073 return 1;
6074 }
6075 else
6076 infile[infileLen++] = '\0'; /* Just to be safe. */
6077 }
6078 rewind(ifp);
6079 return 0;
6080 }
6081
lrs_cache_to_file(char * name,const char * restart)6082 long lrs_cache_to_file(char *name,const char *restart)
6083 {
6084 FILE *ofp = fopen(name, "wb");
6085
6086 if (ofp == NULL)
6087 {
6088 printf("*Error opening output file %s",name);
6089 return 1;
6090 }
6091 fwrite(infile, sizeof(char), infileLen, ofp);
6092
6093
6094 if(lrs_global_list[0]->count[2] > 1L && overflow==2)
6095 fprintf(ofp,"\nrestart %s",restart);
6096
6097 fclose(ofp);
6098 return 0;
6099
6100 }
6101
lrs_setup_R(lrs_dic * P,lrs_dat * Q,lrs_restart_dat * R)6102 void lrs_setup_R(lrs_dic *P, lrs_dat *Q, lrs_restart_dat *R)
6103 {
6104 int i;
6105
6106 R->d = P->d; /* length of R->facets */
6107 R->m = P->m; /* number of input rows*/
6108
6109 Q->startcob = (long int*) CALLOC ((R->d + R->m + 1), sizeof (long));
6110
6111 for(i=0;i<R->d;i++)
6112 Q->startcob[i]=Q->inequality[i];
6113
6114 if (strcmp (Q->fname, "redund") == 0)
6115 {
6116 R->redund=1;
6117 R->lrs=0;
6118 if(R->redineq == NULL)
6119 {
6120 R->redineq = (long int*) CALLOC ((R->m + 1), sizeof (long));
6121 for (i=0; i<= R->m; i++)
6122 R->redineq[i] = Q->redineq[i];
6123 }
6124 }
6125
6126 } /* lrs_setup_R */
6127
lrs_setup(int argc,char * argv[],lrs_dat ** Q,lrs_restart_dat * R)6128 lrs_dic *lrs_setup(int argc, char *argv[], lrs_dat **Q, lrs_restart_dat *R)
6129 /* allocate lrs_dat Q, lrs_dic P, read in the problem data and make a copy of P */
6130 {
6131 lrs_dic *P; /* structure for holding current dictionary and indices */
6132
6133 lrs_ifp = stdin;
6134 lrs_ofp = stdout;
6135
6136
6137 if(strncmp("redund",argv[0],6)==0)
6138 {
6139 if ( !lrs_init ("\n*redund:"))
6140 return NULL;
6141 }
6142 else if ( !lrs_init ("\n*lrs:"))
6143 return NULL;
6144
6145 *Q = lrs_alloc_dat ("LRS globals"); /* allocate and init structure for static problem data */
6146
6147 if (*Q == NULL)
6148 return NULL;
6149
6150 strcpy((*Q)->fname,"lrs");
6151
6152 if(strncmp("redund",argv[0],6)==0)
6153 strcpy((*Q)->fname,"redund");
6154
6155 if((*Q)->mplrs)
6156 (*Q)->messages=R->messages;
6157
6158
6159 if (!lrs_read_dat (*Q, argc, argv)) /* read first part of problem data to get dimensions */
6160 return NULL; /* and problem type: H- or V- input representation */
6161 P = lrs_alloc_dic (*Q); /* allocate and initialize lrs_dic */
6162 if (P == NULL )
6163 return NULL;
6164
6165 if (!lrs_read_dic (P, *Q)) /* read remainder of input to setup P and Q */
6166 return NULL;
6167
6168 return P;
6169 } /* lrs_setup */
6170
6171
lrs_reset(lrs_dic * P_orig,lrs_dat * Q,lrs_restart_dat * R)6172 lrs_dic *lrs_reset(lrs_dic *P_orig, lrs_dat *Q, lrs_restart_dat *R)
6173 {
6174 lrs_dic *P;
6175 long i;
6176
6177 itomp (ZERO, Q->Nvolume);
6178 itomp (ONE, Q->Dvolume);
6179 itomp (ZERO, Q->sumdet);
6180
6181 P=lrs_getdic (Q);
6182 Q->Qhead=P_orig;
6183 Q->Qtail=P_orig;
6184 if( P == P_orig)
6185 {
6186 if(Q->mplrs)
6187 lrs_post_output("warning", "*lrs_reset: copy_dict has dest=src -ignoring copy");
6188 else
6189 fprintf(stderr,"*lrs_reset: copy_dict has dest=src -ignoring copy");
6190 }
6191 copy_dict (Q,P,P_orig); /* restore original input */
6192 Q->Qhead=P;
6193 Q->Qtail=P;
6194
6195 /*if overiding, update Q from R */
6196
6197 if (R->lrs && R->overide == 1)
6198 {
6199 Q->messages=R->messages;
6200 if(R->maxdepth == -1)
6201 Q->maxdepth=MAXD;
6202 else
6203 Q->maxdepth=R->maxdepth;
6204 Q->mindepth=R->mindepth;
6205 Q->maxcobases=R->maxcobases;
6206 if (Q->maxcobases > 0)
6207 Q->maxcobases = Q->maxcobases + R->count[2];
6208 if(R->restart==1)
6209 {
6210 Q->restart=TRUE;
6211 if(!Q->lponly)
6212 Q->giveoutput=FALSE; /* supress first output */
6213
6214 for(i=0;i<R->d;i++)
6215 {
6216 Q->facet[i+Q->nlinearity]=R->facet[i];
6217 Q->inequality[i]=Q->startcob[i];
6218 }
6219 for(i=0;i<5;i++)
6220 {
6221 Q->count[i]=R->count[i];
6222 Q->startcount[i] = Q->count[i]; /* for mplrs subjob counts */
6223 }
6224 }
6225
6226 P->depth = R->depth;
6227 R->maxdepth=MAXD;
6228 }
6229
6230 if (R->redund)
6231 {
6232 for (i=0;i<=Q->m;i++)
6233 Q->redineq[i]=R->redineq[i];
6234
6235 Q->verifyredund=R->verifyredund;
6236 }
6237
6238 return P;
6239 }
6240
update_R(lrs_dic * P,lrs_dat * Q,lrs_restart_dat * R)6241 void update_R(lrs_dic *P, lrs_dat *Q, lrs_restart_dat *R)
6242 {
6243 int i;
6244 for (i=0;i<=4;i++)
6245 R->count[i]=Q->count[i];
6246 R->count[5]=Q->hull;
6247 if(Q->hull)
6248 R->count[6]=Q->nredundcol-Q->homogeneous;
6249 else
6250 R->count[6]=Q->nredundcol;
6251 R->count[7]=Q->deepest;
6252 return;
6253 }
6254
6255
6256 #ifdef GMP
6257 /* compiled with gmp arithmetic */
6258
lrsgmp_main(int argc,char * argv[],lrs_dic ** P_orig,lrs_dat ** Q,long overf,long stage,char * tmp,lrs_restart_dat * R)6259 long lrsgmp_main(int argc, char *argv[],lrs_dic **P_orig, lrs_dat **Q,long overf,long stage,char *tmp, lrs_restart_dat *R)
6260 {
6261 return lrsv2_main(argc,argv,P_orig,Q,overf,stage,tmp,R);
6262 }
6263
6264
6265 #elif defined(LRSLONG)
6266
6267 #ifdef B128
6268
lrs2_main(int argc,char * argv[],lrs_dic ** P_orig,lrs_dat ** Q,long overf,long stage,char * tmp,lrs_restart_dat * R)6269 long lrs2_main(int argc, char *argv[],lrs_dic **P_orig, lrs_dat **Q,long overf,long stage,char *tmp, lrs_restart_dat *R)
6270 {
6271 return lrsv2_main(argc,argv,P_orig,Q,overf,stage,tmp,R);
6272 }
6273
6274
6275 #else
6276
lrs1_main(int argc,char * argv[],lrs_dic ** P_orig,lrs_dat ** Q,long overf,long stage,char * tmp,lrs_restart_dat * R)6277 long lrs1_main(int argc, char *argv[],lrs_dic **P_orig, lrs_dat **Q,long overf,long stage,char *tmp, lrs_restart_dat *R)
6278 {
6279 return lrsv2_main(argc,argv,P_orig,Q,overf,stage,tmp,R);
6280 }
6281
6282 #endif
6283 #endif
6284
6285
lrs_main(int argc,char * argv[])6286 long lrs_main(int argc, char *argv[])
6287 /* legacy version, replaced by lrsv2_main but still maintained */
6288
6289 {
6290 lrs_dic *P;
6291 lrs_dat *Q;
6292 lrs_restart_dat *R;
6293 char* tmp; /* when overflow occurs a new input file name is returned */
6294 long overfl=0; /* =0 no overflow =1 restart overwrite =2 restart append */
6295
6296
6297 P=NULL;
6298 Q=NULL;
6299 R=NULL;
6300 tmp=NULL;
6301
6302 R = lrs_alloc_restart();
6303 if (R == NULL)
6304 exit(1);
6305
6306 overfl=lrsv2_main(argc,argv,&P,&Q,0,0,tmp,R); /* set up, read input, no run */
6307
6308 if(overfl == -1) /* lrs_setup failed due to bad input file etc. - no cleanup*/
6309 return 0;
6310 if(overfl == 0)
6311 lrsv2_main(argc,argv,&P,&Q,0,1,tmp,R); /* standard lrs run - argc, argv, R not used */
6312
6313 lrsv2_main(argc,argv,&P,&Q,0,2,tmp,R); /* free memory and close, does not access argc, argv */
6314
6315 free(R->facet);
6316 if (R->redund)
6317 free(R->redineq);
6318
6319 return 0;
6320 }
6321
lrsv2_main(int argc,char * argv[],lrs_dic ** P_orig,lrs_dat ** Q,long overf,long stage,char * tmp,lrs_restart_dat * R)6322 long lrsv2_main(int argc, char *argv[],lrs_dic **P_orig, lrs_dat **Q,long overf,long stage,char *tmp, lrs_restart_dat *R)
6323
6324 /* compiled independently with all supported arithmetic packages */
6325 /* should be called from one of lrsX_main where X is an arithmetic package */
6326
6327 {
6328 lrs_dic *P; /* structure for holding current dictionary and indices */
6329 int i;
6330
6331 overflow=overf;
6332
6333 if (!setjmp(buf1)) /* normal processing - jump to end if overflow occurs */
6334 {
6335 /* initial call: allocate lrs_dat, lrs_dic and set up the problem - no run */
6336 if(stage==0)
6337 {
6338 *P_orig=lrs_setup(argc,argv,Q,R);
6339 if(*P_orig==NULL)
6340 {
6341 fprintf(stderr,"\n*lrs_setup failed\n");
6342 fflush(stderr);
6343 return -1;
6344 }
6345 lrs_setup_R(*P_orig,*Q,R);
6346 return 0;
6347 }
6348
6349 /* reverse search runs: restore P and update Q as needed from R */
6350 if(stage==1)
6351 {
6352 P=lrs_reset(*P_orig,*Q,R); /* restore P and reset Q from R */
6353 if(P==NULL)
6354 return -1;
6355 if(overf==2)
6356 (*Q)->giveoutput=FALSE; /* suppress first output */
6357
6358 if(R->redund)
6359 {
6360 redund_run(P,*Q);
6361 return 0;
6362 }
6363 lrs_run(P,*Q); /* do reverse search */
6364
6365 update_R(P,*Q,R); /* update counts for mplrs */
6366 return 0;
6367 }
6368
6369 /* final cleanup */
6370 if(stage == 2 )
6371 {
6372
6373 (*Q)->Qhead=*P_orig;
6374 (*Q)->Qtail=*P_orig;
6375
6376 lrs_free_all_memory(*P_orig,*Q);
6377 lrs_close ("lrs:");
6378
6379 return 0;
6380 }
6381 }
6382
6383 /* overflow occurred */
6384
6385 if (R->redund )
6386 {
6387 if(R->redineq != NULL)
6388 {
6389 overflow=2;
6390 for (i=0;i<=R->m;i++) /*i=0 contains next ineq to process */
6391 R->redineq[i]=(*Q)->redineq[i];
6392 }
6393 lrs_clear_mp_matrix((*Q)->Ain,(*P_orig)->m_A,(*P_orig)->d);
6394 }
6395
6396 if (tmp != NULL)
6397 strcpy(tmp,tmpfilename);
6398
6399 (*Q)->Qhead=*P_orig;
6400 (*Q)->Qtail=*P_orig;
6401 lrs_free_all_memory(*P_orig,*Q);
6402 fflush(lrs_ofp);
6403
6404 return overflow; /* overflow */
6405 } /* lrsv2_main */
6406
6407 void
lrs_warning(lrs_dat * Q,char * type,char * ss)6408 lrs_warning(lrs_dat *Q, char* type, char* ss)
6409 {
6410 if(Q->messages)
6411 {
6412 if(Q->mplrs)
6413 lrs_post_output(type,ss);
6414 else
6415 {
6416 fprintf (lrs_ofp, "\n%s",ss);
6417 if(lrs_ofp != stdout)
6418 fprintf (stderr, "\n%s",ss);
6419 }
6420 }
6421 }
6422