isl_lp_piplib.c

   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
  10 #include <isl/map.h>
  11 #include <isl/vec.h>
  12 #include <isl/lp.h>
  13 #include "isl_piplib.h"
  14 #include "isl_map_piplib.h"
  15
  16 static void copy_solution(struct isl_vec *vec, int maximize, isl_int *opt,
  17         isl_int *opt_denom, PipQuast *sol)
  18 {
  19         int i;
  20         PipList *list;
  21         isl_int tmp;
  22
  23         if (opt) {
  24                 if (opt_denom) {
  25                         isl_seq_cpy_from_pip(opt,
  26                                  &sol->list->vector->the_vector[0], 1);
  27                         isl_seq_cpy_from_pip(opt_denom,
  28                                  &sol->list->vector->the_deno[0], 1);
  29                 } else if (maximize)
  30                         mpz_fdiv_q(*opt, sol->list->vector->the_vector[0],
  31                                          sol->list->vector->the_deno[0]);
  32                 else
  33                         mpz_cdiv_q(*opt, sol->list->vector->the_vector[0],
  34                                          sol->list->vector->the_deno[0]);
  35         }
  36
  37         if (!vec)
  38                 return;
  39
  40         isl_int_init(tmp);
  41         isl_int_set_si(vec->el[0], 1);
  42         for (i = 0, list = sol->list->next; list; ++i, list = list->next) {
  43                 isl_seq_cpy_from_pip(&vec->el[1 + i],
  44                         &list->vector->the_deno[0], 1);
  45                 isl_int_lcm(vec->el[0], vec->el[0], vec->el[1 + i]);
  46         }
  47         for (i = 0, list = sol->list->next; list; ++i, list = list->next) {
  48                 isl_seq_cpy_from_pip(&tmp, &list->vector->the_deno[0], 1);
  49                 isl_int_divexact(tmp, vec->el[0], tmp);
  50                 isl_seq_cpy_from_pip(&vec->el[1 + i],
  51                         &list->vector->the_vector[0], 1);
  52                 isl_int_mul(vec->el[1 + i], vec->el[1 + i], tmp);
  53         }
  54         isl_int_clear(tmp);
  55 }
  56
  57 enum isl_lp_result isl_pip_solve_lp(struct isl_basic_map *bmap, int maximize,
  58                                       isl_int *f, isl_int denom, isl_int *opt,
  59                                       isl_int *opt_denom,
  60                                       struct isl_vec **vec)
  61 {
  62         enum isl_lp_result res = isl_lp_ok;
  63         PipMatrix       *domain = NULL;
  64         PipOptions      *options;
  65         PipQuast        *sol;
  66         unsigned         total;
  67
  68         total = isl_basic_map_total_dim(bmap);
  69         domain = isl_basic_map_to_pip(bmap, 0, 1, 0);
  70         if (!domain)
  71                 goto error;
  72         entier_set_si(domain->p[0][1], -1);
  73         isl_int_set(domain->p[0][domain->NbColumns - 1], f[0]);
  74         isl_seq_cpy_to_pip(domain->p[0]+2, f+1, total);
  75
  76         options = pip_options_init();
  77         if (!options)
  78                 goto error;
  79         options->Urs_unknowns = -1;
  80         options->Maximize = maximize;
  81         options->Nq = 0;
  82         sol = pip_solve(domain, NULL, -1, options);
  83         pip_options_free(options);
  84         if (!sol)
  85                 goto error;
  86
  87         if (vec) {
  88                 isl_ctx *ctx = isl_basic_map_get_ctx(bmap);
  89                 *vec = isl_vec_alloc(ctx, 1 + total);
  90         }
  91         if (vec && !*vec)
  92                 res = isl_lp_error;
  93         else if (!sol->list)
  94                 res = isl_lp_empty;
  95         else if (entier_zero_p(sol->list->vector->the_deno[0]))
  96                 res = isl_lp_unbounded;
  97         else
  98                 copy_solution(*vec, maximize, opt, opt_denom, sol);
  99         pip_matrix_free(domain);
 100         pip_quast_free(sol);
 101         return res;
 102 error:
 103         if (domain)
 104                 pip_matrix_free(domain);
 105         return isl_lp_error;
 106 }