source/kernel/alpha.c

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