/* * Copyright 2008-2009 Katholieke Universiteit Leuven * Copyright 2010 INRIA Saclay * Copyright 2016-2017 Sven Verdoolaege * * Use of this software is governed by the MIT license * * Written by Sven Verdoolaege, K.U.Leuven, Departement * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium * and INRIA Saclay - Ile-de-France, Parc Club Orsay Universite, * ZAC des vignes, 4 rue Jacques Monod, 91893 Orsay, France */ #include <isl_ctx_private.h> #include "isl_map_private.h" #include <isl_seq.h> #include "isl_tab.h" #include "isl_sample.h" #include <isl_mat_private.h> #include <isl_vec_private.h> #include <isl_aff_private.h> #include <isl_constraint_private.h> #include <isl_options_private.h> #include <isl_config.h> #include <bset_to_bmap.c> /* * The implementation of parametric integer linear programming in this file * was inspired by the paper "Parametric Integer Programming" and the * report "Solving systems of affine (in)equalities" by Paul Feautrier * (and others). * * The strategy used for obtaining a feasible solution is different * from the one used in isl_tab.c. In particular, in isl_tab.c, * upon finding a constraint that is not yet satisfied, we pivot * in a row that increases the constant term of the row holding the * constraint, making sure the sample solution remains feasible * for all the constraints it already satisfied. * Here, we always pivot in the row holding the constraint, * choosing a column that induces the lexicographically smallest * increment to the sample solution. * * By starting out from a sample value that is lexicographically * smaller than any integer point in the problem space, the first * feasible integer sample point we find will also be the lexicographically * smallest. If all variables can be assumed to be non-negative, * then the initial sample value may be chosen equal to zero. * However, we will not make this assumption. Instead, we apply * the "big parameter" trick. Any variable x is then not directly * used in the tableau, but instead it is represented by another * variable x' = M + x, where M is an arbitrarily large (positive) * value. x' is therefore always non-negative, whatever the value of x. * Taking as initial sample value x' = 0 corresponds to x = -M, * which is always smaller than any possible value of x. * * The big parameter trick is used in the main tableau and * also in the context tableau if isl_context_lex is used. * In this case, each tableaus has its own big parameter. * Before doing any real work, we check if all the parameters * happen to be non-negative. If so, we drop the column corresponding * to M from the initial context tableau. * If isl_context_gbr is used, then the big parameter trick is only * used in the main tableau. */ struct isl_context; struct isl_context_op { … }; /* Shared parts of context representation. * * "n_unknown" is the number of final unknown integer divisions * in the input domain. */ struct isl_context { … }; struct isl_context_lex { … }; /* A stack (linked list) of solutions of subtrees of the search space. * * "ma" describes the solution as a function of "dom". * In particular, the domain space of "ma" is equal to the space of "dom". * * If "ma" is NULL, then there is no solution on "dom". */ struct isl_partial_sol { … }; struct isl_sol; struct isl_sol_callback { … }; /* isl_sol is an interface for constructing a solution to * a parametric integer linear programming problem. * Every time the algorithm reaches a state where a solution * can be read off from the tableau, the function "add" is called * on the isl_sol passed to find_solutions_main. In a state where * the tableau is empty, "add_empty" is called instead. * "free" is called to free the implementation specific fields, if any. * * "error" is set if some error has occurred. This flag invalidates * the remainder of the data structure. * If "rational" is set, then a rational optimization is being performed. * "level" is the current level in the tree with nodes for each * split in the context. * If "max" is set, then a maximization problem is being solved, rather than * a minimization problem, which means that the variables in the * tableau have value "M - x" rather than "M + x". * "n_out" is the number of output dimensions in the input. * "space" is the space in which the solution (and also the input) lives. * * The context tableau is owned by isl_sol and is updated incrementally. * * There are currently two implementations of this interface, * isl_sol_map, which simply collects the solutions in an isl_map * and (optionally) the parts of the context where there is no solution * in an isl_set, and * isl_sol_pma, which collects an isl_pw_multi_aff instead. */ struct isl_sol { … }; static void sol_free(struct isl_sol *sol) { … } /* Push a partial solution represented by a domain and function "ma" * onto the stack of partial solutions. * If "ma" is NULL, then "dom" represents a part of the domain * with no solution. */ static void sol_push_sol(struct isl_sol *sol, __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma) { … } /* Check that the final columns of "M", starting at "first", are zero. */ static isl_stat check_final_columns_are_zero(__isl_keep isl_mat *M, unsigned first) { … } /* Set the affine expressions in "ma" according to the rows in "M", which * are defined over the local space "ls". * The matrix "M" may have extra (zero) columns beyond the number * of variables in "ls". */ static __isl_give isl_multi_aff *set_from_affine_matrix( __isl_take isl_multi_aff *ma, __isl_take isl_local_space *ls, __isl_take isl_mat *M) { … } /* Push a partial solution represented by a domain and mapping M * onto the stack of partial solutions. * * The affine matrix "M" maps the dimensions of the context * to the output variables. Convert it into an isl_multi_aff and * then call sol_push_sol. * * Note that the description of the initial context may have involved * existentially quantified variables, in which case they also appear * in "dom". These need to be removed before creating the affine * expression because an affine expression cannot be defined in terms * of existentially quantified variables without a known representation. * Since newly added integer divisions are inserted before these * existentially quantified variables, they are still in the final * positions and the corresponding final columns of "M" are zero * because align_context_divs adds the existentially quantified * variables of the context to the main tableau without any constraints and * any equality constraints that are added later on can only serve * to eliminate these existentially quantified variables. */ static void sol_push_sol_mat(struct isl_sol *sol, __isl_take isl_basic_set *dom, __isl_take isl_mat *M) { … } /* Pop one partial solution from the partial solution stack and * pass it on to sol->add or sol->add_empty. */ static void sol_pop_one(struct isl_sol *sol) { … } /* Return a fresh copy of the domain represented by the context tableau. */ static struct isl_basic_set *sol_domain(struct isl_sol *sol) { … } /* Check whether two partial solutions have the same affine expressions. */ static isl_bool same_solution(struct isl_partial_sol *s1, struct isl_partial_sol *s2) { … } /* Swap the initial two partial solutions in "sol". * * That is, go from * * sol->partial = p1; p1->next = p2; p2->next = p3 * * to * * sol->partial = p2; p2->next = p1; p1->next = p3 */ static void swap_initial(struct isl_sol *sol) { … } /* Combine the initial two partial solution of "sol" into * a partial solution with the current context domain of "sol" and * the function description of the second partial solution in the list. * The level of the new partial solution is set to the current level. * * That is, the first two partial solutions (D1,M1) and (D2,M2) are * replaced by (D,M2), where D is the domain of "sol", which is assumed * to be the union of D1 and D2, while M1 is assumed to be equal to M2 * (at least on D1). */ static isl_stat combine_initial_into_second(struct isl_sol *sol) { … } /* Are "ma1" and "ma2" equal to each other on "dom"? * * Combine "ma1" and "ma2" with "dom" and check if the results are the same. * "dom" may have existentially quantified variables. Eliminate them first * as otherwise they would have to be eliminated twice, in a more complicated * context. */ static isl_bool equal_on_domain(__isl_keep isl_multi_aff *ma1, __isl_keep isl_multi_aff *ma2, __isl_keep isl_basic_set *dom) { … } /* The initial two partial solutions of "sol" are known to be at * the same level. * If they represent the same solution (on different parts of the domain), * then combine them into a single solution at the current level. * Otherwise, pop them both. * * Even if the two partial solution are not obviously the same, * one may still be a simplification of the other over its own domain. * Also check if the two sets of affine functions are equal when * restricted to one of the domains. If so, combine the two * using the set of affine functions on the other domain. * That is, for two partial solutions (D1,M1) and (D2,M2), * if M1 = M2 on D1, then the pair of partial solutions can * be replaced by (D1+D2,M2) and similarly when M1 = M2 on D2. */ static isl_stat combine_initial_if_equal(struct isl_sol *sol) { … } /* Pop all solutions from the partial solution stack that were pushed onto * the stack at levels that are deeper than the current level. * If the two topmost elements on the stack have the same level * and represent the same solution, then their domains are combined. * This combined domain is the same as the current context domain * as sol_pop is called each time we move back to a higher level. * If the outer level (0) has been reached, then all partial solutions * at the current level are also popped off. */ static void sol_pop(struct isl_sol *sol) { … } static void sol_dec_level(struct isl_sol *sol) { … } static isl_stat sol_dec_level_wrap(struct isl_tab_callback *cb) { … } /* Move down to next level and push callback onto context tableau * to decrease the level again when it gets rolled back across * the current state. That is, dec_level will be called with * the context tableau in the same state as it is when inc_level * is called. */ static void sol_inc_level(struct isl_sol *sol) { … } static void scale_rows(struct isl_mat *mat, isl_int m, int n_row) { … } /* Add the solution identified by the tableau and the context tableau. * * The layout of the variables is as follows. * tab->n_var is equal to the total number of variables in the input * map (including divs that were copied from the context) * + the number of extra divs constructed * Of these, the first tab->n_param and the last tab->n_div variables * correspond to the variables in the context, i.e., * tab->n_param + tab->n_div = context_tab->n_var * tab->n_param is equal to the number of parameters and input * dimensions in the input map * tab->n_div is equal to the number of divs in the context * * If there is no solution, then call add_empty with a basic set * that corresponds to the context tableau. (If add_empty is NULL, * then do nothing). * * If there is a solution, then first construct a matrix that maps * all dimensions of the context to the output variables, i.e., * the output dimensions in the input map. * The divs in the input map (if any) that do not correspond to any * div in the context do not appear in the solution. * The algorithm will make sure that they have an integer value, * but these values themselves are of no interest. * We have to be careful not to drop or rearrange any divs in the * context because that would change the meaning of the matrix. * * To extract the value of the output variables, it should be noted * that we always use a big parameter M in the main tableau and so * the variable stored in this tableau is not an output variable x itself, but * x' = M + x (in case of minimization) * or * x' = M - x (in case of maximization) * If x' appears in a column, then its optimal value is zero, * which means that the optimal value of x is an unbounded number * (-M for minimization and M for maximization). * We currently assume that the output dimensions in the original map * are bounded, so this cannot occur. * Similarly, when x' appears in a row, then the coefficient of M in that * row is necessarily 1. * If the row in the tableau represents * d x' = c + d M + e(y) * then, in case of minimization, the corresponding row in the matrix * will be * a c + a e(y) * with a d = m, the (updated) common denominator of the matrix. * In case of maximization, the row will be * -a c - a e(y) */ static void sol_add(struct isl_sol *sol, struct isl_tab *tab) { … } struct isl_sol_map { … }; static void sol_map_free(struct isl_sol *sol) { … } /* This function is called for parts of the context where there is * no solution, with "bset" corresponding to the context tableau. * Simply add the basic set to the set "empty". */ static void sol_map_add_empty(struct isl_sol_map *sol, struct isl_basic_set *bset) { … } static void sol_map_add_empty_wrap(struct isl_sol *sol, struct isl_basic_set *bset) { … } /* Given a basic set "dom" that represents the context and a tuple of * affine expressions "ma" defined over this domain, construct a basic map * that expresses this function on the domain. */ static void sol_map_add(struct isl_sol_map *sol, __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma) { … } static void sol_map_add_wrap(struct isl_sol *sol, __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma) { … } /* Store the "parametric constant" of row "row" of tableau "tab" in "line", * i.e., the constant term and the coefficients of all variables that * appear in the context tableau. * Note that the coefficient of the big parameter M is NOT copied. * The context tableau may not have a big parameter and even when it * does, it is a different big parameter. */ static void get_row_parameter_line(struct isl_tab *tab, int row, isl_int *line) { … } /* Check if rows "row1" and "row2" have identical "parametric constants", * as explained above. * In this case, we also insist that the coefficients of the big parameter * be the same as the values of the constants will only be the same * if these coefficients are also the same. */ static int identical_parameter_line(struct isl_tab *tab, int row1, int row2) { … } /* Return an inequality that expresses that the "parametric constant" * should be non-negative. * This function is only called when the coefficient of the big parameter * is equal to zero. */ static struct isl_vec *get_row_parameter_ineq(struct isl_tab *tab, int row) { … } /* Normalize a div expression of the form * * [(g*f(x) + c)/(g * m)] * * with c the constant term and f(x) the remaining coefficients, to * * [(f(x) + [c/g])/m] */ static void normalize_div(__isl_keep isl_vec *div) { … } /* Return an integer division for use in a parametric cut based * on the given row. * In particular, let the parametric constant of the row be * * \sum_i a_i y_i * * where y_0 = 1, but none of the y_i corresponds to the big parameter M. * The div returned is equal to * * floor(\sum_i {-a_i} y_i) = floor((\sum_i (-a_i mod d) y_i)/d) */ static struct isl_vec *get_row_parameter_div(struct isl_tab *tab, int row) { … } /* Return an integer division for use in transferring an integrality constraint * to the context. * In particular, let the parametric constant of the row be * * \sum_i a_i y_i * * where y_0 = 1, but none of the y_i corresponds to the big parameter M. * The the returned div is equal to * * floor(\sum_i {a_i} y_i) = floor((\sum_i (a_i mod d) y_i)/d) */ static struct isl_vec *get_row_split_div(struct isl_tab *tab, int row) { … } /* Construct and return an inequality that expresses an upper bound * on the given div. * In particular, if the div is given by * * d = floor(e/m) * * then the inequality expresses * * m d <= e */ static __isl_give isl_vec *ineq_for_div(__isl_keep isl_basic_set *bset, unsigned div) { … } /* Given a row in the tableau and a div that was created * using get_row_split_div and that has been constrained to equality, i.e., * * d = floor(\sum_i {a_i} y_i) = \sum_i {a_i} y_i * * replace the expression "\sum_i {a_i} y_i" in the row by d, * i.e., we subtract "\sum_i {a_i} y_i" and add 1 d. * The coefficients of the non-parameters in the tableau have been * verified to be integral. We can therefore simply replace coefficient b * by floor(b). For the coefficients of the parameters we have * floor(a_i) = a_i - {a_i}, while for the other coefficients, we have * floor(b) = b. */ static struct isl_tab *set_row_cst_to_div(struct isl_tab *tab, int row, int div) { … } /* Check if the (parametric) constant of the given row is obviously * negative, meaning that we don't need to consult the context tableau. * If there is a big parameter and its coefficient is non-zero, * then this coefficient determines the outcome. * Otherwise, we check whether the constant is negative and * all non-zero coefficients of parameters are negative and * belong to non-negative parameters. */ static int is_obviously_neg(struct isl_tab *tab, int row) { … } /* Check if the (parametric) constant of the given row is obviously * non-negative, meaning that we don't need to consult the context tableau. * If there is a big parameter and its coefficient is non-zero, * then this coefficient determines the outcome. * Otherwise, we check whether the constant is non-negative and * all non-zero coefficients of parameters are positive and * belong to non-negative parameters. */ static int is_obviously_nonneg(struct isl_tab *tab, int row) { … } /* Given a row r and two columns, return the column that would * lead to the lexicographically smallest increment in the sample * solution when leaving the basis in favor of the row. * Pivoting with column c will increment the sample value by a non-negative * constant times a_{V,c}/a_{r,c}, with a_{V,c} the elements of column c * corresponding to the non-parametric variables. * If variable v appears in a column c_v, then a_{v,c} = 1 iff c = c_v, * with all other entries in this virtual row equal to zero. * If variable v appears in a row, then a_{v,c} is the element in column c * of that row. * * Let v be the first variable with a_{v,c1}/a_{r,c1} != a_{v,c2}/a_{r,c2}. * Then if a_{v,c1}/a_{r,c1} < a_{v,c2}/a_{r,c2}, i.e., * a_{v,c2} a_{r,c1} - a_{v,c1} a_{r,c2} > 0, c1 results in the minimal * increment. Otherwise, it's c2. */ static int lexmin_col_pair(struct isl_tab *tab, int row, int col1, int col2, isl_int tmp) { … } /* Does the index into the tab->var or tab->con array "index" * correspond to a variable in the context tableau? * In particular, it needs to be an index into the tab->var array and * it needs to refer to either one of the first tab->n_param variables or * one of the last tab->n_div variables. */ static int is_parameter_var(struct isl_tab *tab, int index) { … } /* Does column "col" of "tab" refer to a variable in the context tableau? */ static int col_is_parameter_var(struct isl_tab *tab, int col) { … } /* Does row "row" of "tab" refer to a variable in the context tableau? */ static int row_is_parameter_var(struct isl_tab *tab, int row) { … } /* Given a row in the tableau, find and return the column that would * result in the lexicographically smallest, but positive, increment * in the sample point. * If there is no such column, then return tab->n_col. * If anything goes wrong, return -1. */ static int lexmin_pivot_col(struct isl_tab *tab, int row) { … } /* Return the first known violated constraint, i.e., a non-negative * constraint that currently has an either obviously negative value * or a previously determined to be negative value. * * If any constraint has a negative coefficient for the big parameter, * if any, then we return one of these first. */ static int first_neg(struct isl_tab *tab) { … } /* Check whether the invariant that all columns are lexico-positive * is satisfied. This function is not called from the current code * but is useful during debugging. */ static void check_lexpos(struct isl_tab *tab) __attribute__ ((unused)); static void check_lexpos(struct isl_tab *tab) { … } /* Report to the caller that the given constraint is part of an encountered * conflict. */ static int report_conflicting_constraint(struct isl_tab *tab, int con) { … } /* Given a conflicting row in the tableau, report all constraints * involved in the row to the caller. That is, the row itself * (if it represents a constraint) and all constraint columns with * non-zero (and therefore negative) coefficients. */ static int report_conflict(struct isl_tab *tab, int row) { … } /* Resolve all known or obviously violated constraints through pivoting. * In particular, as long as we can find any violated constraint, we * look for a pivoting column that would result in the lexicographically * smallest increment in the sample point. If there is no such column * then the tableau is infeasible. */ static int restore_lexmin(struct isl_tab *tab) WARN_UNUSED; static int restore_lexmin(struct isl_tab *tab) { … } /* Given a row that represents an equality, look for an appropriate * pivoting column. * In particular, if there are any non-zero coefficients among * the non-parameter variables, then we take the last of these * variables. Eliminating this variable in terms of the other * variables and/or parameters does not influence the property * that all column in the initial tableau are lexicographically * positive. The row corresponding to the eliminated variable * will only have non-zero entries below the diagonal of the * initial tableau. That is, we transform * * I I * 1 into a * I I * * If there is no such non-parameter variable, then we are dealing with * pure parameter equality and we pick any parameter with coefficient 1 or -1 * for elimination. This will ensure that the eliminated parameter * always has an integer value whenever all the other parameters are integral. * If there is no such parameter then we return -1. */ static int last_var_col_or_int_par_col(struct isl_tab *tab, int row) { … } /* Add an equality that is known to be valid to the tableau. * We first check if we can eliminate a variable or a parameter. * If not, we add the equality as two inequalities. * In this case, the equality was a pure parameter equality and there * is no need to resolve any constraint violations. * * This function assumes that at least two more rows and at least * two more elements in the constraint array are available in the tableau. */ static struct isl_tab *add_lexmin_valid_eq(struct isl_tab *tab, isl_int *eq) { … } /* Check if the given row is a pure constant. */ static int is_constant(struct isl_tab *tab, int row) { … } /* Is the given row a parametric constant? * That is, does it only involve variables that also appear in the context? */ static int is_parametric_constant(struct isl_tab *tab, int row) { … } /* Add an equality that may or may not be valid to the tableau. * If the resulting row is a pure constant, then it must be zero. * Otherwise, the resulting tableau is empty. * * If the row is not a pure constant, then we add two inequalities, * each time checking that they can be satisfied. * In the end we try to use one of the two constraints to eliminate * a column. * * This function assumes that at least two more rows and at least * two more elements in the constraint array are available in the tableau. */ static int add_lexmin_eq(struct isl_tab *tab, isl_int *eq) WARN_UNUSED; static int add_lexmin_eq(struct isl_tab *tab, isl_int *eq) { … } /* Add an inequality to the tableau, resolving violations using * restore_lexmin. * * This function assumes that at least one more row and at least * one more element in the constraint array are available in the tableau. */ static struct isl_tab *add_lexmin_ineq(struct isl_tab *tab, isl_int *ineq) { … } /* Check if the coefficients of the parameters are all integral. */ static int integer_parameter(struct isl_tab *tab, int row) { … } /* Check if the coefficients of the non-parameter variables are all integral. */ static int integer_variable(struct isl_tab *tab, int row) { … } /* Check if the constant term is integral. */ static int integer_constant(struct isl_tab *tab, int row) { … } #define I_CST … #define I_PAR … #define I_VAR … /* Check for next (non-parameter) variable after "var" (first if var == -1) * that is non-integer and therefore requires a cut and return * the index of the variable. * For parametric tableaus, there are three parts in a row, * the constant, the coefficients of the parameters and the rest. * For each part, we check whether the coefficients in that part * are all integral and if so, set the corresponding flag in *f. * If the constant and the parameter part are integral, then the * current sample value is integral and no cut is required * (irrespective of whether the variable part is integral). */ static int next_non_integer_var(struct isl_tab *tab, int var, int *f) { … } /* Check for first (non-parameter) variable that is non-integer and * therefore requires a cut and return the corresponding row. * For parametric tableaus, there are three parts in a row, * the constant, the coefficients of the parameters and the rest. * For each part, we check whether the coefficients in that part * are all integral and if so, set the corresponding flag in *f. * If the constant and the parameter part are integral, then the * current sample value is integral and no cut is required * (irrespective of whether the variable part is integral). */ static int first_non_integer_row(struct isl_tab *tab, int *f) { … } /* Add a (non-parametric) cut to cut away the non-integral sample * value of the given row. * * If the row is given by * * m r = f + \sum_i a_i y_i * * then the cut is * * c = - {-f/m} + \sum_i {a_i/m} y_i >= 0 * * The big parameter, if any, is ignored, since it is assumed to be big * enough to be divisible by any integer. * If the tableau is actually a parametric tableau, then this function * is only called when all coefficients of the parameters are integral. * The cut therefore has zero coefficients for the parameters. * * The current value is known to be negative, so row_sign, if it * exists, is set accordingly. * * Return the row of the cut or -1. */ static int add_cut(struct isl_tab *tab, int row) { … } #define CUT_ALL … #define CUT_ONE … /* Given a non-parametric tableau, add cuts until an integer * sample point is obtained or until the tableau is determined * to be integer infeasible. * As long as there is any non-integer value in the sample point, * we add appropriate cuts, if possible, for each of these * non-integer values and then resolve the violated * cut constraints using restore_lexmin. * If one of the corresponding rows is equal to an integral * combination of variables/constraints plus a non-integral constant, * then there is no way to obtain an integer point and we return * a tableau that is marked empty. * The parameter cutting_strategy controls the strategy used when adding cuts * to remove non-integer points. CUT_ALL adds all possible cuts * before continuing the search. CUT_ONE adds only one cut at a time. */ static struct isl_tab *cut_to_integer_lexmin(struct isl_tab *tab, int cutting_strategy) { … } /* Check whether all the currently active samples also satisfy the inequality * "ineq" (treated as an equality if eq is set). * Remove those samples that do not. */ static struct isl_tab *check_samples(struct isl_tab *tab, isl_int *ineq, int eq) { … } /* Check whether the sample value of the tableau is finite, * i.e., either the tableau does not use a big parameter, or * all values of the variables are equal to the big parameter plus * some constant. This constant is the actual sample value. */ static int sample_is_finite(struct isl_tab *tab) { … } /* Check if the context tableau of sol has any integer points. * Leave tab in empty state if no integer point can be found. * If an integer point can be found and if moreover it is finite, * then it is added to the list of sample values. * * This function is only called when none of the currently active sample * values satisfies the most recently added constraint. */ static struct isl_tab *check_integer_feasible(struct isl_tab *tab) { … } /* Check if any of the currently active sample values satisfies * the inequality "ineq" (an equality if eq is set). */ static int tab_has_valid_sample(struct isl_tab *tab, isl_int *ineq, int eq) { … } /* Insert a div specified by "div" to the tableau "tab" at position "pos" and * return isl_bool_true if the div is obviously non-negative. */ static isl_bool context_tab_insert_div(struct isl_tab *tab, int pos, __isl_keep isl_vec *div, isl_stat (*add_ineq)(void *user, isl_int *), void *user) { … } /* Add a div specified by "div" to both the main tableau and * the context tableau. In case of the main tableau, we only * need to add an extra div. In the context tableau, we also * need to express the meaning of the div. * Return the index of the div or -1 if anything went wrong. * * The new integer division is added before any unknown integer * divisions in the context to ensure that it does not get * equated to some linear combination involving unknown integer * divisions. */ static int add_div(struct isl_tab *tab, struct isl_context *context, __isl_keep isl_vec *div) { … } /* Return the position of the integer division that is equal to div/denom * if there is one. Otherwise, return a position beyond the integer divisions. */ static int find_div(struct isl_tab *tab, isl_int *div, isl_int denom) { … } /* Return the index of a div that corresponds to "div". * We first check if we already have such a div and if not, we create one. */ static int get_div(struct isl_tab *tab, struct isl_context *context, struct isl_vec *div) { … } /* Add a parametric cut to cut away the non-integral sample value * of the given row. * Let a_i be the coefficients of the constant term and the parameters * and let b_i be the coefficients of the variables or constraints * in basis of the tableau. * Let q be the div q = floor(\sum_i {-a_i} y_i). * * The cut is expressed as * * c = \sum_i -{-a_i} y_i + \sum_i {b_i} x_i + q >= 0 * * If q did not already exist in the context tableau, then it is added first. * If q is in a column of the main tableau then the "+ q" can be accomplished * by setting the corresponding entry to the denominator of the constraint. * If q happens to be in a row of the main tableau, then the corresponding * row needs to be added instead (taking care of the denominators). * Note that this is very unlikely, but perhaps not entirely impossible. * * The current value of the cut is known to be negative (or at least * non-positive), so row_sign is set accordingly. * * Return the row of the cut or -1. */ static int add_parametric_cut(struct isl_tab *tab, int row, struct isl_context *context) { … } /* Construct a tableau for bmap that can be used for computing * the lexicographic minimum (or maximum) of bmap. * If not NULL, then dom is the domain where the minimum * should be computed. In this case, we set up a parametric * tableau with row signs (initialized to "unknown"). * If M is set, then the tableau will use a big parameter. * If max is set, then a maximum should be computed instead of a minimum. * This means that for each variable x, the tableau will contain the variable * x' = M - x, rather than x' = M + x. This in turn means that the coefficient * of the variables in all constraints are negated prior to adding them * to the tableau. */ static __isl_give struct isl_tab *tab_for_lexmin(__isl_keep isl_basic_map *bmap, __isl_keep isl_basic_set *dom, unsigned M, int max) { … } /* Given a main tableau where more than one row requires a split, * determine and return the "best" row to split on. * * If any of the rows requiring a split only involves * variables that also appear in the context tableau, * then the negative part is guaranteed not to have a solution. * It is therefore best to split on any of these rows first. * * Otherwise, * given two rows in the main tableau, if the inequality corresponding * to the first row is redundant with respect to that of the second row * in the current tableau, then it is better to split on the second row, * since in the positive part, both rows will be positive. * (In the negative part a pivot will have to be performed and just about * anything can happen to the sign of the other row.) * * As a simple heuristic, we therefore select the row that makes the most * of the other rows redundant. * * Perhaps it would also be useful to look at the number of constraints * that conflict with any given constraint. * * best is the best row so far (-1 when we have not found any row yet). * best_r is the number of other rows made redundant by row best. * When best is still -1, bset_r is meaningless, but it is initialized * to some arbitrary value (0) anyway. Without this redundant initialization * valgrind may warn about uninitialized memory accesses when isl * is compiled with some versions of gcc. */ static int best_split(struct isl_tab *tab, struct isl_tab *context_tab) { … } static struct isl_basic_set *context_lex_peek_basic_set( struct isl_context *context) { … } static struct isl_tab *context_lex_peek_tab(struct isl_context *context) { … } static void context_lex_add_eq(struct isl_context *context, isl_int *eq, int check, int update) { … } static void context_lex_add_ineq(struct isl_context *context, isl_int *ineq, int check, int update) { … } static isl_stat context_lex_add_ineq_wrap(void *user, isl_int *ineq) { … } /* Check which signs can be obtained by "ineq" on all the currently * active sample values. See row_sign for more information. */ static enum isl_tab_row_sign tab_ineq_sign(struct isl_tab *tab, isl_int *ineq, int strict) { … } static enum isl_tab_row_sign context_lex_ineq_sign(struct isl_context *context, isl_int *ineq, int strict) { … } /* Check whether "ineq" can be added to the tableau without rendering * it infeasible. */ static int context_lex_test_ineq(struct isl_context *context, isl_int *ineq) { … } static int context_lex_get_div(struct isl_context *context, struct isl_tab *tab, struct isl_vec *div) { … } /* Insert a div specified by "div" to the context tableau at position "pos" and * return isl_bool_true if the div is obviously non-negative. * context_tab_add_div will always return isl_bool_true, because all variables * in a isl_context_lex tableau are non-negative. * However, if we are using a big parameter in the context, then this only * reflects the non-negativity of the variable used to _encode_ the * div, i.e., div' = M + div, so we can't draw any conclusions. */ static isl_bool context_lex_insert_div(struct isl_context *context, int pos, __isl_keep isl_vec *div) { … } static int context_lex_detect_equalities(struct isl_context *context, struct isl_tab *tab) { … } static int context_lex_best_split(struct isl_context *context, struct isl_tab *tab) { … } static int context_lex_is_empty(struct isl_context *context) { … } static void *context_lex_save(struct isl_context *context) { … } static void context_lex_restore(struct isl_context *context, void *save) { … } static void context_lex_discard(void *save) { … } static int context_lex_is_ok(struct isl_context *context) { … } /* For each variable in the context tableau, check if the variable can * only attain non-negative values. If so, mark the parameter as non-negative * in the main tableau. This allows for a more direct identification of some * cases of violated constraints. */ static struct isl_tab *tab_detect_nonnegative_parameters(struct isl_tab *tab, struct isl_tab *context_tab) { … } static struct isl_tab *context_lex_detect_nonnegative_parameters( struct isl_context *context, struct isl_tab *tab) { … } static void context_lex_invalidate(struct isl_context *context) { … } static __isl_null struct isl_context *context_lex_free( struct isl_context *context) { … } struct isl_context_op isl_context_lex_op = …; static struct isl_tab *context_tab_for_lexmin(__isl_take isl_basic_set *bset) { … } static struct isl_context *isl_context_lex_alloc(struct isl_basic_set *dom) { … } /* Representation of the context when using generalized basis reduction. * * "shifted" contains the offsets of the unit hypercubes that lie inside the * context. Any rational point in "shifted" can therefore be rounded * up to an integer point in the context. * If the context is constrained by any equality, then "shifted" is not used * as it would be empty. */ struct isl_context_gbr { … }; static struct isl_tab *context_gbr_detect_nonnegative_parameters( struct isl_context *context, struct isl_tab *tab) { … } static struct isl_basic_set *context_gbr_peek_basic_set( struct isl_context *context) { … } static struct isl_tab *context_gbr_peek_tab(struct isl_context *context) { … } /* Initialize the "shifted" tableau of the context, which * contains the constraints of the original tableau shifted * by the sum of all negative coefficients. This ensures * that any rational point in the shifted tableau can * be rounded up to yield an integer point in the original tableau. */ static void gbr_init_shifted(struct isl_context_gbr *cgbr) { … } /* Check if the shifted tableau is non-empty, and if so * use the sample point to construct an integer point * of the context tableau. */ static struct isl_vec *gbr_get_shifted_sample(struct isl_context_gbr *cgbr) { … } static __isl_give isl_basic_set *drop_constant_terms( __isl_take isl_basic_set *bset) { … } static int use_shifted(struct isl_context_gbr *cgbr) { … } static struct isl_vec *gbr_get_sample(struct isl_context_gbr *cgbr) { … } static void check_gbr_integer_feasible(struct isl_context_gbr *cgbr) { … } static struct isl_tab *add_gbr_eq(struct isl_tab *tab, isl_int *eq) { … } /* Add the equality described by "eq" to the context. * If "check" is set, then we check if the context is empty after * adding the equality. * If "update" is set, then we check if the samples are still valid. * * We do not explicitly add shifted copies of the equality to * cgbr->shifted since they would conflict with each other. * Instead, we directly mark cgbr->shifted empty. */ static void context_gbr_add_eq(struct isl_context *context, isl_int *eq, int check, int update) { … } static void add_gbr_ineq(struct isl_context_gbr *cgbr, isl_int *ineq) { … } static void context_gbr_add_ineq(struct isl_context *context, isl_int *ineq, int check, int update) { … } static isl_stat context_gbr_add_ineq_wrap(void *user, isl_int *ineq) { … } static enum isl_tab_row_sign context_gbr_ineq_sign(struct isl_context *context, isl_int *ineq, int strict) { … } /* Check whether "ineq" can be added to the tableau without rendering * it infeasible. */ static int context_gbr_test_ineq(struct isl_context *context, isl_int *ineq) { … } /* Return the column of the last of the variables associated to * a column that has a non-zero coefficient. * This function is called in a context where only coefficients * of parameters or divs can be non-zero. */ static int last_non_zero_var_col(struct isl_tab *tab, isl_int *p) { … } /* Look through all the recently added equalities in the context * to see if we can propagate any of them to the main tableau. * * The newly added equalities in the context are encoded as pairs * of inequalities starting at inequality "first". * * We tentatively add each of these equalities to the main tableau * and if this happens to result in a row with a final coefficient * that is one or negative one, we use it to kill a column * in the main tableau. Otherwise, we discard the tentatively * added row. * This tentative addition of equality constraints turns * on the undo facility of the tableau. Turn it off again * at the end, assuming it was turned off to begin with. * * Return 0 on success and -1 on failure. */ static int propagate_equalities(struct isl_context_gbr *cgbr, struct isl_tab *tab, unsigned first) { … } static int context_gbr_detect_equalities(struct isl_context *context, struct isl_tab *tab) { … } static int context_gbr_get_div(struct isl_context *context, struct isl_tab *tab, struct isl_vec *div) { … } static isl_bool context_gbr_insert_div(struct isl_context *context, int pos, __isl_keep isl_vec *div) { … } static int context_gbr_best_split(struct isl_context *context, struct isl_tab *tab) { … } static int context_gbr_is_empty(struct isl_context *context) { … } struct isl_gbr_tab_undo { … }; static void *context_gbr_save(struct isl_context *context) { … } static void context_gbr_restore(struct isl_context *context, void *save) { … } static void context_gbr_discard(void *save) { … } static int context_gbr_is_ok(struct isl_context *context) { … } static void context_gbr_invalidate(struct isl_context *context) { … } static __isl_null struct isl_context *context_gbr_free( struct isl_context *context) { … } struct isl_context_op isl_context_gbr_op = …; static struct isl_context *isl_context_gbr_alloc(__isl_keep isl_basic_set *dom) { … } /* Allocate a context corresponding to "dom". * The representation specific fields are initialized by * isl_context_lex_alloc or isl_context_gbr_alloc. * The shared "n_unknown" field is initialized to the number * of final unknown integer divisions in "dom". */ static struct isl_context *isl_context_alloc(__isl_keep isl_basic_set *dom) { … } /* Initialize some common fields of "sol", which keeps track * of the solution of an optimization problem on "bmap" over * the domain "dom". * If "max" is set, then a maximization problem is being solved, rather than * a minimization problem, which means that the variables in the * tableau have value "M - x" rather than "M + x". */ static isl_stat sol_init(struct isl_sol *sol, __isl_keep isl_basic_map *bmap, __isl_keep isl_basic_set *dom, int max) { … } /* Construct an isl_sol_map structure for accumulating the solution. * If track_empty is set, then we also keep track of the parts * of the context where there is no solution. * If max is set, then we are solving a maximization, rather than * a minimization problem, which means that the variables in the * tableau have value "M - x" rather than "M + x". */ static struct isl_sol *sol_map_init(__isl_keep isl_basic_map *bmap, __isl_take isl_basic_set *dom, int track_empty, int max) { … } /* Check whether all coefficients of (non-parameter) variables * are non-positive, meaning that no pivots can be performed on the row. */ static int is_critical(struct isl_tab *tab, int row) { … } /* Check whether the inequality represented by vec is strict over the integers, * i.e., there are no integer values satisfying the constraint with * equality. This happens if the gcd of the coefficients is not a divisor * of the constant term. If so, scale the constraint down by the gcd * of the coefficients. */ static int is_strict(struct isl_vec *vec) { … } /* Determine the sign of the given row of the main tableau. * The result is one of * isl_tab_row_pos: always non-negative; no pivot needed * isl_tab_row_neg: always non-positive; pivot * isl_tab_row_any: can be both positive and negative; split * * We first handle some simple cases * - the row sign may be known already * - the row may be obviously non-negative * - the parametric constant may be equal to that of another row * for which we know the sign. This sign will be either "pos" or * "any". If it had been "neg" then we would have pivoted before. * * If none of these cases hold, we check the value of the row for each * of the currently active samples. Based on the signs of these values * we make an initial determination of the sign of the row. * * all zero -> unk(nown) * all non-negative -> pos * all non-positive -> neg * both negative and positive -> all * * If we end up with "all", we are done. * Otherwise, we perform a check for positive and/or negative * values as follows. * * samples neg unk pos * <0 ? Y N Y N * pos any pos * >0 ? Y N Y N * any neg any neg * * There is no special sign for "zero", because we can usually treat zero * as either non-negative or non-positive, whatever works out best. * However, if the row is "critical", meaning that pivoting is impossible * then we don't want to limp zero with the non-positive case, because * then we we would lose the solution for those values of the parameters * where the value of the row is zero. Instead, we treat 0 as non-negative * ensuring a split if the row can attain both zero and negative values. * The same happens when the original constraint was one that could not * be satisfied with equality by any integer values of the parameters. * In this case, we normalize the constraint, but then a value of zero * for the normalized constraint is actually a positive value for the * original constraint, so again we need to treat zero as non-negative. * In both these cases, we have the following decision tree instead: * * all non-negative -> pos * all negative -> neg * both negative and non-negative -> all * * samples neg pos * <0 ? Y N * any pos * >=0 ? Y N * any neg */ static enum isl_tab_row_sign row_sign(struct isl_tab *tab, struct isl_sol *sol, int row) { … } static void find_solutions(struct isl_sol *sol, struct isl_tab *tab); /* Find solutions for values of the parameters that satisfy the given * inequality. * * We currently take a snapshot of the context tableau that is reset * when we return from this function, while we make a copy of the main * tableau, leaving the original main tableau untouched. * These are fairly arbitrary choices. Making a copy also of the context * tableau would obviate the need to undo any changes made to it later, * while taking a snapshot of the main tableau could reduce memory usage. * If we were to switch to taking a snapshot of the main tableau, * we would have to keep in mind that we need to save the row signs * and that we need to do this before saving the current basis * such that the basis has been restore before we restore the row signs. */ static void find_in_pos(struct isl_sol *sol, struct isl_tab *tab, isl_int *ineq) { … } /* Record the absence of solutions for those values of the parameters * that do not satisfy the given inequality with equality. */ static void no_sol_in_strict(struct isl_sol *sol, struct isl_tab *tab, struct isl_vec *ineq) { … } /* Reset all row variables that are marked to have a sign that may * be both positive and negative to have an unknown sign. */ static void reset_any_to_unknown(struct isl_tab *tab) { … } /* Compute the lexicographic minimum of the set represented by the main * tableau "tab" within the context "sol->context_tab". * On entry the sample value of the main tableau is lexicographically * less than or equal to this lexicographic minimum. * Pivots are performed until a feasible point is found, which is then * necessarily equal to the minimum, or until the tableau is found to * be infeasible. Some pivots may need to be performed for only some * feasible values of the context tableau. If so, the context tableau * is split into a part where the pivot is needed and a part where it is not. * * Whenever we enter the main loop, the main tableau is such that no * "obvious" pivots need to be performed on it, where "obvious" means * that the given row can be seen to be negative without looking at * the context tableau. In particular, for non-parametric problems, * no pivots need to be performed on the main tableau. * The caller of find_solutions is responsible for making this property * hold prior to the first iteration of the loop, while restore_lexmin * is called before every other iteration. * * Inside the main loop, we first examine the signs of the rows of * the main tableau within the context of the context tableau. * If we find a row that is always non-positive for all values of * the parameters satisfying the context tableau and negative for at * least one value of the parameters, we perform the appropriate pivot * and start over. An exception is the case where no pivot can be * performed on the row. In this case, we require that the sign of * the row is negative for all values of the parameters (rather than just * non-positive). This special case is handled inside row_sign, which * will say that the row can have any sign if it determines that it can * attain both negative and zero values. * * If we can't find a row that always requires a pivot, but we can find * one or more rows that require a pivot for some values of the parameters * (i.e., the row can attain both positive and negative signs), then we split * the context tableau into two parts, one where we force the sign to be * non-negative and one where we force is to be negative. * The non-negative part is handled by a recursive call (through find_in_pos). * Upon returning from this call, we continue with the negative part and * perform the required pivot. * * If no such rows can be found, all rows are non-negative and we have * found a (rational) feasible point. If we only wanted a rational point * then we are done. * Otherwise, we check if all values of the sample point of the tableau * are integral for the variables. If so, we have found the minimal * integral point and we are done. * If the sample point is not integral, then we need to make a distinction * based on whether the constant term is non-integral or the coefficients * of the parameters. Furthermore, in order to decide how to handle * the non-integrality, we also need to know whether the coefficients * of the other columns in the tableau are integral. This leads * to the following table. The first two rows do not correspond * to a non-integral sample point and are only mentioned for completeness. * * constant parameters other * * int int int | * int int rat | -> no problem * * rat int int -> fail * * rat int rat -> cut * * int rat rat | * rat rat rat | -> parametric cut * * int rat int | * rat rat int | -> split context * * If the parametric constant is completely integral, then there is nothing * to be done. If the constant term is non-integral, but all the other * coefficient are integral, then there is nothing that can be done * and the tableau has no integral solution. * If, on the other hand, one or more of the other columns have rational * coefficients, but the parameter coefficients are all integral, then * we can perform a regular (non-parametric) cut. * Finally, if there is any parameter coefficient that is non-integral, * then we need to involve the context tableau. There are two cases here. * If at least one other column has a rational coefficient, then we * can perform a parametric cut in the main tableau by adding a new * integer division in the context tableau. * If all other columns have integral coefficients, then we need to * enforce that the rational combination of parameters (c + \sum a_i y_i)/m * is always integral. We do this by introducing an integer division * q = floor((c + \sum a_i y_i)/m) and stipulating that its argument should * always be integral in the context tableau, i.e., m q = c + \sum a_i y_i. * Since q is expressed in the tableau as * c + \sum a_i y_i - m q >= 0 * -c - \sum a_i y_i + m q + m - 1 >= 0 * it is sufficient to add the inequality * -c - \sum a_i y_i + m q >= 0 * In the part of the context where this inequality does not hold, the * main tableau is marked as being empty. */ static void find_solutions(struct isl_sol *sol, struct isl_tab *tab) { … } /* Does "sol" contain a pair of partial solutions that could potentially * be merged? * * We currently only check that "sol" is not in an error state * and that there are at least two partial solutions of which the final two * are defined at the same level. */ static int sol_has_mergeable_solutions(struct isl_sol *sol) { … } /* Compute the lexicographic minimum of the set represented by the main * tableau "tab" within the context "sol->context_tab". * * As a preprocessing step, we first transfer all the purely parametric * equalities from the main tableau to the context tableau, i.e., * parameters that have been pivoted to a row. * These equalities are ignored by the main algorithm, because the * corresponding rows may not be marked as being non-negative. * In parts of the context where the added equality does not hold, * the main tableau is marked as being empty. * * Before we embark on the actual computation, we save a copy * of the context. When we return, we check if there are any * partial solutions that can potentially be merged. If so, * we perform a rollback to the initial state of the context. * The merging of partial solutions happens inside calls to * sol_dec_level that are pushed onto the undo stack of the context. * If there are no partial solutions that can potentially be merged * then the rollback is skipped as it would just be wasted effort. */ static void find_solutions_main(struct isl_sol *sol, struct isl_tab *tab) { … } /* Check if integer division "div" of "dom" also occurs in "bmap". * If so, return its position within the divs. * Otherwise, return a position beyond the integer divisions. */ static int find_context_div(__isl_keep isl_basic_map *bmap, __isl_keep isl_basic_set *dom, unsigned div) { … } /* The correspondence between the variables in the main tableau, * the context tableau, and the input map and domain is as follows. * The first n_param and the last n_div variables of the main tableau * form the variables of the context tableau. * In the basic map, these n_param variables correspond to the * parameters and the input dimensions. In the domain, they correspond * to the parameters and the set dimensions. * The n_div variables correspond to the integer divisions in the domain. * To ensure that everything lines up, we may need to copy some of the * integer divisions of the domain to the map. These have to be placed * in the same order as those in the context and they have to be placed * after any other integer divisions that the map may have. * This function performs the required reordering. */ static __isl_give isl_basic_map *align_context_divs( __isl_take isl_basic_map *bmap, __isl_keep isl_basic_set *dom) { … } /* Base case of isl_tab_basic_map_partial_lexopt, after removing * some obvious symmetries. * * We make sure the divs in the domain are properly ordered, * because they will be added one by one in the given order * during the construction of the solution map. * Furthermore, make sure that the known integer divisions * appear before any unknown integer division because the solution * may depend on the known integer divisions, while anything that * depends on any variable starting from the first unknown integer * division is ignored in sol_pma_add. */ static struct isl_sol *basic_map_partial_lexopt_base_sol( __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom, __isl_give isl_set **empty, int max, struct isl_sol *(*init)(__isl_keep isl_basic_map *bmap, __isl_take isl_basic_set *dom, int track_empty, int max)) { … } /* Base case of isl_tab_basic_map_partial_lexopt, after removing * some obvious symmetries. * * We call basic_map_partial_lexopt_base_sol and extract the results. */ static __isl_give isl_map *basic_map_partial_lexopt_base( __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom, __isl_give isl_set **empty, int max) { … } /* Return a count of the number of occurrences of the "n" first * variables in the inequality constraints of "bmap". */ static __isl_give int *count_occurrences(__isl_keep isl_basic_map *bmap, int n) { … } /* Do all of the "n" variables with non-zero coefficients in "c" * occur in exactly a single constraint. * "occurrences" is an array of length "n" containing the number * of occurrences of each of the variables in the inequality constraints. */ static int single_occurrence(int n, isl_int *c, int *occurrences) { … } /* Do all of the "n" initial variables that occur in inequality constraint * "ineq" of "bmap" only occur in that constraint? */ static int all_single_occurrence(__isl_keep isl_basic_map *bmap, int ineq, int n) { … } /* Structure used during detection of parallel constraints. * n_in: number of "input" variables: isl_dim_param + isl_dim_in * n_out: number of "output" variables: isl_dim_out + isl_dim_div * val: the coefficients of the output variables */ struct isl_constraint_equal_info { … }; /* Check whether the coefficients of the output variables * of the constraint in "entry" are equal to info->val. */ static isl_bool constraint_equal(const void *entry, const void *val) { … } /* Check whether "bmap" has a pair of constraints that have * the same coefficients for the output variables. * Note that the coefficients of the existentially quantified * variables need to be zero since the existentially quantified * of the result are usually not the same as those of the input. * Furthermore, check that each of the input variables that occur * in those constraints does not occur in any other constraint. * If so, return true and return the row indices of the two constraints * in *first and *second. */ static isl_bool parallel_constraints(__isl_keep isl_basic_map *bmap, int *first, int *second) { … } /* Given a set of upper bounds in "var", add constraints to "bset" * that make the i-th bound smallest. * * In particular, if there are n bounds b_i, then add the constraints * * b_i <= b_j for j > i * b_i < b_j for j < i */ static __isl_give isl_basic_set *select_minimum(__isl_take isl_basic_set *bset, __isl_keep isl_mat *var, int i) { … } /* Given a set of upper bounds on the last "input" variable m, * construct a set that assigns the minimal upper bound to m, i.e., * construct a set that divides the space into cells where one * of the upper bounds is smaller than all the others and assign * this upper bound to m. * * In particular, if there are n bounds b_i, then the result * consists of n basic sets, each one of the form * * m = b_i * b_i <= b_j for j > i * b_i < b_j for j < i */ static __isl_give isl_set *set_minimum(__isl_take isl_space *space, __isl_take isl_mat *var) { … } /* Given that the last input variable of "bmap" represents the minimum * of the bounds in "cst", check whether we need to split the domain * based on which bound attains the minimum. * * A split is needed when the minimum appears in an integer division * or in an equality. Otherwise, it is only needed if it appears in * an upper bound that is different from the upper bounds on which it * is defined. */ static isl_bool need_split_basic_map(__isl_keep isl_basic_map *bmap, __isl_keep isl_mat *cst) { … } /* Given that the last set variable of "bset" represents the minimum * of the bounds in "cst", check whether we need to split the domain * based on which bound attains the minimum. * * We simply call need_split_basic_map here. This is safe because * the position of the minimum is computed from "cst" and not * from "bmap". */ static isl_bool need_split_basic_set(__isl_keep isl_basic_set *bset, __isl_keep isl_mat *cst) { … } /* Given that the last set variable of "set" represents the minimum * of the bounds in "cst", check whether we need to split the domain * based on which bound attains the minimum. */ static isl_bool need_split_set(__isl_keep isl_set *set, __isl_keep isl_mat *cst) { … } /* Given a map of which the last input variable is the minimum * of the bounds in "cst", split each basic set in the set * in pieces where one of the bounds is (strictly) smaller than the others. * This subdivision is given in "min_expr". * The variable is subsequently projected out. * * We only do the split when it is needed. * For example if the last input variable m = min(a,b) and the only * constraints in the given basic set are lower bounds on m, * i.e., l <= m = min(a,b), then we can simply project out m * to obtain l <= a and l <= b, without having to split on whether * m is equal to a or b. */ static __isl_give isl_map *split_domain(__isl_take isl_map *opt, __isl_take isl_set *min_expr, __isl_take isl_mat *cst) { … } /* Given a set of which the last set variable is the minimum * of the bounds in "cst", split each basic set in the set * in pieces where one of the bounds is (strictly) smaller than the others. * This subdivision is given in "min_expr". * The variable is subsequently projected out. */ static __isl_give isl_set *split(__isl_take isl_set *empty, __isl_take isl_set *min_expr, __isl_take isl_mat *cst) { … } static __isl_give isl_map *basic_map_partial_lexopt( __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom, __isl_give isl_set **empty, int max); /* This function is called from basic_map_partial_lexopt_symm. * The last variable of "bmap" and "dom" corresponds to the minimum * of the bounds in "cst". "map_space" is the space of the original * input relation (of basic_map_partial_lexopt_symm) and "set_space" * is the space of the original domain. * * We recursively call basic_map_partial_lexopt and then plug in * the definition of the minimum in the result. */ static __isl_give isl_map *basic_map_partial_lexopt_symm_core( __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom, __isl_give isl_set **empty, int max, __isl_take isl_mat *cst, __isl_take isl_space *map_space, __isl_take isl_space *set_space) { … } /* Extract a domain from "bmap" for the purpose of computing * a lexicographic optimum. * * This function is only called when the caller wants to compute a full * lexicographic optimum, i.e., without specifying a domain. In this case, * the caller is not interested in the part of the domain space where * there is no solution and the domain can be initialized to those constraints * of "bmap" that only involve the parameters and the input dimensions. * This relieves the parametric programming engine from detecting those * inequalities and transferring them to the context. More importantly, * it ensures that those inequalities are transferred first and not * intermixed with inequalities that actually split the domain. * * If the caller does not require the absence of existentially quantified * variables in the result (i.e., if ISL_OPT_QE is not set in "flags"), * then the actual domain of "bmap" can be used. This ensures that * the domain does not need to be split at all just to separate out * pieces of the domain that do not have a solution from piece that do. * This domain cannot be used in general because it may involve * (unknown) existentially quantified variables which will then also * appear in the solution. */ static __isl_give isl_basic_set *extract_domain(__isl_keep isl_basic_map *bmap, unsigned flags) { … } #undef TYPE #define TYPE … #undef SUFFIX #define SUFFIX #include "isl_tab_lexopt_templ.c" /* Extract the subsequence of the sample value of "tab" * starting at "pos" and of length "len". */ static __isl_give isl_vec *extract_sample_sequence(struct isl_tab *tab, int pos, int len) { … } /* Check if the sequence of variables starting at "pos" * represents a trivial solution according to "trivial". * That is, is the result of applying "trivial" to this sequence * equal to the zero vector? */ static isl_bool region_is_trivial(struct isl_tab *tab, int pos, __isl_keep isl_mat *trivial) { … } /* Global internal data for isl_tab_basic_set_non_trivial_lexmin. * * "n_op" is the number of initial coordinates to optimize, * as passed to isl_tab_basic_set_non_trivial_lexmin. * "region" is the "n_region"-sized array of regions passed * to isl_tab_basic_set_non_trivial_lexmin. * * "tab" is the tableau that corresponds to the ILP problem. * "local" is an array of local data structure, one for each * (potential) level of the backtracking procedure of * isl_tab_basic_set_non_trivial_lexmin. * "v" is a pre-allocated vector that can be used for adding * constraints to the tableau. * * "sol" contains the best solution found so far. * It is initialized to a vector of size zero. */ struct isl_lexmin_data { … }; /* Return the index of the first trivial region, "n_region" if all regions * are non-trivial or -1 in case of error. */ static int first_trivial_region(struct isl_lexmin_data *data) { … } /* Check if the solution is optimal, i.e., whether the first * n_op entries are zero. */ static int is_optimal(__isl_keep isl_vec *sol, int n_op) { … } /* Add constraints to "tab" that ensure that any solution is significantly * better than that represented by "sol". That is, find the first * relevant (within first n_op) non-zero coefficient and force it (along * with all previous coefficients) to be zero. * If the solution is already optimal (all relevant coefficients are zero), * then just mark the table as empty. * "n_zero" is the number of coefficients that have been forced zero * by previous calls to this function at the same level. * Return the updated number of forced zero coefficients or -1 on error. * * This function assumes that at least 2 * (n_op - n_zero) more rows and * at least 2 * (n_op - n_zero) more elements in the constraint array * are available in the tableau. */ static int force_better_solution(struct isl_tab *tab, __isl_keep isl_vec *sol, int n_op, int n_zero) { … } /* Fix triviality direction "dir" of the given region to zero. * * This function assumes that at least two more rows and at least * two more elements in the constraint array are available in the tableau. */ static isl_stat fix_zero(struct isl_tab *tab, struct isl_trivial_region *region, int dir, struct isl_lexmin_data *data) { … } /* This function selects case "side" for non-triviality region "region", * assuming all the equality constraints have been imposed already. * In particular, the triviality direction side/2 is made positive * if side is even and made negative if side is odd. * * This function assumes that at least one more row and at least * one more element in the constraint array are available in the tableau. */ static struct isl_tab *pos_neg(struct isl_tab *tab, struct isl_trivial_region *region, int side, struct isl_lexmin_data *data) { … } /* Local data at each level of the backtracking procedure of * isl_tab_basic_set_non_trivial_lexmin. * * "update" is set if a solution has been found in the current case * of this level, such that a better solution needs to be enforced * in the next case. * "n_zero" is the number of initial coordinates that have already * been forced to be zero at this level. * "region" is the non-triviality region considered at this level. * "side" is the index of the current case at this level. * "n" is the number of triviality directions. * "snap" is a snapshot of the tableau holding a state that needs * to be satisfied by all subsequent cases. */ struct isl_local_region { … }; /* Initialize the global data structure "data" used while solving * the ILP problem "bset". */ static isl_stat init_lexmin_data(struct isl_lexmin_data *data, __isl_keep isl_basic_set *bset) { … } /* Mark all outer levels as requiring a better solution * in the next cases. */ static void update_outer_levels(struct isl_lexmin_data *data, int level) { … } /* Initialize "local" to refer to region "region" and * to initiate processing at this level. */ static isl_stat init_local_region(struct isl_local_region *local, int region, struct isl_lexmin_data *data) { … } /* What to do next after entering a level of the backtracking procedure. * * error: some error has occurred; abort * done: an optimal solution has been found; stop search * backtrack: backtrack to the previous level * handle: add the constraints for the current level and * move to the next level */ enum isl_next { … }; /* Have all cases of the current region been considered? * If there are n directions, then there are 2n cases. * * The constraints in the current tableau are imposed * in all subsequent cases. This means that if the current * tableau is empty, then none of those cases should be considered * anymore and all cases have effectively been considered. */ static int finished_all_cases(struct isl_local_region *local, struct isl_lexmin_data *data) { … } /* Enter level "level" of the backtracking search and figure out * what to do next. "init" is set if the level was entered * from a higher level and needs to be initialized. * Otherwise, the level is entered as a result of backtracking and * the tableau needs to be restored to a position that can * be used for the next case at this level. * The snapshot is assumed to have been saved in the previous case, * before the constraints specific to that case were added. * * In the initialization case, the local region is initialized * to point to the first violated region. * If the constraints of all regions are satisfied by the current * sample of the tableau, then tell the caller to continue looking * for a better solution or to stop searching if an optimal solution * has been found. * * If the tableau is empty or if all cases at the current level * have been considered, then the caller needs to backtrack as well. */ static enum isl_next enter_level(int level, int init, struct isl_lexmin_data *data) { … } /* If a solution has been found in the previous case at this level * (marked by local->update being set), then add constraints * that enforce a better solution in the present and all following cases. * The constraints only need to be imposed once because they are * included in the snapshot (taken in pick_side) that will be used in * subsequent cases. */ static isl_stat better_next_side(struct isl_local_region *local, struct isl_lexmin_data *data) { … } /* Add constraints to data->tab that select the current case (local->side) * at the current level. * * If the linear combinations v should not be zero, then the cases are * v_0 >= 1 * v_0 <= -1 * v_0 = 0 and v_1 >= 1 * v_0 = 0 and v_1 <= -1 * v_0 = 0 and v_1 = 0 and v_2 >= 1 * v_0 = 0 and v_1 = 0 and v_2 <= -1 * ... * in this order. * * A snapshot is taken after the equality constraint (if any) has been added * such that the next case can start off from this position. * The rollback to this position is performed in enter_level. */ static isl_stat pick_side(struct isl_local_region *local, struct isl_lexmin_data *data) { … } /* Free the memory associated to "data". */ static void clear_lexmin_data(struct isl_lexmin_data *data) { … } /* Return the lexicographically smallest non-trivial solution of the * given ILP problem. * * All variables are assumed to be non-negative. * * n_op is the number of initial coordinates to optimize. * That is, once a solution has been found, we will only continue looking * for solutions that result in significantly better values for those * initial coordinates. That is, we only continue looking for solutions * that increase the number of initial zeros in this sequence. * * A solution is non-trivial, if it is non-trivial on each of the * specified regions. Each region represents a sequence of * triviality directions on a sequence of variables that starts * at a given position. A solution is non-trivial on such a region if * at least one of the triviality directions is non-zero * on that sequence of variables. * * Whenever a conflict is encountered, all constraints involved are * reported to the caller through a call to "conflict". * * We perform a simple branch-and-bound backtracking search. * Each level in the search represents an initially trivial region * that is forced to be non-trivial. * At each level we consider 2 * n cases, where n * is the number of triviality directions. * In terms of those n directions v_i, we consider the cases * v_0 >= 1 * v_0 <= -1 * v_0 = 0 and v_1 >= 1 * v_0 = 0 and v_1 <= -1 * v_0 = 0 and v_1 = 0 and v_2 >= 1 * v_0 = 0 and v_1 = 0 and v_2 <= -1 * ... * in this order. */ __isl_give isl_vec *isl_tab_basic_set_non_trivial_lexmin( __isl_take isl_basic_set *bset, int n_op, int n_region, struct isl_trivial_region *region, int (*conflict)(int con, void *user), void *user) { … } /* Wrapper for a tableau that is used for computing * the lexicographically smallest rational point of a non-negative set. * This point is represented by the sample value of "tab", * unless "tab" is empty. */ struct isl_tab_lexmin { … }; /* Free "tl" and return NULL. */ __isl_null isl_tab_lexmin *isl_tab_lexmin_free(__isl_take isl_tab_lexmin *tl) { … } /* Construct an isl_tab_lexmin for computing * the lexicographically smallest rational point in "bset", * assuming that all variables are non-negative. */ __isl_give isl_tab_lexmin *isl_tab_lexmin_from_basic_set( __isl_take isl_basic_set *bset) { … } /* Return the dimension of the set represented by "tl". */ int isl_tab_lexmin_dim(__isl_keep isl_tab_lexmin *tl) { … } /* Add the equality with coefficients "eq" to "tl", updating the optimal * solution if needed. * The equality is added as two opposite inequality constraints. */ __isl_give isl_tab_lexmin *isl_tab_lexmin_add_eq(__isl_take isl_tab_lexmin *tl, isl_int *eq) { … } /* Add cuts to "tl" until the sample value reaches an integer value or * until the result becomes empty. */ __isl_give isl_tab_lexmin *isl_tab_lexmin_cut_to_integer( __isl_take isl_tab_lexmin *tl) { … } /* Return the lexicographically smallest rational point in the basic set * from which "tl" was constructed. * If the original input was empty, then return a zero-length vector. */ __isl_give isl_vec *isl_tab_lexmin_get_solution(__isl_keep isl_tab_lexmin *tl) { … } struct isl_sol_pma { … }; static void sol_pma_free(struct isl_sol *sol) { … } /* This function is called for parts of the context where there is * no solution, with "bset" corresponding to the context tableau. * Simply add the basic set to the set "empty". */ static void sol_pma_add_empty(struct isl_sol_pma *sol, __isl_take isl_basic_set *bset) { … } /* Given a basic set "dom" that represents the context and a tuple of * affine expressions "maff" defined over this domain, construct * an isl_pw_multi_aff with a single cell corresponding to "dom" and * the affine expressions in "maff". */ static void sol_pma_add(struct isl_sol_pma *sol, __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *maff) { … } static void sol_pma_add_empty_wrap(struct isl_sol *sol, __isl_take isl_basic_set *bset) { … } static void sol_pma_add_wrap(struct isl_sol *sol, __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma) { … } /* Construct an isl_sol_pma structure for accumulating the solution. * If track_empty is set, then we also keep track of the parts * of the context where there is no solution. * If max is set, then we are solving a maximization, rather than * a minimization problem, which means that the variables in the * tableau have value "M - x" rather than "M + x". */ static struct isl_sol *sol_pma_init(__isl_keep isl_basic_map *bmap, __isl_take isl_basic_set *dom, int track_empty, int max) { … } /* Base case of isl_tab_basic_map_partial_lexopt, after removing * some obvious symmetries. * * We call basic_map_partial_lexopt_base_sol and extract the results. */ static __isl_give isl_pw_multi_aff *basic_map_partial_lexopt_base_pw_multi_aff( __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom, __isl_give isl_set **empty, int max) { … } /* Given that the last input variable of "maff" represents the minimum * of some bounds, check whether we need to plug in the expression * of the minimum. * * In particular, check if the last input variable appears in any * of the expressions in "maff". */ static isl_bool need_substitution(__isl_keep isl_multi_aff *maff) { … } /* Given a set of upper bounds on the last "input" variable m, * construct a piecewise affine expression that selects * the minimal upper bound to m, i.e., * divide the space into cells where one * of the upper bounds is smaller than all the others and select * this upper bound on that cell. * * In particular, if there are n bounds b_i, then the result * consists of n cell, each one of the form * * b_i <= b_j for j > i * b_i < b_j for j < i * * The affine expression on this cell is * * b_i */ static __isl_give isl_pw_aff *set_minimum_pa(__isl_take isl_space *space, __isl_take isl_mat *var) { … } /* Given a piecewise multi-affine expression of which the last input variable * is the minimum of the bounds in "cst", plug in the value of the minimum. * This minimum expression is given in "min_expr_pa". * The set "min_expr" contains the same information, but in the form of a set. * The variable is subsequently projected out. * * The implementation is similar to those of "split" and "split_domain". * If the variable appears in a given expression, then minimum expression * is plugged in. Otherwise, if the variable appears in the constraints * and a split is required, then the domain is split. Otherwise, no split * is performed. */ static __isl_give isl_pw_multi_aff *split_domain_pma( __isl_take isl_pw_multi_aff *opt, __isl_take isl_pw_aff *min_expr_pa, __isl_take isl_set *min_expr, __isl_take isl_mat *cst) { … } static __isl_give isl_pw_multi_aff *basic_map_partial_lexopt_pw_multi_aff( __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom, __isl_give isl_set **empty, int max); /* This function is called from basic_map_partial_lexopt_symm. * The last variable of "bmap" and "dom" corresponds to the minimum * of the bounds in "cst". "map_space" is the space of the original * input relation (of basic_map_partial_lexopt_symm) and "set_space" * is the space of the original domain. * * We recursively call basic_map_partial_lexopt and then plug in * the definition of the minimum in the result. */ static __isl_give isl_pw_multi_aff * basic_map_partial_lexopt_symm_core_pw_multi_aff( __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom, __isl_give isl_set **empty, int max, __isl_take isl_mat *cst, __isl_take isl_space *map_space, __isl_take isl_space *set_space) { … } #undef TYPE #define TYPE … #undef SUFFIX #define SUFFIX … #include "isl_tab_lexopt_templ.c"
Codebase Browser
llvm