normalization.c

   1 #include <assert.h>
   2 #include <barvinok/util.h>
   3 #include "normalization.h"
   4
   5 static int is_unit_row(Value *row, int pos, int len)
   6 {
   7     if (!value_one_p(row[pos]) && !value_mone_p(row[pos]))
   8         return 0;
   9     return First_Non_Zero(row+pos+1, len-(pos+1)) == -1;
  10 }
  11
  12 /* Transform the constraints C to "standard form".
  13  * In particular, if C is described by
  14  *              A x + b(p) >= 0
  15  * then this function returns a matrix H = A U, A = H Q, such
  16  * that D x' = D Q x >= -b(p), with D a diagonal matrix with
  17  * positive entries.  The calling function can then construct
  18  * the standard form H' x' - I s + b(p) = 0, with H' the rows of H
  19  * that are not positive multiples of unit vectors
  20  * (since those correspond to D x' >= -b(p)).
  21  * The number of rows in H' is returned in *rows_p.
  22  * Optionally returns the matrix that maps the new variables
  23  * back to the old variables x = U x'.
  24  * Note that the rows of H (and C) may be reordered by this function.
  25  */
  26 Matrix *standard_constraints(Matrix *C, unsigned nparam, int *rows_p,
  27                              Matrix **T)
  28 {
  29     unsigned dim = C->NbColumns - 2;
  30     unsigned nvar = dim - nparam;
  31     int i, j, d;
  32     int rows;
  33     Matrix *M;
  34     Matrix *H, *U, *Q;
  35
  36     /* move constraints only involving parameters down
  37      * and move unit vectors (if there are any) to the right place.
  38      */
  39     for (d = 0, j = C->NbRows; d < j; ++d) {
  40         int index;
  41         index = First_Non_Zero(C->p[d]+1, nvar);
  42         if (index != -1) {
  43             if (index != d &&
  44                 is_unit_row(C->p[d]+1, index, nvar)) {
  45                 Vector_Exchange(C->p[d], C->p[index], dim+2);
  46                 if (index > d)
  47                     --d;
  48             }
  49             continue;
  50         }
  51         while (d < --j && First_Non_Zero(C->p[j]+1, nvar) == -1)
  52             ;
  53         if (d >= j)
  54             break;
  55         Vector_Exchange(C->p[d], C->p[j], dim+2);
  56     }
  57     M = Matrix_Alloc(d, nvar);
  58     for (j = 0; j < d; ++j)
  59         Vector_Copy(C->p[j]+1, M->p[j], nvar);
  60
  61     neg_left_hermite(M, &H, &Q, &U);
  62     Matrix_Free(M);
  63     Matrix_Free(Q);
  64     if (T)
  65         *T = U;
  66     else
  67         Matrix_Free(U);
  68
  69     /* Rearrange rows such that top of H is lower diagonal and
  70      * compute the number of non (multiple of) unit-vector rows.
  71      */
  72     rows = H->NbRows-nvar;
  73     for (i = 0; i < H->NbColumns; ++i) {
  74         for (j = i; j < H->NbRows; ++j)
  75             if (value_notzero_p(H->p[j][i]))
  76                 break;
  77         if (j != i) {
  78             Vector_Exchange(C->p[i], C->p[j], dim+2);
  79             Vector_Exchange(H->p[i], H->p[j], H->NbColumns);
  80         }
  81         if (First_Non_Zero(H->p[i], i) != -1)
  82             rows++;
  83     }
  84     if (rows_p)
  85         *rows_p = rows;
  86
  87     return H;
  88 }
  89
  90 /*
  91  * Transform the constraints of PE involving existentially
  92  * quantified variables, placed in front in PE, to "standard form".
  93  *
  94  * The resulting matrix has the existential variables after
  95  * the unquantified variables again.
  96  */
  97 static Matrix *exist_standard_constraints(Polyhedron *PE, unsigned nvar)
  98 {
  99     int i, i1, i2;
 100     unsigned exist = PE->Dimension - nvar;
 101     Matrix *M = Matrix_Alloc(PE->NbConstraints, 2+PE->Dimension);
 102     Matrix *M2 = Matrix_Alloc(PE->NbConstraints, 2+PE->Dimension);
 103     Matrix *S2;
 104     Value c;
 105
 106     for (i = i1 = i2 = 0; i < PE->NbConstraints; ++i) {
 107         if (First_Non_Zero(PE->Constraint[i]+1, exist) == -1) {
 108             value_assign(M->p[i1][0], PE->Constraint[i][0]);
 109             Vector_Copy(PE->Constraint[i]+1+exist, M->p[i1]+1, nvar);
 110             value_assign(M->p[i1][1+PE->Dimension],
 111                          PE->Constraint[i][1+PE->Dimension]);
 112             ++i1;
 113         } else {
 114             Vector_Copy(PE->Constraint[i], M2->p[i2], M2->NbColumns);
 115             ++i2;
 116         }
 117     }
 118     M->NbRows = i1;
 119     M2->NbRows = i2;
 120
 121     S2 = standard_constraints(M2, nvar, NULL, NULL);
 122
 123     for (i = 0; i < M2->NbRows; ++i) {
 124         value_assign(M->p[M->NbRows+i][0], M2->p[i][0]);
 125         Vector_Copy(M2->p[i]+1+exist, M->p[M->NbRows+i]+1, nvar);
 126         Vector_Copy(S2->p[i], M->p[M->NbRows+i]+1+nvar, exist);
 127         value_assign(M->p[M->NbRows+i][1+PE->Dimension],
 128                         M2->p[i][1+PE->Dimension]);
 129     }
 130
 131     /* Make sure "standard" constraints for existential variables
 132      * have only non-positive coefficients for the variables.
 133      */
 134     assert(M2->NbRows >= exist);
 135     value_init(c);
 136     for (i = 0; i < exist; ++i) {
 137         assert(value_pos_p(S2->p[i][i]));
 138         for (i1 = 0; i1 < nvar; ++i1) {
 139             if (value_negz_p(M->p[M->NbRows+i][1+i1]))
 140                 continue;
 141             mpz_cdiv_q(c, M->p[M->NbRows+i][1+i1], S2->p[i][i]);
 142             value_oppose(c, c);
 143             for (i2 = i; i2 < M2->NbRows; ++i2)
 144                 value_addmul(M->p[M->NbRows+i2][1+i1], c, S2->p[i2][i]);
 145         }
 146     }
 147     value_clear(c);
 148     M->NbRows += M2->NbRows;
 149
 150     Matrix_Free(M2);
 151
 152     return M;
 153 }
 154
 155 /*
 156  * For those variables x that may attain negative values,
 157  * compute a translations t such that x' = x + t is sure
 158  * to attain only non-negative values, i.e., if M is
 159  *
 160  *              A x + b >= 0
 161  *
 162  * then the x' in
 163  *
 164  *              A x' + b - A t >= 0
 165  *
 166  * are sure to attain only non-negative values.
 167  *
 168  * The first nvar shifts have already been set prior to
 169  * calling this function.
 170  */
 171 static void compute_shifts(Matrix *M, unsigned nvar, Vector *shifts)
 172 {
 173     int c, r;
 174     unsigned dim = M->NbColumns-2;
 175     for (c = nvar, r = 0; c < dim; ++c) {
 176         for ( ; r < M->NbRows; ++r) {
 177             if (value_notzero_p(M->p[r][1+c]))
 178                 break;
 179         }
 180         assert(r < M->NbRows);
 181         assert(value_pos_p(M->p[r][1+c]));
 182         if (c > 0)
 183             Inner_Product(M->p[r]+1, shifts->p+1, c, &shifts->p[1+c]);
 184         value_subtract(shifts->p[1+c], M->p[r][1+dim], shifts->p[1+c]);
 185         if (value_pos_p(shifts->p[1+c]))
 186             mpz_cdiv_q(shifts->p[1+c], shifts->p[1+c], M->p[r][1+c]);
 187         else
 188             value_set_si(shifts->p[1+c], 0);
 189     }
 190     return shifts;
 191 }
 192
 193 /*
 194  * Move the existential variables before the unquantified variables.
 195  */
 196 static Polyhedron *move_exists_in_front(Polyhedron *D, unsigned nvar,
 197                                         unsigned MaxRays)
 198 {
 199     Polyhedron *P;
 200     Polyhedron *DE = NULL;
 201     Polyhedron **next = &DE;
 202
 203     for (P = D; P; P = P->next) {
 204         int i;
 205         unsigned exist = P->Dimension - nvar;
 206         Matrix *M = Matrix_Alloc(P->NbConstraints, 2+P->Dimension);
 207
 208         for (i = 0; i < P->NbConstraints; ++i) {
 209             value_assign(M->p[i][0], P->Constraint[i][0]);
 210             Vector_Copy(P->Constraint[i]+1, M->p[i]+1+exist, nvar);
 211             Vector_Copy(P->Constraint[i]+1+nvar, M->p[i]+1, exist);
 212             value_assign(M->p[i][1+P->Dimension],
 213                          P->Constraint[i][1+P->Dimension]);
 214         }
 215
 216         *next = Constraints2Polyhedron(M, MaxRays);
 217         Matrix_Free(M);
 218
 219         next = &(*next)->next;
 220     }
 221
 222     return DE;
 223 }
 224
 225 /*
 226  * Compute the shifts s for the variables x, such that x + s >= 0.
 227  * DE is a union of polyhedra with the existential variables first.
 228  *
 229  * We simply project out the existential variables, compute
 230  * the convex hull and then compute the minimal coordinate
 231  * in each direction of all the vertices.
 232  * The shift is the opposite of this minimal coordinate or zero.
 233  */
 234 static Vector *compute_var_shifts(Polyhedron *DE, unsigned nvar,
 235                                         unsigned MaxRays)
 236 {
 237     Polyhedron *PE;
 238     Polyhedron *D = NULL;
 239     Polyhedron **next = &D;
 240     Polyhedron *C;
 241     Vector *shifts;
 242     int i, j;
 243     Value tmp;
 244
 245     for (PE = DE; PE; PE = PE->next) {
 246         *next = Polyhedron_Project(PE, nvar);
 247         next = &(*next)->next;
 248     }
 249     C = DomainConvex(D, MaxRays);
 250     Domain_Free(D);
 251     POL_ENSURE_VERTICES(C);
 252
 253     shifts = Vector_Alloc(C->Dimension);
 254     value_init(tmp);
 255     for (i = 0; i < C->Dimension; ++i) {
 256         for (j = 0; j < C->NbRays; ++j) {
 257             assert(value_notzero_p(C->Ray[j][1+C->Dimension]));
 258             if (value_posz_p(C->Ray[j][1+i]))
 259                 continue;
 260             mpz_fdiv_q(tmp, C->Ray[j][1+i], C->Ray[j][1+C->Dimension]);
 261             if (value_lt(tmp, shifts->p[i]))
 262                 value_assign(shifts->p[i], tmp);
 263         }
 264         value_oppose(shifts->p[i], shifts->p[i]);
 265     }
 266     value_clear(tmp);
 267     Polyhedron_Free(C);
 268
 269     return shifts;
 270 }
 271
 272 /*
 273  * Given a union of polyhedra, with possibly existentially
 274  * quantified variables, but no parameters, apply a unimodular
 275  * transformation with constant translation to the variables,
 276  * and a unimodular transformation with a translation possibly
 277  * depending on the variables to the existentially quantified variables,
 278  * such that all variables and all existentially quantified variables
 279  * attain only non-negative values.
 280  *
 281  * The projection onto the variables is assumed to be a polytope.
 282  */
 283 Polyhedron *skew_to_positive_orthant(Polyhedron *D, unsigned nvar,
 284                                         unsigned MaxRays)
 285 {
 286     Polyhedron *PE;
 287     Polyhedron *R = NULL;
 288     Polyhedron **next = &R;
 289     Vector *var_shifts;
 290     Polyhedron *DE;
 291
 292     DE = move_exists_in_front(D, nvar, MaxRays);
 293     var_shifts = compute_var_shifts(DE, nvar, MaxRays);
 294
 295     for (PE = DE; PE; PE = PE->next) {
 296         Matrix *M, *C;
 297         Vector *shifts, *offsets;
 298         int i;
 299
 300         M = exist_standard_constraints(PE, nvar);
 301         shifts = Vector_Alloc(M->NbColumns);
 302         Vector_Copy(var_shifts->p, shifts->p+1, nvar);
 303         compute_shifts(M, nvar, shifts);
 304
 305         offsets = Vector_Alloc(M->NbRows);
 306         Matrix_Vector_Product(M, shifts->p, offsets->p);
 307         Vector_Free(shifts);
 308
 309         C = Matrix_Alloc(M->NbRows, M->NbColumns);
 310         for (i = 0; i < M->NbRows; ++i) {
 311             value_assign(C->p[i][0], M->p[i][0]);
 312             Vector_Copy(M->p[i]+1, C->p[i]+1, M->NbColumns-2);
 313             value_subtract(C->p[i][M->NbColumns-1],
 314                             M->p[i][M->NbColumns-1], offsets->p[i]);
 315         }
 316         Vector_Free(offsets);
 317         Matrix_Free(M);
 318
 319         *next = Constraints2Polyhedron(C, MaxRays);
 320         Matrix_Free(C);
 321
 322         next = &(*next)->next;
 323     }
 324     Vector_Free(var_shifts);
 325     Domain_Free(DE);
 326     return R;
 327 }