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