diff -rc octave-old//scripts/optimization/glpk.m octave-diff//scripts/optimization/glpk.m *** octave-old//scripts/optimization/glpk.m 2011-02-08 11:00:51.000000000 +0100 --- octave-diff//scripts/optimization/glpk.m 2011-03-05 13:43:23.894275999 +0100 *************** *** 132,137 **** --- 132,140 ---- ## ## @item "I" ## An integer variable. + ## + ## @item "B" + ## A binary variable. ## @end table ## ## @item sense *************** *** 167,182 **** ## Scaling option: ## @table @asis ## @item 0 ! ## No scaling. ## ## @item 1 ! ## Equilibration scaling. ## ## @item 2 ! ## Geometric mean scaling, then equilibration scaling. ## @end table ## ! ## @item dual (@address@hidden, default: 0) ## Dual simplex option: ## @table @asis ## @item 0 --- 170,191 ---- ## Scaling option: ## @table @asis ## @item 0 ! ## Skip scaling. ## ## @item 1 ! ## Perform geometric mean scaling. ## ## @item 2 ! ## 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) ## Dual simplex option: ## @table @asis ## @item 0 *************** *** 184,189 **** --- 193,201 ---- ## ## @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,202 **** ## Steepest edge pricing. ## @end table ## ! ## @item round (@address@hidden, default: 0) ## Solution rounding option: ## @table @asis ## @item 0 --- 208,223 ---- ## Steepest edge pricing. ## @end table ## ! ## @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,212 **** ## Replace tiny primal and dual values by exact zero. ## @end table ## ! ## @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,233 ---- ## Replace tiny primal and dual values by exact zero. ## @end table ## ! ## @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,232 **** --- 248,256 ---- ## 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,245 **** --- 264,272 ---- ## Breadth first search. ## ## @item 2 + ## Best local bound. + ## + ## @item 3 ## Backtrack using the best projection heuristic. ## @end table ## *************** *** 256,267 **** ## ## @item 2 ## Interior point method. ## @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. ## @end table ## ## Real parameters: --- 283,378 ---- ## ## @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, 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,328 **** --- 434,467 ---- ## 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,403 **** ## @item status ## Status of the optimization. ## - ## Simplex Method: ## @table @asis ! ## @item 180 (@address@hidden) ! ## 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. ## @end table ## @noindent ## If an error occurs, @var{status} will contain one of the following ! ## codes: ## ## @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). ## @end table ## ## @item extra --- 477,582 ---- ## @item status ## Status of the optimization. ## ## @table @asis ! ## @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 6 ! ## Solution is unbounded. ## @end table ## @noindent ## If an error occurs, @var{status} will contain one of the following ! ## codes. ## + ## Simplex method routines: ## @table @asis ! ## @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,460 **** endif if (all (size (c) > 1) || iscomplex (c) || ischar (c)) ! error ("glpk:C must be a real vector"); return; endif nx = length (c); --- 633,639 ---- endif if (all (size (c) > 1) || iscomplex (c) || ischar (c)) ! error ("glpk: c must be a real vector"); return; endif nx = length (c); *************** *** 536,543 **** || 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"); return; endif else --- 715,722 ---- || length (vartype) != nx) error ("glpk: VARTYPE must be a char valued vector of length %d", nx); return; ! elseif (! all (vartype == "C" | vartype == "I" | vartype == "B")) ! error ("glpk: VARTYPE must contain only C, I or B"); return; endif else diff -rc octave-old//src/DLD-FUNCTIONS/__glpk__.cc octave-diff//src/DLD-FUNCTIONS/__glpk__.cc *** octave-old//src/DLD-FUNCTIONS/__glpk__.cc 2011-02-21 02:04:59.441776066 +0100 --- octave-diff//src/DLD-FUNCTIONS/__glpk__.cc 2011-03-05 13:42:54.226275998 +0100 *************** *** 1,5 **** --- 1,6 ---- /* + Copyright (C) 2011 Tommaso Balercia Copyright (C) 2005-2011 Nicolo' Giorgetti This file is part of Octave. *************** *** 27,32 **** --- 28,35 ---- #include #include #include + #include + #include #include "lo-ieee.h" *************** *** 39,416 **** #if defined (HAVE_GLPK) ! extern "C" ! { ! #if defined (HAVE_GLPK_GLPK_H) ! #include ! #else #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 ! 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 --- 42,634 ---- #if defined (HAVE_GLPK) ! extern "C" { #include } ! #define NIntP 21 ! #define NRealP 11 ! // 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_USECUTS, ! LPX_K_PREPROCESS, ! LPX_K_BINARIZE, ! LPX_K_RATIO_TEST }; ! //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_MIPGAP }; ! static jmp_buf mark; /*-- Address for long jump */ ! static int glpk_print_hook (void *info, const char *msg) { ! int msglev = *reinterpret_cast(info); ! if(msglev > 0) ! { ! message(0, "%s", msg); ! } ! return 1; } ! // ! 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) { ! int typx = 0; ! int method; ! 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)); ! ! //-- Create an empty LP/MILP object ! glp_prob *lp = glp_create_prob (); ! ! //-- Set the sense of optimization ! if (sense == 1) ! glp_set_obj_dir (lp, GLP_MIN); ! else ! glp_set_obj_dir (lp, GLP_MAX); ! ! //-- 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]) { ! 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 (! freeLB[i] && freeUB[i]) ! glp_set_col_bnds (lp, i+1, GLP_LO, lb[i], ub[i]); ! else { ! 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. ! glp_set_obj_coef(lp,i+1,c[i]); ! if (isMIP) ! glp_set_col_kind (lp, i+1, vartype[i]); } ! glp_add_rows (lp, m); ! for (int i = 0; i < m; i++) { ! /* 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]) { case 'F': ! typx = GLP_FR; ! break; ! // upper bound case 'U': ! typx = GLP_UP; ! break; ! // lower bound case 'L': ! typx = GLP_LO; ! break; ! // fixed constraint case 'S': ! typx = GLP_FX; ! break; ! // double-bounded variable case 'D': ! typx = GLP_DB; ! break; } ! 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); } ! 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 (lpsolver == 1 || lpsolver == 3) { ! *status = glp_get_status (lp); ! *fmin = glp_get_obj_val (lp); } ! else { ! *status = glp_ipt_status (lp); ! *fmin = glp_ipt_obj_val (lp); } } ! // 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 || 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 || 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 < 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; ! glp_long tpeak; ! glp_mem_usage(NULL, NULL, NULL, &tpeak); ! *mem=static_cast((4294967296.0 * tpeak.hi + tpeak.lo) / (1024)); ! glp_delete_prob (lp); ! ! return 0; ! } ! else { ! // printf("errnum is %d\n", errnum); } ! glp_delete_prob (lp); ! /* this shouldn't be nessiary with glp_deleted_prob, but try it ! if we have weird behavior again... */ ! glp_free_env(); ! ! *status = errnum; ! ! return errnum; } #endif *************** *** 424,430 **** { \ if (! tmp.is_empty ()) \ { \ ! lpxRealParam[IDX] = tmp.scalar_value (); \ \ if (error_state) \ { \ --- 642,648 ---- { \ if (! tmp.is_empty ()) \ { \ ! glpRealParam[IDX] = tmp.scalar_value (); \ \ if (error_state) \ { \ *************** *** 468,856 **** while (0) DEFUN_DLD (__glpk__, args, , ! "-*- 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; #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 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 ()) { ! 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 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 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 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; } --- 686,1170 ---- while (0) DEFUN_DLD (__glpk__, args, , ! "-*- 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; #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 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 ()) { ! 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 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 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 vartype (dim_vector (mrowsc, 1)); ! volatile int isMIP = 0; ! for (int i = 0; i < mrowsc ; i++) { ! 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; } } ! //-- 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", glpIntParam[0]); ! if (glpIntParam[0] < 0 || glpIntParam[0] > 3) { ! 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; } ! //-- scaling option ! OCTAVE_GLPK_GET_INT_PARAM ("scale", glpIntParam[1]); ! if (glpIntParam[1] < 0 || glpIntParam[1] > 4) { ! error ("glpk: the parameter \"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", glpIntParam[2]); ! if (glpIntParam[2] < 0 || glpIntParam[2] > 2) { ! error ("glpk: the parameters \"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", glpIntParam[3]); ! if (glpIntParam[3] < 0 || glpIntParam[3] > 1) { ! error ("glpk: the parameter \"price\" must be 0 (textbook pricing) or 1 (steepest edge pricing [default])"); ! return retval; } ! //-- Ratio test option ! OCTAVE_GLPK_GET_INT_PARAM ("r_test", glpIntParam[20]); ! if (glpIntParam[20] < 0 || glpIntParam[20] > 1) { ! error ("glpk: the parameter \"r_test\" must be 0 (textbook [default]) or 1 (Harris's Two pass ratio test)"); ! return retval; } ! //-- Solution rounding option ! OCTAVE_GLPK_GET_INT_PARAM ("round", glpIntParam[4]); ! if (glpIntParam[4] < 0 || glpIntParam[4] > 1) ! { ! 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; ! } ! //-- 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: 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; ! } ! ! //-- Backtracking heuristic option ! OCTAVE_GLPK_GET_INT_PARAM ("btrack", glpIntParam[15]); ! if (glpIntParam[15] < 0 || glpIntParam[15] > 3) ! { ! 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; ! } ! ! //-- Presolver option ! OCTAVE_GLPK_GET_INT_PARAM ("presol", glpIntParam[16]); ! if (glpIntParam[16] < 0 || glpIntParam[16] > 1) ! { ! error ("glpk: the parameter \"presol\" must be 0 (do NOT use LP presolver [default]) or 1 (use LP presolver [default])"); ! return retval; ! } ! ! //-- 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; } ! //-- 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"); #endif ! return retval; }