isl_lp.c 9.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371
  1. /*
  2. * Copyright 2008-2009 Katholieke Universiteit Leuven
  3. *
  4. * Use of this software is governed by the MIT license
  5. *
  6. * Written by Sven Verdoolaege, K.U.Leuven, Departement
  7. * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
  8. */
  9. #include <isl_ctx_private.h>
  10. #include <isl_map_private.h>
  11. #include <isl/lp.h>
  12. #include <isl_seq.h>
  13. #include "isl_tab.h"
  14. #include <isl_options_private.h>
  15. #include <isl_local_space_private.h>
  16. #include <isl_aff_private.h>
  17. #include <isl_mat_private.h>
  18. #include <isl_val_private.h>
  19. #include <isl_vec_private.h>
  20. #include <bset_to_bmap.c>
  21. #include <set_to_map.c>
  22. enum isl_lp_result isl_tab_solve_lp(__isl_keep isl_basic_map *bmap,
  23. int maximize, isl_int *f, isl_int denom, isl_int *opt,
  24. isl_int *opt_denom, __isl_give isl_vec **sol)
  25. {
  26. struct isl_tab *tab;
  27. enum isl_lp_result res;
  28. isl_size dim = isl_basic_map_dim(bmap, isl_dim_all);
  29. if (dim < 0)
  30. return isl_lp_error;
  31. if (maximize)
  32. isl_seq_neg(f, f, 1 + dim);
  33. bmap = isl_basic_map_gauss(bmap, NULL);
  34. tab = isl_tab_from_basic_map(bmap, 0);
  35. res = isl_tab_min(tab, f, denom, opt, opt_denom, 0);
  36. if (res == isl_lp_ok && sol) {
  37. *sol = isl_tab_get_sample_value(tab);
  38. if (!*sol)
  39. res = isl_lp_error;
  40. }
  41. isl_tab_free(tab);
  42. if (maximize)
  43. isl_seq_neg(f, f, 1 + dim);
  44. if (maximize && opt)
  45. isl_int_neg(*opt, *opt);
  46. return res;
  47. }
  48. /* Given a basic map "bmap" and an affine combination of the variables "f"
  49. * with denominator "denom", set *opt / *opt_denom to the minimal
  50. * (or maximal if "maximize" is true) value attained by f/d over "bmap",
  51. * assuming the basic map is not empty and the expression cannot attain
  52. * arbitrarily small (or large) values.
  53. * If opt_denom is NULL, then *opt is rounded up (or down)
  54. * to the nearest integer.
  55. * The return value reflects the nature of the result (empty, unbounded,
  56. * minimal or maximal value returned in *opt).
  57. */
  58. enum isl_lp_result isl_basic_map_solve_lp(__isl_keep isl_basic_map *bmap,
  59. int max, isl_int *f, isl_int d, isl_int *opt, isl_int *opt_denom,
  60. __isl_give isl_vec **sol)
  61. {
  62. if (sol)
  63. *sol = NULL;
  64. if (!bmap)
  65. return isl_lp_error;
  66. return isl_tab_solve_lp(bmap, max, f, d, opt, opt_denom, sol);
  67. }
  68. enum isl_lp_result isl_basic_set_solve_lp(__isl_keep isl_basic_set *bset,
  69. int max, isl_int *f, isl_int d, isl_int *opt, isl_int *opt_denom,
  70. __isl_give isl_vec **sol)
  71. {
  72. return isl_basic_map_solve_lp(bset_to_bmap(bset), max,
  73. f, d, opt, opt_denom, sol);
  74. }
  75. enum isl_lp_result isl_map_solve_lp(__isl_keep isl_map *map, int max,
  76. isl_int *f, isl_int d, isl_int *opt,
  77. isl_int *opt_denom,
  78. __isl_give isl_vec **sol)
  79. {
  80. int i;
  81. isl_int o;
  82. isl_int t;
  83. isl_int opt_i;
  84. isl_int opt_denom_i;
  85. enum isl_lp_result res;
  86. int max_div;
  87. isl_vec *v = NULL;
  88. if (!map)
  89. return isl_lp_error;
  90. if (map->n == 0)
  91. return isl_lp_empty;
  92. max_div = 0;
  93. for (i = 0; i < map->n; ++i)
  94. if (map->p[i]->n_div > max_div)
  95. max_div = map->p[i]->n_div;
  96. if (max_div > 0) {
  97. isl_size total = isl_map_dim(map, isl_dim_all);
  98. if (total < 0)
  99. return isl_lp_error;
  100. v = isl_vec_alloc(map->ctx, 1 + total + max_div);
  101. if (!v)
  102. return isl_lp_error;
  103. isl_seq_cpy(v->el, f, 1 + total);
  104. isl_seq_clr(v->el + 1 + total, max_div);
  105. f = v->el;
  106. }
  107. if (!opt && map->n > 1 && sol) {
  108. isl_int_init(o);
  109. opt = &o;
  110. }
  111. if (map->n > 0)
  112. isl_int_init(opt_i);
  113. if (map->n > 0 && opt_denom) {
  114. isl_int_init(opt_denom_i);
  115. isl_int_init(t);
  116. }
  117. res = isl_basic_map_solve_lp(map->p[0], max, f, d,
  118. opt, opt_denom, sol);
  119. if (res == isl_lp_error || res == isl_lp_unbounded)
  120. goto done;
  121. if (sol)
  122. *sol = NULL;
  123. for (i = 1; i < map->n; ++i) {
  124. isl_vec *sol_i = NULL;
  125. enum isl_lp_result res_i;
  126. int better;
  127. res_i = isl_basic_map_solve_lp(map->p[i], max, f, d,
  128. &opt_i,
  129. opt_denom ? &opt_denom_i : NULL,
  130. sol ? &sol_i : NULL);
  131. if (res_i == isl_lp_error || res_i == isl_lp_unbounded) {
  132. res = res_i;
  133. goto done;
  134. }
  135. if (res_i == isl_lp_empty)
  136. continue;
  137. if (res == isl_lp_empty) {
  138. better = 1;
  139. } else if (!opt_denom) {
  140. if (max)
  141. better = isl_int_gt(opt_i, *opt);
  142. else
  143. better = isl_int_lt(opt_i, *opt);
  144. } else {
  145. isl_int_mul(t, opt_i, *opt_denom);
  146. isl_int_submul(t, *opt, opt_denom_i);
  147. if (max)
  148. better = isl_int_is_pos(t);
  149. else
  150. better = isl_int_is_neg(t);
  151. }
  152. if (better) {
  153. res = res_i;
  154. if (opt)
  155. isl_int_set(*opt, opt_i);
  156. if (opt_denom)
  157. isl_int_set(*opt_denom, opt_denom_i);
  158. if (sol) {
  159. isl_vec_free(*sol);
  160. *sol = sol_i;
  161. }
  162. } else
  163. isl_vec_free(sol_i);
  164. }
  165. done:
  166. isl_vec_free(v);
  167. if (map->n > 0 && opt_denom) {
  168. isl_int_clear(opt_denom_i);
  169. isl_int_clear(t);
  170. }
  171. if (map->n > 0)
  172. isl_int_clear(opt_i);
  173. if (opt == &o)
  174. isl_int_clear(o);
  175. return res;
  176. }
  177. enum isl_lp_result isl_set_solve_lp(__isl_keep isl_set *set, int max,
  178. isl_int *f, isl_int d, isl_int *opt,
  179. isl_int *opt_denom,
  180. __isl_give isl_vec **sol)
  181. {
  182. return isl_map_solve_lp(set_to_map(set), max,
  183. f, d, opt, opt_denom, sol);
  184. }
  185. /* Return the optimal (rational) value of "obj" over "bset", assuming
  186. * that "obj" and "bset" have aligned parameters and divs.
  187. * If "max" is set, then the maximal value is computed.
  188. * Otherwise, the minimal value is computed.
  189. *
  190. * Return infinity or negative infinity if the optimal value is unbounded and
  191. * NaN if "bset" is empty.
  192. *
  193. * Call isl_basic_set_solve_lp and translate the results.
  194. */
  195. static __isl_give isl_val *basic_set_opt_lp(
  196. __isl_keep isl_basic_set *bset, int max, __isl_keep isl_aff *obj)
  197. {
  198. isl_ctx *ctx;
  199. isl_val *res;
  200. enum isl_lp_result lp_res;
  201. if (!bset || !obj)
  202. return NULL;
  203. ctx = isl_aff_get_ctx(obj);
  204. res = isl_val_alloc(ctx);
  205. if (!res)
  206. return NULL;
  207. lp_res = isl_basic_set_solve_lp(bset, max, obj->v->el + 1,
  208. obj->v->el[0], &res->n, &res->d, NULL);
  209. if (lp_res == isl_lp_ok)
  210. return isl_val_normalize(res);
  211. isl_val_free(res);
  212. if (lp_res == isl_lp_error)
  213. return NULL;
  214. if (lp_res == isl_lp_empty)
  215. return isl_val_nan(ctx);
  216. if (max)
  217. return isl_val_infty(ctx);
  218. else
  219. return isl_val_neginfty(ctx);
  220. }
  221. /* Return the optimal (rational) value of "obj" over "bset", assuming
  222. * that "obj" and "bset" have aligned parameters.
  223. * If "max" is set, then the maximal value is computed.
  224. * Otherwise, the minimal value is computed.
  225. *
  226. * Return infinity or negative infinity if the optimal value is unbounded and
  227. * NaN if "bset" is empty.
  228. *
  229. * Align the divs of "bset" and "obj" and call basic_set_opt_lp.
  230. */
  231. static __isl_give isl_val *isl_basic_set_opt_lp_val_aligned(
  232. __isl_keep isl_basic_set *bset, int max, __isl_keep isl_aff *obj)
  233. {
  234. int *exp1 = NULL;
  235. int *exp2 = NULL;
  236. isl_ctx *ctx;
  237. isl_mat *bset_div = NULL;
  238. isl_mat *div = NULL;
  239. isl_val *res;
  240. isl_size bset_n_div, obj_n_div;
  241. if (!bset || !obj)
  242. return NULL;
  243. ctx = isl_aff_get_ctx(obj);
  244. if (!isl_space_is_equal(bset->dim, obj->ls->dim))
  245. isl_die(ctx, isl_error_invalid,
  246. "spaces don't match", return NULL);
  247. bset_n_div = isl_basic_set_dim(bset, isl_dim_div);
  248. obj_n_div = isl_aff_dim(obj, isl_dim_div);
  249. if (bset_n_div < 0 || obj_n_div < 0)
  250. return NULL;
  251. if (bset_n_div == 0 && obj_n_div == 0)
  252. return basic_set_opt_lp(bset, max, obj);
  253. bset = isl_basic_set_copy(bset);
  254. obj = isl_aff_copy(obj);
  255. bset_div = isl_basic_set_get_divs(bset);
  256. exp1 = isl_alloc_array(ctx, int, bset_n_div);
  257. exp2 = isl_alloc_array(ctx, int, obj_n_div);
  258. if (!bset_div || (bset_n_div && !exp1) || (obj_n_div && !exp2))
  259. goto error;
  260. div = isl_merge_divs(bset_div, obj->ls->div, exp1, exp2);
  261. bset = isl_basic_set_expand_divs(bset, isl_mat_copy(div), exp1);
  262. obj = isl_aff_expand_divs(obj, isl_mat_copy(div), exp2);
  263. res = basic_set_opt_lp(bset, max, obj);
  264. isl_mat_free(bset_div);
  265. isl_mat_free(div);
  266. free(exp1);
  267. free(exp2);
  268. isl_basic_set_free(bset);
  269. isl_aff_free(obj);
  270. return res;
  271. error:
  272. isl_mat_free(div);
  273. isl_mat_free(bset_div);
  274. free(exp1);
  275. free(exp2);
  276. isl_basic_set_free(bset);
  277. isl_aff_free(obj);
  278. return NULL;
  279. }
  280. /* Return the optimal (rational) value of "obj" over "bset".
  281. * If "max" is set, then the maximal value is computed.
  282. * Otherwise, the minimal value is computed.
  283. *
  284. * Return infinity or negative infinity if the optimal value is unbounded and
  285. * NaN if "bset" is empty.
  286. */
  287. static __isl_give isl_val *isl_basic_set_opt_lp_val(
  288. __isl_keep isl_basic_set *bset, int max, __isl_keep isl_aff *obj)
  289. {
  290. isl_bool equal;
  291. isl_val *res;
  292. if (!bset || !obj)
  293. return NULL;
  294. equal = isl_basic_set_space_has_equal_params(bset, obj->ls->dim);
  295. if (equal < 0)
  296. return NULL;
  297. if (equal)
  298. return isl_basic_set_opt_lp_val_aligned(bset, max, obj);
  299. bset = isl_basic_set_copy(bset);
  300. obj = isl_aff_copy(obj);
  301. bset = isl_basic_set_align_params(bset, isl_aff_get_domain_space(obj));
  302. obj = isl_aff_align_params(obj, isl_basic_set_get_space(bset));
  303. res = isl_basic_set_opt_lp_val_aligned(bset, max, obj);
  304. isl_basic_set_free(bset);
  305. isl_aff_free(obj);
  306. return res;
  307. }
  308. /* Return the minimal (rational) value of "obj" over "bset".
  309. *
  310. * Return negative infinity if the minimal value is unbounded and
  311. * NaN if "bset" is empty.
  312. */
  313. __isl_give isl_val *isl_basic_set_min_lp_val(__isl_keep isl_basic_set *bset,
  314. __isl_keep isl_aff *obj)
  315. {
  316. return isl_basic_set_opt_lp_val(bset, 0, obj);
  317. }
  318. /* Return the maximal (rational) value of "obj" over "bset".
  319. *
  320. * Return infinity if the maximal value is unbounded and
  321. * NaN if "bset" is empty.
  322. */
  323. __isl_give isl_val *isl_basic_set_max_lp_val(__isl_keep isl_basic_set *bset,
  324. __isl_keep isl_aff *obj)
  325. {
  326. return isl_basic_set_opt_lp_val(bset, 1, obj);
  327. }