?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
corrected a bug in reading matrices on MacOS (fgets) and cleaned a bit the code
[polylib.git] / source / kernel / polyparam.c
blob46e4c88633a3502ae77441dfb61de9641baf2e0f
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 /***********************************************************************/
19 /* Parametrized polyhedra V4.20 */
20 /* copyright 1995-2000 Vincent Loechner */
21 /* copyright 1996-1997, Doran Wilde */
22 /***********************************************************************/
23
24 /********************* -----------USER #DEFS-------- ***********************/
25 /* These are mainly for debug purposes. You shouldn't need to change */
26 /* anything for daily usage... */
27 /***************************************************************************/
28
29 /* you may define each macro independently */
30 /* #define DEBUGPP */
31 /* #define DEBUGPP3 */ /* initialization of domain, context, ... */
32 /* #define DEBUGPP31 */ /* even more init-domains */
33 /* #define DEBUGPP32 */ /* even even more... (Elim_Columns) */
34 /* #define DEBUGPP4 */ /* m-faces scan */
35 /* #define DEBUGPP41 */ /* inverse Di in scan */
36 /* #define DEBUGPP5 */ /* Compute_PDomains */
37 /********************* ---------END USER #DEFS------ ***********************/
38
39 #include <stdio.h>
40 #include <stdlib.h>
41 #include <string.h>
42 #include <assert.h>
43 #ifdef DEBUGPP
44 #include <time.h>
45 #endif
46
47 #include <polylib/polylib.h>
48
49 static void traite_m_face(Polyhedron *, unsigned int *, unsigned int *);
50 static void scan_m_face(int,int,Polyhedron *,unsigned int *);
51
52 /*
53 * Return the intersection of two polyhedral domains 'Pol1' and 'Pol2' such
54 * that if the intersection is a polyhedron of lower dimension (a degenerate
55 * polyhedron) than the operands, it is discarded from the resulting polyhedra
56 * list.
57 */
58 Polyhedron *PDomainIntersection(Polyhedron *Pol1,Polyhedron *Pol2,unsigned NbMaxRays) {
59
60 Polyhedron *p1, *p2, *p3, *d;
61
62 if (!Pol1 || !Pol2) return (Polyhedron*) 0;
63 if((Pol1->Dimension != Pol2->Dimension) || (Pol1->NbEq != Pol2->NbEq)) {
64 fprintf(stderr,
65 "? PDomainIntersection: operation on different dimensions\n");
66 return (Polyhedron*) 0;
67 }
68
69 POL_ENSURE_FACETS(Pol1);
70 POL_ENSURE_VERTICES(Pol1);
71 POL_ENSURE_FACETS(Pol2);
72 POL_ENSURE_VERTICES(Pol2);
73
74 d = (Polyhedron *)0;
75 for (p1=Pol1; p1; p1=p1->next) {
76 for (p2=Pol2; p2; p2=p2->next) {
77 p3 = AddConstraints(p2->Constraint[0],
78 p2->NbConstraints,p1,NbMaxRays);
79 if (!p3) continue;
80
81 /* If the new polyhedron 'p3' has lower dimension, discard it */
82 if (p3->NbEq!=Pol1->NbEq)
83 Polyhedron_Free(p3) ;
84
85 /* Otherwise add it to the new polyhderal domain 'd'. */
86 else
87 d = AddPolyToDomain(p3,d);
88 }
89 }
90 return d;
91 } /* PDomainIntersection */
92
93 /*
94 * Given polyhderal domains 'Pol1' and 'Pol2', return the difference of the
95 * two domains with a modification that the resulting polyhedra in the new
96 * domain don't have a 1 unit space around cut and the degenerate results
97 * (of smaller dimension) are discarded.
98 */
99 Polyhedron *PDomainDifference(Polyhedron *Pol1,Polyhedron *Pol2,unsigned NbMaxRays) {
100
101 Polyhedron *p1, *p2, *p3, *d;
102 int i;
103
104 if (!Pol1 || !Pol2)
105 return (Polyhedron*) 0;
106 if((Pol1->Dimension != Pol2->Dimension) || (Pol1->NbEq != Pol2->NbEq)) {
107 fprintf(stderr,
108 "? PDomainDifference: operation on different dimensions\n");
109 return (Polyhedron*) 0;
110 }
111
112 POL_ENSURE_FACETS(Pol1);
113 POL_ENSURE_VERTICES(Pol1);
114 POL_ENSURE_FACETS(Pol2);
115 POL_ENSURE_VERTICES(Pol2);
116
117 d = (Polyhedron *)0;
118 for (p2=Pol2; p2; p2=p2->next) {
119 for (p1=Pol1; p1; p1=p1->next) {
120 for (i=0; i<p2->NbConstraints; i++) {
121
122 /* Add the constraint (-p2->Constraint[i]) >= 0 in 'p1' */
123 p3 = SubConstraint(p2->Constraint[i],p1,NbMaxRays,2);
124 if (!p3) continue;
125
126 /* If the new polyhedron 'p3' is empty or is a polyhedron of lower */
127 /* dimension, discard it. */
128 if (emptyQ(p3) || p3->NbEq!=Pol1->NbEq)
129 Polyhedron_Free(p3);
130
131 /* Otherwise add 'p3' to the new polyhderal domain 'd' */
132 else
133 d = AddPolyToDomain(p3,d);
134 }
135 }
136 if (p2 != Pol2)
137 Domain_Free(Pol1);
138 Pol1 = d;
139 d = (Polyhedron *)0;
140 }
141 return Pol1;
142 } /* PDomainDifference */
143
144 /*
145 * Return 1 if matrix 'Mat' is full column ranked, otherwise return 0.
146 */
147 static int TestRank(Matrix *Mat) {
148
149 int i,j,k;
150 Value m1,m2,m3,gcd,tmp;
151
152 /* Initialize all the 'Value' variables */
153 value_init(m1); value_init(m2);
154 value_init(m3); value_init(gcd); value_init(tmp);
155
156 for(k=0;k<Mat->NbColumns;++k) {
157
158 /* If the digonal entry (k,k) is zero, search down the column(k) */
159 /* starting from row(k) to find a non-zero entry */
160 if(value_zero_p(Mat->p[k][k])) {
161 for(j=k+1;j<Mat->NbRows;++j) {
162
163 /* If a non-zero entry (j,k) is found */
164 if(value_notzero_p(Mat->p[j][k])) {
165
166 /* Exchange row(k) and row(j) */
167 for(i=k;i<Mat->NbColumns;++i) {
168 value_assign(tmp,Mat->p[j][i]);
169 value_assign(Mat->p[j][i],Mat->p[k][i]);
170 value_assign(Mat->p[k][i],tmp);
171 }
172 break;
173 }
174 }
175
176 /* If no non-zero entry is found then the matrix 'Mat' is not full */
177 /* ranked. Return zero. */
178 if(j>=Mat->NbRows) {
179
180 /* Clear all the 'Value' variables */
181 value_clear(m1); value_clear(m2);
182 value_clear(m3); value_clear(gcd); value_clear(tmp);
183 return 0;
184 }
185 }
186
187 /* Now Mat[k][k] is the pivot element */
188 for(j=k+1;j<Mat->NbRows;++j) {
189
190 /* Make every other entry (below row(k)) in column(k) zero */
191 value_gcd(gcd, Mat->p[j][k], Mat->p[k][k]);
192 for(i=k+1;i<Mat->NbColumns;++i) {
193
194 /* pour tous les indices i > k */
195 value_multiply(m1,Mat->p[j][i],Mat->p[k][k]);
196 value_multiply(m2,Mat->p[j][k],Mat->p[k][i]);
197 value_subtract(m3,m1,m2);
198 value_division(Mat->p[j][i],m3,gcd);
199 }
200 }
201 }
202
203 /* Clear all the 'Value' variables */
204 value_clear(m1); value_clear(m2);
205 value_clear(m3); value_clear(gcd); value_clear(tmp);
206
207 /* The matrix 'Mat' is full ranked, return 1 */
208 return 1;
209 } /* TestRank */
210
211 /*
212 * The Saturation matrix is defined to be an integer (int type) matrix. It is
213 * a boolean matrix which has a row for every constraint and a column for
214 * every line or ray. The bits in the binary format of each integer in the
215 * saturation matrix stores the information whether the corresponding constr-
216 * aint is saturated by ray(line) or not.
217 */
218 typedef struct {
219 unsigned int NbRows;
220 unsigned int NbColumns;
221 unsigned int **p;
222 unsigned int *p_init;
223 } SatMatrix;
224
225 static SatMatrix *SMAlloc(int rows,int cols) {
226
227 unsigned int **q, *p;
228 int i;
229
230 SatMatrix *result = (SatMatrix *)malloc(sizeof(SatMatrix));
231 assert (result != NULL);
232
233 result->NbRows = rows;
234 result->NbColumns = cols;
235 result->p = q = (unsigned int **)malloc(rows * sizeof(unsigned int *));
236 assert (result->p != NULL);
237 result->p_init = p = (unsigned int *)malloc(rows * cols * sizeof(unsigned int));
238 assert (result->p_init != NULL);
239
240 for (i=0;i<rows;i++) {
241 *q++ = p;
242 p += cols;
243 }
244
245 return result;
246 } /* SMAlloc */
247
248 static void SMPrint (SatMatrix *matrix) {
249
250 unsigned int *p;
251 int i, j;
252 unsigned NbRows, NbColumns;
253
254 fprintf(stderr,"%d %d\n",NbRows=matrix->NbRows, NbColumns=matrix->NbColumns);
255 for (i=0;i<NbRows;i++) {
256 p = *(matrix->p+i);
257 for (j=0;j<NbColumns;j++)
258 fprintf(stderr, " %10X ", *p++);
259 fprintf(stderr, "\n");
260 }
261 } /* SMPrint */
262
263
264 static void SMFree (SatMatrix *matrix) {
265
266 free ((char *) matrix->p_init);
267 free ((char *) matrix->p);
268 free ((char *) matrix);
269 return;
270 } /* SMFree */
271
272 /* -------------------------------------------------------------------------
273 * Shared Global Variables:
274 * Used by procedures: Find_m_face, scan_m_face, Poly2Sat, traite_m_face,
275 * count_sat
276 * -------------------------------------------------------------------------
277 */
278 static int m; /* number of parameters */
279 static int m_dim; /* dimension of m-face */
280 static int n; /* dimension (not including parameters) */
281 static int ws; /* Working Space size */
282 static int nr; /* (NbRays-1)/32 + 1 */
283
284 static Polyhedron *CEqualities;/* Equalities in the context */
285 static SatMatrix *Sat; /* Saturation Matrix (row=constraint, col=ray)*/
286 static unsigned int *egalite; /* Bool vector marking constraints in m-face */
287 static Matrix *Xi, *Pi; /* Xi and Pi */
288 static Matrix *PiTest; /* Matrix used to test if Pi is full ranked? */
289 static Matrix *CTest;
290 static Matrix *PiInv; /* Matrix inverse Pi, with the last col of */
291 /* each line = denominator of the line */
292 static Matrix *RaysDi; /* Constraint matrix for computing Di */
293
294 static int KD; /* Flag : keep the full domains in memory ? */
295 /* 1 = yes; 0 = no, keep constraints only */
296
297 static int nbPV; /* The number of parameterized vertices */
298 static Param_Vertices *PV_Result; /* List of parameterized vertices */
299 static Param_Domain *PDomains; /* List of domains. */
300
301 #ifdef DEBUGPP
302 static int nbfaces;
303 #endif
304
305 /*
306 * Add the constraints from the context polyhedron 'CEqualities' to the
307 * constraints of polyhedra in the polyhedral domain 'D' and return the new
308 * polyhedral domain. Polyhedral domain 'D' is subsequently deleted from memory
309 */
310 static Polyhedron *Add_CEqualities(Polyhedron *D) {
311
312 Polyhedron *d,*r,*tmp;
313
314 if(!CEqualities)
315 return D;
316 else {
317 if(!D || emptyQ(D)) {
318 if(D)
319 Domain_Free(D);
320 return(Polyhedron_Copy(CEqualities));
321 }
322 r = AddConstraints(D->Constraint[0],D->NbConstraints,
323 CEqualities,ws);
324 tmp = r;
325 for(d=D->next;d;d=d->next) {
326 tmp->next = AddConstraints(d->Constraint[0],d->NbConstraints,
327 CEqualities,ws);
328 tmp = tmp->next;
329 }
330 Domain_Free(D);
331 return(r);
332 }
333 } /* Add_CEqualities */
334
335 #define INT_BITS (sizeof(unsigned) * 8)
336
337 unsigned int *int_array2bit_vector(unsigned int *array, int n)
338 {
339 int i, ix;
340 unsigned bx;
341 int words = (n+INT_BITS-1)/INT_BITS;
342 unsigned int *bv = (unsigned int *)calloc(words, sizeof(unsigned));
343 assert(bv);
344 for (i = 0, ix = 0, bx = MSB; i < n; ++i) {
345 if (array[i])
346 bv[ix] |= bx;
347 NEXT(ix, bx);
348 }
349 return bv;
350 }
351
352 /*----------------------------------------------------------------------*/
353 /* traite_m_face */
354 /* Given an m-face, compute the parameterized vertex */
355 /* D - The entire domain */
356 /* mf - Bit vector marking the lines/rays in the m-face */
357 /* egalite - boolean vector marking saturated constraints in m-face */
358 /*----------------------------------------------------------------------*/
359 static void traite_m_face(Polyhedron *D, unsigned int *mf,
360 unsigned int *egalite)
361 {
362 Matrix *Si; /* Solution */
363 Polyhedron *PDi; /* polyhedron Di */
364 Param_Vertices *PV;
365 int j,k,c,r;
366 unsigned kx, bx;
367
368 #ifdef DEBUGPP
369 ++nbfaces;
370 #endif
371
372 /* Extract Xi, Pi, and RaysDi from D */
373 RaysDi->NbRows = 0;
374 for(k=0,c=0,kx=0,bx=MSB;k<D->NbRays;++k) {
375 if(mf[kx]&bx) { /* this ray is in the current m-face */
376 if(c<m+1) {
377 int i;
378
379 /* tester si cette nouvelle colonne est lin. indep. des autres */
380 /* i.e. si gauss ne donne pas de '0' sur la colonne Pi */
381 /* jusqu'a l'indice 'c' */
382
383 /* construit PiTest */
384 for(j=0;j<m+1;++j) {
385 for(i=0;i<c;++i)
386
387 /* les c premieres colonnes */
388 value_assign(PiTest->p[j][i],Pi->p[j][i]);
389
390 /* la nouvelle */
391 value_assign(PiTest->p[j][c],D->Ray[k][j+1+n]);
392 }
393 PiTest->NbColumns = c+1;
394 r = TestRank(PiTest);
395 if(r /* TestRank(PiTest) */) {
396
397 /* Ok, c'est lin. indep. */
398 for (j=0;j<n;j++)
399 value_assign(Xi->p[j][c],D->Ray[k][j+1]); /* Xi */
400 for (j=0;j<m;j++)
401 value_assign(Pi->p[j][c],D->Ray[k][j+1+n]); /* Pi */
402 value_assign(Xi->p[n][c],D->Ray[k][n+m+1]); /* const */
403 value_assign(Pi->p[m][c],D->Ray[k][n+m+1]); /* const */
404 c++;
405 }
406 }
407
408 /* Status bit */
409 value_assign(RaysDi->p[RaysDi->NbRows][0],D->Ray[k][0]);
410 Vector_Copy(&D->Ray[k][n+1],&RaysDi->p[RaysDi->NbRows][1],(m+1));
411 ++RaysDi->NbRows;
412 }
413 NEXT(kx,bx);
414 }
415
416 #ifdef DEBUGPP41
417 fprintf(stderr, "\nRaysDi=\n");
418 Matrix_Print(stderr,P_VALUE_FMT,RaysDi);
419 if(c < m+1)
420 fprintf(stderr, "Invalid ");
421 fprintf(stderr, "Pi=\n");
422 Matrix_Print(stderr,P_VALUE_FMT,Pi);
423 #endif
424
425 #ifdef DEBUGPP4
426 if(c < m+1)
427 fprintf(stderr,"Eliminated because of no vertex\n");
428 #endif
429
430 if(c < m+1)
431 return;
432
433 /* RaysDi->numColumns = m+2; */ /* stays the same */
434
435 /* Xi->NbColumns = m+1;*/ /* VIN100: stays the same. was 'c'! */
436 /* Xi->NbRows = n+1; */ /* stays the same */
437 /* Pi->NbColumns = m+1;*/ /* VIN100: stays the same. was 'c'! */
438 /* Pi->NbRows = m+1; */ /* stays the same */
439
440 #ifdef DEBUGPP4
441 fprintf(stderr,"Xi = ");
442 Matrix_Print(stderr,P_VALUE_FMT,Xi);
443 fprintf(stderr,"Pi = ");
444 Matrix_Print(stderr,P_VALUE_FMT,Pi);
445 #endif
446
447 /* (Right) invert Pi if POSSIBLE, if not then next m-face */
448 /* Pi is destroyed */
449 if(!MatInverse(Pi,PiInv)) {
450
451 #ifdef DEBUGPP4
452 fprintf(stderr, "Eliminated because of no inverse Pi\n");
453 #endif
454
455 return;
456 }
457
458 #ifdef DEBUGPP4
459 fprintf(stderr,"FACE GENERATED!\n");
460 fprintf(stderr,"PiInv = ");
461 Matrix_Print(stderr,P_VALUE_FMT,PiInv);
462 #endif
463
464 /* Compute Si (now called Ti in the paper) */
465 Si = Matrix_Alloc(Xi->NbRows,PiInv->NbColumns);
466 rat_prodmat(Si,Xi,PiInv);
467
468 #ifdef DEBUGPP4
469 fprintf(stderr,"Si = ");
470 Matrix_Print(stderr,P_VALUE_FMT,Si);
471 #endif
472
473 Si->NbRows--; /* throw out the last row = 0 ... 0 1 */
474
475 /* Copy all of that into the PV structure */
476 PV = (Param_Vertices *) malloc(sizeof(Param_Vertices));
477 PV->next = PV_Result;
478 PV->Vertex = Si;
479 PV->Domain = NULL;
480 PV->Facets = int_array2bit_vector(egalite, D->NbConstraints);
481 PV_Result = PV;
482 nbPV++; /* increment vertex count */
483
484 /* Ok... now compute the parameter domain */
485 PDi = Rays2Polyhedron(RaysDi,ws);
486
487 #ifdef DEBUGPP3
488 fprintf(stderr,"RaysDi = ");
489 Matrix_Print(stderr,P_VALUE_FMT,RaysDi);
490 fprintf(stderr,"PDi = ");
491 Polyhedron_Print(stderr,P_VALUE_FMT,PDi);
492 #endif
493
494 if(KD==0) {
495
496 /* Add the equalities again to the domain */
497 PDi = Add_CEqualities(PDi);
498 PV->Domain = Polyhedron2Constraints(PDi);
499 Polyhedron_Free(PDi);
500 }
501 else {
502 Param_Domain *PD;
503 PD = (Param_Domain *) malloc(sizeof(Param_Domain));
504 PD->Domain = PDi;
505 PD->F = NULL;
506 PD->next = PDomains;
507 PDomains = PD;
508 }
509 return;
510 } /* traite_m_face */
511
512 /*----------------------------------------------------------------------*/
513 /* count_sat */
514 /* count the number of saturated rays in the bit vector mf */
515 /* Uses nr from global area */
516 /*----------------------------------------------------------------------*/
517 int cntbit[256] = { /* counts for 8 bits */
518 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4,
519 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
520 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
521 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
522
523 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
524 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
525 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
526 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
527
528 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
529 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
530 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
531 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
532
533 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
534 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
535 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
536 4, 5, 5, 6, 5, 6, 6, 7, 5, 6, 6, 7, 6, 7, 7, 8 };
537
538 static int count_sat (unsigned int *mf) {
539
540 register unsigned int i, tmp, cnt=0;
541
542 for (i=0; i<nr; i++) {
543 tmp = mf[i];
544 cnt = cnt
545 + cntbit[ tmp & 0xff ]
546 + cntbit[ (tmp>>8) & 0xff ]
547 + cntbit[ (tmp>>16) & 0xff ]
548 + cntbit[ (tmp>>24) & 0xff ]
549 ;
550 }
551 return cnt;
552 } /* count_sat */
553
554 /* Returns true if all bits in part are also set in bv */
555 static int bit_vector_includes(unsigned int *bv, int len, unsigned int *part)
556 {
557 int j;
558
559 for (j = 0; j < len; j++) {
560 #ifdef DEBUGPP4
561 fprintf(stderr, "mf=%08X Sat=%08X &=%08X\n",
562 part[j],bv[j], (part[j] & bv[j]));
563 #endif
564 if ((part[j] & bv[j]) != part[j])
565 return 0;
566 }
567 return 1;
568 }
569
570 /*----------------------------------------------------------------------*/
571 /* let D + E + L be the dimension of the polyhedron */
572 /* D = Dimension of polytope (ray space) */
573 /* L = Dimension of Lineality space (number of lines, bid) */
574 /* E = Dimension of Affine hull (number of equations) */
575 /* n = number of data indices */
576 /* m = number of parameters */
577 /* full domain: */
578 /* n + m = D + E + L */
579 /* projected domains: */
580 /* m = D_m + E_m + L_m */
581 /* n = D_n + E_n + L_n */
582 /* What dimension M-face, when projected onto parameter space, */
583 /* will give an L_m-face? */
584 /* What are the conditions? */
585 /* - at least one vertex */
586 /* - number of rays >= D_m+1 after removal of redundants */
587 /* */
588 /* dim of face nb saturated constraints nb saturated lines,rays */
589 /* ----------- ------------------------ ----------------------- */
590 /* (0+L)-face all E eqns + >=D ineq all L lines + 1 ray */
591 /* (M+L)-face all E eqns + >=(D-M) ineq all L lines + >=(M+1) rays */
592 /* (D+L)-face all E eqns + 0 ineq all L lines + >=(D+1) rays */
593 /*----------------------------------------------------------------------*/
594 /*----------------------------------------------------------------------*/
595 /* scan_m_face */
596 /* pos : the next candidate constraint position */
597 /* nb_un : number of saturated constraints needed to finish a face */
598 /* D : the source polyhedron (context included ) */
599 /* mf : bit-array marking rays which are saturated so far */
600 /* From Global area: */
601 /* ---------------- */
602 /* n : number of data indices */
603 /* m : number of parameters */
604 /* egalite : boolean vector marking saturated constraints in m-face */
605 /* Sat : Saturation Matrix (row=constraints, col=rays) */
606 /* ws : working space size */
607 /* nr : (NbRays-1)/32 + 1 */
608 /* */
609 /* Recursive function to find the rays and vertices of each m-face */
610 /*----------------------------------------------------------------------*/
611 static void scan_m_face(int pos,int nb_un,Polyhedron *D,unsigned int *mf) {
612 /* pos - the next candidate constraint position */
613 /* nb_un - the number of constraints needed to finish a face */
614 /* D - the source polyhedron */
615 /* mf - (bit vector) marks rays that are saturated so far */
616
617 unsigned int *new_mf;
618
619 #ifdef DEBUGPP4
620 fprintf(stderr,"Start scan_m_face(pos=%d, nb_un=%d, n=%d, m=%d\n",
621 pos,nb_un,n,m);
622 fprintf(stderr,"mf = ");
623 {
624 int i;
625 for(i=0;i<nr;i++)
626 fprintf(stderr,"%08X", mf[i]);
627 fprintf(stderr,"\nequality = [");
628 for(i=0;i<D->NbConstraints;i++)
629 fprintf(stderr," %1d",egalite[i]);
630 fprintf(stderr,"]\n");
631 }
632 #endif
633
634 if(nb_un == 0) { /* Base case */
635 int i;
636
637 /*********** ELIMINATION OF REDUNDANT FACES ***********/
638 /* if all these vertices also verify a previous constraint */
639 /* which is NOT already selected, we eliminate this face */
640 /* This keeps the lexicographically greatest selection */
641 for(i=0;i<pos-1;i++)
642 {
643 if(egalite[i])
644 continue; /* already selected */
645
646 /* if Sat[i] & mf == mf then it's redundant */
647 #ifdef DEBUGPP4
648 fprintf(stderr, "Sat[%d]\n", i);
649 #endif
650 if (bit_vector_includes(Sat->p[i], nr, mf)) {
651 #ifdef DEBUGPP4
652 fprintf(stderr, "Redundant with constraint %d\n", i);
653 #endif
654 return; /* it is redundant */
655 }
656 }
657 /********* END OF ELIMINATION OF DEGENERATE FACES *********/
658 /* Now check for other constraints that are verified */
659 for (i = pos; i < D->NbConstraints; ++i) {
660 if (bit_vector_includes(Sat->p[i], nr, mf))
661 egalite[i] = 1;
662 }
663 /* if we haven't found a constraint verified by all */
664 /* the rays, its OK, it's a new face. */
665 traite_m_face(D, mf, egalite);
666 for (i = pos; i < D->NbConstraints; ++i)
667 egalite[i] = 0;
668 return;
669 }
670
671 /* See if there are enough constraints left to finish */
672 if((pos+nb_un)>D->NbConstraints) return;
673
674 /* Recurring part of the procedure */
675 /* Add the pos'th constraint, compute new saturation vector */
676 {
677 int k;
678 new_mf = (unsigned int *)malloc(nr*sizeof(unsigned int));
679 for (k=0; k<nr; k++)
680 new_mf[k] = mf[k] & Sat->p[pos][k];
681 }
682 #ifdef DEBUGPP4
683 fprintf(stderr,"new_mf = ");
684 {
685 int i;
686 for(i=0;i<nr;i++) {
687 fprintf(stderr,"%08X", new_mf[i]);
688 }
689 fprintf(stderr,"\ncount(new_mf) = %d\n",count_sat(new_mf));
690 }
691 #endif
692
693 {
694 int c;
695 c = count_sat(new_mf);
696 /* optimization : at least m_dim+1 rays must be saturated to add this constraint */
697 if (c>m_dim )
698 {
699 int redundant = 0;
700
701 egalite[pos]=1; /* Try it with the pos-th constraint */
702
703 /* If this constraint does not change anything,
704 * it is redundant with respect to the selected
705 * equalities and the remaining inequalities.
706 * Check whether it is redundant with respect
707 * to just the selected equalities.
708 */
709 if( c==count_sat(mf) ) {
710 int i, c, j;
711
712 for (i = 0, c = 0; i < D->NbConstraints; ++i) {
713 if (egalite[i] == 0 || egalite[i] == -1)
714 continue;
715 for (j = 0; j < D->Dimension+1; ++j)
716 value_assign(CTest->p[j][c],
717 D->Constraint[i][j+1]);
718 ++c;
719 }
720 CTest->NbColumns = c;
721 #ifdef DEBUGPP41
722 Matrix_Print(stderr,P_VALUE_FMT,CTest);
723 #endif
724 redundant = !TestRank(CTest);
725 }
726
727 /* Do not decrement nb_un if equality is redundant. */
728 if( redundant )
729 {
730 egalite[pos]=-1; /* Don't use in further redundance test
731 */
732 scan_m_face(pos+1,nb_un,D,new_mf);
733 }
734 else
735 {
736 scan_m_face(pos+1,nb_un-1,D,new_mf);
737 }
738 }
739 }
740 free(new_mf);
741 egalite[pos]=0; /* Try it without the pos-th constraint */
742 if ((pos+nb_un)>=D->NbConstraints) return;
743 scan_m_face(pos+1,nb_un,D,mf);
744 return;
745 } /* scan_m_face */
746
747 /*
748 * Create a saturation matrix with rows correspond to the constraints and
749 * columns correspond to the rays of the polyhedron 'Pol'. Global variable
750 * 'nr' is set in the function.
751 */
752 static SatMatrix *Poly2Sat(Polyhedron *Pol,unsigned int **L) {
753
754 SatMatrix *Sat;
755 int i, j, k, kx;
756 unsigned int *Temp;
757 Value *p1, *p2, p3,tmp;
758 unsigned Dimension, NbRay, NbCon, bx;
759
760 /* Initialize all the 'Value' variables */
761 value_init(p3); value_init(tmp);
762
763 NbRay = Pol->NbRays;
764 NbCon = Pol->NbConstraints;
765 Dimension = Pol->Dimension+1; /* Homogeneous Dimension */
766
767 /* Build the Sat matrix */
768 nr = (NbRay - 1)/(sizeof(int)*8) + 1; /* Set globally */
769 Sat = SMAlloc(NbCon,nr);
770 Temp = (unsigned int *)malloc(nr*sizeof(unsigned int));
771 memset(Sat->p_init,0,nr*NbCon*sizeof(int));
772 memset(Temp,0,nr*sizeof(unsigned int));
773 kx=0; bx=MSB;
774 for (k=0; k<NbRay; k++) {
775 for (i=0; i<NbCon; i++) {
776 p1 = &Pol->Constraint[i][1];
777 p2 = &Pol->Ray[k][1];
778 value_set_si(p3,0);
779 for (j=0;j<Dimension;j++) {
780 value_multiply(tmp,*p1,*p2);
781 value_addto(p3,p3,tmp);
782 p1++; p2++;
783 }
784 if (value_zero_p(p3))
785 Sat->p[i][kx]|=bx;
786 }
787 Temp[kx] |= bx;
788 NEXT(kx, bx);
789 }
790
791 /* Set 'L' to an array containing ones in every bit position of its */
792 /* elements. */
793 *L = Temp;
794
795 /* Clear all the 'Value' variables */
796 value_clear(p3); value_clear(tmp);
797
798 return Sat;
799 } /* Poly2Sat */
800
801 /*
802 * Create a parametrized polyhedron with zero parameters. This function was
803 * first written by Xavier Redon, and was later modified by others.
804 */
805 Param_Polyhedron *GenParamPolyhedron(Polyhedron *Pol, Matrix *Rays)
806 {
807 Param_Polyhedron *result;
808 int nbRows, nbColumns;
809 int i, size, rays;
810
811 nbRows=Pol->NbRays;
812 nbColumns=Pol->Dimension+2;
813
814 /* Count the number of rays */
815 for(i=0, rays=0; i<nbRows; i++)
816 if(value_notone_p(Pol->Ray[i][0]) ||
817 value_zero_p(Pol->Ray[i][nbColumns-1]))
818 ++rays;
819
820 /* Initialize the result */
821 result=(Param_Polyhedron *)malloc(sizeof(Param_Polyhedron));
822 result->nbV=nbRows-rays;
823 result->V=NULL;
824 result->Constraints = Polyhedron2Constraints(Pol);
825 result->Rays = Rays;
826
827 /* Build the parametric vertices */
828 for(i=0;i<nbRows;i++) {
829 Matrix *vertex;
830 Param_Vertices *paramVertex;
831 int j;
832
833 if (value_notone_p(Pol->Ray[i][0]) ||
834 value_zero_p(Pol->Ray[i][nbColumns-1]))
835 continue;
836
837 vertex=Matrix_Alloc(nbColumns-2,2);
838 for(j=1;j<nbColumns-1;j++) {
839 value_assign(vertex->p[j-1][0],Pol->Ray[i][j]);
840 value_assign(vertex->p[j-1][1],Pol->Ray[i][nbColumns-1]);
841 }
842 paramVertex=(Param_Vertices *)malloc(sizeof(Param_Vertices));
843 paramVertex->Vertex=vertex;
844
845 /* There is one validity domain : universe of dimension 0 */
846 paramVertex->Domain=Matrix_Alloc(1,2);
847 value_set_si(paramVertex->Domain->p[0][0],1);
848 value_set_si(paramVertex->Domain->p[0][1],1);
849 paramVertex->Facets = NULL;
850 paramVertex->next=result->V;
851 result->V=paramVertex;
852 }
853
854 /* Build the parametric domains (only one here) */
855 if (nbRows > 1)
856 size=(nbRows-1)/(8*sizeof(int))+1;
857 else
858 size = 1;
859 result->D=(Param_Domain *)malloc(sizeof(Param_Domain));
860 result->D->next=NULL;
861 result->D->Domain=Universe_Polyhedron(0);
862 result->D->F=(unsigned int *)malloc(size*sizeof(int));
863 memset(&result->D->F[0],0xFF,size*sizeof(int));
864
865 return result;
866 } /* GenParamPolyhedron */
867
868
869 /*----------------------------------------------------------------------*/
870 /* PreElim_Columns */
871 /* function being called before Elim_Columns */
872 /* Equalities in E are analysed to initialize ref and p. */
873 /* These two vectors are used to construct the new constraint matrix */
874 /* PreElim_Columns returns the transformation matrix to re-convert the */
875 /* resulting domains in the same format (by adding empty columns) */
876 /* in the parameter space */
877 /*----------------------------------------------------------------------*/
878 Matrix *PreElim_Columns(Polyhedron *E,int *p,int *ref,int m) {
879
880 int i,j,l;
881 Matrix *T;
882
883 /* find which columns to eliminate */
884 /* p contains, for each line in E, the column to eliminate */
885 /* (i.e. the corresponding parameter index to eliminate) */
886 /* 0 <= i <= E->NbEq, and 1 <= p[i] <= E->Dimension */
887 memset(p,0,sizeof(int)*(E->NbEq));
888
889 #ifdef DEBUGPP32
890 fprintf(stderr,"\n\nPreElim_Columns starting\n");
891 fprintf(stderr,"E =\n");
892 Polyhedron_Print(stderr,P_VALUE_FMT,E);
893 #endif
894
895 for(l=0;l<E->NbEq;++l) {
896 if(value_notzero_p(E->Constraint[l][0])) {
897 fprintf(stderr,"Internal error: Elim_Columns (polyparam.c), expecting equalities first in E.\n");
898 free(p);
899 return(NULL);
900 }
901 for(i=1;value_zero_p(E->Constraint[l][i]);++i) {
902 if(i==E->Dimension+1) {
903 fprintf(stderr,"Internal error: Elim_Columns (polyparam.c), expecting non-empty constraint in E.\n");
904 free(p);
905 return( NULL );
906 }
907 }
908 p[l] = i;
909
910 #ifdef DEBUGPP32
911 fprintf(stderr,"p[%d] = %d,",l,p[l]);
912 #endif
913 }
914
915 /* Reference vector: column ref[i] in A corresponds to column i in M */
916 for(i=0;i<E->Dimension+2-E->NbEq;++i) {
917 ref[i]=i;
918 for(j=0;j<E->NbEq;++j)
919 if(p[j]<=ref[i])
920 ref[i]++;
921
922 #ifdef DEBUGPP32
923 fprintf(stderr,"ref[%d] = %d,",i,ref[i]);
924 #endif
925 }
926
927 /* Size of T : phdim-nbEq * phdim */
928 T = Matrix_Alloc(m+1-E->NbEq,m+1);
929 for(i=0;i<T->NbColumns;i++)
930 for(j=0;j<T->NbRows;j++)
931 if(ref[E->Dimension-m+j+1] == E->Dimension-m+i+1)
932 value_set_si(T->p[j][i],1);
933 else
934 value_set_si(T->p[j][i],0);
935
936 #ifdef DEBUGPP32
937 fprintf(stderr,"\nTransMatrix =\n");
938 Matrix_Print(stderr,P_VALUE_FMT,T);
939 #endif
940
941 return(T);
942
943 } /* PreElim_Columns */
944
945 /*----------------------------------------------------------------------*/
946 /* Elim_Columns */
947 /* Eliminate columns from A, using the equalities in E. */
948 /* ref and p are precalculated by PreElim_Columns, using E; */
949 /* these two vectors are used to construct the new constraint matrix */
950 /*----------------------------------------------------------------------*/
951 Polyhedron *Elim_Columns(Polyhedron *A,Polyhedron *E,int *p,int *ref) {
952
953 int i,l,c;
954 Matrix *M, *C;
955 Polyhedron *R;
956 Value tmp1,tmp2;
957
958 /* Initialize all the 'Value' variables */
959 value_init(tmp1); value_init(tmp2);
960
961 #ifdef DEBUGPP32
962 fprintf(stderr,"\nElim_Columns starting\n");
963 fprintf(stderr,"A =\n" );
964 Polyhedron_Print(stderr,P_VALUE_FMT,A);
965 #endif
966
967 /* Builds M = constraint matrix of A, useless columns zeroed */
968 M = Polyhedron2Constraints(A);
969 for(l=0;l<E->NbEq;++l) {
970 for(c=0;c<M->NbRows;++c) {
971 if(value_notzero_p(M->p[c][p[l]])) {
972
973 /* A parameter to eliminate here ! */
974 for(i=1;i<M->NbColumns;++i) {
975 value_multiply(tmp1,M->p[c][i],E->Constraint[l][p[l]]);
976 value_multiply(tmp2,E->Constraint[l][i],M->p[c][p[l]]);
977 value_subtract(M->p[c][i],tmp1,tmp2);
978 }
979 }
980 }
981 }
982
983 #ifdef DEBUGPP32
984 fprintf(stderr,"\nElim_Columns after zeroing columns of A.\n");
985 fprintf(stderr,"M =\n");
986 Matrix_Print(stderr,P_VALUE_FMT,M);
987 #endif
988
989 /* Builds C = constraint matrix, useless columns eliminated */
990 C = Matrix_Alloc(M->NbRows,M->NbColumns - E->NbEq);
991 for(l=0;l<C->NbRows;++l)
992 for(c=0;c<C->NbColumns;++c) {
993 value_assign(C->p[l][c],M->p[l][ref[c]]);
994 }
995
996 #ifdef DEBUGPP32
997 fprintf(stderr,"\nElim_Columns after eliminating columns of A.\n");
998 fprintf(stderr,"C =\n");
999 Matrix_Print(stderr,P_VALUE_FMT,C);
1000 #endif
1001
1002 R = Constraints2Polyhedron(C,ws);
1003 Matrix_Free(C);
1004 Matrix_Free(M);
1005 value_clear(tmp1); value_clear(tmp2);
1006 return(R);
1007 } /* Elim_Columns */
1008
1009
1010 static Polyhedron *Recession_Cone(Polyhedron *P, unsigned nvar, unsigned MaxRays)
1011 {
1012 int i;
1013 Matrix *M = Matrix_Alloc(P->NbConstraints, 1 + nvar + 1);
1014 Polyhedron *R;
1015 for (i = 0; i < P->NbConstraints; ++i)
1016 Vector_Copy(P->Constraint[i], M->p[i], 1+nvar);
1017 R = Constraints2Polyhedron(M, MaxRays);
1018 Matrix_Free(M);
1019 return R;
1020 }
1021
1022 /* Compute lines/unidirectional rays of the (non parametric) polyhedron */
1023 /* Input :
1024 * D1 : combined polyhedron,
1025 * Output :
1026 * Rays : non parametric ray matrix
1027 * return value : number of lines
1028 */
1029 static int ComputeNPLinesRays(int n, Polyhedron *D1, Matrix **Rays)
1030 {
1031 int i, j, nbr, dimfaces;
1032 Polyhedron *RC; /* Recession Cone */
1033
1034 RC = Recession_Cone(D1, n, ws);
1035
1036 /* get the rays/lines from RC */
1037 nbr = 0;
1038 for (i = 0; i < RC->NbRays ;i++)
1039 if (value_zero_p(RC->Ray[i][n+1]))
1040 nbr++;
1041 *Rays=Matrix_Alloc(nbr, n+2);
1042 for (i = 0, j = 0; j < nbr ;i++)
1043 if (value_zero_p(RC->Ray[i][n+1]))
1044 Vector_Copy(RC->Ray[i], (*Rays)->p[j++], n+2);
1045
1046 dimfaces = RC->NbBid;
1047 Polyhedron_Free(RC);
1048
1049 #ifdef DEBUGPP31
1050 fprintf(stderr, "Rays = ");
1051 Matrix_Print(stderr, P_VALUE_FMT, *Rays);
1052 fprintf(stderr, "dimfaces = %d\n", dimfaces);
1053 #endif
1054
1055 return dimfaces;
1056 }
1057
1058 /*
1059 * Given a polyhedron 'Di' in combined data and parameter space and a context
1060 * polyhedron 'C' representing the constraints on the parameter space, create
1061 * a list of parameterized vertices and assign values to global variables:
1062 * m,n,ws,Sat,egalite,mf,Xi,Pi,PiInv,RaysDi,CEqualities.
1063 */
1064 Param_Polyhedron *Find_m_faces(Polyhedron **Di,Polyhedron *C,int keep_dom,int working_space,Polyhedron **CEq,Matrix **CT) {
1065
1066 unsigned int *mf;
1067 int i, j, dimfaces;
1068 Polyhedron *D=*Di,
1069 *C1, /* true context */
1070 *D1; /* the combined polyhedron, including context C */
1071 Matrix *Rays, /* lines/rays (non parametric) */
1072 *M;
1073 Param_Polyhedron *res;
1074 int *p, *ref;
1075
1076 CEqualities = NULL;
1077
1078 if(CT) {
1079 *CEq = NULL;
1080 *CT = NULL;
1081 }
1082
1083 if(!D || !C)
1084 return (Param_Polyhedron *) 0;
1085
1086 ws = working_space;
1087 m = C->Dimension;
1088 n = D->Dimension - m;
1089 if(n<0) {
1090 fprintf(stderr,
1091 "Find_m_faces: ?%d parameters of a %d-polyhedron !\n",m,n);
1092 return (Param_Polyhedron *) 0;
1093 }
1094 if (m==0)
1095 return GenParamPolyhedron(D,Matrix_Alloc(0,2));
1096
1097 /* Add constraints from Context to D -> result in D1 */
1098 C1 = align_context(C,D->Dimension,ws);
1099
1100 #ifdef DEBUGPP31
1101 fprintf(stderr,"m = %d\n",m);
1102 fprintf(stderr, "D = ");
1103 Polyhedron_Print(stderr,P_VALUE_FMT,D);
1104 fprintf(stderr,"C1 = ");
1105 Polyhedron_Print(stderr,P_VALUE_FMT,C1);
1106 #endif
1107
1108 D1 = DomainIntersection(D,C1,ws);
1109
1110 #ifdef DEBUGPP31
1111 fprintf(stderr,"D1 = ");
1112 Polyhedron_Print(stderr,P_VALUE_FMT,D1);
1113 #endif
1114
1115 Domain_Free(C1);
1116
1117 if (!D1)
1118 return(NULL);
1119 if (emptyQ(D1)) {
1120 Polyhedron_Free(D1);
1121 return(NULL);
1122 }
1123
1124 /* Compute the true context C1 */
1125 /* M : lines in the direction of the first n indices (index space) */
1126 M = Matrix_Alloc(n, D1->Dimension+2);
1127 for (i=0; i<n; i++)
1128 value_set_si(M->p[i][i+1],1);
1129 C1 = DomainAddRays(D1,M,ws);
1130 Matrix_Free(M);
1131
1132 #ifdef DEBUGPP31
1133 fprintf(stderr,"True context C1 = ");
1134 Polyhedron_Print(stderr,P_VALUE_FMT,C1);
1135 #endif
1136
1137 dimfaces = ComputeNPLinesRays(n, D1, &Rays);
1138
1139 /* CEqualities contains the constraints (to be added again later) */
1140 /* *CT is the transformation matrix to add the removed parameters */
1141 if(!CT) {
1142 if (C1->NbEq == 0) {
1143 Polyhedron_Free(C1);
1144 C1 = NULL;
1145 } else {
1146 Polyhedron *CEq1, /* CEqualities, in homogeneous dim */
1147 *D2; /* D1, (temporary) simplified */
1148
1149 /* Remove equalities from true context C1 and from D1 */
1150 /* Compute CEqualities = matrix of equalities in C1, projected in */
1151 /* the parameter space */
1152 M = Matrix_Alloc(C1->NbEq,m+2);
1153 for(j=0,i=0;i<C1->NbEq;++i,++j) {
1154 while(value_notzero_p(C1->Constraint[j][0]))
1155 ++j;
1156 value_assign(M->p[i][0],C1->Constraint[j][0]);
1157 Vector_Copy(&C1->Constraint[j][D->Dimension-m+1],&M->p[i][1],(m+1));
1158 }
1159 CEqualities = Constraints2Polyhedron(M,ws);
1160 Matrix_Free(M);
1161 CEq1 = align_context(CEqualities,D->Dimension,ws);
1162
1163 /* Simplify D1 and C1 (remove the equalities) */
1164 D2 = DomainSimplify(D1,CEq1,ws);
1165 Polyhedron_Free(D1);
1166 Polyhedron_Free(C1);
1167 Polyhedron_Free(CEq1);
1168 D1 = D2;
1169 C1 = NULL;
1170 }
1171 }
1172 else { /* if( CT ) */
1173 Polyhedron *CEq1, /* CEqualities */
1174 *D2; /* D1, (temporary) simplified */
1175
1176 /* Suppress all useless constraints in parameter domain */
1177 /* when CT is not NULL (ehrhart) */
1178 /* Vin100, march 01 */
1179 CEq1 = C1;
1180 M = Matrix_Alloc(C1->NbConstraints,m+2);
1181 for(i=0;i<C1->NbConstraints;++i) {
1182 value_assign(M->p[i][0],C1->Constraint[i][0]);
1183 Vector_Copy(&C1->Constraint[i][D->Dimension-m+1],&M->p[i][1],(m+1));
1184 }
1185 CEqualities = Constraints2Polyhedron( M, ws );
1186 Matrix_Free(M);
1187
1188 D2 = DomainSimplify(D1,CEq1,ws);
1189 Polyhedron_Free(D1);
1190 D1 = D2;
1191 C1 = Universe_Polyhedron(D2->Dimension);
1192
1193 /* if CT is not NULL, the constraints are eliminated */
1194 /* *CT will contain the transformation matrix to come back to the */
1195 /* original dimension (for a polyhedron, in the parameter space) */
1196 if( CEq1->NbEq )
1197 {
1198 m -= CEq1->NbEq;
1199 p = (int *)malloc(sizeof(int)*(CEq1->NbEq));
1200 }
1201 else
1202 p = NULL;
1203 ref = (int*) malloc(sizeof(int)*
1204 (CEq1->Dimension+2-CEq1->NbEq));
1205 *CT = PreElim_Columns(CEq1,p,ref,CEqualities->Dimension);
1206 D2 = Elim_Columns(D1,CEq1,p,ref);
1207 if (p)
1208 free(p);
1209 free(ref);
1210
1211 #ifdef DEBUGPP3
1212 fprintf(stderr,"D2\t Dim = %3d\tNbEq = %3d\tLines = %3d\n",
1213 D2->Dimension,D2->NbEq,D2->NbBid);
1214 C2 = Elim_Columns(C1,CEq1,p,ref);
1215 fprintf(stderr,"C2\t Dim = %3d\tNbEq = %3d\tLines = %3d\n",
1216 C2->Dimension,C2->NbEq,C2->NbBid);
1217 Polyhedron_Free(C2);
1218 #endif
1219
1220 Polyhedron_Free(D1);
1221 Polyhedron_Free(C1);
1222 D1 = D2;
1223 C1 = NULL;
1224 *CEq = CEqualities;
1225
1226 #ifdef DEBUGPP3
1227 fprintf(stderr,"Polyhedron CEq = ");
1228 Polyhedron_Print(stderr,P_VALUE_FMT,*CEq);
1229 fprintf(stderr,"Matrix CT = ");
1230 Matrix_Print(stderr,P_VALUE_FMT,*CT);
1231 #endif
1232
1233 Polyhedron_Free(CEq1);
1234 CEqualities = NULL; /* don't simplify ! */
1235
1236 /* m changed !!! */
1237 if(m==0) {
1238 /* return the new D1 too */
1239 *Di = D1;
1240
1241 return GenParamPolyhedron(D1, Rays);
1242 }
1243 }
1244
1245 #ifdef DEBUGPP3
1246 fprintf(stderr,"Polyhedron D1 (D AND C) = ");
1247 Polyhedron_Print(stderr,P_VALUE_FMT, D1);
1248 fprintf(stderr,"Polyhedron CEqualities = ");
1249 if(CEqualities) Polyhedron_Print(stderr,P_VALUE_FMT, CEqualities);
1250 else fprintf(stderr,"NULL\n");
1251 #endif
1252
1253 KD = keep_dom;
1254 PDomains = NULL;
1255 PV_Result = NULL;
1256 nbPV = 0;
1257
1258 if (emptyQ(D1)) {
1259 Polyhedron_Free(D1);
1260 Matrix_Free(Rays);
1261 return NULL;
1262 }
1263
1264 /* mf : a bit array indicating which rays are part of the m-face */
1265 /* Poly2Sat initializes mf to all ones */
1266 /* set global variable nr to size (number of words) of mf */
1267 Sat = Poly2Sat(D1,&mf);
1268
1269 #ifdef DEBUGPP4
1270 fprintf(stderr,"Sat = ");
1271 SMPrint(Sat);
1272
1273 fprintf(stderr,"mf = ");
1274 for (i=0; i<nr; i++)
1275 fprintf(stderr,"%08X", mf[i]);
1276 fprintf(stderr, "\n");
1277 #endif
1278
1279 /* A boolean array saying which constraints are part of the m-face */
1280 egalite = (unsigned int *)malloc(sizeof(int)*D1->NbConstraints);
1281 memset(egalite,0, sizeof(int)*D1->NbConstraints);
1282
1283 for (i=0; i<D1->NbEq; i++)
1284 egalite[i] = 1;
1285
1286 Xi = Matrix_Alloc(n+1,m+1);
1287 Pi = Matrix_Alloc(m+1,m+1);
1288 PiTest = Matrix_Alloc(m+1,m+1);
1289 CTest = Matrix_Alloc(D->Dimension+1,D->NbConstraints);
1290 PiInv = Matrix_Alloc(m+1,m+2);
1291 RaysDi = Matrix_Alloc(D1->NbRays,m+2);
1292 m_dim = m;
1293
1294 /* m_dim has to be increased by the dimension of the smallest faces
1295 * of the (non parametric) polyhedron
1296 */
1297 m_dim += dimfaces;
1298
1299 /* if the smallest face is of smaller dimension than m_dim,
1300 * then increase m_dim (I think this should never happen --Vincent)
1301 */
1302 #ifdef DEBUGPP3
1303 if (m_dim < D1->NbBid)
1304 fprintf(stderr, "m_dim (%d) < D1->NbBid (%d)\n", m_dim, D1->NbBid );
1305 #endif
1306 if (m_dim < D1->NbBid)
1307 m_dim = D1->NbBid;
1308
1309 #ifdef DEBUGPP
1310 nbfaces=0;
1311 #endif
1312 #ifdef DEBUGPP3
1313 fprintf(stderr, "m_dim = %d\n", m_dim);
1314 fprintf(stderr,
1315 "Target: find faces that saturate %d constraints and %d rays/lines\n",
1316 D1->Dimension - m_dim,m_dim+1);
1317 #endif
1318
1319 /* D1->NbEq constraints already saturated ! */
1320 scan_m_face(D1->NbEq,(D1->Dimension - m_dim - D1->NbEq),D1,mf);
1321
1322 /* pos, number of constraints needed */
1323
1324 #ifdef DEBUGPP
1325 fprintf( stderr, "Number of m-faces: %d\n", nbfaces );
1326 #endif
1327
1328 Matrix_Free(RaysDi);
1329 Matrix_Free(PiInv);
1330 Matrix_Free(PiTest);
1331 Matrix_Free(CTest);
1332 Matrix_Free(Pi);
1333 Matrix_Free(Xi);
1334 free(egalite);
1335 free(mf);
1336 SMFree(Sat);
1337
1338 /* if(CEqualities && keep_dom==0) {
1339 Domain_Free(CEqualities);
1340 }
1341 */
1342
1343 res = (Param_Polyhedron *) malloc (sizeof(Param_Polyhedron));
1344 res->nbV = nbPV;
1345 res->V = PV_Result;
1346 res->D = PDomains;
1347 res->Constraints = Polyhedron2Constraints(D1);
1348 res->Rays = Rays;
1349
1350 if(CT) /* return the new D1 too ! */
1351 *Di = D1;
1352 else
1353 Domain_Free(D1);
1354
1355 return(res);
1356 } /* Find_m_faces */
1357
1358 /*
1359 * Given parametric domain 'PD' and number of parametric vertices 'nb_domains',
1360 * find the vertices that belong to distinct sub-domains.
1361 */
1362 void Compute_PDomains(Param_Domain *PD,int nb_domains,int working_space) {
1363
1364 unsigned bx;
1365 int i, ix, nv;
1366 Polyhedron *dx, *d1, *d2;
1367 Param_Domain *p1, *p2, *p2prev, *PDNew;
1368
1369 if (nb_domains==0) {
1370
1371 #ifdef DEBUGPP5
1372 fprintf(stderr,"No domains\n");
1373 #endif
1374
1375 return;
1376 }
1377
1378 /* Already filled out by GenParamPolyhedron */
1379 if (!PD->next && PD->F)
1380 return;
1381
1382 /* Initialization */
1383 nv = (nb_domains - 1)/(8*sizeof(int)) + 1;
1384
1385 #ifdef DEBUGPP5
1386 fprintf(stderr,"nv = %d\n",nv);
1387 #endif
1388
1389 for(p1=PD,i=0,ix=0,bx=MSB;p1;p1=p1->next,i++) {
1390
1391 /* Assign a bit array 'p1->F' of suitable size to include the vertices */
1392 p1->F = (unsigned *) malloc (nv * sizeof(unsigned));
1393
1394 /* Set the bit array to zeros */
1395 memset(p1->F,0,nv * sizeof(unsigned));
1396 p1->F[ix] |= bx; /* Set i'th bit to one */
1397 NEXT(ix, bx);
1398 }
1399
1400 #ifdef DEBUGPP5
1401 fprintf(stderr,"nb of vertices=%d\n",i);
1402 #endif
1403
1404 /* Walk the PD list with two pointers */
1405 ix = 0; bx=MSB;
1406 for (p1=PD;p1;p1=p1->next) {
1407 for (p2prev=p1,p2=p1->next;p2;p2prev=p2,p2=p2->next) {
1408
1409 /* Find intersection */
1410 dx = PDomainIntersection(p1->Domain,p2->Domain,working_space);
1411
1412 if (!dx || emptyQ(dx)) {
1413
1414 #ifdef DEBUGPP5
1415 fprintf( stderr, "Empty dx (p1 inter p2). Continuing\n");
1416 #endif
1417 if(dx)
1418 Domain_Free(dx);
1419 continue;
1420 }
1421
1422 #ifdef DEBUGPP5
1423 fprintf(stderr,"Begin PDomainDifference\n");
1424 fprintf(stderr, "p1=");
1425 Polyhedron_Print(stderr,P_VALUE_FMT,p1->Domain);
1426 fprintf(stderr,"p2=");
1427 Polyhedron_Print(stderr,P_VALUE_FMT,p2->Domain);
1428 #endif
1429 d1 = PDomainDifference(p1->Domain,p2->Domain,working_space);
1430 d2 = PDomainDifference(p2->Domain,p1->Domain,working_space);
1431
1432 #ifdef DEBUGPP5
1433 fprintf(stderr,"p1\\p2=");
1434 Polyhedron_Print(stderr,P_VALUE_FMT,d1);
1435 fprintf(stderr,"p2\\p1=");
1436 Polyhedron_Print(stderr,P_VALUE_FMT,d2);
1437 fprintf(stderr,"END PDomainDifference\n\n");
1438 #endif
1439 if (!d1 || emptyQ(d1) || d1->NbEq!=0) {
1440
1441 #ifdef DEBUGPP5
1442 fprintf(stderr,"Empty d1\n");
1443 #endif
1444 if (d1)
1445 Domain_Free(d1);
1446 Domain_Free(dx);
1447
1448 if (!d2 || emptyQ(d2) || d2->NbEq!=0) {
1449
1450 #ifdef DEBUGPP5
1451 fprintf( stderr, "Empty d2 (deleting)\n");
1452 #endif
1453 /* dx = p1->Domain = p2->Domain */
1454 if (d2) Domain_Free(d2);
1455
1456 /* Update p1 */
1457 for (i=0;i<nv;i++)
1458 p1->F[i] |= p2->F[i];
1459
1460 /* Delete p2 */
1461 p2prev->next = p2->next;
1462 Domain_Free(p2->Domain);
1463 free(p2->F);
1464 free(p2);
1465 p2 = p2prev;
1466 }
1467 else { /* d2 is not empty --> dx==p1->domain */
1468
1469 #ifdef DEBUGPP5
1470 fprintf( stderr, "p2 replaced by d2\n");
1471 #endif
1472 /* Update p1 */
1473 for(i=0;i<nv;i++)
1474 p1->F[i] |= p2->F[i];
1475
1476 /* Replace p2 with d2 */
1477 Domain_Free( p2->Domain );
1478 p2->Domain = d2;
1479 }
1480 }
1481 else { /* d1 is not empty */
1482 if (!d2 || emptyQ(d2) || d2->NbEq!=0) {
1483
1484 #ifdef DEBUGPP5
1485 fprintf( stderr, "p1 replaced by d1\n");
1486 #endif
1487 if (d2) Domain_Free(d2);
1488
1489 /* dx = p2->domain */
1490 Domain_Free(dx);
1491
1492 /* Update p2 */
1493 for(i=0;i<nv;i++)
1494 p2->F[i] |= p1->F[i];
1495
1496 /* Replace p1 with d1 */
1497 Domain_Free(p1->Domain);
1498 p1->Domain = d1;
1499 }
1500 else { /*d2 is not empty-->d1,d2,dx are distinct */
1501
1502 #ifdef DEBUGPP5
1503 fprintf(stderr,"Non-empty d1 and d2\nNew node created\n");
1504 #endif
1505 /* Create a new node for dx */
1506 PDNew = (Param_Domain *) malloc( sizeof(Param_Domain) );
1507 PDNew->F = (unsigned int *)malloc( nv*sizeof(int) );
1508 memset(PDNew->F,0,nv*sizeof(int));
1509 PDNew->Domain = dx;
1510
1511 for (i=0;i<nv;i++)
1512 PDNew->F[i] = p1->F[i] | p2->F[i];
1513
1514 /* Replace p1 with d1 */
1515 Domain_Free( p1->Domain );
1516 p1->Domain = d1;
1517
1518 /* Replace p2 with d2 */
1519 Domain_Free( p2->Domain );
1520 p2->Domain = d2;
1521
1522 /* Insert new node after p1 */
1523 PDNew->next = p1->next;
1524 p1->next = PDNew;
1525 }
1526 }
1527 } /* end of p2 scan */
1528 if (p1->Domain->next) {
1529 Polyhedron *C = DomainConvex(p1->Domain, working_space);
1530 Domain_Free(p1->Domain);
1531 p1->Domain = C;
1532 }
1533 } /* end of p1 scan */
1534 } /* Compute_PDomains */
1535
1536 /*
1537 * Given a polyhedron 'Din' in combined data and parametre space, a context
1538 * polyhedron 'Cin' representing the constraints on the parameter space and
1539 * a working space size 'working_space', return a parametric polyhedron with
1540 * a list of parametric vertices and their defining domains.
1541 */
1542 Param_Polyhedron *Polyhedron2Param_Vertices(Polyhedron *Din,Polyhedron *Cin,int working_space) {
1543
1544 Param_Polyhedron *result;
1545
1546 POL_ENSURE_FACETS(Din);
1547 POL_ENSURE_VERTICES(Din);
1548 POL_ENSURE_FACETS(Cin);
1549 POL_ENSURE_VERTICES(Cin);
1550
1551 #ifdef DEBUGPP
1552 fprintf(stderr,"Polyhedron2Param_Vertices algorithm starting at : %.2fs\n",
1553 (float)clock()/CLOCKS_PER_SEC);
1554 #endif
1555
1556 /***************** Scan the m-faces ****************/
1557 result = Find_m_faces(&Din,Cin,0,working_space,NULL,NULL);
1558
1559 #ifdef DEBUGPP
1560 fprintf(stderr, "nb of points : %d\n",result->nbV);
1561 #endif
1562
1563 #ifdef DEBUGPP
1564 fprintf(stderr, "end main loop : %.2fs\n", (float)clock()/CLOCKS_PER_SEC);
1565 #endif
1566
1567 return(result);
1568 } /* Polyhedron2Param_Vertices */
1569
1570 /*
1571 * Free the memory allocated to a list of parametrized vertices
1572 */
1573 void Param_Vertices_Free(Param_Vertices *PV) {
1574
1575 Param_Vertices *next_pv;
1576
1577 while(PV) {
1578 next_pv = PV->next;
1579 if (PV->Vertex) Matrix_Free(PV->Vertex);
1580 if (PV->Domain) Matrix_Free(PV->Domain);
1581 if (PV->Facets) free(PV->Facets);
1582 free(PV);
1583 PV = next_pv;
1584 }
1585 } /* Param_Vertices_Free */
1586
1587 /*
1588 * Print a list of parametrized vertices *
1589 */
1590 void Print_Vertex(FILE *DST, Matrix *V, const char **param_names)
1591 {
1592 int l, v;
1593 int first;
1594 Value gcd,tmp;
1595
1596 value_init(gcd); value_init(tmp);
1597
1598 fprintf(DST, "[" );
1599 for(l=0;l<V->NbRows;++l){
1600
1601 /* Variables */
1602 first=1;
1603 fprintf(DST, " " );
1604 for(v=0;v < V->NbColumns-2;++v) {
1605 if(value_notzero_p(V->p[l][v])) {
1606 value_gcd(gcd, V->p[l][v], V->p[l][V->NbColumns-1]);
1607 value_divexact(tmp, V->p[l][v], gcd);
1608 if(value_posz_p(tmp)) {
1609 if(!first)
1610 fprintf(DST, "+");
1611 if(value_notone_p(tmp)) {
1612 value_print(DST,VALUE_FMT,tmp);
1613 }
1614 }
1615 else { /* V->p[l][v]/gcd<0 */
1616 if(value_mone_p(tmp))
1617 fprintf(DST, "-" );
1618 else {
1619 value_print(DST,VALUE_FMT,tmp);
1620 }
1621 }
1622 value_divexact(tmp, V->p[l][V->NbColumns-1], gcd);
1623 if(value_notone_p(tmp)) {
1624 fprintf(DST, "%s/", param_names[v]);
1625 value_print(DST,VALUE_FMT,tmp);
1626 }
1627 else
1628 fprintf(DST, "%s", param_names[v]);
1629 first=0;
1630 }
1631 }
1632
1633 /* Constant */
1634 if(value_notzero_p(V->p[l][v]) || first) {
1635 if(value_posz_p(V->p[l][v]) && !first)
1636 fprintf(DST,"+");
1637 value_gcd(gcd, V->p[l][v], V->p[l][V->NbColumns-1]);
1638 value_divexact(tmp, V->p[l][v], gcd);
1639 value_print(DST,VALUE_FMT,tmp);
1640 value_divexact(tmp, V->p[l][V->NbColumns-1], gcd);
1641 if(value_notone_p(tmp)) {
1642 fprintf(DST,"/");
1643 value_print(DST,VALUE_FMT,tmp);
1644 fprintf(DST," ");
1645 }
1646 }
1647 if (l<V->NbRows-1)
1648 fprintf(DST, ", ");
1649 }
1650 fprintf(DST, " ]");
1651 value_clear(gcd); value_clear(tmp);
1652 return;
1653 } /* Print_Vertex */
1654
1655 /*----------------------------------------------------------------------*/
1656 /* VertexCT */
1657 /* convert a paramvertex from reduced space to normal m-space */
1658 /*----------------------------------------------------------------------*/
1659 Matrix *VertexCT(Matrix *V,Matrix *CT) {
1660
1661 Matrix *Vt;
1662 int i,j,k;
1663
1664 if(CT) {
1665
1666 /* Have to transform the vertices to original dimension */
1667 Vt = Matrix_Alloc(V->NbRows,CT->NbColumns+1);
1668 for(i=0;i<V->NbRows;++i) {
1669 value_assign(Vt->p[i][CT->NbColumns],V->p[i][V->NbColumns-1]);
1670 for(j=0;j<CT->NbColumns;j++) {
1671 for(k=0;k<CT->NbRows;k++)
1672 if(value_notzero_p(CT->p[k][j]))
1673 break;
1674 if(k<CT->NbRows)
1675 value_assign(Vt->p[i][j],V->p[i][k]);
1676 else
1677 value_set_si(Vt->p[i][j],0);
1678 }
1679 }
1680 return(Vt);
1681 }
1682 else
1683 return(NULL);
1684 } /* VertexCT */
1685
1686 /*
1687 * Print the validity Domain 'D' of a parametric polyhedron
1688 */
1689 void Print_Domain(FILE *DST, Polyhedron *D, const char **pname)
1690 {
1691 int l, v;
1692 int first;
1693
1694 POL_ENSURE_FACETS(D);
1695 POL_ENSURE_VERTICES(D);
1696
1697 for(l=0;l<D->NbConstraints;++l) {
1698 fprintf(DST, " ");
1699 first = 1;
1700 for(v=1;v<=D->Dimension;++v) {
1701 if(value_notzero_p(D->Constraint[l][v])) {
1702 if(value_one_p(D->Constraint[l][v])) {
1703 if(first)
1704 fprintf(DST, "%s ", pname[v-1]);
1705 else
1706 fprintf(DST, "+ %s ", pname[v-1] );
1707 }
1708 else if(value_mone_p(D->Constraint[l][v]))
1709 fprintf(DST, "- %s ", pname[v-1] );
1710 else {
1711 if(value_pos_p(D->Constraint[l][v]) && !first )
1712 fprintf(DST, "+ " );
1713 value_print(DST,VALUE_FMT,D->Constraint[l][v]);
1714 fprintf(DST,"%s ",pname[v-1]);
1715 }
1716 first = 0;
1717 }
1718 }
1719 if(value_notzero_p(D->Constraint[l][v])) {
1720 if(value_pos_p(D->Constraint[l][v]) && !first)
1721 fprintf(DST,"+");
1722 fprintf(DST," ");
1723 value_print(DST,VALUE_FMT,D->Constraint[l][v]);
1724 }
1725 fprintf(DST,(value_notzero_p(D->Constraint[l][0])) ?" >= 0":" = 0");
1726 fprintf(DST, "\n" );
1727 }
1728 fprintf(DST, "\n");
1729
1730 if( D->next )
1731 {
1732 fprintf( DST, "UNION\n" );
1733 Print_Domain( DST, D->next, pname );
1734 }
1735 return;
1736 } /* Print_Domain */
1737
1738 /*
1739 * Given a list of parametrized vertices and an array of parameter names, Print
1740 * a list of parametrized vertices in a comprehensible format.
1741 */
1742 void Param_Vertices_Print(FILE *DST, Param_Vertices *PV, const char **param_names)
1743 {
1744 Polyhedron *poly;
1745
1746 while(PV) {
1747 fprintf(DST, "Vertex :\n" );
1748 Print_Vertex(DST,PV->Vertex,param_names);
1749
1750 /* Pour le domaine : */
1751 fprintf(DST, " If :\n" );
1752 poly = Constraints2Polyhedron(PV->Domain,200);
1753 Print_Domain(DST,poly,param_names);
1754 Domain_Free(poly);
1755 PV = PV->next;
1756 }
1757 return;
1758 } /* Param_Vertices_Print */
1759
1760 /*
1761 * Given a polyhedron 'Din' in combined data and parametre space, a context
1762 * polyhedron 'Cin' representing the constraints on the parameter space and
1763 * a working space size 'working_space', return a parametric polyhedron with
1764 * a list of distinct validity domains and a complete list of valid vertices
1765 * associated to each validity domain.
1766 */
1767 Param_Polyhedron *Polyhedron2Param_Domain(Polyhedron *Din,Polyhedron *Cin,int working_space) {
1768
1769 Param_Polyhedron *result;
1770 Param_Domain *D;
1771
1772 POL_ENSURE_FACETS(Din);
1773 POL_ENSURE_VERTICES(Din);
1774 POL_ENSURE_FACETS(Cin);
1775 POL_ENSURE_VERTICES(Cin);
1776
1777 if (emptyQ(Din) || emptyQ(Cin))
1778 return NULL;
1779
1780 #ifdef DEBUGPP
1781 fprintf(stderr,"Polyhedron2Param_Polyhedron algorithm starting at : %.2fs\n",
1782 (float)clock()/CLOCKS_PER_SEC);
1783 #endif
1784
1785 /* Find the m-faces, keeping the corresponding domains */
1786 /* in the linked list PDomains */
1787 result = Find_m_faces(&Din,Cin,1,working_space,NULL,NULL);
1788
1789 #ifdef DEBUGPP
1790 if(result) fprintf(stderr, "Number of vertices : %d\n",result->nbV);
1791 fprintf(stderr,"Vertices found at : %.2fs\n",(float)clock()/CLOCKS_PER_SEC);
1792 #endif
1793
1794 /* Processing of PVResult and PDomains */
1795 if(result && Cin->Dimension>0) /* at least 1 parameter */
1796 Compute_PDomains(result->D,result->nbV,working_space);
1797 if(result && CEqualities)
1798 for(D=result->D;D;D=D->next)
1799 D->Domain = Add_CEqualities(D->Domain);
1800 Polyhedron_Free(CEqualities);
1801
1802 #ifdef DEBUGPP
1803 fprintf(stderr, "domains found at : %.2fs\n", (float)clock()/CLOCKS_PER_SEC);
1804 #endif
1805
1806 return(result);
1807 } /* Polyhedon2Param_Domain */
1808
1809 /*
1810 *
1811 */
1812 Param_Polyhedron *Polyhedron2Param_SimplifiedDomain(Polyhedron **Din,Polyhedron *Cin,int working_space,Polyhedron **CEq,Matrix **CT) {
1813
1814 Param_Polyhedron *result;
1815
1816 assert(CEq != NULL);
1817 assert(CT != NULL);
1818
1819 POL_ENSURE_FACETS(*Din);
1820 POL_ENSURE_VERTICES(*Din);
1821 POL_ENSURE_FACETS(Cin);
1822 POL_ENSURE_VERTICES(Cin);
1823
1824 #ifdef DEBUGPP
1825 fprintf(stderr,"Polyhedron2Param_Polyhedron algorithm starting at : %.2fs\n",
1826 (float)clock()/CLOCKS_PER_SEC);
1827 #endif
1828
1829 /* Find the m-faces, keeping the corresponding domains */
1830 /* in the linked list PDomains */
1831 result = Find_m_faces(Din,Cin,1,working_space,CEq,CT);
1832
1833 #ifdef DEBUGPP
1834 if(result) fprintf(stderr, "Number of vertices : %d\n",result->nbV);
1835 fprintf(stderr,"Vertices found at : %.2fs\n",(float)clock()/CLOCKS_PER_SEC);
1836 #endif
1837
1838 /* Processing of PVResult and PDomains */
1839 if(result && Cin->Dimension>0) /* at least 1 parameter */
1840 Compute_PDomains(result->D,result->nbV,working_space);
1841
1842 /* Removed this, Vin100, March 01 */
1843 /* if(result && CEqualities )
1844 for(D=result->D;D;D=D->next)
1845 D->Domain = Add_CEqualities(D->Domain);
1846 */
1847
1848 #ifdef DEBUGPP
1849 fprintf(stderr, "domains found at : %.2fs\n", (float)clock()/CLOCKS_PER_SEC);
1850 #endif
1851
1852 return(result);
1853 } /* Polyhedron2Param_SimplifiedDomain */
1854
1855 /*
1856 * Free the memory allocated to a list of validity domain of a parametrized
1857 * polyhedron.
1858 */
1859 void Param_Domain_Free(Param_Domain *PD) {
1860
1861 Param_Domain *next_pd;
1862
1863 while(PD) {
1864 free(PD->F);
1865 Domain_Free(PD->Domain);
1866 next_pd = PD->next;
1867 free(PD);
1868 PD = next_pd;
1869 }
1870 return;
1871 } /* Param_Domain_Free */
1872
1873 /*
1874 * Free the memory allocated to a parametric polyhedron 'P'
1875 */
1876 void Param_Polyhedron_Free(Param_Polyhedron *P) {
1877
1878 if (!P) return;
1879 Param_Vertices_Free(P->V);
1880 Param_Domain_Free(P->D);
1881 if (P->Constraints)
1882 Matrix_Free(P->Constraints);
1883 if (P->Rays)
1884 Matrix_Free(P->Rays);
1885 free(P);
1886 return;
1887 } /* Param_Polyhedron_Free */
1888
1889 /*
1890 * Scales the parametric polyhedron such that all vertices are integer.
1891 */
1892 void Param_Polyhedron_Scale_Integer(Param_Polyhedron *PP, Polyhedron **P,
1893 Value *det, unsigned MaxRays)
1894 {
1895 int i;
1896 int nb_param, nb_vars;
1897 Vector *denoms;
1898 Param_Vertices *V;
1899 Value global_var_lcm;
1900 Matrix *expansion;
1901
1902 value_set_si(*det, 1);
1903 if (!PP->nbV)
1904 return;
1905
1906 nb_param = PP->D->Domain->Dimension;
1907 nb_vars = PP->V->Vertex->NbRows;
1908
1909 /* Scan the vertices and make an orthogonal expansion of the variable
1910 space */
1911 /* a- prepare the array of common denominators */
1912 denoms = Vector_Alloc(nb_vars);
1913 value_init(global_var_lcm);
1914
1915 /* b- scan the vertices and compute the variables' global lcms */
1916 for (V = PP->V; V; V = V->next)
1917 for (i = 0; i < nb_vars; i++)
1918 value_lcm(denoms->p[i], denoms->p[i], V->Vertex->p[i][nb_param+1]);
1919
1920 value_set_si(global_var_lcm, 1);
1921 for (i = 0; i < nb_vars; i++) {
1922 value_multiply(*det, *det, denoms->p[i]);
1923 value_lcm(global_var_lcm, global_var_lcm, denoms->p[i]);
1924 }
1925
1926 /* scale vertices */
1927 for (V = PP->V; V; V = V->next)
1928 for (i = 0; i < nb_vars; i++) {
1929 Vector_Scale(V->Vertex->p[i], V->Vertex->p[i], denoms->p[i], nb_param+1);
1930 Vector_Normalize(V->Vertex->p[i], nb_param+2);
1931 }
1932
1933 /* the expansion can be actually writen as global_var_lcm.L^{-1} */
1934 /* this is equivalent to multiply the rows of P by denoms_det */
1935 for (i = 0; i < nb_vars; i++)
1936 value_division(denoms->p[i], global_var_lcm, denoms->p[i]);
1937
1938 /* OPT : we could use a vector instead of a diagonal matrix here (c- and d-).*/
1939 /* c- make the quick expansion matrix */
1940 expansion = Matrix_Alloc(nb_vars+nb_param+1, nb_vars+nb_param+1);
1941 for (i = 0; i < nb_vars; i++)
1942 value_assign(expansion->p[i][i], denoms->p[i]);
1943 for (i = nb_vars; i < nb_vars+nb_param+1; i++)
1944 value_assign(expansion->p[i][i], global_var_lcm);
1945
1946 /* d- apply the variable expansion to the polyhedron */
1947 if (P)
1948 *P = Polyhedron_Preimage(*P, expansion, MaxRays);
1949
1950 Matrix_Free(expansion);
1951 value_clear(global_var_lcm);
1952 Vector_Free(denoms);
1953 }
Polyhedral library