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     for (i = 0; i < Q->NbEq; ++i)
  51         if (First_Non_Zero(Q->Constraint[i]+1, Q->Dimension-nparam) == -1)
  52             break;
  53
  54     return i;
  55 }
  56
  57 static void remove_parameter_equalities(Polyhedron **Q, Polyhedron **D,
  58                                         Matrix **CP, unsigned *nparam,
  59                                         int free, unsigned MaxRays)
  60 {
  61     int i;
  62
  63     /* We need to loop until we can't find any more equalities
  64      * because the transformation may enable a simplification of the
  65      * constraints resulting in new implicit equalities.
  66      */
  67     while ((i = first_parameter_equality(*Q, *nparam)) < (*Q)->NbEq) {
  68         Matrix *M = Matrix_Alloc((*Q)->NbEq, 1+*nparam+1);
  69         Matrix *T;
  70         int n = 0;
  71
  72         for (; i < (*Q)->NbEq; ++i) {
  73             if (First_Non_Zero((*Q)->Constraint[i]+1, (*Q)->Dimension-*nparam) == -1)
  74                 Vector_Copy((*Q)->Constraint[i]+1+(*Q)->Dimension-*nparam,
  75                             M->p[n++]+1, *nparam+1);
  76         }
  77         M->NbRows = n;
  78         T = compress_variables(M, 0);
  79         Matrix_Free(M);
  80         if (!T) {
  81             *Q = replace_by_empty_polyhedron(*Q, free);
  82             break;
  83         } else {
  84             transform(Q, D, T, free, MaxRays);
  85             *nparam = T->NbColumns-1;
  86             *CP = compose_transformations(*CP, T);
  87             free = 1;
  88         }
  89     }
  90 }
  91
  92 /* Remove all equalities in P and the context C (if not NULL).
  93  * Does not destroy P (or C).
  94  * Returns transformation on the parameters in the Matrix pointed to by CPP
  95  * (unless NULL) and transformation on the variables in the Matrix pointed to
  96  * by CVP (unless NULL).
  97  * Each of these transformation matrices maps the new parameters/variables
  98  * back to the original ones.
  99  * If it can be shown that there is no integer point in P, then
 100  * an empty polyhedron will be returned.
 101  */
 102 int remove_all_equalities(Polyhedron **P, Polyhedron **C, Matrix **CPP, Matrix **CVP,
 103                           unsigned nparam, unsigned MaxRays)
 104 {
 105     Matrix *CV = NULL;
 106     Matrix *CP = NULL;
 107     Polyhedron *Q = *P;
 108     Polyhedron *D = NULL;
 109     Polyhedron *R;
 110     int changed;
 111     Matrix M;
 112
 113     if (C) {
 114         D = *C;
 115         assert(D->Dimension == nparam);
 116     }
 117
 118     if (Q->NbEq == 0 && (!D || D->NbEq == 0))
 119         return 0;
 120
 121     if (D && D->NbEq) {
 122         Polyhedron_Matrix_View(D, &M, D->NbEq);
 123         CP = compress_variables(&M, 0);
 124         transform(&Q, &D, CP, Q != *P, MaxRays);
 125         nparam = CP->NbColumns-1;
 126     }
 127
 128     /* compress_parms doesn't like equalities that only involve parameters */
 129     remove_parameter_equalities(&Q, &D, &CP, &nparam, Q != *P, MaxRays);
 130
 131     if (!emptyQ2(Q) && Q->NbEq) {
 132         Matrix *T;
 133         Polyhedron_Matrix_View(Q, &M, Q->NbEq);
 134         T = compress_parms(&M, nparam);
 135         if (!T) {
 136             Q = replace_by_empty_polyhedron(Q, Q != *P);
 137         } else if (isIdentity(T)) {
 138             Matrix_Free(T);
 139         } else {
 140             transform(&Q, &D, T, Q != *P, MaxRays);
 141             CP = compose_transformations(CP, T);
 142             nparam = CP->NbColumns-1;
 143             remove_parameter_equalities(&Q, &D, &CP, &nparam, Q != *P, MaxRays);
 144         }
 145     }
 146
 147     /* We need to loop until we can't find any more equalities
 148      * because the transformation may enable a simplification of the
 149      * constraints resulting in new implicit equalities.
 150      */
 151     while (!emptyQ2(Q) && Q->NbEq) {
 152         Matrix *T;
 153         Polyhedron_Matrix_View(Q, &M, Q->NbEq);
 154         T = compress_variables(&M, nparam);
 155         if (!T) {
 156             Q = replace_by_empty_polyhedron(Q, Q != *P);
 157         } else if (isIdentity(T)) {
 158             Matrix_Free(T);
 159         } else {
 160             R = Polyhedron_Preimage(Q, T, MaxRays);
 161             CV = compose_transformations(CV, T);
 162             if (Q != *P)
 163                 Polyhedron_Free(Q);
 164             Q = R;
 165         }
 166     }
 167
 168     changed = *P != Q;
 169     *P = Q;
 170     if (C)
 171         *C = D;
 172     if (CPP)
 173         *CPP = CP;
 174     else if (CP)
 175         Matrix_Free(CP);
 176     if (CVP)
 177         *CVP = CV;
 178     else if (CV)
 179         Matrix_Free(CV);
 180     return changed;
 181 }