polylib: alpha.c Source File

00001 /* polytest.c */
00002 #include <stdio.h>
00003 #include <polylib/polylib.h>
00004 
00005 /* alpha.c
00006      COPYRIGHT
00007           Both this software and its documentation are
00008 
00009               Copyright 1993 by IRISA /Universite de Rennes I - France,
00010               Copyright 1995,1996 by BYU, Provo, Utah
00011                          all rights reserved.
00012 
00013           Permission is granted to copy, use, and distribute
00014           for any commercial or noncommercial purpose under the terms
00015           of the GNU General Public license, version 2, June 1991
00016           (see file : LICENSING).
00017 */
00018 
00019 #include <stdio.h>
00020 #include <stdlib.h>
00021 #include <string.h>
00022 
00023 /*---------------------------------------------------------------------*/
00024 /* int exist_points(pos,P,context)                                     */
00025 /*    pos : index position of current loop index (0..hdim-1)           */
00026 /*    P: loop domain                                                   */
00027 /*    context : context values for fixed indices                       */
00028 /* recursive procedure, recurs for each imbriquation                   */
00029 /* returns 1 if there exists any integer points in this polyhedron     */
00030 /* returns 0 if there are none                                         */
00031 /*---------------------------------------------------------------------*/
00032 static int exist_points(int pos,Polyhedron *Pol,Value *context) {
00033   
00034   Value LB, UB, k,tmp;
00035   
00036   value_init(LB); value_init(UB); 
00037   value_init(k);  value_init(tmp);
00038   value_set_si(LB,0);
00039   value_set_si(UB,0);
00040   
00041   /* Problem if UB or LB is INFINITY */
00042   if (lower_upper_bounds(pos,Pol,context,&LB,&UB) !=0) {
00043     errormsg1("exist_points", "infdom", "infinite domain");
00044     value_clear(LB);
00045     value_clear(UB);
00046     value_clear(k);
00047     value_clear(tmp);
00048     return -1;
00049   }
00050   value_set_si(context[pos],0);
00051   if(value_lt(UB,LB)) {
00052     value_clear(LB); 
00053     value_clear(UB);
00054     value_clear(k);
00055     value_clear(tmp);
00056     return 0;
00057   }  
00058   if (!Pol->next) {
00059     value_subtract(tmp,UB,LB);
00060     value_increment(tmp,tmp);
00061     value_clear(UB);
00062     value_clear(LB);
00063     value_clear(k);
00064     return (value_pos_p(tmp));
00065   }
00066   
00067   for (value_assign(k,LB);value_le(k,UB);value_increment(k,k)) {
00068     
00069     /* insert k in context */
00070     value_assign(context[pos],k);    
00071     if (exist_points(pos+1,Pol->next,context) > 0 ) {
00072       value_clear(LB); value_clear(UB);
00073       value_clear(k); value_clear(tmp);
00074       return 1;
00075     }
00076   }   
00077   /* Reset context */
00078   value_set_si(context[pos],0);
00079   value_clear(UB); value_clear(LB);
00080   value_clear(k); value_clear(tmp);
00081   return 0;
00082 }
00083     
00084 /*--------------------------------------------------------------*/
00085 /* Test to see if there are any integral points in a polyhedron */
00086 /*--------------------------------------------------------------*/
00087 int Polyhedron_Not_Empty(Polyhedron *P,Polyhedron *C,int MAXRAYS) {
00088 
00089   int res,i;
00090   Value *context;
00091   Polyhedron *L;
00092   
00093   POL_ENSURE_FACETS(P);
00094   POL_ENSURE_VERTICES(P);
00095   POL_ENSURE_FACETS(C);
00096   POL_ENSURE_VERTICES(C);
00097 
00098   /* Create a context vector size dim+2 and set it to all zeros */
00099   context = (Value *) malloc((P->Dimension+2)*sizeof(Value));
00100   
00101   /* Initialize array 'context' */
00102   for (i=0;i<(P->Dimension+2);i++) 
00103     value_init(context[i]);
00104   
00105   Vector_Set(context,0,(P->Dimension+2));
00106   
00107   /* Set context[P->Dimension+1] = 1  (the constant) */
00108   value_set_si(context[P->Dimension+1],1);
00109   
00110   L = Polyhedron_Scan(P,C,MAXRAYS);
00111   res = exist_points(1,L,context);
00112   Domain_Free(L);
00113   
00114   /* Clear array 'context' */
00115   for (i=0;i<(P->Dimension+2);i++) 
00116     value_clear(context[i]);
00117   free(context);
00118   return res;
00119 }
00120 
00121 /* PolyhedronLTQ ( P1, P2 ) */
00122 /* P1 and P2 must be simple polyhedra */
00123 /* result =  0 --> not comparable */
00124 /* result = -1 --> P1 < P2        */
00125 /* result =  1 --> P1 > P2        */
00126 /* INDEX  = 1 .... Dimension      */
00127 int PolyhedronLTQ (Polyhedron *Pol1,Polyhedron *Pol2,int INDEX, int PDIM, int NbMaxConstrs) { 
00128   
00129   int res, dim, i, j, k;
00130   Polyhedron *Q1, *Q2, *Q3, *Q4, *Q;
00131   Matrix *Mat;
00132 
00133   if (Pol1->next || Pol2->next) {
00134     errormsg1("PolyhedronLTQ", "compoly", "Can only compare polyhedra");
00135     return 0;
00136   }
00137   if (Pol1->Dimension != Pol2->Dimension) {
00138     errormsg1("PolyhedronLTQ","diffdim","Polyhedra are not same dimension");
00139     return 0;
00140   }
00141   dim = Pol1->Dimension+2;
00142 
00143   POL_ENSURE_FACETS(Pol1);
00144   POL_ENSURE_VERTICES(Pol1);
00145   POL_ENSURE_FACETS(Pol2);
00146   POL_ENSURE_VERTICES(Pol2);
00147   
00148 #ifdef DEBUG
00149   fprintf(stdout, "P1\n");
00150   Polyhedron_Print(stdout,P_VALUE_FMT,Pol1);
00151   fprintf(stdout, "P2\n");
00152   Polyhedron_Print(stdout,P_VALUE_FMT,Pol2);
00153 #endif
00154   
00155   /* Create the Line to add */
00156   k = Pol1->Dimension-INDEX+1-PDIM;
00157   Mat = Matrix_Alloc(k,dim);
00158   Vector_Set(Mat->p_Init,0,dim*k);
00159   for(j=0,i=INDEX;j<k;i++,j++)
00160     value_set_si(Mat->p[j][i],1);
00161   
00162   Q1 = AddRays(Mat->p[0],k,Pol1,NbMaxConstrs);
00163   Q2 = AddRays(Mat->p[0],k,Pol2,NbMaxConstrs);
00164 
00165 #ifdef DEBUG
00166   fprintf(stdout, "Q1\n");
00167   Polyhedron_Print(stdout,P_VALUE_FMT,Q1);
00168   fprintf(stdout, "Q2\n");
00169   Polyhedron_Print(stdout,P_VALUE_FMT,Q2);
00170 #endif
00171   
00172   Matrix_Free(Mat);
00173   Q  = DomainIntersection(Q1,Q2,NbMaxConstrs);
00174   
00175 #ifdef DEBUG
00176   fprintf(stdout, "Q\n");
00177   Polyhedron_Print(stdout,P_VALUE_FMT,Q);
00178 #endif
00179   
00180   Domain_Free(Q1);
00181   Domain_Free(Q2);
00182   
00183   if (emptyQ(Q)) res = 0;       /* not comparable */
00184   else {
00185     Q1 = DomainIntersection(Pol1,Q,NbMaxConstrs);
00186     Q2 = DomainIntersection(Pol2,Q,NbMaxConstrs);
00187     
00188 #ifdef DEBUG
00189     fprintf(stdout, "Q1\n");
00190     Polyhedron_Print(stdout,P_VALUE_FMT,Q1);
00191     fprintf(stdout, "Q2\n");
00192     Polyhedron_Print(stdout,P_VALUE_FMT,Q2);
00193 #endif
00194 
00195     k = Q1->NbConstraints + Q2->NbConstraints;
00196     Mat = Matrix_Alloc(k, dim);
00197     Vector_Set(Mat->p_Init,0,k*dim);
00198     
00199     /* First compute surrounding polyhedron */    
00200     j=0;
00201     for (i=0; i<Q1->NbConstraints; i++) {
00202       if ((value_one_p(Q1->Constraint[i][0])) && (value_pos_p(Q1->Constraint[i][INDEX]))) {
00203         
00204         /* keep Q1's lower bounds */
00205         for (k=0; k<dim; k++) 
00206           value_assign(Mat->p[j][k],Q1->Constraint[i][k]);
00207         j++;
00208       }
00209     }
00210     for (i=0; i<Q2->NbConstraints; i++) {
00211       if ((value_one_p(Q2->Constraint[i][0])) && (value_neg_p(Q2->Constraint[i][INDEX]))) {
00212         
00213         /* and keep Q2's upper bounds */
00214         for (k=0; k<dim; k++) 
00215           value_assign(Mat->p[j][k],Q2->Constraint[i][k]);
00216         j++;
00217       }
00218     }
00219     Q4 = AddConstraints(Mat->p[0], j, Q, NbMaxConstrs);
00220     Matrix_Free(Mat);
00221     
00222 #ifdef debug
00223     fprintf(stderr, "Q4 surrounding polyhedron\n");
00224     Polyhderon_Print(stderr,P_VALUE_FMT, Q4);
00225 #endif
00226 
00227     /* if surrounding polyhedron is empty, D1>D2 */
00228     if (emptyQ(Q4)) {
00229       res = 1;
00230       
00231 #ifdef debug
00232       fprintf(stderr, "Surrounding polyhedron is empty\n");
00233 #endif
00234       goto LTQdone2; 
00235     }
00236     
00237     /* Test if Q1 < Q2 */      
00238     /* Build a constraint array for >= Q1 and <= Q2 */
00239     Mat = Matrix_Alloc(2,dim);
00240     Vector_Set(Mat->p_Init,0,2*dim);
00241     
00242     /* Choose a contraint from Q1 */
00243     for (i=0; i<Q1->NbConstraints; i++) {
00244       if (value_zero_p(Q1->Constraint[i][0])) {
00245         
00246         /* Equality */
00247         if (value_zero_p(Q1->Constraint[i][INDEX])) {
00248           
00249           /* Ignore side constraint (they are in Q) */
00250           continue;
00251         }
00252         else if (value_neg_p(Q1->Constraint[i][INDEX])) {
00253           
00254           /* copy -constraint to Mat */
00255           value_set_si(Mat->p[0][0],1);
00256           for (k=1; k<dim; k++)
00257             value_oppose(Mat->p[0][k],Q1->Constraint[i][k]);
00258         }
00259         else {
00260           
00261           /* Copy constraint to Mat */
00262           
00263           value_set_si(Mat->p[0][0],1);
00264           for (k=1; k<dim; k++)
00265             value_assign(Mat->p[0][k],Q1->Constraint[i][k]);
00266         }
00267       }
00268       else if(value_neg_p(Q1->Constraint[i][INDEX])) {
00269         
00270         /* Upper bound -- make a lower bound from it */
00271         value_set_si(Mat->p[0][0],1);
00272         for (k=1; k<dim; k++)
00273           value_oppose(Mat->p[0][k],Q1->Constraint[i][k]);
00274       }
00275       else {    
00276         
00277         /* Lower or side bound -- ignore it */
00278         continue;
00279       }
00280       
00281       /* Choose a constraint from Q2 */
00282       for (j=0; j<Q2->NbConstraints; j++) {
00283         if (value_zero_p(Q2->Constraint[j][0])) {   /* equality */
00284           if (value_zero_p(Q2->Constraint[j][INDEX])) {
00285             
00286             /* Ignore side constraint (they are in Q) */
00287             continue;
00288           }
00289           else if (value_pos_p(Q2->Constraint[j][INDEX])) {
00290             
00291             /* Copy -constraint to Mat */
00292             value_set_si(Mat->p[1][0],1);
00293             for (k=1; k<dim; k++)
00294               value_oppose(Mat->p[1][k],Q2->Constraint[j][k]);
00295           }
00296           else {
00297             
00298             /* Copy constraint to Mat */
00299             value_set_si(Mat->p[1][0],1);
00300             for (k=1; k<dim; k++)
00301               value_assign(Mat->p[1][k],Q2->Constraint[j][k]);
00302           };
00303         }
00304         else if (value_pos_p(Q2->Constraint[j][INDEX])) {
00305           
00306           /* Lower bound -- make an upper bound from it */
00307           value_set_si(Mat->p[1][0],1);
00308           for(k=1;k<dim;k++)
00309             value_oppose(Mat->p[1][k],Q2->Constraint[j][k]);
00310         }
00311         else {
00312           
00313           /* Upper or side bound -- ignore it */
00314           continue;
00315         };
00316         
00317 #ifdef DEBUG
00318         fprintf(stdout, "i=%d j=%d M=\n", i+1, j+1);
00319         Matrix_Print(stdout,P_VALUE_FMT,Mat);
00320 #endif
00321         
00322         /* Add Mat to Q and see if anything is made */
00323         Q3 = AddConstraints(Mat->p[0],2,Q,NbMaxConstrs);
00324 
00325 #ifdef DEBUG
00326         fprintf(stdout, "Q3\n");
00327         Polyhedron_Print(stdout,P_VALUE_FMT,Q3);
00328 #endif
00329         
00330         if (!emptyQ(Q3)) { 
00331           Domain_Free(Q3);
00332           
00333 #ifdef DEBUG
00334           fprintf(stdout, "not empty\n");
00335 #endif
00336           res = -1;
00337           goto LTQdone;
00338         }
00339 #ifdef DEBUG
00340         fprintf(stdout,"empty\n");      
00341 #endif
00342         Domain_Free(Q3);
00343       } /* end for j */
00344     } /* end for i */
00345     res = 1;
00346 LTQdone:
00347     Matrix_Free(Mat);
00348 LTQdone2: 
00349     Domain_Free(Q4);
00350     Domain_Free(Q1);
00351     Domain_Free(Q2);
00352   }
00353   Domain_Free(Q);
00354   
00355 #ifdef DEBUG
00356   fprintf(stdout, "res = %d\n", res);
00357 #endif
00358   
00359   return res;
00360 } /* PolyhedronLTQ */
00361 
00362 /* GaussSimplify --
00363    Given Mat1, a matrix of equalities, performs Gaussian elimination.
00364    Find a minimum basis, Returns the rank.
00365    Mat1 is context, Mat2 is reduced in context of Mat1
00366 */
00367 int GaussSimplify(Matrix *Mat1,Matrix *Mat2) {
00368   
00369   int NbRows = Mat1->NbRows;
00370   int NbCols = Mat1->NbColumns;
00371   int *column_index;
00372   int i, j, k, n, t, pivot, Rank; 
00373   Value gcd, tmp, *cp; 
00374   
00375   column_index=(int *)malloc(NbCols * sizeof(int));
00376   if (!column_index) {
00377     errormsg1("GaussSimplify", "outofmem", "out of memory space\n");
00378     Pol_status = 1;
00379     return 0;
00380   }
00381   
00382   /* Initialize all the 'Value' variables */
00383   value_init(gcd); value_init(tmp);
00384   
00385   Rank=0;
00386   for (j=0; j<NbCols; j++) {              /* for each column starting at */ 
00387     for (i=Rank; i<NbRows; i++)           /* diagonal, look down to find */
00388       if (value_notzero_p(Mat1->p[i][j])) /* the first non-zero entry    */
00389         break;                           
00390     if (i!=NbRows) {                      /* was one found ? */
00391       if (i!=Rank)                        /* was it found below the diagonal?*/
00392         Vector_Exchange(Mat1->p[Rank],Mat1->p[i],NbCols);
00393       
00394       /* Normalize the pivot row */
00395       Vector_Gcd(Mat1->p[Rank],NbCols,&gcd);
00396       
00397       /* If (gcd >= 2) */
00398       value_set_si(tmp,2);
00399       if (value_ge(gcd,tmp)) {
00400         cp = Mat1->p[Rank];
00401         for (k=0; k<NbCols; k++,cp++)
00402           value_division(*cp,*cp,gcd);          
00403       }
00404       if (value_neg_p(Mat1->p[Rank][j])) {
00405         cp = Mat1->p[Rank];
00406         for (k=0; k<NbCols; k++,cp++)
00407           value_oppose(*cp,*cp);
00408       }
00409       /* End of normalize */
00410       pivot=i;
00411       for (i=0;i<NbRows;i++)    /* Zero out the rest of the column */
00412         if (i!=Rank) {
00413           if (value_notzero_p(Mat1->p[i][j])) {
00414             Value a, a1, a2, a1abs, a2abs;
00415             value_init(a); value_init(a1); value_init(a2);
00416             value_init(a1abs); value_init(a2abs);
00417             value_assign(a1,Mat1->p[i][j]);
00418             value_absolute(a1abs,a1);
00419             value_assign(a2,Mat1->p[Rank][j]); 
00420             value_absolute(a2abs,a2);
00421             value_gcd(a, a1abs, a2abs);
00422             value_divexact(a1, a1, a);
00423             value_divexact(a2, a2, a);
00424             value_oppose(a1,a1);
00425             Vector_Combine(Mat1->p[i],Mat1->p[Rank],Mat1->p[i],a2, 
00426                            a1,NbCols);
00427             Vector_Normalize(Mat1->p[i],NbCols);
00428             value_clear(a); value_clear(a1); value_clear(a2);
00429             value_clear(a1abs); value_clear(a2abs);
00430           }
00431         }
00432       column_index[Rank]=j;
00433       Rank++;
00434     }
00435   } /* end of Gauss elimination */
00436 
00437 
00438   if (Mat2) {  /* Mat2 is a transformation matrix  (i,j->f(i,j))....
00439                   can't scale it because can't scale both sides of -> */
00440     /* normalizes an affine transformation        */
00441     /* priority of forms                          */
00442     /*    1. i' -> i                (identity)    */
00443     /*    2. i' -> i + constant     (uniform)     */
00444     /*    3. i' -> constant         (broadcast)   */
00445     /*    4. i' -> j                (permutation) */
00446     /*    5. i' -> j + constant     (      )      */
00447     /*    6. i' -> i + j + constant (non-uniform) */
00448     for (k=0; k<Rank; k++) {
00449       j = column_index[k];
00450       for (i=0; i<(Mat2->NbRows-1);i++) {   /* all but the last row 0...0 1 */
00451         if ((i!=j) && value_notzero_p(Mat2->p[i][j])) {
00452           
00453           /* Remove dependency of i' on j */
00454           Value a, a1, a1abs, a2, a2abs;
00455           value_init(a); value_init(a1); value_init(a2);
00456           value_init(a1abs); value_init(a2abs);
00457           value_assign(a1,Mat2->p[i][j]);
00458           value_absolute(a1abs,a1);
00459           value_assign(a2,Mat1->p[k][j]);
00460           value_absolute(a2abs,a2);
00461           value_gcd(a, a1abs, a2abs);
00462           value_divexact(a1, a1, a);
00463           value_divexact(a2, a2, a);
00464           value_oppose(a1,a1);
00465           if (value_one_p(a2)) {
00466             Vector_Combine(Mat2->p[i],Mat1->p[k],Mat2->p[i],a2,
00467                            a1,NbCols);
00468             
00469             /* Vector_Normalize(Mat2->p[i],NbCols); -- can't do T        */
00470           } /* otherwise, can't do it without mult lhs prod (2i,3j->...) */
00471           value_clear(a); value_clear(a1); value_clear(a2);
00472           value_clear(a1abs); value_clear(a2abs);
00473                 
00474         }
00475         else if ((i==j) && value_zero_p(Mat2->p[i][j])) {
00476           
00477           /* 'i' does not depend on j */
00478           for (n=j+1; n < (NbCols-1); n++) {
00479             if (value_notzero_p(Mat2->p[i][n])) { /* i' depends on some n */
00480               value_set_si(tmp,1);
00481               Vector_Combine(Mat2->p[i],Mat1->p[k],Mat2->p[i],tmp,
00482                              tmp,NbCols);
00483               break;
00484             }  /* if 'i' depends on just a constant, then leave it alone.*/
00485           }
00486         }
00487       }
00488     }
00489     
00490     /* Check last row of transformation Mat2 */
00491     for (j=0; j<(NbCols-1); j++)
00492       if (value_notzero_p(Mat2->p[Mat2->NbRows-1][j])) {
00493         errormsg1("GaussSimplify", "corrtrans", "Corrupted transformation\n");
00494         break;
00495       }
00496     
00497     if (value_notone_p(Mat2->p[Mat2->NbRows-1][NbCols-1])) {
00498       errormsg1("GaussSimplify", "corrtrans", "Corrupted transformation\n");
00499     }
00500   }
00501   value_clear(gcd); value_clear(tmp);
00502   free(column_index);
00503   return Rank;
00504 } /* GaussSimplify */
00505 
00506 /* 
00507  * Topologically sort 'n' polyhdera and return 0 on failure, otherwise return 
00508  * 1 on success. Here 'L' is a an array of pointers to polyhedra, 'n' is the 
00509  * number of polyhedra, 'index' is the level to consider for partial ordering
00510  * 'pdim' is the parameter space dimension, 'time' is an array of 'n' integers
00511  * to store logical time values, 'pvect', if not NULL, is an array of 'n' 
00512  * integers that contains a permutation specification after call and MAXRAYS is
00513  * the workspace size for polyhedra operations. 
00514  */
00515 int PolyhedronTSort (Polyhedron **L,unsigned int n,unsigned int index,unsigned int pdim,int *time,int *pvect,unsigned int MAXRAYS) {
00516  
00517   unsigned int const nbcells = ((n*(n-1))>>1);    /* Number of memory cells 
00518                                                      to allocate, see below */
00519   int *dag;  /* The upper triangular matrix */
00520   int **p;   /* Array of matrix row addresses */
00521   unsigned int i, j, k;
00522   unsigned int t, nb, isroot;
00523 
00524   if (n<2) return 0;
00525 
00526   /* we need an upper triangular matrix (example with n=4):
00527 
00528      . o o o
00529      . . o o     . are unuseful cells, o are useful cells
00530      . . . o
00531      . . . .
00532      
00533      so we need to allocate (n)(n-1)/2 cells
00534      - dag will point to this memory.
00535      - p[i] will point to row i of the matrix
00536      p[0] = dag - 1 (such that p[0][1] == dag[0])
00537      p[1] = dag - 1 + (n-1)
00538      p[2] = dag - 1 + (n-1) + (n-2)
00539      ...
00540      p[i] = p[i-1] + (n-1-i)
00541   */
00542 
00543   /* malloc n(n-1)/2 for dag and n-1 for p (node n does not have any row) */
00544   dag = (int *) malloc(nbcells*sizeof(int));
00545   if (!dag) return 0;
00546   p = (int **) malloc ((n-1) * sizeof(int *));
00547   if (!p) { 
00548     free(dag); return 0; 
00549   }
00550 
00551   /* Initialize 'p' and 'dag' */
00552   p[0] = dag-1;
00553   for (i=1; i<n-1; i++)
00554     p[i] = p[i-1] + (n-1)-i;
00555   for (i=0; i<nbcells; i++)
00556     dag[i] = -2;      /* -2 means 'not computed yet' */
00557   for (i=0; i<n; i++) time[i] = -1;
00558 
00559   /* Compute the dag using transitivity to reduce the number of */
00560   /*   PolyhedronLTQ calls.                                     */
00561   for (i=0; i<n-1; i++) {
00562     POL_ENSURE_FACETS(L[i]);
00563     POL_ENSURE_VERTICES(L[i]);
00564     for (j=i+1; j<n; j++) {
00565       if (p[i][j] == -2) /* not computed yes */
00566         p[i][j] = PolyhedronLTQ(L[i], L[j], index, pdim, MAXRAYS);
00567       if (p[i][j] != 0) {
00568         
00569         /* if p[i][j] is 1 or -1, look for -p[i][j] on the same row:
00570            p[i][j] == -p[i][k] ==> p[j][k] = p[i][k] (transitivity)
00571            note: p[r][c] == -p[c][r], use this to avoid reading or writing
00572            under the matrix diagonal
00573         */  
00574            
00575         /* first, k<i so look for p[i][j] == p[k][i] (i.e. -p[i][k]) */
00576         for (k=0; k<i; k++)
00577           if (p[k][i] == p[i][j]) p[k][j] = p[k][i];
00578         
00579         /* then, i<k<j so look for p[i][j] == -p[i][k] */
00580         for (k=i+1; k<j; k++)
00581           if (p[i][k] == -p[i][j]) p[k][j] = -p[i][k];
00582         
00583         /* last, k>j same search but */
00584         for (k=j+1; k<n; k++)
00585           if (p[i][k] == -p[i][j]) p[j][k] = p[i][k];
00586       }
00587     }
00588   }
00589   
00590   /* walk thru the dag to compute the partial order (and optionally
00591      the permutation)
00592      Note: this is not the fastest way to do it but it takes
00593      negligible time compared to a single call of PolyhedronLTQ !
00594      Each macro-step (while loop) gives the same logical time to all
00595      found roots and optionally add these nodes in the permutation vector
00596   */ 
00597   
00598   t = 0; /* current logical time, assigned to current roots and
00599             increased by 1 at the end of each step */
00600   nb = 0; /* number of processed nodes (have been given a time) */
00601   while (nb<n) {
00602     for (i=0; i<n; i++) { /* search for roots */
00603       /* for any node, if it's not already been given a logical time
00604          then walk thru the node row; if we find a 1 at some column j,
00605          it means that node j preceeds the current node, so it is not
00606          a root */
00607       if (time[i]<0) {
00608         isroot = 1; /* assume that it is until it is definitely not */
00609         /* first search on a column */
00610         for (j=0; j<i; j++) {
00611           if (p[j][i]==-1) { /* found a node lower than it */
00612             isroot = 0; break;
00613           }
00614         }
00615         if /*still*/ (isroot)
00616           for (j=i+1; j<n; j++) {
00617             if (p[i][j]==1) { /* greater than this node */
00618               isroot = 0; break;
00619             }
00620           }
00621         if (isroot) { /* then it definitely is */
00622           time[i] = t; /* this node gets current logical time */
00623           if (pvect)
00624             pvect[nb] = i+1; /* nodes will be numbered from 1 to n */
00625           nb++; /* one more node processed */
00626         }
00627       }
00628     }
00629     /* now make roots become neutral, i.e. non comparable with all other nodes */
00630     for (i=0; i<n; i++) {
00631       if (time[i]==t) {
00632         for (j=0; j<i; j++)
00633           p[j][i] = 0;
00634         for (j=i+1; j<n; j++)
00635           p[i][j] = 0;
00636       }
00637     }
00638     t++; /* ready for next set of root nodes */
00639   }
00640   
00641   free (p);   /* let's be clean, collect the garbage */
00642   free (dag);
00643   return 1;
00644 } /* PolyhedronTSort */
00645 
00646 
00647 
00648