# HG changeset patch
# User Tommaso Balercia <tommaso.balercia@googlemail.com>
# Date 1301826798 -7200
# Node ID f31c421e9e7305e46d9d0a5892d277784e2ba397
# Parent  6ef23b4a34026cf72145ae20a7b7474c91c3f013
Added support for all the new features provided by GLPK 4.36+

diff -r 6ef23b4a3402 -r f31c421e9e73 scripts/optimization/glpk.m
--- a/scripts/optimization/glpk.m	Sat Apr 02 13:37:52 2011 -0700
+++ b/scripts/optimization/glpk.m	Sun Apr 03 12:33:18 2011 +0200
@@ -132,6 +132,9 @@
 ##
 ## @item "I"
 ## An integer variable.
+##
+## @item "B"
+## A binary variable.
 ## @end table
 ##
 ## @item sense
@@ -167,16 +170,22 @@
 ## Scaling option:
 ## @table @asis
 ## @item 0
-## No scaling.
+## Skip scaling.
 ##
 ## @item 1
-## Equilibration scaling.
+## Perform geometric mean scaling.
 ##
 ## @item 2
-## Geometric mean scaling, then equilibration scaling.
+## Perform equilibration scaling.
+##
+## @item 3
+## Perform automatic scaling
+##
+## @item 4
+## Round scale factors to nearest power of two
 ## @end table
 ##
-## @item dual    (@w{@code{LPX_K_DUAL}}, default: 0)
+## @item dual (@w{@code{LPX_K_DUAL}}, default: 0)
 ## Dual simplex option:
 ## @table @asis
 ## @item 0
@@ -184,6 +193,9 @@
 ##
 ## @item 1
 ## If initial basic solution is dual feasible, use the dual simplex.
+## 
+## @item 2
+## Use two phase dual simplex, or primal simplex if dual fails
 ## @end table
 ##
 ## @item price   (@w{@code{LPX_K_PRICE}}, default: 1)
@@ -196,7 +208,16 @@
 ## Steepest edge pricing.
 ## @end table
 ##
-## @item round   (@w{@code{LPX_K_ROUND}}, default: 0)
+## @item r_test   (default: 0)
+## @table @asis
+## @item 0
+## Standard (textbook).
+## 
+## @item 1
+## Harris's two-pass ratio test.
+## @end table
+##
+## @item round (@w{@code{LPX_K_ROUND}}, default: 0)
 ## Solution rounding option:
 ## @table @asis
 ## @item 0
@@ -206,7 +227,7 @@
 ## Replace tiny primal and dual values by exact zero.
 ## @end table
 ##
-## @item itlim   (@w{@code{LPX_K_ITLIM}}, default: -1)
+## @item itlim (@w{@code{LPX_K_ITLIM}}, default: -1)
 ## Simplex iterations limit.  If this value is positive, it is decreased by
 ## one each time when one simplex iteration has been performed, and
 ## reaching zero value signals the solver to stop the search.  Negative
@@ -227,6 +248,9 @@
 ## Branch on the last variable.
 ##
 ## @item 2
+## Branch on the most fractional variable.
+##
+## @item 3
 ## Branch using a heuristic by Driebeck and Tomlin.
 ## @end table
 ##
@@ -240,6 +264,9 @@
 ## Breadth first search.
 ##
 ## @item 2
+## Best local bound.
+##
+## @item 3
 ## Backtrack using the best projection heuristic.
 ## @end table
 ##
@@ -256,12 +283,96 @@
 ##
 ## @item 2
 ## Interior point method.
+##
+## @item 3
+## Simplex method with exact arithmetic
+## @end table
+##
+## @item pprocess (default: 2)
+## Pre-processing technique option (for MIP only).
+## @table @asis
+## @item 0
+## Disable preprocessing.
+##
+## @item 1
+## Perform preprocessing for root only.
+##
+## @item 2
+## Perform preprocessing for all levels.
+## @end table
+##
+## @item usecuts (default: 1)
+## Adds cutting planes in order to improve the LP relaxation before
+## applying the branch and bound method (for MIP only)
+## @table @asis
+## @item 0
+## Do not use cuts.
+##
+## @item 1
+## Use Gomoy's mixed integer cuts.
+##
+## @item 2
+## Use mixed integer rounding cuts.
+##
+## @item 3
+## Use mixed cover cuts.
+##
+## @item 4
+## Use clique cuts.
+##
+## @item 5
+## Use all cuts.
+## @end table
+##
+## @item binarize (default: 0)
+## Binarization option (for MIP only and used only if the presolver is
+## enabled).
+## @table @asis
+## @item 0
+## Do not use binarization.
+##
+## @item 1
+## Replace general integer variables with binary ones.
 ## @end table
 ##
 ## @item save (default: 0)
-## If this parameter is nonzero, save a copy of the problem in
-## CPLEX LP format to the file @file{"outpb.lp"}.  There is currently no
-## way to change the name of the output file.
+## If this parameter is nonzero, a copy of the problem is saved to file. 
+## Desired filename and format can be specified using the string parameters
+## described below. If these are not defined, the original problem is 
+## saved in CPLEX LP format to the file @file{outpb.lp}.  
+##
+## @item mpsinfo (default: 1)
+## If this is set, the interface writes to file several comment cards,
+## which contains some information about the problem. Otherwise the routine
+## writes no comment cards.
+##
+## @item mpsobj (default: 2)
+## This parameter tells the routing how to output the objective function
+## row. 
+## @table @asis
+## @item 0
+## Never output objective function row.
+##
+## @item 1
+## Always output objective function row.
+##
+## @item 2
+## Output objective function row if the problem has no free rows.
+## @end table
+##
+## @item mpsorig (default: 0)
+## If this is set, the routine uses the original symbolic names of rows
+## and columns. Otherwise the routine generates plain names using ordinal
+## numbers of rows and columns.
+##
+## @item mpswide (default: 1)
+## If this is set, the routine uses all data fields. Otherwise the routine
+## keeps fields 5 and 6 empty.
+##
+## @item mpsfree (default: 0)
+## If this is set, the routine omits column and vector names every time
+## when possible (free style). Otherwise the routine never omits these
+## names (pedantic style).
 ## @end table
 ##
 ## Real parameters:
@@ -323,6 +434,34 @@
 ## is not better than in the best known integer feasible solution.  It is
 ## not recommended that you change this parameter unless you have a
 ## detailed understanding of its purpose.
+##
+## @item mipgap (@w{@code{LPX_K_MIPGAPJ}}, default: 0.0)
+## The relative MIP gap tolerance. Inf the relative mip gap for currently
+## known best integer feasible solution falls below this tolerance, the
+## solver terminates the search. This allows obtaining suboptimal integer
+## feasible solutions if solving the problem to optimality takes too long.
+## @end table
+##
+## String Parameters:
+## @table @code
+## @item savefilename (default: @file{outpb})
+## The parameter allows saving the original problem in a file of desired
+## filename.
+## @item savefiletype (default: @file{cplex})
+## The parameter allows to specify the format type used to save the problem.
+## @table @asis
+## @item @file{fixedmps}
+## Fixed MPS format.
+##
+## @item @file{freemps}
+## Free MPS format.
+##
+## @item @file{cplex}
+## CPLEX LP format.
+##
+## @item @file{plain}
+## Plain text.
+## @end table
 ## @end table
 ## @end table
 ##
@@ -338,66 +477,106 @@
 ## @item status
 ## Status of the optimization.
 ##
-## Simplex Method:
 ## @table @asis
-## @item 180 (@w{@code{LPX_OPT}})
+## @item 1
+## Solution is undefined.
+## @item 2
+## Solution is feasible.
+## @item 3
+## Solution is infeasible.
+## @item 4
+## Problem has no feasible solution.
+## @item 5
 ## Solution is optimal.
-## @item 181 (@w{@code{LPX_FEAS}})
-## Solution is feasible.
-## @item 182 (@w{@code{LPX_INFEAS}})
-## Solution is infeasible.
-## @item 183 (@w{@code{LPX_NOFEAS}})
-## Problem has no feasible solution.
-## @item 184 (@w{@code{LPX_UNBND}})
-## Problem has no unbounded solution.
-## @item 185 (@w{@code{LPX_UNDEF}})
-## Solution status is undefined.
-## @end table
-## Interior Point Method:
-## @table @asis
-## @item 150 (@w{@code{LPX_T_UNDEF}})
-## The interior point method is undefined.
-## @item 151 (@w{@code{LPX_T_OPT}})
-## The interior point method is optimal.
-## @end table
-## Mixed Integer Method:
-## @table @asis
-## @item 170 (@w{@code{LPX_I_UNDEF}})
-## The status is undefined.
-## @item 171 (@w{@code{LPX_I_OPT}})
-## The solution is integer optimal.
-## @item 172 (@w{@code{LPX_I_FEAS}})
-## Solution integer feasible but its optimality has not been proven
-## @item 173 (@w{@code{LPX_I_NOFEAS}})
-## No integer feasible solution.
+## @item 6
+## Solution is unbounded.
 ## @end table
 ## @noindent
 ## If an error occurs, @var{status} will contain one of the following
-## codes:
+## codes.
 ##
+## Simplex method routines:
 ## @table @asis
