diff --git a/src/lib/glpk.ml b/src/lib/glpk.ml
new file mode 100644
index 0000000000000000000000000000000000000000..1dee50fadefd530c0ef8717aba8823fdf12a6834
--- /dev/null
+++ b/src/lib/glpk.ml
@@ -0,0 +1,188 @@
+(*
+ * ocaml-glpk - OCaml bindings to glpk
+ * Copyright (C) 2004-2006 Samuel Mimram, 2014 Dominik Paulus
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place - Suite 330,
+ * Boston, MA 02111-1307, USA.
+ *)
+
+(* $Id$ *)
+
+type lp
+
+type direction = Minimize | Maximize
+
+type aux_var_type = Free_var | Lower_bounded_var | Upper_bounded_var | Double_bounded_var | Fixed_var
+
+type prob_class = Linear_prog | Mixed_integer_prog
+
+type var_kind = Continuous_var | Integer_var
+
+exception Fault
+exception Lower_limit
+exception Upper_limit
+exception No_primal_feasible_solution
+exception No_dual_feasible_solution
+exception Iteration_limit
+exception Time_limit
+exception Solver_failure
+exception Empty
+exception Bad_basis
+exception No_convergence
+exception Unknown_error
+
+let _ =
+  Callback.register_exception "ocaml_glpk_exn_fault" Fault;
+  Callback.register_exception "ocaml_glpk_exn_objll" Lower_limit;
+  Callback.register_exception "ocaml_glpk_exn_objul" Upper_limit;
+  Callback.register_exception "ocaml_glpk_exn_nopfs" No_primal_feasible_solution;
+  Callback.register_exception "ocaml_glpk_exn_nodfs" No_dual_feasible_solution;
+  Callback.register_exception "ocaml_glpk_exn_itlim" Iteration_limit;
+  Callback.register_exception "ocaml_glpk_exn_tmlim" Time_limit;
+  Callback.register_exception "ocaml_glpk_exn_sing" Solver_failure;
+  Callback.register_exception "ocaml_glpk_exn_empty" Empty;
+  Callback.register_exception "ocaml_glpk_exn_badb" Bad_basis;
+  Callback.register_exception "ocaml_glpk_exn_noconv" No_convergence;
+  Callback.register_exception "ocaml_glpk_exn_unkown" Unknown_error;
+
+external new_problem : unit -> lp = "ocaml_glpk_new_prob"
+
+external set_prob_name : lp -> string -> unit = "ocaml_glpk_set_prob_name"
+
+external get_prob_name : lp -> string = "ocaml_glpk_get_prob_name"
+
+external set_obj_name : lp -> string -> unit = "ocaml_glpk_set_obj_name"
+
+external get_obj_name : lp -> string = "ocaml_glpk_get_obj_name"
+
+external set_direction : lp -> direction -> unit = "ocaml_glpk_set_direction"
+
+external get_direction : lp -> direction = "ocaml_glpk_get_direction"
+
+external add_rows : lp -> int -> unit = "ocaml_glpk_add_rows"
+
+external set_row_name : lp -> int -> string -> unit = "ocaml_glpk_set_row_name"
+
+external get_row_name : lp -> int -> string = "ocaml_glpk_get_row_name"
+
+external set_row_bounds : lp -> int -> aux_var_type -> float -> float -> unit = "ocaml_glpk_set_row_bounds"
+
+external add_columns : lp -> int -> unit = "ocaml_glpk_add_cols"
+
+external set_col_name : lp -> int -> string -> unit = "ocaml_glpk_set_col_name"
+
+external get_col_name : lp -> int -> string = "ocaml_glpk_get_col_name"
+
+external set_col_bounds : lp -> int -> aux_var_type -> float -> float -> unit = "ocaml_glpk_set_col_bounds"
+
+external set_obj_coef : lp -> int -> float -> unit = "ocaml_glpk_set_obj_coef"
+
+external load_matrix : lp -> float array array -> unit = "ocaml_glpk_load_matrix"
+
+external load_sparse_matrix : lp -> ((int * int) * float) array -> unit = "ocaml_glpk_load_sparse_matrix"
+
+external simplex : lp -> unit = "ocaml_glpk_simplex"
+
+external get_obj_val : lp -> float = "ocaml_glpk_get_obj_val"
+
+external get_col_primal : lp -> int -> float = "ocaml_glpk_get_col_prim"
+
+external get_row_primal : lp -> int -> float = "ocaml_glpk_get_row_prim"
+
+external get_row_dual : lp -> int -> float = "ocaml_glpk_get_row_dual"
+
+let make_problem dir zcoefs constr pbounds xbounds =
+  let lp = new_problem () in
+    set_direction lp dir;
+    add_rows lp (Array.length pbounds);
+    for i = 0 to (Array.length pbounds) - 1 do
+      match pbounds.(i) with
+        | lb, ub when lb = -.infinity && ub = infinity -> set_row_bounds lp i Free_var 0. 0.
+        | lb, ub when ub = infinity -> set_row_bounds lp i Lower_bounded_var lb 0.
+        | lb, ub when lb = -.infinity -> set_row_bounds lp i Upper_bounded_var 0. ub
+        | lb, ub when lb = ub -> set_row_bounds lp i Fixed_var lb ub
+        | lb, ub -> set_row_bounds lp i Double_bounded_var lb ub
+    done;
+    add_columns lp (Array.length xbounds);
+    for i = 0 to (Array.length xbounds) - 1 do
+      set_obj_coef lp i zcoefs.(i);
+      match xbounds.(i) with
+        | lb, ub when lb = -.infinity && ub = infinity -> set_col_bounds lp i Free_var 0. 0.
+        | lb, ub when ub = infinity -> set_col_bounds lp i Lower_bounded_var lb 0.
+        | lb, ub when lb = -.infinity -> set_col_bounds lp i Upper_bounded_var 0. ub
+        | lb, ub when lb = ub -> set_col_bounds lp i Fixed_var lb ub
+        | lb, ub -> set_col_bounds lp i Double_bounded_var lb ub
+    done;
+    load_matrix lp constr;
+    lp
+
+external get_num_rows : lp -> int = "ocaml_glpk_get_num_rows"
+
+external get_num_cols : lp -> int = "ocaml_glpk_get_num_cols"
+
+let get_col_primals lp =
+  let n = get_num_cols lp in
+  let ans = Array.make n 0. in
+    for i = 0 to (n - 1)
+    do
+      ans.(i) <- get_col_primal lp i
+    done;
+    ans
+
+external scale_problem : lp -> unit = "ocaml_glpk_scale_problem"
+
+external unscale_problem : lp -> unit = "ocaml_glpk_unscale_problem"
+
+external interior : lp -> unit = "ocaml_glpk_interior"
+
+external set_class : lp -> prob_class -> unit = "ocaml_glpk_set_class"
+
+external get_class : lp -> prob_class = "ocaml_glpk_get_class"
+
+external set_col_kind : lp -> int -> var_kind -> unit = "ocaml_glpk_set_col_kind"
+
+external branch_and_bound : lp -> unit = "ocaml_glpk_integer"
+
+external branch_and_bound_opt : lp -> unit = "ocaml_glpk_integer"
+
+external warm_up : lp -> unit = "ocaml_glpk_warm_up"
+
+external use_presolver : lp -> bool -> unit = "ocaml_glpk_set_use_presolver"
+
+external read_cplex : string -> lp = "ocaml_glpk_read_cplex"
+
+external write_cplex : lp -> string -> unit = "ocaml_glpk_write_cplex"
+
+external set_simplex_iteration_count : lp -> int -> unit = "ocaml_glpk_set_iteration_count"
+
+let reset_simplex_iteration_count lp =
+  set_simplex_iteration_count lp 0
+
+external get_simplex_iteration_count : lp -> int = "ocaml_glpk_get_iteration_count"
+
+external _set_message_level : lp -> int -> unit = "ocaml_glpk_set_message_level"
+
+let set_message_level lp n =
+    if (n < 0 && n > 3) then
+        raise (Invalid_argument "set_message_level");
+    _set_message_level lp n
+
+external set_simplex_iteration_limit : lp -> int -> unit = "ocaml_glpk_set_iteration_limit"
+
+external get_simplex_iteration_limit : lp -> int = "ocaml_glpk_get_iteration_limit"
+
+external set_simplex_time_limit : lp -> float -> unit = "ocaml_glpk_set_time_limit"
+
+external get_simplex_time_limit : lp -> float = "ocaml_glpk_get_time_limit"
diff --git a/src/lib/glpk.mli b/src/lib/glpk.mli
new file mode 100644
index 0000000000000000000000000000000000000000..84153e7efadb235c950d56f21e45c6ef03e69a91
--- /dev/null
+++ b/src/lib/glpk.mli
@@ -0,0 +1,234 @@
+(**
+  * OCaml bindings to glpk. Please see the glpk manual for further explanations 
+  * on the semantics of functions.
+  *
+  * Warning: contrarily to the C version of glpk, all indexes are 0-based.
+  *
+  * @author Samuel Mimram
+  *)
+
+
+(* $Id$ *)
+
+
+(** {1 Types} *) (* TODO: better comment! *)
+
+(** A linear programmation problem. *)
+type lp
+
+(** Direction of the optimization. *)
+type direction = Minimize | Maximize
+
+(** Type of bounds of an auxiliary variable. *)
+type aux_var_type = Free_var | Lower_bounded_var | Upper_bounded_var | Double_bounded_var | Fixed_var
+
+(** Class of a problem. *)
+type prob_class =
+  | Linear_prog (** linear programming *)
+  | Mixed_integer_prog (** mixed integer programming *)
+
+(** Kind of a variable. *)
+type var_kind =
+  | Continuous_var (** continuous variable *)
+  | Integer_var (** integer variable *)
+
+(** {1 Exceptions} *)
+
+(** The problem has no rows/columns, or the initial basis is invalid, or the initial basis matrix is singular or ill-conditionned. *)
+exception Fault
+
+(** The problem has no rows and/or column. *)
+exception Empty
+
+(** The LP basis is invalid beacause the number of basic variables is not the same as the number of rows. *)
+exception Bad_basis
+
+(** The objective function being minimized has reached its lower limit and continues decreasing. *)
+exception Lower_limit
+
+(** The objective function being maximized has reached its upper limit and continues increasing. *)
+exception Upper_limit
+
+(** The problem has no primal feasible solution. *)
+exception No_primal_feasible_solution
+
+(** The problem has no dual feasible solution. *)
+exception No_dual_feasible_solution
+
+(** Iterations limit exceeded. *)
+exception Iteration_limit
+
+(** Time limit exceeded. *)
+exception Time_limit
+
+(** Very slow convergence or divergence. *)
+exception No_convergence
+
+(** Failure of the solver (the current basis matrix got singular or ill-conditionned). *)
+exception Solver_failure
+
+(** Unknown error (this exception should disappear in future versions). *)
+exception Unknown_error
+
+
+(** {1 Functions} *)
+
+(** {2 Creating, reading and saving problems} *)
+
+(** Create a new linear programmation problem. *)
+val new_problem : unit -> lp
+
+(** [make_problem dir zcoefs constrs pbounds xbounds] creates the new linear programmation problem where Z = Sum_i [zcoefs.(i)] * x_ i should be optimized in the direction [dir] under the constraints [fst pbounds.(i)] <= p_i <= [snd pbounds.(i)] and [fst xbounds.(i)] <= x_i <= [snd xbounds.(i)] where p_i = Sum_j [constrs.(i).(j)] * x_j. The bounds may be [+] / [- infinity]. *)
+val make_problem : direction -> float array -> float array array -> (float * float) array -> (float * float) array -> lp
+
+(** Read problem data in CPLEX LP format from a file. *)
+val read_cplex : string -> lp
+
+(** Write prblem data in CPLEX LP format into a file. *)
+val write_cplex : lp -> string -> unit
+
+
+(** {2 Setting and retreiving paramters of a problem} *)
+
+(** Set the problem name. *)
+val set_prob_name : lp -> string -> unit
+
+(** Retrieve the problem name. *)
+val get_prob_name : lp -> string
+
+(** Set the problem class. *)
+val set_class : lp -> prob_class -> unit
+
+(** Retrieve the problem class. *)
+val get_class : lp -> prob_class
+
+(** Set the direction of the optimization. *)
+val set_direction : lp -> direction -> unit
+
+(** Retrieve the direction of the optimization. *)
+val get_direction : lp -> direction
+
+(** Set the objective name. *)
+val set_obj_name : lp -> string -> unit
+
+(** Retrieve the objective name. *)
+val get_obj_name : lp -> string
+
+(** Add rows. *)
+val add_rows : lp -> int -> unit
+
+(** Retreive the number of rows. *)
+val get_num_rows : lp -> int
+
+(** Set the name of a row. *)
+val set_row_name : lp -> int -> string -> unit
+
+(** Retrieve the name of a row. *)
+val get_row_name : lp -> int -> string
+
+(** Set a row bound. *)
+val set_row_bounds : lp -> int -> aux_var_type -> float -> float -> unit
+
+(** Add columns. *)
+val add_columns : lp -> int -> unit
+
+(** Retreive the number of columns. *)
+val get_num_cols : lp -> int
+
+(** Set the name of a column. *)
+val set_col_name : lp -> int -> string -> unit
+
+(** Retrieve the name of a column. *)
+val get_col_name : lp -> int -> string
+
+(** Set column kind. *)
+val set_col_kind : lp -> int -> var_kind -> unit
+
+(** Set a column boudaries. *)
+val set_col_bounds : lp -> int -> aux_var_type -> float -> float -> unit
+
+(** Set an objective coefficient. *)
+val set_obj_coef : lp -> int -> float -> unit
+
+(** Load a constraint matrix. *)
+val load_matrix : lp -> float array array -> unit
+
+(** Load a sparse constraint matrix stored as an array whose elements are of the
+  * form ((row, column), value) indicating non-null elements of the matrix. *)
+val load_sparse_matrix : lp -> ((int * int) * float) array -> unit
+
+
+(** {2 Solving problems and retreiving solutions} *)
+
+(** Scale problem data. *)
+val scale_problem : lp -> unit
+
+(** Unscale problem data. *)
+val unscale_problem : lp -> unit
+
+(** Warm up the LP basis for the specified problem object using current statuses assigned to rows and columns. *)
+val warm_up : lp -> unit
+
+(** Solve an LP problem using the simplex method. You must use builtin presolver
+  * (see [use_presolver]) to get an exception if the problem has no feasible
+  * solution. *)
+val simplex : lp -> unit
+
+(** Solve an LP problem using the primal-dual interior point method. *)
+val interior : lp -> unit
+
+(** Solve a MIP proble using the branch-and-bound method. *)
+val branch_and_bound : lp -> unit
+
+(** Solve a MIP proble using and optimized version of the branch-and-bound method. *)
+val branch_and_bound_opt : lp -> unit
+
+(** Retrieve objective value. *)
+val get_obj_val : lp -> float
+
+(** Get the primal value of the structural variable associated with a column. *)
+val get_col_primal : lp -> int -> float
+
+(** Get the primal values of the structural variables associated with each column. *)
+val get_col_primals : lp -> float array
+
+(** Get the primal value of the structural variable associated with a row. *)
+val get_row_primal : lp -> int -> float
+
+(** Get the dual value of the structural variable associated with a row. *)
+val get_row_dual : lp -> int -> float
+
+
+(** {2 Setting parameters of the solver} *)
+
+(** Set the level of messages output by sover routines. The second argument might be:
+  - 0: no output
+  - 1: error message only
+  - 2: normal output
+  - 3: full output (includes informational messages)
+*)
+val set_message_level : lp -> int -> unit
+
+(** Use the builtin LP-presolver in [simplex]? *)
+val use_presolver : lp -> bool -> unit
+
+(** Initialize the simplex iteration counter. *)
+val set_simplex_iteration_count : lp -> int -> unit
+
+(** Reset the simplex iteration counter. *)
+val reset_simplex_iteration_count : lp -> unit
+
+(** This number is incremented after each simplex iteration. *)
+val get_simplex_iteration_count : lp -> int
+
+(** Set the maximum number of iterations that [simplex] should do. *)
+val set_simplex_iteration_limit : lp -> int -> unit
+
+(** Retrieve the maximum number of iterations that [simplex] should do. *)
+val get_simplex_iteration_limit : lp -> int
+
+(** Set the maximum amount of time that [simplex] should take. *)
+val set_simplex_time_limit : lp -> float -> unit
+
+(** Retrieve the maximum amount of time that [simplex] should take. *)
+val get_simplex_time_limit : lp -> float
diff --git a/src/lib/glpk_stub.c b/src/lib/glpk_stub.c
new file mode 100644
index 0000000000000000000000000000000000000000..6cd8501cba1418dc36bbbb5b9c0a55e169631b90
--- /dev/null
+++ b/src/lib/glpk_stub.c
@@ -0,0 +1,526 @@
+/*
+ * ocaml-glpk - OCaml bindings to glpk
+ * Copyright (C) 2004 Samuel Mimram
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place - Suite 330,
+ * Boston, MA 02111-1307, USA.
+ */
+
+/* $Id$ */
+
+#include <caml/alloc.h>
+#include <caml/callback.h>
+#include <caml/custom.h>
+#include <caml/fail.h>
+#include <caml/memory.h>
+#include <caml/misc.h>
+#include <caml/mlvalues.h>
+#include <caml/signals.h>
+
+#include <assert.h>
+
+#include <glpk.h>
+
+static void raise_on_error(int ret)
+{
+  switch(ret)
+  {
+    case LPX_E_OK:
+      return;
+
+    case LPX_E_FAULT:
+      caml_raise_constant(*caml_named_value("ocaml_glpk_exn_fault"));
+
+    case LPX_E_OBJLL:
+      caml_raise_constant(*caml_named_value("ocaml_glpk_exn_objll"));
+
+    case LPX_E_OBJUL:
+      caml_raise_constant(*caml_named_value("ocaml_glpk_exn_objul"));
+
+    case LPX_E_NOPFS:
+      caml_raise_constant(*caml_named_value("ocaml_glpk_exn_nopfs"));
+
+    case LPX_E_NODFS:
+      caml_raise_constant(*caml_named_value("ocaml_glpk_exn_nodfs"));
+
+    case LPX_E_ITLIM:
+      caml_raise_constant(*caml_named_value("ocaml_glpk_exn_itlim"));
+
+    case LPX_E_TMLIM:
+      caml_raise_constant(*caml_named_value("ocaml_glpk_exn_tmlim"));
+
+    case LPX_E_SING:
+      caml_raise_constant(*caml_named_value("ocaml_glpk_exn_sing"));
+
+    case LPX_E_EMPTY:
+      caml_raise_constant(*caml_named_value("ocaml_glpk_exn_empty"));
+
+    case LPX_E_BADB:
+      caml_raise_constant(*caml_named_value("ocaml_glpk_exn_badb"));
+
+    case LPX_E_NOCONV:
+      caml_raise_constant(*caml_named_value("ocaml_glpk_exn_noconv"));
+
+    default:
+      caml_raise_constant(*caml_named_value("ocaml_glpk_exn_unknown"));
+  }
+  assert(0); /* TODO */
+}
+
+#define Lpx_val(v) (*((LPX**)Data_custom_val(v)))
+
+static void finalize_lpx(value block)
+{
+  lpx_delete_prob(Lpx_val(block));
+}
+
+static struct custom_operations lpx_ops =
+{
+  "ocaml_glpk_lpx",
+  finalize_lpx,
+  custom_compare_default,
+  custom_hash_default,
+  custom_serialize_default,
+  custom_deserialize_default
+};
+
+static value new_blp(LPX* lp)
+{
+  value block = caml_alloc_custom(&lpx_ops, sizeof(LPX*), 0, 1);
+  Lpx_val(block) = lp;
+  return block;
+}
+
+CAMLprim value ocaml_glpk_new_prob(value unit)
+{
+  LPX *lp = lpx_create_prob();
+  return new_blp(lp);
+}
+
+CAMLprim value ocaml_glpk_set_prob_name(value blp, value name)
+{
+  LPX *lp = Lpx_val(blp);
+  lpx_set_prob_name(lp, String_val(name));
+  return Val_unit;
+}
+
+CAMLprim value ocaml_glpk_get_prob_name(value blp)
+{
+  CAMLparam1(blp);
+  LPX *lp = Lpx_val(blp);
+  CAMLreturn(caml_copy_string(lpx_get_prob_name(lp)));
+}
+
+CAMLprim value ocaml_glpk_set_obj_name(value blp, value name)
+{
+  LPX *lp = Lpx_val(blp);
+  lpx_set_obj_name(lp, String_val(name));
+  return Val_unit;
+}
+
+CAMLprim value ocaml_glpk_get_obj_name(value blp)
+{
+  CAMLparam1(blp);
+  LPX *lp = Lpx_val(blp);
+  CAMLreturn(caml_copy_string(lpx_get_obj_name(lp)));
+}
+
+static int direction_table[] = {LPX_MIN, LPX_MAX};
+
+CAMLprim value ocaml_glpk_set_direction(value blp, value direction)
+{
+  LPX *lp = Lpx_val(blp);
+  lpx_set_obj_dir(lp, direction_table[Int_val(direction)]);
+  return Val_unit;
+}
+
+CAMLprim value ocaml_glpk_get_direction(value blp)
+{
+  LPX *lp = Lpx_val(blp);
+  switch(lpx_get_obj_dir(lp))
+  {
+    case LPX_MIN:
+      return Val_int(0);
+
+    case LPX_MAX:
+      return Val_int(1);
+
+    default:
+      assert(0);
+  }
+}
+
+CAMLprim value ocaml_glpk_add_rows(value blp, value n)
+{
+  LPX *lp = Lpx_val(blp);
+  lpx_add_rows(lp, Int_val(n));
+  return Val_unit;
+}
+
+CAMLprim value ocaml_glpk_set_row_name(value blp, value n, value name)
+{
+  LPX *lp = Lpx_val(blp);
+  lpx_set_row_name(lp, Int_val(n) + 1, String_val(name));
+  return Val_unit;
+}
+
+CAMLprim value ocaml_glpk_get_row_name(value blp, value n)
+{
+  CAMLparam1(blp);
+  LPX *lp = Lpx_val(blp);
+  CAMLreturn(caml_copy_string(lpx_get_row_name(lp, Int_val(n) + 1)));
+}
+
+static int auxvartype_table[] = {LPX_FR, LPX_LO, LPX_UP, LPX_DB, LPX_FX};
+
+CAMLprim value ocaml_glpk_set_row_bounds(value blp, value n, value type, value lb, value ub)
+{
+  LPX *lp = Lpx_val(blp);
+  lpx_set_row_bnds(lp, Int_val(n) + 1, auxvartype_table[Int_val(type)], Double_val(lb), Double_val(ub));
+  return Val_unit;
+}
+
+CAMLprim value ocaml_glpk_add_cols(value blp, value n)
+{
+  LPX *lp = Lpx_val(blp);
+  lpx_add_cols(lp, Int_val(n));
+  return Val_unit;
+}
+
+CAMLprim value ocaml_glpk_set_col_name(value blp, value n, value name)
+{
+  LPX *lp = Lpx_val(blp);
+  lpx_set_col_name(lp, Int_val(n) + 1, String_val(name));
+  return Val_unit;
+}
+
+CAMLprim value ocaml_glpk_get_col_name(value blp, value n)
+{
+  CAMLparam1(blp);
+  LPX *lp = Lpx_val(blp);
+  CAMLreturn(caml_copy_string(lpx_get_col_name(lp, Int_val(n) + 1)));
+}
+
+CAMLprim value ocaml_glpk_set_col_bounds(value blp, value n, value type, value lb, value ub)
+{
+  LPX *lp = Lpx_val(blp);
+  lpx_set_col_bnds(lp, Int_val(n) + 1, auxvartype_table[Int_val(type)], Double_val(lb), Double_val(ub));
+  return Val_unit;
+}
+
+CAMLprim value ocaml_glpk_set_obj_coef(value blp, value n, value coef)
+{
+  LPX *lp = Lpx_val(blp);
+  lpx_set_obj_coef(lp, Int_val(n) + 1, Double_val(coef));
+  return Val_unit;
+}
+
+CAMLprim value ocaml_glpk_load_matrix(value blp, value matrix)
+{
+  LPX *lp = Lpx_val(blp);
+  int i_dim = Wosize_val(matrix), j_dim;
+  int *ia, *ja;
+  double *ar;
+  double x;
+  int i, j, n;
+
+  if (i_dim <= 0)
+    return Val_unit;
+
+  j_dim = Wosize_val(Field(matrix, 0)) / Double_wosize;
+  ia = (int*)malloc((i_dim * j_dim + 1) * sizeof(int));
+  ja = (int*)malloc((i_dim * j_dim + 1) * sizeof(int));
+  ar = (double*)malloc((i_dim * j_dim + 1) * sizeof(double));
+  n = 1;
+
+  for(i = 0; i < i_dim; i++)
+  {
+    /* TODO: raise an error */
+    assert(Wosize_val(Field(matrix, i)) / Double_wosize == j_dim);
+    for(j = 0; j < j_dim; j++)
+    {
+      x = Double_field(Field(matrix, i), j);
+      /* We only want non null elements. */
+      if (x != 0)
+      {
+        ia[n] = i + 1;
+        ja[n] = j + 1;
+        ar[n] = x;
+        n++;
+      }
+    }
+  }
+  lpx_load_matrix(lp, n - 1, ia, ja, ar);
+
+  free(ia);
+  free(ja);
+  free(ar);
+
+  return Val_unit;
+}
+
+CAMLprim value ocaml_glpk_load_sparse_matrix(value blp, value matrix)
+{
+  LPX *lp = Lpx_val(blp);
+  int len = Wosize_val(matrix);
+  int *ia, *ja;
+  double *ar;
+  int i;
+  value e;
+
+  ia = (int*)malloc((len + 1) * sizeof(int));
+  ja = (int*)malloc((len + 1) * sizeof(int));
+  ar = (double*)malloc((len + 1) * sizeof(double));
+
+  for(i = 0; i < len; i++)
+  {
+    e = Field(matrix, i);
+    ia[i+1] = Int_val(Field(Field(e, 0), 0)) + 1;
+    ja[i+1] = Int_val(Field(Field(e, 0), 1)) + 1;
+    ar[i+1] = Double_val(Field(e, 1));
+  }
+  lpx_load_matrix(lp, len, ia, ja, ar);
+
+  free(ia);
+  free(ja);
+  free(ar);
+
+  return Val_unit;
+}
+
+CAMLprim value ocaml_glpk_simplex(value blp)
+{
+  CAMLparam1(blp);
+  LPX *lp = Lpx_val(blp);
+  int ret;
+
+  caml_enter_blocking_section();
+  ret = lpx_simplex(lp);
+  caml_leave_blocking_section();
+
+  raise_on_error(ret);
+  CAMLreturn(Val_unit);
+}
+
+CAMLprim value ocaml_glpk_get_obj_val(value blp)
+{
+  LPX *lp = Lpx_val(blp);
+  double ans;
+  if (lpx_get_class(lp) == LPX_MIP)
+    ans = lpx_mip_obj_val(lp);
+  else
+    ans = lpx_get_obj_val(lp);
+  return caml_copy_double(ans);
+}
+
+CAMLprim value ocaml_glpk_get_col_prim(value blp, value n)
+{
+  LPX *lp = Lpx_val(blp);
+  double ans;
+  /* TODO: is it the right thing to do? */
+  if (lpx_get_class(lp) == LPX_MIP)
+    ans = lpx_mip_col_val(lp, Int_val(n) + 1);
+  else
+    ans = lpx_get_col_prim(lp, Int_val(n) + 1);
+  return caml_copy_double(ans);
+}
+
+CAMLprim value ocaml_glpk_get_row_prim(value blp, value n)
+{
+  LPX *lp = Lpx_val(blp);
+  return caml_copy_double(lpx_get_row_prim(lp, Int_val(n) + 1));
+}
+
+CAMLprim value ocaml_glpk_get_row_dual(value blp, value n)
+{
+  LPX *lp = Lpx_val(blp);
+  return caml_copy_double(lpx_get_row_dual(lp, Int_val(n) + 1));
+}
+
+CAMLprim value ocaml_glpk_get_num_rows(value blp)
+{
+  LPX *lp = Lpx_val(blp);
+  return Val_int(lpx_get_num_rows(lp));
+}
+
+CAMLprim value ocaml_glpk_get_num_cols(value blp)
+{
+  LPX *lp = Lpx_val(blp);
+  return Val_int(lpx_get_num_cols(lp));
+}
+
+CAMLprim value ocaml_glpk_scale_problem(value blp)
+{
+  LPX *lp = Lpx_val(blp);
+  lpx_scale_prob(lp);
+  return Val_unit;
+}
+
+CAMLprim value ocaml_glpk_unscale_problem(value blp)
+{
+  LPX *lp = Lpx_val(blp);
+  lpx_unscale_prob(lp);
+  return Val_unit;
+}
+
+/* TODO */
+/*
+CAMLprim value ocaml_glpk_check_kkt(value blp, value scaled, value vkkt)
+{
+
+}
+*/
+
+CAMLprim value ocaml_glpk_interior(value blp)
+{
+  CAMLparam1(blp);
+  LPX *lp = Lpx_val(blp);
+  int ret;
+
+  caml_enter_blocking_section();
+  ret = lpx_interior(lp);
+  caml_leave_blocking_section();
+
+  raise_on_error(ret);
+  CAMLreturn(Val_unit);
+}
+
+static int class_table[] = {LPX_LP, LPX_MIP};
+
+CAMLprim value ocaml_glpk_set_class(value blp, value class)
+{
+  LPX *lp = Lpx_val(blp);
+  lpx_set_class(lp, class_table[Int_val(class)]);
+  return Val_unit;
+}
+
+CAMLprim value ocaml_glpk_get_class(value blp)
+{
+  LPX *lp = Lpx_val(blp);
+  switch(lpx_get_class(lp))
+  {
+    case LPX_LP:
+      return Val_int(0);
+
+    case LPX_MIP:
+      return Val_int(1);
+
+    default:
+      assert(0);
+  }
+}
+
+static int kind_table[] = {LPX_CV, LPX_IV};
+
+CAMLprim value ocaml_glpk_set_col_kind(value blp, value n, value kind)
+{
+  LPX *lp = Lpx_val(blp);
+  lpx_set_col_kind(lp, Int_val(n) + 1, kind_table[Int_val(kind)]);
+  return Val_unit;
+}
+
+CAMLprim value ocaml_glpk_integer(value blp)
+{
+  CAMLparam1(blp);
+  LPX *lp = Lpx_val(blp);
+  int ret;
+
+  caml_enter_blocking_section();
+  ret = lpx_integer(lp);
+  caml_leave_blocking_section();
+
+  raise_on_error(ret);
+  CAMLreturn(Val_unit);
+}
+
+CAMLprim value ocaml_glpk_intopt(value blp)
+{
+  CAMLparam1(blp);
+  LPX *lp = Lpx_val(blp);
+  int ret;
+
+  caml_enter_blocking_section();
+  ret = lpx_intopt(lp);
+  caml_leave_blocking_section();
+
+  raise_on_error(ret);
+  CAMLreturn(Val_unit);
+}
+
+CAMLprim value ocaml_glpk_warm_up(value blp)
+{
+  LPX *lp = Lpx_val(blp);
+  raise_on_error(lpx_warm_up(lp));
+  return Val_unit;
+}
+
+#define BIND_INT_PARAM(name, param) \
+CAMLprim value ocaml_glpk_get_##name(value blp) \
+{ \
+  LPX *lp = Lpx_val(blp); \
+  return Val_int(lpx_get_int_parm(lp, param)); \
+} \
+CAMLprim value ocaml_glpk_set_##name(value blp, value n) \
+{ \
+  LPX *lp = Lpx_val(blp); \
+  lpx_set_int_parm(lp, param, Int_val(n)); \
+  return Val_unit; \
+}
+
+#define BIND_REAL_PARAM(name, param) \
+CAMLprim value ocaml_glpk_get_##name(value blp) \
+{ \
+  LPX *lp = Lpx_val(blp); \
+  double ans = lpx_get_real_parm(lp, param); \
+  return caml_copy_double(ans); \
+} \
+CAMLprim value ocaml_glpk_set_##name(value blp, value n) \
+{ \
+  LPX *lp = Lpx_val(blp); \
+  lpx_set_real_parm(lp, param, Double_val(n)); \
+  return Val_unit; \
+}
+
+BIND_INT_PARAM(message_level, LPX_K_MSGLEV);
+BIND_INT_PARAM(scaling, LPX_K_SCALE);
+BIND_INT_PARAM(use_dual_simplex, LPX_K_DUAL);
+BIND_INT_PARAM(pricing, LPX_K_PRICE);
+BIND_REAL_PARAM(relaxation, LPX_K_RELAX);
+/*
+BIND_REAL_PARAM(relative_tolerance, LPX_K_TOLBND);
+BIND_REAL_PARAM(absolute_tolerance, LPX_K_TOLDJ);
+*/
+BIND_INT_PARAM(solution_rounding, LPX_K_ROUND);
+BIND_INT_PARAM(iteration_limit, LPX_K_ITLIM);
+BIND_INT_PARAM(iteration_count, LPX_K_ITCNT);
+BIND_REAL_PARAM(time_limit, LPX_K_TMLIM);
+BIND_INT_PARAM(branching_heuristic, LPX_K_BRANCH);
+BIND_INT_PARAM(backtracking_heuristic, LPX_K_BTRACK);
+BIND_INT_PARAM(use_presolver, LPX_K_PRESOL);
+
+CAMLprim value ocaml_glpk_read_cplex(value fname)
+{
+  LPX *lp = lpx_read_cpxlp(String_val(fname));
+  if (!lp)
+    caml_failwith("Error while reading data in CPLEX LP format.");
+  return new_blp(lp);
+}
+
+CAMLprim value ocaml_glpk_write_cplex(value blp, value fname)
+{
+  if (lpx_write_cpxlp(Lpx_val(blp), String_val(fname)))
+    caml_failwith("Error while writing data in CPLEX LP format.");
+  return Val_unit;
+}