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 * $Id: compress_parms.c,v 1.32 2006/11/03 17:34:26 skimo Exp $ 20 * 21 * The integer points in a parametric linear subspace of Q^n are generally 22 * lying on a sub-lattice of Z^n. 23 * Functions here use and compute validity lattices, i.e. lattices induced on a 24 * set of variables by such equalities involving another set of integer 25 * variables. 26 * @author B. Meister 12/2003-2006 meister@icps.u-strasbg.fr 27 * LSIIT -ICPS 28 * UMR 7005 CNRS 29 * Louis Pasteur University (ULP), Strasbourg, France 30 */ 31 32 #include <stdlib.h> 33 #include <polylib/polylib.h> 34 35 /** 36 * debug flags (2 levels) 37 */ 38 #define dbgCompParm 0 39 #define dbgCompParmMore 0 40 41 #define dbgStart(a) if (dbgCompParmMore) { printf(" -- begin "); \ 42 printf(#a); \ 43 printf(" --\n"); } \ 44 while(0) 45 #define dbgEnd(a) if (dbgCompParmMore) { printf(" -- end "); \ 46 printf(#a); \ 47 printf(" --\n"); } \ 48 while(0) 49 50 51 /** 52 * Given a full-row-rank nxm matrix M made of m row-vectors), computes the 53 * basis K (made of n-m column-vectors) of the integer kernel of the rows of M 54 * so we have: M.K = 0 55 */ 56 Matrix * int_ker(Matrix * M) { 57 Matrix *U, *Q, *H, *H2, *K=NULL; 58 int i, j, rk; 59 60 if (dbgCompParm) 61 show_matrix(M); 62 /* eliminate redundant rows : UM = H*/ 63 right_hermite(M, &H, &Q, &U); 64 for (rk=H->NbRows-1; (rk>=0) && Vector_IsZero(H->p[rk], H->NbColumns); rk--); 65 rk++; 66 if (dbgCompParmMore) { 67 printf("rank = %d\n", rk); 68 } 69 70 /* there is a non-null kernel if and only if the dimension m of 71 the space spanned by the rows 72 is inferior to the number n of variables */ 73 if (M->NbColumns <= rk) { 74 Matrix_Free(H); 75 Matrix_Free(Q); 76 Matrix_Free(U); 77 K = Matrix_Alloc(M->NbColumns, 0); 78 return K; 79 } 80 Matrix_Free(U); 81 Matrix_Free(Q); 82 /* fool left_hermite by giving NbRows =rank of M*/ 83 H->NbRows=rk; 84 /* computes MU = [H 0] */ 85 left_hermite(H, &H2, &Q, &U); 86 if (dbgCompParmMore) { 87 printf("-- Int. Kernel -- \n"); 88 show_matrix(M); 89 printf(" = \n"); 90 show_matrix(H2); 91 show_matrix(U); 92 } 93 H->NbRows==M->NbRows; 94 Matrix_Free(H); 95 /* the Integer Kernel is made of the last n-rk columns of U */ 96 Matrix_subMatrix(U, 0, rk, U->NbRows, U->NbColumns, &K); 97 98 /* clean up */ 99 Matrix_Free(H2); 100 Matrix_Free(U); 101 Matrix_Free(Q); 102 return K; 103 } /* int_ker */ 104 105 106 /** 107 * Computes the intersection of two linear lattices, whose base vectors are 108 * respectively represented in A and B. 109 * If I and/or Lb is set to NULL, then the matrix is allocated. 110 * Else, the matrix is assumed to be allocated already. 111 * I and Lb are rk x rk, where rk is the rank of A (or B). 112 * @param A the full-row rank matrix whose column-vectors are the basis for the 113 * first linear lattice. 114 * @param B the matrix whose column-vectors are the basis for the second linear 115 * lattice. 116 * @param Lb the matrix such that B.Lb = I, where I is the intersection. 117 * @return their intersection. 118 */ 119 static void linearInter(Matrix * A, Matrix * B, Matrix ** I, Matrix **Lb) { 120 Matrix * AB=NULL; 121 int rk = A->NbRows; 122 int a = A->NbColumns; 123 int b = B->NbColumns; 124 int i,j, z=0; 125 126 Matrix * H, *U, *Q; 127 /* ensure that the spanning vectors are in the same space */ 128 assert(B->NbRows==rk); 129 /* 1- build the matrix 130 * (A 0 1) 131 * (0 B 1) 132 */ 133 AB = Matrix_Alloc(2*rk, a+b+rk); 134 Matrix_copySubMatrix(A, 0, 0, rk, a, AB, 0, 0); 135 Matrix_copySubMatrix(B, 0, 0, rk, b, AB, rk, a); 136 for (i=0; i< rk; i++) { 137 value_set_si(AB->p[i][a+b+i], 1); 138 value_set_si(AB->p[i+rk][a+b+i], 1); 139 } 140 if (dbgCompParm) { 141 show_matrix(AB); 142 } 143 144 /* 2- Compute its left Hermite normal form. AB.U = [H 0] */ 145 left_hermite(AB, &H, &Q, &U); 146 Matrix_Free(AB); 147 Matrix_Free(Q); 148 /* count the number of non-zero colums in H */ 149 for (z=H->NbColumns-1; value_zero_p(H->p[H->NbRows-1][z]); z--); 150 z++; 151 if (dbgCompParm) { 152 show_matrix(H); 153 printf("z=%d\n", z); 154 } 155 Matrix_Free(H); 156 /* if you split U in 9 submatrices, you have: 157 * A.U_13 = -U_33 158 * B.U_23 = -U_33, 159 * where the nb of cols of U_{*3} equals the nb of zero-cols of H 160 * U_33 is a (the smallest) combination of col-vectors of A and B at the same 161 * time: their intersection. 162 */ 163 Matrix_subMatrix(U, a+b, z, U->NbColumns, U->NbColumns, I); 164 Matrix_subMatrix(U, a, z, a+b, U->NbColumns, Lb); 165 if (dbgCompParm) { 166 show_matrix(U); 167 } 168 Matrix_Free(U); 169 } /* linearInter */ 170 171 172 /** 173 * Given a system of equalities, looks if it has an integer solution in the 174 * combined space, and if yes, returns one solution. 175 * <p>pre-condition: the equalities are full-row rank (without the constant 176 * part)</p> 177 * @param Eqs the system of equations (as constraints) 178 * @param I a feasible integer solution if it exists, else NULL. Allocated if 179 * initially set to NULL, else reused. 180 */ 181 void Equalities_integerSolution(Matrix * Eqs, Matrix **I) { 182 Matrix * Hm, *H=NULL, *U, *Q, *M=NULL, *C=NULL, *Hi; 183 Matrix *Ip; 184 int i; 185 Value mod; 186 unsigned int rk; 187 if (Eqs==NULL){ 188 if ((*I)!=NULL) Matrix_Free(*I); 189 I = NULL; 190 return; 191 } 192 /* we use: AI = C = (Ha 0).Q.I = (Ha 0)(I' 0)^T */ 193 /* with I = Qinv.I' = U.I'*/ 194 /* 1- compute I' = Hainv.(-C) */ 195 /* HYP: the equalities are full-row rank */ 196 rk = Eqs->NbRows; 197 Matrix_subMatrix(Eqs, 0, 1, rk, Eqs->NbColumns-1, &M); 198 left_hermite(M, &Hm, &Q, &U); 199 Matrix_Free(M); 200 Matrix_subMatrix(Hm, 0, 0, rk, rk, &H); 201 if (dbgCompParmMore) { 202 show_matrix(Hm); 203 show_matrix(H); 204 show_matrix(U); 205 } 206 Matrix_Free(Q); 207 Matrix_Free(Hm); 208 Matrix_subMatrix(Eqs, 0, Eqs->NbColumns-1, rk, Eqs->NbColumns, &C); 209 Matrix_oppose(C); 210 Hi = Matrix_Alloc(rk, rk+1); 211 MatInverse(H, Hi); 212 if (dbgCompParmMore) { 213 show_matrix(C); 214 show_matrix(Hi); 215 } 216 /* put the numerator of Hinv back into H */ 217 Matrix_subMatrix(Hi, 0, 0, rk, rk, &H); 218 Ip = Matrix_Alloc(Eqs->NbColumns-2, 1); 219 /* fool Matrix_Product on the size of Ip */ 220 Ip->NbRows = rk; 221 Matrix_Product(H, C, Ip); 222 Ip->NbRows = Eqs->NbColumns-2; 223 Matrix_Free(H); 224 Matrix_Free(C); 225 value_init(mod); 226 for (i=0; i< rk; i++) { 227 /* if Hinv.C is not integer, return NULL (no solution) */ 228 value_pmodulus(mod, Ip->p[i][0], Hi->p[i][rk]); 229 if (value_notzero_p(mod)) { 230 if ((*I)!=NULL) Matrix_Free(*I); 231 value_clear(mod); 232 Matrix_Free(U); 233 Matrix_Free(Ip); 234 Matrix_Free(Hi); 235 I = NULL; 236 return; 237 } 238 else { 239 value_pdivision(Ip->p[i][0], Ip->p[i][0], Hi->p[i][rk]); 240 } 241 } 242 /* fill the rest of I' with zeros */ 243 for (i=rk; i< Eqs->NbColumns-2; i++) { 244 value_set_si(Ip->p[i][0], 0); 245 } 246 value_clear(mod); 247 Matrix_Free(Hi); 248 /* 2 - Compute the particular solution I = U.(I' 0) */ 249 ensureMatrix((*I), Eqs->NbColumns-2, 1); 250 Matrix_Product(U, Ip, (*I)); 251 Matrix_Free(U); 252 Matrix_Free(Ip); 253 if (dbgCompParm) { 254 show_matrix(*I); 255 } 256 } 257 258 259 /** 260 * Computes the validity lattice of a set of equalities. I.e., the lattice 261 * induced on the last <tt>b</tt> variables by the equalities involving the 262 * first <tt>a</tt> integer existential variables. The submatrix of Eqs that 263 * concerns only the existential variables (so the first a columns) is assumed 264 * to be full-row rank. 265 * @param Eqs the equalities 266 * @param a the number of existential integer variables, placed as first 267 * variables 268 * @param vl the (returned) validity lattice, in homogeneous form. It is 269 * allocated if initially set to null, or reused if already allocated. 270 */ 271 void Equalities_validityLattice(Matrix * Eqs, int a, Matrix** vl) { 272 unsigned int b = Eqs->NbColumns-2-a; 273 unsigned int r = Eqs->NbRows; 274 Matrix * A=NULL, * B=NULL, *I = NULL, *Lb=NULL, *sol=NULL; 275 Matrix *H, *U, *Q; 276 unsigned int i; 277 278 if (dbgCompParm) { 279 printf("Computing validity lattice induced by the %d first variables of:" 280 ,a); 281 show_matrix(Eqs); 282 } 283 if (b==0) { 284 ensureMatrix((*vl), 1, 1); 285 value_set_si((*vl)->p[0][0], 1); 286 return; 287 } 288 289 /* 1- check that there is an integer solution to the equalities */ 290 /* OPT: could change integerSolution's profile to allocate or not*/ 291 Equalities_integerSolution(Eqs, &sol); 292 /* if there is no integer solution, there is no validity lattice */ 293 if (sol==NULL) { 294 if ((*vl)!=NULL) Matrix_Free(*vl); 295 return; 296 } 297 Matrix_subMatrix(Eqs, 0, 1, r, 1+a, &A); 298 Matrix_subMatrix(Eqs, 0, 1+a, r, 1+a+b, &B); 299 linearInter(A, B, &I, &Lb); 300 Matrix_Free(A); 301 Matrix_Free(B); 302 Matrix_Free(I); 303 if (dbgCompParm) { 304 show_matrix(Lb); 305 } 306 307 /* 2- The linear part of the validity lattice is the left HNF of Lb */ 308 left_hermite(Lb, &H, &Q, &U); 309 Matrix_Free(Lb); 310 Matrix_Free(Q); 311 Matrix_Free(U); 312 313 /* 3- build the validity lattice */ 314 ensureMatrix((*vl), b+1, b+1); 315 Matrix_copySubMatrix(H, 0, 0, b, b, (*vl), 0,0); 316 Matrix_Free(H); 317 for (i=0; i< b; i++) { 318 value_assign((*vl)->p[i][b], sol->p[0][a+i]); 319 } 320 Matrix_Free(sol); 321 Vector_Set((*vl)->p[b],0, b); 322 value_set_si((*vl)->p[b][b], 1); 323 324 } /* validityLattice */ 325 326 327 /** 328 * Eliminate the columns corresponding to a list of eliminated parameters. 329 * @param M the constraints matrix whose columns are to be removed 330 * @param nbVars an offset to be added to the ranks of the variables to be 331 * removed 332 * @param elimParms the list of ranks of the variables to be removed 333 * @param newM (output) the matrix without the removed columns 334 */ 335 void Constraints_removeElimCols(Matrix * M, unsigned int nbVars, 336 unsigned int *elimParms, Matrix ** newM) { 337 unsigned int i, j, k; 338 if (elimParms[0]==0) { 339 Matrix_clone(M, newM); 340 return; 341 } 342 if ((*newM)==NULL) { 343 (*newM) = Matrix_Alloc(M->NbRows, M->NbColumns - elimParms[0]); 344 } 345 else { 346 assert ((*newM)->NbColumns==M->NbColumns - elimParms[0]); 347 } 348 for (i=0; i< M->NbRows; i++) { 349 value_assign((*newM)->p[i][0], M->p[i][0]); /* kind of cstr */ 350 k=0; 351 Vector_Copy(&(M->p[i][1]), &((*newM)->p[i][1]), nbVars); 352 for (j=0; j< M->NbColumns-2-nbVars; j++) { 353 if (j!=elimParms[k+1]) { 354 value_assign((*newM)->p[i][j-k+nbVars+1], M->p[i][j+nbVars+1]); 355 } 356 else { 357 k++; 358 } 359 } 360 value_assign((*newM)->p[i][(*newM)->NbColumns-1], 361 M->p[i][M->NbColumns-1]); /* cst part */ 362 } 363 } /* Constraints_removeElimCols */ 364 365 366 /** 367 * Eliminates all the equalities in a set of constraints and returns the set of 368 * constraints defining a full-dimensional polyhedron, such that there is a 369 * bijection between integer points of the original polyhedron and these of the 370 * resulting (projected) polyhedron). 371 * If VL is set to NULL, this funciton allocates it. Else, it assumes that 372 * (*VL) points to a matrix of the right size. 373 * <p> The following things are done: 374 * <ol> 375 * <li> remove equalities involving only parameters, and remove as many 376 * parameters as there are such equalities. From that, the list of 377 * eliminated parameters <i>elimParms</i> is built. 378 * <li> remove equalities that involve variables. This requires a compression 379 * of the parameters and of the other variables that are not eliminated. 380 * The affine compresson is represented by matrix VL (for <i>validity 381 * lattice</i>) and is such that (N I 1)^T = VL.(N' I' 1), where N', I' 382 * are integer (they are the parameters and variables after compression). 383 *</ol> 384 *</p> 385 */ 386 void Constraints_fullDimensionize(Matrix ** M, Matrix ** C, Matrix ** VL, 387 Matrix ** Eqs, Matrix ** ParmEqs, 388 unsigned int ** elimVars, 389 unsigned int ** elimParms, 390 int maxRays) { 391 unsigned int i, j; 392 Matrix * A=NULL, *B=NULL; 393 Matrix * Ineqs=NULL; 394 unsigned int nbVars = (*M)->NbColumns - (*C)->NbColumns; 395 unsigned int nbParms; 396 int nbElimVars; 397 Matrix * fullDim = NULL; 398 399 /* variables for permutations */ 400 unsigned int * permutation; 401 Matrix * permutedEqs=NULL, * permutedIneqs=NULL; 402 403 /* 1- Eliminate the equalities involving only parameters. */ 404 (*ParmEqs) = Constraints_removeParmEqs(M, C, 0, elimParms); 405 /* if the polyehdron is empty, return now. */ 406 if ((*M)->NbColumns==0) return; 407 /* eliminate the columns corresponding to the eliminated parameters */ 408 if (elimParms[0]!=0) { 409 Constraints_removeElimCols(*M, nbVars, (*elimParms), &A); 410 Matrix_Free(*M); 411 (*M) = A; 412 Constraints_removeElimCols(*C, 0, (*elimParms), &B); 413 Matrix_Free(*C); 414 (*C) = B; 415 if (dbgCompParm) { 416 printf("After false parameter elimination: \n"); 417 show_matrix(*M); 418 show_matrix(*C); 419 } 420 } 421 nbParms = (*C)->NbColumns-2; 422 423 /* 2- Eliminate the equalities involving variables */ 424 /* a- extract the (remaining) equalities from the poyhedron */ 425 split_constraints((*M), Eqs, &Ineqs); 426 nbElimVars = (*Eqs)->NbRows; 427 /* if the polyhedron is already full-dimensional, return */ 428 if ((*Eqs)->NbRows==0) { 429 Matrix_identity(nbParms+1, VL); 430 return; 431 } 432 /* b- choose variables to be eliminated */ 433 permutation = find_a_permutation((*Eqs), nbParms); 434 435 if (dbgCompParm) { 436 printf("Permuting the vars/parms this way: [ "); 437 for (i=0; i< (*Eqs)->NbColumns-2; i++) { 438 printf("%d ", permutation[i]); 439 } 440 printf("]\n"); 441 } 442 443 Constraints_permute((*Eqs), permutation, &permutedEqs); 444 Equalities_validityLattice(permutedEqs, (*Eqs)->NbRows, VL); 445 446 if (dbgCompParm) { 447 printf("Validity lattice: "); 448 show_matrix(*VL); 449 } 450 Constraints_compressLastVars(permutedEqs, (*VL)); 451 Constraints_permute(Ineqs, permutation, &permutedIneqs); 452 if (dbgCompParmMore) { 453 show_matrix(permutedIneqs); 454 show_matrix(permutedEqs); 455 } 456 Matrix_Free(*Eqs); 457 Matrix_Free(Ineqs); 458 Constraints_compressLastVars(permutedIneqs, (*VL)); 459 if (dbgCompParm) { 460 printf("After compression: "); 461 show_matrix(permutedIneqs); 462 } 463 /* c- eliminate the first variables */ 464 assert(Constraints_eliminateFirstVars(permutedEqs, permutedIneqs)); 465 if (dbgCompParmMore) { 466 printf("After elimination of the variables: "); 467 show_matrix(permutedIneqs); 468 } 469 470 /* d- get rid of the first (zero) columns, 471 which are now useless, and put the parameters back at the end */ 472 fullDim = Matrix_Alloc(permutedIneqs->NbRows, 473 permutedIneqs->NbColumns-nbElimVars); 474 for (i=0; i< permutedIneqs->NbRows; i++) { 475 value_set_si(fullDim->p[i][0], 1); 476 for (j=0; j< nbParms; j++) { 477 value_assign(fullDim->p[i][j+fullDim->NbColumns-nbParms-1], 478 permutedIneqs->p[i][j+nbElimVars+1]); 479 } 480 for (j=0; j< permutedIneqs->NbColumns-nbParms-2-nbElimVars; j++) { 481 value_assign(fullDim->p[i][j+1], 482 permutedIneqs->p[i][nbElimVars+nbParms+j+1]); 483 } 484 value_assign(fullDim->p[i][fullDim->NbColumns-1], 485 permutedIneqs->p[i][permutedIneqs->NbColumns-1]); 486 } 487 Matrix_Free(permutedIneqs); 488 489 } /* Constraints_fullDimensionize */ 490 491 492 /** 493 * Given a matrix that defines a full-dimensional affine lattice, returns the 494 * affine sub-lattice spanned in the k first dimensions. 495 * Useful for instance when you only look for the parameters' validity lattice. 496 * @param lat the original full-dimensional lattice 497 * @param subLat the sublattice 498 */ 499 void Lattice_extractSubLattice(Matrix * lat, unsigned int k, Matrix ** subLat) { 500 Matrix * H, *Q, *U, *linLat = NULL; 501 unsigned int i; 502 dbgStart(Lattice_extractSubLattice); 503 /* if the dimension is already good, just copy the initial lattice */ 504 if (k==lat->NbRows-1) { 505 if (*subLat==NULL) { 506 (*subLat) = Matrix_Copy(lat); 507 } 508 else { 509 Matrix_copySubMatrix(lat, 0, 0, lat->NbRows, lat->NbColumns, (*subLat), 0, 0); 510 } 511 return; 512 } 513 assert(k<lat->NbRows-1); 514 /* 1- Make the linear part of the lattice triangular to eliminate terms from 515 other dimensions */ 516 Matrix_subMatrix(lat, 0, 0, lat->NbRows, lat->NbColumns-1, &linLat); 517 /* OPT: any integer column-vector elimination is ok indeed. */ 518 /* OPT: could test if the lattice is already in triangular form. */ 519 left_hermite(linLat, &H, &Q, &U); 520 if (dbgCompParmMore) { 521 show_matrix(H); 522 } 523 Matrix_Free(Q); 524 Matrix_Free(U); 525 Matrix_Free(linLat); 526 /* if not allocated yet, allocate it */ 527 if (*subLat==NULL) { 528 (*subLat) = Matrix_Alloc(k+1, k+1); 529 } 530 Matrix_copySubMatrix(H, 0, 0, k, k, (*subLat), 0, 0); 531 Matrix_Free(H); 532 Matrix_copySubMatrix(lat, 0, lat->NbColumns-1, k, 1, (*subLat), 0, k); 533 for (i=0; i<k; i++) { 534 value_set_si((*subLat)->p[k][i], 0); 535 } 536 value_set_si((*subLat)->p[k][k], 1); 537 dbgEnd(Lattice_extractSubLattice); 538 } /* Lattice_extractSubLattice */ 539 540 541 /** 542 * Computes the overall period of the variables I for (MI) mod |d|, where M is 543 * a matrix and |d| a vector. Produce a diagonal matrix S = (s_k) where s_k is 544 * the overall period of i_k 545 * @param M the set of affine functions of I (row-vectors) 546 * @param d the column-vector representing the modulos 547 */ 548 Matrix * affine_periods(Matrix * M, Matrix * d) { 549 Matrix * S; 550 unsigned int i,j; 551 Value tmp; 552 Value * periods = (Value *)malloc(sizeof(Value) * M->NbColumns); 553 value_init(tmp); 554 for(i=0; i< M->NbColumns; i++) { 555 value_init(periods[i]); 556 value_set_si(periods[i], 1); 557 } 558 for (i=0; i<M->NbRows; i++) { 559 for (j=0; j< M->NbColumns; j++) { 560 value_gcd(tmp, d->p[i][0], M->p[i][j]); 561 value_divexact(tmp, d->p[i][0], tmp); 562 value_lcm(periods[j], periods[j], tmp); 563 } 564 } 565 value_clear(tmp); 566 567 /* 2- build S */ 568 S = Matrix_Alloc(M->NbColumns, M->NbColumns); 569 for (i=0; i< M->NbColumns; i++) 570 for (j=0; j< M->NbColumns; j++) 571 if (i==j) value_assign(S->p[i][j],periods[j]); 572 else value_set_si(S->p[i][j], 0); 573 574 /* 3- clean up */ 575 for(i=0; i< M->NbColumns; i++) value_clear(periods[i]); 576 free(periods); 577 return S; 578 } /* affine_periods */ 579 580 581 /** 582 * Given an integer matrix B with m rows and integer m-vectors C and d, 583 * computes the basis of the integer solutions to (BN+C) mod d = 0 (1). 584 * This is an affine lattice (G): (N 1)^T= G(N' 1)^T, forall N' in Z^b. 585 * If there is no solution, returns NULL. 586 * @param B B, a (m x b) matrix 587 * @param C C, a (m x 1) integer matrix 588 * @param d d, a (1 x m) integer matrix 589 * @param imb the affine (b+1)x(b+1) basis of solutions, in the homogeneous 590 * form. Allocated if initially set to NULL, reused if not. 591 */ 592 void Equalities_intModBasis(Matrix * B, Matrix * C, Matrix * d, Matrix ** imb) { 593 int b = B->NbColumns; 594 /* FIXME: treat the case d=0 as a regular equality B_kN+C_k = 0: */ 595 /* OPT: could keep only equalities for which d>1 */ 596 int nbEqs = B->NbRows; 597 unsigned int i; 598 599 /* 1- buid the problem DI+BN+C = 0 */ 600 Matrix * eqs = Matrix_Alloc(nbEqs, nbEqs+b+1); 601 for (i=0; i< nbEqs; i++) { 602 value_assign(eqs->p[i][i], d->p[0][i]); 603 } 604 Matrix_copySubMatrix(B, 0, 0, nbEqs, b, eqs, 0, nbEqs); 605 Matrix_copySubMatrix(C, 0, 0, nbEqs, 1, eqs, 0, nbEqs+b); 606 607 /* 2- the solution is the validity lattice of the equalities */ 608 Equalities_validityLattice(eqs, nbEqs, imb); 609 Matrix_Free(eqs); 610 } /* Equalities_intModBasis */ 611 612 613 /** kept here for backwards compatiblity. Wrapper to Equalities_intModBasis() */ 614 Matrix * int_mod_basis(Matrix * B, Matrix * C, Matrix * d) { 615 Matrix * imb = NULL; 616 Equalities_intModBasis(B, C, d, &imb); 617 return imb; 618 } /* int_mod_basis */ 619 620 621 /** 622 * Given a parameterized constraints matrix with m equalities, computes the 623 * compression matrix G such that there is an integer solution in the variables 624 * space for each value of N', with N = G N' (N are the "nb_parms" parameters) 625 * @param E a matrix of parametric equalities @param nb_parms the number of 626 * parameters 627 * <b>Note: </b>this function is mostly here for backwards 628 * compatibility. Prefer the use of <tt>Equalities_validityLattice</tt>. 629 */ 630 Matrix * compress_parms(Matrix * E, int nbParms) { 631 Matrix * vl=NULL; 632 Equalities_validityLattice(E, E->NbColumns-2-nbParms, &vl); 633 return vl; 634 }/* compress_parms */ 635 636 637 /** Removes the equalities that involve only parameters, by eliminating some 638 * parameters in the polyhedron's constraints and in the context.<p> 639 * <b>Updates M and Ctxt.</b> 640 * @param M1 the polyhedron's constraints 641 * @param Ctxt1 the constraints of the polyhedron's context 642 * @param renderSpace tells if the returned equalities must be expressed in the 643 * parameters space (renderSpace=0) or in the combined var/parms space 644 * (renderSpace = 1) 645 * @param elimParms the list of parameters that have been removed: an array 646 * whose 1st element is the number of elements in the list. (returned) 647 * @return the system of equalities that involve only parameters. 648 */ 649 Matrix * Constraints_Remove_parm_eqs(Matrix ** M1, Matrix ** Ctxt1, 650 int renderSpace, 651 unsigned int ** elimParms) { 652 int i, j, k, nbEqsParms =0; 653 int nbEqsM, nbEqsCtxt, allZeros, nbTautoM = 0, nbTautoCtxt = 0; 654 Matrix * M = (*M1); 655 Matrix * Ctxt = (*Ctxt1); 656 int nbVars = M->NbColumns-Ctxt->NbColumns; 657 Matrix * Eqs; 658 Matrix * EqsMTmp; 659 660 /* 1- build the equality matrix(ces) */ 661 nbEqsM = 0; 662 for (i=0; i< M->NbRows; i++) { 663 k = First_Non_Zero(M->p[i], M->NbColumns); 664 /* if it is a tautology, count it as such */ 665 if (k==-1) { 666 nbTautoM++; 667 } 668 else { 669 /* if it only involves parameters, count it */ 670 if (k>= nbVars+1) nbEqsM++; 671 } 672 } 673 674 nbEqsCtxt = 0; 675 for (i=0; i< Ctxt->NbRows; i++) { 676 if (value_zero_p(Ctxt->p[i][0])) { 677 if (First_Non_Zero(Ctxt->p[i], Ctxt->NbColumns)==-1) { 678 nbTautoCtxt++; 679 } 680 else { 681 nbEqsCtxt ++; 682 } 683 } 684 } 685 nbEqsParms = nbEqsM + nbEqsCtxt; 686 687 /* nothing to do in this case */ 688 if (nbEqsParms+nbTautoM+nbTautoCtxt==0) { 689 (*elimParms) = (unsigned int*) malloc(sizeof(int)); 690 (*elimParms)[0] = 0; 691 if (renderSpace==0) { 692 return Matrix_Alloc(0,Ctxt->NbColumns); 693 } 694 else { 695 return Matrix_Alloc(0,M->NbColumns); 696 } 697 } 698 699 Eqs= Matrix_Alloc(nbEqsParms, Ctxt->NbColumns); 700 EqsMTmp= Matrix_Alloc(nbEqsParms, M->NbColumns); 701 702 /* copy equalities from the context */ 703 k = 0; 704 for (i=0; i< Ctxt->NbRows; i++) { 705 if (value_zero_p(Ctxt->p[i][0]) 706 && First_Non_Zero(Ctxt->p[i], Ctxt->NbColumns)!=-1) { 707 Vector_Copy(Ctxt->p[i], Eqs->p[k], Ctxt->NbColumns); 708 Vector_Copy(Ctxt->p[i]+1, EqsMTmp->p[k]+nbVars+1, 709 Ctxt->NbColumns-1); 710 k++; 711 } 712 } 713 for (i=0; i< M->NbRows; i++) { 714 j=First_Non_Zero(M->p[i], M->NbColumns); 715 /* copy equalities that involve only parameters from M */ 716 if (j>=nbVars+1) { 717 Vector_Copy(M->p[i]+nbVars+1, Eqs->p[k]+1, Ctxt->NbColumns-1); 718 Vector_Copy(M->p[i]+nbVars+1, EqsMTmp->p[k]+nbVars+1, 719 Ctxt->NbColumns-1); 720 /* mark these equalities for removal */ 721 value_set_si(M->p[i][0], 2); 722 k++; 723 } 724 /* mark the all-zero equalities for removal */ 725 if (j==-1) { 726 value_set_si(M->p[i][0], 2); 727 } 728 } 729 730 /* 2- eliminate parameters until all equalities are used or until we find a 731 contradiction (overconstrained system) */ 732 (*elimParms) = (unsigned int *) malloc((Eqs->NbRows+1) * sizeof(int)); 733 (*elimParms)[0] = 0; 734 allZeros = 0; 735 for (i=0; i< Eqs->NbRows; i++) { 736 /* find a variable that can be eliminated */ 737 k = First_Non_Zero(Eqs->p[i], Eqs->NbColumns); 738 if (k!=-1) { /* nothing special to do for tautologies */ 739 740 /* if there is a contradiction, return empty matrices */ 741 if (k==Eqs->NbColumns-1) { 742 printf("Contradiction in %dth row of Eqs: ",k); 743 show_matrix(Eqs); 744 Matrix_Free(Eqs); 745 Matrix_Free(EqsMTmp); 746 (*M1) = Matrix_Alloc(0, M->NbColumns); 747 Matrix_Free(M); 748 (*Ctxt1) = Matrix_Alloc(0,Ctxt->NbColumns); 749 Matrix_Free(Ctxt); 750 free(*elimParms); 751 (*elimParms) = (unsigned int *) malloc(sizeof(int)); 752 (*elimParms)[0] = 0; 753 if (renderSpace==1) { 754 return Matrix_Alloc(0,(*M1)->NbColumns); 755 } 756 else { 757 return Matrix_Alloc(0,(*Ctxt1)->NbColumns); 758 } 759 } 760 /* if we have something we can eliminate, do it in 3 places: 761 Eqs, Ctxt, and M */ 762 else { 763 k--; /* k is the rank of the variable, now */ 764 (*elimParms)[0]++; 765 (*elimParms)[(*elimParms[0])]=k; 766 for (j=0; j< Eqs->NbRows; j++) { 767 if (i!=j) { 768 eliminate_var_with_constr(Eqs, i, Eqs, j, k); 769 eliminate_var_with_constr(EqsMTmp, i, EqsMTmp, j, k+nbVars); 770 } 771 } 772 for (j=0; j< Ctxt->NbRows; j++) { 773 if (value_notzero_p(Ctxt->p[i][0])) { 774 eliminate_var_with_constr(Eqs, i, Ctxt, j, k); 775 } 776 } 777 for (j=0; j< M->NbRows; j++) { 778 if (value_cmp_si(M->p[i][0], 2)) { 779 eliminate_var_with_constr(EqsMTmp, i, M, j, k+nbVars); 780 } 781 } 782 } 783 } 784 /* if (k==-1): count the tautologies in Eqs to remove them later */ 785 else { 786 allZeros++; 787 } 788 } 789 790 /* elimParms may have been overallocated. Now we know how many parms have 791 been eliminated so we can reallocate the right amount of memory. */ 792 if (!realloc((*elimParms), ((*elimParms)[0]+1)*sizeof(int))) { 793 fprintf(stderr, "Constraints_Remove_parm_eqs > cannot realloc()"); 794 } 795 796 Matrix_Free(EqsMTmp); 797 798 /* 3- remove the "bad" equalities from the input matrices 799 and copy the equalities involving only parameters */ 800 EqsMTmp = Matrix_Alloc(M->NbRows-nbEqsM-nbTautoM, M->NbColumns); 801 k=0; 802 for (i=0; i< M->NbRows; i++) { 803 if (value_cmp_si(M->p[i][0], 2)) { 804 Vector_Copy(M->p[i], EqsMTmp->p[k], M->NbColumns); 805 k++; 806 } 807 } 808 Matrix_Free(M); 809 (*M1) = EqsMTmp; 810 811 EqsMTmp = Matrix_Alloc(Ctxt->NbRows-nbEqsCtxt-nbTautoCtxt, Ctxt->NbColumns); 812 k=0; 813 for (i=0; i< Ctxt->NbRows; i++) { 814 if (value_notzero_p(Ctxt->p[i][0])) { 815 Vector_Copy(Ctxt->p[i], EqsMTmp->p[k], Ctxt->NbColumns); 816 k++; 817 } 818 } 819 Matrix_Free(Ctxt); 820 (*Ctxt1) = EqsMTmp; 821 822 if (renderSpace==0) {/* renderSpace=0: equalities in the parameter space */ 823 EqsMTmp = Matrix_Alloc(Eqs->NbRows-allZeros, Eqs->NbColumns); 824 k=0; 825 for (i=0; i<Eqs->NbRows; i++) { 826 if (First_Non_Zero(Eqs->p[i], Eqs->NbColumns)!=-1) { 827 Vector_Copy(Eqs->p[i], EqsMTmp->p[k], Eqs->NbColumns); 828 k++; 829 } 830 } 831 } 832 else {/* renderSpace=1: equalities rendered in the combined space */ 833 EqsMTmp = Matrix_Alloc(Eqs->NbRows-allZeros, (*M1)->NbColumns); 834 k=0; 835 for (i=0; i<Eqs->NbRows; i++) { 836 if (First_Non_Zero(Eqs->p[i], Eqs->NbColumns)!=-1) { 837 Vector_Copy(Eqs->p[i], &(EqsMTmp->p[k][nbVars]), Eqs->NbColumns); 838 k++; 839 } 840 } 841 } 842 Matrix_Free(Eqs); 843 Eqs = EqsMTmp; 844 845 return Eqs; 846 } /* Constraints_Remove_parm_eqs */ 847 848 849 /** Removes equalities involving only parameters, but starting from a 850 * Polyhedron and its context. 851 * @param P the polyhedron 852 * @param C P's context 853 * @param renderSpace: 0 for the parameter space, =1 for the combined space. 854 * @maxRays Polylib's usual <i>workspace</i>. 855 */ 856 Polyhedron * Polyhedron_Remove_parm_eqs(Polyhedron ** P, Polyhedron ** C, 857 int renderSpace, 858 unsigned int ** elimParms, 859 int maxRays) { 860 Matrix * Eqs; 861 Polyhedron * Peqs; 862 Matrix * M = Polyhedron2Constraints((*P)); 863 Matrix * Ct = Polyhedron2Constraints((*C)); 864 865 /* if the Minkowski representation is not computed yet, do not compute it in 866 Constraints2Polyhedron */ 867 if (F_ISSET((*P), POL_VALID | POL_INEQUALITIES) && 868 (F_ISSET((*C), POL_VALID | POL_INEQUALITIES))) { 869 FL_INIT(maxRays, POL_NO_DUAL); 870 } 871 872 Eqs = Constraints_Remove_parm_eqs(&M, &Ct, renderSpace, elimParms); 873 Peqs = Constraints2Polyhedron(Eqs, maxRays); 874 Matrix_Free(Eqs); 875 876 /* particular case: no equality involving only parms is found */ 877 if (Eqs->NbRows==0) { 878 Matrix_Free(M); 879 Matrix_Free(Ct); 880 return Peqs; 881 } 882 Polyhedron_Free(*P); 883 Polyhedron_Free(*C); 884 (*P) = Constraints2Polyhedron(M, maxRays); 885 (*C) = Constraints2Polyhedron(Ct, maxRays); 886 Matrix_Free(M); 887 Matrix_Free(Ct); 888 return Peqs; 889 } /* Polyhedron_Remove_parm_eqs */ 890 891 892 /** 893 * Given a matrix with m parameterized equations, compress the nb_parms 894 * parameters and n-m variables so that m variables are integer, and transform 895 * the variable space into a n-m space by eliminating the m variables (using 896 * the equalities) the variables to be eliminated are chosen automatically by 897 * the function. 898 * <b>Deprecated.</b> Try to use Constraints_fullDimensionize instead. 899 * @param M the constraints 900 * @param the number of parameters 901 * @param validityLattice the the integer lattice underlying the integer 902 * solutions. 903 */ 904 Matrix * full_dimensionize(Matrix const * M, int nbParms, 905 Matrix ** validityLattice) { 906 Matrix * Eqs, * Ineqs; 907 Matrix * permutedEqs, * permutedIneqs; 908 Matrix * Full_Dim; 909 Matrix * WVL; /* The Whole Validity Lattice (vars+parms) */ 910 unsigned int i,j; 911 int nbElimVars; 912 unsigned int * permutation, * permutationInv; 913 /* 0- Split the equalities and inequalities from each other */ 914 split_constraints(M, &Eqs, &Ineqs); 915 916 /* 1- if the polyhedron is already full-dimensional, return it */ 917 if (Eqs->NbRows==0) { 918 Matrix_Free(Eqs); 919 (*validityLattice) = Identity_Matrix(nbParms+1); 920 return Ineqs; 921 } 922 nbElimVars = Eqs->NbRows; 923 924 /* 2- put the vars to be eliminated at the first positions, 925 and compress the other vars/parms 926 -> [ variables to eliminate / parameters / variables to keep ] */ 927 permutation = find_a_permutation(Eqs, nbParms); 928 if (dbgCompParm) { 929 printf("Permuting the vars/parms this way: [ "); 930 for (i=0; i< Eqs->NbColumns; i++) { 931 printf("%d ", permutation[i]); 932 } 933 printf("]\n"); 934 } 935 permutedEqs = mpolyhedron_permute(Eqs, permutation); 936 WVL = compress_parms(permutedEqs, Eqs->NbColumns-2-Eqs->NbRows); 937 if (dbgCompParm) { 938 printf("Whole validity lattice: "); 939 show_matrix(WVL); 940 } 941 mpolyhedron_compress_last_vars(permutedEqs, WVL); 942 permutedIneqs = mpolyhedron_permute(Ineqs, permutation); 943 if (dbgCompParm) { 944 show_matrix(permutedEqs); 945 } 946 Matrix_Free(Eqs); 947 Matrix_Free(Ineqs); 948 mpolyhedron_compress_last_vars(permutedIneqs, WVL); 949 if (dbgCompParm) { 950 printf("After compression: "); 951 show_matrix(permutedIneqs); 952 } 953 /* 3- eliminate the first variables */ 954 if (!mpolyhedron_eliminate_first_variables(permutedEqs, permutedIneqs)) { 955 fprintf(stderr,"full-dimensionize > variable elimination failed. \n"); 956 return NULL; 957 } 958 if (dbgCompParm) { 959 printf("After elimination of the variables: "); 960 show_matrix(permutedIneqs); 961 } 962 963 /* 4- get rid of the first (zero) columns, 964 which are now useless, and put the parameters back at the end */ 965 Full_Dim = Matrix_Alloc(permutedIneqs->NbRows, 966 permutedIneqs->NbColumns-nbElimVars); 967 for (i=0; i< permutedIneqs->NbRows; i++) { 968 value_set_si(Full_Dim->p[i][0], 1); 969 for (j=0; j< nbParms; j++) 970 value_assign(Full_Dim->p[i][j+Full_Dim->NbColumns-nbParms-1], 971 permutedIneqs->p[i][j+nbElimVars+1]); 972 for (j=0; j< permutedIneqs->NbColumns-nbParms-2-nbElimVars; j++) 973 value_assign(Full_Dim->p[i][j+1], 974 permutedIneqs->p[i][nbElimVars+nbParms+j+1]); 975 value_assign(Full_Dim->p[i][Full_Dim->NbColumns-1], 976 permutedIneqs->p[i][permutedIneqs->NbColumns-1]); 977 } 978 Matrix_Free(permutedIneqs); 979 980 /* 5- Keep only the the validity lattice restricted to the parameters */ 981 *validityLattice = Matrix_Alloc(nbParms+1, nbParms+1); 982 for (i=0; i< nbParms; i++) { 983 for (j=0; j< nbParms; j++) 984 value_assign((*validityLattice)->p[i][j], 985 WVL->p[i][j]); 986 value_assign((*validityLattice)->p[i][nbParms], 987 WVL->p[i][WVL->NbColumns-1]); 988 } 989 for (j=0; j< nbParms; j++) 990 value_set_si((*validityLattice)->p[nbParms][j], 0); 991 value_assign((*validityLattice)->p[nbParms][nbParms], 992 WVL->p[WVL->NbColumns-1][WVL->NbColumns-1]); 993 994 /* 6- Clean up */ 995 Matrix_Free(WVL); 996 return Full_Dim; 997 } /* full_dimensionize */ 998 999 #undef dbgCompParm 1000 #undef dbgCompParmMore 1001