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