-## @item 204 (@w{@code{LPX_E_FAULT}})
-## Unable to start the search.
-## @item 205 (@w{@code{LPX_E_OBJLL}})
-## Objective function lower limit reached.
-## @item 206 (@w{@code{LPX_E_OBJUL}})
-## Objective function upper limit reached.
-## @item 207 (@w{@code{LPX_E_ITLIM}})
-## Iterations limit exhausted.
-## @item 208 (@w{@code{LPX_E_TMLIM}})
-## Time limit exhausted.
-## @item 209 (@w{@code{LPX_E_NOFEAS}})
-## No feasible solution.
-## @item 210 (@w{@code{LPX_E_INSTAB}})
-## Numerical instability.
-## @item 211 (@w{@code{LPX_E_SING}})
-## Problems with basis matrix.
-## @item 212 (@w{@code{LPX_E_NOCONV}})
-## No convergence (interior).
-## @item 213 (@w{@code{LPX_E_NOPFS}})
-## No primal feasible solution (LP presolver).
-## @item 214 (@w{@code{LPX_E_NODFS}})
-## No dual feasible solution (LP presolver).
+## @item 101 (@w{@code{GLP_EBADB}})
+## Unable to start the search, because the initial basis specifed in the 
+## problem object is invalid|the number of basic (auxiliary and structural)
+## variables is not the same as the number of rows in the problem object.
+## @item 102 (@w{@code{GLP_ESING}})
+## Unable to start the search, because the basis matrix corresponding
+## to the initial basis is singular within the working precision.
+## @item 103 (@w{@code{GLP_ECOND}})
+## Unable to start the search, because the basis matrix corresponding
+## to the initial basis is ill-conditioned, i.e. its condition number is 
+## too large.
+## @item 104 (@w{@code{GLP_EBOUND}})
+## Unable to start the search, because some double-bounded (auxiliary or 
+## structural) variables have incorrect bounds.
+## @item 105 (@w{@code{GLP_EFAIL}})
+## The search was prematurely terminated due to the solver failure.
+## @item 106 (@w{@code{GLP_EOBJLL}})
+## The search was prematurely terminated, because the objective function 
+## being maximized has reached its lower limit and continues decreasing
+## (the dual simplex only).
+## @item 107 (@w{@code{GLP_EOBJUL}})
+## The search was prematurely terminated, because the objective function 
+## being minimized has reached its upper limit and continues increasing 
+## (the dual simplex only).
+## @item 108 (@w{@code{GLP_EITLIM}})
+## The search was prematurely terminated, because the simplex iteration 
+## limit has been exceeded.
+## @item 109 (@w{@code{GLP_ETMLIM}})
+## The search was prematurely terminated, because the time limit has been 
+## exceeded.
+## @item 110 (@w{@code{GLP_ENOPFS}})
+## The LP problem instance has no primal feasible solution (only if the LP
+## presolver is used).
+## @item 111 (@w{@code{GLP_ENODFS}})
+## The LP problem instance has no dual feasible solution (only if the LP 
+## presolver is used).
+## @end table
+##
+## Mixed integer programming routines:
+## @table @asis
+## @item 204 (@w{@code{GLP_EBOUND}})
+## Unable to start the search, because some double-bounded
+## variables have incorrect bounds or some integer variables
+## have non-integer (fractional) bounds.
+## @item 205 (@w{@code{GLP_EFAIL}})
+## The search was prematurely terminated due to the solver failure.
+## @item 209 (@w{@code{GLP_ETMLIM}})
+## The search was prematurely terminated, because the time limit has been 
+## exceeded.
+## @item 210 (@w{@code{GLP_ENOPFS}})
+## Unable to start the search, because LP relaxation of the
+## MIP problem instance has no primal feasible solution.
+## (This code may appear only if the presolver is enabled.)
+## @item 211 (@w{@code{GLP_ENODFS}})
+## Unable to start the search, because LP relaxation of the
+## MIP problem instance has no dual feasible solution. In
+## other word, this code means that if the LP relaxation has
+## at least one primal feasible solution, its optimal solution is
+## unbounded, so if the MIP problem has at least one integer
+## feasible solution, its (integer) optimal solution is also unbounded.
+## (This code may appear only if the presolver is enabled.)
+## @item 212 (@w{@code{GLP_EROOT}})
+## Unable to start the search, because optimal basis for initial
+## LP relaxation is not provided. (This code may appear only
+## if the presolver is disabled.)
+## @item 214 (@w{@code{GLP_EMIPGAP}})
+## The search was prematurely terminated, because the relative MIP gap 
+## tolerance has been reached.
+## @end table
+##
+## Interior point method routines:
+## @table @asis
+## @item 305 (@w{@code{GLP_EFAIL}})
+## The problem has no rows/columns.
+## @item 308 (@w{@code{GLP_EITLIM}})
+## Iteration limit exceeded.
+## @item 316 (@w{@code{GLP_ENOCVG}})
+## Very slow convergence or divergence.
+## @item 317 (@w{@code{GLP_EINSTAB}})
+## Numerical instability on solving Newtonian system.
 ## @end table
 ##
 ## @item extra
@@ -454,7 +633,7 @@
   endif
 
   if (all (size (c) > 1) || iscomplex (c) || ischar (c))
-    error ("glpk:C must be a real vector");
+    error ("glpk: c must be a real vector");
     return;
   endif
   nx = length (c);
@@ -536,8 +715,8 @@
             || length (vartype) != nx)
       error ("glpk: VARTYPE must be a char valued vector of length %d", nx);
       return;
-    elseif (! all (vartype == "C" | vartype == "I"))
-      error ("glpk: VARTYPE must contain only C or I");
+    elseif (! all (vartype == "C" | vartype == "I" | vartype == "B"))
+      error ("glpk: VARTYPE must contain only C, I or B");
       return;
     endif
   else
diff -r 6ef23b4a3402 -r f31c421e9e73 src/DLD-FUNCTIONS/__glpk__.cc
--- a/src/DLD-FUNCTIONS/__glpk__.cc	Sat Apr 02 13:37:52 2011 -0700
+++ b/src/DLD-FUNCTIONS/__glpk__.cc	Sun Apr 03 12:33:18 2011 +0200
@@ -1,5 +1,6 @@
 /*
 
+Copyright (C) 2011 Tommaso Balercia
 Copyright (C) 2005-2011 Nicolo' Giorgetti
 
 This file is part of Octave.
@@ -27,6 +28,8 @@
 #include <cfloat>
 #include <csetjmp>
 #include <ctime>
+#include <cstring>
+#include <string>
 
 #include "lo-ieee.h"
 
@@ -39,378 +42,593 @@
 
 #if defined (HAVE_GLPK)
 
-extern "C"
-{
-#if defined (HAVE_GLPK_GLPK_H)
-#include <glpk/glpk.h>
-#else
+extern "C" {
 #include <glpk.h>
-#endif
-
-#if 0
-#ifdef GLPK_PRE_4_14
-
-#ifndef _GLPLIB_H
-#include <glplib.h>
-#endif
-#ifndef lib_set_fault_hook
-#define lib_set_fault_hook lib_fault_hook
-#endif
-#ifndef lib_set_print_hook
-#define lib_set_print_hook lib_print_hook
-#endif
-
-#else
-
-void _glp_lib_print_hook (int (*func)(void *info, char *buf), void *info);
-void _glp_lib_fault_hook (int (*func)(void *info, char *buf), void *info);
-
-#endif
-#endif
 }
 
-#define NIntP 17
-#define NRealP 10
+#define NIntP 21
+#define NRealP 11
 
-int lpxIntParam[NIntP] = {
-  0,
-  1,
-  0,
-  1,
-  0,
-  -1,
-  0,
-  200,
-  1,
-  2,
-  0,
-  1,
-  0,
-  0,
-  2,
-  2,
-  1
+// control parameters not defined in glpk.h
+#define LPX_K_PREPROCESS    401   /* preprocessing  */
+#define LPX_K_RATIO_TEST    402   /* ratio test */
+
+
+// Integer Param Defaluts
+int glpIntParam[NIntP] = {
+    0,
+    1,
+    0,
+    1,
+    0,
+    INT_MAX,
+    INT_MAX,
+    200,
+    1,
+    2,
+    0,
+    1,
+    0,
+    0,
+    3,
+    2,
+    1,
+    0,
+    2,
+    0,
+    1
 };
 
+// Integer Param Names
 int IParam[NIntP] = {
-  LPX_K_MSGLEV,
-  LPX_K_SCALE,
-  LPX_K_DUAL,
-  LPX_K_PRICE,
-  LPX_K_ROUND,
-  LPX_K_ITLIM,
-  LPX_K_ITCNT,
-  LPX_K_OUTFRQ,
-  LPX_K_MPSINFO,
-  LPX_K_MPSOBJ,
-  LPX_K_MPSORIG,
-  LPX_K_MPSWIDE,
-  LPX_K_MPSFREE,
-  LPX_K_MPSSKIP,
-  LPX_K_BRANCH,
-  LPX_K_BTRACK,
-  LPX_K_PRESOL
+    LPX_K_MSGLEV,
+    LPX_K_SCALE,
+    LPX_K_DUAL,
+    LPX_K_PRICE,
+    LPX_K_ROUND,
+    LPX_K_ITLIM,
+    LPX_K_ITCNT,
+    LPX_K_OUTFRQ,
+    LPX_K_MPSINFO,
+    LPX_K_MPSOBJ,
+    LPX_K_MPSORIG,
+    LPX_K_MPSWIDE,
+    LPX_K_MPSFREE,
+    LPX_K_MPSSKIP,
+    LPX_K_BRANCH,
+    LPX_K_BTRACK,
+    LPX_K_PRESOL,
+    LPX_K_USECUTS,
+    LPX_K_PREPROCESS,
+    LPX_K_BINARIZE,
+    LPX_K_RATIO_TEST
 };
 
