remove_equalities.c

   1 #include <barvinok/util.h>
   2 #include "remove_equalities.h"
   3
   4 static void transform(Polyhedron **P, Polyhedron **C, Matrix *CP, int free,
   5                       unsigned MaxRays)
   6 {
   7     Polyhedron *Q = *P;
   8     Polyhedron *D = *C;
   9     Matrix *T;
  10
  11     T = align_matrix(CP, Q->Dimension+1);
  12     *P = Polyhedron_Preimage(Q, T, MaxRays);
  13     if (free)
  14         Polyhedron_Free(Q);
  15     Matrix_Free(T);
  16     if (!D)
  17         return;
  18     *C = Polyhedron_Preimage(D, CP, MaxRays);
  19     if (free)
  20         Polyhedron_Free(D);
  21 }
  22
  23 /* Remove all equalities in P and the context C (if not NULL).
  24  * Does not destroy P (or C).
  25  * Returns transformation on the parameters in the Matrix pointed to by CPP
  26  * (unless NULL) and transformation on the variables in the Matrix pointed to
  27  * by CVP (unless NULL).
  28  * Each of these transformation matrices maps the new parameters/variables
  29  * back to the original ones.
  30  */
  31 int remove_all_equalities(Polyhedron **P, Polyhedron **C, Matrix **CPP, Matrix **CVP,
  32                           unsigned nparam, unsigned MaxRays)
  33 {
  34     Matrix *CV = NULL;
  35     Matrix *CP = NULL;
  36     Polyhedron *Q = *P;
  37     Polyhedron *D = NULL;
  38     Polyhedron *R;
  39     int i;
  40     Matrix M;
  41
  42     if (Q->NbEq == 0)
  43         return 0;
  44
  45     if (C) {
  46         D = *C;
  47         assert(D->Dimension == nparam);
  48     }
  49
  50     /* compress_parms doesn't like equalities that only involve parameters */
  51     for (i = 0; i < Q->NbEq; ++i)
  52         if (First_Non_Zero(Q->Constraint[i]+1, Q->Dimension-nparam) == -1)
  53             break;
  54
  55     if (i < Q->NbEq) {
  56         Matrix *M = Matrix_Alloc(Q->NbEq, 1+nparam+1);
  57         int n = 0;
  58         for (; i < Q->NbEq; ++i) {
  59             if (First_Non_Zero(Q->Constraint[i]+1, Q->Dimension-nparam) == -1)
  60                 Vector_Copy(Q->Constraint[i]+1+Q->Dimension-nparam,
  61                             M->p[n++]+1, nparam+1);
  62         }
  63         M->NbRows = n;
  64         CV = compress_variables(M, 0);
  65         Matrix_Free(M);
  66         transform(&Q, &D, CV, Q != *P, MaxRays);
  67         nparam = CV->NbColumns-1;
  68     }
  69
  70     if (emptyQ2(Q)) {
  71         if (CV)
  72             Matrix_Free(CV);
  73         *P = Q;
  74         if (C)
  75             *C = D;
  76         return 0;
  77     }
  78
  79     if (Q->NbEq == 0) {
  80         CP = CV;
  81         CV = NULL;
  82     } else {
  83         /* Matrix "view" of equalities */
  84         M.NbRows = Q->NbEq;
  85         M.NbColumns = Q->Dimension+2;
  86         M.p_Init = Q->p_Init;
  87         M.p = Q->Constraint;
  88         CP = compress_parms(&M, nparam);
  89
  90         if (isIdentity(CP)) {
  91             Matrix_Free(CP);
  92             CP = CV;
  93             CV = NULL;
  94         } else {
  95             transform(&Q, &D, CP, Q != *P, MaxRays);
  96             if (CV) {
  97                 Matrix *T = CP;
  98                 CP = Matrix_Alloc(CV->NbRows, T->NbColumns);
  99                 Matrix_Product(CV, T, CP);
 100                 Matrix_Free(T);
 101                 Matrix_Free(CV);
 102                 CV = NULL;
 103             }
 104             nparam = CP->NbColumns-1;
 105         }
 106
 107         /* Matrix "view" of equalities */
 108         M.NbRows = Q->NbEq;
 109         M.NbColumns = Q->Dimension+2;
 110         M.p_Init = Q->p_Init;
 111         M.p = Q->Constraint;
 112         CV = compress_variables(&M, nparam);
 113         if (!CV) {
 114             if (Q != *P)
 115                 Polyhedron_Free(Q);
 116             Q = NULL;
 117         } else if (isIdentity(CV)) {
 118             Matrix_Free(CV);
 119             CV = NULL;
 120         } else {
 121             R = Polyhedron_Preimage(Q, CV, MaxRays);
 122             if (Q != *P)
 123                 Polyhedron_Free(Q);
 124             Q = R;
 125         }
 126     }
 127
 128     *P = Q;
 129     if (C)
 130         *C = D;
 131     if (CPP)
 132         *CPP = CP;
 133     else if (CP)
 134         Matrix_Free(CP);
 135     if (CVP)
 136         *CVP = CV;
 137     else if (CV)
 138         Matrix_Free(CV);
 139     return 1;
 140 }