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