first experimental version of gurobi ILP solver
[libfirm] / ir / lpp / lpp_gurobi.c
1 /**
2  * Author:      Matthias Braun
3  * Copyright:   (c) Universitaet Karlsruhe
4  * Licence:     This file protected by GPL -  GNU GENERAL PUBLIC LICENSE.
5  */
6 #include "config.h"
7
8 #include "lpp_gurobi.h"
9
10 #include <stdio.h>
11 #include <stdlib.h>
12
13 #ifdef WITH_GUROBI
14
15
16 #ifdef _WIN32
17 #include <malloc.h>
18 #else
19 #include <sys/time.h>
20 #include <alloca.h>
21 #endif
22 #include <assert.h>
23
24 #include "obst.h"
25
26 #include <gurobi_c.h>
27
28 #include "error.h"
29 #include "sp_matrix.h"
30
31 static char gurobi_cst_encoding[4] = { 0, GRB_EQUAL, GRB_LESS_EQUAL, GRB_GREATER_EQUAL };
32 static char gurobi_var_encoding[4] = { 0, 0, GRB_CONTINUOUS, GRB_BINARY };
33
34 #define my_timersub(tvp, uvp, vvp)                     \
35     do {                                \
36         (vvp)->tv_sec = (tvp)->tv_sec - (uvp)->tv_sec;      \
37         (vvp)->tv_usec = (tvp)->tv_usec - (uvp)->tv_usec;   \
38         if ((vvp)->tv_usec < 0) {               \
39             (vvp)->tv_sec--;                \
40             (vvp)->tv_usec += 1000000;          \
41         }                           \
42     } while (0)
43
44 typedef struct _gurobi_t {
45         lpp_t *lpp;
46         GRBenv *env;
47         GRBmodel *model;
48 } gurobi_t;
49
50 static void check_gurobi_error(gurobi_t *grb, int error)
51 {
52         if (error != 0) {
53                 panic("gurobi error: %s", GRBgeterrormsg(grb->env));
54         }
55 }
56
57 static gurobi_t *new_gurobi(lpp_t *lpp)
58 {
59         int error;
60
61         gurobi_t *grb = XMALLOCZ(gurobi_t);
62         grb->lpp = lpp;
63         error = GRBloadenv(&grb->env, NULL);
64         check_gurobi_error(grb, error);
65         error = GRBsetlogfile(grb->env, lpp->log);
66         check_gurobi_error(grb, error);
67
68         return grb;
69 }
70
71 static void free_gurobi(gurobi_t *grb)
72 {
73         GRBfreeenv(grb->env);
74         free(grb);
75 }
76
77 /**
78  * Build CPLEX data structure from LPP matrix.
79  * @note: The LPP matrix is freed after this step, to save memory.
80  */
81 static void gurobi_construct(gurobi_t *grb)
82 {
83         const matrix_elem_t *elem;
84         int                  i, o, sv_cnt;
85         int                  numcols, numrows, numentries;
86         int                  objsen, *matbeg, *matcnt, *matind, *indices;
87         double               *obj, *rhs, *matval, *lb, *startv;
88         char                 *sense, *vartype;
89         char                 **colname, **rowname;
90         struct obstack       obst;
91         lpp_t                *lpp = grb->lpp;
92         int                  error;
93
94         numcols    = lpp->var_next-1;
95         numrows    = lpp->cst_next-1;
96         numentries = matrix_get_entries(lpp->m);
97         objsen     = lpp->opt_type == lpp_minimize ? 1 : -1;
98         obstack_init(&obst);
99
100         obj     = obstack_alloc(&obst, numcols * sizeof(*obj));
101         lb      = obstack_alloc(&obst, numcols * sizeof(*lb));
102         colname = obstack_alloc(&obst, numcols * sizeof(*colname));
103         rowname = obstack_alloc(&obst, numrows * sizeof(*rowname));
104         vartype = obstack_alloc(&obst, numcols * sizeof(*vartype));
105         indices = obstack_alloc(&obst, numcols * sizeof(*indices));
106         startv  = obstack_alloc(&obst, numcols * sizeof(*startv));
107         matbeg  = obstack_alloc(&obst, numcols * sizeof(*matbeg));
108         matcnt  = obstack_alloc(&obst, numcols * sizeof(*matcnt));
109         matind  = obstack_alloc(&obst, numentries * sizeof(*matind));
110         matval  = obstack_alloc(&obst, numentries * sizeof(*matval));
111         rhs     = obstack_alloc(&obst, numrows * sizeof(*rhs));
112         sense   = obstack_alloc(&obst, numrows * sizeof(*sense));
113
114         o      = 0;
115         sv_cnt = 0;
116         /* fill the CPLEX matrix*/
117         for (i = 0; i < numcols; ++i) {
118                 lpp_name_t *curr_var = lpp->vars[1+i];
119
120                 obj[i] = matrix_get(lpp->m, 0, 1+i);
121                 lb[i]  = 0.0;
122
123                 colname[i] = (char*) curr_var->name;
124                 vartype[i] = gurobi_var_encoding[curr_var->type.var_type];
125
126 #if 0
127                 if (curr_var->value_kind == lpp_value_start) {
128                         panic("start values not supported in gurobi yet");
129                         indices[sv_cnt]  = i;
130                         startv[sv_cnt++] = curr_var->value;
131                 }
132 #endif
133
134                 matbeg[i] = o;
135                 matcnt[i] = 0;
136                 matrix_foreach_in_col(lpp->m, 1 + i, elem) {
137                         if (elem->row == 0)
138                                 continue;
139                         matind[o] = elem->row-1;
140                         matval[o] = elem->val;
141                         matcnt[i]++;
142                         o++;
143                 }
144         }
145
146         /* get constraint stuff (right hand side, type, name) */
147         for (i = 0; i < numrows; ++i) {
148                 lpp_name_t *curr_cst = lpp->csts[1 + i];
149
150                 rhs[i]     = matrix_get(lpp->m, 1 + i, 0);
151                 sense[i]   = gurobi_cst_encoding[curr_cst->type.cst_type];
152                 rowname[i] = (char*) curr_cst->name;
153         }
154
155         error = GRBloadmodel(grb->env, &grb->model, lpp->name, numcols, numrows,
156                              objsen, 0, obj, sense, rhs, matbeg, matcnt, matind,
157                              matval, lb, NULL, vartype, colname, rowname);
158         check_gurobi_error(grb, error);
159
160         obstack_free(&obst, NULL);
161         free_lpp_matrix(lpp);
162 }
163
164 static void gurobi_solve(gurobi_t *grb)
165 {
166         lpp_t *lpp = grb->lpp;
167         int i;
168         int optimstatus;
169         int error;
170         int numcols = lpp->var_next-1;
171         double *values;
172         struct timeval tvb, tva, tvdiff;
173         double iterations;
174
175         /* set performance parameters */
176         // CPXsetintparam(grb->env, CPX_PARAM_MIPSTART, CPX_ON);
177         //CPXsetintparam(grb->env, CPX_PARAM_MIPORDTYPE, CPX_MIPORDER_COST);
178         /* output every search tree node */
179         // CPXsetintparam(grb->env, CPX_PARAM_MIPINTERVAL, 1);
180
181         /* experimental switches */
182         // CPXsetintparam(grb->env, CPX_PARAM_VARSEL, CPX_VARSEL_STRONG);
183         // CPXsetdblparam(grb->env, CPX_PARAM_BTTOL, 1.0);
184         // CPXsetintparam(grb->env, CPX_PARAM_BRDIR, CPX_BRDIR_UP);
185
186         /* Set the time limit appropriately */
187         if(lpp->time_limit_secs > 0.0) {
188                 error = GRBsetdblparam(grb->env, GRB_DBL_PAR_TIMELIMIT, lpp->time_limit_secs);
189                 check_gurobi_error(grb, error);
190         }
191
192         /*
193          * If we have enough time, we instruct cplex to imply some
194          * of its higher order magic to pursue the best solution
195          */
196         if(lpp->emphasis) {
197                 /* not implemented */
198         }
199
200         /*
201          * If a bound of the objective function is supplied,
202          * set it accordingly, dependign on minimization or maximization.
203          */
204         if(lpp->set_bound) {
205                 //panic("bound not implemented yet");
206                 fprintf(stderr, "Warning: gurobi bound not implemented yet\n");
207         }
208
209         /* solve */
210         gettimeofday(&tvb, NULL);
211         error = GRBoptimize(grb->model);
212         check_gurobi_error(grb, error);
213         gettimeofday(&tva, NULL);
214
215         /* get solution status */
216         error = GRBgetintattr(grb->model, GRB_INT_ATTR_STATUS, &optimstatus);
217         check_gurobi_error(grb, error);
218
219         switch (optimstatus) {
220         case GRB_OPTIMAL:           lpp->sol_state = lpp_optimal; break;
221         case GRB_INFEASIBLE:        lpp->sol_state = lpp_infeasible; break;
222         case GRB_INF_OR_UNBD:       lpp->sol_state = lpp_inforunb; break;
223         case GRB_UNBOUNDED:         lpp->sol_state = lpp_unbounded; break;
224         /* TODO: is this correct? */
225         default:                    lpp->sol_state = lpp_feasible; break;
226         }
227
228         /* get variable solution values */
229         values = alloca(numcols * sizeof(*values));
230         error = GRBgetdblattrarray(grb->model, GRB_DBL_ATTR_X, 0, numcols, values);
231         check_gurobi_error(grb, error);
232         for(i=0; i<numcols; ++i) {
233                 lpp->vars[1+i]->value      = values[i];
234                 lpp->vars[1+i]->value_kind = lpp_value_solution;
235         }
236
237         /* Get the value of the objective function. */
238         error = GRBgetdblattr(grb->model, GRB_DBL_ATTR_OBJVAL, &lpp->objval);
239         check_gurobi_error(grb, error);
240         error = GRBgetdblattr(grb->model , GRB_DBL_ATTR_OBJBOUND, &lpp->best_bound);
241         check_gurobi_error(grb, error);
242
243         /* get some statistics */
244         my_timersub(&tva, &tvb, &tvdiff);
245         lpp->sol_time = tvdiff.tv_sec + tvdiff.tv_usec / 1e6;
246         error = GRBgetdblattr(grb->model, GRB_DBL_ATTR_ITERCOUNT, &iterations);
247         check_gurobi_error(grb, error);
248         lpp->iterations = (unsigned) iterations;
249 }
250
251 void lpp_solve_gurobi(lpp_t *lpp)
252 {
253         gurobi_t *grb = new_gurobi(lpp);
254         gurobi_construct(grb);
255         gurobi_solve(grb);
256         free_gurobi(grb);
257 }
258
259 #endif