remove_equalities.c

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