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