1 /* cddlp.c: dual simplex method c-code
2 written by Komei Fukuda, fukuda@math.ethz.ch
3 Version 0.94i, March 9, 2018
4 */
5
6 /* cddlp.c : C-Implementation of the dual simplex method for
7 solving an LP: max/min A_(m-1).x subject to x in P, where
8 P= {x : A_i.x >= 0, i=0,...,m-2, and x_0=1}, and
9 A_i is the i-th row of an m x n matrix A.
10 Please read COPYING (GNU General Public Licence) and
11 the manual cddlibman.tex for detail.
12 */
13
14 #include "setoper.h" /* set operation library header (Ver. May 18, 2000 or later) */
15 #include "cdd.h"
16 #include "splitmix64.h"
17 #include <stdio.h>
18 #include <stdlib.h>
19 #include <time.h>
20 #include <math.h>
21 #include <string.h>
22
23 #if defined GMPRATIONAL
24 #include "cdd_f.h"
25 #endif
26
27 #define dd_CDDLPVERSION "Version 0.94b (August 25, 2005)"
28
29 #define dd_FALSE 0
30 #define dd_TRUE 1
31
32 typedef set_type rowset; /* set_type defined in setoper.h */
33 typedef set_type colset;
34
35 void dd_CrissCrossSolve(dd_LPPtr lp,dd_ErrorType *);
36 void dd_DualSimplexSolve(dd_LPPtr lp,dd_ErrorType *);
37 void dd_CrissCrossMinimize(dd_LPPtr,dd_ErrorType *);
38 void dd_CrissCrossMaximize(dd_LPPtr,dd_ErrorType *);
39 void dd_DualSimplexMinimize(dd_LPPtr,dd_ErrorType *);
40 void dd_DualSimplexMaximize(dd_LPPtr,dd_ErrorType *);
41 void dd_FindLPBasis(dd_rowrange,dd_colrange,dd_Amatrix,dd_Bmatrix,dd_rowindex,dd_rowset,
42 dd_colindex,dd_rowindex,dd_rowrange,dd_colrange,
43 dd_colrange *,int *,dd_LPStatusType *,long *);
44 void dd_FindDualFeasibleBasis(dd_rowrange,dd_colrange,dd_Amatrix,dd_Bmatrix,dd_rowindex,
45 dd_colindex,long *,dd_rowrange,dd_colrange,dd_boolean,
46 dd_colrange *,dd_ErrorType *,dd_LPStatusType *,long *, long maxpivots);
47
48
49 #ifdef GMPRATIONAL
50 void dd_BasisStatus(ddf_LPPtr lpf, dd_LPPtr lp, dd_boolean*);
51 void dd_BasisStatusMinimize(dd_rowrange,dd_colrange, dd_Amatrix,dd_Bmatrix,dd_rowset,
52 dd_rowrange,dd_colrange,ddf_LPStatusType,mytype *,dd_Arow,dd_Arow,dd_rowset,ddf_colindex,
53 ddf_rowrange,ddf_colrange,dd_colrange *,long *, int *, int *);
54 void dd_BasisStatusMaximize(dd_rowrange,dd_colrange,dd_Amatrix,dd_Bmatrix,dd_rowset,
55 dd_rowrange,dd_colrange,ddf_LPStatusType,mytype *,dd_Arow,dd_Arow,dd_rowset,ddf_colindex,
56 ddf_rowrange,ddf_colrange,dd_colrange *,long *, int *, int *);
57 #endif
58
59 void dd_WriteBmatrix(FILE *f,dd_colrange d_size,dd_Bmatrix T);
60 void dd_SetNumberType(char *line,dd_NumberType *number,dd_ErrorType *Error);
61 void dd_ComputeRowOrderVector2(dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,
62 dd_rowindex OV,dd_RowOrderType ho,unsigned int rseed);
63 void dd_SelectPreorderedNext2(dd_rowrange m_size,dd_colrange d_size,
64 rowset excluded,dd_rowindex OV,dd_rowrange *hnext);
65 void dd_SetSolutions(dd_rowrange,dd_colrange,
66 dd_Amatrix,dd_Bmatrix,dd_rowrange,dd_colrange,dd_LPStatusType,
67 mytype *,dd_Arow,dd_Arow,dd_rowset,dd_colindex,dd_rowrange,dd_colrange,dd_rowindex);
68
69 void dd_WriteTableau(FILE *,dd_rowrange,dd_colrange,dd_Amatrix,dd_Bmatrix,
70 dd_colindex,dd_rowindex);
71
72 void dd_WriteSignTableau(FILE *,dd_rowrange,dd_colrange,dd_Amatrix,dd_Bmatrix,
73 dd_colindex,dd_rowindex);
74
75
dd_CopyLPSolution(dd_LPPtr lp)76 dd_LPSolutionPtr dd_CopyLPSolution(dd_LPPtr lp)
77 {
78 dd_LPSolutionPtr lps;
79 dd_colrange j;
80 long i;
81
82 lps=(dd_LPSolutionPtr) calloc(1,sizeof(dd_LPSolutionType));
83 for (i=1; i<=dd_filenamelen; i++) lps->filename[i-1]=lp->filename[i-1];
84 lps->objective=lp->objective;
85 lps->solver=lp->solver;
86 lps->m=lp->m;
87 lps->d=lp->d;
88 lps->numbtype=lp->numbtype;
89
90 lps->LPS=lp->LPS; /* the current solution status */
91 dd_init(lps->optvalue);
92 dd_set(lps->optvalue,lp->optvalue); /* optimal value */
93 dd_InitializeArow(lp->d+1,&(lps->sol));
94 dd_InitializeArow(lp->d+1,&(lps->dsol));
95 lps->nbindex=(long*) calloc((lp->d)+1,sizeof(long)); /* dual solution */
96 for (j=0; j<=lp->d; j++){
97 dd_set(lps->sol[j],lp->sol[j]);
98 dd_set(lps->dsol[j],lp->dsol[j]);
99 lps->nbindex[j]=lp->nbindex[j];
100 }
101 lps->pivots[0]=lp->pivots[0];
102 lps->pivots[1]=lp->pivots[1];
103 lps->pivots[2]=lp->pivots[2];
104 lps->pivots[3]=lp->pivots[3];
105 lps->pivots[4]=lp->pivots[4];
106 lps->total_pivots=lp->total_pivots;
107
108 return lps;
109 }
110
111
dd_CreateLPData(dd_LPObjectiveType obj,dd_NumberType nt,dd_rowrange m,dd_colrange d)112 dd_LPPtr dd_CreateLPData(dd_LPObjectiveType obj,
113 dd_NumberType nt,dd_rowrange m,dd_colrange d)
114 {
115 dd_LPType *lp;
116
117 lp=(dd_LPPtr) calloc(1,sizeof(dd_LPType));
118 lp->solver=dd_choiceLPSolverDefault; /* set the default lp solver */
119 lp->d=d;
120 lp->m=m;
121 lp->numbtype=nt;
122 lp->objrow=m;
123 lp->rhscol=1L;
124 lp->objective=dd_LPnone;
125 lp->LPS=dd_LPSundecided;
126 lp->eqnumber=0; /* the number of equalities */
127
128 lp->nbindex=(long*) calloc(d+1,sizeof(long));
129 lp->given_nbindex=(long*) calloc(d+1,sizeof(long));
130 set_initialize(&(lp->equalityset),m);
131 /* i must be in the set iff i-th row is equality . */
132
133 lp->redcheck_extensive=dd_FALSE; /* this is on only for RedundantExtensive */
134 lp->ired=0;
135 set_initialize(&(lp->redset_extra),m);
136 /* i is in the set if i-th row is newly recognized redundant (during the checking the row ired). */
137 set_initialize(&(lp->redset_accum),m);
138 /* i is in the set if i-th row is recognized redundant (during the checking the row ired). */
139 set_initialize(&(lp->posset_extra),m);
140 /* i is in the set if i-th row is recognized non-linearity (during the course of computation). */
141 lp->lexicopivot=dd_choiceLexicoPivotQ; /* dd_choice... is set in dd_set_global_constants() */
142
143 lp->m_alloc=lp->m+2;
144 lp->d_alloc=lp->d+2;
145 lp->objective=obj;
146 dd_InitializeBmatrix(lp->d_alloc,&(lp->B));
147 dd_InitializeAmatrix(lp->m_alloc,lp->d_alloc,&(lp->A));
148 dd_InitializeArow(lp->d_alloc,&(lp->sol));
149 dd_InitializeArow(lp->d_alloc,&(lp->dsol));
150 dd_init(lp->optvalue);
151 return lp;
152 }
153
154
dd_Matrix2LP(dd_MatrixPtr M,dd_ErrorType * err)155 dd_LPPtr dd_Matrix2LP(dd_MatrixPtr M, dd_ErrorType *err)
156 {
157 dd_rowrange m, i, irev, linc;
158 dd_colrange d, j;
159 dd_LPType *lp;
160 dd_boolean localdebug=dd_FALSE;
161
162 *err=dd_NoError;
163 linc=set_card(M->linset);
164 m=M->rowsize+1+linc;
165 /* We represent each equation by two inequalities.
166 This is not the best way but makes the code simple. */
167 d=M->colsize;
168 if (localdebug) fprintf(stderr,"number of equalities = %ld\n", linc);
169
170 lp=dd_CreateLPData(M->objective, M->numbtype, m, d);
171 lp->Homogeneous = dd_TRUE;
172 lp->eqnumber=linc; /* this records the number of equations */
173
174 irev=M->rowsize; /* the first row of the linc reversed inequalities. */
175 for (i = 1; i <= M->rowsize; i++) {
176 if (set_member(i, M->linset)) {
177 irev=irev+1;
178 set_addelem(lp->equalityset,i); /* it is equality. */
179 /* the reversed row irev is not in the equality set. */
180 for (j = 1; j <= M->colsize; j++) {
181 dd_neg(lp->A[irev-1][j-1],M->matrix[i-1][j-1]);
182 } /*of j*/
183 if (localdebug) fprintf(stderr,"equality row %ld generates the reverse row %ld.\n",i,irev);
184 }
185 for (j = 1; j <= M->colsize; j++) {
186 dd_set(lp->A[i-1][j-1],M->matrix[i-1][j-1]);
187 if (j==1 && i<M->rowsize && dd_Nonzero(M->matrix[i-1][j-1])) lp->Homogeneous = dd_FALSE;
188 } /*of j*/
189 } /*of i*/
190 for (j = 1; j <= M->colsize; j++) {
191 dd_set(lp->A[m-1][j-1],M->rowvec[j-1]); /* objective row */
192 } /*of j*/
193
194 return lp;
195 }
196
dd_Matrix2Feasibility(dd_MatrixPtr M,dd_ErrorType * err)197 dd_LPPtr dd_Matrix2Feasibility(dd_MatrixPtr M, dd_ErrorType *err)
198 /* Load a matrix to create an LP object for feasibility. It is
199 essentially the dd_Matrix2LP except that the objject function
200 is set to identically ZERO (maximization).
201
202 */
203 /* 094 */
204 {
205 dd_rowrange m, linc;
206 dd_colrange j;
207 dd_LPType *lp;
208
209 *err=dd_NoError;
210 linc=set_card(M->linset);
211 m=M->rowsize+1+linc;
212 /* We represent each equation by two inequalities.
213 This is not the best way but makes the code simple. */
214
215 lp=dd_Matrix2LP(M, err);
216 lp->objective = dd_LPmax; /* since the objective is zero, this is not important */
217 for (j = 1; j <= M->colsize; j++) {
218 dd_set(lp->A[m-1][j-1],dd_purezero); /* set the objective to zero. */
219 } /*of j*/
220
221 return lp;
222 }
223
dd_Matrix2Feasibility2(dd_MatrixPtr M,dd_rowset R,dd_rowset S,dd_ErrorType * err)224 dd_LPPtr dd_Matrix2Feasibility2(dd_MatrixPtr M, dd_rowset R, dd_rowset S, dd_ErrorType *err)
225 /* Load a matrix to create an LP object for feasibility with additional equality and
226 strict inequality constraints given by R and S. There are three types of inequalities:
227
228 b_r + A_r x = 0 Linearity (Equations) specified by M
229 b_s + A_s x > 0 Strict Inequalities specified by row index set S
230 b_t + A_t x >= 0 The rest inequalities in M
231
232 Where the linearity is considered here as the union of linearity specified by
233 M and the additional set R. When S contains any linearity rows, those
234 rows are considered linearity (equation). Thus S does not overlide linearity.
235 To find a feasible solution, we set an LP
236
237 maximize z
238 subject to
239 b_r + A_r x = 0 all r in Linearity
240 b_s + A_s x - z >= 0 for all s in S
241 b_t + A_t x >= 0 for all the rest rows t
242 1 - z >= 0 to make the LP bounded.
243
244 Clearly, the feasibility problem has a solution iff the LP has a positive optimal value.
245 The variable z will be the last variable x_{d+1}.
246
247 */
248 /* 094 */
249 {
250 dd_rowrange m, i, irev, linc;
251 dd_colrange d, j;
252 dd_LPType *lp;
253 dd_rowset L;
254 dd_boolean localdebug=dd_FALSE;
255
256 *err=dd_NoError;
257 set_initialize(&L, M->rowsize);
258 set_uni(L,M->linset,R);
259 linc=set_card(L);
260 m=M->rowsize+1+linc+1;
261 /* We represent each equation by two inequalities.
262 This is not the best way but makes the code simple. */
263 d=M->colsize+1;
264 if (localdebug) fprintf(stderr,"number of equalities = %ld\n", linc);
265
266 lp=dd_CreateLPData(dd_LPmax, M->numbtype, m, d);
267 lp->Homogeneous = dd_TRUE;
268 lp->eqnumber=linc; /* this records the number of equations */
269
270 irev=M->rowsize; /* the first row of the linc reversed inequalities. */
271 for (i = 1; i <= M->rowsize; i++) {
272 if (set_member(i, L)) {
273 irev=irev+1;
274 set_addelem(lp->equalityset,i); /* it is equality. */
275 /* the reversed row irev is not in the equality set. */
276 for (j = 1; j <= M->colsize; j++) {
277 dd_neg(lp->A[irev-1][j-1],M->matrix[i-1][j-1]);
278 } /*of j*/
279 if (localdebug) fprintf(stderr,"equality row %ld generates the reverse row %ld.\n",i,irev);
280 } else if (set_member(i, S)) {
281 dd_set(lp->A[i-1][M->colsize],dd_minusone);
282 }
283 for (j = 1; j <= M->colsize; j++) {
284 dd_set(lp->A[i-1][j-1],M->matrix[i-1][j-1]);
285 if (j==1 && i<M->rowsize && dd_Nonzero(M->matrix[i-1][j-1])) lp->Homogeneous = dd_FALSE;
286 } /*of j*/
287 } /*of i*/
288 for (j = 1; j <= d; j++) {
289 dd_set(lp->A[m-2][j-1],dd_purezero); /* initialize */
290 } /*of j*/
291 dd_set(lp->A[m-2][0],dd_one); /* the bounding constraint. */
292 dd_set(lp->A[m-2][M->colsize],dd_minusone); /* the bounding constraint. */
293 for (j = 1; j <= d; j++) {
294 dd_set(lp->A[m-1][j-1],dd_purezero); /* initialize */
295 } /*of j*/
296 dd_set(lp->A[m-1][M->colsize],dd_one); /* maximize z */
297
298 set_free(L);
299 return lp;
300 }
301
302
303
dd_FreeLPData(dd_LPPtr lp)304 void dd_FreeLPData(dd_LPPtr lp)
305 {
306 if ((lp)!=NULL){
307 dd_clear(lp->optvalue);
308 dd_FreeArow(lp->d_alloc,lp->dsol);
309 dd_FreeArow(lp->d_alloc,lp->sol);
310 dd_FreeBmatrix(lp->d_alloc,lp->B);
311 dd_FreeAmatrix(lp->m_alloc,lp->d_alloc,lp->A);
312 set_free(lp->equalityset);
313 set_free(lp->redset_extra);
314 set_free(lp->redset_accum);
315 set_free(lp->posset_extra);
316 free(lp->nbindex);
317 free(lp->given_nbindex);
318 free(lp);
319 }
320 }
321
dd_FreeLPSolution(dd_LPSolutionPtr lps)322 void dd_FreeLPSolution(dd_LPSolutionPtr lps)
323 {
324 if (lps!=NULL){
325 free(lps->nbindex);
326 dd_FreeArow(lps->d+1,lps->dsol);
327 dd_FreeArow(lps->d+1,lps->sol);
328 dd_clear(lps->optvalue);
329
330 free(lps);
331 }
332 }
333
dd_LPReverseRow(dd_LPPtr lp,dd_rowrange i)334 int dd_LPReverseRow(dd_LPPtr lp, dd_rowrange i)
335 {
336 dd_colrange j;
337 int success=0;
338
339 if (i>=1 && i<=lp->m){
340 lp->LPS=dd_LPSundecided;
341 for (j=1; j<=lp->d; j++) {
342 dd_neg(lp->A[i-1][j-1],lp->A[i-1][j-1]);
343 /* negating the i-th constraint of A */
344 }
345 success=1;
346 }
347 return success;
348 }
349
dd_LPReplaceRow(dd_LPPtr lp,dd_rowrange i,dd_Arow a)350 int dd_LPReplaceRow(dd_LPPtr lp, dd_rowrange i, dd_Arow a)
351 {
352 dd_colrange j;
353 int success=0;
354
355 if (i>=1 && i<=lp->m){
356 lp->LPS=dd_LPSundecided;
357 for (j=1; j<=lp->d; j++) {
358 dd_set(lp->A[i-1][j-1],a[j-1]);
359 /* replacing the i-th constraint by a */
360 }
361 success=1;
362 }
363 return success;
364 }
365
dd_LPCopyRow(dd_LPPtr lp,dd_rowrange i)366 dd_Arow dd_LPCopyRow(dd_LPPtr lp, dd_rowrange i)
367 {
368 dd_colrange j;
369 dd_Arow a;
370
371 if (i>=1 && i<=lp->m){
372 dd_InitializeArow(lp->d, &a);
373 for (j=1; j<=lp->d; j++) {
374 dd_set(a[j-1],lp->A[i-1][j-1]);
375 /* copying the i-th row to a */
376 }
377 }
378 return a;
379 }
380
381
dd_SetNumberType(char * line,dd_NumberType * number,dd_ErrorType * Error)382 void dd_SetNumberType(char *line,dd_NumberType *number,dd_ErrorType *Error)
383 {
384 if (strncmp(line,"integer",7)==0) {
385 *number = dd_Integer;
386 return;
387 }
388 else if (strncmp(line,"rational",8)==0) {
389 *number = dd_Rational;
390 return;
391 }
392 else if (strncmp(line,"real",4)==0) {
393 *number = dd_Real;
394 return;
395 }
396 else {
397 *number=dd_Unknown;
398 *Error=dd_ImproperInputFormat;
399 }
400 }
401
402
dd_WriteTableau(FILE * f,dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_Bmatrix T,dd_colindex nbindex,dd_rowindex bflag)403 void dd_WriteTableau(FILE *f,dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_Bmatrix T,
404 dd_colindex nbindex,dd_rowindex bflag)
405 /* Write the tableau A.T */
406 {
407 dd_colrange j;
408 dd_rowrange i;
409 mytype x;
410
411 dd_init(x);
412 fprintf(f," %ld %ld real\n",m_size,d_size);
413 fprintf(f," |");
414 for (j=1; j<= d_size; j++) {
415 fprintf(f," %ld",nbindex[j]);
416 } fprintf(f,"\n");
417 for (j=1; j<= d_size+1; j++) {
418 fprintf(f," ----");
419 } fprintf(f,"\n");
420 for (i=1; i<= m_size; i++) {
421 fprintf(f," %3ld(%3ld) |",i,bflag[i]);
422 for (j=1; j<= d_size; j++) {
423 dd_TableauEntry(&x,m_size,d_size,A,T,i,j);
424 dd_WriteNumber(f,x);
425 }
426 fprintf(f,"\n");
427 }
428 fprintf(f,"end\n");
429 dd_clear(x);
430 }
431
dd_WriteSignTableau(FILE * f,dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_Bmatrix T,dd_colindex nbindex,dd_rowindex bflag)432 void dd_WriteSignTableau(FILE *f,dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_Bmatrix T,
433 dd_colindex nbindex,dd_rowindex bflag)
434 /* Write the sign tableau A.T */
435 {
436 dd_colrange j;
437 dd_rowrange i;
438 mytype x;
439
440 dd_init(x);
441 fprintf(f," %ld %ld real\n",m_size,d_size);
442 fprintf(f," |");
443 for (j=1; j<= d_size; j++) {
444 fprintf(f,"%3ld",nbindex[j]);
445 } fprintf(f,"\n ------- | ");
446 for (j=1; j<= d_size; j++) {
447 fprintf(f,"---");
448 } fprintf(f,"\n");
449 for (i=1; i<= m_size; i++) {
450 fprintf(f," %3ld(%3ld) |",i,bflag[i]);
451 for (j=1; j<= d_size; j++) {
452 dd_TableauEntry(&x,m_size,d_size,A,T,i,j);
453 if (dd_Positive(x)) fprintf(f, " +");
454 else if (dd_Negative(x)) fprintf(f, " -");
455 else fprintf(f, " 0");
456 }
457 fprintf(f,"\n");
458 }
459 fprintf(f,"end\n");
460 dd_clear(x);
461 }
462
dd_WriteSignTableau2(FILE * f,dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_Bmatrix T,dd_colindex nbindex_ref,dd_colindex nbindex,dd_rowindex bflag)463 void dd_WriteSignTableau2(FILE *f,dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_Bmatrix T,
464 dd_colindex nbindex_ref, dd_colindex nbindex,dd_rowindex bflag)
465 /* Write the sign tableau A.T and the reference basis */
466 {
467 dd_colrange j;
468 dd_rowrange i;
469 mytype x;
470
471 dd_init(x);
472 fprintf(f," %ld %ld real\n",m_size,d_size);
473 fprintf(f," |");
474 for (j=1; j<= d_size; j++) fprintf(f,"%3ld",nbindex_ref[j]);
475 fprintf(f,"\n |");
476 for (j=1; j<= d_size; j++) {
477 fprintf(f,"%3ld",nbindex[j]);
478 } fprintf(f,"\n ------- | ");
479 for (j=1; j<= d_size; j++) {
480 fprintf(f,"---");
481 } fprintf(f,"\n");
482 for (i=1; i<= m_size; i++) {
483 fprintf(f," %3ld(%3ld) |",i,bflag[i]);
484 for (j=1; j<= d_size; j++) {
485 dd_TableauEntry(&x,m_size,d_size,A,T,i,j);
486 if (dd_Positive(x)) fprintf(f, " +");
487 else if (dd_Negative(x)) fprintf(f, " -");
488 else fprintf(f, " 0");
489 }
490 fprintf(f,"\n");
491 }
492 fprintf(f,"end\n");
493 dd_clear(x);
494 }
495
496
dd_GetRedundancyInformation(dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_Bmatrix T,dd_colindex nbindex,dd_rowindex bflag,dd_rowset redset)497 void dd_GetRedundancyInformation(dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_Bmatrix T,
498 dd_colindex nbindex,dd_rowindex bflag, dd_rowset redset)
499 /* Some basic variables that are forced to be nonnegative will be output. These are
500 variables whose dictionary row components are all nonnegative. */
501 {
502 dd_colrange j;
503 dd_rowrange i;
504 mytype x;
505 dd_boolean red=dd_FALSE,localdebug=dd_FALSE;
506 long numbred=0;
507
508 dd_init(x);
509 for (i=1; i<= m_size; i++) {
510 red=dd_TRUE;
511 for (j=1; j<= d_size; j++) {
512 dd_TableauEntry(&x,m_size,d_size,A,T,i,j);
513 if (red && dd_Negative(x)) red=dd_FALSE;
514 }
515 if (bflag[i]<0 && red) {
516 numbred+=1;
517 set_addelem(redset,i);
518 }
519 }
520 if (localdebug) fprintf(stderr,"\ndd_GetRedundancyInformation: %ld redundant rows over %ld\n",numbred, m_size);
521 dd_clear(x);
522 }
523
524
dd_SelectDualSimplexPivot(dd_rowrange m_size,dd_colrange d_size,int Phase1,dd_Amatrix A,dd_Bmatrix T,dd_rowindex OV,dd_colindex nbindex_ref,dd_colindex nbindex,dd_rowindex bflag,dd_rowrange objrow,dd_colrange rhscol,dd_boolean lexicopivot,dd_rowrange * r,dd_colrange * s,int * selected,dd_LPStatusType * lps)525 void dd_SelectDualSimplexPivot(dd_rowrange m_size,dd_colrange d_size,
526 int Phase1,dd_Amatrix A,dd_Bmatrix T,dd_rowindex OV,
527 dd_colindex nbindex_ref, dd_colindex nbindex,dd_rowindex bflag,
528 dd_rowrange objrow,dd_colrange rhscol, dd_boolean lexicopivot,
529 dd_rowrange *r,dd_colrange *s,int *selected,dd_LPStatusType *lps)
530 {
531 /* selects a dual simplex pivot (*r,*s) if the current
532 basis is dual feasible and not optimal. If not dual feasible,
533 the procedure returns *selected=dd_FALSE and *lps=LPSundecided.
534 If Phase1=dd_TRUE, the RHS column will be considered as the negative
535 of the column of the largest variable (==m_size). For this case, it is assumed
536 that the caller used the auxiliary row (with variable m_size) to make the current
537 dictionary dual feasible before calling this routine so that the nonbasic
538 column for m_size corresponds to the auxiliary variable.
539 */
540 dd_boolean colselected=dd_FALSE,rowselected=dd_FALSE,
541 dualfeasible=dd_TRUE,localdebug=dd_FALSE;
542 dd_rowrange i,iref;
543 dd_colrange j,k;
544 mytype val,valn, minval,rat,minrat;
545 static dd_Arow rcost;
546 static dd_colrange d_last=0;
547 static dd_colset tieset,stieset; /* store the column indices with tie */
548
549 dd_init(val); dd_init(valn); dd_init(minval); dd_init(rat); dd_init(minrat);
550 if (d_last<d_size) {
551 if (d_last>0) {
552 for (j=1; j<=d_last; j++){ dd_clear(rcost[j-1]);}
553 free(rcost);
554 set_free(tieset);
555 set_free(stieset);
556 }
557 rcost=(mytype*) calloc(d_size,sizeof(mytype));
558 for (j=1; j<=d_size; j++){ dd_init(rcost[j-1]);}
559 set_initialize(&tieset,d_size);
560 set_initialize(&stieset,d_size);
561 d_last=d_size;
562 }
563
564 *r=0; *s=0;
565 *selected=dd_FALSE;
566 *lps=dd_LPSundecided;
567 for (j=1; j<=d_size; j++){
568 if (j!=rhscol){
569 dd_TableauEntry(&(rcost[j-1]),m_size,d_size,A,T,objrow,j);
570 if (dd_Positive(rcost[j-1])) {
571 dualfeasible=dd_FALSE;
572 }
573 }
574 }
575 if (dualfeasible){
576 while ((*lps==dd_LPSundecided) && (!rowselected) && (!colselected)) {
577 for (i=1; i<=m_size; i++) {
578 if (i!=objrow && bflag[i]==-1) { /* i is a basic variable */
579 if (Phase1){
580 dd_TableauEntry(&val, m_size,d_size,A,T,i,bflag[m_size]);
581 dd_neg(val,val);
582 /* for dual Phase I. The RHS (dual objective) is the negative of the auxiliary variable column. */
583 }
584 else {dd_TableauEntry(&val,m_size,d_size,A,T,i,rhscol);}
585 if (dd_Smaller(val,minval)) {
586 *r=i;
587 dd_set(minval,val);
588 }
589 }
590 }
591 if (dd_Nonnegative(minval)) {
592 *lps=dd_Optimal;
593 }
594 else {
595 rowselected=dd_TRUE;
596 set_emptyset(tieset);
597 for (j=1; j<=d_size; j++){
598 dd_TableauEntry(&val,m_size,d_size,A,T,*r,j);
599 if (j!=rhscol && dd_Positive(val)) {
600 dd_div(rat,rcost[j-1],val);
601 dd_neg(rat,rat);
602 if (*s==0 || dd_Smaller(rat,minrat)){
603 dd_set(minrat,rat);
604 *s=j;
605 set_emptyset(tieset);
606 set_addelem(tieset, j);
607 } else if (dd_Equal(rat,minrat)){
608 set_addelem(tieset,j);
609 }
610 }
611 }
612 if (*s>0) {
613 if (!lexicopivot || set_card(tieset)==1){
614 colselected=dd_TRUE; *selected=dd_TRUE;
615 } else { /* lexicographic rule with respect to the given reference cobasis. */
616 if (localdebug) {printf("Tie occurred at:"); set_write(tieset); printf("\n");
617 dd_WriteTableau(stderr,m_size,d_size,A,T,nbindex,bflag);
618 }
619 *s=0;
620 k=2; /* k runs through the column indices except RHS. */
621 do {
622 iref=nbindex_ref[k]; /* iref runs though the reference basic indices */
623 if (iref>0) {
624 j=bflag[iref];
625 if (j>0) {
626 if (set_member(j,tieset) && set_card(tieset)==1) {
627 *s=j;
628 colselected=dd_TRUE;
629 } else {
630 set_delelem(tieset, j);
631 /* iref is cobasic, and the corresponding col is not the pivot column except it is the last one. */
632 }
633 } else {
634 *s=0;
635 for (j=1; j<=d_size; j++){
636 if (set_member(j,tieset)) {
637 dd_TableauEntry(&val,m_size,d_size,A,T,*r,j);
638 dd_TableauEntry(&valn,m_size,d_size,A,T,iref,j);
639 if (j!=rhscol && dd_Positive(val)) {
640 dd_div(rat,valn,val);
641 if (*s==0 || dd_Smaller(rat,minrat)){
642 dd_set(minrat,rat);
643 *s=j;
644 set_emptyset(stieset);
645 set_addelem(stieset, j);
646 } else if (dd_Equal(rat,minrat)){
647 set_addelem(stieset,j);
648 }
649 }
650 }
651 }
652 set_copy(tieset,stieset);
653 if (set_card(tieset)==1) colselected=dd_TRUE;
654 }
655 }
656 k+=1;
657 } while (!colselected && k<=d_size);
658 *selected=dd_TRUE;
659 }
660 } else *lps=dd_Inconsistent;
661 }
662 } /* end of while */
663 }
664 if (localdebug) {
665 if (Phase1) fprintf(stderr,"Phase 1 : select %ld,%ld\n",*r,*s);
666 else fprintf(stderr,"Phase 2 : select %ld,%ld\n",*r,*s);
667 }
668 dd_clear(val); dd_clear(valn); dd_clear(minval); dd_clear(rat); dd_clear(minrat);
669 }
670
dd_TableauEntry(mytype * x,dd_rowrange m_size,dd_colrange d_size,dd_Amatrix X,dd_Bmatrix T,dd_rowrange r,dd_colrange s)671 void dd_TableauEntry(mytype *x,dd_rowrange m_size, dd_colrange d_size, dd_Amatrix X, dd_Bmatrix T,
672 dd_rowrange r, dd_colrange s)
673 /* Compute the (r,s) entry of X.T */
674 {
675 dd_colrange j;
676 mytype temp;
677
678 dd_init(temp);
679 dd_set(*x,dd_purezero);
680 for (j=0; j< d_size; j++) {
681 dd_mul(temp,X[r-1][j], T[j][s-1]);
682 dd_add(*x, *x, temp);
683 }
684 dd_clear(temp);
685 }
686
dd_SelectPivot2(dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_Bmatrix T,dd_RowOrderType roworder,dd_rowindex ordervec,rowset equalityset,dd_rowrange rowmax,rowset NopivotRow,colset NopivotCol,dd_rowrange * r,dd_colrange * s,dd_boolean * selected)687 void dd_SelectPivot2(dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_Bmatrix T,
688 dd_RowOrderType roworder,dd_rowindex ordervec, rowset equalityset,
689 dd_rowrange rowmax,rowset NopivotRow,
690 colset NopivotCol,dd_rowrange *r,dd_colrange *s,
691 dd_boolean *selected)
692 /* Select a position (*r,*s) in the matrix A.T such that (A.T)[*r][*s] is nonzero
693 The choice is feasible, i.e., not on NopivotRow and NopivotCol, and
694 best with respect to the specified roworder
695 */
696 {
697 int stop;
698 dd_rowrange i,rtemp;
699 rowset rowexcluded;
700 mytype Xtemp;
701 dd_boolean localdebug=dd_FALSE;
702
703 stop = dd_FALSE;
704 localdebug=dd_debug;
705 dd_init(Xtemp);
706 set_initialize(&rowexcluded,m_size);
707 set_copy(rowexcluded,NopivotRow);
708 for (i=rowmax+1;i<=m_size;i++) {
709 set_addelem(rowexcluded,i); /* cannot pivot on any row > rmax */
710 }
711 *selected = dd_FALSE;
712 do {
713 rtemp=0; i=1;
714 while (i<=m_size && rtemp==0) { /* equalityset vars have highest priorities */
715 if (set_member(i,equalityset) && !set_member(i,rowexcluded)){
716 if (localdebug) fprintf(stderr,"marked set %ld chosen as a candidate\n",i);
717 rtemp=i;
718 }
719 i++;
720 }
721 if (rtemp==0) dd_SelectPreorderedNext2(m_size,d_size,rowexcluded,ordervec,&rtemp);;
722 if (rtemp>=1) {
723 *r=rtemp;
724 *s=1;
725 while (*s <= d_size && !*selected) {
726 dd_TableauEntry(&Xtemp,m_size,d_size,A,T,*r,*s);
727 if (!set_member(*s,NopivotCol) && dd_Nonzero(Xtemp)) {
728 *selected = dd_TRUE;
729 stop = dd_TRUE;
730 } else {
731 (*s)++;
732 }
733 }
734 if (!*selected) {
735 set_addelem(rowexcluded,rtemp);
736 }
737 }
738 else {
739 *r = 0;
740 *s = 0;
741 stop = dd_TRUE;
742 }
743 } while (!stop);
744 set_free(rowexcluded); dd_clear(Xtemp);
745 }
746
dd_GaussianColumnPivot(dd_rowrange m_size,dd_colrange d_size,dd_Amatrix X,dd_Bmatrix T,dd_rowrange r,dd_colrange s)747 void dd_GaussianColumnPivot(dd_rowrange m_size, dd_colrange d_size,
748 dd_Amatrix X, dd_Bmatrix T, dd_rowrange r, dd_colrange s)
749 /* Update the Transformation matrix T with the pivot operation on (r,s)
750 This procedure performs a implicit pivot operation on the matrix X by
751 updating the dual basis inverse T.
752 */
753 {
754 dd_colrange j, j1;
755 mytype Xtemp0, Xtemp1, Xtemp;
756 static dd_Arow Rtemp;
757 static dd_colrange last_d=0;
758
759 dd_init(Xtemp0); dd_init(Xtemp1); dd_init(Xtemp);
760 if (last_d!=d_size){
761 if (last_d>0) {
762 for (j=1; j<=last_d; j++) dd_clear(Rtemp[j-1]);
763 free(Rtemp);
764 }
765 Rtemp=(mytype*)calloc(d_size,sizeof(mytype));
766 for (j=1; j<=d_size; j++) dd_init(Rtemp[j-1]);
767 last_d=d_size;
768 }
769
770 for (j=1; j<=d_size; j++) {
771 dd_TableauEntry(&(Rtemp[j-1]), m_size, d_size, X, T, r,j);
772 }
773 dd_set(Xtemp0,Rtemp[s-1]);
774 for (j = 1; j <= d_size; j++) {
775 if (j != s) {
776 dd_div(Xtemp,Rtemp[j-1],Xtemp0);
777 dd_set(Xtemp1,dd_purezero);
778 for (j1 = 1; j1 <= d_size; j1++){
779 dd_mul(Xtemp1,Xtemp,T[j1-1][s - 1]);
780 dd_sub(T[j1-1][j-1],T[j1-1][j-1],Xtemp1);
781 /* T[j1-1][j-1] -= T[j1-1][s - 1] * Xtemp / Xtemp0; */
782 }
783 }
784 }
785 for (j = 1; j <= d_size; j++)
786 dd_div(T[j-1][s - 1],T[j-1][s - 1],Xtemp0);
787
788 dd_clear(Xtemp0); dd_clear(Xtemp1); dd_clear(Xtemp);
789 }
790
dd_GaussianColumnPivot2(dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_Bmatrix T,dd_colindex nbindex,dd_rowindex bflag,dd_rowrange r,dd_colrange s)791 void dd_GaussianColumnPivot2(dd_rowrange m_size,dd_colrange d_size,
792 dd_Amatrix A,dd_Bmatrix T,dd_colindex nbindex,dd_rowindex bflag,dd_rowrange r,dd_colrange s)
793 /* Update the Transformation matrix T with the pivot operation on (r,s)
794 This procedure performs a implicit pivot operation on the matrix A by
795 updating the dual basis inverse T.
796 */
797 {
798 int localdebug=dd_FALSE;
799 long entering;
800
801 if (dd_debug) localdebug=dd_debug;
802 dd_GaussianColumnPivot(m_size,d_size,A,T,r,s);
803 entering=nbindex[s];
804 bflag[r]=s; /* the nonbasic variable r corresponds to column s */
805 nbindex[s]=r; /* the nonbasic variable on s column is r */
806
807 if (entering>0) bflag[entering]=-1;
808 /* original variables have negative index and should not affect the row index */
809
810 if (localdebug) {
811 fprintf(stderr,"dd_GaussianColumnPivot2\n");
812 fprintf(stderr," pivot: (leaving, entering) = (%ld, %ld)\n", r,entering);
813 fprintf(stderr, " bflag[%ld] is set to %ld\n", r, s);
814 }
815 }
816
817
dd_ResetTableau(dd_rowrange m_size,dd_colrange d_size,dd_Bmatrix T,dd_colindex nbindex,dd_rowindex bflag,dd_rowrange objrow,dd_colrange rhscol)818 void dd_ResetTableau(dd_rowrange m_size,dd_colrange d_size,dd_Bmatrix T,
819 dd_colindex nbindex,dd_rowindex bflag,dd_rowrange objrow,dd_colrange rhscol)
820 {
821 dd_rowrange i;
822 dd_colrange j;
823
824 /* Initialize T and nbindex */
825 for (j=1; j<=d_size; j++) nbindex[j]=-j;
826 nbindex[rhscol]=0;
827 /* RHS is already in nonbasis and is considered to be associated
828 with the zero-th row of input. */
829 dd_SetToIdentity(d_size,T);
830
831 /* Set the bflag according to nbindex */
832 for (i=1; i<=m_size; i++) bflag[i]=-1;
833 /* all basic variables have index -1 */
834 bflag[objrow]= 0;
835 /* bflag of the objective variable is 0,
836 different from other basic variables which have -1 */
837 for (j=1; j<=d_size; j++) if (nbindex[j]>0) bflag[nbindex[j]]=j;
838 /* bflag of a nonbasic variable is its column number */
839
840 }
841
dd_SelectCrissCrossPivot(dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_Bmatrix T,dd_rowindex bflag,dd_rowrange objrow,dd_colrange rhscol,dd_rowrange * r,dd_colrange * s,int * selected,dd_LPStatusType * lps)842 void dd_SelectCrissCrossPivot(dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_Bmatrix T,
843 dd_rowindex bflag,dd_rowrange objrow,dd_colrange rhscol,
844 dd_rowrange *r,dd_colrange *s,
845 int *selected,dd_LPStatusType *lps)
846 {
847 int colselected=dd_FALSE,rowselected=dd_FALSE;
848 dd_rowrange i;
849 mytype val;
850
851 dd_init(val);
852 *selected=dd_FALSE;
853 *lps=dd_LPSundecided;
854 while ((*lps==dd_LPSundecided) && (!rowselected) && (!colselected)) {
855 for (i=1; i<=m_size; i++) {
856 if (i!=objrow && bflag[i]==-1) { /* i is a basic variable */
857 dd_TableauEntry(&val,m_size,d_size,A,T,i,rhscol);
858 if (dd_Negative(val)) {
859 rowselected=dd_TRUE;
860 *r=i;
861 break;
862 }
863 }
864 else if (bflag[i] >0) { /* i is nonbasic variable */
865 dd_TableauEntry(&val,m_size,d_size,A,T,objrow,bflag[i]);
866 if (dd_Positive(val)) {
867 colselected=dd_TRUE;
868 *s=bflag[i];
869 break;
870 }
871 }
872 }
873 if ((!rowselected) && (!colselected)) {
874 *lps=dd_Optimal;
875 return;
876 }
877 else if (rowselected) {
878 for (i=1; i<=m_size; i++) {
879 if (bflag[i] >0) { /* i is nonbasic variable */
880 dd_TableauEntry(&val,m_size,d_size,A,T,*r,bflag[i]);
881 if (dd_Positive(val)) {
882 colselected=dd_TRUE;
883 *s=bflag[i];
884 *selected=dd_TRUE;
885 break;
886 }
887 }
888 }
889 }
890 else if (colselected) {
891 for (i=1; i<=m_size; i++) {
892 if (i!=objrow && bflag[i]==-1) { /* i is a basic variable */
893 dd_TableauEntry(&val,m_size,d_size,A,T,i,*s);
894 if (dd_Negative(val)) {
895 rowselected=dd_TRUE;
896 *r=i;
897 *selected=dd_TRUE;
898 break;
899 }
900 }
901 }
902 }
903 if (!rowselected) {
904 *lps=dd_DualInconsistent;
905 }
906 else if (!colselected) {
907 *lps=dd_Inconsistent;
908 }
909 }
910 dd_clear(val);
911 }
912
dd_CrissCrossSolve(dd_LPPtr lp,dd_ErrorType * err)913 void dd_CrissCrossSolve(dd_LPPtr lp, dd_ErrorType *err)
914 {
915 switch (lp->objective) {
916 case dd_LPmax:
917 dd_CrissCrossMaximize(lp,err);
918 break;
919
920 case dd_LPmin:
921 dd_CrissCrossMinimize(lp,err);
922 break;
923
924 case dd_LPnone: *err=dd_NoLPObjective; break;
925 }
926
927 }
928
dd_DualSimplexSolve(dd_LPPtr lp,dd_ErrorType * err)929 void dd_DualSimplexSolve(dd_LPPtr lp, dd_ErrorType *err)
930 {
931 switch (lp->objective) {
932 case dd_LPmax:
933 dd_DualSimplexMaximize(lp,err);
934 break;
935
936 case dd_LPmin:
937 dd_DualSimplexMinimize(lp,err);
938 break;
939
940 case dd_LPnone: *err=dd_NoLPObjective; break;
941 }
942 }
943
944 #ifdef GMPRATIONAL
945
LPSf2LPS(ddf_LPStatusType lpsf)946 dd_LPStatusType LPSf2LPS(ddf_LPStatusType lpsf)
947 {
948 dd_LPStatusType lps=dd_LPSundecided;
949
950 switch (lpsf) {
951 case ddf_LPSundecided: lps=dd_LPSundecided; break;
952 case ddf_Optimal: lps=dd_Optimal; break;
953 case ddf_Inconsistent: lps=dd_Inconsistent; break;
954 case ddf_DualInconsistent: lps=dd_DualInconsistent; break;
955 case ddf_StrucInconsistent: lps=dd_StrucInconsistent; break;
956 case ddf_StrucDualInconsistent: lps=dd_StrucDualInconsistent; break;
957 case ddf_Unbounded: lps=dd_Unbounded; break;
958 case ddf_DualUnbounded: lps=dd_DualUnbounded; break;
959 }
960 return lps;
961 }
962
963
dd_BasisStatus(ddf_LPPtr lpf,dd_LPPtr lp,dd_boolean * LPScorrect)964 void dd_BasisStatus(ddf_LPPtr lpf, dd_LPPtr lp, dd_boolean *LPScorrect)
965 {
966 int i;
967 dd_colrange se, j;
968 dd_boolean basisfound;
969
970 switch (lp->objective) {
971 case dd_LPmax:
972 dd_BasisStatusMaximize(lp->m,lp->d,lp->A,lp->B,lp->equalityset,lp->objrow,lp->rhscol,
973 lpf->LPS,&(lp->optvalue),lp->sol,lp->dsol,lp->posset_extra,lpf->nbindex,lpf->re,lpf->se,&se,lp->pivots,
974 &basisfound, LPScorrect);
975 if (*LPScorrect) {
976 /* printf("BasisStatus Check: the current basis is verified with GMP\n"); */
977 lp->LPS=LPSf2LPS(lpf->LPS);
978 lp->re=lpf->re;
979 lp->se=se;
980 for (j=1; j<=lp->d; j++) lp->nbindex[j]=lpf->nbindex[j];
981 }
982 for (i=1; i<=5; i++) lp->pivots[i-1]+=lpf->pivots[i-1];
983 break;
984 case dd_LPmin:
985 dd_BasisStatusMinimize(lp->m,lp->d,lp->A,lp->B,lp->equalityset,lp->objrow,lp->rhscol,
986 lpf->LPS,&(lp->optvalue),lp->sol,lp->dsol,lp->posset_extra,lpf->nbindex,lpf->re,lpf->se,&se,lp->pivots,
987 &basisfound, LPScorrect);
988 if (*LPScorrect) {
989 /* printf("BasisStatus Check: the current basis is verified with GMP\n"); */
990 lp->LPS=LPSf2LPS(lpf->LPS);
991 lp->re=lpf->re;
992 lp->se=se;
993 for (j=1; j<=lp->d; j++) lp->nbindex[j]=lpf->nbindex[j];
994 }
995 for (i=1; i<=5; i++) lp->pivots[i-1]+=lpf->pivots[i-1];
996 break;
997 case dd_LPnone: break;
998 }
999 }
1000 #endif
1001
dd_FindLPBasis(dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_Bmatrix T,dd_rowindex OV,dd_rowset equalityset,dd_colindex nbindex,dd_rowindex bflag,dd_rowrange objrow,dd_colrange rhscol,dd_colrange * cs,int * found,dd_LPStatusType * lps,long * pivot_no)1002 void dd_FindLPBasis(dd_rowrange m_size,dd_colrange d_size,
1003 dd_Amatrix A, dd_Bmatrix T,dd_rowindex OV,dd_rowset equalityset, dd_colindex nbindex,
1004 dd_rowindex bflag,dd_rowrange objrow,dd_colrange rhscol,
1005 dd_colrange *cs,int *found,dd_LPStatusType *lps,long *pivot_no)
1006 {
1007 /* Find a LP basis using Gaussian pivots.
1008 If the problem has an LP basis,
1009 the procedure returns *found=dd_TRUE,*lps=LPSundecided and an LP basis.
1010 If the constraint matrix A (excluding the rhs and objective) is not
1011 column independent, there are two cases. If the dependency gives a dual
1012 inconsistency, this returns *found=dd_FALSE, *lps=dd_StrucDualInconsistent and
1013 the evidence column *s. Otherwise, this returns *found=dd_TRUE,
1014 *lps=LPSundecided and an LP basis of size less than d_size. Columns j
1015 that do not belong to the basis (i.e. cannot be chosen as pivot because
1016 they are all zero) will be indicated in nbindex vector: nbindex[j] will
1017 be negative and set to -j.
1018 */
1019 int chosen,stop;
1020 long pivots_p0=0,rank;
1021 colset ColSelected;
1022 rowset RowSelected;
1023 mytype val;
1024
1025 dd_rowrange r;
1026 dd_colrange j,s;
1027
1028 dd_init(val);
1029 *found=dd_FALSE; *cs=0; rank=0;
1030 stop=dd_FALSE;
1031 *lps=dd_LPSundecided;
1032
1033 set_initialize(&RowSelected,m_size);
1034 set_initialize(&ColSelected,d_size);
1035 set_addelem(RowSelected,objrow);
1036 set_addelem(ColSelected,rhscol);
1037
1038 stop=dd_FALSE;
1039 do { /* Find a LP basis */
1040 dd_SelectPivot2(m_size,d_size,A,T,dd_MinIndex,OV,equalityset,
1041 m_size,RowSelected,ColSelected,&r,&s,&chosen);
1042 if (chosen) {
1043 set_addelem(RowSelected,r);
1044 set_addelem(ColSelected,s);
1045 dd_GaussianColumnPivot2(m_size,d_size,A,T,nbindex,bflag,r,s);
1046 pivots_p0++;
1047 rank++;
1048 } else {
1049 for (j=1;j<=d_size && *lps==dd_LPSundecided; j++) {
1050 if (j!=rhscol && nbindex[j]<0){
1051 dd_TableauEntry(&val,m_size,d_size,A,T,objrow,j);
1052 if (dd_Nonzero(val)){ /* dual inconsistent */
1053 *lps=dd_StrucDualInconsistent;
1054 *cs=j;
1055 /* dual inconsistent because the nonzero reduced cost */
1056 }
1057 }
1058 }
1059 if (*lps==dd_LPSundecided) *found=dd_TRUE;
1060 /* dependent columns but not dual inconsistent. */
1061 stop=dd_TRUE;
1062 }
1063 /* printf("d_size=%ld, rank=%ld\n",d_size,rank); */
1064 if (rank==d_size-1) {
1065 stop = dd_TRUE;
1066 *found=dd_TRUE;
1067 }
1068 } while (!stop);
1069
1070 *pivot_no=pivots_p0;
1071 dd_statBApivots+=pivots_p0;
1072 set_free(RowSelected);
1073 set_free(ColSelected);
1074 dd_clear(val);
1075 }
1076
1077
dd_FindLPBasis2(dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_Bmatrix T,dd_rowindex OV,dd_rowset equalityset,dd_colindex nbindex,dd_rowindex bflag,dd_rowrange objrow,dd_colrange rhscol,dd_colrange * cs,int * found,long * pivot_no)1078 void dd_FindLPBasis2(dd_rowrange m_size,dd_colrange d_size,
1079 dd_Amatrix A, dd_Bmatrix T,dd_rowindex OV,dd_rowset equalityset, dd_colindex nbindex,
1080 dd_rowindex bflag,dd_rowrange objrow,dd_colrange rhscol,
1081 dd_colrange *cs,int *found,long *pivot_no)
1082 {
1083 /* Similar to dd_FindLPBasis but it is much simpler. This tries to recompute T for
1084 the specified basis given by nbindex. It will return *found=dd_FALSE if the specified
1085 basis is not a basis.
1086 */
1087 int chosen,stop;
1088 long pivots_p0=0,rank;
1089 dd_colset ColSelected,DependentCols;
1090 dd_rowset RowSelected, NopivotRow;
1091 mytype val;
1092 dd_boolean localdebug=dd_FALSE;
1093
1094 dd_rowrange r,negcount=0;
1095 dd_colrange j,s;
1096
1097 dd_init(val);
1098 *found=dd_FALSE; *cs=0; rank=0;
1099
1100 set_initialize(&RowSelected,m_size);
1101 set_initialize(&DependentCols,d_size);
1102 set_initialize(&ColSelected,d_size);
1103 set_initialize(&NopivotRow,m_size);
1104 set_addelem(RowSelected,objrow);
1105 set_addelem(ColSelected,rhscol);
1106 set_compl(NopivotRow, NopivotRow); /* set NopivotRow to be the groundset */
1107
1108 for (j=2; j<=d_size; j++)
1109 if (nbindex[j]>0)
1110 set_delelem(NopivotRow, nbindex[j]);
1111 else if (nbindex[j]<0){
1112 negcount++;
1113 set_addelem(DependentCols, -nbindex[j]);
1114 set_addelem(ColSelected, -nbindex[j]);
1115 }
1116
1117 set_uni(RowSelected, RowSelected, NopivotRow); /* RowSelected is the set of rows not allowed to poviot on */
1118
1119 stop=dd_FALSE;
1120 do { /* Find a LP basis */
1121 dd_SelectPivot2(m_size,d_size,A,T,dd_MinIndex,OV,equalityset, m_size,RowSelected,ColSelected,&r,&s,&chosen);
1122 if (chosen) {
1123 set_addelem(RowSelected,r);
1124 set_addelem(ColSelected,s);
1125
1126 dd_GaussianColumnPivot2(m_size,d_size,A,T,nbindex,bflag,r,s);
1127 if (localdebug && m_size <=10){
1128 dd_WriteBmatrix(stderr,d_size,T);
1129 dd_WriteTableau(stderr,m_size,d_size,A,T,nbindex,bflag);
1130 }
1131 pivots_p0++;
1132 rank++;
1133 } else{
1134 *found=dd_FALSE; /* cannot pivot on any of the spacified positions. */
1135 stop=dd_TRUE;
1136 }
1137 if (rank==d_size-1-negcount) {
1138 if (negcount){
1139 /* Now it tries to pivot on rows that are supposed to be dependent. */
1140 set_diff(ColSelected, ColSelected, DependentCols);
1141 dd_SelectPivot2(m_size,d_size,A,T,dd_MinIndex,OV,equalityset, m_size,RowSelected,ColSelected,&r,&s,&chosen);
1142 if (chosen) *found=dd_FALSE; /* not supposed to be independent */
1143 else *found=dd_TRUE;
1144 if (localdebug){
1145 printf("Try to check the dependent cols:");
1146 set_write(DependentCols);
1147 if (chosen) printf("They are not dependent. Can still pivot on (%ld, %ld)\n",r, s);
1148 else printf("They are indeed dependent.\n");
1149 }
1150 } else {
1151 *found=dd_TRUE;
1152 }
1153 stop = dd_TRUE;
1154 }
1155 } while (!stop);
1156
1157 for (j=1; j<=d_size; j++) if (nbindex[j]>0) bflag[nbindex[j]]=j;
1158 *pivot_no=pivots_p0;
1159 set_free(RowSelected);
1160 set_free(ColSelected);
1161 set_free(NopivotRow);
1162 set_free(DependentCols);
1163 dd_clear(val);
1164 }
1165
dd_FindDualFeasibleBasis(dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_Bmatrix T,dd_rowindex OV,dd_colindex nbindex,dd_rowindex bflag,dd_rowrange objrow,dd_colrange rhscol,dd_boolean lexicopivot,dd_colrange * s,dd_ErrorType * err,dd_LPStatusType * lps,long * pivot_no,long maxpivots)1166 void dd_FindDualFeasibleBasis(dd_rowrange m_size,dd_colrange d_size,
1167 dd_Amatrix A,dd_Bmatrix T,dd_rowindex OV,
1168 dd_colindex nbindex,dd_rowindex bflag,dd_rowrange objrow,dd_colrange rhscol, dd_boolean lexicopivot,
1169 dd_colrange *s,dd_ErrorType *err,dd_LPStatusType *lps,long *pivot_no, long maxpivots)
1170 {
1171 /* Find a dual feasible basis using Phase I of Dual Simplex method.
1172 If the problem is dual feasible,
1173 the procedure returns *err=NoError, *lps=LPSundecided and a dual feasible
1174 basis. If the problem is dual infeasible, this returns
1175 *err=NoError, *lps=DualInconsistent and the evidence column *s.
1176 Caution: matrix A must have at least one extra row: the row space A[m_size] must
1177 have been allocated.
1178 */
1179 dd_boolean phase1,dualfeasible=dd_TRUE;
1180 dd_boolean localdebug=dd_FALSE,chosen,stop;
1181 dd_LPStatusType LPSphase1;
1182 long pivots_p1=0;
1183 dd_rowrange i,r_val;
1184 dd_colrange j,l,ms=0,s_val,local_m_size;
1185 mytype x,val,maxcost,axvalue,maxratio;
1186 static dd_colrange d_last=0;
1187 static dd_Arow rcost;
1188 static dd_colindex nbindex_ref; /* to be used to store the initial feasible basis for lexico rule */
1189
1190 mytype scaling,svalue; /* random scaling mytype value */
1191 mytype minval;
1192
1193 if (dd_debug) localdebug=dd_debug;
1194 dd_init(x); dd_init(val); dd_init(scaling); dd_init(svalue); dd_init(axvalue);
1195 dd_init(maxcost); dd_set(maxcost,dd_minuszero);
1196 dd_init(maxratio); dd_set(maxratio,dd_minuszero);
1197 if (d_last<d_size) {
1198 if (d_last>0) {
1199 for (j=1; j<=d_last; j++){ dd_clear(rcost[j-1]);}
1200 free(rcost);
1201 free(nbindex_ref);
1202 }
1203 rcost=(mytype*) calloc(d_size,sizeof(mytype));
1204 nbindex_ref=(long*) calloc(d_size+1,sizeof(long));
1205 for (j=1; j<=d_size; j++){ dd_init(rcost[j-1]);}
1206 d_last=d_size;
1207 }
1208
1209 *err=dd_NoError; *lps=dd_LPSundecided; *s=0;
1210 local_m_size=m_size+1; /* increase m_size by 1 */
1211
1212 ms=0; /* ms will be the index of column which has the largest reduced cost */
1213 for (j=1; j<=d_size; j++){
1214 if (j!=rhscol){
1215 dd_TableauEntry(&(rcost[j-1]),local_m_size,d_size,A,T,objrow,j);
1216 if (dd_Larger(rcost[j-1],maxcost)) {dd_set(maxcost,rcost[j-1]); ms = j;}
1217 }
1218 }
1219 if (dd_Positive(maxcost)) dualfeasible=dd_FALSE;
1220
1221 if (!dualfeasible){
1222 for (j=1; j<=d_size; j++){
1223 dd_set(A[local_m_size-1][j-1], dd_purezero);
1224 for (l=1; l<=d_size; l++){
1225 if (nbindex[l]>0) {
1226 dd_set_si(scaling,l+10);
1227 dd_mul(svalue,A[nbindex[l]-1][j-1],scaling);
1228 dd_sub(A[local_m_size-1][j-1],A[local_m_size-1][j-1],svalue);
1229 /* To make the auxiliary row (0,-11,-12,...,-d-10).
1230 It is likely to be better than (0, -1, -1, ..., -1)
1231 to avoid a degenerate LP. Version 093c. */
1232 }
1233 }
1234 }
1235
1236 if (localdebug){
1237 fprintf(stderr,"\ndd_FindDualFeasibleBasis: curruent basis is not dual feasible.\n");
1238 fprintf(stderr,"because of the column %ld assoc. with var %ld dual cost =",
1239 ms,nbindex[ms]);
1240 dd_WriteNumber(stderr, maxcost);
1241 if (localdebug) {
1242 if (m_size <=100 && d_size <=30){
1243 printf("\ndd_FindDualFeasibleBasis: the starting dictionary.\n");
1244 dd_WriteTableau(stdout,m_size+1,d_size,A,T,nbindex,bflag);
1245 }
1246 }
1247 }
1248
1249 ms=0;
1250 /* Ratio Test: ms will be now the index of column which has the largest reduced cost
1251 over the auxiliary row entry */
1252 for (j=1; j<=d_size; j++){
1253 if ((j!=rhscol) && dd_Positive(rcost[j-1])){
1254 dd_TableauEntry(&axvalue,local_m_size,d_size,A,T,local_m_size,j);
1255 if (dd_Nonnegative(axvalue)) {
1256 *err=dd_NumericallyInconsistent;
1257 /* This should not happen as they are set negative above. Quit the phase I.*/
1258 if (localdebug) fprintf(stderr,"dd_FindDualFeasibleBasis: Numerical Inconsistency detected.\n");
1259 goto _L99;
1260 }
1261 dd_neg(axvalue,axvalue);
1262 dd_div(axvalue,rcost[j-1],axvalue); /* axvalue is the negative of ratio that is to be maximized. */
1263 if (dd_Larger(axvalue,maxratio)) {
1264 dd_set(maxratio,axvalue);
1265 ms = j;
1266 }
1267 }
1268 }
1269
1270 if (ms==0) {
1271 *err=dd_NumericallyInconsistent; /* This should not happen. Quit the phase I.*/
1272 if (localdebug) fprintf(stderr,"dd_FindDualFeasibleBasis: Numerical Inconsistency detected.\n");
1273 goto _L99;
1274 }
1275
1276 /* Pivot on (local_m_size,ms) so that the dual basic solution becomes feasible */
1277 dd_GaussianColumnPivot2(local_m_size,d_size,A,T,nbindex,bflag,local_m_size,ms);
1278 pivots_p1=pivots_p1+1;
1279 if (localdebug) {
1280 printf("\ndd_FindDualFeasibleBasis: Pivot on %ld %ld.\n",local_m_size,ms);
1281 }
1282
1283 for (j=1; j<=d_size; j++) nbindex_ref[j]=nbindex[j];
1284 /* set the reference basis to be the current feasible basis. */
1285 if (localdebug){
1286 fprintf(stderr, "Store the current feasible basis:");
1287 for (j=1; j<=d_size; j++) fprintf(stderr, " %ld", nbindex_ref[j]);
1288 fprintf(stderr, "\n");
1289 if (m_size <=100 && d_size <=30)
1290 dd_WriteSignTableau2(stdout,m_size+1,d_size,A,T,nbindex_ref,nbindex,bflag);
1291 }
1292
1293 phase1=dd_TRUE; stop=dd_FALSE;
1294 do { /* Dual Simplex Phase I */
1295 chosen=dd_FALSE; LPSphase1=dd_LPSundecided;
1296 if (pivots_p1>maxpivots) {
1297 *err=dd_LPCycling;
1298 fprintf(stderr,"max number %ld of pivots performed in Phase I. Switch to the anticycling phase.\n", maxpivots);
1299 goto _L99; /* failure due to max no. of pivots performed */
1300 }
1301 dd_SelectDualSimplexPivot(local_m_size,d_size,phase1,A,T,OV,nbindex_ref,nbindex,bflag,
1302 objrow,rhscol,lexicopivot,&r_val,&s_val,&chosen,&LPSphase1);
1303 if (!chosen) {
1304 /* The current dictionary is terminal. There are two cases:
1305 dd_TableauEntry(local_m_size,d_size,A,T,objrow,ms) is negative or zero.
1306 The first case implies dual infeasible,
1307 and the latter implies dual feasible but local_m_size is still in nonbasis.
1308 We must pivot in the auxiliary variable local_m_size.
1309 */
1310 dd_TableauEntry(&x,local_m_size,d_size,A,T,objrow,ms);
1311 if (dd_Negative(x)){
1312 *err=dd_NoError; *lps=dd_DualInconsistent; *s=ms;
1313 }
1314 if (localdebug) {
1315 fprintf(stderr,"\ndd_FindDualFeasibleBasis: the auxiliary variable was forced to enter the basis (# pivots = %ld).\n",pivots_p1);
1316 fprintf(stderr," -- objrow %ld, ms %ld entry: ",objrow,ms);
1317 dd_WriteNumber(stderr, x); fprintf(stderr,"\n");
1318 if (dd_Negative(x)){
1319 fprintf(stderr,"->The basis is dual inconsistent. Terminate.\n");
1320 } else {
1321 fprintf(stderr,"->The basis is feasible. Go to phase II.\n");
1322 }
1323 }
1324
1325 dd_init(minval);
1326 r_val=0;
1327 for (i=1; i<=local_m_size; i++){
1328 if (bflag[i]<0) {
1329 /* i is basic and not the objective variable */
1330 dd_TableauEntry(&val,local_m_size,d_size,A,T,i,ms); /* auxiliary column*/
1331 if (dd_Smaller(val, minval)) {
1332 r_val=i;
1333 dd_set(minval,val);
1334 }
1335 }
1336 }
1337 dd_clear(minval);
1338
1339 if (r_val==0) {
1340 *err=dd_NumericallyInconsistent; /* This should not happen. Quit the phase I.*/
1341 if (localdebug) fprintf(stderr,"dd_FindDualFeasibleBasis: Numerical Inconsistency detected (r_val is 0).\n");
1342 goto _L99;
1343 }
1344
1345 dd_GaussianColumnPivot2(local_m_size,d_size,A,T,nbindex,bflag,r_val,ms);
1346 pivots_p1=pivots_p1+1;
1347 if (localdebug) {
1348 printf("\ndd_FindDualFeasibleBasis: make the %ld-th pivot on %ld %ld to force the auxiliary variable to enter the basis.\n",pivots_p1,r_val,ms);
1349 if (m_size <=100 && d_size <=30)
1350 dd_WriteSignTableau2(stdout,m_size+1,d_size,A,T,nbindex_ref,nbindex,bflag);
1351 }
1352
1353 stop=dd_TRUE;
1354
1355 } else {
1356 dd_GaussianColumnPivot2(local_m_size,d_size,A,T,nbindex,bflag,r_val,s_val);
1357 pivots_p1=pivots_p1+1;
1358 if (localdebug) {
1359 printf("\ndd_FindDualFeasibleBasis: make a %ld-th pivot on %ld %ld\n",pivots_p1,r_val,s_val);
1360 if (m_size <=100 && d_size <=30)
1361 dd_WriteSignTableau2(stdout,local_m_size,d_size,A,T,nbindex_ref,nbindex,bflag);
1362 }
1363
1364
1365 if (bflag[local_m_size]<0) {
1366 stop=dd_TRUE;
1367 if (localdebug)
1368 fprintf(stderr,"\nDualSimplex Phase I: the auxiliary variable entered the basis (# pivots = %ld).\nGo to phase II\n",pivots_p1);
1369 }
1370 }
1371 } while(!stop);
1372 }
1373 _L99:
1374 *pivot_no=pivots_p1;
1375 dd_statDS1pivots+=pivots_p1;
1376 dd_clear(x); dd_clear(val); dd_clear(maxcost); dd_clear(maxratio);
1377 dd_clear(scaling); dd_clear(svalue); dd_clear(axvalue);
1378 }
1379
dd_DualSimplexMinimize(dd_LPPtr lp,dd_ErrorType * err)1380 void dd_DualSimplexMinimize(dd_LPPtr lp,dd_ErrorType *err)
1381 {
1382 dd_colrange j;
1383
1384 *err=dd_NoError;
1385 for (j=1; j<=lp->d; j++)
1386 dd_neg(lp->A[lp->objrow-1][j-1],lp->A[lp->objrow-1][j-1]);
1387 dd_DualSimplexMaximize(lp,err);
1388 dd_neg(lp->optvalue,lp->optvalue);
1389 for (j=1; j<=lp->d; j++){
1390 if (lp->LPS!=dd_Inconsistent) {
1391 /* Inconsistent certificate stays valid for minimization, 0.94e */
1392 dd_neg(lp->dsol[j-1],lp->dsol[j-1]);
1393 }
1394 dd_neg(lp->A[lp->objrow-1][j-1],lp->A[lp->objrow-1][j-1]);
1395 }
1396 }
1397
dd_DualSimplexMaximize(dd_LPPtr lp,dd_ErrorType * err)1398 void dd_DualSimplexMaximize(dd_LPPtr lp,dd_ErrorType *err)
1399 /*
1400 When LP is inconsistent then lp->re returns the evidence row.
1401 When LP is dual-inconsistent then lp->se returns the evidence column.
1402 */
1403 {
1404 int stop,chosen,phase1,found;
1405 long pivots_ds=0,pivots_p0=0,pivots_p1=0,pivots_pc=0,maxpivots,maxpivfactor=20;
1406 dd_boolean localdebug=dd_FALSE,localdebug1=dd_FALSE;
1407
1408 #if !defined GMPRATIONAL
1409 long maxccpivots,maxccpivfactor=100;
1410 /* criss-cross should not cycle, but with floating-point arithmetics, it happens
1411 (very rarely). Jorg Rambau reported such an LP, in August 2003. Thanks Jorg!
1412 */
1413 #endif
1414
1415 dd_rowrange i,r;
1416 dd_colrange j,s;
1417 static dd_rowindex bflag;
1418 static long mlast=0,nlast=0;
1419 static dd_rowindex OrderVector; /* the permutation vector to store a preordered row indeces */
1420 static dd_colindex nbindex_ref; /* to be used to store the initial feasible basis for lexico rule */
1421
1422 double redpercent=0,redpercent_prev=0,redgain=0;
1423 unsigned int rseed=1;
1424
1425 /* *err=dd_NoError; */
1426 if (dd_debug) localdebug=dd_debug;
1427 set_emptyset(lp->redset_extra);
1428 for (i=0; i<= 4; i++) lp->pivots[i]=0;
1429 maxpivots=maxpivfactor*lp->d; /* maximum pivots to be performed before cc pivot is applied. */
1430 #if !defined GMPRATIONAL
1431 maxccpivots=maxccpivfactor*lp->d; /* maximum pivots to be performed with emergency cc pivots. */
1432 #endif
1433 if (mlast!=lp->m || nlast!=lp->d){
1434 if (mlast>0) { /* called previously with different lp->m */
1435 free(OrderVector);
1436 free(bflag);
1437 free(nbindex_ref);
1438 }
1439 OrderVector=(long *)calloc(lp->m+1,sizeof(*OrderVector));
1440 bflag=(long *) calloc(lp->m+2,sizeof(*bflag)); /* one more element for an auxiliary variable */
1441 nbindex_ref=(long*) calloc(lp->d+1,sizeof(long));
1442 mlast=lp->m;nlast=lp->d;
1443 }
1444 /* Initializing control variables. */
1445 dd_ComputeRowOrderVector2(lp->m,lp->d,lp->A,OrderVector,dd_MinIndex,rseed);
1446
1447 lp->re=0; lp->se=0;
1448
1449 dd_ResetTableau(lp->m,lp->d,lp->B,lp->nbindex,bflag,lp->objrow,lp->rhscol);
1450
1451 dd_FindLPBasis(lp->m,lp->d,lp->A,lp->B,OrderVector,lp->equalityset,lp->nbindex,bflag,
1452 lp->objrow,lp->rhscol,&s,&found,&(lp->LPS),&pivots_p0);
1453 lp->pivots[0]=pivots_p0;
1454
1455 if (!found){
1456 lp->se=s;
1457 goto _L99;
1458 /* No LP basis is found, and thus Inconsistent.
1459 Output the evidence column. */
1460 }
1461
1462 dd_FindDualFeasibleBasis(lp->m,lp->d,lp->A,lp->B,OrderVector,lp->nbindex,bflag,
1463 lp->objrow,lp->rhscol,lp->lexicopivot,&s, err,&(lp->LPS),&pivots_p1, maxpivots);
1464 lp->pivots[1]=pivots_p1;
1465
1466 for (j=1; j<=lp->d; j++) nbindex_ref[j]=lp->nbindex[j];
1467 /* set the reference basis to be the current feasible basis. */
1468 if (localdebug){
1469 fprintf(stderr, "dd_DualSimplexMaximize: Store the current feasible basis:");
1470 for (j=1; j<=lp->d; j++) fprintf(stderr, " %ld", nbindex_ref[j]);
1471 fprintf(stderr, "\n");
1472 if (lp->m <=100 && lp->d <=30)
1473 dd_WriteSignTableau2(stdout,lp->m+1,lp->d,lp->A,lp->B,nbindex_ref,lp->nbindex,bflag);
1474 }
1475
1476 if (*err==dd_LPCycling || *err==dd_NumericallyInconsistent){
1477 if (localdebug) fprintf(stderr, "Phase I failed and thus switch to the Criss-Cross method\n");
1478 dd_CrissCrossMaximize(lp,err);
1479 return;
1480 }
1481
1482 if (lp->LPS==dd_DualInconsistent){
1483 lp->se=s;
1484 goto _L99;
1485 /* No dual feasible basis is found, and thus DualInconsistent.
1486 Output the evidence column. */
1487 }
1488
1489 /* Dual Simplex Method */
1490 stop=dd_FALSE;
1491 do {
1492 chosen=dd_FALSE; lp->LPS=dd_LPSundecided; phase1=dd_FALSE;
1493 if (pivots_ds<maxpivots) {
1494 dd_SelectDualSimplexPivot(lp->m,lp->d,phase1,lp->A,lp->B,OrderVector,nbindex_ref,lp->nbindex,bflag,
1495 lp->objrow,lp->rhscol,lp->lexicopivot,&r,&s,&chosen,&(lp->LPS));
1496 }
1497 if (chosen) {
1498 pivots_ds=pivots_ds+1;
1499 if (lp->redcheck_extensive) {
1500 dd_GetRedundancyInformation(lp->m,lp->d,lp->A,lp->B,lp->nbindex, bflag, lp->redset_extra);
1501 set_uni(lp->redset_accum, lp->redset_accum,lp->redset_extra);
1502 redpercent=100*(double)set_card(lp->redset_extra)/(double)lp->m;
1503 redgain=redpercent-redpercent_prev;
1504 redpercent_prev=redpercent;
1505 if (localdebug1){
1506 fprintf(stderr,"\ndd_DualSimplexMaximize: Phase II pivot %ld on (%ld, %ld).\n",pivots_ds,r,s);
1507 fprintf(stderr," redundancy %f percent: redset size = %ld\n",redpercent,set_card(lp->redset_extra));
1508 }
1509 }
1510 }
1511 if (!chosen && lp->LPS==dd_LPSundecided) {
1512 if (localdebug1){
1513 fprintf(stderr,"Warning: an emergency CC pivot in Phase II is performed\n");
1514 /* In principle this should not be executed because we already have dual feasibility
1515 attained and dual simplex pivot should have been chosen. This might occur
1516 under floating point computation, or the case of cycling.
1517 */
1518 if (localdebug && lp->m <=100 && lp->d <=30){
1519 fprintf(stderr,"\ndd_DualSimplexMaximize: The current dictionary.\n");
1520 dd_WriteSignTableau2(stdout,lp->m,lp->d,lp->A,lp->B,nbindex_ref,lp->nbindex,bflag);
1521 }
1522 }
1523
1524 #if !defined GMPRATIONAL
1525 if (pivots_pc>maxccpivots) {
1526 *err=dd_LPCycling;
1527 stop=dd_TRUE;
1528 goto _L99;
1529 }
1530 #endif
1531
1532 dd_SelectCrissCrossPivot(lp->m,lp->d,lp->A,lp->B,bflag,
1533 lp->objrow,lp->rhscol,&r,&s,&chosen,&(lp->LPS));
1534 if (chosen) pivots_pc=pivots_pc+1;
1535 }
1536 if (chosen) {
1537 dd_GaussianColumnPivot2(lp->m,lp->d,lp->A,lp->B,lp->nbindex,bflag,r,s);
1538 if (localdebug && lp->m <=100 && lp->d <=30){
1539 fprintf(stderr,"\ndd_DualSimplexMaximize: The current dictionary.\n");
1540 dd_WriteSignTableau2(stdout,lp->m,lp->d,lp->A,lp->B,nbindex_ref,lp->nbindex,bflag);
1541 }
1542 } else {
1543 switch (lp->LPS){
1544 case dd_Inconsistent: lp->re=r;
1545 case dd_DualInconsistent: lp->se=s;
1546
1547 default: break;
1548 }
1549 stop=dd_TRUE;
1550 }
1551 } while(!stop);
1552 _L99:
1553 lp->pivots[2]=pivots_ds;
1554 lp->pivots[3]=pivots_pc;
1555 dd_statDS2pivots+=pivots_ds;
1556 dd_statACpivots+=pivots_pc;
1557
1558 dd_SetSolutions(lp->m,lp->d,lp->A,lp->B,lp->objrow,lp->rhscol,lp->LPS,&(lp->optvalue),lp->sol,lp->dsol,lp->posset_extra,lp->nbindex,lp->re,lp->se,bflag);
1559
1560 }
1561
1562
1563
dd_CrissCrossMinimize(dd_LPPtr lp,dd_ErrorType * err)1564 void dd_CrissCrossMinimize(dd_LPPtr lp,dd_ErrorType *err)
1565 {
1566 dd_colrange j;
1567
1568 *err=dd_NoError;
1569 for (j=1; j<=lp->d; j++)
1570 dd_neg(lp->A[lp->objrow-1][j-1],lp->A[lp->objrow-1][j-1]);
1571 dd_CrissCrossMaximize(lp,err);
1572 dd_neg(lp->optvalue,lp->optvalue);
1573 for (j=1; j<=lp->d; j++){
1574 if (lp->LPS!=dd_Inconsistent) {
1575 /* Inconsistent certificate stays valid for minimization, 0.94e */
1576 dd_neg(lp->dsol[j-1],lp->dsol[j-1]);
1577 }
1578 dd_neg(lp->A[lp->objrow-1][j-1],lp->A[lp->objrow-1][j-1]);
1579 }
1580 }
1581
dd_CrissCrossMaximize(dd_LPPtr lp,dd_ErrorType * err)1582 void dd_CrissCrossMaximize(dd_LPPtr lp,dd_ErrorType *err)
1583 /*
1584 When LP is inconsistent then lp->re returns the evidence row.
1585 When LP is dual-inconsistent then lp->se returns the evidence column.
1586 */
1587 {
1588 int stop,chosen,found;
1589 long pivots0,pivots1;
1590 #if !defined GMPRATIONAL
1591 long maxpivots,maxpivfactor=1000;
1592 /* criss-cross should not cycle, but with floating-point arithmetics, it happens
1593 (very rarely). Jorg Rambau reported such an LP, in August 2003. Thanks Jorg!
1594 */
1595 #endif
1596
1597 dd_rowrange i,r;
1598 dd_colrange s;
1599 static dd_rowindex bflag;
1600 static long mlast=0;
1601 static dd_rowindex OrderVector; /* the permutation vector to store a preordered row indeces */
1602 unsigned int rseed=1;
1603 dd_colindex nbtemp;
1604
1605 *err=dd_NoError;
1606 #if !defined GMPRATIONAL
1607 maxpivots=maxpivfactor*lp->d; /* maximum pivots to be performed when floating-point arithmetics is used. */
1608 #endif
1609 nbtemp=(long *) calloc(lp->d+1,sizeof(long));
1610 for (i=0; i<= 4; i++) lp->pivots[i]=0;
1611 if (bflag==NULL || mlast!=lp->m){
1612 if (mlast!=lp->m && mlast>0) {
1613 free(bflag); /* called previously with different lp->m */
1614 free(OrderVector);
1615 }
1616 bflag=(long *) calloc(lp->m+1,sizeof(long));
1617 OrderVector=(long *)calloc(lp->m+1,sizeof(long));
1618 /* initialize only for the first time or when a larger space is needed */
1619
1620 mlast=lp->m;
1621 }
1622 /* Initializing control variables. */
1623 dd_ComputeRowOrderVector2(lp->m,lp->d,lp->A,OrderVector,dd_MinIndex,rseed);
1624
1625 lp->re=0; lp->se=0; pivots1=0;
1626
1627 dd_ResetTableau(lp->m,lp->d,lp->B,lp->nbindex,bflag,lp->objrow,lp->rhscol);
1628
1629 dd_FindLPBasis(lp->m,lp->d,lp->A,lp->B,OrderVector,lp->equalityset,
1630 lp->nbindex,bflag,lp->objrow,lp->rhscol,&s,&found,&(lp->LPS),&pivots0);
1631 lp->pivots[0]+=pivots0;
1632
1633 if (!found){
1634 lp->se=s;
1635 goto _L99;
1636 /* No LP basis is found, and thus Inconsistent.
1637 Output the evidence column. */
1638 }
1639
1640 stop=dd_FALSE;
1641 do { /* Criss-Cross Method */
1642 #if !defined GMPRATIONAL
1643 if (pivots1>maxpivots) {
1644 *err=dd_LPCycling;
1645 fprintf(stderr,"max number %ld of pivots performed by the criss-cross method. Most likely due to the floating-point arithmetics error.\n", maxpivots);
1646 goto _L99; /* failure due to max no. of pivots performed */
1647 }
1648 #endif
1649
1650 dd_SelectCrissCrossPivot(lp->m,lp->d,lp->A,lp->B,bflag,
1651 lp->objrow,lp->rhscol,&r,&s,&chosen,&(lp->LPS));
1652 if (chosen) {
1653 dd_GaussianColumnPivot2(lp->m,lp->d,lp->A,lp->B,lp->nbindex,bflag,r,s);
1654 pivots1++;
1655 } else {
1656 switch (lp->LPS){
1657 case dd_Inconsistent: lp->re=r;
1658 case dd_DualInconsistent: lp->se=s;
1659
1660 default: break;
1661 }
1662 stop=dd_TRUE;
1663 }
1664 } while(!stop);
1665
1666 _L99:
1667 lp->pivots[1]+=pivots1;
1668 dd_statCCpivots+=pivots1;
1669 dd_SetSolutions(lp->m,lp->d,lp->A,lp->B,
1670 lp->objrow,lp->rhscol,lp->LPS,&(lp->optvalue),lp->sol,lp->dsol,lp->posset_extra,lp->nbindex,lp->re,lp->se,bflag);
1671 free(nbtemp);
1672 }
1673
dd_SetSolutions(dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_Bmatrix T,dd_rowrange objrow,dd_colrange rhscol,dd_LPStatusType LPS,mytype * optvalue,dd_Arow sol,dd_Arow dsol,dd_rowset posset,dd_colindex nbindex,dd_rowrange re,dd_colrange se,dd_rowindex bflag)1674 void dd_SetSolutions(dd_rowrange m_size,dd_colrange d_size,
1675 dd_Amatrix A,dd_Bmatrix T,
1676 dd_rowrange objrow,dd_colrange rhscol,dd_LPStatusType LPS,
1677 mytype *optvalue,dd_Arow sol,dd_Arow dsol,dd_rowset posset, dd_colindex nbindex,
1678 dd_rowrange re,dd_colrange se,dd_rowindex bflag)
1679 /*
1680 Assign the solution vectors to sol,dsol,*optvalue after solving
1681 the LP.
1682 */
1683 {
1684 dd_rowrange i;
1685 dd_colrange j;
1686 mytype x,sw;
1687 int localdebug=dd_FALSE;
1688
1689 dd_init(x); dd_init(sw);
1690 if (localdebug) fprintf(stderr,"SetSolutions:\n");
1691 switch (LPS){
1692 case dd_Optimal:
1693 for (j=1;j<=d_size; j++) {
1694 dd_set(sol[j-1],T[j-1][rhscol-1]);
1695 dd_TableauEntry(&x,m_size,d_size,A,T,objrow,j);
1696 dd_neg(dsol[j-1],x);
1697 dd_TableauEntry(optvalue,m_size,d_size,A,T,objrow,rhscol);
1698 if (localdebug) {fprintf(stderr,"dsol[%ld]= ",nbindex[j]); dd_WriteNumber(stderr, dsol[j-1]); }
1699 }
1700 for (i=1; i<=m_size; i++) {
1701 if (bflag[i]==-1) { /* i is a basic variable */
1702 dd_TableauEntry(&x,m_size,d_size,A,T,i,rhscol);
1703 if (dd_Positive(x)) set_addelem(posset, i);
1704 }
1705 }
1706
1707 break;
1708 case dd_Inconsistent:
1709 if (localdebug) fprintf(stderr,"SetSolutions: LP is inconsistent.\n");
1710 for (j=1;j<=d_size; j++) {
1711 dd_set(sol[j-1],T[j-1][rhscol-1]);
1712 dd_TableauEntry(&x,m_size,d_size,A,T,re,j);
1713 dd_neg(dsol[j-1],x);
1714 if (localdebug) {fprintf(stderr,"dsol[%ld]=",nbindex[j]);
1715 dd_WriteNumber(stderr,dsol[j-1]);
1716 fprintf(stderr,"\n");
1717 }
1718 }
1719 break;
1720
1721 case dd_DualInconsistent:
1722 if (localdebug) printf( "SetSolutions: LP is dual inconsistent.\n");
1723 for (j=1;j<=d_size; j++) {
1724 dd_set(sol[j-1],T[j-1][se-1]);
1725 dd_TableauEntry(&x,m_size,d_size,A,T,objrow,j);
1726 dd_neg(dsol[j-1],x);
1727 if (localdebug) {fprintf(stderr,"dsol[%ld]=",nbindex[j]);
1728 dd_WriteNumber(stderr,dsol[j-1]);
1729 fprintf(stderr,"\n");
1730 }
1731 }
1732 break;
1733
1734 case dd_StrucDualInconsistent:
1735 dd_TableauEntry(&x,m_size,d_size,A,T,objrow,se);
1736 if (dd_Positive(x)) dd_set(sw,dd_one);
1737 else dd_neg(sw,dd_one);
1738 for (j=1;j<=d_size; j++) {
1739 dd_mul(sol[j-1],sw,T[j-1][se-1]);
1740 dd_TableauEntry(&x,m_size,d_size,A,T,objrow,j);
1741 dd_neg(dsol[j-1],x);
1742 if (localdebug) {fprintf(stderr,"dsol[%ld]= ",nbindex[j]);dd_WriteNumber(stderr,dsol[j-1]);}
1743 }
1744 if (localdebug) fprintf(stderr,"SetSolutions: LP is dual inconsistent.\n");
1745 break;
1746
1747 default:break;
1748 }
1749 dd_clear(x); dd_clear(sw);
1750 }
1751
1752
dd_RandomPermutation2(dd_rowindex OV,long t,unsigned int seed)1753 void dd_RandomPermutation2(dd_rowindex OV,long t,unsigned int seed)
1754 {
1755 long k,j,ovj;
1756 double u,xk,r,rand_max=(double) UINT64_MAX;
1757 int localdebug=dd_FALSE;
1758
1759 srand_splitmix64(seed);
1760 for (j=t; j>1 ; j--) {
1761 r=rand_splitmix64();
1762 u=r/rand_max;
1763 xk=(double)(j*u +1);
1764 k=(long)xk;
1765 if (localdebug) fprintf(stderr,"u=%g, k=%ld, r=%g, randmax= %g\n",u,k,r,rand_max);
1766 ovj=OV[j];
1767 OV[j]=OV[k];
1768 OV[k]=ovj;
1769 if (localdebug) fprintf(stderr,"row %ld is exchanged with %ld\n",j,k);
1770 }
1771 }
1772
dd_ComputeRowOrderVector2(dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_rowindex OV,dd_RowOrderType ho,unsigned int rseed)1773 void dd_ComputeRowOrderVector2(dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,
1774 dd_rowindex OV,dd_RowOrderType ho,unsigned int rseed)
1775 {
1776 long i,itemp;
1777
1778 OV[0]=0;
1779 switch (ho){
1780 case dd_MaxIndex:
1781 for(i=1; i<=m_size; i++) OV[i]=m_size-i+1;
1782 break;
1783
1784 case dd_LexMin:
1785 for(i=1; i<=m_size; i++) OV[i]=i;
1786 dd_QuickSort(OV,1,m_size,A,d_size);
1787 break;
1788
1789 case dd_LexMax:
1790 for(i=1; i<=m_size; i++) OV[i]=i;
1791 dd_QuickSort(OV,1,m_size,A,d_size);
1792 for(i=1; i<=m_size/2;i++){ /* just reverse the order */
1793 itemp=OV[i];
1794 OV[i]=OV[m_size-i+1];
1795 OV[m_size-i+1]=itemp;
1796 }
1797 break;
1798
1799 case dd_RandomRow:
1800 for(i=1; i<=m_size; i++) OV[i]=i;
1801 if (rseed<=0) rseed=1;
1802 dd_RandomPermutation2(OV,m_size,rseed);
1803 break;
1804
1805 case dd_MinIndex:
1806 for(i=1; i<=m_size; i++) OV[i]=i;
1807 break;
1808
1809 default:
1810 for(i=1; i<=m_size; i++) OV[i]=i;
1811 break;
1812 }
1813 }
1814
dd_SelectPreorderedNext2(dd_rowrange m_size,dd_colrange d_size,rowset excluded,dd_rowindex OV,dd_rowrange * hnext)1815 void dd_SelectPreorderedNext2(dd_rowrange m_size,dd_colrange d_size,
1816 rowset excluded,dd_rowindex OV,dd_rowrange *hnext)
1817 {
1818 dd_rowrange i,k;
1819
1820 *hnext=0;
1821 for (i=1; i<=m_size && *hnext==0; i++){
1822 k=OV[i];
1823 if (!set_member(k,excluded)) *hnext=k ;
1824 }
1825 }
1826
1827 #ifdef GMPRATIONAL
1828
Obj2Obj(dd_LPObjectiveType obj)1829 ddf_LPObjectiveType Obj2Obj(dd_LPObjectiveType obj)
1830 {
1831 ddf_LPObjectiveType objf=ddf_LPnone;
1832
1833 switch (obj) {
1834 case dd_LPnone: objf=ddf_LPnone; break;
1835 case dd_LPmax: objf=ddf_LPmax; break;
1836 case dd_LPmin: objf=ddf_LPmin; break;
1837 }
1838 return objf;
1839 }
1840
dd_LPgmp2LPf(dd_LPPtr lp)1841 ddf_LPPtr dd_LPgmp2LPf(dd_LPPtr lp)
1842 {
1843 dd_rowrange i;
1844 dd_colrange j;
1845 ddf_LPType *lpf;
1846 double val;
1847 dd_boolean localdebug=dd_FALSE;
1848
1849 if (localdebug) fprintf(stderr,"Converting a GMP-LP to a float-LP.\n");
1850
1851 lpf=ddf_CreateLPData(Obj2Obj(lp->objective), ddf_Real, lp->m, lp->d);
1852 lpf->Homogeneous = lp->Homogeneous;
1853 lpf->eqnumber=lp->eqnumber; /* this records the number of equations */
1854
1855 for (i = 1; i <= lp->m; i++) {
1856 if (set_member(i, lp->equalityset)) set_addelem(lpf->equalityset,i);
1857 /* it is equality. Its reversed row will not be in this set */
1858 for (j = 1; j <= lp->d; j++) {
1859 val=mpq_get_d(lp->A[i-1][j-1]);
1860 ddf_set_d(lpf->A[i-1][j-1],val);
1861 } /*of j*/
1862 } /*of i*/
1863
1864 return lpf;
1865 }
1866
1867
1868 #endif
1869
1870
dd_LPSolve(dd_LPPtr lp,dd_LPSolverType solver,dd_ErrorType * err)1871 dd_boolean dd_LPSolve(dd_LPPtr lp,dd_LPSolverType solver,dd_ErrorType *err)
1872 /*
1873 The current version of dd_LPSolve that solves an LP with floating-arithmetics first
1874 and then with the specified arithimetics if it is GMP.
1875
1876 When LP is inconsistent then *re returns the evidence row.
1877 When LP is dual-inconsistent then *se returns the evidence column.
1878 */
1879 {
1880 int i;
1881 dd_boolean found=dd_FALSE;
1882 #ifdef GMPRATIONAL
1883 ddf_LPPtr lpf;
1884 ddf_ErrorType errf;
1885 dd_boolean LPScorrect=dd_FALSE;
1886 dd_boolean localdebug=dd_FALSE;
1887 if (dd_debug) localdebug=dd_debug;
1888 #endif
1889
1890 *err=dd_NoError;
1891 lp->solver=solver;
1892
1893 time(&lp->starttime);
1894
1895 #ifndef GMPRATIONAL
1896 switch (lp->solver) {
1897 case dd_CrissCross:
1898 dd_CrissCrossSolve(lp,err);
1899 break;
1900 case dd_DualSimplex:
1901 dd_DualSimplexSolve(lp,err);
1902 break;
1903 }
1904 #else
1905 lpf=dd_LPgmp2LPf(lp);
1906 switch (lp->solver) {
1907 case dd_CrissCross:
1908 ddf_CrissCrossSolve(lpf,&errf); /* First, run with double float. */
1909 if (errf==ddf_NoError){ /* 094a: fix for a bug reported by Dima Pasechnik */
1910 dd_BasisStatus(lpf,lp, &LPScorrect); /* Check the basis. */
1911 } else {LPScorrect=dd_FALSE;}
1912 if (!LPScorrect) {
1913 if (localdebug) printf("BasisStatus: the current basis is NOT verified with GMP. Rerun with GMP.\n");
1914 dd_CrissCrossSolve(lp,err); /* Rerun with GMP if fails. */
1915 } else {
1916 if (localdebug) printf("BasisStatus: the current basis is verified with GMP. The LP Solved.\n");
1917 }
1918 break;
1919 case dd_DualSimplex:
1920 ddf_DualSimplexSolve(lpf,&errf); /* First, run with double float. */
1921 if (errf==ddf_NoError){ /* 094a: fix for a bug reported by Dima Pasechnik */
1922 dd_BasisStatus(lpf,lp, &LPScorrect); /* Check the basis. */
1923 } else {LPScorrect=dd_FALSE;}
1924 if (!LPScorrect){
1925 if (localdebug) printf("BasisStatus: the current basis is NOT verified with GMP. Rerun with GMP.\n");
1926 dd_DualSimplexSolve(lp,err); /* Rerun with GMP if fails. */
1927 if (localdebug){
1928 printf("*total number pivots = %ld (ph0 = %ld, ph1 = %ld, ph2 = %ld, ph3 = %ld, ph4 = %ld)\n",
1929 lp->total_pivots,lp->pivots[0],lp->pivots[1],lp->pivots[2],lp->pivots[3],lp->pivots[4]);
1930 ddf_WriteLPResult(stdout, lpf, errf);
1931 dd_WriteLP(stdout, lp);
1932 }
1933 } else {
1934 if (localdebug) printf("BasisStatus: the current basis is verified with GMP. The LP Solved.\n");
1935 }
1936 break;
1937 }
1938 ddf_FreeLPData(lpf);
1939 #endif
1940
1941 time(&lp->endtime);
1942 lp->total_pivots=0;
1943 for (i=0; i<=4; i++) lp->total_pivots+=lp->pivots[i];
1944 if (*err==dd_NoError) found=dd_TRUE;
1945 return found;
1946 }
1947
1948
dd_LPSolve0(dd_LPPtr lp,dd_LPSolverType solver,dd_ErrorType * err)1949 dd_boolean dd_LPSolve0(dd_LPPtr lp,dd_LPSolverType solver,dd_ErrorType *err)
1950 /*
1951 The original version of dd_LPSolve that solves an LP with specified arithimetics.
1952
1953 When LP is inconsistent then *re returns the evidence row.
1954 When LP is dual-inconsistent then *se returns the evidence column.
1955 */
1956 {
1957 int i;
1958 dd_boolean found=dd_FALSE;
1959
1960 *err=dd_NoError;
1961 lp->solver=solver;
1962 time(&lp->starttime);
1963
1964 switch (lp->solver) {
1965 case dd_CrissCross:
1966 dd_CrissCrossSolve(lp,err);
1967 break;
1968 case dd_DualSimplex:
1969 dd_DualSimplexSolve(lp,err);
1970 break;
1971 }
1972
1973 time(&lp->endtime);
1974 lp->total_pivots=0;
1975 for (i=0; i<=4; i++) lp->total_pivots+=lp->pivots[i];
1976 if (*err==dd_NoError) found=dd_TRUE;
1977 return found;
1978 }
1979
1980
dd_MakeLPforInteriorFinding(dd_LPPtr lp)1981 dd_LPPtr dd_MakeLPforInteriorFinding(dd_LPPtr lp)
1982 /* Delete the objective row,
1983 add an extra column with -1's to the matrix A,
1984 add an extra row with (bceil, 0,...,0,-1),
1985 add an objective row with (0,...,0,1), and
1986 rows & columns, and change m_size and d_size accordingly, to output new_A.
1987 This sets up the LP:
1988 maximize x_{d+1}
1989 s.t. A x + x_{d+1} <= b
1990 x_{d+1} <= bm * bmax,
1991 where bm is set to 2 by default, and bmax=max{1, b[1],...,b[m_size]}.
1992 Note that the equalitions (linearity) in the input lp will be ignored.
1993 */
1994 {
1995 dd_rowrange m;
1996 dd_colrange d;
1997 dd_NumberType numbtype;
1998 dd_LPObjectiveType obj;
1999 dd_LPType *lpnew;
2000 dd_rowrange i;
2001 dd_colrange j;
2002 mytype bm,bmax,bceil;
2003 int localdebug=dd_FALSE;
2004
2005 dd_init(bm); dd_init(bmax); dd_init(bceil);
2006 dd_add(bm,dd_one,dd_one); dd_set(bmax,dd_one);
2007 numbtype=lp->numbtype;
2008 m=lp->m+1;
2009 d=lp->d+1;
2010 obj=dd_LPmax;
2011
2012 lpnew=dd_CreateLPData(obj, numbtype, m, d);
2013
2014 for (i=1; i<=lp->m; i++) {
2015 if (dd_Larger(lp->A[i-1][lp->rhscol-1],bmax))
2016 dd_set(bmax,lp->A[i-1][lp->rhscol-1]);
2017 }
2018 dd_mul(bceil,bm,bmax);
2019 if (localdebug) {fprintf(stderr,"bceil is set to "); dd_WriteNumber(stderr, bceil);}
2020
2021 for (i=1; i <= lp->m; i++) {
2022 for (j=1; j <= lp->d; j++) {
2023 dd_set(lpnew->A[i-1][j-1],lp->A[i-1][j-1]);
2024 }
2025 }
2026
2027 for (i=1;i<=lp->m; i++){
2028 dd_neg(lpnew->A[i-1][lp->d],dd_one); /* new column with all minus one's */
2029 }
2030
2031 for (j=1;j<=lp->d;j++){
2032 dd_set(lpnew->A[m-2][j-1],dd_purezero); /* new row (bceil, 0,...,0,-1) */
2033 }
2034 dd_set(lpnew->A[m-2][0],bceil); /* new row (bceil, 0,...,0,-1) */
2035
2036 for (j=1;j<= d-1;j++) {
2037 dd_set(lpnew->A[m-1][j-1],dd_purezero); /* new obj row with (0,...,0,1) */
2038 }
2039 dd_set(lpnew->A[m-1][d-1],dd_one); /* new obj row with (0,...,0,1) */
2040
2041 if (localdebug) dd_WriteAmatrix(stderr, lp->A, lp->m, lp->d);
2042 if (localdebug) dd_WriteAmatrix(stderr, lpnew->A, lpnew->m, lpnew->d);
2043 dd_clear(bm); dd_clear(bmax); dd_clear(bceil);
2044
2045 return lpnew;
2046 }
2047
2048
dd_WriteLPResult(FILE * f,dd_LPPtr lp,dd_ErrorType err)2049 void dd_WriteLPResult(FILE *f,dd_LPPtr lp,dd_ErrorType err)
2050 {
2051 long j;
2052
2053 fprintf(f,"* cdd LP solver result\n");
2054
2055 if (err!=dd_NoError) {
2056 dd_WriteErrorMessages(f,err);
2057 goto _L99;
2058 }
2059
2060 dd_WriteProgramDescription(f);
2061
2062 fprintf(f,"* #constraints = %ld\n",lp->m-1);
2063 fprintf(f,"* #variables = %ld\n",lp->d-1);
2064
2065 switch (lp->solver) {
2066 case dd_DualSimplex:
2067 fprintf(f,"* Algorithm: dual simplex algorithm\n");break;
2068 case dd_CrissCross:
2069 fprintf(f,"* Algorithm: criss-cross method\n");break;
2070 }
2071
2072 switch (lp->objective) {
2073 case dd_LPmax:
2074 fprintf(f,"* maximization is chosen\n");break;
2075 case dd_LPmin:
2076 fprintf(f,"* minimization is chosen\n");break;
2077 case dd_LPnone:
2078 fprintf(f,"* no objective type (max or min) is chosen\n");break;
2079 }
2080
2081 if (lp->objective==dd_LPmax||lp->objective==dd_LPmin){
2082 fprintf(f,"* Objective function is\n");
2083 for (j=0; j<lp->d; j++){
2084 if (j>0 && dd_Nonnegative(lp->A[lp->objrow-1][j]) ) fprintf(f," +");
2085 if (j>0 && (j % 5) == 0) fprintf(f,"\n");
2086 dd_WriteNumber(f,lp->A[lp->objrow-1][j]);
2087 if (j>0) fprintf(f," X[%3ld]",j);
2088 }
2089 fprintf(f,"\n");
2090 }
2091
2092 switch (lp->LPS){
2093 case dd_Optimal:
2094 fprintf(f,"* LP status: a dual pair (x,y) of optimal solutions found.\n");
2095 fprintf(f,"begin\n");
2096 fprintf(f," primal_solution\n");
2097 for (j=1; j<lp->d; j++) {
2098 fprintf(f," %3ld : ",j);
2099 dd_WriteNumber(f,lp->sol[j]);
2100 fprintf(f,"\n");
2101 }
2102 fprintf(f," dual_solution\n");
2103 for (j=1; j<lp->d; j++){
2104 if (lp->nbindex[j+1]>0) {
2105 fprintf(f," %3ld : ",lp->nbindex[j+1]);
2106 dd_WriteNumber(f,lp->dsol[j]); fprintf(f,"\n");
2107 }
2108 }
2109 fprintf(f," optimal_value : "); dd_WriteNumber(f,lp->optvalue);
2110 fprintf(f,"\nend\n");
2111 break;
2112
2113 case dd_Inconsistent:
2114 fprintf(f,"* LP status: LP is inconsistent.\n");
2115 fprintf(f,"* The positive combination of original inequalities with\n");
2116 fprintf(f,"* the following coefficients will prove the inconsistency.\n");
2117 fprintf(f,"begin\n");
2118 fprintf(f," dual_direction\n");
2119 fprintf(f," %3ld : ",lp->re);
2120 dd_WriteNumber(f,dd_one); fprintf(f,"\n");
2121 for (j=1; j<lp->d; j++){
2122 if (lp->nbindex[j+1]>0) {
2123 fprintf(f," %3ld : ",lp->nbindex[j+1]);
2124 dd_WriteNumber(f,lp->dsol[j]); fprintf(f,"\n");
2125 }
2126 }
2127 fprintf(f,"end\n");
2128 break;
2129
2130 case dd_DualInconsistent: case dd_StrucDualInconsistent:
2131 fprintf(f,"* LP status: LP is dual inconsistent.\n");
2132 fprintf(f,"* The linear combination of columns with\n");
2133 fprintf(f,"* the following coefficients will prove the dual inconsistency.\n");
2134 fprintf(f,"* (It is also an unbounded direction for the primal LP.)\n");
2135 fprintf(f,"begin\n");
2136 fprintf(f," primal_direction\n");
2137 for (j=1; j<lp->d; j++) {
2138 fprintf(f," %3ld : ",j);
2139 dd_WriteNumber(f,lp->sol[j]);
2140 fprintf(f,"\n");
2141 }
2142 fprintf(f,"end\n");
2143 break;
2144
2145 default:
2146 break;
2147 }
2148 fprintf(f,"* number of pivot operations = %ld (ph0 = %ld, ph1 = %ld, ph2 = %ld, ph3 = %ld, ph4 = %ld)\n",lp->total_pivots,lp->pivots[0],lp->pivots[1],lp->pivots[2],lp->pivots[3],lp->pivots[4]);
2149 dd_WriteLPTimes(f, lp);
2150 _L99:;
2151 }
2152
dd_CreateLP_H_ImplicitLinearity(dd_MatrixPtr M)2153 dd_LPPtr dd_CreateLP_H_ImplicitLinearity(dd_MatrixPtr M)
2154 {
2155 dd_rowrange m, i, irev, linc;
2156 dd_colrange d, j;
2157 dd_LPPtr lp;
2158 dd_boolean localdebug=dd_FALSE;
2159
2160 linc=set_card(M->linset);
2161 m=M->rowsize+1+linc+1;
2162 /* We represent each equation by two inequalities.
2163 This is not the best way but makes the code simple. */
2164 d=M->colsize+1;
2165
2166 lp=dd_CreateLPData(M->objective, M->numbtype, m, d);
2167 lp->Homogeneous = dd_TRUE;
2168 lp->objective = dd_LPmax;
2169 lp->eqnumber=linc; /* this records the number of equations */
2170 lp->redcheck_extensive=dd_FALSE; /* this is default */
2171
2172 irev=M->rowsize; /* the first row of the linc reversed inequalities. */
2173 for (i = 1; i <= M->rowsize; i++) {
2174 if (set_member(i, M->linset)) {
2175 irev=irev+1;
2176 set_addelem(lp->equalityset,i); /* it is equality. */
2177 /* the reversed row irev is not in the equality set. */
2178 for (j = 1; j <= M->colsize; j++) {
2179 dd_neg(lp->A[irev-1][j-1],M->matrix[i-1][j-1]);
2180 } /*of j*/
2181 } else {
2182 dd_set(lp->A[i-1][d-1],dd_minusone); /* b_I + A_I x - 1 z >= 0 (z=x_d) */
2183 }
2184 for (j = 1; j <= M->colsize; j++) {
2185 dd_set(lp->A[i-1][j-1],M->matrix[i-1][j-1]);
2186 if (j==1 && i<M->rowsize && dd_Nonzero(M->matrix[i-1][j-1])) lp->Homogeneous = dd_FALSE;
2187 } /*of j*/
2188 } /*of i*/
2189 dd_set(lp->A[m-2][0],dd_one); dd_set(lp->A[m-2][d-1],dd_minusone);
2190 /* make the LP bounded. */
2191
2192 dd_set(lp->A[m-1][d-1],dd_one);
2193 /* objective is to maximize z. */
2194
2195 if (localdebug) {
2196 fprintf(stderr,"dd_CreateLP_H_ImplicitLinearity: an new lp is\n");
2197 dd_WriteLP(stderr,lp);
2198 }
2199
2200 return lp;
2201 }
2202
dd_CreateLP_V_ImplicitLinearity(dd_MatrixPtr M)2203 dd_LPPtr dd_CreateLP_V_ImplicitLinearity(dd_MatrixPtr M)
2204 {
2205 dd_rowrange m, i, irev, linc;
2206 dd_colrange d, j;
2207 dd_LPPtr lp;
2208 dd_boolean localdebug=dd_FALSE;
2209
2210 linc=set_card(M->linset);
2211 m=M->rowsize+1+linc+1;
2212 /* We represent each equation by two inequalities.
2213 This is not the best way but makes the code simple. */
2214 d=(M->colsize)+2;
2215 /* Two more columns. This is different from the H-reprentation case */
2216
2217 /* The below must be modified for V-representation!!! */
2218
2219 lp=dd_CreateLPData(M->objective, M->numbtype, m, d);
2220 lp->Homogeneous = dd_FALSE;
2221 lp->objective = dd_LPmax;
2222 lp->eqnumber=linc; /* this records the number of equations */
2223 lp->redcheck_extensive=dd_FALSE; /* this is default */
2224
2225 irev=M->rowsize; /* the first row of the linc reversed inequalities. */
2226 for (i = 1; i <= M->rowsize; i++) {
2227 dd_set(lp->A[i-1][0],dd_purezero); /* It is almost completely degerate LP */
2228 if (set_member(i, M->linset)) {
2229 irev=irev+1;
2230 set_addelem(lp->equalityset,i); /* it is equality. */
2231 /* the reversed row irev is not in the equality set. */
2232 for (j = 2; j <= (M->colsize)+1; j++) {
2233 dd_neg(lp->A[irev-1][j-1],M->matrix[i-1][j-2]);
2234 } /*of j*/
2235 if (localdebug) fprintf(stderr,"equality row %ld generates the reverse row %ld.\n",i,irev);
2236 } else {
2237 dd_set(lp->A[i-1][d-1],dd_minusone); /* b_I x_0 + A_I x - 1 z >= 0 (z=x_d) */
2238 }
2239 for (j = 2; j <= (M->colsize)+1; j++) {
2240 dd_set(lp->A[i-1][j-1],M->matrix[i-1][j-2]);
2241 } /*of j*/
2242 } /*of i*/
2243 dd_set(lp->A[m-2][0],dd_one); dd_set(lp->A[m-2][d-1],dd_minusone);
2244 /* make the LP bounded. */
2245 dd_set(lp->A[m-1][d-1],dd_one);
2246 /* objective is to maximize z. */
2247
2248 if (localdebug) {
2249 fprintf(stderr,"dd_CreateLP_V_ImplicitLinearity: an new lp is\n");
2250 dd_WriteLP(stderr,lp);
2251 }
2252
2253 return lp;
2254 }
2255
2256
dd_CreateLP_H_Redundancy(dd_MatrixPtr M,dd_rowrange itest)2257 dd_LPPtr dd_CreateLP_H_Redundancy(dd_MatrixPtr M, dd_rowrange itest)
2258 {
2259 dd_rowrange m, i, irev, linc;
2260 dd_colrange d, j;
2261 dd_LPPtr lp;
2262 dd_boolean localdebug=dd_FALSE;
2263
2264 linc=set_card(M->linset);
2265 m=M->rowsize+1+linc;
2266 /* We represent each equation by two inequalities.
2267 This is not the best way but makes the code simple. */
2268 d=M->colsize;
2269
2270 lp=dd_CreateLPData(M->objective, M->numbtype, m, d);
2271 lp->Homogeneous = dd_TRUE;
2272 lp->objective = dd_LPmin;
2273 lp->eqnumber=linc; /* this records the number of equations */
2274 lp->redcheck_extensive=dd_FALSE; /* this is default */
2275
2276 irev=M->rowsize; /* the first row of the linc reversed inequalities. */
2277 for (i = 1; i <= M->rowsize; i++) {
2278 if (set_member(i, M->linset)) {
2279 irev=irev+1;
2280 set_addelem(lp->equalityset,i); /* it is equality. */
2281 /* the reversed row irev is not in the equality set. */
2282 for (j = 1; j <= M->colsize; j++) {
2283 dd_neg(lp->A[irev-1][j-1],M->matrix[i-1][j-1]);
2284 } /*of j*/
2285 if (localdebug) fprintf(stderr,"equality row %ld generates the reverse row %ld.\n",i,irev);
2286 }
2287 for (j = 1; j <= M->colsize; j++) {
2288 dd_set(lp->A[i-1][j-1],M->matrix[i-1][j-1]);
2289 if (j==1 && i<M->rowsize && dd_Nonzero(M->matrix[i-1][j-1])) lp->Homogeneous = dd_FALSE;
2290 } /*of j*/
2291 } /*of i*/
2292 for (j = 1; j <= M->colsize; j++) {
2293 dd_set(lp->A[m-1][j-1],M->matrix[itest-1][j-1]);
2294 /* objective is to violate the inequality in question. */
2295 } /*of j*/
2296 dd_add(lp->A[itest-1][0],lp->A[itest-1][0],dd_one); /* relax the original inequality by one */
2297
2298 return lp;
2299 }
2300
2301
dd_CreateLP_V_Redundancy(dd_MatrixPtr M,dd_rowrange itest)2302 dd_LPPtr dd_CreateLP_V_Redundancy(dd_MatrixPtr M, dd_rowrange itest)
2303 {
2304 dd_rowrange m, i, irev, linc;
2305 dd_colrange d, j;
2306 dd_LPPtr lp;
2307 dd_boolean localdebug=dd_FALSE;
2308
2309 linc=set_card(M->linset);
2310 m=M->rowsize+1+linc;
2311 /* We represent each equation by two inequalities.
2312 This is not the best way but makes the code simple. */
2313 d=(M->colsize)+1;
2314 /* One more column. This is different from the H-reprentation case */
2315
2316 /* The below must be modified for V-representation!!! */
2317
2318 lp=dd_CreateLPData(M->objective, M->numbtype, m, d);
2319 lp->Homogeneous = dd_FALSE;
2320 lp->objective = dd_LPmin;
2321 lp->eqnumber=linc; /* this records the number of equations */
2322 lp->redcheck_extensive=dd_FALSE; /* this is default */
2323
2324 irev=M->rowsize; /* the first row of the linc reversed inequalities. */
2325 for (i = 1; i <= M->rowsize; i++) {
2326 if (i==itest){
2327 dd_set(lp->A[i-1][0],dd_one); /* this is to make the LP bounded, ie. the min >= -1 */
2328 } else {
2329 dd_set(lp->A[i-1][0],dd_purezero); /* It is almost completely degerate LP */
2330 }
2331 if (set_member(i, M->linset)) {
2332 irev=irev+1;
2333 set_addelem(lp->equalityset,i); /* it is equality. */
2334 /* the reversed row irev is not in the equality set. */
2335 for (j = 2; j <= (M->colsize)+1; j++) {
2336 dd_neg(lp->A[irev-1][j-1],M->matrix[i-1][j-2]);
2337 } /*of j*/
2338 if (localdebug) fprintf(stderr,"equality row %ld generates the reverse row %ld.\n",i,irev);
2339 }
2340 for (j = 2; j <= (M->colsize)+1; j++) {
2341 dd_set(lp->A[i-1][j-1],M->matrix[i-1][j-2]);
2342 } /*of j*/
2343 } /*of i*/
2344 for (j = 2; j <= (M->colsize)+1; j++) {
2345 dd_set(lp->A[m-1][j-1],M->matrix[itest-1][j-2]);
2346 /* objective is to violate the inequality in question. */
2347 } /*of j*/
2348 dd_set(lp->A[m-1][0],dd_purezero); /* the constant term for the objective is zero */
2349
2350 if (localdebug) dd_WriteLP(stdout, lp);
2351
2352 return lp;
2353 }
2354
2355
dd_CreateLP_V_SRedundancy(dd_MatrixPtr M,dd_rowrange itest)2356 dd_LPPtr dd_CreateLP_V_SRedundancy(dd_MatrixPtr M, dd_rowrange itest)
2357 {
2358 /*
2359 V-representation (=boundary problem)
2360 g* = maximize
2361 1^T b_{I-itest} x_0 + 1^T A_{I-itest} (the sum of slacks)
2362 subject to
2363 b_itest x_0 + A_itest x = 0 (the point has to lie on the boundary)
2364 b_{I-itest} x_0 + A_{I-itest} x >= 0 (all nonlinearity generators in one side)
2365 1^T b_{I-itest} x_0 + 1^T A_{I-itest} x <= 1 (to make an LP bounded)
2366 b_L x_0 + A_L x = 0. (linearity generators)
2367
2368 The redundant row is strongly redundant if and only if g* is zero.
2369 */
2370
2371 dd_rowrange m, i, irev, linc;
2372 dd_colrange d, j;
2373 dd_LPPtr lp;
2374 dd_boolean localdebug=dd_FALSE;
2375
2376 linc=set_card(M->linset);
2377 m=M->rowsize+1+linc+2;
2378 /* We represent each equation by two inequalities.
2379 This is not the best way but makes the code simple.
2380 Two extra constraints are for the first equation and the bouding inequality.
2381 */
2382 d=(M->colsize)+1;
2383 /* One more column. This is different from the H-reprentation case */
2384
2385 /* The below must be modified for V-representation!!! */
2386
2387 lp=dd_CreateLPData(M->objective, M->numbtype, m, d);
2388 lp->Homogeneous = dd_FALSE;
2389 lp->objective = dd_LPmax;
2390 lp->eqnumber=linc; /* this records the number of equations */
2391
2392 irev=M->rowsize; /* the first row of the linc reversed inequalities. */
2393 for (i = 1; i <= M->rowsize; i++) {
2394 if (i==itest){
2395 dd_set(lp->A[i-1][0],dd_purezero); /* this is a half of the boundary constraint. */
2396 } else {
2397 dd_set(lp->A[i-1][0],dd_purezero); /* It is almost completely degerate LP */
2398 }
2399 if (set_member(i, M->linset) || i==itest) {
2400 irev=irev+1;
2401 set_addelem(lp->equalityset,i); /* it is equality. */
2402 /* the reversed row irev is not in the equality set. */
2403 for (j = 2; j <= (M->colsize)+1; j++) {
2404 dd_neg(lp->A[irev-1][j-1],M->matrix[i-1][j-2]);
2405 } /*of j*/
2406 if (localdebug) fprintf(stderr,"equality row %ld generates the reverse row %ld.\n",i,irev);
2407 }
2408 for (j = 2; j <= (M->colsize)+1; j++) {
2409 dd_set(lp->A[i-1][j-1],M->matrix[i-1][j-2]);
2410 dd_add(lp->A[m-1][j-1],lp->A[m-1][j-1],lp->A[i-1][j-1]); /* the objective is the sum of all ineqalities */
2411 } /*of j*/
2412 } /*of i*/
2413 for (j = 2; j <= (M->colsize)+1; j++) {
2414 dd_neg(lp->A[m-2][j-1],lp->A[m-1][j-1]);
2415 /* to make an LP bounded. */
2416 } /*of j*/
2417 dd_set(lp->A[m-2][0],dd_one); /* the constant term for the bounding constraint is 1 */
2418
2419 if (localdebug) dd_WriteLP(stdout, lp);
2420
2421 return lp;
2422 }
2423
dd_Redundant(dd_MatrixPtr M,dd_rowrange itest,dd_Arow certificate,dd_ErrorType * error)2424 dd_boolean dd_Redundant(dd_MatrixPtr M, dd_rowrange itest, dd_Arow certificate, dd_ErrorType *error)
2425 /* 092 */
2426 {
2427 /* Checks whether the row itest is redundant for the representation.
2428 All linearity rows are not checked and considered NONredundant.
2429 This code works for both H- and V-representations. A certificate is
2430 given in the case of non-redundancy, showing a solution x violating only the itest
2431 inequality for H-representation, a hyperplane RHS and normal (x_0, x) that
2432 separates the itest from the rest. More explicitly, the LP to be setup is
2433
2434 H-representation
2435 f* = minimize
2436 b_itest + A_itest x
2437 subject to
2438 b_itest + 1 + A_itest x >= 0 (relaxed inequality to make an LP bounded)
2439 b_{I-itest} + A_{I-itest} x >= 0 (all inequalities except for itest)
2440 b_L + A_L x = 0. (linearity)
2441
2442 V-representation (=separation problem)
2443 f* = minimize
2444 b_itest x_0 + A_itest x
2445 subject to
2446 b_itest x_0 + A_itest x >= -1 (to make an LP bounded)
2447 b_{I-itest} x_0 + A_{I-itest} x >= 0 (all nonlinearity generators except for itest in one side)
2448 b_L x_0 + A_L x = 0. (linearity generators)
2449
2450 Here, the input matrix is considered as (b, A), i.e. b corresponds to the first column of input
2451 and the row indices of input is partitioned into I and L where L is the set of linearity.
2452 In both cases, the itest data is nonredundant if and only if the optimal value f* is negative.
2453 The certificate has dimension one more for V-representation case.
2454 */
2455
2456 dd_colrange j;
2457 dd_LPPtr lp;
2458 dd_LPSolutionPtr lps;
2459 dd_ErrorType err=dd_NoError;
2460 dd_boolean answer=dd_FALSE,localdebug=dd_FALSE;
2461
2462 *error=dd_NoError;
2463 if (set_member(itest, M->linset)){
2464 if (localdebug) printf("The %ld th row is linearity and redundancy checking is skipped.\n",itest);
2465 goto _L99;
2466 }
2467
2468 /* Create an LP data for redundancy checking */
2469 if (M->representation==dd_Generator){
2470 lp=dd_CreateLP_V_Redundancy(M, itest);
2471 } else {
2472 lp=dd_CreateLP_H_Redundancy(M, itest);
2473 }
2474
2475 dd_LPSolve(lp,dd_choiceRedcheckAlgorithm,&err);
2476 if (err!=dd_NoError){
2477 *error=err;
2478 goto _L999;
2479 } else {
2480 lps=dd_CopyLPSolution(lp);
2481
2482 for (j=0; j<lps->d; j++) {
2483 dd_set(certificate[j], lps->sol[j]);
2484 }
2485
2486 if (dd_Negative(lps->optvalue)){
2487 answer=dd_FALSE;
2488 if (localdebug) fprintf(stderr,"==> %ld th row is nonredundant.\n",itest);
2489 } else {
2490 answer=dd_TRUE;
2491 if (localdebug) fprintf(stderr,"==> %ld th row is redundant.\n",itest);
2492 }
2493 dd_FreeLPSolution(lps);
2494 }
2495 _L999:
2496 dd_FreeLPData(lp);
2497 _L99:
2498 return answer;
2499 }
2500
dd_RedundantExtensive(dd_MatrixPtr M,dd_rowrange itest,dd_Arow certificate,dd_rowset * redset,dd_ErrorType * error)2501 dd_boolean dd_RedundantExtensive(dd_MatrixPtr M, dd_rowrange itest, dd_Arow certificate,
2502 dd_rowset *redset,dd_ErrorType *error)
2503 /* 094 */
2504 {
2505 /* This uses the same LP construction as dd_Reduandant. But, while it is checking
2506 the redundancy of itest, it also tries to find some other variable that are
2507 redundant (i.e. forced to be nonnegative). This is expensive as it used
2508 the complete tableau information at each DualSimplex pivot. The redset must
2509 be initialized before this function is called.
2510 */
2511
2512 dd_colrange j;
2513 dd_LPPtr lp;
2514 dd_LPSolutionPtr lps;
2515 dd_ErrorType err=dd_NoError;
2516 dd_boolean answer=dd_FALSE,localdebug=dd_FALSE;
2517
2518 *error=dd_NoError;
2519 if (set_member(itest, M->linset)){
2520 if (localdebug) printf("The %ld th row is linearity and redundancy checking is skipped.\n",itest);
2521 goto _L99;
2522 }
2523
2524 /* Create an LP data for redundancy checking */
2525 if (M->representation==dd_Generator){
2526 lp=dd_CreateLP_V_Redundancy(M, itest);
2527 } else {
2528 lp=dd_CreateLP_H_Redundancy(M, itest);
2529 }
2530
2531 lp->redcheck_extensive=dd_TRUE;
2532
2533 dd_LPSolve0(lp,dd_DualSimplex,&err);
2534 if (err!=dd_NoError){
2535 *error=err;
2536 goto _L999;
2537 } else {
2538 set_copy(*redset,lp->redset_extra);
2539 set_delelem(*redset, itest);
2540 /* itest row might be redundant in the lp but this has nothing to do with its redundancy
2541 in the original system M. Thus we must delete it. */
2542 if (localdebug){
2543 fprintf(stderr, "dd_RedundantExtensive: checking for %ld, extra redset with cardinality %ld (%ld)\n",itest,set_card(*redset),set_card(lp->redset_extra));
2544 set_fwrite(stderr, *redset); fprintf(stderr, "\n");
2545 }
2546 lps=dd_CopyLPSolution(lp);
2547
2548 for (j=0; j<lps->d; j++) {
2549 dd_set(certificate[j], lps->sol[j]);
2550 }
2551
2552 if (dd_Negative(lps->optvalue)){
2553 answer=dd_FALSE;
2554 if (localdebug) fprintf(stderr,"==> %ld th row is nonredundant.\n",itest);
2555 } else {
2556 answer=dd_TRUE;
2557 if (localdebug) fprintf(stderr,"==> %ld th row is redundant.\n",itest);
2558 }
2559 dd_FreeLPSolution(lps);
2560 }
2561 _L999:
2562 dd_FreeLPData(lp);
2563 _L99:
2564 return answer;
2565 }
2566
dd_RedundantRows(dd_MatrixPtr M,dd_ErrorType * error)2567 dd_rowset dd_RedundantRows(dd_MatrixPtr M, dd_ErrorType *error) /* 092 */
2568 {
2569 dd_rowrange i,m;
2570 dd_colrange d;
2571 dd_rowset redset;
2572 dd_MatrixPtr Mcopy;
2573 dd_Arow cvec; /* certificate */
2574 dd_boolean localdebug=dd_FALSE;
2575
2576 m=M->rowsize;
2577 if (M->representation==dd_Generator){
2578 d=(M->colsize)+1;
2579 } else {
2580 d=M->colsize;
2581 }
2582 Mcopy=dd_MatrixCopy(M);
2583 dd_InitializeArow(d,&cvec);
2584 set_initialize(&redset, m);
2585 for (i=m; i>=1; i--) {
2586 if (dd_Redundant(Mcopy, i, cvec, error)) {
2587 if (localdebug) printf("Iteration %ld: the row %ld is redundant.\n",m-i+1,i);
2588 set_addelem(redset, i);
2589 dd_MatrixRowRemove(&Mcopy, i);
2590 } else {
2591 if (localdebug) printf("Iteration %ld: the row %ld is essential.\n",m-i+1, i);
2592 }
2593 if (*error!=dd_NoError) goto _L99;
2594 }
2595 _L99:
2596 dd_FreeMatrix(Mcopy);
2597 dd_FreeArow(d, cvec);
2598 return redset;
2599 }
2600
2601
dd_MatrixRedundancyRemove(dd_MatrixPtr * M,dd_rowset * redset,dd_rowindex * newpos,dd_ErrorType * error)2602 dd_boolean dd_MatrixRedundancyRemove(dd_MatrixPtr *M, dd_rowset *redset,dd_rowindex *newpos, dd_ErrorType *error) /* 094 */
2603 {
2604 /* It returns the set of all redundant rows. This should be called after all
2605 implicit linearity are recognized with dd_MatrixCanonicalizeLinearity.
2606 */
2607
2608
2609 dd_rowrange i,k,m,m1;
2610 dd_colrange d;
2611 dd_rowset redset1;
2612 dd_rowindex newpos1;
2613 dd_MatrixPtr M1=NULL;
2614 dd_Arow cvec; /* certificate */
2615 dd_boolean success=dd_FALSE, localdebug=dd_FALSE;
2616
2617 m=(*M)->rowsize;
2618 set_initialize(redset, m);
2619 M1=dd_MatrixSortedUniqueCopy(*M,newpos);
2620 for (i=1; i<=m; i++){
2621 if ((*newpos)[i]<=0) set_addelem(*redset,i);
2622 if (localdebug) printf(" %ld:%ld",i,(*newpos)[i]);
2623 }
2624 if (localdebug) printf("\n");
2625
2626 if ((*M)->representation==dd_Generator){
2627 d=((*M)->colsize)+1;
2628 } else {
2629 d=(*M)->colsize;
2630 }
2631 m1=M1->rowsize;
2632 if (localdebug){
2633 fprintf(stderr,"dd_MatrixRedundancyRemove: By sorting, %ld rows have been removed. The remaining has %ld rows.\n",m-m1,m1);
2634 /* dd_WriteMatrix(stdout,M1); */
2635 }
2636 dd_InitializeArow(d,&cvec);
2637 set_initialize(&redset1, M1->rowsize);
2638 k=1;
2639 do {
2640 if (dd_RedundantExtensive(M1, k, cvec, &redset1,error)) {
2641 set_addelem(redset1, k);
2642 dd_MatrixRowsRemove2(&M1,redset1,&newpos1);
2643 for (i=1; i<=m; i++){
2644 if ((*newpos)[i]>0){
2645 if (set_member((*newpos)[i],redset1)){
2646 set_addelem(*redset,i);
2647 (*newpos)[i]=0; /* now the original row i is recognized redundant and removed from M1 */
2648 } else {
2649 (*newpos)[i]=newpos1[(*newpos)[i]]; /* update the new pos vector */
2650 }
2651 }
2652 }
2653 set_free(redset1);
2654 set_initialize(&redset1, M1->rowsize);
2655 if (localdebug) {
2656 printf("dd_MatrixRedundancyRemove: the row %ld is redundant. The new matrix has %ld rows.\n", k, M1->rowsize);
2657 /* dd_WriteMatrix(stderr, M1); */
2658 }
2659 free(newpos1);
2660 } else {
2661 if (set_card(redset1)>0) {
2662 dd_MatrixRowsRemove2(&M1,redset1,&newpos1);
2663 for (i=1; i<=m; i++){
2664 if ((*newpos)[i]>0){
2665 if (set_member((*newpos)[i],redset1)){
2666 set_addelem(*redset,i);
2667 (*newpos)[i]=0; /* now the original row i is recognized redundant and removed from M1 */
2668 } else {
2669 (*newpos)[i]=newpos1[(*newpos)[i]]; /* update the new pos vector */
2670 }
2671 }
2672 }
2673 set_free(redset1);
2674 set_initialize(&redset1, M1->rowsize);
2675 free(newpos1);
2676 }
2677 if (localdebug) {
2678 printf("dd_MatrixRedundancyRemove: the row %ld is essential. The new matrix has %ld rows.\n", k, M1->rowsize);
2679 /* dd_WriteMatrix(stderr, M1); */
2680 }
2681 k=k+1;
2682 }
2683 if (*error!=dd_NoError) goto _L99;
2684 } while (k<=M1->rowsize);
2685 if (localdebug) dd_WriteMatrix(stderr, M1);
2686 success=dd_TRUE;
2687
2688 _L99:
2689 dd_FreeMatrix(*M);
2690 *M=M1;
2691 dd_FreeArow(d, cvec);
2692 set_free(redset1);
2693 return success;
2694 }
2695
2696
dd_SRedundant(dd_MatrixPtr M,dd_rowrange itest,dd_Arow certificate,dd_ErrorType * error)2697 dd_boolean dd_SRedundant(dd_MatrixPtr M, dd_rowrange itest, dd_Arow certificate, dd_ErrorType *error)
2698 /* 093a */
2699 {
2700 /* Checks whether the row itest is strongly redundant for the representation.
2701 A row is strongly redundant in H-representation if every point in
2702 the polyhedron satisfies it with strict inequality.
2703 A row is strongly redundant in V-representation if this point is in
2704 the interior of the polyhedron.
2705
2706 All linearity rows are not checked and considered NOT strongly redundant.
2707 This code works for both H- and V-representations. A certificate is
2708 given in the case of non-redundancy, showing a solution x violating only the itest
2709 inequality for H-representation, a hyperplane RHS and normal (x_0, x) that
2710 separates the itest from the rest. More explicitly, the LP to be setup is
2711
2712 H-representation
2713 f* = minimize
2714 b_itest + A_itest x
2715 subject to
2716 b_itest + 1 + A_itest x >= 0 (relaxed inequality to make an LP bounded)
2717 b_{I-itest} + A_{I-itest} x >= 0 (all inequalities except for itest)
2718 b_L + A_L x = 0. (linearity)
2719
2720 V-representation (=separation problem)
2721 f* = minimize
2722 b_itest x_0 + A_itest x
2723 subject to
2724 b_itest x_0 + A_itest x >= -1 (to make an LP bounded)
2725 b_{I-itest} x_0 + A_{I-itest} x >= 0 (all nonlinearity generators except for itest in one side)
2726 b_L x_0 + A_L x = 0. (linearity generators)
2727
2728 Here, the input matrix is considered as (b, A), i.e. b corresponds to the first column of input
2729 and the row indices of input is partitioned into I and L where L is the set of linearity.
2730 In H-representation, the itest data is strongly redundant if and only if the optimal value f* is positive.
2731 In V-representation, the itest data is redundant if and only if the optimal value f* is zero (as the LP
2732 is homogeneous and the optimal value is always non-positive). To recognize strong redundancy, one
2733 can set up a second LP
2734
2735 V-representation (=boundary problem)
2736 g* = maximize
2737 1^T b_{I-itest} x_0 + 1^T A_{I-itest} (the sum of slacks)
2738 subject to
2739 b_itest x_0 + A_itest x = 0 (the point has to lie on the boundary)
2740 b_{I-itest} x_0 + A_{I-itest} x >= 0 (all nonlinearity generators in one side)
2741 1^T b_{I-itest} x_0 + 1^T A_{I-itest} x <= 1 (to make an LP bounded)
2742 b_L x_0 + A_L x = 0. (linearity generators)
2743
2744 The redundant row is strongly redundant if and only if g* is zero.
2745
2746 The certificate has dimension one more for V-representation case.
2747 */
2748
2749 dd_colrange j;
2750 dd_LPPtr lp;
2751 dd_LPSolutionPtr lps;
2752 dd_ErrorType err=dd_NoError;
2753 dd_boolean answer=dd_FALSE,localdebug=dd_FALSE;
2754
2755 *error=dd_NoError;
2756 if (set_member(itest, M->linset)){
2757 if (localdebug) printf("The %ld th row is linearity and strong redundancy checking is skipped.\n",itest);
2758 goto _L99;
2759 }
2760
2761 /* Create an LP data for redundancy checking */
2762 if (M->representation==dd_Generator){
2763 lp=dd_CreateLP_V_Redundancy(M, itest);
2764 } else {
2765 lp=dd_CreateLP_H_Redundancy(M, itest);
2766 }
2767
2768 dd_LPSolve(lp,dd_choiceRedcheckAlgorithm,&err);
2769 if (err!=dd_NoError){
2770 *error=err;
2771 goto _L999;
2772 } else {
2773 lps=dd_CopyLPSolution(lp);
2774
2775 for (j=0; j<lps->d; j++) {
2776 dd_set(certificate[j], lps->sol[j]);
2777 }
2778
2779 if (localdebug){
2780 printf("Optimum value:");
2781 dd_WriteNumber(stdout, lps->optvalue);
2782 printf("\n");
2783 }
2784
2785 if (M->representation==dd_Inequality){
2786 if (dd_Positive(lps->optvalue)){
2787 answer=dd_TRUE;
2788 if (localdebug) fprintf(stderr,"==> %ld th inequality is strongly redundant.\n",itest);
2789 } else {
2790 answer=dd_FALSE;
2791 if (localdebug) fprintf(stderr,"==> %ld th inequality is not strongly redundant.\n",itest);
2792 }
2793 } else {
2794 if (dd_Negative(lps->optvalue)){
2795 answer=dd_FALSE;
2796 if (localdebug) fprintf(stderr,"==> %ld th point is not strongly redundant.\n",itest);
2797 } else {
2798 /* for V-representation, we have to solve another LP */
2799 dd_FreeLPData(lp);
2800 dd_FreeLPSolution(lps);
2801 lp=dd_CreateLP_V_SRedundancy(M, itest);
2802 dd_LPSolve(lp,dd_DualSimplex,&err);
2803 lps=dd_CopyLPSolution(lp);
2804 if (localdebug) dd_WriteLPResult(stdout,lp,err);
2805 if (dd_Positive(lps->optvalue)){
2806 answer=dd_FALSE;
2807 if (localdebug) fprintf(stderr,"==> %ld th point is not strongly redundant.\n",itest);
2808 } else {
2809 answer=dd_TRUE;
2810 if (localdebug) fprintf(stderr,"==> %ld th point is strongly redundant.\n",itest);
2811 }
2812 }
2813 }
2814 dd_FreeLPSolution(lps);
2815 }
2816 _L999:
2817 dd_FreeLPData(lp);
2818 _L99:
2819 return answer;
2820 }
2821
dd_SRedundantRows(dd_MatrixPtr M,dd_ErrorType * error)2822 dd_rowset dd_SRedundantRows(dd_MatrixPtr M, dd_ErrorType *error) /* 093a */
2823 {
2824 dd_rowrange i,m;
2825 dd_colrange d;
2826 dd_rowset redset;
2827 dd_MatrixPtr Mcopy;
2828 dd_Arow cvec; /* certificate */
2829 dd_boolean localdebug=dd_FALSE;
2830
2831 m=M->rowsize;
2832 if (M->representation==dd_Generator){
2833 d=(M->colsize)+1;
2834 } else {
2835 d=M->colsize;
2836 }
2837 Mcopy=dd_MatrixCopy(M);
2838 dd_InitializeArow(d,&cvec);
2839 set_initialize(&redset, m);
2840 for (i=m; i>=1; i--) {
2841 if (dd_SRedundant(Mcopy, i, cvec, error)) {
2842 if (localdebug) printf("dd_SRedundantRows: the row %ld is strongly redundant.\n", i);
2843 set_addelem(redset, i);
2844 dd_MatrixRowRemove(&Mcopy, i);
2845 } else {
2846 if (localdebug) printf("dd_SRedundantRows: the row %ld is not strongly redundant.\n", i);
2847 }
2848 if (*error!=dd_NoError) goto _L99;
2849 }
2850 _L99:
2851 dd_FreeMatrix(Mcopy);
2852 dd_FreeArow(d, cvec);
2853 return redset;
2854 }
2855
dd_RedundantRowsViaShooting(dd_MatrixPtr M,dd_ErrorType * error)2856 dd_rowset dd_RedundantRowsViaShooting(dd_MatrixPtr M, dd_ErrorType *error) /* 092 */
2857 {
2858 /*
2859 For H-representation only and not quite reliable,
2860 especially when floating-point arithmetic is used.
2861 Use the ordinary (slower) method dd_RedundantRows.
2862 */
2863
2864 dd_rowrange i,m, ired, irow=0;
2865 dd_colrange j,k,d;
2866 dd_rowset redset;
2867 dd_rowindex rowflag;
2868 /* ith comp is negative if the ith inequality (i-1 st row) is redundant.
2869 zero if it is not decided.
2870 k > 0 if it is nonredundant and assigned to the (k-1)th row of M1.
2871 */
2872 dd_MatrixPtr M1;
2873 dd_Arow shootdir, cvec=NULL;
2874 dd_LPPtr lp0, lp;
2875 dd_LPSolutionPtr lps;
2876 dd_ErrorType err;
2877 dd_LPSolverType solver=dd_DualSimplex;
2878 dd_boolean localdebug=dd_FALSE;
2879
2880 m=M->rowsize;
2881 d=M->colsize;
2882 M1=dd_CreateMatrix(m,d);
2883 M1->rowsize=0; /* cheat the rowsize so that smaller matrix can be stored */
2884 set_initialize(&redset, m);
2885 dd_InitializeArow(d, &shootdir);
2886 dd_InitializeArow(d, &cvec);
2887
2888 rowflag=(long *)calloc(m+1, sizeof(long));
2889
2890 /* First find some (likely) nonredundant inequalities by Interior Point Find. */
2891 lp0=dd_Matrix2LP(M, &err);
2892 lp=dd_MakeLPforInteriorFinding(lp0);
2893 dd_FreeLPData(lp0);
2894 dd_LPSolve(lp, solver, &err); /* Solve the LP */
2895 lps=dd_CopyLPSolution(lp);
2896
2897 if (dd_Positive(lps->optvalue)){
2898 /* An interior point is found. Use rayshooting to find some nonredundant
2899 inequalities. */
2900 for (j=1; j<d; j++){
2901 for (k=1; k<=d; k++) dd_set(shootdir[k-1], dd_purezero);
2902 dd_set(shootdir[j], dd_one); /* j-th unit vector */
2903 ired=dd_RayShooting(M, lps->sol, shootdir);
2904 if (localdebug) printf("nonredundant row %3ld found by shooting.\n", ired);
2905 if (ired>0 && rowflag[ired]<=0) {
2906 irow++;
2907 rowflag[ired]=irow;
2908 for (k=1; k<=d; k++) dd_set(M1->matrix[irow-1][k-1], M->matrix[ired-1][k-1]);
2909 }
2910
2911 dd_neg(shootdir[j], dd_one); /* negative of the j-th unit vector */
2912 ired=dd_RayShooting(M, lps->sol, shootdir);
2913 if (localdebug) printf("nonredundant row %3ld found by shooting.\n", ired);
2914 if (ired>0 && rowflag[ired]<=0) {
2915 irow++;
2916 rowflag[ired]=irow;
2917 for (k=1; k<=d; k++) dd_set(M1->matrix[irow-1][k-1], M->matrix[ired-1][k-1]);
2918 }
2919 }
2920
2921 M1->rowsize=irow;
2922 if (localdebug) {
2923 printf("The initial nonredundant set is:");
2924 for (i=1; i<=m; i++) if (rowflag[i]>0) printf(" %ld", i);
2925 printf("\n");
2926 }
2927
2928 i=1;
2929 while(i<=m){
2930 if (rowflag[i]==0){ /* the ith inequality is not yet checked */
2931 if (localdebug) fprintf(stderr, "Checking redundancy of %ld th inequality\n", i);
2932 irow++; M1->rowsize=irow;
2933 for (k=1; k<=d; k++) dd_set(M1->matrix[irow-1][k-1], M->matrix[i-1][k-1]);
2934 if (!dd_Redundant(M1, irow, cvec, &err)){
2935 for (k=1; k<=d; k++) dd_sub(shootdir[k-1], cvec[k-1], lps->sol[k-1]);
2936 ired=dd_RayShooting(M, lps->sol, shootdir);
2937 rowflag[ired]=irow;
2938 for (k=1; k<=d; k++) dd_set(M1->matrix[irow-1][k-1], M->matrix[ired-1][k-1]);
2939 if (localdebug) {
2940 fprintf(stderr, "The %ld th inequality is nonredundant for the subsystem\n", i);
2941 fprintf(stderr, "The nonredundancy of %ld th inequality is found by shooting.\n", ired);
2942 }
2943 } else {
2944 if (localdebug) fprintf(stderr, "The %ld th inequality is redundant for the subsystem and thus for the whole.\n", i);
2945 rowflag[i]=-1;
2946 set_addelem(redset, i);
2947 i++;
2948 }
2949 } else {
2950 i++;
2951 }
2952 } /* endwhile */
2953 } else {
2954 /* No interior point is found. Apply the standard LP technique. */
2955 if (localdebug) printf("No interior-point is found and thus the standard LP technique will be used.\n", ired);
2956 redset=dd_RedundantRows(M, error);
2957 }
2958
2959 dd_FreeLPData(lp);
2960 dd_FreeLPSolution(lps);
2961
2962 M1->rowsize=m; M1->colsize=d; /* recover the original sizes */
2963 dd_FreeMatrix(M1);
2964 dd_FreeArow(d, shootdir);
2965 dd_FreeArow(d, cvec);
2966 free(rowflag);
2967 return redset;
2968 }
2969
dd_Matrix2Adjacency(dd_MatrixPtr M,dd_ErrorType * error)2970 dd_SetFamilyPtr dd_Matrix2Adjacency(dd_MatrixPtr M, dd_ErrorType *error) /* 093 */
2971 {
2972 /* This is to generate the (facet) graph of a polyheron (H) V-represented by M using LPs.
2973 Since it does not use the representation conversion, it should work for a large
2974 scale problem.
2975 */
2976 dd_rowrange i,m;
2977 dd_colrange d;
2978 dd_rowset redset;
2979 dd_MatrixPtr Mcopy;
2980 dd_SetFamilyPtr F=NULL;
2981
2982 m=M->rowsize;
2983 d=M->colsize;
2984 if (m<=0 ||d<=0) {
2985 *error=dd_EmptyRepresentation;
2986 goto _L999;
2987 }
2988 Mcopy=dd_MatrixCopy(M);
2989 F=dd_CreateSetFamily(m, m);
2990 for (i=1; i<=m; i++) {
2991 if (!set_member(i, M->linset)){
2992 set_addelem(Mcopy->linset, i);
2993 redset=dd_RedundantRows(Mcopy, error); /* redset should contain all nonadjacent ones */
2994 set_uni(redset, redset, Mcopy->linset); /* all linearity elements should be nonadjacent */
2995 set_compl(F->set[i-1], redset); /* set the adjacency list of vertex i */
2996 set_delelem(Mcopy->linset, i);
2997 set_free(redset);
2998 if (*error!=dd_NoError) goto _L99;
2999 }
3000 }
3001 _L99:
3002 dd_FreeMatrix(Mcopy);
3003 _L999:
3004 return F;
3005 }
3006
dd_Matrix2WeakAdjacency(dd_MatrixPtr M,dd_ErrorType * error)3007 dd_SetFamilyPtr dd_Matrix2WeakAdjacency(dd_MatrixPtr M, dd_ErrorType *error) /* 093a */
3008 {
3009 /* This is to generate the weak-adjacency (facet) graph of a polyheron (H) V-represented by M using LPs.
3010 Since it does not use the representation conversion, it should work for a large
3011 scale problem.
3012 */
3013 dd_rowrange i,m;
3014 dd_colrange d;
3015 dd_rowset redset;
3016 dd_MatrixPtr Mcopy;
3017 dd_SetFamilyPtr F=NULL;
3018
3019 m=M->rowsize;
3020 d=M->colsize;
3021 if (m<=0 ||d<=0) {
3022 *error=dd_EmptyRepresentation;
3023 goto _L999;
3024 }
3025 Mcopy=dd_MatrixCopy(M);
3026 F=dd_CreateSetFamily(m, m);
3027 for (i=1; i<=m; i++) {
3028 if (!set_member(i, M->linset)){
3029 set_addelem(Mcopy->linset, i);
3030 redset=dd_SRedundantRows(Mcopy, error); /* redset should contain all weakly nonadjacent ones */
3031 set_uni(redset, redset, Mcopy->linset); /* all linearity elements should be nonadjacent */
3032 set_compl(F->set[i-1], redset); /* set the adjacency list of vertex i */
3033 set_delelem(Mcopy->linset, i);
3034 set_free(redset);
3035 if (*error!=dd_NoError) goto _L99;
3036 }
3037 }
3038 _L99:
3039 dd_FreeMatrix(Mcopy);
3040 _L999:
3041 return F;
3042 }
3043
3044
dd_ImplicitLinearity(dd_MatrixPtr M,dd_rowrange itest,dd_Arow certificate,dd_ErrorType * error)3045 dd_boolean dd_ImplicitLinearity(dd_MatrixPtr M, dd_rowrange itest, dd_Arow certificate, dd_ErrorType *error)
3046 /* 092 */
3047 {
3048 /* Checks whether the row itest is implicit linearity for the representation.
3049 All linearity rows are not checked and considered non implicit linearity (dd_FALSE).
3050 This code works for both H- and V-representations. A certificate is
3051 given in the case of dd_FALSE, showing a feasible solution x satisfying the itest
3052 strict inequality for H-representation, a hyperplane RHS and normal (x_0, x) that
3053 separates the itest from the rest. More explicitly, the LP to be setup is
3054 the same thing as redundancy case but with maximization:
3055
3056 H-representation
3057 f* = maximize
3058 b_itest + A_itest x
3059 subject to
3060 b_itest + 1 + A_itest x >= 0 (relaxed inequality. This is not necessary but kept for simplicity of the code)
3061 b_{I-itest} + A_{I-itest} x >= 0 (all inequalities except for itest)
3062 b_L + A_L x = 0. (linearity)
3063
3064 V-representation (=separation problem)
3065 f* = maximize
3066 b_itest x_0 + A_itest x
3067 subject to
3068 b_itest x_0 + A_itest x >= -1 (again, this is not necessary but kept for simplicity.)
3069 b_{I-itest} x_0 + A_{I-itest} x >= 0 (all nonlinearity generators except for itest in one side)
3070 b_L x_0 + A_L x = 0. (linearity generators)
3071
3072 Here, the input matrix is considered as (b, A), i.e. b corresponds to the first column of input
3073 and the row indices of input is partitioned into I and L where L is the set of linearity.
3074 In both cases, the itest data is implicit linearity if and only if the optimal value f* is nonpositive.
3075 The certificate has dimension one more for V-representation case.
3076 */
3077
3078 dd_colrange j;
3079 dd_LPPtr lp;
3080 dd_LPSolutionPtr lps;
3081 dd_ErrorType err=dd_NoError;
3082 dd_boolean answer=dd_FALSE,localdebug=dd_FALSE;
3083
3084 *error=dd_NoError;
3085 if (set_member(itest, M->linset)){
3086 if (localdebug) printf("The %ld th row is linearity and redundancy checking is skipped.\n",itest);
3087 goto _L99;
3088 }
3089
3090 /* Create an LP data for redundancy checking */
3091 if (M->representation==dd_Generator){
3092 lp=dd_CreateLP_V_Redundancy(M, itest);
3093 } else {
3094 lp=dd_CreateLP_H_Redundancy(M, itest);
3095 }
3096
3097 lp->objective = dd_LPmax; /* the lp->objective is set by CreateLP* to LPmin */
3098 dd_LPSolve(lp,dd_choiceRedcheckAlgorithm,&err);
3099 if (err!=dd_NoError){
3100 *error=err;
3101 goto _L999;
3102 } else {
3103 lps=dd_CopyLPSolution(lp);
3104
3105 for (j=0; j<lps->d; j++) {
3106 dd_set(certificate[j], lps->sol[j]);
3107 }
3108
3109 if (lps->LPS==dd_Optimal && dd_EqualToZero(lps->optvalue)){
3110 answer=dd_TRUE;
3111 if (localdebug) fprintf(stderr,"==> %ld th data is an implicit linearity.\n",itest);
3112 } else {
3113 answer=dd_FALSE;
3114 if (localdebug) fprintf(stderr,"==> %ld th data is not an implicit linearity.\n",itest);
3115 }
3116 dd_FreeLPSolution(lps);
3117 }
3118 _L999:
3119 dd_FreeLPData(lp);
3120 _L99:
3121 return answer;
3122 }
3123
3124
dd_FreeOfImplicitLinearity(dd_MatrixPtr M,dd_Arow certificate,dd_rowset * imp_linrows,dd_ErrorType * error)3125 int dd_FreeOfImplicitLinearity(dd_MatrixPtr M, dd_Arow certificate, dd_rowset *imp_linrows, dd_ErrorType *error)
3126 /* 092 */
3127 {
3128 /* Checks whether the matrix M constains any implicit linearity at all.
3129 It returns 1 if it is free of any implicit linearity. This means that
3130 the present linearity rows define the linearity correctly. It returns
3131 nonpositive values otherwise.
3132
3133
3134 H-representation
3135 f* = maximize z
3136 subject to
3137 b_I + A_I x - 1 z >= 0
3138 b_L + A_L x = 0 (linearity)
3139 z <= 1.
3140
3141 V-representation (=separation problem)
3142 f* = maximize z
3143 subject to
3144 b_I x_0 + A_I x - 1 z >= 0 (all nonlinearity generators in one side)
3145 b_L x_0 + A_L x = 0 (linearity generators)
3146 z <= 1.
3147
3148 Here, the input matrix is considered as (b, A), i.e. b corresponds to the first column of input
3149 and the row indices of input is partitioned into I and L where L is the set of linearity.
3150 In both cases, any implicit linearity exists if and only if the optimal value f* is nonpositive.
3151 The certificate has dimension one more for V-representation case.
3152 */
3153
3154 dd_LPPtr lp;
3155 dd_rowrange i,m;
3156 dd_colrange j,d1;
3157 dd_ErrorType err=dd_NoError;
3158 dd_Arow cvec; /* certificate for implicit linearity */
3159
3160 int answer=0,localdebug=dd_FALSE;
3161
3162 *error=dd_NoError;
3163 /* Create an LP data for redundancy checking */
3164 if (M->representation==dd_Generator){
3165 lp=dd_CreateLP_V_ImplicitLinearity(M);
3166 } else {
3167 lp=dd_CreateLP_H_ImplicitLinearity(M);
3168 }
3169
3170 dd_LPSolve(lp,dd_choiceRedcheckAlgorithm,&err);
3171 if (err!=dd_NoError){
3172 *error=err;
3173 goto _L999;
3174 } else {
3175
3176 for (j=0; j<lp->d; j++) {
3177 dd_set(certificate[j], lp->sol[j]);
3178 }
3179
3180 if (localdebug) dd_WriteLPResult(stderr,lp,err);
3181
3182 /* *posset contains a set of row indices that are recognized as nonlinearity. */
3183 if (localdebug) {
3184 fprintf(stderr,"==> The following variables are not implicit linearity:\n");
3185 set_fwrite(stderr, lp->posset_extra);
3186 fprintf(stderr,"\n");
3187 }
3188
3189 if (M->representation==dd_Generator){
3190 d1=(M->colsize)+1;
3191 } else {
3192 d1=M->colsize;
3193 }
3194 m=M->rowsize;
3195 dd_InitializeArow(d1,&cvec);
3196 set_initialize(imp_linrows,m);
3197
3198 if (lp->LPS==dd_Optimal){
3199 if (dd_Positive(lp->optvalue)){
3200 answer=1;
3201 if (localdebug) fprintf(stderr,"==> The matrix has no implicit linearity.\n");
3202 } else if (dd_Negative(lp->optvalue)) {
3203 answer=-1;
3204 if (localdebug) fprintf(stderr,"==> The matrix defines the trivial system.\n");
3205 } else {
3206 answer=0;
3207 if (localdebug) fprintf(stderr,"==> The matrix has some implicit linearity.\n");
3208 }
3209 } else {
3210 answer=-2;
3211 if (localdebug) fprintf(stderr,"==> The LP fails.\n");
3212 }
3213 if (answer==0){
3214 /* List the implicit linearity rows */
3215 for (i=m; i>=1; i--) {
3216 if (!set_member(i,lp->posset_extra)) {
3217 if (dd_ImplicitLinearity(M, i, cvec, error)) {
3218 set_addelem(*imp_linrows, i);
3219 if (localdebug) {
3220 fprintf(stderr," row %ld is implicit linearity\n",i);
3221 fprintf(stderr,"\n");
3222 }
3223 }
3224 if (*error!=dd_NoError) goto _L999;
3225 }
3226 }
3227 } /* end of if (answer==0) */
3228 if (answer==-1) {
3229 for (i=m; i>=1; i--) set_addelem(*imp_linrows, i);
3230 } /* all rows are considered implicit linearity */
3231
3232 dd_FreeArow(d1,cvec);
3233 }
3234 _L999:
3235 dd_FreeLPData(lp);
3236
3237 return answer;
3238 }
3239
3240
dd_ImplicitLinearityRows(dd_MatrixPtr M,dd_ErrorType * error)3241 dd_rowset dd_ImplicitLinearityRows(dd_MatrixPtr M, dd_ErrorType *error) /* 092 */
3242 {
3243 dd_colrange d;
3244 dd_rowset imp_linset;
3245 dd_Arow cvec; /* certificate */
3246 int foi;
3247 dd_boolean localdebug=dd_FALSE;
3248
3249 if (M->representation==dd_Generator){
3250 d=(M->colsize)+2;
3251 } else {
3252 d=M->colsize+1;
3253 }
3254
3255 dd_InitializeArow(d,&cvec);
3256 if (localdebug) fprintf(stdout, "\ndd_ImplicitLinearityRows: Check whether the system contains any implicit linearity.\n");
3257 foi=dd_FreeOfImplicitLinearity(M, cvec, &imp_linset, error);
3258 if (localdebug){
3259 switch (foi) {
3260 case 1:
3261 fprintf(stdout, " It is free of implicit linearity.\n");
3262 break;
3263
3264 case 0:
3265 fprintf(stdout, " It is not free of implicit linearity.\n");
3266 break;
3267
3268 case -1:
3269 fprintf(stdout, " The input system is trivial (i.e. the empty H-polytope or the V-rep of the whole space.\n");
3270 break;
3271
3272 default:
3273 fprintf(stdout, " The LP was not solved correctly.\n");
3274 break;
3275
3276 }
3277 }
3278
3279 if (localdebug){
3280 fprintf(stderr, " Implicit linearity rows are:\n");
3281 set_fwrite(stderr,imp_linset);
3282 fprintf(stderr, "\n");
3283 }
3284 dd_FreeArow(d, cvec);
3285 return imp_linset;
3286 }
3287
dd_MatrixCanonicalizeLinearity(dd_MatrixPtr * M,dd_rowset * impl_linset,dd_rowindex * newpos,dd_ErrorType * error)3288 dd_boolean dd_MatrixCanonicalizeLinearity(dd_MatrixPtr *M, dd_rowset *impl_linset,dd_rowindex *newpos,
3289 dd_ErrorType *error) /* 094 */
3290 {
3291 /* This is to recongnize all implicit linearities, and put all linearities at the top of
3292 the matrix. All implicit linearities will be returned by *impl_linset.
3293 */
3294 dd_rowrange rank;
3295 dd_rowset linrows,ignoredrows,basisrows;
3296 dd_colset ignoredcols,basiscols;
3297 dd_rowrange i,k,m;
3298 dd_rowindex newpos1;
3299 dd_boolean success=dd_FALSE;
3300
3301 linrows=dd_ImplicitLinearityRows(*M, error);
3302 if (*error!=dd_NoError) goto _L99;
3303
3304 m=(*M)->rowsize;
3305
3306 set_uni((*M)->linset, (*M)->linset, linrows);
3307 /* add the implicit linrows to the explicit linearity rows */
3308
3309 /* To remove redundancy of the linearity part,
3310 we need to compute the rank and a basis of the linearity part. */
3311 set_initialize(&ignoredrows, (*M)->rowsize);
3312 set_initialize(&ignoredcols, (*M)->colsize);
3313 set_compl(ignoredrows, (*M)->linset);
3314 rank=dd_MatrixRank(*M,ignoredrows,ignoredcols,&basisrows,&basiscols);
3315 set_diff(ignoredrows, (*M)->linset, basisrows);
3316 dd_MatrixRowsRemove2(M,ignoredrows,newpos);
3317
3318 dd_MatrixShiftupLinearity(M,&newpos1);
3319
3320 for (i=1; i<=m; i++){
3321 k=(*newpos)[i];
3322 if (k>0) {
3323 (*newpos)[i]=newpos1[k];
3324 }
3325 }
3326
3327 *impl_linset=linrows;
3328 success=dd_TRUE;
3329 free(newpos1);
3330 set_free(basisrows);
3331 set_free(basiscols);
3332 set_free(ignoredrows);
3333 set_free(ignoredcols);
3334 _L99:
3335 return success;
3336 }
3337
dd_MatrixCanonicalize(dd_MatrixPtr * M,dd_rowset * impl_linset,dd_rowset * redset,dd_rowindex * newpos,dd_ErrorType * error)3338 dd_boolean dd_MatrixCanonicalize(dd_MatrixPtr *M, dd_rowset *impl_linset, dd_rowset *redset,
3339 dd_rowindex *newpos, dd_ErrorType *error) /* 094 */
3340 {
3341 /* This is to find a canonical representation of a matrix *M by
3342 recognizing all implicit linearities and all redundancies.
3343 All implicit linearities will be returned by *impl_linset and
3344 redundancies will be returned by *redset.
3345 */
3346 dd_rowrange i,k,m;
3347 dd_rowindex newpos1,revpos;
3348 dd_rowset redset1;
3349 dd_boolean success=dd_TRUE;
3350
3351 m=(*M)->rowsize;
3352 set_initialize(redset, m);
3353 revpos=(long *)calloc(m+1,sizeof(long));
3354
3355 success=dd_MatrixCanonicalizeLinearity(M, impl_linset, newpos, error);
3356
3357 if (!success) goto _L99;
3358
3359 for (i=1; i<=m; i++){
3360 k=(*newpos)[i];
3361 if (k>0) revpos[k]=i; /* inverse of *newpos[] */
3362 }
3363
3364 success=dd_MatrixRedundancyRemove(M, &redset1, &newpos1, error); /* 094 */
3365
3366 if (!success) goto _L99;
3367
3368 for (i=1; i<=m; i++){
3369 k=(*newpos)[i];
3370 if (k>0) {
3371 (*newpos)[i]=newpos1[k];
3372 if (newpos1[k]<0) (*newpos)[i]=-revpos[-newpos1[k]]; /* update the certificate of its duplicate removal. */
3373 if (set_member(k,redset1)) set_addelem(*redset, i);
3374 }
3375 }
3376
3377 _L99:
3378 set_free(redset1);
3379 free(newpos1);
3380 free(revpos);
3381 return success;
3382 }
3383
3384
dd_ExistsRestrictedFace(dd_MatrixPtr M,dd_rowset R,dd_rowset S,dd_ErrorType * err)3385 dd_boolean dd_ExistsRestrictedFace(dd_MatrixPtr M, dd_rowset R, dd_rowset S, dd_ErrorType *err)
3386 /* 0.94 */
3387 {
3388 /* This function checkes if there is a point that satifies all the constraints of
3389 the matrix M (interpreted as an H-representation) with additional equality contraints
3390 specified by R and additional strict inequality constraints specified by S.
3391 The set S is supposed to be disjoint from both R and M->linset. When it is not,
3392 the set S will be considered as S\(R U M->linset).
3393 */
3394 dd_boolean answer=dd_FALSE;
3395 dd_LPPtr lp=NULL;
3396
3397 /*
3398 printf("\n--- ERF ---\n");
3399 printf("R = "); set_write(R);
3400 printf("S = "); set_write(S);
3401 */
3402
3403 lp=dd_Matrix2Feasibility2(M, R, S, err);
3404
3405 if (*err!=dd_NoError) goto _L99;
3406
3407 /* Solve the LP by cdd LP solver. */
3408 dd_LPSolve(lp, dd_DualSimplex, err); /* Solve the LP */
3409 if (*err!=dd_NoError) goto _L99;
3410 if (lp->LPS==dd_Optimal && dd_Positive(lp->optvalue)) {
3411 answer=dd_TRUE;
3412 }
3413
3414 dd_FreeLPData(lp);
3415 _L99:
3416 return answer;
3417 }
3418
dd_ExistsRestrictedFace2(dd_MatrixPtr M,dd_rowset R,dd_rowset S,dd_LPSolutionPtr * lps,dd_ErrorType * err)3419 dd_boolean dd_ExistsRestrictedFace2(dd_MatrixPtr M, dd_rowset R, dd_rowset S, dd_LPSolutionPtr *lps, dd_ErrorType *err)
3420 /* 0.94 */
3421 {
3422 /* This function checkes if there is a point that satifies all the constraints of
3423 the matrix M (interpreted as an H-representation) with additional equality contraints
3424 specified by R and additional strict inequality constraints specified by S.
3425 The set S is supposed to be disjoint from both R and M->linset. When it is not,
3426 the set S will be considered as S\(R U M->linset).
3427
3428 This function returns a certificate of the answer in terms of the associated LP solutions.
3429 */
3430 dd_boolean answer=dd_FALSE;
3431 dd_LPPtr lp=NULL;
3432
3433 /*
3434 printf("\n--- ERF ---\n");
3435 printf("R = "); set_write(R);
3436 printf("S = "); set_write(S);
3437 */
3438
3439 lp=dd_Matrix2Feasibility2(M, R, S, err);
3440
3441 if (*err!=dd_NoError) goto _L99;
3442
3443 /* Solve the LP by cdd LP solver. */
3444 dd_LPSolve(lp, dd_DualSimplex, err); /* Solve the LP */
3445 if (*err!=dd_NoError) goto _L99;
3446 if (lp->LPS==dd_Optimal && dd_Positive(lp->optvalue)) {
3447 answer=dd_TRUE;
3448 }
3449
3450
3451 (*lps)=dd_CopyLPSolution(lp);
3452 dd_FreeLPData(lp);
3453 _L99:
3454 return answer;
3455 }
3456
dd_FindRelativeInterior(dd_MatrixPtr M,dd_rowset * ImL,dd_rowset * Lbasis,dd_LPSolutionPtr * lps,dd_ErrorType * err)3457 dd_boolean dd_FindRelativeInterior(dd_MatrixPtr M, dd_rowset *ImL, dd_rowset *Lbasis, dd_LPSolutionPtr *lps, dd_ErrorType *err)
3458 /* 0.94 */
3459 {
3460 /* This function computes a point in the relative interior of the H-polyhedron given by M.
3461 Even the representation is V-representation, it simply interprete M as H-representation.
3462 lps returns the result of solving an LP whose solution is a relative interior point.
3463 ImL returns all row indices of M that are implicit linearities, i.e. their inqualities
3464 are satisfied by equality by all points in the polyhedron. Lbasis returns a row basis
3465 of the submatrix of M consisting of all linearities and implicit linearities. This means
3466 that the dimension of the polyhedron is M->colsize - set_card(Lbasis) -1.
3467 */
3468
3469 dd_rowset S;
3470 dd_colset T, Lbasiscols;
3471 dd_boolean success=dd_FALSE;
3472 dd_rowrange i;
3473 dd_colrange rank;
3474
3475
3476 *ImL=dd_ImplicitLinearityRows(M, err);
3477
3478 if (*err!=dd_NoError) goto _L99;
3479
3480 set_initialize(&S, M->rowsize); /* the empty set */
3481 for (i=1; i <=M->rowsize; i++) {
3482 if (!set_member(i, M->linset) && !set_member(i, *ImL)){
3483 set_addelem(S, i); /* all nonlinearity rows go to S */
3484 }
3485 }
3486 if (dd_ExistsRestrictedFace2(M, *ImL, S, lps, err)){
3487 /* printf("a relative interior point found\n"); */
3488 success=dd_TRUE;
3489 }
3490
3491 set_initialize(&T, M->colsize); /* empty set */
3492 rank=dd_MatrixRank(M,S,T,Lbasis,&Lbasiscols); /* the rank of the linearity submatrix of M. */
3493
3494 set_free(S);
3495 set_free(T);
3496 set_free(Lbasiscols);
3497
3498 _L99:
3499 return success;
3500 }
3501
3502
dd_RayShooting(dd_MatrixPtr M,dd_Arow p,dd_Arow r)3503 dd_rowrange dd_RayShooting(dd_MatrixPtr M, dd_Arow p, dd_Arow r)
3504 {
3505 /* 092, find the first inequality "hit" by a ray from an intpt. */
3506 dd_rowrange imin=-1,i,m;
3507 dd_colrange j, d;
3508 dd_Arow vecmin, vec;
3509 mytype min,t1,t2,alpha, t1min;
3510 dd_boolean started=dd_FALSE;
3511 dd_boolean localdebug=dd_FALSE;
3512
3513 m=M->rowsize;
3514 d=M->colsize;
3515 if (!dd_Equal(dd_one, p[0])){
3516 fprintf(stderr, "Warning: RayShooting is called with a point with first coordinate not 1.\n");
3517 dd_set(p[0],dd_one);
3518 }
3519 if (!dd_EqualToZero(r[0])){
3520 fprintf(stderr, "Warning: RayShooting is called with a direction with first coordinate not 0.\n");
3521 dd_set(r[0],dd_purezero);
3522 }
3523
3524 dd_init(alpha); dd_init(min); dd_init(t1); dd_init(t2); dd_init(t1min);
3525 dd_InitializeArow(d,&vecmin);
3526 dd_InitializeArow(d,&vec);
3527
3528 for (i=1; i<=m; i++){
3529 dd_InnerProduct(t1, d, M->matrix[i-1], p);
3530 if (dd_Positive(t1)) {
3531 dd_InnerProduct(t2, d, M->matrix[i-1], r);
3532 dd_div(alpha, t2, t1);
3533 if (!started){
3534 imin=i; dd_set(min, alpha);
3535 dd_set(t1min, t1); /* store the denominator. */
3536 started=dd_TRUE;
3537 if (localdebug) {
3538 fprintf(stderr," Level 1: imin = %ld and min = ", imin);
3539 dd_WriteNumber(stderr, min);
3540 fprintf(stderr,"\n");
3541 }
3542 } else {
3543 if (dd_Smaller(alpha, min)){
3544 imin=i; dd_set(min, alpha);
3545 dd_set(t1min, t1); /* store the denominator. */
3546 if (localdebug) {
3547 fprintf(stderr," Level 2: imin = %ld and min = ", imin);
3548 dd_WriteNumber(stderr, min);
3549 fprintf(stderr,"\n");
3550 }
3551 } else {
3552 if (dd_Equal(alpha, min)) { /* tie break */
3553 for (j=1; j<= d; j++){
3554 dd_div(vecmin[j-1], M->matrix[imin-1][j-1], t1min);
3555 dd_div(vec[j-1], M->matrix[i-1][j-1], t1);
3556 }
3557 if (dd_LexSmaller(vec,vecmin, d)){
3558 imin=i; dd_set(min, alpha);
3559 dd_set(t1min, t1); /* store the denominator. */
3560 if (localdebug) {
3561 fprintf(stderr," Level 3: imin = %ld and min = ", imin);
3562 dd_WriteNumber(stderr, min);
3563 fprintf(stderr,"\n");
3564 }
3565 }
3566 }
3567 }
3568 }
3569 }
3570 }
3571
3572 dd_clear(alpha); dd_clear(min); dd_clear(t1); dd_clear(t2); dd_clear(t1min);
3573 dd_FreeArow(d, vecmin);
3574 dd_FreeArow(d, vec);
3575 return imin;
3576 }
3577
3578 #ifdef GMPRATIONAL
dd_BasisStatusMaximize(dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_Bmatrix T,dd_rowset equalityset,dd_rowrange objrow,dd_colrange rhscol,ddf_LPStatusType LPS,mytype * optvalue,dd_Arow sol,dd_Arow dsol,dd_rowset posset,ddf_colindex nbindex,ddf_rowrange re,ddf_colrange se,dd_colrange * nse,long * pivots,int * found,int * LPScorrect)3579 void dd_BasisStatusMaximize(dd_rowrange m_size,dd_colrange d_size,
3580 dd_Amatrix A,dd_Bmatrix T,dd_rowset equalityset,
3581 dd_rowrange objrow,dd_colrange rhscol,ddf_LPStatusType LPS,
3582 mytype *optvalue,dd_Arow sol,dd_Arow dsol,dd_rowset posset, ddf_colindex nbindex,
3583 ddf_rowrange re,ddf_colrange se, dd_colrange *nse, long *pivots, int *found, int *LPScorrect)
3584 /* This is just to check whether the status LPS of the basis given by
3585 nbindex with extra certificates se or re is correct. It is done
3586 by recomputing the basis inverse matrix T. It does not solve the LP
3587 when the status *LPS is undecided. Thus the input is
3588 m_size, d_size, A, equalityset, LPS, nbindex, re and se.
3589 Other values will be recomputed from scratch.
3590
3591 The main purpose of the function is to verify the correctness
3592 of the result of floating point computation with the GMP rational
3593 arithmetics.
3594 */
3595 {
3596 long pivots0,pivots1,fbasisrank;
3597 dd_rowrange i,is;
3598 dd_colrange s,senew,j;
3599 static dd_rowindex bflag;
3600 static long mlast=0;
3601 static dd_rowindex OrderVector; /* the permutation vector to store a preordered row indices */
3602 unsigned int rseed=1;
3603 mytype val;
3604 dd_colindex nbtemp;
3605 dd_LPStatusType ddlps;
3606 dd_boolean localdebug=dd_FALSE;
3607
3608 if (dd_debug) localdebug=dd_debug;
3609 if (localdebug){
3610 printf("\nEvaluating dd_BasisStatusMaximize:\n");
3611 }
3612 dd_init(val);
3613 nbtemp=(long *) calloc(d_size+1,sizeof(long));
3614 for (i=0; i<= 4; i++) pivots[i]=0;
3615 if (bflag==NULL || mlast!=m_size){
3616 if (mlast!=m_size && mlast>0) {
3617 free(bflag); /* called previously with different m_size */
3618 free(OrderVector);
3619 }
3620 bflag=(long *) calloc(m_size+1,sizeof(long));
3621 OrderVector=(long *)calloc(m_size+1,sizeof(long));
3622 /* initialize only for the first time or when a larger space is needed */
3623 mlast=m_size;
3624 }
3625
3626 /* Initializing control variables. */
3627 dd_ComputeRowOrderVector2(m_size,d_size,A,OrderVector,dd_MinIndex,rseed);
3628
3629 pivots1=0;
3630
3631 dd_ResetTableau(m_size,d_size,T,nbtemp,bflag,objrow,rhscol);
3632
3633 if (localdebug){
3634 printf("\nnbindex:");
3635 for (j=1; j<=d_size; j++) printf(" %ld", nbindex[j]);
3636 printf("\n");
3637 printf("re = %ld, se=%ld\n", re, se);
3638 }
3639
3640 is=nbindex[se];
3641 if (localdebug) printf("se=%ld, is=%ld\n", se, is);
3642
3643 fbasisrank=d_size-1;
3644 for (j=1; j<=d_size; j++){
3645 if (nbindex[j]<0) fbasisrank=fbasisrank-1;
3646 /* fbasisrank=the basis rank computed by floating-point */
3647 }
3648
3649 if (fbasisrank<d_size-1) {
3650 if (localdebug) {
3651 printf("d_size = %ld, the size of basis = %ld\n", d_size, fbasisrank);
3652 printf("dd_BasisStatusMaximize: the size of basis is smaller than d-1.\nIt is safer to run the LP solver with GMP\n");
3653 }
3654 *found=dd_FALSE;
3655 goto _L99;
3656 /* Suspicious case. Rerun the LP solver with GMP. */
3657 }
3658
3659
3660
3661 dd_FindLPBasis2(m_size,d_size,A,T,OrderVector, equalityset,nbindex,bflag,
3662 objrow,rhscol,&s,found,&pivots0);
3663
3664 /* set up the new se column and corresponding variable */
3665 senew=bflag[is];
3666 is=nbindex[senew];
3667 if (localdebug) printf("new se=%ld, is=%ld\n", senew, is);
3668
3669 pivots[4]=pivots0; /*GMP postopt pivots */
3670 dd_statBSpivots+=pivots0;
3671
3672 if (!(*found)){
3673 if (localdebug) {
3674 printf("dd_BasisStatusMaximize: a specified basis DOES NOT exist.\n");
3675 }
3676
3677 goto _L99;
3678 /* No speficied LP basis is found. */
3679 }
3680
3681 if (localdebug) {
3682 printf("dd_BasisStatusMaximize: a specified basis exists.\n");
3683 if (m_size <=100 && d_size <=30)
3684 dd_WriteTableau(stdout,m_size,d_size,A,T,nbindex,bflag);
3685 }
3686
3687 /* Check whether a recomputed basis is of the type specified by LPS */
3688 *LPScorrect=dd_TRUE;
3689 switch (LPS){
3690 case dd_Optimal:
3691 for (i=1; i<=m_size; i++) {
3692 if (i!=objrow && bflag[i]==-1) { /* i is a basic variable */
3693 dd_TableauEntry(&val,m_size,d_size,A,T,i,rhscol);
3694 if (dd_Negative(val)) {
3695 if (localdebug) printf("RHS entry for %ld is negative\n", i);
3696 *LPScorrect=dd_FALSE;
3697 break;
3698 }
3699 } else if (bflag[i] >0) { /* i is nonbasic variable */
3700 dd_TableauEntry(&val,m_size,d_size,A,T,objrow,bflag[i]);
3701 if (dd_Positive(val)) {
3702 if (localdebug) printf("Reduced cost entry for %ld is positive\n", i);
3703 *LPScorrect=dd_FALSE;
3704 break;
3705 }
3706 }
3707 };
3708 break;
3709 case dd_Inconsistent:
3710 for (j=1; j<=d_size; j++){
3711 dd_TableauEntry(&val,m_size,d_size,A,T,re,j);
3712 if (j==rhscol){
3713 if (dd_Nonnegative(val)){
3714 if (localdebug) printf("RHS entry for %ld is nonnegative\n", re);
3715 *LPScorrect=dd_FALSE;
3716 break;
3717 }
3718 } else if (dd_Positive(val)){
3719 if (localdebug) printf("the row entry for(%ld, %ld) is positive\n", re, j);
3720 *LPScorrect=dd_FALSE;
3721 break;
3722 }
3723 };
3724 break;
3725 case dd_DualInconsistent:
3726 for (i=1; i<=m_size; i++){
3727 dd_TableauEntry(&val,m_size,d_size,A,T,i,bflag[is]);
3728 if (i==objrow){
3729 if (dd_Nonpositive(val)){
3730 if (localdebug) printf("Reduced cost entry for %ld is nonpositive\n", bflag[is]);
3731 *LPScorrect=dd_FALSE;
3732 break;
3733 }
3734 } else if (dd_Negative(val)){
3735 if (localdebug) printf("the column entry for(%ld, %ld) is positive\n", i, bflag[is]);
3736 *LPScorrect=dd_FALSE;
3737 break;
3738 }
3739 };
3740 break;
3741 ;
3742 default: break;
3743 }
3744
3745 ddlps=LPSf2LPS(LPS);
3746
3747 dd_SetSolutions(m_size,d_size,A,T,
3748 objrow,rhscol,ddlps,optvalue,sol,dsol,posset,nbindex,re,senew,bflag);
3749 *nse=senew;
3750
3751
3752 _L99:
3753 dd_clear(val);
3754 free(nbtemp);
3755 }
3756
dd_BasisStatusMinimize(dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_Bmatrix T,dd_rowset equalityset,dd_rowrange objrow,dd_colrange rhscol,ddf_LPStatusType LPS,mytype * optvalue,dd_Arow sol,dd_Arow dsol,dd_rowset posset,ddf_colindex nbindex,ddf_rowrange re,ddf_colrange se,dd_colrange * nse,long * pivots,int * found,int * LPScorrect)3757 void dd_BasisStatusMinimize(dd_rowrange m_size,dd_colrange d_size,
3758 dd_Amatrix A,dd_Bmatrix T,dd_rowset equalityset,
3759 dd_rowrange objrow,dd_colrange rhscol,ddf_LPStatusType LPS,
3760 mytype *optvalue,dd_Arow sol,dd_Arow dsol, dd_rowset posset, ddf_colindex nbindex,
3761 ddf_rowrange re,ddf_colrange se,dd_colrange *nse,long *pivots, int *found, int *LPScorrect)
3762 {
3763 dd_colrange j;
3764
3765 for (j=1; j<=d_size; j++) dd_neg(A[objrow-1][j-1],A[objrow-1][j-1]);
3766 dd_BasisStatusMaximize(m_size,d_size,A,T,equalityset, objrow,rhscol,
3767 LPS,optvalue,sol,dsol,posset,nbindex,re,se,nse,pivots,found,LPScorrect);
3768 dd_neg(*optvalue,*optvalue);
3769 for (j=1; j<=d_size; j++){
3770 if (LPS!=dd_Inconsistent) {
3771 /* Inconsistent certificate stays valid for minimization, 0.94e */
3772 dd_neg(dsol[j-1],dsol[j-1]);
3773 }
3774 dd_neg(A[objrow-1][j-1],A[objrow-1][j-1]);
3775 }
3776 }
3777 #endif
3778
3779 /* end of cddlp.c */
3780
3781