remove_equalities.c

   1 #include <assert.h>
   2 #include <barvinok/util.h>
   3 #include "remove_equalities.h"
   4
   5 static void transform(Polyhedron **P, Polyhedron **C, Matrix *CP, int free,
   6                       unsigned MaxRays)
   7 {
   8     Polyhedron *Q = *P;
   9     Polyhedron *D = *C;
  10     Matrix *T;
  11
  12     T = align_matrix(CP, Q->Dimension+1);
  13     *P = Polyhedron_Preimage(Q, T, MaxRays);
  14     if (free)
  15         Polyhedron_Free(Q);
  16     Matrix_Free(T);
  17     if (!D)
  18         return;
  19     *C = Polyhedron_Preimage(D, CP, MaxRays);
  20     if (free)
  21         Polyhedron_Free(D);
  22 }
  23
  24 static Matrix *compose_transformations(Matrix *first, Matrix *second)
  25 {
  26     Matrix *combined;
  27
  28     if (!first)
  29         return second;
  30
  31     combined = Matrix_Alloc(first->NbRows, second->NbColumns);
  32     Matrix_Product(first, second, combined);
  33     Matrix_Free(first);
  34     Matrix_Free(second);
  35     return combined;
  36 }
  37
  38 static Polyhedron *replace_by_empty_polyhedron(Polyhedron *Q, int free)
  39 {
  40     unsigned dim = Q->Dimension;
  41     if (free)
  42         Polyhedron_Free(Q);
  43     return Empty_Polyhedron(dim);
  44 }
  45
  46 static int first_parameter_equality(Polyhedron *Q, unsigned nparam)
  47 {
  48     int i;
  49
  50     if (emptyQ(Q))
  51         return Q->NbEq;
  52
  53     for (i = 0; i < Q->NbEq; ++i)
  54         if (First_Non_Zero(Q->Constraint[i]+1, Q->Dimension-nparam) == -1)
  55             break;
  56
  57     return i;
  58 }
  59
  60 static void remove_parameter_equalities(Polyhedron **Q, Polyhedron **D,
  61                                         Matrix **CP, unsigned *nparam,
  62                                         int free, unsigned MaxRays)
  63 {
  64     int i;
  65
  66     /* We need to loop until we can't find any more equalities
  67      * because the transformation may enable a simplification of the
  68      * constraints resulting in new implicit equalities.
  69      */
  70     while ((i = first_parameter_equality(*Q, *nparam)) < (*Q)->NbEq) {
  71         Matrix *M = Matrix_Alloc((*Q)->NbEq, 1+*nparam+1);
  72         Matrix *T;
  73         int n = 0;
  74
  75         for (; i < (*Q)->NbEq; ++i) {
  76             if (First_Non_Zero((*Q)->Constraint[i]+1, (*Q)->Dimension-*nparam) == -1)
  77                 Vector_Copy((*Q)->Constraint[i]+1+(*Q)->Dimension-*nparam,
  78                             M->p[n++]+1, *nparam+1);
  79         }
  80         M->NbRows = n;
  81         T = compress_variables(M, 0);
  82         Matrix_Free(M);
  83         if (!T) {
  84             *Q = replace_by_empty_polyhedron(*Q, free);
  85             break;
  86         } else {
  87             transform(Q, D, T, free, MaxRays);
  88             *nparam = T->NbColumns-1;
  89             *CP = compose_transformations(*CP, T);
  90             free = 1;
  91         }
  92     }
  93 }
  94
  95 static int is_translation(Matrix *M)
  96 {
  97     unsigned i, j;
  98
  99     if (M->NbRows != M->NbColumns)
 100         return 0;
 101
 102     for (i = 0; i < M->NbRows; ++i)
 103         for (j = 0; j < M->NbColumns-1; ++j)
 104             if (i == j ? value_notone_p(M->p[i][j])
 105                        : value_notzero_p(M->p[i][j]))
 106                 return 0;
 107
 108     return 1;
 109 }
 110
 111 /* Remove all equalities in P and the context C (if not NULL).
 112  * Does not destroy P (or C).
 113  * Returns transformation on the parameters in the Matrix pointed to by CPP
 114  * (unless NULL) and transformation on the variables in the Matrix pointed to
 115  * by CVP (unless NULL).
 116  * Each of these transformation matrices maps the new parameters/variables
 117  * back to the original ones.
 118  * If it can be shown that there is no integer point in P, then
 119  * an empty polyhedron will be returned.
 120  */
 121 int remove_all_equalities(Polyhedron **P, Polyhedron **C, Matrix **CPP, Matrix **CVP,
 122                           unsigned nparam, unsigned MaxRays)
 123 {
 124     Matrix *CV = NULL;
 125     Matrix *CP = NULL;
 126     Polyhedron *Q = *P;
 127     Polyhedron *D = NULL;
 128     Polyhedron *R;
 129     int changed;
 130     Matrix M;
 131
 132     if (C) {
 133         D = *C;
 134         assert(D->Dimension == nparam);
 135     }
 136
 137     if (Q->NbEq == 0 && (!D || D->NbEq == 0))
 138         return 0;
 139
 140     if (D && D->NbEq) {
 141         Polyhedron_Matrix_View(D, &M, D->NbEq);
 142         CP = compress_variables(&M, 0);
 143         transform(&Q, &D, CP, Q != *P, MaxRays);
 144         nparam = CP->NbColumns-1;
 145     }
 146
 147     /* compress_parms doesn't like equalities that only involve parameters */
 148     remove_parameter_equalities(&Q, &D, &CP, &nparam, Q != *P, MaxRays);
 149
 150     /* We need to loop until we can't find any more equalities
 151      * because the transformation may enable a simplification of the
 152      * constraints resulting in new implicit equalities.
 153      */
 154     while (!emptyQ2(Q) && Q->NbEq) {
 155         Matrix *T;
 156
 157         do {
 158             Polyhedron_Matrix_View(Q, &M, Q->NbEq);
 159             T = compress_parms(&M, nparam);
 160             if (!T) {
 161                 Q = replace_by_empty_polyhedron(Q, Q != *P);
 162             } else if (is_translation(T)) {
 163                 Matrix_Free(T);
 164                 T = NULL;
 165             } else {
 166                 transform(&Q, &D, T, Q != *P, MaxRays);
 167                 CP = compose_transformations(CP, T);
 168                 nparam = CP->NbColumns-1;
 169                 remove_parameter_equalities(&Q, &D, &CP, &nparam, Q != *P, MaxRays);
 170             }
 171         } while (!emptyQ(Q) && Q->NbEq && T);
 172
 173         if (emptyQ(Q) || !Q->NbEq)
 174             break;
 175
 176         Polyhedron_Matrix_View(Q, &M, Q->NbEq);
 177         T = compress_variables(&M, nparam);
 178         if (!T) {
 179             Q = replace_by_empty_polyhedron(Q, Q != *P);
 180         } else if (isIdentity(T)) {
 181             Matrix_Free(T);
 182         } else {
 183             R = Polyhedron_Preimage(Q, T, MaxRays);
 184             CV = compose_transformations(CV, T);
 185             if (Q != *P)
 186                 Polyhedron_Free(Q);
 187             Q = R;
 188         }
 189     }
 190
 191     changed = *P != Q;
 192     *P = Q;
 193     if (C)
 194         *C = D;
 195     if (CPP)
 196         *CPP = CP;
 197     else if (CP)
 198         Matrix_Free(CP);
 199     if (CVP)
 200         *CVP = CV;
 201     else if (CV)
 202         Matrix_Free(CV);
 203     return changed;
 204 }