← Back to branch summary

~ubuntu-branches/ubuntu/saucy/gfan/saucy-proposed

« back to all changes in this revision

Viewing changes to .pc/use_system_cdd.patch/lp_cdd.cpp

  • Committer: Package Import Robot
  • Author(s): Cédric Boutillier
  • Date: 2013-07-09 10:44:01 UTC
  • mfrom: (2.1.2 experimental)
  • Revision ID: package-import@ubuntu.com-20130709104401-5q66ozz5j5af0dak
Tags: 0.5+dfsg-3
* Upload to unstable.
* modify remove_failing_tests_on_32bits.patch to replace command of
  0009RenderStairCase test with an empty one instead of deleting it.
* remove lintian override about spelling error

Show diffs side-by-side

added added

removed removed

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