topcom.c

   1 #include <barvinok/util.h>
   2 #include <barvinok/options.h>
   3 #include <unistd.h>
   4 #include "topcom.h"
   5 #include "config.h"
   6
   7 #define ALLOC(type) (type*)malloc(sizeof(type))
   8
   9 void run_points2triangs(pid_t *child, int *in, int *out)
  10 {
  11     int in_fd[2], out_fd[2];
  12
  13     if (pipe(in_fd))
  14         assert(0);
  15     if (pipe(out_fd))
  16         assert(0);
  17     *child = fork();
  18     assert(*child >= 0);
  19     if (!*child) {
  20         int rc;
  21
  22         dup2(in_fd[0], 0);
  23         dup2(out_fd[1], 1);
  24         close(in_fd[0]);
  25         close(out_fd[1]);
  26         close(in_fd[1]);
  27         close(out_fd[0]);
  28
  29         rc = execl(POINTS2TRIANGS_PATH, "points2triangs", "--regular", NULL);
  30         assert(0);
  31     }
  32     close(in_fd[0]);
  33     close(out_fd[1]);
  34     *in = in_fd[1];
  35     *out = out_fd[0];
  36 }
  37
  38 struct domain {
  39     Param_Domain    domain;
  40     int             F_len;
  41 };
  42
  43 static Param_Vertices *construct_vertex(unsigned *vertex_facets,
  44                                        Matrix *Constraints,
  45                                        int d, unsigned nparam, unsigned MaxRays)
  46 {
  47     unsigned nvar = Constraints->NbColumns-2 - nparam;
  48     Matrix *A = Matrix_Alloc(nvar+1, nvar+1);
  49     Matrix *inv = Matrix_Alloc(nvar+1, nvar+1);
  50     Matrix *B = Matrix_Alloc(nvar, nparam+2);
  51     Matrix *V = Matrix_Alloc(nvar, nparam+2);
  52     Matrix *Domain = Matrix_Alloc(d-nvar, nparam+2);
  53     Polyhedron *AD;
  54     unsigned bx;
  55     int i, j, ix;
  56     int ok;
  57     Param_Vertices *vertex;
  58
  59     for (j = 0, i = 0, ix = 0, bx = MSB; i < d; ++i) {
  60         if ((vertex_facets[ix] & bx) == bx) {
  61             Vector_Copy(Constraints->p[i]+1, A->p[j], nvar);
  62             Vector_Oppose(Constraints->p[i]+1+nvar, B->p[j++], nparam+1);
  63         }
  64         NEXT(ix, bx);
  65     }
  66     assert(j == nvar);
  67     value_set_si(A->p[nvar][nvar], 1);
  68     ok = Matrix_Inverse(A, inv);
  69     assert(ok);
  70     Matrix_Free(A);
  71     inv->NbRows = nvar;
  72     inv->NbColumns = nvar;
  73     Matrix_Product(inv, B, V);
  74     Matrix_Free(B);
  75     for (i = 0; i < nvar; ++i) {
  76         value_assign(V->p[i][nparam+1], inv->p[nvar][nvar]);
  77         Vector_Normalize(V->p[i], V->NbColumns);
  78     }
  79     Matrix_Free(inv);
  80     for (j = 0, i = 0, ix = 0, bx = MSB; i < d; ++i) {
  81         if ((vertex_facets[ix] & bx) == bx) {
  82             NEXT(ix, bx);
  83             continue;
  84         }
  85         Param_Inner_Product(Constraints->p[i], V, Domain->p[j]);
  86         if (First_Non_Zero(Domain->p[j]+1, nparam+1) == -1)
  87             vertex_facets[ix] |= bx;
  88         else
  89             value_set_si(Domain->p[j++][0], 1);
  90         NEXT(ix, bx);
  91     }
  92     Domain->NbRows = j;
  93     A = Matrix_Copy(Domain);
  94     AD = Constraints2Polyhedron(A, MaxRays);
  95     Matrix_Free(A);
  96     POL_ENSURE_VERTICES(AD);
  97     /* A vertex with a lower-dimensional activity domain
  98      * saturates more facets than derived above and is actually
  99      * the superimposition of two or more vertices.
 100      * We therefore discard the domain and (ultimately)
 101      * the chamber containing it.
 102      * We keep the vertex, though, since it may reappear
 103      * in other chambers, which will then likewise be discarded.
 104      * The same holds if the activity domain is empty.
 105      */
 106     if (AD->NbEq > 0) {
 107         Matrix_Free(Domain);
 108         Domain = NULL;
 109     }
 110     Polyhedron_Free(AD);
 111     vertex = calloc(1, sizeof(Param_Vertices));
 112     vertex->Facets = vertex_facets;
 113     vertex->Vertex = V;
 114     vertex->Domain = Domain;
 115     return vertex;
 116 }
 117
 118 static int add_vertex_to_domain(Param_Vertices **vertices, int words,
 119                                 unsigned *vertex_facets,
 120                                 Matrix *Constraints, int d, unsigned nparam,
 121                                 struct domain *domain,
 122                                 unsigned MaxRays)
 123 {
 124     Param_Vertices *vertex;
 125     unsigned vbx;
 126     int vi, vix;
 127
 128     for (vi = 0, vix = 0, vbx = MSB;
 129             *vertices;
 130             vertices = &(*vertices)->next, ++vi) {
 131         int i;
 132         for (i = 0; i < words; ++i)
 133             if (((*vertices)->Facets[i] & vertex_facets[i]) != vertex_facets[i])
 134                 break;
 135         if (i >= words) {
 136             if (!(*vertices)->Domain)
 137                 domain->F_len = 0;
 138             else
 139                 domain->domain.F[vix] |= vbx;
 140             free(vertex_facets);
 141             return 0;
 142         }
 143         NEXT(vix, vbx);
 144     }
 145     if (domain->F_len <= vix) {
 146         domain->F_len++;
 147         domain->domain.F = realloc(domain->domain.F,
 148                                    domain->F_len * sizeof(unsigned));
 149         domain->domain.F[domain->F_len-1] = 0;
 150     }
 151     vertex = construct_vertex(vertex_facets, Constraints, d, nparam, MaxRays);
 152     if (!vertex->Domain)
 153         domain->F_len = 0;
 154     else
 155         domain->domain.F[vix] |= vbx;
 156     vertex->next = *vertices;
 157     *vertices = vertex;
 158     return 1;
 159 }
 160
 161 static void compute_domain(struct domain *domain, Param_Vertices *vertices,
 162                            Polyhedron *C, unsigned MaxRays)
 163 {
 164     unsigned bx;
 165     int i, ix, j;
 166     int nbV = bit_vector_count(domain->domain.F, domain->F_len);
 167     unsigned cols = vertices->Domain->NbColumns;
 168     unsigned rows = vertices->Domain->NbRows;
 169     Matrix *Constraints = Matrix_Alloc(nbV * rows + C->NbConstraints, cols);
 170
 171     for (i = 0, j = 0, ix = 0, bx = MSB;
 172             vertices;
 173             vertices = vertices->next, ++i) {
 174         if ((domain->domain.F[ix] & bx) == bx)
 175             Vector_Copy(vertices->Domain->p[0],
 176                         Constraints->p[(j++)*rows], rows * cols);
 177         NEXT(ix, bx);
 178     }
 179     Vector_Copy(C->Constraint[0], Constraints->p[j*rows], C->NbConstraints * cols);
 180     domain->domain.Domain = Constraints2Polyhedron(Constraints, MaxRays);
 181     Matrix_Free(Constraints);
 182 }
 183
 184 static void add_domain(struct domain **domains, struct domain *domain,
 185                        Param_Vertices *vertices, Polyhedron *C,
 186                        struct barvinok_options *options)
 187 {
 188     options->stats->topcom_chambers++;
 189
 190     for (; *domains; domains = (struct domain **)&(*domains)->domain.next) {
 191         int i;
 192         for (i = 0; i < (*domains)->F_len; ++i)
 193             if (((*domains)->domain.F[i] & domain->domain.F[i])
 194                         != domain->domain.F[i])
 195                 break;
 196         if (i < (*domains)->F_len)
 197             continue;
 198         for (; i < domain->F_len; ++i)
 199             if (domain->domain.F[i])
 200                 break;
 201         if (i >= domain->F_len) {
 202             Param_Domain_Free(&domain->domain);
 203             return;
 204         }
 205     }
 206     options->stats->topcom_distinct_chambers++;
 207     compute_domain(domain, vertices, C, options->MaxRays);
 208     *domains = domain;
 209 }
 210
 211 #define INT_BITS (sizeof(unsigned) * 8)
 212
 213 /* Remove any empty or lower-dimensional chamber.  The latter
 214  * lie on the boundary of the context and are facets of other chambers.
 215  *
 216  * While we are examining each chamber, also extend the F vector
 217  * of each chamber to the maximum.
 218  */
 219 static void remove_empty_chambers(Param_Domain **PD, unsigned vertex_words)
 220 {
 221     while (*PD) {
 222         int remove = 0;
 223         int i;
 224
 225         if ((*PD)->Domain->NbEq > 0)
 226             remove = 1;
 227         else {
 228             POL_ENSURE_FACETS((*PD)->Domain);
 229             if ((*PD)->Domain->NbEq > 0)
 230                 remove = 1;
 231         }
 232         if (remove) {
 233             Param_Domain *D = *PD;
 234             *PD = (*PD)->next;
 235             D->next = NULL;
 236             Param_Domain_Free(D);
 237             continue;
 238         }
 239         if ((i = ((struct domain*)(*PD))->F_len) < vertex_words)
 240             (*PD)->F = realloc((*PD)->F, vertex_words * sizeof(unsigned));
 241         for (; i < vertex_words; ++i)
 242             (*PD)->F[i] = 0;
 243         PD = &(*PD)->next;
 244     }
 245 }
 246
 247 static Param_Polyhedron *points2triangs(Matrix *K, Polyhedron *P, Polyhedron *C,
 248                                         struct barvinok_options *options)
 249 {
 250     int in, out;
 251     int i, j;
 252     pid_t child;
 253     FILE *fin, *fout;
 254     int d = K->NbRows;
 255     int words = (d+INT_BITS-1)/INT_BITS;
 256     Param_Vertices *vertices = NULL;
 257     struct domain *domains = NULL;
 258     int vertex_words = 1;
 259     Param_Polyhedron *PP = ALLOC(Param_Polyhedron);
 260     unsigned MaxRays = options->MaxRays;
 261
 262     PP->nbV = 0;
 263     PP->Constraints = Polyhedron2Constraints(P);
 264     /* We need the exact facets, because we may make some of them open later */
 265     POL_UNSET(options->MaxRays, POL_INTEGER);
 266
 267     run_points2triangs(&child, &in, &out);
 268
 269     fin = fdopen(in, "w");
 270     fprintf(fin, "[\n");
 271     for (i = 0; i < K->NbRows; ++i) {
 272         fprintf(fin, "[");
 273         for (j = 0; j < K->NbColumns; ++j)
 274             value_print(fin, P_VALUE_FMT, K->p[i][j]);
 275         fprintf(fin, "]");
 276     }
 277     fprintf(fin, "]\n");
 278     fclose(fin);
 279
 280     fout = fdopen(out, "r");
 281     while (fscanf(fout, "T[%d]:={", &i) == 1) {
 282         struct domain *domain = ALLOC(struct domain);
 283         memset(domain, 0, sizeof(*domain));
 284         domain->domain.F = calloc(vertex_words, sizeof(unsigned));
 285         domain->F_len = vertex_words;
 286
 287         while (fgetc(fout) == '{') {    /* '{' or closing '}' */
 288             int c;
 289             unsigned *F = calloc(words, sizeof(unsigned));
 290
 291             for (j = 0; j < K->NbColumns; ++j) {
 292                 unsigned v, shift;
 293                 fscanf(fout, "%d", &v);
 294                 shift = INT_BITS - (v % INT_BITS) - 1;
 295                 F[v / INT_BITS] |= 1u << shift;
 296                 fgetc(fout); /* , or } */
 297             }
 298             if (!domain->F_len)
 299                 free(F);
 300             else if (add_vertex_to_domain(&vertices, words, F, PP->Constraints,
 301                                           d, C->Dimension,
 302                                           domain, options->MaxRays))
 303                 ++PP->nbV;
 304             if ((c = fgetc(fout)) != ',')       /* , or } */
 305                 ungetc(c, fout);
 306         }
 307         if (domain->F_len)
 308             vertex_words = domain->F_len;
 309         fgetc(fout); /* ; */
 310         fgetc(fout); /* \n */
 311         if (bit_vector_count(domain->domain.F, domain->F_len) > 0)
 312             add_domain(&domains, domain, vertices, C, options);
 313         else {
 314             options->stats->topcom_empty_chambers++;
 315             Param_Domain_Free(&domain->domain);
 316         }
 317     }
 318     fclose(fout);
 319
 320     PP->V = vertices;
 321     PP->D = &domains->domain;
 322
 323     remove_empty_chambers(&PP->D, vertex_words);
 324
 325     options->MaxRays = MaxRays;
 326
 327     return PP;
 328 }
 329
 330 /* Assumes M is of full row rank */
 331 static Matrix *null_space(Matrix *M)
 332 {
 333     Matrix *H, *Q, *U;
 334     Matrix *N;
 335     int i;
 336
 337     left_hermite(M, &H, &Q, &U);
 338     N = Matrix_Alloc(M->NbColumns, M->NbColumns - M->NbRows);
 339     for (i = 0; i < N->NbRows; ++i)
 340         Vector_Copy(U->p[i] + M->NbRows, N->p[i], N->NbColumns);
 341     Matrix_Free(H);
 342     Matrix_Free(Q);
 343     Matrix_Free(U);
 344     return N;
 345 }
 346
 347 static void SwapColumns(Value **V, int n, int i, int j)
 348 {
 349     int r;
 350
 351     for (r = 0; r < n; ++r)
 352         value_swap(V[r][i], V[r][j]);
 353 }
 354
 355 static int is_unit_row(Value *row, int pos, int len)
 356 {
 357     if (!value_one_p(row[pos]) && !value_mone_p(row[pos]))
 358         return 0;
 359     return First_Non_Zero(row+pos+1, len-(pos+1)) == -1;
 360 }
 361
 362 /* Transform the constraints of P to "standard form".
 363  * In particular, if P is described by
 364  *              A x + b(p) >= 0
 365  * then this function returns a matrix H = A U, A = H Q, such
 366  * that D x' = D Q x >= -b(p), with D a diagonal matrix with
 367  * positive entries.  The calling function can then construct
 368  * the standard form H' x' - I s + b(p) = 0, with H' the rows of H
 369  * that are not positive multiples of unit vectors
 370  * (since those correspond to D x' >= -b(p)).
 371  * The number of rows in H' is returned in *rows_p.
 372  * Optionally returns the matrix that maps the new variables
 373  * back to the old variables x = U x'.
 374  * Note that the rows of H (and P) may be reordered by this function.
 375  */
 376 Matrix *standard_constraints(Polyhedron *P, unsigned nparam, int *rows_p,
 377                              Matrix **T)
 378 {
 379     unsigned nvar = P->Dimension - nparam;
 380     int i, j, d;
 381     int rows;
 382     Matrix *M;
 383     Matrix *H, *U, *Q;
 384
 385     assert(P->NbEq == 0);
 386
 387     /* move constraints only involving parameters down
 388      * and move unit vectors (if there are any) to the right place.
 389      */
 390     for (d = 0, j = P->NbConstraints; d < j; ++d) {
 391         int index;
 392         index = First_Non_Zero(P->Constraint[d]+1, nvar);
 393         if (index != -1) {
 394             if (index != d &&
 395                 is_unit_row(P->Constraint[d]+1, index, nvar)) {
 396                 Vector_Exchange(P->Constraint[d], P->Constraint[index],
 397                                 P->Dimension+2);
 398                 if (index > d)
 399                     --d;
 400             }
 401             continue;
 402         }
 403         while (d < --j && First_Non_Zero(P->Constraint[j]+1, nvar) == -1)
 404             ;
 405         if (d >= j)
 406             break;
 407         Vector_Exchange(P->Constraint[d], P->Constraint[j], P->Dimension+2);
 408     }
 409     M = Matrix_Alloc(d, nvar);
 410     for (j = 0; j < d; ++j)
 411         Vector_Copy(P->Constraint[j]+1, M->p[j], nvar);
 412
 413     neg_left_hermite(M, &H, &Q, &U);
 414     Matrix_Free(M);
 415     Matrix_Free(Q);
 416     if (T)
 417         *T = U;
 418     else
 419         Matrix_Free(U);
 420
 421     /* Rearrange rows such that top of H is lower diagonal and
 422      * compute the number of non (multiple of) unit-vector rows.
 423      */
 424     rows = H->NbRows-nvar;
 425     for (i = 0; i < H->NbColumns; ++i) {
 426         for (j = i; j < H->NbRows; ++j)
 427             if (value_notzero_p(H->p[j][i]))
 428                 break;
 429         if (j != i) {
 430             Vector_Exchange(P->Constraint[i], P->Constraint[j], P->Dimension+2);
 431             Vector_Exchange(H->p[i], H->p[j], H->NbColumns);
 432         }
 433         if (First_Non_Zero(H->p[i], i) != -1)
 434             rows++;
 435     }
 436     *rows_p = rows;
 437
 438     return H;
 439 }
 440
 441 /* C is assumed to be the "true" context, i.e., it has been intersected
 442  * with the projection of P onto the parameter space.
 443  * Furthermore, P and C are assumed to be full-dimensional.
 444  */
 445 Param_Polyhedron *TC_P2PP(Polyhedron *P, Polyhedron *C,
 446                           struct barvinok_options *options)
 447 {
 448     unsigned nparam = C->Dimension;
 449     unsigned nvar = P->Dimension - C->Dimension;
 450     int i, j;
 451     Matrix *H;
 452     Matrix *A;
 453     Matrix *K;
 454     int rows;
 455     Param_Polyhedron *PP;
 456
 457     assert(C->NbEq == 0);
 458
 459     H = standard_constraints(P, nparam, &rows, NULL);
 460
 461     A = Matrix_Alloc(rows, nvar+rows);
 462     for (i = nvar; i < H->NbRows; ++i) {
 463         Vector_Oppose(H->p[i], A->p[i-nvar], H->NbColumns);
 464         value_set_si(A->p[i-nvar][i], 1);
 465     }
 466     for (i = 0, j = H->NbRows-nvar; i < nvar; ++i) {
 467         if (First_Non_Zero(H->p[i], i) == -1)
 468             continue;
 469         Vector_Oppose(H->p[i], A->p[j], H->NbColumns);
 470         value_set_si(A->p[j][j+nvar], 1);
 471         SwapColumns(A->p, A->NbRows, j+nvar, i);
 472         ++j;
 473     }
 474     K = null_space(A);
 475     Matrix_Free(A);
 476     /* Ignore extra constraints */
 477     K->NbRows = H->NbRows;
 478     Matrix_Free(H);
 479     PP = points2triangs(K, P, C, options);
 480     Matrix_Free(K);
 481     return PP;
 482 }