1 #include "lp_cdd.h"
2 //extern "C"{
3 #ifdef NOCDDPREFIX
4 #include "setoper.h"
5 #include "cdd.h"
6 #include "cdd_f.h"
7 #else
8 #include "cdd/setoper.h"
9 #include "cdd/cdd.h"
10 #include "cdd/cdd_f.h"
11 #endif
12 //}
13 #include "termorder.h"
14 #include "printer.h"
15 #include "log.h"
16
17 //--------------------------------------------------
18 // LpSolverCdd (double precision)
19 //--------------------------------------------------
20
vectorList2Matrix(int n,const IntegerVectorList & g,ddf_ErrorType * Error)21 static ddf_MatrixPtr vectorList2Matrix(int n, const IntegerVectorList &g, ddf_ErrorType *Error)
22 {
23 ddf_MatrixPtr M=NULL;
24 ddf_rowrange m_input,i;
25 ddf_colrange d_input,j;
26 ddf_RepresentationType rep=ddf_Inequality;
27 ddf_boolean found=ddf_FALSE, newformat=ddf_FALSE, successful=ddf_FALSE;
28 char command[ddf_linelenmax], comsave[ddf_linelenmax];
29 ddf_NumberType NT;
30
31 (*Error)=ddf_NoError;
32
33 rep=ddf_Inequality; newformat=ddf_TRUE;
34
35 m_input=g.size();
36 d_input=n+1;//g.begin()->size()+1;
37
38 NT=ddf_Rational;
39
40 M=ddf_CreateMatrix(m_input, d_input);
41 M->representation=rep;
42 M->numbtype=NT;
43
44 IntegerVectorList::const_iterator I=g.begin();
45 for (i = 0; i < m_input; i++) {
46 ddf_set_si(M->matrix[i][0],0);
47 for (j = 1; j < d_input; j++) {
48 ddf_set_si(M->matrix[i][j],(*I)[j-1]);
49 }
50 I++;
51 }
52
53 successful=ddf_TRUE;
54
55 return M;
56 }
57
cddinit()58 static void cddinit()
59 {
60 static bool initialized;
61 if(!initialized)
62 {
63 debug<<"Initialising cdd.\n";
64 ddf_set_global_constants(); /* First, this must be called. */
65 initialized=true;
66 }
67 }
68
isFacet(const IntegerVectorList & g,IntegerVectorList::const_iterator i)69 bool LpSolverCdd::isFacet(const IntegerVectorList &g, IntegerVectorList::const_iterator i)
70 {
71 bool ret;
72 int index;
73 ddf_MatrixPtr M=NULL,M2=NULL,M3=NULL;
74 ddf_colrange d;
75 ddf_ErrorType err=ddf_NoError;
76 ddf_rowset redrows,linrows,ignoredrows, basisrows;
77 ddf_colset ignoredcols, basiscols;
78 // long rank;
79 mytype val;
80 ddf_DataFileType inputfile;
81 FILE *reading=NULL;
82
83
84 cddinit();
85
86 // ddf_init(val);
87
88 M=vectorList2Matrix(i->size(), g, &err);
89
90 if (err!=ddf_NoError) goto _L99;
91
92 d=M->colsize;
93
94 // redrows=ddf_RedundantRows(M, &err);
95 // set_fwrite(Stderr, redrows);
96
97
98 index=0;
99 for(IntegerVectorList::const_iterator J=g.begin();J!=g.end()&&J!=i;J++){index++;}
100 if(index==g.size())assert(0);
101
102
103 static ddf_Arow temp;
104
105 ddf_InitializeArow(i->size()+1,&temp);
106
107 ret= !ddf_Redundant(M,index+1,temp,&err);
108
109 ddf_FreeMatrix(M);
110 ddf_FreeArow(i->size()+1,temp);
111 // ddf_FreeArow(i->size(),temp);
112
113 if (err!=ddf_NoError) goto _L99;
114 return ret;
115 _L99:
116 assert(0);
117 return false;
118
119 }
120
121 static LpSolverCdd theLpSolverCdd;
122
123
124 //--------------------------------------------------
125 // LpSolverCddGmp (exact arithmetics)
126 //--------------------------------------------------
127
128
cddinitGmp()129 static void cddinitGmp()
130 {
131 static bool initialized;
132 if(!initialized)
133 {
134 dd_set_global_constants(); /* First, this must be called. */
135 initialized=true;
136 }
137 }
138
139
vectorList2MatrixGmp(int n,const IntegerVectorList & g,dd_ErrorType * Error)140 static dd_MatrixPtr vectorList2MatrixGmp(int n, const IntegerVectorList &g, dd_ErrorType *Error)
141 {
142 dd_MatrixPtr M=NULL;
143 dd_rowrange m_input,i;
144 dd_colrange d_input,j;
145 dd_RepresentationType rep=dd_Inequality;
146 dd_boolean found=dd_FALSE, newformat=dd_FALSE, successful=dd_FALSE;
147 char command[dd_linelenmax], comsave[dd_linelenmax];
148 dd_NumberType NT;
149
150 (*Error)=dd_NoError;
151
152 rep=dd_Inequality; newformat=dd_TRUE;
153
154 if(n==-1)
155 {
156 assert(g.size());
157 n=g.begin()->size();
158 }
159 m_input=g.size();
160 // d_input=g.begin()->size()+1;
161 d_input=n+1;
162 if(m_input)
163 {
164 if(g.begin()->size()!=n)
165 {
166 AsciiPrinter(Stderr).printVectorList(g);
167 }
168
169
170 assert(g.begin()->size()==n);
171 }
172
173 NT=dd_Rational;
174
175 M=dd_CreateMatrix(m_input, d_input);
176 M->representation=rep;
177 M->numbtype=NT;
178
179 IntegerVectorList::const_iterator I=g.begin();
180 for (i = 0; i < m_input; i++) {
181 dd_set_si(M->matrix[i][0],0);
182 for (j = 1; j < d_input; j++) {
183 dd_set_si(M->matrix[i][j],(*I)[j-1]);
184 }
185 I++;
186 }
187
188 successful=dd_TRUE;
189
190 return M;
191 }
192
193
vectorList2MatrixIncludingFirstColumnGmp(int n,const IntegerVectorList & inequalities,const IntegerVectorList & equations,dd_ErrorType * Error)194 static dd_MatrixPtr vectorList2MatrixIncludingFirstColumnGmp(int n, const IntegerVectorList &inequalities, const IntegerVectorList &equations, dd_ErrorType *Error)
195 {
196 dd_MatrixPtr M=NULL;
197 dd_rowrange m_input,i;
198 dd_colrange d_input,j;
199 dd_RepresentationType rep=dd_Inequality;
200 dd_boolean found=dd_FALSE, newformat=dd_FALSE, successful=dd_FALSE;
201 // char command[dd_linelenmax], comsave[dd_linelenmax];
202 dd_NumberType NT;
203
204 (*Error)=dd_NoError;
205
206 int numberOfEquations=equations.size();
207 int numberOfInequalities=inequalities.size();
208 m_input=numberOfEquations+numberOfInequalities;
209 assert(m_input>0);
210
211 rep=dd_Inequality; newformat=dd_TRUE;
212
213 d_input=n;
214 assert(d_input>0);
215
216 NT=dd_Rational;
217
218 M=dd_CreateMatrix(m_input, d_input);
219 M->representation=rep;
220 M->numbtype=NT;
221
222 IntegerVectorList::const_iterator I=inequalities.begin();
223 for (i = 0; i < numberOfInequalities; i++) {
224 for (j = 0; j < d_input; j++)dd_set_si(M->matrix[i][j],(*I)[j]);
225 I++;
226 }
227 I=equations.begin();
228 for (; i < m_input; i++) {
229 set_addelem(M->linset, i+1);
230 for (j = 0; j < d_input; j++)dd_set_si(M->matrix[i][j],(*I)[j]);
231 I++;
232 }
233
234 successful=dd_TRUE;
235
236 return M;
237 }
238
hasHomogeneousSolution(int n,const IntegerVectorList & inequalities,const IntegerVectorList & equations)239 bool LpSolverCddGmp::hasHomogeneousSolution(int n, const IntegerVectorList &inequalities, const IntegerVectorList &equations)
240 {
241 if(n==0)return false;
242 int nrows=inequalities.size()+equations.size();
243 if(nrows==0)return true;
244
245 dd_LPSolverType solver=dd_DualSimplex;
246 dd_MatrixPtr A=NULL;
247 dd_LPSolutionPtr lps1;
248 dd_ErrorType err=dd_NoError;
249 int ret=0;
250 cddinitGmp();
251
252 dd_LPPtr lp,lp1;
253
254 A=vectorList2MatrixIncludingFirstColumnGmp(n, inequalities, equations, &err);
255
256 if (err!=dd_NoError) goto _L99;
257
258 A->objective=dd_LPmax;
259 lp=dd_Matrix2LP(A, &err);
260 if (err!=dd_NoError) goto _L99;
261
262
263 // dd_WriteMatrix(Stderr,A);
264
265 /* Find an interior point with cdd LP library. */
266
267 lp1=dd_MakeLPforInteriorFinding(lp);
268 dd_LPSolve(lp1,solver,&err);
269 if (err!=dd_NoError) goto _L99;
270
271 // dd_WriteMatrix(Stderr,A);
272 // dd_WriteLP(Stderr,lp);
273 // dd_WriteLP(Stderr,lp1);
274
275 /* Write an interior point. */
276 lps1=dd_CopyLPSolution(lp1);
277 // dd_WriteLPSolution(Stderr,lps1);
278 /* {
279 fprintf(Stderr,"(");
280 for (int j=1; j <(lps1->d); j++) {
281 if(j!=1)fprintf(Stderr,", ");
282 dd_WriteNumber(Stderr,lps1->sol[j]);
283 }
284 fprintf(Stderr,")\n");
285 }*/
286
287 // dd_WriteNumber(Stderr,lps1->optvalue);
288
289 if(dd_Positive(lps1->optvalue))ret=1;
290 if(dd_Negative(lps1->optvalue))ret=-1;
291
292 // fprintf(Stderr,"ret=%i",ret);
293
294 dd_FreeLPData(lp);
295 dd_FreeLPSolution(lps1);
296 dd_FreeLPData(lp1);
297 dd_FreeMatrix(A);
298
299 return ret>=0;
300 _L99:
301 assert(0);
302 return 0;
303 }
304
305
isFacet(const IntegerVectorList & g,IntegerVectorList::const_iterator i)306 bool LpSolverCddGmp::isFacet(const IntegerVectorList &g, IntegerVectorList::const_iterator i)
307 {
308 bool ret;
309 int index;
310 dd_MatrixPtr M=NULL,M2=NULL,M3=NULL;
311 dd_colrange d;
312 dd_ErrorType err=dd_NoError;
313 dd_rowset redrows,linrows,ignoredrows, basisrows;
314 dd_colset ignoredcols, basiscols;
315 // long rank;
316 mytype val;
317 dd_DataFileType inputfile;
318 FILE *reading=NULL;
319
320
321 cddinitGmp();
322
323
324 // dd_init(val);
325
326 M=vectorList2MatrixGmp(i->size(),g, &err);
327
328 if (err!=dd_NoError) goto _L99;
329
330 d=M->colsize;
331
332 // redrows=dd_RedundantRows(M, &err);
333 // set_fwrite(Stderr, redrows);
334
335
336 index=0;
337 for(IntegerVectorList::const_iterator J=g.begin();J!=g.end()&&J!=i;J++){index++;}
338 if(index==g.size())assert(0);
339
340
341 static dd_Arow temp;
342
343 dd_InitializeArow(i->size()+1,&temp);
344
345 ret= !dd_Redundant(M,index+1,temp,&err);
346
347
348 dd_FreeMatrix(M);
349 dd_FreeArow(i->size()+1,temp);
350
351 if (err!=dd_NoError) goto _L99;
352 return ret;
353 _L99:
354 assert(0);
355 return false;
356
357 }
358
staticInteriorPoint(int n,mpq_t * point,const IntegerVectorList & g,bool strictlyPositive,IntegerVector const * equalitySet=0)359 int staticInteriorPoint(int n, mpq_t *point, const IntegerVectorList &g, bool strictlyPositive, IntegerVector const *equalitySet=0)
360 {// copy-paste from testshoot.c in cdd
361 // AsciiPrinter(Stderr).printVectorList(g);
362 // fprintf(Stderr,"strictlyPositive:%i\n",strictlyPositive);
363
364 //if(equalitySet) AsciiPrinter(Stderr).printVector(*equalitySet);
365 dd_LPSolverType solver=dd_DualSimplex;
366 dd_MatrixPtr A=NULL;
367 dd_LPSolutionPtr lps1;
368 dd_ErrorType err=dd_NoError;
369 int ret=0;
370 cddinitGmp();
371
372 dd_LPPtr lp,lp1;
373
374 assert(g.begin()!=g.end());
375 if(strictlyPositive)
376 {
377 int n=g.begin()->size();
378 IntegerVectorList G=g;
379 for(int i=0;i<n;i++)
380 G.push_back(IntegerVector::standardVector(n,i));
381 A=vectorList2MatrixGmp(n,G, &err);
382 }
383 else
384 A=vectorList2MatrixGmp(n, g, &err);
385 if (err!=dd_NoError) goto _L99;
386
387 if(equalitySet)
388 {
389 for(int i=0;i<g.size();i++)
390 if(!(*equalitySet)[i])
391 dd_set_si(A->matrix[i][0],-1);
392
393 assert(g.size()>=equalitySet->size());
394
395 for(int i=0;i<equalitySet->size();i++)
396 if((*equalitySet)[i])set_addelem(A->linset, i+1);
397 }
398
399 A->objective=dd_LPmax;
400 lp=dd_Matrix2LP(A, &err);
401 if (err!=dd_NoError) goto _L99;
402
403
404 // dd_WriteMatrix(Stderr,A);
405
406 /* Find an interior point with cdd LP library. */
407
408 lp1=dd_MakeLPforInteriorFinding(lp);
409 dd_LPSolve(lp1,solver,&err);
410 if (err!=dd_NoError) goto _L99;
411
412 // dd_WriteMatrix(Stderr,A);
413 // dd_WriteLP(Stderr,lp1);
414
415 /* Write an interior point. */
416 lps1=dd_CopyLPSolution(lp1);
417
418 if(dd_Positive(lps1->optvalue))ret=1;
419 if(dd_Negative(lps1->optvalue))ret=-1;
420
421 // fprintf(Stderr,"ret=%i",ret);
422
423 if (ret>=0)
424 // if (ret>0)
425 for (int j=1; j <(lps1->d)-1; j++)
426 mpq_set(point[j-1],lps1->sol[j]);
427
428 dd_FreeLPData(lp);
429 dd_FreeLPSolution(lps1);
430 dd_FreeLPData(lp1);
431 dd_FreeMatrix(A);
432 return ret;
433 _L99:
434 assert(0);
435 return 0;
436 }
437
staticRelativeInteriorPoint(int n,mpq_t * point,const IntegerVectorList & g,bool strictlyPositive,IntegerVector const * equalitySet=0)438 int staticRelativeInteriorPoint(int n, mpq_t *point, const IntegerVectorList &g, bool strictlyPositive, IntegerVector const *equalitySet=0)
439 {// copy-paste from testshoot.c in cdd
440 // AsciiPrinter(Stderr).printVectorList(g);
441 // fprintf(Stderr,"strictlyPositive:%i\n",strictlyPositive);
442
443 //if(equalitySet) AsciiPrinter(Stderr).printVector(*equalitySet);
444 dd_LPSolverType solver=dd_DualSimplex;
445 dd_MatrixPtr A=NULL;
446 dd_LPSolutionPtr lps1;
447 dd_ErrorType err=dd_NoError;
448 int ret=0;
449 cddinitGmp();
450
451 dd_LPPtr lp,lp1;
452
453 // assert(g.begin()!=g.end());
454 if(strictlyPositive)
455 {
456 // int n=g.begin()->size();
457 IntegerVectorList G=g;
458 for(int i=0;i<n;i++)
459 G.push_back(IntegerVector::standardVector(n,i));
460 A=vectorList2MatrixGmp(n, G, &err);
461 }
462 else
463 A=vectorList2MatrixGmp(n, g, &err);
464 if (err!=dd_NoError) goto _L99;
465
466 if(equalitySet)
467 {
468 for(int i=0;i<g.size();i++)
469 if(!(*equalitySet)[i])
470 dd_set_si(A->matrix[i][0],-1);
471
472 assert(g.size()>=equalitySet->size());
473
474 for(int i=0;i<equalitySet->size();i++)
475 if((*equalitySet)[i])set_addelem(A->linset, i+1);
476 }
477
478 A->objective=dd_LPmax;
479 lp=dd_Matrix2LP(A, &err);
480 if (err!=dd_NoError) goto _L99;
481
482
483 // dd_WriteMatrix(Stderr,A);
484
485 /* Find an interior point with cdd LP library. */
486
487 lp1=dd_MakeLPforInteriorFinding(lp);
488 dd_LPSolve(lp1,solver,&err);
489 if (err!=dd_NoError) goto _L99;
490
491 // dd_WriteMatrix(Stderr,A);
492 // dd_WriteLP(Stderr,lp1);
493
494 /* Write an interior point. */
495 lps1=dd_CopyLPSolution(lp1);
496
497 if(dd_Positive(lps1->optvalue))ret=1;
498 if(dd_Negative(lps1->optvalue))ret=-1;
499
500 // fprintf(Stderr,"ret=%i",ret);
501
502 if (ret>=0)
503 // if (ret>0)
504 for (int j=1; j <(lps1->d)-1; j++)
505 mpq_set(point[j-1],lps1->sol[j]);
506
507 dd_FreeLPData(lp);
508 dd_FreeLPSolution(lps1);
509 dd_FreeLPData(lp1);
510 dd_FreeMatrix(A);
511 return ret;
512 _L99:
513 assert(0);
514 return 0;
515 }
516
hasInteriorPoint(const IntegerVectorList & g,bool strictlyPositive,IntegerVector const * equalitySet)517 bool LpSolverCddGmp::hasInteriorPoint(const IntegerVectorList &g, bool strictlyPositive, IntegerVector const *equalitySet)
518 {
519 assert(!g.empty());
520 int n=g.begin()->size();
521 mpq_t *point = new mpq_t [n];
522 for(int i=0;i<n;i++)mpq_init(point[i]);
523
524 int ret=staticInteriorPoint(n, point,g,strictlyPositive,equalitySet);
525
526 // fprintf(Stderr,"%i\n",ret);
527
528 for(int i=0;i<n;i++)mpq_clear(point[i]);
529 delete [] point;
530
531 if(equalitySet)return ret>=0; //THIS NEEDS TO BE FIXED
532 return ret>0;
533 }
534
535
arrayToIntegerVector(mpq_t * point,int n)536 IntegerVector arrayToIntegerVector(mpq_t *point, int n)
537 {
538 IntegerVector ret(n);
539
540 for (int j=0; j <n; j++)
541 {
542 int den;
543 int num;
544
545 if((!mpz_fits_sint_p(mpq_denref(point[j])))||(!mpz_fits_sint_p(mpq_numref(point[j]))))
546 {
547 fprintf(stderr,"INTEGER OVERFLOW IN POLYHEDRAL COMPUTATION\n");
548 // for (int j=0; j <n; j++)
549 // gmp_printf("%40Qd,\n",point[j]);
550 assert(0);
551 }
552 den=mpz_get_si(mpq_denref(point[j]));
553 num=mpz_get_si(mpq_numref(point[j]));
554
555 assert(den==1);
556 ret[j]=num;
557 }
558
559 return ret;
560 }
561
scaleToIntegerVector(mpq_t * point,int n)562 void scaleToIntegerVector(mpq_t *point, int n)
563 {
564 mpz_t lcm;
565 mpz_t gcd;
566 mpz_init_set_ui(lcm, 1);
567 mpz_init_set_ui(gcd, 0);
568
569 for(int j=0;j<n;j++)
570 {
571 mpz_lcm(lcm,lcm,mpq_denref(point[j]));
572 mpz_gcd(gcd,gcd,mpq_numref(point[j]));
573 }
574
575 if(mpz_sgn(gcd)!=0)
576 if(mpz_cmp(gcd,lcm)!=0)
577 {
578 assert(mpz_sgn(gcd)>0);
579 mpq_t scale;
580 mpq_init(scale);
581 mpq_set_den(scale,gcd);
582 mpq_set_num(scale,lcm);
583 mpq_canonicalize(scale); //according to the gmp manual this call is needed, although it slows things down a bit
584 for(int j=0;j<n;j++)
585 {
586 mpq_mul(point[j],point[j],scale);
587 }
588 mpq_clear(scale);
589 }
590 mpz_clear(lcm);
591 mpz_clear(gcd);
592 }
593
594
relativeInteriorPoint(int n,const IntegerVectorList & g,IntegerVector const * equalitySet)595 IntegerVector LpSolverCddGmp::relativeInteriorPoint(int n, const IntegerVectorList &g, IntegerVector const *equalitySet)
596 {
597 // assert(!g.empty());
598 // int n=g.begin()->size();
599 mpq_t *point = new mpq_t [n];
600 for(int i=0;i<n;i++)mpq_init(point[i]);
601
602 int ret=staticRelativeInteriorPoint(n, point,g,false,equalitySet);
603
604 assert(ret>=0);//-- any cone has a relative interior point
605 // if (ret>0){
606 if (ret>=0)
607 {
608 // fprintf(stderr,"TEST1\n");
609 scaleToIntegerVector(point,n);
610 // fprintf(stderr,"TEST2\n");
611 }
612
613
614 IntegerVector result=arrayToIntegerVector(point,n);
615
616
617 for(int i=0;i<n;i++)mpq_clear(point[i]);
618 delete [] point;
619
620 return result;
621 }
622
623
interiorPoint(const IntegerVectorList & g,IntegerVector & result,bool strictlyPositive,IntegerVector const * equalitySet)624 bool LpSolverCddGmp::interiorPoint(const IntegerVectorList &g, IntegerVector &result, bool strictlyPositive, IntegerVector const *equalitySet)
625 {
626 int n=g.begin()->size();
627 mpq_t *point = new mpq_t [n];
628 for(int i=0;i<n;i++)mpq_init(point[i]);
629
630 int ret=staticInteriorPoint(n, point,g,strictlyPositive,equalitySet);
631
632 // if (ret>0){
633 if (ret>=0)
634 {
635 scaleToIntegerVector(point,n);
636 }
637
638
639 result=arrayToIntegerVector(point,n);
640 // if (ret>0){
641 /* if (ret>=0){
642 fprintf(f,"(");
643 for (int j=0; j <n; j++) {
644 if(j!=0)fprintf(f,", ");
645 dd_WriteNumber(f,point[j]);
646 }
647 fprintf(f,")\n");
648 }
649 if (ret<0)
650 fprintf(f,"The feasible region is empty.\n");
651 if (ret==0)
652 fprintf(f,"The feasible region is nonempty but has no interior point.\n");
653 */
654
655
656 for(int i=0;i<n;i++)mpq_clear(point[i]);
657 delete [] point;
658
659 return (ret>=0);
660 }
661
662 // the following two routines are static to avoid including gmp.h in the header file. Maybe that should be changed...
663
lexicographicShootCompare(IntegerVector const & a,IntegerVector const & b,mpq_t const & aDot,mpq_t const & bDot,mpq_t & aTemp,mpq_t & bTemp)664 static bool lexicographicShootCompare(IntegerVector const &a, IntegerVector const &b, mpq_t const &aDot, mpq_t const &bDot, mpq_t &aTemp, mpq_t &bTemp)
665 {
666 int n=a.size();
667 assert(b.size()==n);
668 for(int i=0;i<n;i++)
669 {
670 mpq_set_si(aTemp,a[i],1);
671 mpq_set_si(bTemp,b[i],1);
672 mpq_mul(aTemp,aTemp,bDot);
673 mpq_mul(bTemp,bTemp,aDot);
674 int cmp=mpq_cmp(aTemp,bTemp);
675 if(cmp>0)return true;
676 if(cmp<0)return false;
677 }
678 assert(0);
679 return false;
680 }
681
computeDotProduct(mpq_t & ret,IntegerVector const & iv,mpq_t const * qv,mpq_t & temp)682 static void computeDotProduct(mpq_t &ret, IntegerVector const &iv, mpq_t const *qv, mpq_t &temp)
683 {
684 mpq_set_si(ret,0,1);
685 int n=iv.size();
686 for(int i=0;i<n;i++)
687 {
688 mpq_set_si(temp,iv[i],1);
689 mpq_mul(temp,temp,qv[i]);
690 mpq_add(ret,ret,temp);
691 }
692 }
693
shoot(const IntegerVectorList & g)694 IntegerVectorList::const_iterator LpSolverCddGmp::shoot(const IntegerVectorList &g)
695 {
696 // fprintf(Stderr,"shoot begin\n");
697
698 // AsciiPrinter(Stderr).printVectorList(g);
699
700 int n=g.begin()->size();
701 mpq_t *point = new mpq_t [n];
702 for(int i=0;i<n;i++)mpq_init(point[i]);
703
704 // fprintf(Stderr,"\nVektor list with no interior point:\n");
705 // AsciiPrinter(Stderr).printVectorList(g);
706
707
708 int ret=staticInteriorPoint(n, point,g,true);
709
710 //fprintf(Stderr,"Interior point\n");
711 // AsciiPrinter(Stderr).printIntegerVector(point);
712
713 assert(ret>0);
714
715 IntegerVectorList::const_iterator bestIterator=g.end();
716 mpq_t tempA;
717 mpq_init(tempA);
718 mpq_t tempB;
719 mpq_init(tempB);
720 mpq_t bestDotProduct;
721 mpq_init(bestDotProduct);
722 mpq_t currentDotProduct;
723 mpq_init(currentDotProduct);
724
725 IntegerVectorList::const_iterator i=g.begin();
726 for(;i!=g.end();i++)
727 if(LexicographicTermOrder()(*i,*i-*i))
728 {
729 computeDotProduct(bestDotProduct,*i,point,tempA);
730 bestIterator=i;
731 break;
732 }
733 // fprintf(Stderr,"A\n");
734 //if(bestIterator!=g.end())
735 // AsciiPrinter(Stderr).printVector(*bestIterator);
736 if(i!=g.end())
737 for(i++;i!=g.end();i++)
738 {
739 computeDotProduct(currentDotProduct,*i,point,tempA);
740 if(lexicographicShootCompare(*bestIterator,*i,bestDotProduct,currentDotProduct,tempA,tempB))
741 {
742 bestIterator=i;
743 mpq_set(bestDotProduct,currentDotProduct);
744 }
745 }
746 mpq_clear(currentDotProduct);
747 mpq_clear(bestDotProduct);
748 mpq_clear(tempB);
749 mpq_clear(tempA);
750
751 for(int i=0;i<n;i++)mpq_clear(point[i]);
752 delete [] point;
753
754 // fprintf(Stderr,"shoot end\n");
755 //if(bestIterator!=g.end())
756 // AsciiPrinter(Stderr).printVector(*bestIterator);
757 return bestIterator;
758 }
759
760
761 //-----------------------------------------
762 // Positive vector in kernel
763 //-----------------------------------------
764
positiveVectorInKernel(const IntegerVectorList & g_,IntegerVector * result)765 bool LpSolverCddGmp::positiveVectorInKernel(const IntegerVectorList &g_, IntegerVector *result)
766 {// copy-paste from testshoot.c in cdd
767 // AsciiPrinter(Stderr).printVectorList(g);
768 // fprintf(Stderr,"strictlyPositive:%i\n",strictlyPositive);
769
770
771 IntegerVectorList g=g_;
772
773
774 // AsciiPrinter(Stderr).printVectorList(g);
775
776 int dim2=g.size();
777
778 assert(g.size()!=0);
779 int dimension=g.begin()->size();
780
781 for(int i=0;i<dimension;i++)
782 g.push_back(IntegerVector::standardVector(dimension,i));
783
784 dd_LPSolverType solver=dd_DualSimplex;
785 dd_MatrixPtr A=NULL;
786 dd_LPSolutionPtr lps1;
787 dd_ErrorType err=dd_NoError;
788 // int ret=0;
789 cddinitGmp();
790
791 dd_LPPtr lp;
792
793 A=vectorList2MatrixGmp(dimension, g, &err);
794
795 for(int i=0;i<dim2;i++)set_addelem(A->linset, i+1);//???
796 // set objective function
797 for (int i = 0; i < dimension; i++)
798 dd_set_si(A->rowvec[i+1],-1);
799
800
801 // set right hand side
802 for (int i = 0; i < dimension; i++)
803 dd_set_si(A->matrix[i+dim2][0],-1);
804
805
806 if (err!=dd_NoError) goto _L99;
807 A->objective=dd_LPmax;
808
809
810 /* for(int i=0;i<dimension;i++)
811 {
812 for(int j=0;j<dimension;j++)
813 dd_WriteNumber(f,point[j]);
814
815 fprintf(Stderr,"\n");
816 }
817 */
818
819 // dd_WriteMatrix(Stderr,A);
820 // assert(0);
821
822 lp=dd_Matrix2LP(A, &err);
823 if (err!=dd_NoError) goto _L99;
824
825 // dd_WriteLP(Stderr,lp);
826 /* Find an interior point with cdd LP library. */
827
828 dd_LPSolve(lp,solver,&err);
829
830
831 if (err!=dd_NoError) goto _L99;
832
833 // IF INFEASIBLE RETURN FALSE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
834 if(lp->LPS==dd_Inconsistent)
835 {
836 dd_FreeLPData(lp);
837 dd_FreeMatrix(A);
838
839 return false;
840 }
841 // dd_WriteLPResult(stdout,lp,err);
842
843 /* Write an interior point. */
844 lps1=dd_CopyLPSolution(lp);
845
846 if (dd_Positive(lps1->optvalue))
847 {
848 // fprintf(Stderr,"Returning false\n");
849 return false;
850 }
851
852
853
854
855 /* {
856 fprintf(Stderr,"(");
857 for (int j=1; j <(lps1->d); j++) {
858 if(j!=1)fprintf(Stderr,", ");
859 dd_WriteNumber(Stderr,lps1->sol[j]);
860 }
861 fprintf(Stderr,")\n");
862 }
863 */
864
865 //transform into integer vectors
866 {
867 int n=dimension;
868 dd_Arow point=lps1->sol+1;
869
870 mpz_t lcm;
871 mpz_t gcd;
872 mpz_init_set_ui(lcm, 1);
873 mpz_init_set_ui(gcd, 0);
874
875 for(int j=0;j<n;j++)
876 {
877 mpz_lcm(lcm,lcm,mpq_denref(point[j]));
878 mpz_gcd(gcd,gcd,mpq_numref(point[j]));
879 }
880
881 mpq_t scale;
882 mpq_init(scale);
883 mpq_set_den(scale,gcd);
884 mpq_set_num(scale,lcm);
885 for(int j=0;j<n;j++)
886 {
887 mpq_mul(point[j],point[j],scale);
888 }
889
890 mpz_clear(lcm);
891 mpz_clear(gcd);
892 }
893
894
895
896 /* {
897 fprintf(Stderr,"(");
898 for (int j=1; j <(lps1->d); j++) {
899 if(j!=1)fprintf(Stderr,", ");
900 dd_WriteNumber(Stderr,lps1->sol[j]);
901 }
902 fprintf(Stderr,")\n");
903 }
904 */
905 for (int j=1; j <(lps1->d); j++)
906 {
907 int den;
908 int num;
909
910 den=mpz_get_si(mpq_denref(lps1->sol[j]));
911 num=mpz_get_si(mpq_numref(lps1->sol[j]));
912
913 assert(den==1);
914 if(result)
915 {
916 assert(j-1<result->size());
917 (*result)[j-1]=num;
918 }
919 }
920
921 // dd_WriteLPSolution(lps1);
922
923 dd_FreeLPData(lp);
924 dd_FreeLPSolution(lps1);
925 dd_FreeMatrix(A);
926 return true;
927 _L99:
928 assert(0);
929 return 0;
930 }
931
932 //-----------------------------------------
933 // Rank of matrix
934 //-----------------------------------------
935
936 #include "subspace.h"
rankOfMatrix(const IntegerVectorList & g_)937 int LpSolverCddGmp::rankOfMatrix(const IntegerVectorList &g_)
938 {
939 if(g_.size()==0)return 0;
940 FieldMatrix temp=integerMatrixToFieldMatrix(rowsToIntegerMatrix(g_),Q);
941 return temp.reduceAndComputeRank();
942 return Subspace(g_).dimension();
943 dd_rowset r=NULL;
944 int result;
945 IntegerVectorList g=g_;
946
947 int dim2=g.size();
948
949 assert(g.size()!=0);
950 int dimension=g.begin()->size();
951
952 /* for(int i=0;i<dimension;i++)
953 g.push_back(IntegerVector::standardVector(dimension,i));
954 */
955 dd_LPSolverType solver=dd_DualSimplex;
956 dd_MatrixPtr A=NULL;
957 dd_ErrorType err=dd_NoError;
958 // int ret=0;
959 cddinitGmp();
960
961
962 A=vectorList2MatrixGmp(dimension, g, &err);
963
964 // dd_WriteMatrix(Stderr,A);
965
966 for(int i=0;i<g.size();i++)
967 set_addelem(A->linset,i+1);
968 // for(int i=0;i<dim2;i++)set_addelem(A->linset, i+1);//???
969 // set objective function
970 for (int i = 0; i < dimension; i++)
971 dd_set_si(A->rowvec[i+1],-1);
972
973
974 // set right hand side
975 /*for (int i = 0; i < dimension; i++)
976 dd_set_si(A->matrix[i+dim2][0],-1);
977 */
978
979 if (err!=dd_NoError) goto _L99;
980 A->objective=dd_LPmax;
981
982 //fprintf(Stderr,"rank of matrix matrix:\n");
983
984 // dd_WriteMatrix(Stderr,A);
985
986 r=dd_RedundantRows(A,&err);
987 if (err!=dd_NoError) goto _L99;
988
989 result=dim2-set_card(r);
990
991 // fprintf(Stderr,"dim2==%i set_card(r)==%i\n",dim2,set_card(r));
992
993 set_free(r);
994
995 dd_FreeMatrix(A);
996 return result;
997 _L99:
998 assert(0);
999 return 0;
1000 }
1001
1002 //-----------------------------------------
1003 // Extreme Rays' Inequality Indices
1004 //-----------------------------------------
1005
1006 // this procedure is take from cddio.c.
dd_ComputeAinc(dd_PolyhedraPtr poly)1007 static void dd_ComputeAinc(dd_PolyhedraPtr poly)
1008 {
1009 /* This generates the input incidence array poly->Ainc, and
1010 two sets: poly->Ared, poly->Adom.
1011 */
1012 dd_bigrange k;
1013 dd_rowrange i,m1;
1014 dd_colrange j;
1015 dd_boolean redundant;
1016 dd_MatrixPtr M=NULL;
1017 mytype sum,temp;
1018
1019 dd_init(sum); dd_init(temp);
1020 if (poly->AincGenerated==dd_TRUE) goto _L99;
1021
1022 M=dd_CopyOutput(poly);
1023 poly->n=M->rowsize;
1024 m1=poly->m1;
1025
1026
1027 /* fprintf(Stderr,"m1:%i\n",m1);
1028 fprintf(Stderr,"n:%i\n",poly->n);
1029 fprintf(Stderr,"m:%i\n",poly->m);
1030 */
1031 /* this number is same as poly->m, except when
1032 poly is given by nonhomogeneous inequalty:
1033 !(poly->homogeneous) && poly->representation==Inequality,
1034 it is poly->m+1. See dd_ConeDataLoad.
1035 */
1036 poly->Ainc=(set_type*)calloc(m1, sizeof(set_type));
1037 for(i=1; i<=m1; i++) set_initialize(&(poly->Ainc[i-1]),poly->n);
1038 set_initialize(&(poly->Ared), m1);
1039 set_initialize(&(poly->Adom), m1);
1040
1041 for (k=1; k<=poly->n; k++){
1042 for (i=1; i<=poly->m; i++){
1043 dd_set(sum,dd_purezero);
1044 for (j=1; j<=poly->d; j++){
1045 dd_mul(temp,poly->A[i-1][j-1],M->matrix[k-1][j-1]);
1046 dd_add(sum,sum,temp);
1047 }
1048 if (dd_EqualToZero(sum)) {
1049 set_addelem(poly->Ainc[i-1], k);
1050 }
1051 }
1052 if (!(poly->homogeneous) && poly->representation==dd_Inequality){
1053 if (dd_EqualToZero(M->matrix[k-1][0])) {
1054 set_addelem(poly->Ainc[m1-1], k); /* added infinity inequality (1,0,0,...,0) */
1055 }
1056 }
1057 }
1058
1059 for (i=1; i<=m1; i++){
1060 if (set_card(poly->Ainc[i-1])==M->rowsize){
1061 set_addelem(poly->Adom, i);
1062 }
1063 }
1064 for (i=m1; i>=1; i--){
1065 if (set_card(poly->Ainc[i-1])==0){
1066 redundant=dd_TRUE;
1067 set_addelem(poly->Ared, i);
1068 }else {
1069 redundant=dd_FALSE;
1070 for (k=1; k<=m1; k++) {
1071 if (k!=i && !set_member(k, poly->Ared) && !set_member(k, poly->Adom) &&
1072 set_subset(poly->Ainc[i-1], poly->Ainc[k-1])){
1073 if (!redundant){
1074 redundant=dd_TRUE;
1075 }
1076 set_addelem(poly->Ared, i);
1077 }
1078 }
1079 }
1080 }
1081 dd_FreeMatrix(M);
1082 poly->AincGenerated=dd_TRUE;
1083 _L99:;
1084 dd_clear(sum); dd_clear(temp);
1085 }
1086
1087
extremeRaysInequalityIndices(const IntegerVectorList & inequalityList)1088 IntegerVectorList LpSolverCddGmp::extremeRaysInequalityIndices(const IntegerVectorList &inequalityList)
1089 {
1090 int result;
1091 IntegerVectorList g=inequalityList;
1092
1093 int dim2=g.size();
1094
1095 // assert(g.size()!=0);
1096 if(g.size()==0)return IntegerVectorList();
1097
1098 int dimension=g.begin()->size();
1099
1100 dd_MatrixPtr A=NULL;
1101 dd_ErrorType err=dd_NoError;
1102
1103 cddinitGmp();
1104
1105 A=vectorList2MatrixGmp(dimension, g, &err);
1106
1107 // dd_WriteMatrix(Stderr,A);
1108
1109 dd_PolyhedraPtr poly;
1110 poly=dd_DDMatrix2Poly2(A, dd_LexMin, &err);
1111
1112
1113 // dd_WritePolyFile(Stderr,poly);
1114
1115 // dd_WriteIncidence(Stderr,poly);
1116 if (poly->child==NULL || poly->child->CompStatus!=dd_AllFound) assert(0);
1117 if (poly->AincGenerated==dd_FALSE) dd_ComputeAinc(poly);
1118
1119
1120 // dd_WriteIncidence(Stderr,poly);
1121
1122 IntegerVectorList ret;
1123
1124 // fprintf(Stderr,"%i %i\n",poly->n,poly->m1);
1125
1126 /*
1127 How do we interpret the cddlib output? For a long ting gfan has
1128 been using poly->n as the number of rays of the cone and thus
1129 returned sets of indices that actually gave the lineality
1130 space. The mistake was then caught later in PolyhedralCone. On Feb
1131 17 2009 gfan was changed to check the length of each set to make
1132 sure that it does not specify the lineality space and only return
1133 those sets giving rise to rays. This does not seem to be the best
1134 strategy and might even be wrong.
1135 */
1136
1137
1138 for (int k=1; k<=poly->n; k++)
1139 //for (int k=1; k<=poly->m1; k++)
1140 {
1141 int length=0;
1142 for (int i=1; i<=poly->m1; i++)
1143 if(set_member(k,poly->Ainc[i-1]))length++;
1144 if(length!=dim2)
1145 {
1146 IntegerVector v(length);
1147 int j=0;
1148 for (int i=1; i<=poly->m1; i++)
1149 if(set_member(k,poly->Ainc[i-1]))v[j++]=i-1;
1150 ret.push_back(v);
1151 }
1152 }
1153
1154 // AsciiPrinter(Stderr).printVectorList(ret);
1155 // dd_WriteIncidence(Stderr,poly);
1156
1157 dd_FreeMatrix(A);
1158 dd_FreePolyhedra(poly);
1159
1160 return ret;
1161 _L99:
1162 assert(0);
1163 return IntegerVectorList();
1164 }
1165
1166
1167
1168
1169 //-----------------------------------------
1170 // Remove Redundant Rows
1171 //-----------------------------------------
1172
removeRedundantRows(IntegerVectorList * inequalities,IntegerVectorList * equalities,bool removeInequalityRedundancies)1173 void LpSolverCddGmp::removeRedundantRows(IntegerVectorList *inequalities, IntegerVectorList *equalities, bool removeInequalityRedundancies)
1174 {
1175 int numberOfEqualities=equalities->size();
1176 int numberOfInequalities=inequalities->size();
1177 int numberOfRows=numberOfEqualities+numberOfInequalities;
1178
1179
1180 // AsciiPrinter(Stderr).printVectorList(*inequalities);
1181 // AsciiPrinter(Stderr).printVectorList(*equalities);
1182
1183 dd_rowset r=NULL;
1184 IntegerVectorList g=*inequalities;
1185 for(IntegerVectorList::const_iterator i=equalities->begin();i!=equalities->end();i++)
1186 g.push_back(*i);
1187
1188 if(numberOfRows==0)return;//the full space, so description is alredy irredundant
1189 assert(numberOfRows>0);
1190
1191 dd_LPSolverType solver=dd_DualSimplex;
1192 dd_MatrixPtr A=NULL;
1193 dd_ErrorType err=dd_NoError;
1194
1195 cddinitGmp();
1196
1197 A=vectorList2MatrixGmp(g.begin()->size(),g, &err);
1198
1199 for(int i=numberOfInequalities;i<numberOfRows;i++)
1200 set_addelem(A->linset,i+1);
1201
1202 // dd_WriteMatrix(Stderr,A);
1203
1204 // for(int i=0;i<dim2;i++)set_addelem(A->linset, i+1);//???
1205 // set objective function
1206 // for (int i = 0; i < dimension; i++)
1207 // dd_set_si(A->rowvec[i+1],-1);
1208
1209
1210 // set right hand side
1211 /*for (int i = 0; i < dimension; i++)
1212 dd_set_si(A->matrix[i+dim2][0],-1);
1213 */
1214
1215
1216 IntegerVectorList newLin;
1217 IntegerVectorList newIn;
1218 int index=0;
1219 if (err!=dd_NoError) goto _L99;
1220 A->objective=dd_LPmax;
1221
1222 // dd_WriteMatrix(Stderr,A);
1223
1224 // r=dd_RedundantRows(A,&err);
1225
1226 dd_rowset impl_linset;
1227 dd_rowset redset;
1228 dd_rowindex newpos;
1229
1230 if(removeInequalityRedundancies)
1231 dd_MatrixCanonicalize(&A, &impl_linset, &redset, &newpos, &err);
1232 else
1233 dd_MatrixCanonicalizeLinearity(&A, &impl_linset, &newpos, &err);
1234
1235 if (err!=dd_NoError) goto _L99;
1236
1237 // set_fwrite(stderr,redset);
1238 // set_fwrite(stderr,impl_linset);
1239
1240 // Maybe the following should be changed... what if cononicalize generates new rows?
1241
1242 if(1)
1243 {
1244 // dd_WriteMatrix(Stderr,A);
1245 int rowsize=A->rowsize;
1246 int n=A->colsize-1;
1247 // fprintf(Stderr,"rowsize: %i ColSize:%i\n",rowsize,n);
1248 for(int i=0;i<rowsize;i++)
1249 {
1250 mpq_t *point = new mpq_t [n];
1251 for(int j=0;j<n;j++)mpq_init(point[j]);
1252
1253 for(int j=0;j<n;j++)mpq_set(point[j],A->matrix[i][j+1]);
1254
1255 IntegerVector v=arrayToIntegerVector(point, n);
1256 // AsciiPrinter(Stderr).printVector(v);
1257
1258 for(int j=0;j<n;j++)mpq_clear(point[j]);
1259 delete [] point;
1260 if(set_member(i+1,A->linset))
1261 newLin.push_back(v);
1262 else
1263 newIn.push_back(v);
1264 }
1265 }
1266 else
1267 {
1268 for(IntegerVectorList::iterator i=inequalities->begin();i!=inequalities->end();index++)
1269 {
1270 int i2=newpos[index+1];
1271 if(i2)
1272 {
1273 if(set_member(i2,A->linset))
1274 newLin.push_back(*i);
1275 else
1276 newIn.push_back(*i);
1277 }
1278 i++;
1279 }
1280
1281 for(IntegerVectorList::iterator i=equalities->begin();i!=equalities->end();index++)
1282 {
1283 int i2=newpos[index+1];
1284 if(i2)
1285 {
1286 if(set_member(i2,A->linset))
1287 newLin.push_back(*i);
1288 else
1289 newIn.push_back(*i);
1290 }
1291 i++;
1292 }
1293 assert(index==numberOfRows);
1294 }
1295
1296 assert(set_card(A->linset)==newLin.size());
1297
1298 if(A->rowsize!=newLin.size()+newIn.size())
1299 {
1300 fprintf(stderr,"A->rowsize: %i\n",(int)A->rowsize);
1301 fprintf(stderr,"newLin.size(): %i\n",(int)newLin.size());
1302 fprintf(stderr,"newIn.size(): %i\n",(int)newIn.size());
1303
1304 dd_WriteMatrix(Stderr,A);
1305
1306 AsciiPrinter(Stderr).printVectorList(newLin);
1307 AsciiPrinter(Stderr).printVectorList(newIn);
1308
1309 }
1310 assert(A->rowsize==newLin.size()+newIn.size());
1311
1312
1313 set_free(impl_linset);
1314 if(removeInequalityRedundancies)
1315 set_free(redset);
1316 free(newpos);
1317 // set_free();//HOW DO WE FREE newpos?
1318
1319 *equalities=newLin;
1320 *inequalities=newIn;
1321
1322 dd_FreeMatrix(A);
1323 return;
1324 _L99:
1325 assert(0);
1326 }
1327
vectorLists2MatrixGmp(int n,IntegerVectorList const & inequalities,IntegerVectorList const & equations,dd_ErrorType * err)1328 static dd_MatrixPtr vectorLists2MatrixGmp(int n, IntegerVectorList const &inequalities, IntegerVectorList const &equations, dd_ErrorType *err)
1329 {
1330 IntegerVectorList g=inequalities;
1331 int numberOfInequalities=inequalities.size();
1332 for(IntegerVectorList::const_iterator i=equations.begin();i!=equations.end();i++)
1333 g.push_back(*i);
1334
1335 // assert(g.size()); // If this restriction turns out to be a problem it should be fixed in vectorList2MatrixGmp()
1336 int numberOfRows=g.size();
1337
1338 dd_MatrixPtr A=NULL;
1339
1340 cddinitGmp();
1341
1342 A=vectorList2MatrixGmp(n, g, err);
1343
1344 for(int i=numberOfInequalities;i<numberOfRows;i++)
1345 set_addelem(A->linset,i+1);
1346 return A;
1347 }
1348
1349
getConstraints(dd_MatrixPtr A,bool returnEquations)1350 static IntegerVectorList getConstraints(dd_MatrixPtr A, bool returnEquations)
1351 {
1352 IntegerVectorList ret;
1353
1354 int rowsize=A->rowsize;
1355 int n=A->colsize-1;
1356 // fprintf(Stderr,"rowsize: %i ColSize:%i\n",rowsize,n);
1357
1358 mpq_t *point = new mpq_t [n];
1359 for(int j=0;j<n;j++)mpq_init(point[j]);
1360
1361 for(int i=0;i<rowsize;i++)
1362 {
1363 bool isEquation=set_member(i+1,A->linset);
1364 if(isEquation==returnEquations)
1365 {
1366 for(int j=0;j<n;j++)mpq_set(point[j],A->matrix[i][j+1]);
1367
1368 scaleToIntegerVector(point,n);
1369 IntegerVector v=arrayToIntegerVector(point, n);
1370 // AsciiPrinter(Stderr).printVector(v);
1371
1372 ret.push_back(v);
1373 }
1374 }
1375
1376 for(int j=0;j<n;j++)mpq_clear(point[j]);
1377 delete [] point;
1378
1379 return ret;
1380 }
1381
1382
dual(int n,const IntegerVectorList & inequalities,const IntegerVectorList & equations,IntegerVectorList * dualInequalities,IntegerVectorList * dualEquations)1383 void LpSolverCddGmp::dual(int n, const IntegerVectorList &inequalities, const IntegerVectorList &equations, IntegerVectorList *dualInequalities, IntegerVectorList *dualEquations)
1384 {
1385 IntegerVectorList dummy1;
1386 IntegerVectorList dummy2;
1387 if(!dualInequalities)dualInequalities=&dummy1;
1388 if(!dualEquations)dualEquations=&dummy2;
1389
1390 int result;
1391
1392 dd_MatrixPtr A=NULL;
1393 dd_ErrorType err=dd_NoError;
1394
1395 cddinitGmp();
1396
1397 A=vectorLists2MatrixGmp(n, inequalities, equations, &err);
1398 // A=vectorList2MatrixGmp(g, &err);
1399
1400 // dd_WriteMatrix(Stderr,A);
1401
1402 dd_PolyhedraPtr poly;
1403 poly=dd_DDMatrix2Poly2(A, dd_LexMin, &err);
1404
1405 // dd_WriteIncidence(Stderr,poly);
1406 if (poly->child==NULL || poly->child->CompStatus!=dd_AllFound) assert(0);
1407
1408 // dd_WriteIncidence(Stderr,poly);
1409 // dd_WritePolyFile(Stderr,poly);
1410 // dd_WriteMatrix(Stderr,poly->child->A);
1411 dd_MatrixPtr A2=dd_CopyGenerators(poly);
1412 // dd_WriteMatrix(Stderr,A2);
1413
1414 *dualInequalities=getConstraints(A2,false);
1415 *dualEquations=getConstraints(A2,true);
1416
1417 dd_FreeMatrix(A2);
1418
1419 // AsciiPrinter(Stderr).printVectorList(*dualInequalities);
1420 // AsciiPrinter(Stderr).printVectorList(*dualEquations);
1421
1422
1423 // assert(0);
1424
1425 // AsciiPrinter(Stderr).printVectorList(ret);
1426 // dd_WriteIncidence(Stderr,poly);
1427
1428 dd_FreeMatrix(A);
1429 dd_FreePolyhedra(poly);
1430
1431 return;
1432 _L99:
1433 assert(0);
1434 }
1435
1436 static LpSolverCddGmp theLpSolverCddGmp;
1437