--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/dldfcn/__glpk__.cc	Tue Jul 31 21:57:58 2012 -0400
@@ -0,0 +1,861 @@
+/*
+
+Copyright (C) 2005-2012 Nicolo' Giorgetti
+
+This file is part of Octave.
+
+Octave 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 3 of the License, or (at your
+option) any later version.
+
+Octave 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 Octave; see the file COPYING.  If not, see
+<http://www.gnu.org/licenses/>.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <cfloat>
+#include <csetjmp>
+#include <ctime>
+
+#include "lo-ieee.h"
+
+#include "defun-dld.h"
+#include "error.h"
+#include "gripes.h"
+#include "oct-map.h"
+#include "oct-obj.h"
+#include "pager.h"
+
+#if defined (HAVE_GLPK)
+
+extern "C"
+{
+#if defined (HAVE_GLPK_GLPK_H)
+#include <glpk/glpk.h>
+#else
+#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
+
+int lpxIntParam[NIntP] = {
+  0,
+  1,
+  0,
+  1,
+  0,
+  -1,
+  0,
+  200,
+  1,
+  2,
+  0,
+  1,
+  0,
+  0,
+  2,
+  2,
+  1
+};
+
+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
+};
+
+
+double lpxRealParam[NRealP] = {
+  0.07,
+  1e-7,
+  1e-7,
+  1e-9,
+  -DBL_MAX,
+  DBL_MAX,
+  -1.0,
+  0.0,
+  1e-6,
+  1e-7
+};
+
+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
+};
+
+static jmp_buf mark;  //-- Address for long jump to jump to
+
+#if 0
+int
+glpk_fault_hook (void * /* info */, char *msg)
+{
+  error ("CRITICAL ERROR in GLPK: %s", msg);
+  longjmp (mark, -1);
+}
+
+int
+glpk_print_hook (void * /* info */, char *msg)
+{
+  message (0, "%s", msg);
+  return 1;
+}
+#endif
+
+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 ();
+
+#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
+
+  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
+
+  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 type of the structural variables
+      if (! freeLB[i] && ! freeUB[i])
+        {
+          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])
+                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]);
+            }
+        }
+
+      // -- Set the objective coefficient of the corresponding
+      // -- structural variable. No constant term is assumed.
+      lpx_set_obj_coef(lp,i+1,c[i]);
+
+      if (isMIP)
+        lpx_set_col_kind (lp, i+1, vartype[i]);
+    }
+
+  lpx_add_rows (lp, m);
+
+  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.
+      */
+
+      switch (ctype[i])
+        {
+        case 'F':
+          typx = LPX_FR;
+          break;
+
+        case 'U':
+          typx = LPX_UP;
+          break;
+
+        case 'L':
+          typx = LPX_LO;
+          break;
+
+        case 'S':
+          typx = LPX_FX;
+          break;
+
+        case 'D':
+          typx = LPX_DB;
+          break;
+        }
+
+      lpx_set_row_bnds (lp, i+1, typx, b[i], b[i]);
+
+    }
+
+  lpx_load_matrix (lp, nz, rn, cn, a);
+
+  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);
+            }
+          else
+            {
+              *status = lpx_ipt_status (lp);
+              *fmin = lpx_ipt_obj_val (lp);
+            }
+        }
+
+      if (isMIP)
+        {
+          for (int i = 0; i < n; i++)
+            xmin[i] = lpx_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);
+            }
+
+          /* 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);
+            }
+
+          /* 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);
+            }
+        }
+
+      *time = (clock () - t_start) / CLOCKS_PER_SEC;
+
+#ifdef GLPK_PRE_4_14
+      *mem = (lib_env_ptr () -> mem_tpeak);
+#else
+      *mem = 0;
+#endif
+
+      lpx_delete_prob (lp);
+      return 0;
+    }
+
+   lpx_delete_prob (lp);
+
+   *status = errnum;
+
+   return errnum;
+}
+
+#endif
+
+#define OCTAVE_GLPK_GET_REAL_PARAM(NAME, IDX) \
+  do \
+    { \
+      octave_value tmp = PARAM.getfield (NAME); \
+ \
+      if (tmp.is_defined ()) \
+        { \
+          if (! tmp.is_empty ()) \
+            { \
+              lpxRealParam[IDX] = tmp.scalar_value (); \
+ \
+              if (error_state) \
+                { \
+                  error ("glpk: invalid value in PARAM." NAME); \
+                  return retval; \
+                } \
+            } \
+          else \
+            { \
+              error ("glpk: invalid value in PARAM." NAME); \
+              return retval; \
+            } \
+        } \
+    } \
+  while (0)
+
+#define OCTAVE_GLPK_GET_INT_PARAM(NAME, VAL) \
+  do \
+    { \
+      octave_value tmp = PARAM.getfield (NAME); \
+ \
+      if (tmp.is_defined ()) \
+        { \
+          if (! tmp.is_empty ()) \
+            { \
+              VAL = tmp.int_value (); \
+ \
+              if (error_state) \
+                { \
+                  error ("glpk: invalid value in PARAM." NAME); \
+                  return retval; \
+                } \
+            } \
+          else \
+            { \
+              error ("glpk: invalid value in PARAM." NAME); \
+              return retval; \
+            } \
+        } \
+    } \
+  while (0)
+
+DEFUN_DLD (__glpk__, args, ,
+  "-*- 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;
+
+#if defined (HAVE_GLPK)
+
+  int nrhs = args.length ();
+
+  if (nrhs != 9)
+    {
+      print_usage ();
+      return retval;
+    }
+
+  //-- 1nd Input. A column array containing the objective function
+  //--            coefficients.
+  volatile int mrowsc = args(0).rows ();
+
+  Matrix C (args(0).matrix_value ());
+
+  if (error_state)
+    {
+      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;
+
+  //-- 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
+
+      if (error_state)
+        {
+          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);
+
+      if (Anc != mrowsc)
+        {
+          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);
+          }
+    }
+  else
+    {
+      Matrix A (args(1).matrix_value ()); // get the matrix
+
+      if (error_state)
+        {
+          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);
+
+      for (int i = 0; i < mrowsA; i++)
+        {
+          for (int j = 0; j < mrowsc; j++)
+            {
+              if (A(i,j) != 0)
+                {
+                  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 ());
+
+  if (error_state)
+    {
+      error ("__glpk__: invalid value of B");
+      return retval;
+    }
+
+  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 ());
+
+  if (error_state || LB.length () < mrowsc)
+    {
+      error ("__glpk__: invalid value of LB");
+      return retval;
+    }
+
+  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)
+    {
+      error ("__glpk__: invalid value of UB");
+      return retval;
+    }
+
+  double *ub = UB.fortran_vec ();
+
+  Array<int> freeUB (dim_vector (mrowsc, 1));
+  for (int i = 0; i < mrowsc; i++)
+    {
+      if (xisinf (ub[i]))
+        {
+          freeUB(i) = 1;
+          ub[i] = octave_Inf;
+        }
+      else
+        freeUB(i) = 0;
+    }
+
+  //-- 6th Input. A column array containing the sense of each constraint
+  //--            in the constraint matrix.
+  charMatrix CTYPE (args(5).char_matrix_value ());
+
+  if (error_state)
+    {
+      error ("__glpk__: invalid value of CTYPE");
+      return retval;
+    }
+
+  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)
+    {
+      error ("__glpk__: invalid value of VARTYPE");
+      return retval;
+    }
+
+  Array<int> vartype (dim_vector (mrowsc, 1));
+  volatile int isMIP = 0;
+  for (int i = 0; i < mrowsc ; i++)
+    {
+      if (VTYPE(i,0) == 'I')
+        {
+          isMIP = 1;
+          vartype(i) = LPX_IV;
+        }
+      else
+        vartype(i) = LPX_CV;
+    }
+
+  //-- 8th Input. Sense of optimization.
+  volatile int sense;
+  double SENSE = args(7).scalar_value ();
+
+  if (error_state)
+    {
+      error ("__glpk__: invalid value of SENSE");
+      return retval;
+    }
+
+  if (SENSE >= 0)
+    sense = 1;
+  else
+    sense = -1;
+
+  //-- 9th Input. A structure containing the control parameters.
+  octave_scalar_map PARAM = args(8).scalar_map_value ();
+
+  if (error_state)
+    {
+      error ("__glpk__: invalid value of PARAM");
+      return retval;
+    }
+
+  //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+  //-- Integer parameters
+  //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+  //-- Level of messages output by the solver
+  OCTAVE_GLPK_GET_INT_PARAM ("msglev", lpxIntParam[0]);
+  if (lpxIntParam[0] < 0 || lpxIntParam[0] > 3)
+    {
+      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;
+    }
+
+  //-- scaling option
+  OCTAVE_GLPK_GET_INT_PARAM ("scale", lpxIntParam[1]);
+  if (lpxIntParam[1] < 0 || lpxIntParam[1] > 2)
+    {
+      error ("__glpk__: PARAM.scale must be 0 (no scaling) or 1 (equilibration scaling [default]) or 2 (geometric mean scaling)");
+      return retval;
+    }
+
+  //-- Dual dimplex option
+  OCTAVE_GLPK_GET_INT_PARAM ("dual", lpxIntParam[2]);
+  if (lpxIntParam[2] < 0 || lpxIntParam[2] > 1)
+    {
+      error ("__glpk__: PARAM.dual must be 0 (do NOT use dual simplex [default]) or 1 (use dual simplex)");
+      return retval;
+    }
+
+  //-- Pricing option
+  OCTAVE_GLPK_GET_INT_PARAM ("price", lpxIntParam[3]);
+  if (lpxIntParam[3] < 0 || lpxIntParam[3] > 1)
+    {
+      error ("__glpk__: PARAM.price must be 0 (textbook pricing) or 1 (steepest edge pricing [default])");
+      return retval;
+    }
+
+  //-- Solution rounding option
+  OCTAVE_GLPK_GET_INT_PARAM ("round", lpxIntParam[4]);
+  if (lpxIntParam[4] < 0 || lpxIntParam[4] > 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;
+    }
+
+  //-- 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)
+    {
+      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;
+    }
+
+  //-- Backtracking heuristic option
+  OCTAVE_GLPK_GET_INT_PARAM ("btrack", lpxIntParam[15]);
+  if (lpxIntParam[15] < 0 || lpxIntParam[15] > 2)
+    {
+      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;
+    }
+
+  //-- Presolver option
+  OCTAVE_GLPK_GET_INT_PARAM ("presol", lpxIntParam[16]);
+  if (lpxIntParam[16] < 0 || lpxIntParam[16] > 1)
+    {
+      error ("__glpk__: PARAM.presol must be 0 (do NOT use LP presolver) or 1 (use LP presolver [default])");
+      return retval;
+    }
+
+  //-- LPsolver option
+  volatile int lpsolver = 1;
+  OCTAVE_GLPK_GET_INT_PARAM ("lpsolver", lpsolver);
+  if (lpsolver < 1 || lpsolver > 2)
+    {
+      error ("__glpk__: PARAM.lpsolver must be 1 (simplex method) or 2 (interior point method)");
+      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)
+    {
+      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");
+
+#endif
+
+  return retval;
+}
+
+/*
+## No test needed for internal helper function.
+%!assert (1)
+*/