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