remove_equalities.c

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