Coverage Report

Created: 2017-11-21 16:49

/Users/buildslave/jenkins/workspace/clang-stage2-coverage-R/llvm/tools/polly/lib/External/isl/isl_equalities.c
 Line Count Source (jump to first uncovered line) 1 /* 2  * Copyright 2008-2009 Katholieke Universiteit Leuven 3  * Copyright 2010 INRIA Saclay 4  * 5  * Use of this software is governed by the MIT license 6  * 7  * Written by Sven Verdoolaege, K.U.Leuven, Departement 8  * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium 9  * and INRIA Saclay - Ile-de-France, Parc Club Orsay Universite, 10  * ZAC des vignes, 4 rue Jacques Monod, 91893 Orsay, France 11  */ 12 13 #include  14 #include  15 #include  16 #include "isl_map_private.h" 17 #include "isl_equalities.h" 18 #include  19 20 /* Given a set of modulo constraints 21  * 22  *    c + A y = 0 mod d 23  * 24  * this function computes a particular solution y_0 25  * 26  * The input is given as a matrix B = [ c A ] and a vector d. 27  * 28  * The output is matrix containing the solution y_0 or 29  * a zero-column matrix if the constraints admit no integer solution. 30  * 31  * The given set of constrains is equivalent to 32  * 33  *    c + A y = -D x 34  * 35  * with D = diag d and x a fresh set of variables. 36  * Reducing both c and A modulo d does not change the 37  * value of y in the solution and may lead to smaller coefficients. 38  * Let M = [ D A ] and [ H 0 ] = M U, the Hermite normal form of M. 39  * Then 40  *     [ x ] 41  *    M [ y ] = - c 42  * and so 43  *     [ x ] 44  *    [ H 0 ] U^{-1} [ y ] = - c 45  * Let 46  *    [ A ] [ x ] 47  *    [ B ] = U^{-1} [ y ] 48  * then 49  *    H A + 0 B = -c 50  * 51  * so B may be chosen arbitrarily, e.g., B = 0, and then 52  * 53  *     [ x ] = [ -c ] 54  *    U^{-1} [ y ] = [ 0 ] 55  * or 56  *    [ x ] [ -c ] 57  *    [ y ] = U [ 0 ] 58  * specifically, 59  * 60  *    y = U_{2,1} (-c) 61  * 62  * If any of the coordinates of this y are non-integer 63  * then the constraints admit no integer solution and 64  * a zero-column matrix is returned. 65  */ 66 static struct isl_mat *particular_solution(struct isl_mat *B, struct isl_vec *d) 67 6.55k { 68 6.55k   int i, j; 69 6.55k   struct isl_mat *M = NULL; 70 6.55k   struct isl_mat *C = NULL; 71 6.55k   struct isl_mat *U = NULL; 72 6.55k   struct isl_mat *H = NULL; 73 6.55k   struct isl_mat *cst = NULL; 74 6.55k   struct isl_mat *T = NULL; 75 6.55k 76 6.55k   M = isl_mat_alloc(B->ctx, B->n_row, B->n_row + B->n_col - 1); 77 6.55k   C = isl_mat_alloc(B->ctx, 1 + B->n_row, 1); 78 6.55k   if (!M || !C) 79 0     goto error; 80 6.55k   isl_int_set_si(C->row[0][0], 1); 81 14.8k   for (i = 0; i < B->n_row; ++i8.32k) { 82 8.32k     isl_seq_clr(M->row[i], B->n_row); 83 8.32k     isl_int_set(M->row[i][i], d->block.data[i]); 84 8.32k     isl_int_neg(C->row[1 + i][0], B->row[i][0]); 85 8.32k     isl_int_fdiv_r(C->row[1+i][0], C->row[1+i][0], M->row[i][i]); 86 46.7k     for (j = 0; j < B->n_col - 1; ++j38.4k) 87 38.4k       isl_int_fdiv_r(M->row[i][B->n_row + j], 88 8.32k           B->row[i][1 + j], M->row[i][i]); 89 8.32k   } 90 6.55k   M = isl_mat_left_hermite(M, 0, &U, NULL); 91 6.55k   if (!M || !U) 92 0     goto error; 93 6.55k   H = isl_mat_sub_alloc(M, 0, B->n_row, 0, B->n_row); 94 6.55k   H = isl_mat_lin_to_aff(H); 95 6.55k   C = isl_mat_inverse_product(H, C); 96 6.55k   if (!C) 97 0     goto error; 98 14.2k   for (i = 0; 6.55ki < B->n_row; ++i7.70k) { 99 8.30k     if (!isl_int_is_divisible_by(C->row[1+i][0], C->row[0][0])) 100 8.30k       break599; 101 7.70k     isl_int_divexact(C->row[1+i][0], C->row[1+i][0], C->row[0][0]); 102 7.70k   } 103 6.55k   if (i < B->n_row) 104 599     cst = isl_mat_alloc(B->ctx, B->n_row, 0); 105 5.96k   else 106 5.96k     cst = isl_mat_sub_alloc(C, 1, B->n_row, 0, 1); 107 6.55k   T = isl_mat_sub_alloc(U, B->n_row, B->n_col - 1, 0, B->n_row); 108 6.55k   cst = isl_mat_product(T, cst); 109 6.55k   isl_mat_free(M); 110 6.55k   isl_mat_free(C); 111 6.55k   isl_mat_free(U); 112 6.55k   return cst; 113 0 error: 114 0   isl_mat_free(M); 115 0   isl_mat_free(C); 116 0   isl_mat_free(U); 117 0   return NULL; 118 6.55k } 119 120 /* Compute and return the matrix 121  * 122  *    U_1^{-1} diag(d_1, 1, ..., 1) 123  * 124  * with U_1 the unimodular completion of the first (and only) row of B. 125  * The columns of this matrix generate the lattice that satisfies 126  * the single (linear) modulo constraint. 127  */ 128 static struct isl_mat *parameter_compression_1( 129       struct isl_mat *B, struct isl_vec *d) 130 4.65k { 131 4.65k   struct isl_mat *U; 132 4.65k 133 4.65k   U = isl_mat_alloc(B->ctx, B->n_col - 1, B->n_col - 1); 134 4.65k   if (!U) 135 0     return NULL; 136 4.65k   isl_seq_cpy(U->row[0], B->row[0] + 1, B->n_col - 1); 137 4.65k   U = isl_mat_unimodular_complete(U, 1); 138 4.65k   U = isl_mat_right_inverse(U); 139 4.65k   if (!U) 140 0     return NULL; 141 4.65k   isl_mat_col_mul(U, 0, d->block.data[0], 0); 142 4.65k   U = isl_mat_lin_to_aff(U); 143 4.65k   return U; 144 4.65k } 145 146 /* Compute a common lattice of solutions to the linear modulo 147  * constraints specified by B and d. 148  * See also the documentation of isl_mat_parameter_compression. 149  * We put the matrix 150  *  151  *    A = [ L_1^{-T} L_2^{-T} ... L_k^{-T} ] 152  * 153  * on a common denominator. This denominator D is the lcm of modulos d. 154  * Since L_i = U_i^{-1} diag(d_i, 1, ... 1), we have 155  * L_i^{-T} = U_i^T diag(d_i, 1, ... 1)^{-T} = U_i^T diag(1/d_i, 1, ..., 1). 156  * Putting this on the common denominator, we have 157  * D * L_i^{-T} = U_i^T diag(D/d_i, D, ..., D). 158  */ 159 static struct isl_mat *parameter_compression_multi( 160       struct isl_mat *B, struct isl_vec *d) 161 1.26k { 162 1.26k   int i, j, k; 163 1.26k   isl_int D; 164 1.26k   struct isl_mat *A = NULL, *U = NULL; 165 1.26k   struct isl_mat *T; 166 1.26k   unsigned size; 167 1.26k 168 1.26k   isl_int_init(D); 169 1.26k 170 1.26k   isl_vec_lcm(d, &D); 171 1.26k 172 1.26k   size = B->n_col - 1; 173 1.26k   A = isl_mat_alloc(B->ctx, size, B->n_row * size); 174 1.26k   U = isl_mat_alloc(B->ctx, size, size); 175 1.26k   if (!U || !A) 176 0     goto error; 177 3.93k   for (i = 0; 1.26ki < B->n_row; ++i2.67k) { 178 2.67k     isl_seq_cpy(U->row[0], B->row[i] + 1, size); 179 2.67k     U = isl_mat_unimodular_complete(U, 1); 180 2.67k     if (!U) 181 0       goto error; 182 2.67k     isl_int_divexact(D, D, d->block.data[i]); 183 17.0k     for (k = 0; k < U->n_col; ++k14.4k) 184 14.4k       isl_int_mul(A->row[k][i*size+0], D, U->row[0][k]); 185 2.67k     isl_int_mul(D, D, d->block.data[i]); 186 14.4k     for (j = 1; j < U->n_row; ++j11.7k) 187 90.4k       for (k = 0; 11.7kk < U->n_col; ++k78.7k) 188 78.7k         isl_int_mul(A->row[k][i*size+j], 189 2.67k             D, U->row[j][k]); 190 2.67k   } 191 1.26k   A = isl_mat_left_hermite(A, 0, NULL, NULL); 192 1.26k   T = isl_mat_sub_alloc(A, 0, A->n_row, 0, A->n_row); 193 1.26k   T = isl_mat_lin_to_aff(T); 194 1.26k   if (!T) 195 0     goto error; 196 1.26k   isl_int_set(T->row[0][0], D); 197 1.26k   T = isl_mat_right_inverse(T); 198 1.26k   if (!T) 199 0     goto error; 200 1.26k   isl_assert(T->ctx, isl_int_is_one(T->row[0][0]), goto error); 201 1.26k   T = isl_mat_transpose(T); 202 1.26k   isl_mat_free(A); 203 1.26k   isl_mat_free(U); 204 1.26k 205 1.26k   isl_int_clear(D); 206 1.26k   return T; 207 0 error: 208 0   isl_mat_free(A); 209 0   isl_mat_free(U); 210 0   isl_int_clear(D); 211 0   return NULL; 212 1.26k } 213 214 /* Given a set of modulo constraints 215  * 216  *    c + A y = 0 mod d 217  * 218  * this function returns an affine transformation T, 219  * 220  *    y = T y' 221  * 222  * that bijectively maps the integer vectors y' to integer 223  * vectors y that satisfy the modulo constraints. 224  * 225  * This function is inspired by Section 2.5.3 226  * of B. Meister, "Stating and Manipulating Periodicity in the Polytope 227  * Model. Applications to Program Analysis and Optimization". 228  * However, the implementation only follows the algorithm of that 229  * section for computing a particular solution and not for computing 230  * a general homogeneous solution. The latter is incomplete and 231  * may remove some valid solutions. 232  * Instead, we use an adaptation of the algorithm in Section 7 of 233  * B. Meister, S. Verdoolaege, "Polynomial Approximations in the Polytope 234  * Model: Bringing the Power of Quasi-Polynomials to the Masses". 235  * 236  * The input is given as a matrix B = [ c A ] and a vector d. 237  * Each element of the vector d corresponds to a row in B. 238  * The output is a lower triangular matrix. 239  * If no integer vector y satisfies the given constraints then 240  * a matrix with zero columns is returned. 241  * 242  * We first compute a particular solution y_0 to the given set of 243  * modulo constraints in particular_solution. If no such solution 244  * exists, then we return a zero-columned transformation matrix. 245  * Otherwise, we compute the generic solution to 246  * 247  *    A y = 0 mod d 248  * 249  * That is we want to compute G such that 250  * 251  *    y = G y'' 252  * 253  * with y'' integer, describes the set of solutions. 254  * 255  * We first remove the common factors of each row. 256  * In particular if gcd(A_i,d_i) != 1, then we divide the whole 257  * row i (including d_i) by this common factor. If afterwards gcd(A_i) != 1, 258  * then we divide this row of A by the common factor, unless gcd(A_i) = 0. 259  * In the later case, we simply drop the row (in both A and d). 260  * 261  * If there are no rows left in A, then G is the identity matrix. Otherwise, 262  * for each row i, we now determine the lattice of integer vectors 263  * that satisfies this row. Let U_i be the unimodular extension of the 264  * row A_i. This unimodular extension exists because gcd(A_i) = 1. 265  * The first component of 266  * 267  *    y' = U_i y 268  * 269  * needs to be a multiple of d_i. Let y' = diag(d_i, 1, ..., 1) y''. 270  * Then, 271  * 272  *    y = U_i^{-1} diag(d_i, 1, ..., 1) y'' 273  * 274  * for arbitrary integer vectors y''. That is, y belongs to the lattice 275  * generated by the columns of L_i = U_i^{-1} diag(d_i, 1, ..., 1). 276  * If there is only one row, then G = L_1. 277  * 278  * If there is more than one row left, we need to compute the intersection 279  * of the lattices. That is, we need to compute an L such that 280  * 281  *    L = L_i L_i'  for all i 282  * 283  * with L_i' some integer matrices. Let A be constructed as follows 284  * 285  *    A = [ L_1^{-T} L_2^{-T} ... L_k^{-T} ] 286  * 287  * and computed the Hermite Normal Form of A = [ H 0 ] U 288  * Then, 289  * 290  *    L_i^{-T} = H U_{1,i} 291  * 292  * or 293  * 294  *    H^{-T} = L_i U_{1,i}^T 295  * 296  * In other words G = L = H^{-T}. 297  * To ensure that G is lower triangular, we compute and use its Hermite 298  * normal form. 299  * 300  * The affine transformation matrix returned is then 301  * 302  *    [ 1 0 ] 303  *    [ y_0 G ] 304  * 305  * as any y = y_0 + G y' with y' integer is a solution to the original 306  * modulo constraints. 307  */ 308 __isl_give isl_mat *isl_mat_parameter_compression(__isl_take isl_mat *B, 309   __isl_take isl_vec *d) 310 6.55k { 311 6.55k   int i; 312 6.55k   struct isl_mat *cst = NULL; 313 6.55k   struct isl_mat *T = NULL; 314 6.55k   isl_int D; 315 6.55k 316 6.55k   if (!B || !d) 317 0     goto error; 318 6.55k   isl_assert(B->ctx, B->n_row == d->size, goto error); 319 6.55k   cst = particular_solution(B, d); 320 6.55k   if (!cst) 321 0     goto error; 322 6.55k   if (cst->n_col == 0) { 323 599     T = isl_mat_alloc(B->ctx, B->n_col, 0); 324 599     isl_mat_free(cst); 325 599     isl_mat_free(B); 326 599     isl_vec_free(d); 327 599     return T; 328 599   } 329 5.96k   isl_int_init(D); 330 5.96k   /* Replace a*g*row = 0 mod g*m by row = 0 mod m */ 331 13.3k   for (i = 0; i < B->n_row; ++i7.39k) { 332 7.39k     isl_seq_gcd(B->row[i] + 1, B->n_col - 1, &D); 333 7.39k     if (isl_int_is_one(D)) 334 7.39k       continue7.06k; 335 333     if (isl_int_is_zero(D)) { 336 68       B = isl_mat_drop_rows(B, i, 1); 337 68       d = isl_vec_cow(d); 338 68       if (!B || !d) 339 0         goto error2; 340 68       isl_seq_cpy(d->block.data+i, d->block.data+i+1, 341 68               d->size - (i+1)); 342 68       d->size--; 343 68       i--; 344 68       continue; 345 68     } 346 265     B = isl_mat_cow(B); 347 265     if (!B) 348 0       goto error2; 349 265     isl_seq_scale_down(B->row[i] + 1, B->row[i] + 1, D, B->n_col-1); 350 265     isl_int_gcd(D, D, d->block.data[i]); 351 265     d = isl_vec_cow(d); 352 265     if (!d) 353 0       goto error2; 354 265     isl_int_divexact(d->block.data[i], d->block.data[i], D); 355 265   } 356 5.96k   isl_int_clear(D); 357 5.96k   if (B->n_row == 0) 358 38     T = isl_mat_identity(B->ctx, B->n_col); 359 5.92k   else if (B->n_row == 1) 360 4.65k     T = parameter_compression_1(B, d); 361 1.26k   else 362 1.26k     T = parameter_compression_multi(B, d); 363 5.96k   T = isl_mat_left_hermite(T, 0, NULL, NULL); 364 5.96k   if (!T) 365 0     goto error; 366 5.96k   isl_mat_sub_copy(T->ctx, T->row + 1, cst->row, cst->n_row, 0, 0, 1); 367 5.96k   isl_mat_free(cst); 368 5.96k   isl_mat_free(B); 369 5.96k   isl_vec_free(d); 370 5.96k   return T; 371 0 error2: 372 0   isl_int_clear(D); 373 0 error: 374 0   isl_mat_free(cst); 375 0   isl_mat_free(B); 376 0   isl_vec_free(d); 377 0   return NULL; 378 6.55k } 379 380 /* Given a set of equalities 381  * 382  *    B(y) + A x = 0            (*) 383  * 384  * compute and return an affine transformation T, 385  * 386  *    y = T y' 387  * 388  * that bijectively maps the integer vectors y' to integer 389  * vectors y that satisfy the modulo constraints for some value of x. 390  * 391  * Let [H 0] be the Hermite Normal Form of A, i.e., 392  * 393  *    A = [H 0] Q 394  * 395  * Then y is a solution of (*) iff 396  * 397  *    H^-1 B(y) (= - [I 0] Q x) 398  * 399  * is an integer vector. Let d be the common denominator of H^-1. 400  * We impose 401  * 402  *    d H^-1 B(y) = 0 mod d 403  * 404  * and compute the solution using isl_mat_parameter_compression. 405  */ 406 __isl_give isl_mat *isl_mat_parameter_compression_ext(__isl_take isl_mat *B, 407   __isl_take isl_mat *A) 408 0 { 409 0   isl_ctx *ctx; 410 0   isl_vec *d; 411 0   int n_row, n_col; 412 0 413 0   if (!A) 414 0     return isl_mat_free(B); 415 0 416 0   ctx = isl_mat_get_ctx(A); 417 0   n_row = A->n_row; 418 0   n_col = A->n_col; 419 0   A = isl_mat_left_hermite(A, 0, NULL, NULL); 420 0   A = isl_mat_drop_cols(A, n_row, n_col - n_row); 421 0   A = isl_mat_lin_to_aff(A); 422 0   A = isl_mat_right_inverse(A); 423 0   d = isl_vec_alloc(ctx, n_row); 424 0   if (A) 425 0     d = isl_vec_set(d, A->row[0][0]); 426 0   A = isl_mat_drop_rows(A, 0, 1); 427 0   A = isl_mat_drop_cols(A, 0, 1); 428 0   B = isl_mat_product(A, B); 429 0 430 0   return isl_mat_parameter_compression(B, d); 431 0 } 432 433 /* Return a compression matrix that indicates that there are no solutions 434  * to the original constraints. In particular, return a zero-column 435  * matrix with 1 + dim rows. If "T2" is not NULL, then assign *T2 436  * the inverse of this matrix. *T2 may already have been assigned 437  * matrix, so free it first. 438  * "free1", "free2" and "free3" are temporary matrices that are 439  * not useful when an empty compression is returned. They are 440  * simply freed. 441  */ 442 static __isl_give isl_mat *empty_compression(isl_ctx *ctx, unsigned dim, 443   __isl_give isl_mat **T2, __isl_take isl_mat *free1, 444   __isl_take isl_mat *free2, __isl_take isl_mat *free3) 445 0 { 446 0   isl_mat_free(free1); 447 0   isl_mat_free(free2); 448 0   isl_mat_free(free3); 449 0   if (T2) { 450 0     isl_mat_free(*T2); 451 0     *T2 = isl_mat_alloc(ctx, 0, 1 + dim); 452 0   } 453 0   return isl_mat_alloc(ctx, 1 + dim, 0); 454 0 } 455 456 /* Given a matrix that maps a (possibly) parametric domain to 457  * a parametric domain, add in rows that map the "nparam" parameters onto 458  * themselves. 459  */ 460 static __isl_give isl_mat *insert_parameter_rows(__isl_take isl_mat *mat, 461   unsigned nparam) 462 60.6k { 463 60.6k   int i; 464 60.6k 465 60.6k   if (nparam == 0) 466 60.6k     return mat; 467 0   if (!mat) 468 0     return NULL; 469 0 470 0   mat = isl_mat_insert_rows(mat, 1, nparam); 471 0   if (!mat) 472 0     return NULL; 473 0 474 0   for (i = 0; i < nparam; ++i) { 475 0     isl_seq_clr(mat->row[1 + i], mat->n_col); 476 0     isl_int_set(mat->row[1 + i][1 + i], mat->row[0][0]); 477 0   } 478 60.6k 479 60.6k   return mat; 480 60.6k } 481 482 /* Given a set of equalities 483  * 484  *    -C(y) + M x = 0 485  * 486  * this function computes a unimodular transformation from a lower-dimensional 487  * space to the original space that bijectively maps the integer points x' 488  * in the lower-dimensional space to the integer points x in the original 489  * space that satisfy the equalities. 490  * 491  * The input is given as a matrix B = [ -C M ] and the output is a 492  * matrix that maps [1 x'] to [1 x]. 493  * The number of equality constraints in B is assumed to be smaller than 494  * or equal to the number of variables x. 495  * "first" is the position of the first x variable. 496  * The preceding variables are considered to be y-variables. 497  * If T2 is not NULL, then *T2 is set to a matrix mapping [1 x] to [1 x']. 498  * 499  * First compute the (left) Hermite normal form of M, 500  * 501  *    M [U1 U2] = M U = H = [H1 0] 502  * or 503  *     M = H Q = [H1 0] [Q1] 504  * [Q2] 505  * 506  * with U, Q unimodular, Q = U^{-1} (and H lower triangular). 507  * Define the transformed variables as 508  * 509  *    x = [U1 U2] [ x1' ] = [U1 U2] [Q1] x 510  *     [ x2' ] [Q2] 511  * 512  * The equalities then become 513  * 514  *    -C(y) + H1 x1' = 0 or x1' = H1^{-1} C(y) = C'(y) 515  * 516  * If the denominator of the constant term does not divide the 517  * the common denominator of the coefficients of y, then every 518  * integer point is mapped to a non-integer point and then the original set 519  * has no integer solutions (since the x' are a unimodular transformation 520  * of the x). In this case, a zero-column matrix is returned. 521  * Otherwise, the transformation is given by 522  * 523  *    x = U1 H1^{-1} C(y) + U2 x2' 524  * 525  * The inverse transformation is simply 526  * 527  *    x2' = Q2 x 528  */ 529 __isl_give isl_mat *isl_mat_final_variable_compression(__isl_take isl_mat *B, 530   int first, __isl_give isl_mat **T2) 531 60.6k { 532 60.6k   int i, n; 533 60.6k   isl_ctx *ctx; 534 60.6k   isl_mat *H = NULL, *C, *H1, *U = NULL, *U1, *U2; 535 60.6k   unsigned dim; 536 60.6k 537 60.6k   if (T2) 538 19.0k     *T2 = NULL; 539 60.6k   if (!B) 540 0     goto error; 541 60.6k 542 60.6k   ctx = isl_mat_get_ctx(B); 543 60.6k   dim = B->n_col - 1; 544 60.6k   n = dim - first; 545 60.6k   if (n < B->n_row) 546 60.6k     isl_die0(ctx, isl_error_invalid, "too many equality constraints", 547 60.6k       goto error); 548 60.6k   H = isl_mat_sub_alloc(B, 0, B->n_row, 1 + first, n); 549 60.6k   H = isl_mat_left_hermite(H, 0, &U, T2); 550 60.6k   if (!H || !U || (T2 && !*T219.0k)) 551 0     goto error; 552 60.6k   if (T2) { 553 19.0k     *T2 = isl_mat_drop_rows(*T2, 0, B->n_row); 554 19.0k     *T2 = isl_mat_diagonal(isl_mat_identity(ctx, 1 + first), *T2); 555 19.0k     if (!*T2) 556 0       goto error; 557 60.6k   } 558 60.6k   C = isl_mat_alloc(ctx, 1 + B->n_row, 1 + first); 559 60.6k   if (!C) 560 0     goto error; 561 60.6k   isl_int_set_si(C->row[0][0], 1); 562 60.6k   isl_seq_clr(C->row[0] + 1, first); 563 60.6k   isl_mat_sub_neg(ctx, C->row + 1, B->row, B->n_row, 0, 0, 1 + first); 564 60.6k   H1 = isl_mat_sub_alloc(H, 0, H->n_row, 0, H->n_row); 565 60.6k   H1 = isl_mat_lin_to_aff(H1); 566 60.6k   C = isl_mat_inverse_product(H1, C); 567 60.6k   if (!C) 568 0     goto error; 569 60.6k   isl_mat_free(H); 570 60.6k   if (!isl_int_is_one(C->row[0][0])) { 571 379     isl_int g; 572 379 573 379     isl_int_init(g); 574 2.11k     for (i = 0; i < B->n_row; ++i1.73k) { 575 1.73k       isl_seq_gcd(C->row[1 + i] + 1, first, &g); 576 1.73k       isl_int_gcd(g, g, C->row[0][0]); 577 1.73k       if (!isl_int_is_divisible_by(C->row[1 + i][0], g)) 578 1.73k         break0; 579 1.73k     } 580 379     isl_int_clear(g); 581 379 582 379     if (i < B->n_row) 583 0       return empty_compression(ctx, dim, T2, B, C, U); 584 379     C = isl_mat_normalize(C); 585 379   } 586 60.6k   U1 = isl_mat_sub_alloc(U, 0, U->n_row, 0, B->n_row); 587 60.6k   U1 = isl_mat_lin_to_aff(U1); 588 60.6k   U2 = isl_mat_sub_alloc(U, 0, U->n_row, B->n_row, U->n_row - B->n_row); 589 60.6k   U2 = isl_mat_lin_to_aff(U2); 590 60.6k   isl_mat_free(U); 591 60.6k   C = isl_mat_product(U1, C); 592 60.6k   C = isl_mat_aff_direct_sum(C, U2); 593 60.6k   C = insert_parameter_rows(C, first); 594 60.6k 595 60.6k   isl_mat_free(B); 596 60.6k 597 60.6k   return C; 598 0 error: 599 0   isl_mat_free(B); 600 0   isl_mat_free(H); 601 0   isl_mat_free(U); 602 0   if (T2) { 603 0     isl_mat_free(*T2); 604 0     *T2 = NULL; 605 0   } 606 0   return NULL; 607 60.6k } 608 609 /* Given a set of equalities 610  * 611  *    M x - c = 0 612  * 613  * this function computes a unimodular transformation from a lower-dimensional 614  * space to the original space that bijectively maps the integer points x' 615  * in the lower-dimensional space to the integer points x in the original 616  * space that satisfy the equalities. 617  * 618  * The input is given as a matrix B = [ -c M ] and the output is a 619  * matrix that maps [1 x'] to [1 x]. 620  * The number of equality constraints in B is assumed to be smaller than 621  * or equal to the number of variables x. 622  * If T2 is not NULL, then *T2 is set to a matrix mapping [1 x] to [1 x']. 623  */ 624 __isl_give isl_mat *isl_mat_variable_compression(__isl_take isl_mat *B, 625   __isl_give isl_mat **T2) 626 60.6k { 627 60.6k   return isl_mat_final_variable_compression(B, 0, T2); 628 60.6k } 629 630 /* Return "bset" and set *T and *T2 to the identity transformation 631  * on "bset" (provided T and T2 are not NULL). 632  */ 633 static __isl_give isl_basic_set *return_with_identity( 634   __isl_take isl_basic_set *bset, __isl_give isl_mat **T, 635   __isl_give isl_mat **T2) 636 0 { 637 0   unsigned dim; 638 0   isl_mat *id; 639 0 640 0   if (!bset) 641 0     return NULL; 642 0   if (!T && !T2) 643 0     return bset; 644 0 645 0   dim = isl_basic_set_dim(bset, isl_dim_set); 646 0   id = isl_mat_identity(isl_basic_map_get_ctx(bset), 1 + dim); 647 0   if (T) 648 0     *T = isl_mat_copy(id); 649 0   if (T2) 650 0     *T2 = isl_mat_copy(id); 651 0   isl_mat_free(id); 652 0 653 0   return bset; 654 0 } 655 656 /* Use the n equalities of bset to unimodularly transform the 657  * variables x such that n transformed variables x1' have a constant value 658  * and rewrite the constraints of bset in terms of the remaining 659  * transformed variables x2'. The matrix pointed to by T maps 660  * the new variables x2' back to the original variables x, while T2 661  * maps the original variables to the new variables. 662  */ 663 static struct isl_basic_set *compress_variables( 664   struct isl_basic_set *bset, struct isl_mat **T, struct isl_mat **T2) 665 52.2k { 666 52.2k   struct isl_mat *B, *TC; 667 52.2k   unsigned dim; 668 52.2k 669 52.2k   if (T) 670 52.2k     *T = NULL; 671 52.2k   if (T2) 672 14.7k     *T2 = NULL; 673 52.2k   if (!bset) 674 0     goto error; 675 52.2k   isl_assert(bset->ctx, isl_basic_set_n_param(bset) == 0, goto error); 676 52.2k   isl_assert(bset->ctx, bset->n_div == 0, goto error); 677 52.2k   dim = isl_basic_set_n_dim(bset); 678 52.2k   isl_assert(bset->ctx, bset->n_eq <= dim, goto error); 679 52.2k   if (bset->n_eq == 0) 680 0     return return_with_identity(bset, T, T2); 681 52.2k 682 52.2k   B = isl_mat_sub_alloc6(bset->ctx, bset->eq, 0, bset->n_eq, 0, 1 + dim); 683 52.2k   TC = isl_mat_variable_compression(B, T2); 684 52.2k   if (!TC) 685 0     goto error; 686 52.2k   if (TC->n_col == 0) { 687 0     isl_mat_free(TC); 688 0     if (T2) { 689 0       isl_mat_free(*T2); 690 0       *T2 = NULL; 691 0     } 692 0     bset = isl_basic_set_set_to_empty(bset); 693 0     return return_with_identity(bset, T, T2); 694 0   } 695 52.2k 696 52.2k   bset = isl_basic_set_preimage(bset, T ? isl_mat_copy(TC) : TC0); 697 52.2k   if (T) 698 52.2k     *T = TC; 699 52.2k   return bset; 700 0 error: 701 0   isl_basic_set_free(bset); 702 0   return NULL; 703 52.2k } 704 705 struct isl_basic_set *isl_basic_set_remove_equalities( 706   struct isl_basic_set *bset, struct isl_mat **T, struct isl_mat **T2) 707 52.2k { 708 52.2k   if (T) 709 52.2k     *T = NULL; 710 52.2k   if (T2) 711 14.7k     *T2 = NULL; 712 52.2k   if (!bset) 713 0     return NULL; 714 52.2k   isl_assert(bset->ctx, isl_basic_set_n_param(bset) == 0, goto error); 715 52.2k   bset = isl_basic_set_gauss(bset, NULL); 716 52.2k   if (ISL_F_ISSET(bset, ISL_BASIC_SET_EMPTY)) 717 52.2k     return return_with_identity(bset, T, T2)0; 718 52.2k   bset = compress_variables(bset, T, T2); 719 52.2k   return bset; 720 0 error: 721 0   isl_basic_set_free(bset); 722 0   *T = NULL; 723 0   return NULL; 724 52.2k } 725 726 /* Check if dimension dim belongs to a residue class 727  *    i_dim \equiv r mod m 728  * with m != 1 and if so return m in *modulo and r in *residue. 729  * As a special case, when i_dim has a fixed value v, then 730  * *modulo is set to 0 and *residue to v. 731  * 732  * If i_dim does not belong to such a residue class, then *modulo 733  * is set to 1 and *residue is set to 0. 734  */ 735 isl_stat isl_basic_set_dim_residue_class(__isl_keep isl_basic_set *bset, 736   int pos, isl_int *modulo, isl_int *residue) 737 21 { 738 21   isl_bool fixed; 739 21   struct isl_ctx *ctx; 740 21   struct isl_mat *H = NULL, *U = NULL, *C, *H1, *U1; 741 21   unsigned total; 742 21   unsigned nparam; 743 21 744 21   if (!bset || !modulo || !residue) 745 0     return isl_stat_error; 746 21 747 21   fixed = isl_basic_set_plain_dim_is_fixed(bset, pos, residue); 748 21   if (fixed < 0) 749 0     return isl_stat_error; 750 21   if (fixed) { 751 0     isl_int_set_si(*modulo, 0); 752 0     return isl_stat_ok; 753 0   } 754 21 755 21   ctx = isl_basic_set_get_ctx(bset); 756 21   total = isl_basic_set_total_dim(bset); 757 21   nparam = isl_basic_set_n_param(bset); 758 21   H = isl_mat_sub_alloc6(ctx, bset->eq, 0, bset->n_eq, 1, total); 759 21   H = isl_mat_left_hermite(H, 0, &U, NULL); 760 21   if (!H) 761 0     return isl_stat_error; 762 21 763 21   isl_seq_gcd(U->row[nparam + pos]+bset->n_eq, 764 21       total-bset->n_eq, modulo); 765 21   if (isl_int_is_zero(*modulo)) 766 21     isl_int_set_si0(*modulo, 1); 767 21   if (isl_int_is_one(*modulo)) { 768 16     isl_int_set_si(*residue, 0); 769 16     isl_mat_free(H); 770 16     isl_mat_free(U); 771 16     return isl_stat_ok; 772 16   } 773 5 774 5   C = isl_mat_alloc(ctx, 1 + bset->n_eq, 1); 775 5   if (!C) 776 0     goto error; 777 5   isl_int_set_si(C->row[0][0], 1); 778 5   isl_mat_sub_neg(ctx, C->row + 1, bset->eq, bset->n_eq, 0, 0, 1); 779 5   H1 = isl_mat_sub_alloc(H, 0, H->n_row, 0, H->n_row); 780 5   H1 = isl_mat_lin_to_aff(H1); 781 5   C = isl_mat_inverse_product(H1, C); 782 5   isl_mat_free(H); 783 5   U1 = isl_mat_sub_alloc(U, nparam+pos, 1, 0, bset->n_eq); 784 5   U1 = isl_mat_lin_to_aff(U1); 785 5   isl_mat_free(U); 786 5   C = isl_mat_product(U1, C); 787 5   if (!C) 788 0     return isl_stat_error; 789 5   if (!isl_int_is_divisible_by(C->row[1][0], C->row[0][0])) { 790 0     bset = isl_basic_set_copy(bset); 791 0     bset = isl_basic_set_set_to_empty(bset); 792 0     isl_basic_set_free(bset); 793 0     isl_int_set_si(*modulo, 1); 794 0     isl_int_set_si(*residue, 0); 795 0     return isl_stat_ok; 796 0   } 797 5   isl_int_divexact(*residue, C->row[1][0], C->row[0][0]); 798 5   isl_int_fdiv_r(*residue, *residue, *modulo); 799 5   isl_mat_free(C); 800 5   return isl_stat_ok; 801 0 error: 802 0   isl_mat_free(H); 803 0   isl_mat_free(U); 804 0   return isl_stat_error; 805 21 } 806 807 /* Check if dimension dim belongs to a residue class 808  *    i_dim \equiv r mod m 809  * with m != 1 and if so return m in *modulo and r in *residue. 810  * As a special case, when i_dim has a fixed value v, then 811  * *modulo is set to 0 and *residue to v. 812  * 813  * If i_dim does not belong to such a residue class, then *modulo 814  * is set to 1 and *residue is set to 0. 815  */ 816 isl_stat isl_set_dim_residue_class(__isl_keep isl_set *set, 817   int pos, isl_int *modulo, isl_int *residue) 818 21 { 819 21   isl_int m; 820 21   isl_int r; 821 21   int i; 822 21 823 21   if (!set || !modulo || !residue) 824 0     return isl_stat_error; 825 21 826 21   if (set->n == 0) { 827 0     isl_int_set_si(*modulo, 0); 828 0     isl_int_set_si(*residue, 0); 829 0     return isl_stat_ok; 830 0   } 831 21 832 21   if (isl_basic_set_dim_residue_class(set->p[0], pos, modulo, residue)<0) 833 0     return isl_stat_error; 834 21 835 21   if (set->n == 1) 836 18     return isl_stat_ok; 837 3 838 3   if (isl_int_is_one(*modulo)) 839 3     return isl_stat_ok; 840 0 841 0   isl_int_init(m); 842 0   isl_int_init(r); 843 0 844 0   for (i = 1; i < set->n; ++i) { 845 0     if (isl_basic_set_dim_residue_class(set->p[i], pos, &m, &r) < 0) 846 0       goto error; 847 0     isl_int_gcd(*modulo, *modulo, m); 848 0     isl_int_sub(m, *residue, r); 849 0     isl_int_gcd(*modulo, *modulo, m); 850 0     if (!isl_int_is_zero(*modulo)) 851 0       isl_int_fdiv_r(*residue, *residue, *modulo); 852 0     if (isl_int_is_one(*modulo)) 853 0       break; 854 0   } 855 0 856 0   isl_int_clear(m); 857 0   isl_int_clear(r); 858 0 859 0   return isl_stat_ok; 860 0 error: 861 0   isl_int_clear(m); 862 0   isl_int_clear(r); 863 0   return isl_stat_error; 864 21 } 865 866 /* Check if dimension "dim" belongs to a residue class 867  *    i_dim \equiv r mod m 868  * with m != 1 and if so return m in *modulo and r in *residue. 869  * As a special case, when i_dim has a fixed value v, then 870  * *modulo is set to 0 and *residue to v. 871  * 872  * If i_dim does not belong to such a residue class, then *modulo 873  * is set to 1 and *residue is set to 0. 874  */ 875 isl_stat isl_set_dim_residue_class_val(__isl_keep isl_set *set, 876   int pos, __isl_give isl_val **modulo, __isl_give isl_val **residue) 877 21 { 878 21   *modulo = NULL; 879 21   *residue = NULL; 880 21   if (!set) 881 0     return isl_stat_error; 882 21   *modulo = isl_val_alloc(isl_set_get_ctx(set)); 883 21   *residue = isl_val_alloc(isl_set_get_ctx(set)); 884 21   if (!*modulo || !*residue) 885 0     goto error; 886 21   if (isl_set_dim_residue_class(set, pos, 887 21           &(*modulo)->n, &(*residue)->n) < 0) 888 0     goto error; 889 21   isl_int_set_si((*modulo)->d, 1); 890 21   isl_int_set_si((*residue)->d, 1); 891 21   return isl_stat_ok; 892 0 error: 893 0   isl_val_free(*modulo); 894 0   isl_val_free(*residue); 895 0   return isl_stat_error; 896 21 }