polysign.c

   1 #include <assert.h>
   2 #include <barvinok/options.h>
   3 #include <barvinok/util.h>
   4 #include "polysign.h"
   5 #include "config.h"
   6
   7 static int is_unit_row(Value *row, int pos, int len)
   8 {
   9     if (!value_one_p(row[pos]) && !value_mone_p(row[pos]))
  10         return 0;
  11     return First_Non_Zero(row+pos+1, len-(pos+1)) == -1;
  12 }
  13
  14 /* Transform the constraints of P to "standard form".
  15  * In particular, if P is described by
  16  *              A x + b(p) >= 0
  17  * then this function returns a matrix H = A U, A = H Q, such
  18  * that D x' = D Q x >= -b(p), with D a diagonal matrix with
  19  * positive entries.  The calling function can then construct
  20  * the standard form H' x' - I s + b(p) = 0, with H' the rows of H
  21  * that are not positive multiples of unit vectors
  22  * (since those correspond to D x' >= -b(p)).
  23  * The number of rows in H' is returned in *rows_p.
  24  * Optionally returns the matrix that maps the new variables
  25  * back to the old variables x = U x'.
  26  * Note that the rows of H (and P) may be reordered by this function.
  27  */
  28 Matrix *standard_constraints(Polyhedron *P, unsigned nparam, int *rows_p,
  29                              Matrix **T)
  30 {
  31     unsigned nvar = P->Dimension - nparam;
  32     int i, j, d;
  33     int rows;
  34     Matrix *M;
  35     Matrix *H, *U, *Q;
  36
  37     assert(P->NbEq == 0);
  38
  39     /* move constraints only involving parameters down
  40      * and move unit vectors (if there are any) to the right place.
  41      */
  42     for (d = 0, j = P->NbConstraints; d < j; ++d) {
  43         int index;
  44         index = First_Non_Zero(P->Constraint[d]+1, nvar);
  45         if (index != -1) {
  46             if (index != d &&
  47                 is_unit_row(P->Constraint[d]+1, index, nvar)) {
  48                 Vector_Exchange(P->Constraint[d], P->Constraint[index],
  49                                 P->Dimension+2);
  50                 if (index > d)
  51                     --d;
  52             }
  53             continue;
  54         }
  55         while (d < --j && First_Non_Zero(P->Constraint[j]+1, nvar) == -1)
  56             ;
  57         if (d >= j)
  58             break;
  59         Vector_Exchange(P->Constraint[d], P->Constraint[j], P->Dimension+2);
  60     }
  61     M = Matrix_Alloc(d, nvar);
  62     for (j = 0; j < d; ++j)
  63         Vector_Copy(P->Constraint[j]+1, M->p[j], nvar);
  64
  65     neg_left_hermite(M, &H, &Q, &U);
  66     Matrix_Free(M);
  67     Matrix_Free(Q);
  68     if (T)
  69         *T = U;
  70     else
  71         Matrix_Free(U);
  72
  73     /* Rearrange rows such that top of H is lower diagonal and
  74      * compute the number of non (multiple of) unit-vector rows.
  75      */
  76     rows = H->NbRows-nvar;
  77     for (i = 0; i < H->NbColumns; ++i) {
  78         for (j = i; j < H->NbRows; ++j)
  79             if (value_notzero_p(H->p[j][i]))
  80                 break;
  81         if (j != i) {
  82             Vector_Exchange(P->Constraint[i], P->Constraint[j], P->Dimension+2);
  83             Vector_Exchange(H->p[i], H->p[j], H->NbColumns);
  84         }
  85         if (First_Non_Zero(H->p[i], i) != -1)
  86             rows++;
  87     }
  88     *rows_p = rows;
  89
  90     return H;
  91 }
  92
  93 #ifndef HAVE_LIBGLPK
  94 enum order_sign glpk_polyhedron_affine_sign(Polyhedron *D, Matrix *T,
  95                                             struct barvinok_options *options)
  96 {
  97     assert(0);
  98 }
  99
 100 enum lp_result glpk_constraints_opt(Matrix *C, Value *obj, Value denom,
 101                                     enum lp_dir dir, Value *opt)
 102 {
 103     assert(0);
 104 }
 105
 106 enum lp_result glpk_polyhedron_range(Polyhedron *D, Value *obj, Value denom,
 107                                 Value *min, Value *max,
 108                                 struct barvinok_options *options)
 109 {
 110     assert(0);
 111 }
 112 #endif
 113
 114 #ifndef HAVE_LIBCDDGMP
 115 enum order_sign cdd_polyhedron_affine_sign(Polyhedron *D, Matrix *T,
 116                                             struct barvinok_options *options)
 117 {
 118     assert(0);
 119 }
 120
 121 enum order_sign cddf_polyhedron_affine_sign(Polyhedron *D, Matrix *T,
 122                                             struct barvinok_options *options)
 123 {
 124     assert(0);
 125 }
 126
 127 enum lp_result cdd_constraints_opt(Matrix *C, Value *obj, Value denom,
 128                                     enum lp_dir dir, Value *opt)
 129 {
 130     assert(0);
 131 }
 132
 133 enum lp_result cddf_constraints_opt(Matrix *C, Value *obj, Value denom,
 134                                     enum lp_dir dir, Value *opt)
 135 {
 136     assert(0);
 137 }
 138
 139 enum lp_result cdd_polyhedron_range(Polyhedron *D, Value *obj, Value denom,
 140                                 Value *min, Value *max,
 141                                 struct barvinok_options *options)
 142 {
 143     assert(0);
 144 }
 145
 146 enum lp_result cddf_polyhedron_range(Polyhedron *D, Value *obj, Value denom,
 147                                 Value *min, Value *max,
 148                                 struct barvinok_options *options)
 149 {
 150     assert(0);
 151 }
 152 #endif
 153
 154 enum order_sign polyhedron_affine_sign(Polyhedron *D, Matrix *T,
 155                                             struct barvinok_options *options)
 156 {
 157     if (options->lp_solver == BV_LP_POLYLIB)
 158         return PL_polyhedron_affine_sign(D, T, options);
 159     else if (options->lp_solver == BV_LP_GLPK)
 160         return glpk_polyhedron_affine_sign(D, T, options);
 161     else if (options->lp_solver == BV_LP_CDD)
 162         return cdd_polyhedron_affine_sign(D, T, options);
 163     else if (options->lp_solver == BV_LP_CDDF)
 164         return cddf_polyhedron_affine_sign(D, T, options);
 165     else
 166         assert(0);
 167 }
 168
 169 /*
 170  * Optimize (minimize or maximize depending on dir) the affine
 171  * objective function obj (of length dimension+1), with denominator
 172  * denom over the polyhedron specified by the constraints C.
 173  * The result is returned in opt.
 174  */
 175 enum lp_result constraints_opt(Matrix *C, Value *obj, Value denom,
 176                                 enum lp_dir dir, Value *opt,
 177                                 struct barvinok_options *options)
 178 {
 179     if (options->lp_solver == BV_LP_POLYLIB)
 180         return PL_constraints_opt(C, obj, denom, dir, opt, options->MaxRays);
 181     else if (options->lp_solver == BV_LP_GLPK)
 182         return glpk_constraints_opt(C, obj, denom, dir, opt);
 183     else if (options->lp_solver == BV_LP_CDD)
 184         return cdd_constraints_opt(C, obj, denom, dir, opt);
 185     else if (options->lp_solver == BV_LP_CDDF)
 186         return cddf_constraints_opt(C, obj, denom, dir, opt);
 187     else if (options->lp_solver == BV_LP_PIP)
 188         return pip_constraints_opt(C, obj, denom, dir, opt);
 189     else
 190         assert(0);
 191 }
 192
 193 /*
 194  * Optimize (minimize or maximize depending on dir) the affine
 195  * objective function obj (of length dimension+1), with denominator
 196  * denom over the polyhedron specified by P.
 197  * The result is returned in opt.
 198  */
 199 enum lp_result polyhedron_opt(Polyhedron *P, Value *obj, Value denom,
 200                                 enum lp_dir dir, Value *opt,
 201                                 struct barvinok_options *options)
 202 {
 203     if (options->lp_solver == BV_LP_POLYLIB)
 204         return PL_polyhedron_opt(P, obj, denom, dir, opt);
 205     else {
 206         Matrix M;
 207         Polyhedron_Matrix_View(P, &M, P->NbConstraints);
 208         return constraints_opt(&M, obj, denom, dir, opt, options);
 209     }
 210 }
 211
 212 enum lp_result polyhedron_range(Polyhedron *D, Value *obj, Value denom,
 213                                 Value *min, Value *max,
 214                                 struct barvinok_options *options)
 215 {
 216     if (options->lp_solver == BV_LP_POLYLIB)
 217         return PL_polyhedron_range(D, obj, denom, min, max, options);
 218     else if (options->lp_solver == BV_LP_GLPK)
 219         return glpk_polyhedron_range(D, obj, denom, min, max, options);
 220     else if (options->lp_solver == BV_LP_CDD)
 221         return cdd_polyhedron_range(D, obj, denom, min, max, options);
 222     else if (options->lp_solver == BV_LP_CDDF)
 223         return cddf_polyhedron_range(D, obj, denom, min, max, options);
 224     else if (options->lp_solver == BV_LP_PIP)
 225         return pip_polyhedron_range(D, obj, denom, min, max, options);
 226     else
 227         assert(0);
 228 }