source/kernel/alpha.c

   1 /*
   2     This file is part of PolyLib.
   3
   4     PolyLib is free software: you can redistribute it and/or modify
   5     it under the terms of the GNU General Public License as published by
   6     the Free Software Foundation, either version 3 of the License, or
   7     (at your option) any later version.
   8
   9     PolyLib is distributed in the hope that it will be useful,
  10     but WITHOUT ANY WARRANTY; without even the implied warranty of
  11     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12     GNU General Public License for more details.
  13
  14     You should have received a copy of the GNU General Public License
  15     along with PolyLib.  If not, see <http://www.gnu.org/licenses/>.
  16 */
  17
  18 /* polytest.c */
  19 #include <stdio.h>
  20 #include <polylib/polylib.h>
  21
  22 /* alpha.c
  23      COPYRIGHT
  24           Both this software and its documentation are
  25
  26               Copyright 1993 by IRISA /Universite de Rennes I - France,
  27               Copyright 1995,1996 by BYU, Provo, Utah
  28                          all rights reserved.
  29
  30 */
  31
  32 #include <stdio.h>
  33 #include <stdlib.h>
  34 #include <string.h>
  35
  36 /*---------------------------------------------------------------------*/
  37 /* int exist_points(pos,P,context)                                     */
  38 /*    pos : index position of current loop index (0..hdim-1)           */
  39 /*    P: loop domain                                                   */
  40 /*    context : context values for fixed indices                       */
  41 /* recursive procedure, recurs for each imbriquation                   */
  42 /* returns 1 if there exists any integer points in this polyhedron     */
  43 /* returns 0 if there are none                                         */
  44 /*---------------------------------------------------------------------*/
  45 static int exist_points(int pos,Polyhedron *Pol,Value *context) {
  46
  47   Value LB, UB, k,tmp;
  48
  49   value_init(LB); value_init(UB);
  50   value_init(k);  value_init(tmp);
  51   value_set_si(LB,0);
  52   value_set_si(UB,0);
  53
  54   /* Problem if UB or LB is INFINITY */
  55   if (lower_upper_bounds(pos,Pol,context,&LB,&UB) !=0) {
  56     errormsg1("exist_points", "infdom", "infinite domain");
  57     value_clear(LB);
  58     value_clear(UB);
  59     value_clear(k);
  60     value_clear(tmp);
  61     return -1;
  62   }
  63   value_set_si(context[pos],0);
  64   if(value_lt(UB,LB)) {
  65     value_clear(LB);
  66     value_clear(UB);
  67     value_clear(k);
  68     value_clear(tmp);
  69     return 0;
  70   }
  71   if (!Pol->next) {
  72     value_subtract(tmp,UB,LB);
  73     value_increment(tmp,tmp);
  74     value_clear(UB);
  75     value_clear(LB);
  76     value_clear(k);
  77     return (value_pos_p(tmp));
  78   }
  79
  80   for (value_assign(k,LB);value_le(k,UB);value_increment(k,k)) {
  81
  82     /* insert k in context */
  83     value_assign(context[pos],k);
  84     if (exist_points(pos+1,Pol->next,context) > 0 ) {
  85       value_clear(LB); value_clear(UB);
  86       value_clear(k); value_clear(tmp);
  87       return 1;
  88     }
  89   }
  90   /* Reset context */
  91   value_set_si(context[pos],0);
  92   value_clear(UB); value_clear(LB);
  93   value_clear(k); value_clear(tmp);
  94   return 0;
  95 }
  96
  97 /*--------------------------------------------------------------*/
  98 /* Test to see if there are any integral points in a polyhedron */
  99 /*--------------------------------------------------------------*/
 100 int Polyhedron_Not_Empty(Polyhedron *P,Polyhedron *C,int MAXRAYS) {
 101
 102   int res,i;
 103   Value *context;
 104   Polyhedron *L;
 105
 106   POL_ENSURE_FACETS(P);
 107   POL_ENSURE_VERTICES(P);
 108   POL_ENSURE_FACETS(C);
 109   POL_ENSURE_VERTICES(C);
 110
 111   /* Create a context vector size dim+2 and set it to all zeros */
 112   context = (Value *) malloc((P->Dimension+2)*sizeof(Value));
 113
 114   /* Initialize array 'context' */
 115   for (i=0;i<(P->Dimension+2);i++)
 116     value_init(context[i]);
 117
 118   Vector_Set(context,0,(P->Dimension+2));
 119
 120   /* Set context[P->Dimension+1] = 1  (the constant) */
 121   value_set_si(context[P->Dimension+1],1);
 122
 123   L = Polyhedron_Scan(P,C,MAXRAYS);
 124   res = exist_points(1,L,context);
 125   Domain_Free(L);
 126
 127   /* Clear array 'context' */
 128   for (i=0;i<(P->Dimension+2);i++)
 129     value_clear(context[i]);
 130   free(context);
 131   return res;
 132 }
 133
 134 /* PolyhedronLTQ ( P1, P2 ) */
 135 /* P1 and P2 must be simple polyhedra */
 136 /* result =  0 --> not comparable */
 137 /* result = -1 --> P1 < P2        */
 138 /* result =  1 --> P1 > P2        */
 139 /* INDEX  = 1 .... Dimension      */
 140 int PolyhedronLTQ (Polyhedron *Pol1,Polyhedron *Pol2,int INDEX, int PDIM, int NbMaxConstrs) {
 141
 142   int res, dim, i, j, k;
 143   Polyhedron *Q1, *Q2, *Q3, *Q4, *Q;
 144   Matrix *Mat;
 145
 146   if (Pol1->next || Pol2->next) {
 147     errormsg1("PolyhedronLTQ", "compoly", "Can only compare polyhedra");
 148     return 0;
 149   }
 150   if (Pol1->Dimension != Pol2->Dimension) {
 151     errormsg1("PolyhedronLTQ","diffdim","Polyhedra are not same dimension");
 152     return 0;
 153   }
 154   dim = Pol1->Dimension+2;
 155
 156   POL_ENSURE_FACETS(Pol1);
 157   POL_ENSURE_VERTICES(Pol1);
 158   POL_ENSURE_FACETS(Pol2);
 159   POL_ENSURE_VERTICES(Pol2);
 160
 161 #ifdef DEBUG
 162   fprintf(stdout, "P1\n");
 163   Polyhedron_Print(stdout,P_VALUE_FMT,Pol1);
 164   fprintf(stdout, "P2\n");
 165   Polyhedron_Print(stdout,P_VALUE_FMT,Pol2);
 166 #endif
 167
 168   /* Create the Line to add */
 169   k = Pol1->Dimension-INDEX+1-PDIM;
 170   Mat = Matrix_Alloc(k,dim);
 171   Vector_Set(Mat->p_Init,0,dim*k);
 172   for(j=0,i=INDEX;j<k;i++,j++)
 173     value_set_si(Mat->p[j][i],1);
 174
 175   Q1 = AddRays(Mat->p[0],k,Pol1,NbMaxConstrs);
 176   Q2 = AddRays(Mat->p[0],k,Pol2,NbMaxConstrs);
 177
 178 #ifdef DEBUG
 179   fprintf(stdout, "Q1\n");
 180   Polyhedron_Print(stdout,P_VALUE_FMT,Q1);
 181   fprintf(stdout, "Q2\n");
 182   Polyhedron_Print(stdout,P_VALUE_FMT,Q2);
 183 #endif
 184
 185   Matrix_Free(Mat);
 186   Q  = DomainIntersection(Q1,Q2,NbMaxConstrs);
 187
 188 #ifdef DEBUG
 189   fprintf(stdout, "Q\n");
 190   Polyhedron_Print(stdout,P_VALUE_FMT,Q);
 191 #endif
 192
 193   Domain_Free(Q1);
 194   Domain_Free(Q2);
 195
 196   if (emptyQ(Q)) res = 0;       /* not comparable */
 197   else {
 198     Q1 = DomainIntersection(Pol1,Q,NbMaxConstrs);
 199     Q2 = DomainIntersection(Pol2,Q,NbMaxConstrs);
 200
 201 #ifdef DEBUG
 202     fprintf(stdout, "Q1\n");
 203     Polyhedron_Print(stdout,P_VALUE_FMT,Q1);
 204     fprintf(stdout, "Q2\n");
 205     Polyhedron_Print(stdout,P_VALUE_FMT,Q2);
 206 #endif
 207
 208     k = Q1->NbConstraints + Q2->NbConstraints;
 209     Mat = Matrix_Alloc(k, dim);
 210     Vector_Set(Mat->p_Init,0,k*dim);
 211
 212     /* First compute surrounding polyhedron */
 213     j=0;
 214     for (i=0; i<Q1->NbConstraints; i++) {
 215       if ((value_one_p(Q1->Constraint[i][0])) && (value_pos_p(Q1->Constraint[i][INDEX]))) {
 216
 217         /* keep Q1's lower bounds */
 218         for (k=0; k<dim; k++)
 219           value_assign(Mat->p[j][k],Q1->Constraint[i][k]);
 220         j++;
 221       }
 222     }
 223     for (i=0; i<Q2->NbConstraints; i++) {
 224       if ((value_one_p(Q2->Constraint[i][0])) && (value_neg_p(Q2->Constraint[i][INDEX]))) {
 225
 226         /* and keep Q2's upper bounds */
 227         for (k=0; k<dim; k++)
 228           value_assign(Mat->p[j][k],Q2->Constraint[i][k]);
 229         j++;
 230       }
 231     }
 232     Q4 = AddConstraints(Mat->p[0], j, Q, NbMaxConstrs);
 233     Matrix_Free(Mat);
 234
 235 #ifdef debug
 236     fprintf(stderr, "Q4 surrounding polyhedron\n");
 237     Polyhderon_Print(stderr,P_VALUE_FMT, Q4);
 238 #endif
 239
 240     /* if surrounding polyhedron is empty, D1>D2 */
 241     if (emptyQ(Q4)) {
 242       res = 1;
 243
 244 #ifdef debug
 245       fprintf(stderr, "Surrounding polyhedron is empty\n");
 246 #endif
 247       goto LTQdone2;
 248     }
 249
 250     /* Test if Q1 < Q2 */
 251     /* Build a constraint array for >= Q1 and <= Q2 */
 252     Mat = Matrix_Alloc(2,dim);
 253     Vector_Set(Mat->p_Init,0,2*dim);
 254
 255     /* Choose a contraint from Q1 */
 256     for (i=0; i<Q1->NbConstraints; i++) {
 257       if (value_zero_p(Q1->Constraint[i][0])) {
 258
 259         /* Equality */
 260         if (value_zero_p(Q1->Constraint[i][INDEX])) {
 261
 262           /* Ignore side constraint (they are in Q) */
 263           continue;
 264         }
 265         else if (value_neg_p(Q1->Constraint[i][INDEX])) {
 266
 267           /* copy -constraint to Mat */
 268           value_set_si(Mat->p[0][0],1);
 269           for (k=1; k<dim; k++)
 270             value_oppose(Mat->p[0][k],Q1->Constraint[i][k]);
 271         }
 272         else {
 273
 274           /* Copy constraint to Mat */
 275
 276           value_set_si(Mat->p[0][0],1);
 277           for (k=1; k<dim; k++)
 278             value_assign(Mat->p[0][k],Q1->Constraint[i][k]);
 279         }
 280       }
 281       else if(value_neg_p(Q1->Constraint[i][INDEX])) {
 282
 283         /* Upper bound -- make a lower bound from it */
 284         value_set_si(Mat->p[0][0],1);
 285         for (k=1; k<dim; k++)
 286           value_oppose(Mat->p[0][k],Q1->Constraint[i][k]);
 287       }
 288       else {
 289
 290         /* Lower or side bound -- ignore it */
 291         continue;
 292       }
 293
 294       /* Choose a constraint from Q2 */
 295       for (j=0; j<Q2->NbConstraints; j++) {
 296         if (value_zero_p(Q2->Constraint[j][0])) {   /* equality */
 297           if (value_zero_p(Q2->Constraint[j][INDEX])) {
 298
 299             /* Ignore side constraint (they are in Q) */
 300             continue;
 301           }
 302           else if (value_pos_p(Q2->Constraint[j][INDEX])) {
 303
 304             /* Copy -constraint to Mat */
 305             value_set_si(Mat->p[1][0],1);
 306             for (k=1; k<dim; k++)
 307               value_oppose(Mat->p[1][k],Q2->Constraint[j][k]);
 308           }
 309           else {
 310
 311             /* Copy constraint to Mat */
 312             value_set_si(Mat->p[1][0],1);
 313             for (k=1; k<dim; k++)
 314               value_assign(Mat->p[1][k],Q2->Constraint[j][k]);
 315           };
 316         }
 317         else if (value_pos_p(Q2->Constraint[j][INDEX])) {
 318
 319           /* Lower bound -- make an upper bound from it */
 320           value_set_si(Mat->p[1][0],1);
 321           for(k=1;k<dim;k++)
 322             value_oppose(Mat->p[1][k],Q2->Constraint[j][k]);
 323         }
 324         else {
 325
 326           /* Upper or side bound -- ignore it */
 327           continue;
 328         };
 329
 330 #ifdef DEBUG
 331         fprintf(stdout, "i=%d j=%d M=\n", i+1, j+1);
 332         Matrix_Print(stdout,P_VALUE_FMT,Mat);
 333 #endif
 334
 335         /* Add Mat to Q and see if anything is made */
 336         Q3 = AddConstraints(Mat->p[0],2,Q,NbMaxConstrs);
 337
 338 #ifdef DEBUG
 339         fprintf(stdout, "Q3\n");
 340         Polyhedron_Print(stdout,P_VALUE_FMT,Q3);
 341 #endif
 342
 343         if (!emptyQ(Q3)) {
 344           Domain_Free(Q3);
 345
 346 #ifdef DEBUG
 347           fprintf(stdout, "not empty\n");
 348 #endif
 349           res = -1;
 350           goto LTQdone;
 351         }
 352 #ifdef DEBUG
 353         fprintf(stdout,"empty\n");
 354 #endif
 355         Domain_Free(Q3);
 356       } /* end for j */
 357     } /* end for i */
 358     res = 1;
 359 LTQdone:
 360     Matrix_Free(Mat);
 361 LTQdone2:
 362     Domain_Free(Q4);
 363     Domain_Free(Q1);
 364     Domain_Free(Q2);
 365   }
 366   Domain_Free(Q);
 367
 368 #ifdef DEBUG
 369   fprintf(stdout, "res = %d\n", res);
 370 #endif
 371
 372   return res;
 373 } /* PolyhedronLTQ */
 374
 375 /* GaussSimplify --
 376    Given Mat1, a matrix of equalities, performs Gaussian elimination.
 377    Find a minimum basis, Returns the rank.
 378    Mat1 is context, Mat2 is reduced in context of Mat1
 379 */
 380 int GaussSimplify(Matrix *Mat1,Matrix *Mat2) {
 381
 382   int NbRows = Mat1->NbRows;
 383   int NbCols = Mat1->NbColumns;
 384   int *column_index;
 385   int i, j, k, n, t, pivot, Rank;
 386   Value gcd, tmp, *cp;
 387
 388   column_index=(int *)malloc(NbCols * sizeof(int));
 389   if (!column_index) {
 390     errormsg1("GaussSimplify", "outofmem", "out of memory space\n");
 391     Pol_status = 1;
 392     return 0;
 393   }
 394
 395   /* Initialize all the 'Value' variables */
 396   value_init(gcd); value_init(tmp);
 397
 398   Rank=0;
 399   for (j=0; j<NbCols; j++) {              /* for each column starting at */
 400     for (i=Rank; i<NbRows; i++)           /* diagonal, look down to find */
 401       if (value_notzero_p(Mat1->p[i][j])) /* the first non-zero entry    */
 402         break;
 403     if (i!=NbRows) {                      /* was one found ? */
 404       if (i!=Rank)                        /* was it found below the diagonal?*/
 405         Vector_Exchange(Mat1->p[Rank],Mat1->p[i],NbCols);
 406
 407       /* Normalize the pivot row */
 408       Vector_Gcd(Mat1->p[Rank],NbCols,&gcd);
 409
 410       /* If (gcd >= 2) */
 411       value_set_si(tmp,2);
 412       if (value_ge(gcd,tmp)) {
 413         cp = Mat1->p[Rank];
 414         for (k=0; k<NbCols; k++,cp++)
 415           value_division(*cp,*cp,gcd);
 416       }
 417       if (value_neg_p(Mat1->p[Rank][j])) {
 418         cp = Mat1->p[Rank];
 419         for (k=0; k<NbCols; k++,cp++)
 420           value_oppose(*cp,*cp);
 421       }
 422       /* End of normalize */
 423       pivot=i;
 424       for (i=0;i<NbRows;i++)    /* Zero out the rest of the column */
 425         if (i!=Rank) {
 426           if (value_notzero_p(Mat1->p[i][j])) {
 427             Value a, a1, a2, a1abs, a2abs;
 428             value_init(a); value_init(a1); value_init(a2);
 429             value_init(a1abs); value_init(a2abs);
 430             value_assign(a1,Mat1->p[i][j]);
 431             value_absolute(a1abs,a1);
 432             value_assign(a2,Mat1->p[Rank][j]);
 433             value_absolute(a2abs,a2);
 434             value_gcd(a, a1abs, a2abs);
 435             value_divexact(a1, a1, a);
 436             value_divexact(a2, a2, a);
 437             value_oppose(a1,a1);
 438             Vector_Combine(Mat1->p[i],Mat1->p[Rank],Mat1->p[i],a2,
 439                            a1,NbCols);
 440             Vector_Normalize(Mat1->p[i],NbCols);
 441             value_clear(a); value_clear(a1); value_clear(a2);
 442             value_clear(a1abs); value_clear(a2abs);
 443           }
 444         }
 445       column_index[Rank]=j;
 446       Rank++;
 447     }
 448   } /* end of Gauss elimination */
 449
 450
 451   if (Mat2) {  /* Mat2 is a transformation matrix  (i,j->f(i,j))....
 452                   can't scale it because can't scale both sides of -> */
 453     /* normalizes an affine transformation        */
 454     /* priority of forms                          */
 455     /*    1. i' -> i                (identity)    */
 456     /*    2. i' -> i + constant     (uniform)     */
 457     /*    3. i' -> constant         (broadcast)   */
 458     /*    4. i' -> j                (permutation) */
 459     /*    5. i' -> j + constant     (      )      */
 460     /*    6. i' -> i + j + constant (non-uniform) */
 461     for (k=0; k<Rank; k++) {
 462       j = column_index[k];
 463       for (i=0; i<(Mat2->NbRows-1);i++) {   /* all but the last row 0...0 1 */
 464         if ((i!=j) && value_notzero_p(Mat2->p[i][j])) {
 465
 466           /* Remove dependency of i' on j */
 467           Value a, a1, a1abs, a2, a2abs;
 468           value_init(a); value_init(a1); value_init(a2);
 469           value_init(a1abs); value_init(a2abs);
 470           value_assign(a1,Mat2->p[i][j]);
 471           value_absolute(a1abs,a1);
 472           value_assign(a2,Mat1->p[k][j]);
 473           value_absolute(a2abs,a2);
 474           value_gcd(a, a1abs, a2abs);
 475           value_divexact(a1, a1, a);
 476           value_divexact(a2, a2, a);
 477           value_oppose(a1,a1);
 478           if (value_one_p(a2)) {
 479             Vector_Combine(Mat2->p[i],Mat1->p[k],Mat2->p[i],a2,
 480                            a1,NbCols);
 481
 482             /* Vector_Normalize(Mat2->p[i],NbCols); -- can't do T        */
 483           } /* otherwise, can't do it without mult lhs prod (2i,3j->...) */
 484           value_clear(a); value_clear(a1); value_clear(a2);
 485           value_clear(a1abs); value_clear(a2abs);
 486
 487         }
 488         else if ((i==j) && value_zero_p(Mat2->p[i][j])) {
 489
 490           /* 'i' does not depend on j */
 491           for (n=j+1; n < (NbCols-1); n++) {
 492             if (value_notzero_p(Mat2->p[i][n])) { /* i' depends on some n */
 493               value_set_si(tmp,1);
 494               Vector_Combine(Mat2->p[i],Mat1->p[k],Mat2->p[i],tmp,
 495                              tmp,NbCols);
 496               break;
 497             }  /* if 'i' depends on just a constant, then leave it alone.*/
 498           }
 499         }
 500       }
 501     }
 502
 503     /* Check last row of transformation Mat2 */
 504     for (j=0; j<(NbCols-1); j++)
 505       if (value_notzero_p(Mat2->p[Mat2->NbRows-1][j])) {
 506         errormsg1("GaussSimplify", "corrtrans", "Corrupted transformation\n");
 507         break;
 508       }
 509
 510     if (value_notone_p(Mat2->p[Mat2->NbRows-1][NbCols-1])) {
 511       errormsg1("GaussSimplify", "corrtrans", "Corrupted transformation\n");
 512     }
 513   }
 514   value_clear(gcd); value_clear(tmp);
 515   free(column_index);
 516   return Rank;
 517 } /* GaussSimplify */
 518
 519 /*
 520  * Topologically sort 'n' polyhdera and return 0 on failure, otherwise return
 521  * 1 on success. Here 'L' is a an array of pointers to polyhedra, 'n' is the
 522  * number of polyhedra, 'index' is the level to consider for partial ordering
 523  * 'pdim' is the parameter space dimension, 'time' is an array of 'n' integers
 524  * to store logical time values, 'pvect', if not NULL, is an array of 'n'
 525  * integers that contains a permutation specification after call and MAXRAYS is
 526  * the workspace size for polyhedra operations.
 527  */
 528 int PolyhedronTSort (Polyhedron **L,unsigned int n,unsigned int index,unsigned int pdim,int *time,int *pvect,unsigned int MAXRAYS) {
 529
 530   unsigned int const nbcells = ((n*(n-1))>>1);    /* Number of memory cells
 531                                                      to allocate, see below */
 532   int *dag;  /* The upper triangular matrix */
 533   int **p;   /* Array of matrix row addresses */
 534   unsigned int i, j, k;
 535   unsigned int t, nb, isroot;
 536
 537   if (n<2) return 0;
 538
 539   /* we need an upper triangular matrix (example with n=4):
 540
 541      . o o o
 542      . . o o     . are unuseful cells, o are useful cells
 543      . . . o
 544      . . . .
 545
 546      so we need to allocate (n)(n-1)/2 cells
 547      - dag will point to this memory.
 548      - p[i] will point to row i of the matrix
 549      p[0] = dag - 1 (such that p[0][1] == dag[0])
 550      p[1] = dag - 1 + (n-1)
 551      p[2] = dag - 1 + (n-1) + (n-2)
 552      ...
 553      p[i] = p[i-1] + (n-1-i)
 554   */
 555
 556   /* malloc n(n-1)/2 for dag and n-1 for p (node n does not have any row) */
 557   dag = (int *) malloc(nbcells*sizeof(int));
 558   if (!dag) return 0;
 559   p = (int **) malloc ((n-1) * sizeof(int *));
 560   if (!p) {
 561     free(dag); return 0;
 562   }
 563
 564   /* Initialize 'p' and 'dag' */
 565   p[0] = dag-1;
 566   for (i=1; i<n-1; i++)
 567     p[i] = p[i-1] + (n-1)-i;
 568   for (i=0; i<nbcells; i++)
 569     dag[i] = -2;      /* -2 means 'not computed yet' */
 570   for (i=0; i<n; i++) time[i] = -1;
 571
 572   /* Compute the dag using transitivity to reduce the number of */
 573   /*   PolyhedronLTQ calls.                                     */
 574   for (i=0; i<n-1; i++) {
 575     POL_ENSURE_FACETS(L[i]);
 576     POL_ENSURE_VERTICES(L[i]);
 577     for (j=i+1; j<n; j++) {
 578       if (p[i][j] == -2) /* not computed yes */
 579         p[i][j] = PolyhedronLTQ(L[i], L[j], index, pdim, MAXRAYS);
 580       if (p[i][j] != 0) {
 581
 582         /* if p[i][j] is 1 or -1, look for -p[i][j] on the same row:
 583            p[i][j] == -p[i][k] ==> p[j][k] = p[i][k] (transitivity)
 584            note: p[r][c] == -p[c][r], use this to avoid reading or writing
 585            under the matrix diagonal
 586         */
 587
 588         /* first, k<i so look for p[i][j] == p[k][i] (i.e. -p[i][k]) */
 589         for (k=0; k<i; k++)
 590           if (p[k][i] == p[i][j]) p[k][j] = p[k][i];
 591
 592         /* then, i<k<j so look for p[i][j] == -p[i][k] */
 593         for (k=i+1; k<j; k++)
 594           if (p[i][k] == -p[i][j]) p[k][j] = -p[i][k];
 595
 596         /* last, k>j same search but */
 597         for (k=j+1; k<n; k++)
 598           if (p[i][k] == -p[i][j]) p[j][k] = p[i][k];
 599       }
 600     }
 601   }
 602
 603   /* walk thru the dag to compute the partial order (and optionally
 604      the permutation)
 605      Note: this is not the fastest way to do it but it takes
 606      negligible time compared to a single call of PolyhedronLTQ !
 607      Each macro-step (while loop) gives the same logical time to all
 608      found roots and optionally add these nodes in the permutation vector
 609   */
 610
 611   t = 0; /* current logical time, assigned to current roots and
 612             increased by 1 at the end of each step */
 613   nb = 0; /* number of processed nodes (have been given a time) */
 614   while (nb<n) {
 615     for (i=0; i<n; i++) { /* search for roots */
 616       /* for any node, if it's not already been given a logical time
 617          then walk thru the node row; if we find a 1 at some column j,
 618          it means that node j preceeds the current node, so it is not
 619          a root */
 620       if (time[i]<0) {
 621         isroot = 1; /* assume that it is until it is definitely not */
 622         /* first search on a column */
 623         for (j=0; j<i; j++) {
 624           if (p[j][i]==-1) { /* found a node lower than it */
 625             isroot = 0; break;
 626           }
 627         }
 628         if /*still*/ (isroot)
 629           for (j=i+1; j<n; j++) {
 630             if (p[i][j]==1) { /* greater than this node */
 631               isroot = 0; break;
 632             }
 633           }
 634         if (isroot) { /* then it definitely is */
 635           time[i] = t; /* this node gets current logical time */
 636           if (pvect)
 637             pvect[nb] = i+1; /* nodes will be numbered from 1 to n */
 638           nb++; /* one more node processed */
 639         }
 640       }
 641     }
 642     /* now make roots become neutral, i.e. non comparable with all other nodes */
 643     for (i=0; i<n; i++) {
 644       if (time[i]==t) {
 645         for (j=0; j<i; j++)
 646           p[j][i] = 0;
 647         for (j=i+1; j<n; j++)
 648           p[i][j] = 0;
 649       }
 650     }
 651     t++; /* ready for next set of root nodes */
 652   }
 653
 654   free (p);   /* let's be clean, collect the garbage */
 655   free (dag);
 656   return 1;
 657 } /* PolyhedronTSort */
 658
 659
 660
 661