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 /* Remove all equalities in P and the context C (if not NULL).
  96  * Does not destroy P (or C).
  97  * Returns transformation on the parameters in the Matrix pointed to by CPP
  98  * (unless NULL) and transformation on the variables in the Matrix pointed to
  99  * by CVP (unless NULL).
 100  * Each of these transformation matrices maps the new parameters/variables
 101  * back to the original ones.
 102  * If it can be shown that there is no integer point in P, then
 103  * an empty polyhedron will be returned.
 104  */
 105 int remove_all_equalities(Polyhedron **P, Polyhedron **C, Matrix **CPP, Matrix **CVP,
 106                           unsigned nparam, unsigned MaxRays)
 107 {
 108     Matrix *CV = NULL;
 109     Matrix *CP = NULL;
 110     Polyhedron *Q = *P;
 111     Polyhedron *D = NULL;
 112     Polyhedron *R;
 113     int changed;
 114     Matrix M;
 115
 116     if (C) {
 117         D = *C;
 118         assert(D->Dimension == nparam);
 119     }
 120
 121     if (Q->NbEq == 0 && (!D || D->NbEq == 0))
 122         return 0;
 123
 124     if (D && D->NbEq) {
 125         Polyhedron_Matrix_View(D, &M, D->NbEq);
 126         CP = compress_variables(&M, 0);
 127         transform(&Q, &D, CP, Q != *P, MaxRays);
 128         nparam = CP->NbColumns-1;
 129     }
 130
 131     /* compress_parms doesn't like equalities that only involve parameters */
 132     remove_parameter_equalities(&Q, &D, &CP, &nparam, Q != *P, MaxRays);
 133
 134     if (!emptyQ2(Q) && Q->NbEq) {
 135         Matrix *T;
 136         Polyhedron_Matrix_View(Q, &M, Q->NbEq);
 137         T = compress_parms(&M, nparam);
 138         if (!T) {
 139             Q = replace_by_empty_polyhedron(Q, Q != *P);
 140         } else if (isIdentity(T)) {
 141             Matrix_Free(T);
 142         } else {
 143             transform(&Q, &D, T, Q != *P, MaxRays);
 144             CP = compose_transformations(CP, T);
 145             nparam = CP->NbColumns-1;
 146             remove_parameter_equalities(&Q, &D, &CP, &nparam, Q != *P, MaxRays);
 147         }
 148     }
 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         Polyhedron_Matrix_View(Q, &M, Q->NbEq);
 157         T = compress_variables(&M, nparam);
 158         if (!T) {
 159             Q = replace_by_empty_polyhedron(Q, Q != *P);
 160         } else if (isIdentity(T)) {
 161             Matrix_Free(T);
 162         } else {
 163             R = Polyhedron_Preimage(Q, T, MaxRays);
 164             CV = compose_transformations(CV, T);
 165             if (Q != *P)
 166                 Polyhedron_Free(Q);
 167             Q = R;
 168         }
 169     }
 170
 171     changed = *P != Q;
 172     *P = Q;
 173     if (C)
 174         *C = D;
 175     if (CPP)
 176         *CPP = CP;
 177     else if (CP)
 178         Matrix_Free(CP);
 179     if (CVP)
 180         *CVP = CV;
 181     else if (CV)
 182         Matrix_Free(CV);
 183     return changed;
 184 }