polysign_cdd_template.cc

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