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         if (!CP) {
 144             *C = replace_by_empty_polyhedron(D, D != *C);
 145             return 1;
 146         }
 147         transform(&Q, &D, CP, Q != *P, MaxRays);
 148         nparam = CP->NbColumns-1;
 149     }
 150
 151     /* compress_parms doesn't like equalities that only involve parameters */
 152     remove_parameter_equalities(&Q, &D, &CP, &nparam, Q != *P, MaxRays);
 153
 154     /* We need to loop until we can't find any more equalities
 155      * because the transformation may enable a simplification of the
 156      * constraints resulting in new implicit equalities.
 157      */
 158     while (!emptyQ2(Q) && Q->NbEq) {
 159         Matrix *T;
 160
 161         do {
 162             Polyhedron_Matrix_View(Q, &M, Q->NbEq);
 163             T = compress_parms(&M, nparam);
 164             if (!T) {
 165                 Q = replace_by_empty_polyhedron(Q, Q != *P);
 166             } else if (is_translation(T)) {
 167                 Matrix_Free(T);
 168                 T = NULL;
 169             } else {
 170                 transform(&Q, &D, T, Q != *P, MaxRays);
 171                 CP = compose_transformations(CP, T);
 172                 nparam = CP->NbColumns-1;
 173                 remove_parameter_equalities(&Q, &D, &CP, &nparam, Q != *P, MaxRays);
 174             }
 175         } while (!emptyQ(Q) && Q->NbEq && T);
 176
 177         if (emptyQ(Q) || !Q->NbEq)
 178             break;
 179
 180         Polyhedron_Matrix_View(Q, &M, Q->NbEq);
 181         T = compress_variables(&M, nparam);
 182         if (!T) {
 183             Q = replace_by_empty_polyhedron(Q, Q != *P);
 184         } else if (isIdentity(T)) {
 185             Matrix_Free(T);
 186         } else {
 187             R = Polyhedron_Preimage(Q, T, MaxRays);
 188             CV = compose_transformations(CV, T);
 189             if (Q != *P)
 190                 Polyhedron_Free(Q);
 191             Q = R;
 192         }
 193     }
 194
 195     changed = *P != Q;
 196     *P = Q;
 197     if (C)
 198         *C = D;
 199     if (CPP)
 200         *CPP = CP;
 201     else if (CP)
 202         Matrix_Free(CP);
 203     if (CVP)
 204         *CVP = CV;
 205     else if (CV)
 206         Matrix_Free(CV);
 207     return changed;
 208 }