-
-double lpxRealParam[NRealP] = {
-  0.07,
-  1e-7,
-  1e-7,
-  1e-9,
-  -DBL_MAX,
-  DBL_MAX,
-  -1.0,
-  0.0,
-  1e-6,
-  1e-7
+//Real Param Values
+double glpRealParam[NRealP] = {
+    0.07,
+    1e-7,
+    1e-7,
+    1e-10,
+    -DBL_MAX,
+    DBL_MAX,
+    INT_MAX,
+    0.0,
+    1e-5,
+    1e-7,
+    0.0
 };
 
+//Real Param Names
 int RParam[NRealP] = {
-  LPX_K_RELAX,
-  LPX_K_TOLBND,
-  LPX_K_TOLDJ,
-  LPX_K_TOLPIV,
-  LPX_K_OBJLL,
-  LPX_K_OBJUL,
-  LPX_K_TMLIM,
-  LPX_K_OUTDLY,
-  LPX_K_TOLINT,
-  LPX_K_TOLOBJ
+    LPX_K_RELAX,
+    LPX_K_TOLBND,
+    LPX_K_TOLDJ,
+    LPX_K_TOLPIV,
+    LPX_K_OBJLL,
+    LPX_K_OBJUL,
+    LPX_K_TMLIM,
+    LPX_K_OUTDLY,
+    LPX_K_TOLINT,
+    LPX_K_TOLOBJ,
+    LPX_K_MIPGAP
 };
 
-static jmp_buf mark;  //-- Address for long jump to jump to
+static jmp_buf mark;  /*-- Address for long jump */
 
-#if 0
-int
-glpk_fault_hook (void * /* info */, char *msg)
+static int glpk_print_hook (void *info, const char *msg)
 {
-  error ("CRITICAL ERROR in GLPK: %s", msg);
-  longjmp (mark, -1);
+    int msglev = *reinterpret_cast<int*>(info);
+
+    if(msglev > 0)
+    {
+        message(0, "%s", msg);
+    }
+    return 1;
 }
 
-int
-glpk_print_hook (void * /* info */, char *msg)
+//
+int glpk (int sense, int n, int m, const double *c, int nz, const int *rn, const int *cn,
+          const double *a, const double *b, const char *ctype, const int *freeLB, const double *lb,
+          const int *freeUB, const double *ub, const int *vartype, int isMIP, int lpsolver,
+          int save_pb, const char * filename, const char * filetype, double *xmin, double *fmin, double *status,
+          double *lambda, double *redcosts, double *time, double *mem)
 {
-  message (0, "%s", msg);
-  return 1;
-}
-#endif
+    int typx = 0;
+    int method;
 
-int
-glpk (int sense, int n, int m, double *c, int nz, int *rn, int *cn,
-      double *a, double *b, char *ctype, int *freeLB, double *lb,
-      int *freeUB, double *ub, int *vartype, int isMIP, int lpsolver,
-      int save_pb, double *xmin, double *fmin, double *status,
-      double *lambda, double *redcosts, double *time, double *mem)
-{
-  int errnum;
-  int typx = 0;
-  int method;
+    clock_t t_start = clock();
 
-  clock_t t_start = clock();
+    //Redirect standard output
+    // Unfortunately at the moment parts of GLPK still print out messages
+    // without checking the specified options (e.g. the conflict graph part). 
+    // Passing msglev directly to the callback function is a workaround 
+    // for this problem, when the user desires no output.
+    int msglev=glpIntParam[0];
+    glp_term_hook (glpk_print_hook, reinterpret_cast<void*>(&msglev));
 
-#if 0
-#ifdef GLPK_PRE_4_14
-  lib_set_fault_hook (0, glpk_fault_hook);
-#else
-  _glp_lib_fault_hook (glpk_fault_hook, 0);
-#endif
+    //-- Create an empty LP/MILP object
+    glp_prob *lp = glp_create_prob ();
 
-  if (lpxIntParam[0] > 1)
-#ifdef GLPK_PRE_4_14
-    lib_set_print_hook (0, glpk_print_hook);
-#else
-    _glp_lib_print_hook (glpk_print_hook, 0);
-#endif
-#endif
+    //-- Set the sense of optimization
+    if (sense == 1)
+        glp_set_obj_dir (lp, GLP_MIN);
+    else
+        glp_set_obj_dir (lp, GLP_MAX);
 
-  LPX *lp = lpx_create_prob ();
-
-
-  //-- Set the sense of optimization
-  if (sense == 1)
-    lpx_set_obj_dir (lp, LPX_MIN);
-  else
-    lpx_set_obj_dir (lp, LPX_MAX);
-
-  //-- If the problem has integer structural variables switch to MIP
-  if (isMIP)
-    lpx_set_class (lp, LPX_MIP);
-
-  lpx_add_cols (lp, n);
-  for (int i = 0; i < n; i++)
+    //-- Define the number of unknowns and their domains.
+    glp_add_cols (lp, n);
+    for (int i = 0; i < n; i++)
     {
-      //-- Define type of the structural variables
-      if (! freeLB[i] && ! freeUB[i])
+        //-- Define type of the structural variables
+        if (! freeLB[i] && ! freeUB[i]) {
+            if ( lb[i] == ub[i] )
+                glp_set_col_bnds (lp, i+1, GLP_FX, lb[i], ub[i]);
+            else
+                glp_set_col_bnds (lp, i+1, GLP_DB, lb[i], ub[i]);
+        }
+        else
         {
-          if (lb[i] != ub[i])
-            lpx_set_col_bnds (lp, i+1, LPX_DB, lb[i], ub[i]);
-          else
-            lpx_set_col_bnds (lp, i+1, LPX_FX, lb[i], ub[i]);
-        }
-      else
-        {
-          if (! freeLB[i] && freeUB[i])
-            lpx_set_col_bnds (lp, i+1, LPX_LO, lb[i], ub[i]);
-          else
+            if (! freeLB[i] && freeUB[i])
+                glp_set_col_bnds (lp, i+1, GLP_LO, lb[i], ub[i]);
+            else
             {
-              if (freeLB[i] && ! freeUB[i])
-                lpx_set_col_bnds (lp, i+1, LPX_UP, lb[i], ub[i]);
-              else
-                lpx_set_col_bnds (lp, i+1, LPX_FR, lb[i], ub[i]);
+                if (freeLB[i] && ! freeUB[i])
+                    glp_set_col_bnds (lp, i+1, GLP_UP, lb[i], ub[i]);
+                else
+                    glp_set_col_bnds (lp, i+1, GLP_FR, lb[i], ub[i]);
             }
         }
 
-      // -- Set the objective coefficient of the corresponding
-      // -- structural variable. No constant term is assumed.
-      lpx_set_obj_coef(lp,i+1,c[i]);
+        // -- Set the objective coefficient of the corresponding
+        // -- structural variable. No constant term is assumed.
+        glp_set_obj_coef(lp,i+1,c[i]);
 
-      if (isMIP)
-        lpx_set_col_kind (lp, i+1, vartype[i]);
+        if (isMIP)
+            glp_set_col_kind (lp, i+1, vartype[i]);
     }
 
-  lpx_add_rows (lp, m);
+    glp_add_rows (lp, m);
 
-  for (int i = 0; i < m; i++)
+    for (int i = 0; i < m; i++)
     {
-      /* If the i-th row has no lower bound (types F,U), the
-         corrispondent parameter will be ignored.
-         If the i-th row has no upper bound (types F,L), the corrispondent
-         parameter will be ignored.
-         If the i-th row is of S type, the i-th LB is used, but
-         the i-th UB is ignored.
-      */
+        /*  If the i-th row has no lower bound (types F,U), the
+            correspondent parameter will be ignored.
+            If the i-th row has no upper bound (types F,L), the correspondent
+            parameter will be ignored.
+            If the i-th row is of S type, the i-th LB is used, but
+            the i-th UB is ignored.
+        */
 
-      switch (ctype[i])
+        switch (ctype[i])
         {
         case 'F':
-          typx = LPX_FR;
-          break;
-
+            typx = GLP_FR;
+            break;
+            // upper bound
         case 'U':
-          typx = LPX_UP;
-          break;
-
+            typx = GLP_UP;
+            break;
+            // lower bound
         case 'L':
-          typx = LPX_LO;
-          break;
-
+            typx = GLP_LO;
+            break;
+            // fixed constraint
         case 'S':
-          typx = LPX_FX;
-          break;
-
+            typx = GLP_FX;
+            break;
+            // double-bounded variable
         case 'D':
-          typx = LPX_DB;
-          break;
+            typx = GLP_DB;
+            break;
         }
 
-      lpx_set_row_bnds (lp, i+1, typx, b[i], b[i]);
+        if ( typx == GLP_DB && -b[i] < b[i]) {
+            glp_set_row_bnds (lp, i+1, typx, -b[i], b[i]);
+        }
+        else if(typx == GLP_DB && -b[i] == b[i]) {
+            glp_set_row_bnds (lp, i+1, GLP_FX, b[i], b[i]);
+        }
+        else {
+            glp_set_row_bnds (lp, i+1, typx, b[i], b[i]);
+        }
+
+    }
+    // Load constraint matrix A
+    glp_load_matrix (lp, nz, rn, cn, a);
+
+    // Save problem
+    if (save_pb) {
+        if (!strcmp(filetype,"cplex")) {
+            if (glp_write_lp (lp, NULL, filename) != 0) {
+                error("glpk: unable to write the problem to file");
+                longjmp (mark, -1);
+            }
+        } else {
+            if (!strcmp(filetype,"fixedmps")) {
+                if (glp_write_mps (lp, GLP_MPS_DECK, NULL, filename) != 0) {
+                    error("glpk: unable to write the problem to file");
+                    longjmp (mark, -1);
+                }
+            } else {
+                if (!strcmp(filetype,"freemps")) {
+                    if (glp_write_mps (lp, GLP_MPS_FILE, NULL, filename) != 0) {
+                        error("glpk: unable to write the problem to file");
+                        longjmp (mark, -1);
+                    }
+                } else { // plain text
+                    if (lpx_print_prob (lp, filename) != 0) {
+                        error("glpk: unable to write the problem to file");
+                        longjmp (mark, -1);
+                    }
+                }
+            }
+        }
+    }
+    //-- scale the problem data (if required)
+    if (! glpIntParam[16] || lpsolver != 1) {
+        switch ( glpIntParam[1] ) {
+        case ( 0 ):
+            glp_scale_prob( lp, GLP_SF_SKIP );
+            break;
+        case ( 1 ):
+            glp_scale_prob( lp, GLP_SF_GM );
+            break;
+        case ( 2 ):
+            glp_scale_prob( lp, GLP_SF_EQ );
+            break;
+        case ( 3 ):
+            glp_scale_prob( lp, GLP_SF_AUTO  );
+            break;
+        case ( 4 ):
+            glp_scale_prob( lp, GLP_SF_2N );
+            break;
+        default :
+            error("glpk: unrecognized scaling option");
+            longjmp (mark, -1);
+        }
+    }
+
+    //-- build advanced initial basis (if required)
+    if (lpsolver == 1 && ! glpIntParam[16])
+        glp_adv_basis (lp, 0);
+
+    glp_smcp sParam;
+    glp_init_smcp(&sParam);
+
+    //-- set control parameters for simplex/exact method
+    if (lpsolver == 1 || lpsolver == 3) {
+        //remap of control parameters for simplex method
+        sParam.msg_lev=glpIntParam[0];	// message level
+
+        // simplex method: primal/dual
+        switch ( glpIntParam[2] ) {
+        case 0:
+            sParam.meth=GLP_PRIMAL;
+            break;
+        case 1:
+            sParam.meth=GLP_DUAL;
+            break;
+        case 2:
+            sParam.meth=GLP_DUALP;
+            break;
+        default:
+            error("glpk: unrecognized primal/dual method");
+            longjmp (mark, -1);
+        }
+
+        // pricing technique
+        if (glpIntParam[3]==0) sParam.pricing=GLP_PT_STD;
+        else sParam.pricing=GLP_PT_PSE;
+
+        // ratio test
+        if (glpIntParam[20]==0) sParam.r_test = GLP_RT_STD;
+        else sParam.r_test=GLP_RT_HAR;
+
+        //tollerances
+        sParam.tol_bnd=glpRealParam[1];	// primal feasible tollerance
+        sParam.tol_dj=glpRealParam[2];	// dual feasible tollerance
+        sParam.tol_piv=glpRealParam[3];	// pivot tollerance
+        sParam.obj_ll=glpRealParam[4];	// lower limit
+        sParam.obj_ul=glpRealParam[5];	// upper limit
+
+        // iteration limit
+        if (glpIntParam[5]==-1) sParam.it_lim=INT_MAX;
+        else sParam.it_lim=glpIntParam[5];
+
+        // time limit
+        if (glpRealParam[6]==-1) sParam.tm_lim=INT_MAX;
+        else sParam.tm_lim=static_cast<int>(glpRealParam[6]);
+        sParam.out_frq=glpIntParam[7];	// output frequency
+        sParam.out_dly=static_cast<int>(glpRealParam[7]);	// output delay
+        // presolver
+        if (glpIntParam[16]) sParam.presolve=GLP_ON;
+        else sParam.presolve=GLP_OFF;
+    } else {
+        for(int i = 0; i < NIntP; i++) {
+            // skip assinging ratio test or
+            if ( i == 18 || i == 20) continue;
+            lpx_set_int_parm (lp, IParam[i], glpIntParam[i]);
+        }
+
+        for (int i = 0; i < NRealP; i++) {
+            lpx_set_real_parm (lp, RParam[i], glpRealParam[i]);
+        }
+    }
+
+    //set MIP params if MIP....
+    glp_iocp iParam;
+    glp_init_iocp(&iParam);
+
+    if ( isMIP ) {
+        method = 'I';
+
+        switch (glpIntParam[0]) { //message level
+        case 0:
+            iParam.msg_lev = GLP_MSG_OFF;
+            break;
+        case 1:
+            iParam.msg_lev = GLP_MSG_ERR;
+            break;
+        case 2:
+            iParam.msg_lev = GLP_MSG_ON;
+            break;
+        case 3:
+            iParam.msg_lev = GLP_MSG_ALL;
+            break;
+        default:
+            error("glpk: unrecognized msg_lev option");
+        }
+        switch (glpIntParam[14]) { //branching param
+        case 0:
+            iParam.br_tech = GLP_BR_FFV;
+            break;
+        case 1:
+            iParam.br_tech = GLP_BR_LFV;
+            break;
+        case 2:
+            iParam.br_tech = GLP_BR_MFV;
+            break;
+        case 3:
+            iParam.br_tech = GLP_BR_DTH;
+            break;
+        default:
+            error("glpk: unrecognized branch option");
+        }
+        switch (glpIntParam[15]) { //backtracking heuristic
+        case 0:
+            iParam.bt_tech = GLP_BT_DFS;
+            break;
+        case 1:
+            iParam.bt_tech = GLP_BT_BFS;
+            break;
+        case 2:
+            iParam.bt_tech = GLP_BT_BLB;
+            break;
+        case 3:
+            iParam.bt_tech = GLP_BT_BPH;
+            break;
+        default:
+            error("glpk: unrecognized backtrack option");
+        }
+
+        if (  glpRealParam[8] > 0.0 && glpRealParam[8] < 1.0 )
+            iParam.tol_int = glpRealParam[8];  // absolute tolorence
+        else
+            error("glpk: tolint must be between 0 and 1");
+
+        iParam.tol_obj = glpRealParam[9];  // relative tolarence
+        iParam.mip_gap = glpRealParam[10]; // realative gap tolerance
+
+        // set time limit for mip
+        if ( glpRealParam[6] < 0.0 || glpRealParam[6] > 1e6 )
+            iParam.tm_lim = INT_MAX;
+        else
+            iParam.tm_lim = static_cast<int>(1000.0 * glpRealParam[6] );
+
+        // Choose Cutsets for mip
+        // shut all cuts off, then start over....
+        iParam.gmi_cuts = GLP_OFF;
+        iParam.mir_cuts = GLP_OFF;
+        iParam.cov_cuts = GLP_OFF;
+        iParam.clq_cuts = GLP_OFF;
+
+        switch( glpIntParam[17] ) {
+        case 0:
+            break;
+        case 1:
+            iParam.gmi_cuts = GLP_ON;
+            break;
+        case 2:
+            iParam.mir_cuts = GLP_ON;
+            break;
+        case 3:
+            iParam.cov_cuts = GLP_ON;
+            break;
+        case 4:
+            iParam.clq_cuts = GLP_ON;
+            break;
+        case 5:
+            iParam.clq_cuts = GLP_ON;
+            iParam.gmi_cuts = GLP_ON;
+            iParam.mir_cuts = GLP_ON;
+            iParam.cov_cuts = GLP_ON;
+            break;
+        default:
+            error("glpk: unrecognized cut option");
+        }
+
+        switch( glpIntParam[18] ) { // pre-processing for mip
+        case 0:
+            iParam.pp_tech = GLP_PP_NONE;
+            break;
+        case 1:
+            iParam.pp_tech = GLP_PP_ROOT;
+            break;
+        case 2:
+            iParam.pp_tech = GLP_PP_ALL;
+            break;
+        default:
+            error("glpk: unrecognized pprocess option");
+        }
+
+        if (glpIntParam[16])  iParam.presolve=GLP_ON;
+        else                  iParam.presolve=GLP_OFF;
+
+        if (glpIntParam[19])  iParam.binarize = GLP_ON;
+        else                  iParam.binarize = GLP_OFF;
+
+    }
+    else {
+        /* Choose simplex method ('S')
+        or interior point method ('T')
+        or Exact method          ('E')
+        to solve the problem  */
+        switch (lpsolver) {
+        case 1:
+            method = 'S';
+            break;
+        case 2:
+            method = 'T';
+            break;
+        case 3:
+            method = 'E';
+            break;
+        default:
+            error("glpk:  unrecognized lpsolver option");
+            longjmp (mark, -1);
+        }
+    }
+
+    // now run the problem...
+    int errnum = 0;
+
+    switch (method) {
+    case 'I':
+        errnum = glp_intopt( lp, &iParam );
+        errnum += 200; //this is to avoid ambiguity in the return codes.
+        break;
+
+    case 'S':
+        errnum = glp_simplex(lp, &sParam);
+        errnum += 100; //this is to avoid ambiguity in the return codes.
+        break;
+
+    case 'T':
+        errnum = glp_interior(lp, NULL );
+        errnum += 300; //this is to avoid ambiguity in the return codes.
+        break;
+
+    case 'E':
+        errnum = glp_exact(lp, &sParam);
+        errnum += 100; //this is to avoid ambiguity in the return codes.
+        break;
+
+    default:  /*xassert (method != method); */
+        error("glpk: unrecognized problem");
+        longjmp (mark, -1);
 
     }
 
-  lpx_load_matrix (lp, nz, rn, cn, a);
+    if (errnum==100 || errnum==200 || errnum==300 || errnum==106 || errnum==107 || errnum==108 || errnum==109 || errnum==209 || errnum==214 || errnum==308) {
+        // Get status and object value
+        if (isMIP) {
+            *status = glp_mip_status (lp);
+            *fmin = glp_mip_obj_val (lp);
+        }
+        else {
 
-  if (save_pb)
-    {
-      static char tmp[] = "outpb.lp";
-      if (lpx_write_cpxlp (lp, tmp) != 0)
-        {
-          error ("__glpk__: unable to write problem");
-          longjmp (mark, -1);
-        }
-    }
-
-  //-- scale the problem data (if required)
-  //-- if (scale && (!presol || method == 1)) lpx_scale_prob(lp);
-  //-- LPX_K_SCALE=IParam[1]  LPX_K_PRESOL=IParam[16]
-  if (lpxIntParam[1] && (! lpxIntParam[16] || lpsolver != 1))
-    lpx_scale_prob (lp);
-
-  //-- build advanced initial basis (if required)
-  if (lpsolver == 1 && ! lpxIntParam[16])
-    lpx_adv_basis (lp);
-
-  for(int i = 0; i < NIntP; i++)
-    lpx_set_int_parm (lp, IParam[i], lpxIntParam[i]);
-
-  for (int i = 0; i < NRealP; i++)
-    lpx_set_real_parm (lp, RParam[i], lpxRealParam[i]);
-
-  if (lpsolver == 1)
-    method = 'S';
-  else
-    method = 'T';
-
-  switch (method)
-    {
-    case 'S':
-      {
-        if (isMIP)
-          {
-            method = 'I';
-            errnum = lpx_simplex (lp);
-            errnum = lpx_integer (lp);
-          }
-        else
-          errnum = lpx_simplex(lp);
-      }
-     break;
-
-    case 'T':
-      errnum = lpx_interior(lp);
-      break;
-
-    default:
-      break;
-#if 0
-#ifdef GLPK_PRE_4_14
-      insist (method != method);
-#else
-      static char tmp[] = "method != method";
-      glpk_fault_hook (0, tmp);
-#endif
-#endif
-    }
-
-  /*  errnum assumes the following results:
-      errnum = 0 <=> No errors
-      errnum = 1 <=> Iteration limit exceeded.
-      errnum = 2 <=> Numerical problems with basis matrix.
-  */
-  if (errnum == LPX_E_OK)
-    {
-      if (isMIP)
-        {
-          *status = lpx_mip_status (lp);
-          *fmin = lpx_mip_obj_val (lp);
-        }
-      else
-        {
-          if (lpsolver == 1)
-            {
-              *status = lpx_get_status (lp);
-              *fmin = lpx_get_obj_val (lp);
+            if (lpsolver == 1 || lpsolver == 3) {
+                *status = glp_get_status (lp);
+                *fmin = glp_get_obj_val (lp);
             }
-          else
-            {
-              *status = lpx_ipt_status (lp);
-              *fmin = lpx_ipt_obj_val (lp);
+            else {
+                *status = glp_ipt_status (lp);
+                *fmin = glp_ipt_obj_val (lp);
             }
         }
 
-      if (isMIP)
-        {
-          for (int i = 0; i < n; i++)
-            xmin[i] = lpx_mip_col_val (lp, i+1);
+        // Get optimal solution (if exists)
+        if (isMIP) {
+
+            for (int i = 0; i < n; i++)
+                xmin[i] = glp_mip_col_val (lp, i+1);
         }
-      else
-        {
-          /* Primal values */
-          for (int i = 0; i < n; i++)
-            {
-              if (lpsolver == 1)
-                xmin[i] = lpx_get_col_prim (lp, i+1);
-              else
-                xmin[i] = lpx_ipt_col_prim (lp, i+1);
+        else {
+
+            /* Primal values */
+            for (int i = 0; i < n; i++) {
+
+                if (lpsolver == 1 || lpsolver == 3)
+                    xmin[i] = glp_get_col_prim (lp, i+1);
+                else
+                    xmin[i] = glp_ipt_col_prim (lp, i+1);
             }
 
-          /* Dual values */
-          for (int i = 0; i < m; i++)
-            {
-              if (lpsolver == 1)
-                lambda[i] = lpx_get_row_dual (lp, i+1);
-              else
-                lambda[i] = lpx_ipt_row_dual (lp, i+1);
+            /* Dual values */
+            for (int i = 0; i < m; i++) {
+
+                if (lpsolver == 1 || lpsolver == 3)
+                    lambda[i] = glp_get_row_dual (lp, i+1);
+                else
+                    lambda[i] = glp_ipt_row_dual (lp, i+1);
             }
 
-          /* Reduced costs */
-          for (int i = 0; i < lpx_get_num_cols (lp); i++)
-            {
-              if (lpsolver == 1)
-                redcosts[i] = lpx_get_col_dual (lp, i+1);
-              else
-                redcosts[i] = lpx_ipt_col_dual (lp, i+1);
+            /* Reduced costs */
+            for (int i = 0; i < glp_get_num_cols (lp); i++) {
+
+                if (lpsolver == 1 || lpsolver == 3)
+                    redcosts[i] = glp_get_col_dual (lp, i+1);
+                else
+                    redcosts[i] = glp_ipt_col_dual (lp, i+1);
             }
+
         }
 
-      *time = (clock () - t_start) / CLOCKS_PER_SEC;
+        *time = (clock () - t_start) / CLOCKS_PER_SEC;
 
-#ifdef GLPK_PRE_4_14
-      *mem = (lib_env_ptr () -> mem_tpeak);
-#else
-      *mem = 0;
-#endif
+        glp_long tpeak;
+        glp_mem_usage(NULL, NULL, NULL, &tpeak);
+        *mem=static_cast<double>((4294967296.0 * tpeak.hi + tpeak.lo) / (1024));
 
-      lpx_delete_prob (lp);
-      return 0;
+        glp_delete_prob (lp);
+
+        return 0;
+    }
+    else {
+        // printf("errnum is %d\n", errnum);
     }
 
-   lpx_delete_prob (lp);
+    glp_delete_prob (lp);
 
-   *status = errnum;
+    /* this shouldn't be nessiary with glp_deleted_prob, but try it
+    if we have weird behavior again... */
+    glp_free_env();
 
-   return errnum;
+
+    *status = errnum;
+
+    return errnum;
 }
 
 #endif
@@ -424,7 +642,7 @@
         { \
           if (! tmp.is_empty ()) \
             { \
-              lpxRealParam[IDX] = tmp.scalar_value (); \
+              glpRealParam[IDX] = tmp.scalar_value (); \
  \
               if (error_state) \
                 { \
@@ -468,389 +686,485 @@
   while (0)
 
 DEFUN_DLD (__glpk__, args, ,
-  "-*- texinfo -*-\n\
+           "-*- texinfo -*-\n\
 @deftypefn {Loadable Function} {[@var{values}] =} __glpk__ (@var{args})\n\
 Undocumented internal function.\n\
 @end deftypefn")
 {
-  // The list of values to return.  See the declaration in oct-obj.h
-  octave_value_list retval;
+    // The list of values to return.  See the declaration in oct-obj.h
+    octave_value_list retval;
 
 #if defined (HAVE_GLPK)
 
-  int nrhs = args.length ();
+    int nrhs = args.length ();
 
-  if (nrhs != 9)
+    if (nrhs != 9)
     {
-      print_usage ();
-      return retval;
+        print_usage ();
+        return retval;
     }
 
-  //-- 1nd Input. A column array containing the objective function
-  //--            coefficients.
-  volatile int mrowsc = args(0).rows();
+    //-- 1nd Input. A column array containing the objective function
+    //--            coefficients.
+    volatile int mrowsc = args(0).rows();
 
-  Matrix C (args(0).matrix_value ());
+    Matrix C (args(0).matrix_value ());
 
-  if (error_state)
+    if (error_state)
     {
-      error ("__glpk__: invalid value of C");
-      return retval;
+        error ("glpk: invalid value of c");
+        return retval;
     }
 
-  double *c = C.fortran_vec ();
-  Array<int> rn;
-  Array<int> cn;
-  ColumnVector a;
-  volatile int mrowsA;
-  volatile int nz = 0;
+    double *c = C.fortran_vec ();
+    Array<int> rn;
+    Array<int> cn;
+    ColumnVector a;
+    volatile int mrowsA;
+    volatile int nz = 0;
 
-  //-- 2nd Input. A matrix containing the constraints coefficients.
-  // If matrix A is NOT a sparse matrix
-  if (args(1).is_sparse_type ())
+    //-- 2nd Input. A matrix containing the constraints coefficients.
+    // If matrix A is NOT a sparse matrix
+    if (args(1).is_sparse_type ())
     {
-      SparseMatrix A = args(1).sparse_matrix_value (); // get the sparse matrix
+        SparseMatrix A = args(1).sparse_matrix_value (); // get the sparse matrix
 
-      if (error_state)
+        if (error_state)
         {
-          error ("__glpk__: invalid value of A");
-          return retval;
+            error ("glpk: invalid value of A");
+            return retval;
         }
 
-      mrowsA = A.rows ();
-      octave_idx_type Anc = A.cols ();
-      octave_idx_type Anz = A.nnz ();
-      rn.resize (dim_vector (Anz+1, 1));
-      cn.resize (dim_vector (Anz+1, 1));
-      a.resize (Anz+1, 0.0);
+        mrowsA = A.rows ();
+        octave_idx_type Anc = A.cols ();
+        octave_idx_type Anz = A.nnz ();
+        rn.resize (dim_vector (Anz+1, 1));
+        cn.resize (dim_vector (Anz+1, 1));
+        a.resize (Anz+1, 0.0);
 
-      if (Anc != mrowsc)
+        if (Anc != mrowsc)
         {
-          error ("__glpk__: invalid value of A");
-          return retval;
+            error ("glpk: invalid value of A");
+            return retval;
         }
 
-      for (octave_idx_type j = 0; j < Anc; j++)
-        for (octave_idx_type i = A.cidx(j); i < A.cidx(j+1); i++)
-          {
-            nz++;
-            rn(nz) = A.ridx(i) + 1;
-            cn(nz) = j + 1;
-            a(nz) = A.data(i);
-          }
+        for (octave_idx_type j = 0; j < Anc; j++)
+            for (octave_idx_type i = A.cidx(j); i < A.cidx(j+1); i++)
+            {
+                nz++;
+                rn(nz) = A.ridx(i) + 1;
+                cn(nz) = j + 1;
+                a(nz) = A.data(i);
+            }
     }
-  else
+    else
     {
-      Matrix A (args(1).matrix_value ()); // get the matrix
+        Matrix A (args(1).matrix_value ()); // get the matrix
 
-      if (error_state)
+        if (error_state)
         {
-          error ("__glpk__: invalid value of A");
-          return retval;
+            error ("glpk: invalid value of A");
+            return retval;
         }
 
-      mrowsA = A.rows ();
-      rn.resize (dim_vector (mrowsA*mrowsc+1, 1));
-      cn.resize (dim_vector (mrowsA*mrowsc+1, 1));
-      a.resize (mrowsA*mrowsc+1, 0.0);
+        mrowsA = A.rows ();
+        rn.resize (dim_vector (mrowsA*mrowsc+1, 1));
+        cn.resize (dim_vector (mrowsA*mrowsc+1, 1));
+        a.resize (mrowsA*mrowsc+1, 0.0);
 
-      for (int i = 0; i < mrowsA; i++)
+        for (int i = 0; i < mrowsA; i++)
         {
-          for (int j = 0; j < mrowsc; j++)
+            for (int j = 0; j < mrowsc; j++)
             {
-              if (A(i,j) != 0)
+                if (A(i,j) != 0)
                 {
-                  nz++;
-                  rn(nz) = i + 1;
-                  cn(nz) = j + 1;
-                  a(nz) = A(i,j);
+                    nz++;
+                    rn(nz) = i + 1;
+                    cn(nz) = j + 1;
+                    a(nz) = A(i,j);
                 }
             }
         }
 
     }
 
-  //-- 3rd Input. A column array containing the right-hand side value
-  //               for each constraint in the constraint matrix.
-  Matrix B (args(2).matrix_value ());
+    //-- 3rd Input. A column array containing the right-hand side value
+    //               for each constraint in the constraint matrix.
+    Matrix B (args(2).matrix_value ());
 
-  if (error_state)
+    if (error_state)
     {
-      error ("__glpk__: invalid value of B");
-      return retval;
+        error ("glpk: invalid value of B");
+        return retval;
     }
 
-  double *b = B.fortran_vec ();
+    double *b = B.fortran_vec ();
 
-  //-- 4th Input. An array of length mrowsc containing the lower
-  //--            bound on each of the variables.
-  Matrix LB (args(3).matrix_value ());
+    //-- 4th Input. An array of length mrowsc containing the lower
+    //--            bound on each of the variables.
+    Matrix LB (args(3).matrix_value ());
 
-  if (error_state || LB.length () < mrowsc)
+    if (error_state || LB.length () < mrowsc)
     {
-      error ("__glpk__: invalid value of LB");
-      return retval;
+        error ("glpk: invalid value of LB");
+        return retval;
     }
 
-  double *lb = LB.fortran_vec ();
+    double *lb = LB.fortran_vec ();
 
-  //-- LB argument, default: Free
-  Array<int> freeLB (dim_vector (mrowsc, 1));
-  for (int i = 0; i < mrowsc; i++)
-     {
-       if (xisinf (lb[i]))
-         {
-           freeLB(i) = 1;
-           lb[i] = -octave_Inf;
-         }
-       else
-         freeLB(i) = 0;
-     }
-
-  //-- 5th Input. An array of at least length numcols containing the upper
-  //--            bound on each of the variables.
-  Matrix UB (args(4).matrix_value ());
-
-  if (error_state || UB.length () < mrowsc)
+    //-- LB argument, default: Free
+    Array<int> freeLB (dim_vector (mrowsc, 1));
+    for (int i = 0; i < mrowsc; i++)
     {
-      error ("__glpk__: invalid value of UB");
-      return retval;
+        if (xisinf (lb[i]))
+        {
+            freeLB(i) = 1;
+            lb[i] = -octave_Inf;
+        }
+        else
+            freeLB(i) = 0;
     }
 
-  double *ub = UB.fortran_vec ();
+    //-- 5th Input. An array of at least length numcols containing the upper
+    //--            bound on each of the variables.
+    Matrix UB (args(4).matrix_value ());
 
-  Array<int> freeUB (dim_vector (mrowsc, 1));
-  for (int i = 0; i < mrowsc; i++)
+    if (error_state || UB.length () < mrowsc)
     {
-      if (xisinf (ub[i]))
-        {
-          freeUB(i) = 1;
-          ub[i] = octave_Inf;
-        }
-      else
-        freeUB(i) = 0;
+        error ("glpk: invalid value of UB");
+        return retval;
     }
 
-  //-- 6th Input. A column array containing the sense of each constraint
-  //--            in the constraint matrix.
-  charMatrix CTYPE (args(5).char_matrix_value ());
+    double *ub = UB.fortran_vec ();
 
-  if (error_state)
+    Array<int> freeUB (dim_vector (mrowsc, 1));
+    for (int i = 0; i < mrowsc; i++)
     {
-      error ("__glpk__: invalid value of CTYPE");
-      return retval;
+        if (xisinf (ub[i]))
+        {
+            freeUB(i) = 1;
+            ub[i] = octave_Inf;
+        }
+        else
+            freeUB(i) = 0;
     }
 
-  char *ctype = CTYPE.fortran_vec ();
+    //-- 6th Input. A column array containing the sense of each constraint
+    //--            in the constraint matrix.
+    charMatrix CTYPE (args(5).char_matrix_value ());
 
-  //-- 7th Input. A column array containing the types of the variables.
-  charMatrix VTYPE (args(6).char_matrix_value ());
-
-  if (error_state)
+    if (error_state)
     {
-      error ("__glpk__: invalid value of VARTYPE");
-      return retval;
+        error ("glpk: invalid value of CTYPE");
+        return retval;
     }
 
-  Array<int> vartype (dim_vector (mrowsc, 1));
-  volatile int isMIP = 0;
-  for (int i = 0; i < mrowsc ; i++)
+    char *ctype = CTYPE.fortran_vec ();
+
+    //-- 7th Input. A column array containing the types of the variables.
+    charMatrix VTYPE (args(6).char_matrix_value ());
+
+    if (error_state)
     {
-      if (VTYPE(i,0) == 'I')
-        {
-          isMIP = 1;
-          vartype(i) = LPX_IV;
-        }
-      else
-        vartype(i) = LPX_CV;
+        error ("glpk: invalid value of VARTYPE");
+        return retval;
     }
 
-  //-- 8th Input. Sense of optimization.
-  volatile int sense;
-  double SENSE = args(7).scalar_value ();
-
-  if (error_state)
+    Array<int> vartype (dim_vector (mrowsc, 1));
+    volatile int isMIP = 0;
+    for (int i = 0; i < mrowsc ; i++)
     {
-      error ("__glpk__: invalid value of SENSE");
-      return retval;
+        switch(VTYPE(i,0)) {
+        case 'I':
+            vartype(i) = GLP_IV;
+            isMIP = 1;
+            break;
+        case 'B':
+            vartype(i) = GLP_BV;
+            isMIP = 1;
+            break;
+        default:
+            vartype(i) = GLP_CV;
+            isMIP = 0;
+        }
     }
 
-  if (SENSE >= 0)
-    sense = 1;
-  else
-    sense = -1;
+    //-- 8th Input. Sense of optimization.
+    volatile int sense;
+    double SENSE = args(7).scalar_value ();
 
-  //-- 9th Input. A structure containing the control parameters.
-  octave_scalar_map PARAM = args(8).scalar_map_value ();
-
-  if (error_state)
+    if (error_state)
     {
-      error ("__glpk__: invalid value of PARAM");
-      return retval;
+        error ("glpk: invalid value of SENSE");
+        return retval;
     }
 
-  //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-  //-- Integer parameters
-  //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+    if (SENSE >= 0)
+        sense = 1;
+    else
+        sense = -1;
 
-  //-- Level of messages output by the solver
-  OCTAVE_GLPK_GET_INT_PARAM ("msglev", lpxIntParam[0]);
-  if (lpxIntParam[0] < 0 || lpxIntParam[0] > 3)
+    //-- 9th Input. A structure containing the control parameters.
+    octave_scalar_map PARAM = args(8).scalar_map_value ();
+
+    if (error_state)
     {
-      error ("__glpk__: PARAM.msglev must be 0 (no output [default]) or 1 (error messages only) or 2 (normal output) or 3 (full output)");
-      return retval;
+        error ("glpk: invalid value of PARAM");
+        return retval;
     }
 
-  //-- scaling option
-  OCTAVE_GLPK_GET_INT_PARAM ("scale", lpxIntParam[1]);
-  if (lpxIntParam[1] < 0 || lpxIntParam[1] > 2)
+    //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+    //-- Integer parameters
+    //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+    //-- Level of messages output by the solver
+    OCTAVE_GLPK_GET_INT_PARAM ("msglev", glpIntParam[0]);
+    if (glpIntParam[0] < 0 || glpIntParam[0] > 3)
     {
-      error ("__glpk__: PARAM.scale must be 0 (no scaling) or 1 (equilibration scaling [default]) or 2 (geometric mean scaling)");
-      return retval;
+        error ("glpk: the parameter \"msglev\" must be 0 (no output [default]) or 1 (error messages only) or 2 (normal output) or 3 (full output)");
+        return retval;
     }
 
-  //-- Dual dimplex option
-  OCTAVE_GLPK_GET_INT_PARAM ("dual", lpxIntParam[2]);
-  if (lpxIntParam[2] < 0 || lpxIntParam[2] > 1)
+    //-- scaling option
+    OCTAVE_GLPK_GET_INT_PARAM ("scale", glpIntParam[1]);
+    if (glpIntParam[1] < 0 || glpIntParam[1] > 4)
     {
-      error ("__glpk__: PARAM.dual must be 0 (do NOT use dual simplex [default]) or 1 (use dual simplex)");
-      return retval;
+        error ("glpk: the parameter \"scale\" must be 0 (no scaling) or 1 (equilibration scaling [default]) or 2 (geometric mean scaling)");
+        return retval;
     }
 
-  //-- Pricing option
-  OCTAVE_GLPK_GET_INT_PARAM ("price", lpxIntParam[3]);
-  if (lpxIntParam[3] < 0 || lpxIntParam[3] > 1)
+    //-- Dual dimplex option
+    OCTAVE_GLPK_GET_INT_PARAM ("dual", glpIntParam[2]);
+    if (glpIntParam[2] < 0 || glpIntParam[2] > 2)
     {
-      error ("__glpk__: PARAM.price must be 0 (textbook pricing) or 1 (steepest edge pricing [default])");
-      return retval;
+        error ("glpk: the parameters \"dual\" must be 0 (do NOT use dual simplex [default]) or 1 (use dual simplex)");
+        return retval;
     }
 
-  //-- Solution rounding option
-  OCTAVE_GLPK_GET_INT_PARAM ("round", lpxIntParam[4]);
-  if (lpxIntParam[4] < 0 || lpxIntParam[4] > 1)
+    //-- Pricing option
+    OCTAVE_GLPK_GET_INT_PARAM ("price", glpIntParam[3]);
+    if (glpIntParam[3] < 0 || glpIntParam[3] > 1)
     {
-      error ("__glpk__: PARAM.round must be 0 (report all primal and dual values [default]) or 1 (replace tiny primal and dual values by exact zero)");
-      return retval;
+        error ("glpk: the parameter \"price\" must be 0 (textbook pricing) or 1 (steepest edge pricing [default])");
+        return retval;
     }
 
-  //-- Simplex iterations limit
-  OCTAVE_GLPK_GET_INT_PARAM ("itlim", lpxIntParam[5]);
-
-  //-- Simplex iterations count
-  OCTAVE_GLPK_GET_INT_PARAM ("itcnt", lpxIntParam[6]);
-
-  //-- Output frequency, in iterations
-  OCTAVE_GLPK_GET_INT_PARAM ("outfrq", lpxIntParam[7]);
-
-  //-- Branching heuristic option
-  OCTAVE_GLPK_GET_INT_PARAM ("branch", lpxIntParam[14]);
-  if (lpxIntParam[14] < 0 || lpxIntParam[14] > 2)
+    //-- Ratio test option
+    OCTAVE_GLPK_GET_INT_PARAM ("r_test", glpIntParam[20]);
+    if (glpIntParam[20] < 0 || glpIntParam[20] > 1)
     {
-      error ("__glpk__: PARAM.branch must be (MIP only) 0 (branch on first variable) or 1 (branch on last variable) or 2 (branch using a heuristic by Driebeck and Tomlin [default]");
-      return retval;
+        error ("glpk: the parameter \"r_test\" must be 0 (textbook [default]) or 1 (Harris's Two pass ratio test)");
+        return retval;
     }
 
-  //-- Backtracking heuristic option
-  OCTAVE_GLPK_GET_INT_PARAM ("btrack", lpxIntParam[15]);
-  if (lpxIntParam[15] < 0 || lpxIntParam[15] > 2)
+    //-- Solution rounding option
+    OCTAVE_GLPK_GET_INT_PARAM ("round", glpIntParam[4]);
+    if (glpIntParam[4] < 0 || glpIntParam[4] > 1)
     {
-      error ("__glpk__: PARAM.btrack must be (MIP only) 0 (depth first search) or 1 (breadth first search) or 2 (backtrack using the best projection heuristic [default]");
-      return retval;
+        error ("glpk: the parameter \"round\" must be 0 (report all primal and dual values [default]) or 1 (replace tiny primal and dual values by exact zero)");
+        return retval;
     }
 
-  //-- Presolver option
-  OCTAVE_GLPK_GET_INT_PARAM ("presol", lpxIntParam[16]);
-  if (lpxIntParam[16] < 0 || lpxIntParam[16] > 1)
+    //-- Simplex iterations limit
+    OCTAVE_GLPK_GET_INT_PARAM ("itlim", glpIntParam[5]);
+
+    //-- Simplex iterations count
+    OCTAVE_GLPK_GET_INT_PARAM ("itcnt", glpIntParam[6]);
+
+    //-- Output frequency, in iterations
+    OCTAVE_GLPK_GET_INT_PARAM ("outfrq", glpIntParam[7]);
+
+    //-- Branching heuristic option
+    OCTAVE_GLPK_GET_INT_PARAM ("branch", glpIntParam[14]);
+    if (glpIntParam[14] < 0 || glpIntParam[14] > 3)
     {
-      error ("__glpk__: PARAM.presol must be 0 (do NOT use LP presolver) or 1 (use LP presolver [default])");
-      return retval;
+        error ("glpk: the parameter \"branch\" must be 0 (branch on first variable) or 1 (branch on last variable) or 2 (most fractional variable) or 3 (branch using a heuristic by Driebeck and Tomlin [default]");
+        return retval;
     }
 
-  //-- LPsolver option
-  volatile int lpsolver = 1;
-  OCTAVE_GLPK_GET_INT_PARAM ("lpsolver", lpsolver);
-  if (lpsolver < 1 || lpsolver > 2)
+    //-- Backtracking heuristic option
+    OCTAVE_GLPK_GET_INT_PARAM ("btrack", glpIntParam[15]);
+    if (glpIntParam[15] < 0 || glpIntParam[15] > 3)
     {
-      error ("__glpk__: PARAM.lpsolver must be 1 (simplex method) or 2 (interior point method)");
-      return retval;
+        error ("glpk: the parameter \"btrack\" must be 0 (depth first search) or 1 (breadth first search) or 2 ( best local bound ) or 3 (backtrack using the best projection heuristic [default]");
+        return retval;
     }
 
-  //-- Save option
-  volatile int save_pb = 0;
-  OCTAVE_GLPK_GET_INT_PARAM ("save", save_pb);
-  save_pb = save_pb != 0;
-
-  //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-  //-- Real parameters
-  //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-  //-- Ratio test option
-  OCTAVE_GLPK_GET_REAL_PARAM ("relax", 0);
-
-  //-- Relative tolerance used to check if the current basic solution
-  //-- is primal feasible
-  OCTAVE_GLPK_GET_REAL_PARAM ("tolbnd", 1);
-
-  //-- Absolute tolerance used to check if the current basic solution
-  //-- is dual feasible
-  OCTAVE_GLPK_GET_REAL_PARAM ("toldj", 2);
-
-  //-- Relative tolerance used to choose eligible pivotal elements of
-  //--  the simplex table in the ratio test
-  OCTAVE_GLPK_GET_REAL_PARAM ("tolpiv", 3);
-
-  OCTAVE_GLPK_GET_REAL_PARAM ("objll", 4);
-
-  OCTAVE_GLPK_GET_REAL_PARAM ("objul", 5);
-
-  OCTAVE_GLPK_GET_REAL_PARAM ("tmlim", 6);
-
-  OCTAVE_GLPK_GET_REAL_PARAM ("outdly", 7);
-
-  OCTAVE_GLPK_GET_REAL_PARAM ("tolint", 8);
-
-  OCTAVE_GLPK_GET_REAL_PARAM ("tolobj", 9);
-
-  //-- Assign pointers to the output parameters
-  ColumnVector xmin (mrowsc, octave_NA);
-  double fmin = octave_NA;
-  double status;
-  ColumnVector lambda (mrowsA, octave_NA);
-  ColumnVector redcosts (mrowsc, octave_NA);
-  double time;
-  double mem;
-
-  int jmpret = setjmp (mark);
-
-  if (jmpret == 0)
-    glpk (sense, mrowsc, mrowsA, c, nz, rn.fortran_vec (),
-          cn.fortran_vec (), a.fortran_vec (), b, ctype,
-          freeLB.fortran_vec (), lb, freeUB.fortran_vec (), ub,
-          vartype.fortran_vec (), isMIP, lpsolver, save_pb,
-          xmin.fortran_vec (), &fmin, &status, lambda.fortran_vec (),
-          redcosts.fortran_vec (), &time, &mem);
-
-  octave_scalar_map extra;
-
-  if (! isMIP)
+    //-- Presolver option
+    OCTAVE_GLPK_GET_INT_PARAM ("presol", glpIntParam[16]);
+    if (glpIntParam[16] < 0 || glpIntParam[16] > 1)
     {
-      extra.assign ("lambda", lambda);
-      extra.assign ("redcosts", redcosts);
+        error ("glpk: the parameter \"presol\" must be 0 (do NOT use LP presolver [default]) or 1 (use LP presolver [default])");
+        return retval;
     }
 
-  extra.assign ("time", time);
-  extra.assign ("mem", mem);
+    //-- Generating cuts
+    OCTAVE_GLPK_GET_INT_PARAM ("usecuts", glpIntParam[17]);
+    if (glpIntParam[17] < 0 || glpIntParam[17] > 5)
+    {
+        error ("glpk: the parameter \"usecuts\" must be 0 (do NOT generate cuts) or 1 (generate Gomory's cuts [default]) or 2 (generate mixed integer rounding cuts) or 3 (generate cover cuts) or 4 (generate clique cuts) or 5 (generate all cuts)");
+        return retval;
+    }
 
-  retval(3) = extra;
-  retval(2) = status;
-  retval(1) = fmin;
-  retval(0) = xmin;
+    //-- PreProcessing
+    OCTAVE_GLPK_GET_INT_PARAM ("pprocess", glpIntParam[18]);
+    if (glpIntParam[18] < 0 || glpIntParam[18] > 2)
+    {
+        error ("glpk: the parameter \"pprocess\" must be 0 (disable preprocessing) or 1 (preprocessing for root only) or 2 (preprocessing for all levels [default])");
+        return retval;
+    }
+
+    //-- Binarize
+    OCTAVE_GLPK_GET_INT_PARAM ("binarize", glpIntParam[19]);
+    if (glpIntParam[19] < 0 || glpIntParam[19] > 1)
+    {
+        error ("glpk: the parameter \"binarize\" must be 0 (do use binarization [default]) or 1 (replace general integer variable by binary ones)");
+        return retval;
+    }
+
+    //-- LPsolver option
+    volatile int lpsolver = 1;
+    OCTAVE_GLPK_GET_INT_PARAM ("lpsolver", lpsolver);
+    if (lpsolver < 1 || lpsolver > 3)
+    {
+        error ("glpk: the parameter \"lpsolver\" must be 1 (simplex method [default]) or 2 (interior point method) or 3 (LP in exact arithmetic)");
+        return retval;
+    }
+
+    //-- Save option
+    volatile int save_pb = 0;
+    std::string filename("outpb");
+    std::string filetype("cplex");
+    OCTAVE_GLPK_GET_INT_PARAM ("save", save_pb);
+    save_pb = save_pb != 0;
+    if (save_pb) {
+        octave_value name = PARAM.getfield("savefilename");
+
+        if (name.is_defined() )
+        {
+            if (!name.is_empty())
+            {
+                filename = name.string_value();
+            }
+        }
+
+        octave_value type = PARAM.getfield("savefiletype");
+
+        if(type.is_defined())
+        {
+            if (!name.is_empty())
+            {
+                filetype = type.string_value();
+            }
+        }
+
+        if(filetype.compare("cplex") == 0)
+        {
+            filename += ".lp";
+        }
+        else
+        {
+            if(filetype.compare("fixedmps") == 0 || filetype.compare("freemps") == 0)
+            {
+                filename += std::string(".mps");
+            }
+            else
+            {
+                filename += std::string(".txt");
+            }
+        }
+    }
+
+    // MPS parameters
+    //-- mpsinfo
+    OCTAVE_GLPK_GET_INT_PARAM ("mpsinfo", glpIntParam[8]);
+    //-- mpsobj
+    OCTAVE_GLPK_GET_INT_PARAM ("mpsobj", glpIntParam[9]);
+    if (glpIntParam[9] < 0 || glpIntParam[9] > 2)
+    {
+        error("glpk: the parameter \"mpsobj\" must be 0 (never output objective function row) or 1 (always output objective function row ) or 2 (output objective function row if the problem has no free rows [default])");
+        return retval;
+    }
+    //-- mpsorig
+    OCTAVE_GLPK_GET_INT_PARAM ("mpsorig", glpIntParam[10]);
+    //-- mpswide
+    OCTAVE_GLPK_GET_INT_PARAM ("mpswide", glpIntParam[11]);
+    //-- mpsfree
+    OCTAVE_GLPK_GET_INT_PARAM ("mpsfree", glpIntParam[12]);
+
+    //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+    //-- Real parameters
+    //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+    //-- Ratio test option
+    OCTAVE_GLPK_GET_REAL_PARAM ("relax", 0);
+
+    //-- Relative tolerance used to check if the current basic solution
+    //-- is primal feasible
+    OCTAVE_GLPK_GET_REAL_PARAM ("tolbnd", 1);
+
+    //-- Absolute tolerance used to check if the current basic solution
+    //-- is dual feasible
+    OCTAVE_GLPK_GET_REAL_PARAM ("toldj", 2);
+
+    //-- Relative tolerance used to choose eligible pivotal elements of
+    //--  the simplex table in the ratio test
+    OCTAVE_GLPK_GET_REAL_PARAM ("tolpiv", 3);
+
+    OCTAVE_GLPK_GET_REAL_PARAM ("objll", 4);
+
+    OCTAVE_GLPK_GET_REAL_PARAM ("objul", 5);
+
+    OCTAVE_GLPK_GET_REAL_PARAM ("tmlim", 6);
+
+    OCTAVE_GLPK_GET_REAL_PARAM ("outdly", 7);
+
+    OCTAVE_GLPK_GET_REAL_PARAM ("tolint", 8);
+
+    OCTAVE_GLPK_GET_REAL_PARAM ("tolobj", 9);
+
+    OCTAVE_GLPK_GET_REAL_PARAM ("mipgap", 10);
+
+    //-- Assign pointers to the output parameters
+    ColumnVector xmin (mrowsc, octave_NA);
+    double fmin = octave_NA;
+    double status;
+    ColumnVector lambda (mrowsA, octave_NA);
+    ColumnVector redcosts (mrowsc, octave_NA);
+    double time;
+    double mem;
+
+    int jmpret = setjmp (mark);
+
+    if (jmpret == 0)
+        glpk (sense, mrowsc, mrowsA, c, nz, rn.fortran_vec (),
+              cn.fortran_vec (), a.fortran_vec (), b, ctype,
+              freeLB.fortran_vec (), lb, freeUB.fortran_vec (), ub,
+              vartype.fortran_vec (), isMIP, lpsolver, save_pb, filename.c_str(), filetype.c_str(),
+              xmin.fortran_vec (), &fmin, &status, lambda.fortran_vec (),
+              redcosts.fortran_vec (), &time, &mem);
+
+    octave_scalar_map extra;
+
+    if (! isMIP)
+    {
+        extra.assign ("lambda", lambda);
+        extra.assign ("redcosts", redcosts);
+    }
+
+    extra.assign ("time", time);
+    extra.assign ("mem", mem);
+
+    retval(3) = extra;
+    retval(2) = status;
+    retval(1) = fmin;
+    retval(0) = xmin;
 
 #else
 
-  gripe_not_supported ("glpk");
+    gripe_not_supported ("glpk");
 
 #endif
 
-  return retval;
+    return retval;
 }
