diff -r 61cd086e9aeb -r 94d2a3d746c9 scripts/optimization/glpk.m --- a/scripts/optimization/glpk.m Fri Apr 01 03:52:05 2011 -0600 +++ b/scripts/optimization/glpk.m Sat Apr 02 13:53:58 2011 -0600 @@ -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 (@address@hidden, default: 0) +## @item dual (@address@hidden, 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 (@address@hidden, default: 1) @@ -196,7 +208,16 @@ ## Steepest edge pricing. ## @end table ## -## @item round (@address@hidden, 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 (@address@hidden, 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 (@address@hidden, default: -1) +## @item itlim (@address@hidden, 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 (@address@hidden, 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 (@address@hidden) +## @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 (@address@hidden) -## Solution is feasible. -## @item 182 (@address@hidden) -## Solution is infeasible. -## @item 183 (@address@hidden) -## Problem has no feasible solution. -## @item 184 (@address@hidden) -## Problem has no unbounded solution. -## @item 185 (@address@hidden) -## Solution status is undefined. -## @end table -## Interior Point Method: -## @table @asis -## @item 150 (@address@hidden) -## The interior point method is undefined. -## @item 151 (@address@hidden) -## The interior point method is optimal. -## @end table -## Mixed Integer Method: -## @table @asis -## @item 170 (@address@hidden) -## The status is undefined. -## @item 171 (@address@hidden) -## The solution is integer optimal. -## @item 172 (@address@hidden) -## Solution integer feasible but its optimality has not been proven -## @item 173 (@address@hidden) -## 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 (@address@hidden) -## Unable to start the search. -## @item 205 (@address@hidden) -## Objective function lower limit reached. -## @item 206 (@address@hidden) -## Objective function upper limit reached. -## @item 207 (@address@hidden) -## Iterations limit exhausted. -## @item 208 (@address@hidden) -## Time limit exhausted. -## @item 209 (@address@hidden) -## No feasible solution. -## @item 210 (@address@hidden) -## Numerical instability. -## @item 211 (@address@hidden) -## Problems with basis matrix. -## @item 212 (@address@hidden) -## No convergence (interior). -## @item 213 (@address@hidden) -## No primal feasible solution (LP presolver). -## @item 214 (@address@hidden) -## No dual feasible solution (LP presolver). +## @item 101 (@address@hidden) +## 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 (@address@hidden) +## Unable to start the search, because the basis matrix corresponding +## to the initial basis is singular within the working precision. +## @item 103 (@address@hidden) +## 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 (@address@hidden) +## Unable to start the search, because some double-bounded (auxiliary or +## structural) variables have incorrect bounds. +## @item 105 (@address@hidden) +## The search was prematurely terminated due to the solver failure. +## @item 106 (@address@hidden) +## 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 (@address@hidden) +## 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 (@address@hidden) +## The search was prematurely terminated, because the simplex iteration +## limit has been exceeded. +## @item 109 (@address@hidden) +## The search was prematurely terminated, because the time limit has been +## exceeded. +## @item 110 (@address@hidden) +## The LP problem instance has no primal feasible solution (only if the LP +## presolver is used). +## @item 111 (@address@hidden) +## 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 (@address@hidden) +## Unable to start the search, because some double-bounded +## variables have incorrect bounds or some integer variables +## have non-integer (fractional) bounds. +## @item 205 (@address@hidden) +## The search was prematurely terminated due to the solver failure. +## @item 209 (@address@hidden) +## The search was prematurely terminated, because the time limit has been +## exceeded. +## @item 210 (@address@hidden) +## 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 (@address@hidden) +## 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 (@address@hidden) +## 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 (@address@hidden) +## The search was prematurely terminated, because the relative MIP gap +## tolerance has been reached. +## @end table +## +## Interior point method routines: +## @table @asis +## @item 305 (@address@hidden) +## The problem has no rows/columns. +## @item 308 (@address@hidden) +## Iteration limit exceeded. +## @item 316 (@address@hidden) +## Very slow convergence or divergence. +## @item 317 (@address@hidden) +## 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 61cd086e9aeb -r 94d2a3d746c9 src/DLD-FUNCTIONS/__glpk__.cc --- a/src/DLD-FUNCTIONS/__glpk__.cc Fri Apr 01 03:52:05 2011 -0600 +++ b/src/DLD-FUNCTIONS/__glpk__.cc Sat Apr 02 13:53:58 2011 -0600 @@ -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 #include #include +#include +#include #include "lo-ieee.h" @@ -39,378 +42,593 @@ #if defined (HAVE_GLPK) -extern "C" -{ -#if defined (HAVE_GLPK_GLPK_H) -#include -#else +extern "C" { #include -#endif - -#if 0 -#ifdef GLPK_PRE_4_14 - -#ifndef _GLPLIB_H -#include -#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(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(&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(glpRealParam[6]); + sParam.out_frq=glpIntParam[7]; // output frequency + sParam.out_dly=static_cast(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(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((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} address@hidden =} __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 rn; - Array cn; - ColumnVector a; - volatile int mrowsA; - volatile int nz = 0; + double *c = C.fortran_vec (); + Array rn; + Array 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 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 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 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 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 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 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; }