isl_equalities.c

   1 #include "isl_mat.h"
   2 #include "isl_map_private.h"
   3 #include "isl_equalities.h"
   4
   5 /* Use the n equalities of bset to unimodularly transform the
   6  * variables x such that n transformed variables x1' have a constant value
   7  * and rewrite the constraints of bset in terms of the remaining
   8  * transformed variables x2'.  The matrix pointed to by T maps
   9  * the new variables x2' back to the original variables x, while T2
  10  * maps the original variables to the new variables.
  11  *
  12  * Let the equalities of bset be
  13  *
  14  *              M x - c = 0
  15  *
  16  * Compute the (left) Hermite normal form of M,
  17  *
  18  *              M [U1 U2] = M U = H = [H1 0]
  19  * or
  20  *                            M = H Q = [H1 0] [Q1]
  21  *                                             [Q2]
  22  *
  23  * with U, Q unimodular, Q = U^{-1} (and H lower triangular).
  24  * Define the transformed variables as
  25  *
  26  *              x = [U1 U2] [ x1' ] = [U1 U2] [Q1] x
  27  *                          [ x2' ]           [Q2]
  28  *
  29  * The equalities then become
  30  *
  31  *              H1 x1' - c = 0   or   x1' = H1^{-1} c = c'
  32  *
  33  * If any of the c' is non-integer, then the original set has no
  34  * integer solutions (since the x' are a unimodular transformation
  35  * of the x).
  36  * Otherwise, the transformation is given by
  37  *
  38  *              x = U1 H1^{-1} c + U2 x2'
  39  *
  40  * The inverse transformation is simply
  41  *
  42  *              x2' = Q2 x
  43  */
  44 static struct isl_basic_set *compress_variables(struct isl_ctx *ctx,
  45         struct isl_basic_set *bset, struct isl_mat **T, struct isl_mat **T2)
  46 {
  47         int i;
  48         struct isl_mat *H = NULL, *C = NULL, *H1, *U = NULL, *U1, *U2, *TC;
  49
  50         if (T)
  51                 *T = NULL;
  52         if (T2)
  53                 *T2 = NULL;
  54         if (!bset)
  55                 goto error;
  56         isl_assert(ctx, bset->nparam == 0, goto error);
  57         isl_assert(ctx, bset->n_div == 0, goto error);
  58         isl_assert(ctx, bset->n_eq <= bset->dim, goto error);
  59         if (bset->n_eq == 0)
  60                 return bset;
  61
  62         H = isl_mat_sub_alloc(ctx, bset->eq, 0, bset->n_eq, 1, bset->dim);
  63         H = isl_mat_left_hermite(ctx, H, 0, &U, T2);
  64         if (!H || !U || (T2 && !*T2))
  65                 goto error;
  66         if (T2) {
  67                 *T2 = isl_mat_drop_rows(ctx, *T2, 0, bset->n_eq);
  68                 *T2 = isl_mat_lin_to_aff(ctx, *T2);
  69                 if (!*T2)
  70                         goto error;
  71         }
  72         C = isl_mat_alloc(ctx, 1+bset->n_eq, 1);
  73         if (!C)
  74                 goto error;
  75         isl_int_set_si(C->row[0][0], 1);
  76         isl_mat_sub_neg(ctx, C->row+1, bset->eq, bset->n_eq, 0, 0, 1);
  77         H1 = isl_mat_sub_alloc(ctx, H->row, 0, H->n_row, 0, H->n_row);
  78         H1 = isl_mat_lin_to_aff(ctx, H1);
  79         TC = isl_mat_inverse_product(ctx, H1, C);
  80         if (!TC)
  81                 goto error;
  82         isl_mat_free(ctx, H);
  83         if (!isl_int_is_one(TC->row[0][0])) {
  84                 for (i = 0; i < bset->n_eq; ++i) {
  85                         if (!isl_int_is_divisible_by(TC->row[1+i][0], TC->row[0][0])) {
  86                                 isl_mat_free(ctx, TC);
  87                                 isl_mat_free(ctx, U);
  88                                 if (T2) {
  89                                         isl_mat_free(ctx, *T2);
  90                                         *T2 = NULL;
  91                                 }
  92                                 return isl_basic_set_set_to_empty(bset);
  93                         }
  94                         isl_seq_scale_down(TC->row[1+i], TC->row[1+i], TC->row[0][0], 1);
  95                 }
  96                 isl_int_set_si(TC->row[0][0], 1);
  97         }
  98         U1 = isl_mat_sub_alloc(ctx, U->row, 0, U->n_row, 0, bset->n_eq);
  99         U1 = isl_mat_lin_to_aff(ctx, U1);
 100         U2 = isl_mat_sub_alloc(ctx, U->row, 0, U->n_row,
 101                                 bset->n_eq, U->n_row - bset->n_eq);
 102         U2 = isl_mat_lin_to_aff(ctx, U2);
 103         isl_mat_free(ctx, U);
 104         TC = isl_mat_product(ctx, U1, TC);
 105         TC = isl_mat_aff_direct_sum(ctx, TC, U2);
 106         bset = isl_basic_set_preimage(ctx, bset, T ? isl_mat_copy(ctx, TC) : TC);
 107         if (T)
 108                 *T = TC;
 109         return bset;
 110 error:
 111         isl_mat_free(ctx, H);
 112         isl_mat_free(ctx, U);
 113         if (T2)
 114                 isl_mat_free(ctx, *T2);
 115         isl_basic_set_free(bset);
 116         if (T)
 117                 *T = NULL;
 118         if (T2)
 119                 *T2 = NULL;
 120         return NULL;
 121 }
 122
 123 struct isl_basic_set *isl_basic_set_remove_equalities(
 124         struct isl_basic_set *bset, struct isl_mat **T, struct isl_mat **T2)
 125 {
 126         if (T)
 127                 *T = NULL;
 128         if (T2)
 129                 *T2 = NULL;
 130         if (!bset)
 131                 return NULL;
 132         isl_assert(bset->ctx, bset->nparam == 0, goto error);
 133         bset = isl_basic_set_gauss(bset, NULL);
 134         if (F_ISSET(bset, ISL_BASIC_SET_EMPTY))
 135                 return bset;
 136         bset = compress_variables(bset->ctx, bset, T, T2);
 137         return bset;
 138 error:
 139         isl_basic_set_free(bset);
 140         *T = NULL;
 141         return NULL;
 142 }