beifg: Simplify the quite complicated way to divide a number by 2 in be_ifg_stat().
[libfirm] / ir / lpp / lpp_gurobi.c
1 /*
2  * This file is part of libFirm.
3  * Copyright (C) 2012 University of Karlsruhe.
4  */
5
6 /**
7  * @file
8  * @author  Matthias Braun
9  */
10 #include "config.h"
11
12 #ifdef WITH_GUROBI
13 #include "lpp_gurobi.h"
14
15 #include <stdio.h>
16 #include <stdlib.h>
17 #include <math.h>
18
19 #include "obst.h"
20
21 #include <gurobi_c.h>
22
23 #include "error.h"
24 #include "sp_matrix.h"
25
26 static char gurobi_cst_encoding[4] = { 0, GRB_EQUAL, GRB_LESS_EQUAL, GRB_GREATER_EQUAL };
27 static char gurobi_var_encoding[4] = { 0, 0, GRB_CONTINUOUS, GRB_BINARY };
28
29 typedef struct _gurobi_t {
30         lpp_t *lpp;
31         GRBenv *env;
32         GRBenv *modelenv;
33         GRBmodel *model;
34 } gurobi_t;
35
36 static void check_gurobi_error(gurobi_t *grb, int error)
37 {
38         if (error != 0) {
39                 panic("gurobi error: %s", GRBgeterrormsg(grb->env));
40         }
41 }
42
43 static gurobi_t *new_gurobi(lpp_t *lpp)
44 {
45         int error;
46
47         gurobi_t *grb = XMALLOCZ(gurobi_t);
48         grb->lpp = lpp;
49         /* /tmp/firm_gurobi.log is a hack (see below) */
50         error = GRBloadenv(&grb->env, "/tmp/firm_gurobi.log");
51         check_gurobi_error(grb, error);
52         /* Matze: do not set the FILE* for logging output. Because:
53          *  a) the function is deprecated
54          *  b) gurobi closes the FILE handle when it is done, which leads to
55          *     very unexpected effects when you pass stdout or stderr as logging
56          *     output.
57          * The only thing gurobi sanely supports is giving a string with a filename
58          * :-( ...so we use /tmp/firm_gurobi.log as a temporary measure...
59          */
60         if (lpp->log != stdout && lpp->log != stderr) {
61                 error = GRBsetintparam(grb->env, GRB_INT_PAR_OUTPUTFLAG, 0);
62                 check_gurobi_error(grb, error);
63         }
64
65         return grb;
66 }
67
68 static void free_gurobi(gurobi_t *grb)
69 {
70         GRBfreemodel(grb->model);
71         GRBfreeenv(grb->env);
72         free(grb);
73 }
74
75 /**
76  * Build CPLEX data structure from LPP matrix.
77  * @note: The LPP matrix is freed after this step, to save memory.
78  */
79 static void gurobi_construct(gurobi_t *grb)
80 {
81         int            i, o;
82         //int            sv_cnt;
83         //int           *indices;
84         //double        *startv;
85         int            numcols, numrows, numentries;
86         int            objsen, *matbeg, *matcnt, *matind;
87         double        *obj, *rhs, *matval, *lb;
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                 matbeg[i] = o;
127                 matcnt[i] = 0;
128                 matrix_foreach_in_col(lpp->m, 1 + i, elem) {
129                         if (elem->row == 0)
130                                 continue;
131                         matind[o] = elem->row-1;
132                         matval[o] = elem->val;
133                         matcnt[i]++;
134                         o++;
135                 }
136         }
137
138         /* get constraint stuff (right hand side, type, name) */
139         for (i = 0; i < numrows; ++i) {
140                 lpp_name_t *curr_cst = lpp->csts[1 + i];
141
142                 rhs[i]     = matrix_get(lpp->m, 1 + i, 0);
143                 sense[i]   = gurobi_cst_encoding[curr_cst->type.cst_type];
144                 rowname[i] = (char*) curr_cst->name;
145         }
146
147         error = GRBloadmodel(grb->env, &grb->model, lpp->name, numcols, numrows,
148                              objsen, 0, obj, sense, rhs, matbeg, matcnt, matind,
149                              matval, lb, NULL, vartype, colname, rowname);
150         check_gurobi_error(grb, error);
151         grb->modelenv = GRBgetenv(grb->model);
152
153         obstack_free(&obst, NULL);
154         lpp_free_matrix(lpp);
155 }
156
157 static void gurobi_solve(gurobi_t *grb)
158 {
159         lpp_t *lpp = grb->lpp;
160         int i;
161         int optimstatus;
162         int error;
163         int numcols = lpp->var_next-1;
164         double *values;
165         double  iterations;
166
167         /* Set the time limit appropriately */
168         if(lpp->time_limit_secs > 0.0) {
169                 error = GRBsetdblparam(grb->modelenv, GRB_DBL_PAR_TIMELIMIT, lpp->time_limit_secs);
170                 check_gurobi_error(grb, error);
171         }
172
173         /* solve */
174         error = GRBoptimize(grb->model);
175         check_gurobi_error(grb, error);
176
177         /* get solution status */
178         error = GRBgetintattr(grb->model, GRB_INT_ATTR_STATUS, &optimstatus);
179         check_gurobi_error(grb, error);
180
181         switch (optimstatus) {
182         case GRB_OPTIMAL:           lpp->sol_state = lpp_optimal; break;
183         case GRB_INFEASIBLE:        lpp->sol_state = lpp_infeasible; break;
184         case GRB_INF_OR_UNBD:       lpp->sol_state = lpp_inforunb; break;
185         case GRB_UNBOUNDED:         lpp->sol_state = lpp_unbounded; break;
186         /* TODO: is this correct? */
187         default:                    lpp->sol_state = lpp_feasible; break;
188         }
189
190         if (lpp->sol_state >= lpp_feasible) {
191                 /* get variable solution values */
192                 values = alloca(numcols * sizeof(*values));
193                 error = GRBgetdblattrarray(grb->model, GRB_DBL_ATTR_X, 0, numcols,
194                                            values);
195                 check_gurobi_error(grb, error);
196                 for(i=0; i<numcols; ++i) {
197                         lpp->vars[1+i]->value      = values[i];
198                         lpp->vars[1+i]->value_kind = lpp_value_solution;
199                 }
200
201                 /* Get the value of the objective function. */
202                 error = GRBgetdblattr(grb->model, GRB_DBL_ATTR_OBJVAL, &lpp->objval);
203                 check_gurobi_error(grb, error);
204                 error = GRBgetdblattr(grb->model , GRB_DBL_ATTR_OBJBOUND,
205                                       &lpp->best_bound);
206                 if (error != 0) {
207                         lpp->best_bound = FP_NAN;
208                 }
209         }
210
211         /* get some statistics */
212         error = GRBgetdblattr(grb->model, GRB_DBL_ATTR_ITERCOUNT, &iterations);
213         check_gurobi_error(grb, error);
214         lpp->iterations = (unsigned) iterations;
215
216         error = GRBgetdblattr(grb->model, GRB_DBL_ATTR_RUNTIME, &lpp->sol_time);
217         check_gurobi_error(grb, error);
218 }
219
220 void lpp_solve_gurobi(lpp_t *lpp)
221 {
222         gurobi_t *grb = new_gurobi(lpp);
223         gurobi_construct(grb);
224         gurobi_solve(grb);
225         free_gurobi(grb);
226 }
227
228 #endif