polysign_cdd_template.cc

   1 #define GMPRATIONAL
   2 #include <setoper.h>
   3 #include <cdd.h>
   4 #include <barvinok/util.h>
   5 #include <barvinok/options.h>
   6 #include "polysign.h"
   7 #include "initcdd.h"
   8
   9 #define EMPTY_DOMAIN    -2
  10
  11 static DD_LPType *solve_lp(DD_LPObjectiveType obj, Matrix *C,
  12                            Value *f)
  13 {
  14     unsigned dim = C->NbColumns - 2;
  15     DD_LPType *lp;
  16     DD_rowrange irev = C->NbRows;
  17     DD_rowrange rows, cols;
  18     unsigned NbEq = 0;
  19
  20     for (int i = 0; i < C->NbRows; ++i)
  21         if (value_zero_p(C->p[i][0]))
  22             NbEq++;
  23
  24     rows = irev + NbEq + 1;
  25     cols = 1 + dim;
  26     lp = DD_CreateLPData(obj, DD_Rational, rows, cols);
  27     lp->Homogeneous = DD_FALSE;
  28     lp->objective = obj;
  29
  30     for (DD_rowrange j = 0; j < C->NbRows; ++j) {
  31         for (DD_colrange k = 0; k < dim; ++k)
  32             DD_set_z(lp->A[j][1+k], C->p[j][1+k]);
  33         DD_set_z(lp->A[j][0], C->p[j][1+dim]);
  34         if (value_zero_p(C->p[j][0])) {
  35             set_addelem(lp->equalityset, j+1);
  36             for (DD_colrange k = 0; k < dim; ++k)
  37                 DD_neg(lp->A[irev][1+k], lp->A[j][1+k]);
  38             DD_neg(lp->A[irev][0], lp->A[j][0]);
  39             ++irev;
  40         }
  41     }
  42     /* objective function */
  43     for (DD_colrange k = 0; k < dim; ++k)
  44         DD_set_z(lp->A[rows-1][1+k], f[k]);
  45     DD_set_z(lp->A[rows-1][0], f[dim]);
  46
  47     DD_ErrorType err = DD_NoError;
  48     DD_LPSolve(lp, DD_DualSimplex, &err);
  49     assert(err == DD_NoError);
  50
  51     return lp;
  52 }
  53
  54 static lp_result constraints_affine_minmax(DD_LPObjectiveType obj, Matrix *C,
  55                                           Value *f, Value denom, Value *opt)
  56 {
  57     lp_result res = lp_ok;
  58     DD_LPType *lp = solve_lp(obj, C, f);
  59     assert(value_one_p(denom));
  60
  61     switch(lp->LPS) {
  62     case DD_Optimal:
  63         if (obj == DD_LPmin)
  64             DD_ceil(*opt, lp->optvalue);
  65         else
  66             DD_floor(*opt, lp->optvalue);
  67         break;
  68     case DD_DualInconsistent:
  69         res = lp_unbounded;
  70         break;
  71     case DD_Inconsistent:
  72         res = lp_empty;
  73         break;
  74     default:
  75         assert(0);
  76     }
  77     DD_FreeLPData(lp);
  78     return res;
  79 }
  80
  81 static int constraints_affine_minmax_sign(DD_LPObjectiveType obj, Matrix *C,
  82                                          Matrix *T, bool rational)
  83 {
  84     unsigned dim = C->NbColumns-2;
  85     assert(dim == T->NbColumns-1);
  86     assert(T->NbRows == 2);
  87     DD_LPType *lp = solve_lp(obj, C, T->p[0]);
  88
  89     int sign;
  90     if (lp->LPS == DD_Optimal) {
  91         if (rational)
  92             DD_rat_sign(sign, obj, lp->optvalue);
  93         else
  94             /* The objective function has integer coefficients,
  95              * so the optimal should be an integer (over the integer points)
  96              */
  97             DD_int_sign(sign, obj, lp->optvalue);
  98     } else if (lp->LPS == DD_DualInconsistent) {
  99         if (obj == DD_LPmin)
 100             sign = -1;
 101         else
 102             sign = 1;
 103     } else if (lp->LPS == DD_Inconsistent) {
 104         sign = EMPTY_DOMAIN;
 105     } else
 106         assert(0);
 107
 108     DD_FreeLPData(lp);
 109     return sign;
 110 }
 111
 112 enum order_sign cdd_polyhedron_affine_sign(Polyhedron *D, Matrix *T,
 113                                             struct barvinok_options *options)
 114 {
 115     Matrix M;
 116
 117     if (emptyQ2(D))
 118         return order_undefined;
 119
 120     INIT_CDD;
 121     bool rational = !POL_ISSET(options->MaxRays, POL_INTEGER);
 122     Polyhedron_Matrix_View(D, &M, D->NbConstraints);
 123     int min = constraints_affine_minmax_sign(DD_LPmin, &M, T, rational);
 124     if (min == EMPTY_DOMAIN)
 125         return order_undefined;
 126     if (min > 0)
 127         return order_gt;
 128     int max = constraints_affine_minmax_sign(DD_LPmax, &M, T, rational);
 129     assert(max != EMPTY_DOMAIN);
 130     if (max < 0)
 131         return order_lt;
 132     if (min == max)
 133         return order_eq;
 134     if (max == 0)
 135         return order_le;
 136     if (min == 0)
 137         return order_ge;
 138     return order_unknown;
 139 }
 140
 141 enum lp_result cdd_constraints_opt(Matrix *C, Value *obj, Value denom,
 142                                     enum lp_dir dir, Value *opt)
 143 {
 144     DD_LPObjectiveType cdd_dir = dir == lp_min ? DD_LPmin : DD_LPmax;
 145
 146     INIT_CDD;
 147     return constraints_affine_minmax(cdd_dir, C, obj, denom, opt);
 148 }
 149
 150 enum lp_result cdd_polyhedron_range(Polyhedron *D, Value *obj, Value denom,
 151                                 Value *min, Value *max,
 152                                 struct barvinok_options *options)
 153 {
 154     lp_result res;
 155     Matrix M;
 156
 157     if (emptyQ2(D))
 158         return lp_empty;
 159
 160     INIT_CDD;
 161     Polyhedron_Matrix_View(D, &M, D->NbConstraints);
 162     res = constraints_affine_minmax(DD_LPmin, &M, obj, denom, min);
 163     if (res != lp_ok)
 164         return res;
 165     res = constraints_affine_minmax(DD_LPmax, &M, obj, denom, max);
 166     return res;
 167 }