hilbert.c

   1 #define Vector ZSolveVector
   2 #define Matrix ZSolveMatrix
   3 #include "zsolve/libzsolve.h"
   4 #undef Vector
   5 #undef Matrix
   6 #include "hilbert.h"
   7
   8 static ZSolveMatrix Matrix2zsolve(Matrix *M)
   9 {
  10     int i, j;
  11     ZSolveMatrix zmatrix;
  12
  13     zmatrix = createMatrix(M->NbColumns-2, M->NbRows);
  14     for (i = 0; i < M->NbRows; ++i)
  15         for (j = 0; j < M->NbColumns-2; ++j) {
  16             assert(mpz_cmp_si(M->p[i][1+j], -MAXINT) > 0);
  17             assert(mpz_cmp_si(M->p[i][1+j], MAXINT) < 0);
  18             zmatrix->Data[i*zmatrix->Width+j] = mpz_get_si(M->p[i][1+j]);
  19         }
  20
  21     return zmatrix;
  22 }
  23
  24 static Matrix *VectorArray2Matrix(VectorArray array, unsigned cols)
  25 {
  26     int i, j;
  27     Matrix *M = Matrix_Alloc(array->Size, cols+1);
  28
  29     for (i = 0; i < array->Size; ++i) {
  30         for (j = 0; j < cols; ++j)
  31             value_set_si(M->p[i][j], array->Data[i][j]);
  32         value_set_si(M->p[i][cols], 1);
  33     }
  34     return M;
  35 }
  36
  37 static void Polyhedron_Remove_Positivity_Constraint(Polyhedron *P)
  38 {
  39     int i;
  40
  41     for (i = 0; i < P->NbConstraints; ++i) {
  42         if (First_Non_Zero(P->Constraint[i]+1, P->Dimension) != -1)
  43             continue;
  44         if (i < P->NbConstraints-1)
  45             Vector_Exchange(P->Constraint[i],
  46                             P->Constraint[P->NbConstraints-1],
  47                             P->Dimension+2);
  48         P->NbConstraints--;
  49         --i;
  50     }
  51 }
  52
  53 static Matrix *Polyhedron2standard_form(Polyhedron *P, Matrix **T)
  54 {
  55     int i;
  56     unsigned dim = P->Dimension;
  57     Matrix *M2 = Matrix_Alloc(P->NbConstraints+1, dim+1);
  58     Matrix *H, *Q;
  59
  60     Polyhedron_Remove_Positivity_Constraint(P);
  61     for (i = 0; i < P->NbConstraints; ++i) {
  62         assert(value_zero_p(P->Constraint[i][1+dim]));
  63         Vector_Copy(P->Constraint[i]+1, M2->p[i], dim);
  64     }
  65     value_set_si(M2->p[P->NbConstraints][dim], 1);
  66     neg_left_hermite(M2, &H, &Q, T);
  67     Matrix_Free(Q);
  68     Matrix_Free(M2);
  69
  70     M2 = Matrix_Alloc(P->NbConstraints, 2+dim+(P->NbConstraints-P->NbEq));
  71
  72     for (i = 0; i < P->NbEq; ++i)
  73         Vector_Copy(H->p[i], M2->p[i]+1, dim);
  74     for (i = P->NbEq; i < P->NbConstraints; ++i) {
  75         Vector_Copy(H->p[i], M2->p[i]+1, dim);
  76         value_set_si(M2->p[i][1+dim+i-P->NbEq], -1);
  77     }
  78     Matrix_Free(H);
  79     return M2;
  80 }
  81
  82 /* Assumes C is a linear cone (i.e. with apex zero).
  83  * All equalities are removed first to speed up the computation
  84  * in zsolve.
  85  */
  86 Matrix *Cone_Hilbert_Basis(Polyhedron *C, unsigned MaxRays)
  87 {
  88     unsigned dim;
  89     int i;
  90     Matrix *M2, *M3, *T;
  91     Matrix *CV = NULL;
  92     LinearSystem initialsystem;
  93     ZSolveMatrix matrix;
  94     ZSolveVector rhs;
  95     ZSolveContext ctx;
  96
  97     remove_all_equalities(&C, NULL, NULL, &CV, 0, MaxRays);
  98     dim = C->Dimension;
  99
 100     for (i = 0; i < C->NbConstraints; ++i)
 101         assert(value_zero_p(C->Constraint[i][1+dim]) ||
 102                First_Non_Zero(C->Constraint[i]+1, dim) == -1);
 103
 104     M2 = Polyhedron2standard_form(C, &T);
 105     matrix = Matrix2zsolve(M2);
 106     Matrix_Free(M2);
 107
 108     rhs = createVector(matrix->Height);
 109     for (i = 0; i < matrix->Height; i++)
 110         rhs[i] = 0;
 111
 112     initialsystem = createLinearSystem();
 113     setLinearSystemMatrix(initialsystem, matrix);
 114     deleteMatrix(matrix);
 115
 116     setLinearSystemRHS(initialsystem, rhs);
 117     deleteVector(rhs);
 118
 119     setLinearSystemLimit(initialsystem, -1, 0, MAXINT, 0);
 120     setLinearSystemEquationType(initialsystem, -1, EQUATION_EQUAL, 0);
 121
 122     ctx = createZSolveContextFromSystem(initialsystem, NULL, 0, 0, NULL, NULL);
 123     zsolveSystem(ctx, 0);
 124
 125     M2 = VectorArray2Matrix(ctx->Homs, C->Dimension);
 126     deleteZSolveContext(ctx, 1);
 127     Matrix_Transposition(T);
 128     M3 = Matrix_Alloc(M2->NbRows, M2->NbColumns);
 129     Matrix_Product(M2, T, M3);
 130     Matrix_Free(M2);
 131     Matrix_Free(T);
 132
 133     if (CV) {
 134         Matrix *T, *M;
 135         T = Transpose(CV);
 136         M = Matrix_Alloc(M3->NbRows, T->NbColumns);
 137         Matrix_Product(M3, T, M);
 138         Matrix_Free(M3);
 139         Matrix_Free(CV);
 140         Matrix_Free(T);
 141         Polyhedron_Free(C);
 142         M3 = M;
 143     }
 144
 145     return M3;
 146 }