24 #if defined (HAVE_CONFIG_H)
40 #if defined (HAVE_GLPK)
44 #if defined (HAVE_GLPK_GLPK_H)
45 # include <glpk/glpk.h>
74 glpk (
int sense,
int n,
int m,
double *
c,
int nz,
int *rn,
int *cn,
75 double *
a,
double *
b,
char *ctype,
int *freeLB,
double *lb,
76 int *freeUB,
double *ub,
int *vartype,
int isMIP,
int lpsolver,
78 double *
xmin,
double *fmin,
int *status,
79 double *lambda,
double *redcosts,
double *time)
84 clock_t t_start = clock ();
86 glp_prob *lp = glp_create_prob ();
90 glp_set_obj_dir (lp, GLP_MIN);
92 glp_set_obj_dir (lp, GLP_MAX);
95 for (
int i = 0;
i < n;
i++)
98 if (! freeLB[
i] && ! freeUB[
i])
101 glp_set_col_bnds (lp, i+1, GLP_DB, lb[i], ub[i]);
103 glp_set_col_bnds (lp, i+1, GLP_FX, lb[i], ub[i]);
107 if (! freeLB[i] && freeUB[i])
108 glp_set_col_bnds (lp, i+1, GLP_LO, lb[i], ub[i]);
111 if (freeLB[i] && ! freeUB[i])
112 glp_set_col_bnds (lp, i+1, GLP_UP, lb[i], ub[i]);
114 glp_set_col_bnds (lp, i+1, GLP_FR, lb[i], ub[i]);
120 glp_set_obj_coef(lp,i+1,c[i]);
123 glp_set_col_kind (lp, i+1, vartype[i]);
126 glp_add_rows (lp, m);
128 for (
int i = 0;
i <
m;
i++)
159 glp_set_row_bnds (lp, i+1, typx, b[i], b[i]);
163 glp_load_matrix (lp, nz, rn, cn, a);
167 static char tmp[] =
"outpb.lp";
168 if (glp_write_lp (lp, 0, tmp) != 0)
169 error (
"__glpk__: unable to write problem");
173 if (! par->
presol || lpsolver != 1)
174 glp_scale_prob (lp, scale);
177 if (lpsolver == 1 && ! par->
presol)
178 glp_adv_basis (lp, 0);
182 if ((! isMIP && lpsolver == 1)
183 || (isMIP && ! par->
presol))
186 glp_init_smcp (&smcp);
187 smcp.msg_lev = par->
msglev;
188 smcp.meth = par->
dual;
189 smcp.pricing = par->
price;
190 smcp.r_test = par->
rtest;
191 smcp.tol_bnd = par->
tolbnd;
192 smcp.tol_dj = par->
toldj;
193 smcp.tol_piv = par->
tolpiv;
194 smcp.obj_ll = par->
objll;
195 smcp.obj_ul = par->
objul;
196 smcp.it_lim = par->
itlim;
197 smcp.tm_lim = par->
tmlim;
198 smcp.out_frq = par->
outfrq;
199 smcp.out_dly = par->
outdly;
200 smcp.presolve = par->
presol;
201 errnum = glp_simplex (lp, &smcp);
207 glp_init_iocp (&iocp);
208 iocp.msg_lev = par->
msglev;
209 iocp.br_tech = par->
branch;
210 iocp.bt_tech = par->
btrack;
211 iocp.tol_int = par->
tolint;
212 iocp.tol_obj = par->
tolobj;
213 iocp.tm_lim = par->
tmlim;
214 iocp.out_frq = par->
outfrq;
215 iocp.out_dly = par->
outdly;
216 iocp.presolve = par->
presol;
217 errnum = glp_intopt (lp, &iocp);
220 if (! isMIP && lpsolver == 2)
223 glp_init_iptcp (&iptcp);
224 iptcp.msg_lev = par->
msglev;
225 errnum = glp_interior (lp, &iptcp);
232 *status = glp_mip_status (lp);
233 *fmin = glp_mip_obj_val (lp);
239 *status = glp_get_status (lp);
240 *fmin = glp_get_obj_val (lp);
244 *status = glp_ipt_status (lp);
245 *fmin = glp_ipt_obj_val (lp);
251 for (
int i = 0;
i < n;
i++)
252 xmin[
i] = glp_mip_col_val (lp,
i+1);
257 for (
int i = 0;
i < n;
i++)
260 xmin[
i] = glp_get_col_prim (lp,
i+1);
262 xmin[
i] = glp_ipt_col_prim (lp,
i+1);
266 for (
int i = 0;
i <
m;
i++)
269 lambda[
i] = glp_get_row_dual (lp,
i+1);
271 lambda[
i] = glp_ipt_row_dual (lp,
i+1);
275 for (
int i = 0;
i < glp_get_num_cols (lp);
i++)
278 redcosts[
i] = glp_get_col_dual (lp,
i+1);
280 redcosts[
i] = glp_ipt_col_dual (lp,
i+1);
284 *time = (clock () - t_start) / CLOCKS_PER_SEC;
287 glp_delete_prob (lp);
299 #define OCTAVE_GLPK_GET_REAL_PARAM(NAME, VAL) \
302 octave_value tmp = PARAM.getfield (NAME); \
304 if (tmp.is_defined ()) \
306 if (! tmp.is_empty ()) \
307 VAL = tmp.xscalar_value ("glpk: invalid value in PARAM" NAME); \
309 error ("glpk: invalid value in PARAM" NAME); \
314 #define OCTAVE_GLPK_GET_INT_PARAM(NAME, VAL) \
317 octave_value tmp = PARAM.getfield (NAME); \
319 if (tmp.is_defined ()) \
321 if (! tmp.is_empty ()) \
322 VAL = tmp.xint_value ("glpk: invalid value in PARAM" NAME); \
324 error ("glpk: invalid value in PARAM" NAME); \
335 #if defined (HAVE_GLPK)
338 if (
args.length () != 9)
342 int mrowsc =
args(0).rows ();
344 Matrix C =
args(0).xmatrix_value (
"__glpk__: invalid value of C");
355 if (
args(1).is_sparse_type ())
367 error (
"__glpk__: invalid value of A");
373 rn(nz) = A.
ridx (
i) + 1;
380 Matrix A =
args(1).xmatrix_value (
"__glpk__: invalid value of A");
385 a.
resize (mrowsA*mrowsc+1, 0.0);
387 for (
int i = 0;
i < mrowsA;
i++)
389 for (
int j = 0; j < mrowsc; j++)
405 Matrix B =
args(2).xmatrix_value (
"__glpk__: invalid value of B");
411 Matrix LB =
args(3).xmatrix_value (
"__glpk__: invalid value of LB");
413 if (LB.
numel () < mrowsc)
414 error (
"__glpk__: invalid dimensions for LB");
420 for (
int i = 0;
i < mrowsc;
i++)
433 Matrix UB =
args(4).xmatrix_value (
"__glpk__: invalid value of UB");
435 if (UB.
numel () < mrowsc)
436 error (
"__glpk__: invalid dimensions for UB");
441 for (
int i = 0;
i < mrowsc;
i++)
454 charMatrix CTYPE =
args(5).xchar_matrix_value (
"__glpk__: invalid value of CTYPE");
459 charMatrix VTYPE =
args(6).xchar_matrix_value (
"__glpk__: invalid value of VARTYPE");
463 for (
int i = 0;
i < mrowsc ;
i++)
465 if (VTYPE(
i,0) ==
'I')
476 double SENSE =
args(7).xscalar_value (
"__glpk__: invalid value of SENSE");
494 error (
"__glpk__: PARAM.msglev must be 0 (no output) or 1 (error and warning messages only [default]) or 2 (normal output) or 3 (full output)");
499 if (scale < 0 || scale > 128)
500 error (
"__glpk__: PARAM.scale must either be 128 (automatic selection of scaling options), or a bitwise or of: 1 (geometric mean scaling), 16 (equilibration scaling), 32 (round scale factors to power of two), 64 (skip if problem is well scaled");
506 error (
"__glpk__: PARAM.dual must be 1 (use two-phase primal simplex [default]) or 2 (use two-phase dual simplex) or 3 (use two-phase dual simplex, and if it fails, switch to the primal simplex)");
512 error (
"__glpk__: PARAM.price must be 17 (textbook pricing) or 34 (steepest edge pricing [default])");
526 error (
"__glpk__: PARAM.branch must be 1 (first fractional variable) or 2 (last fractional variable) or 3 (most fractional variable) or 4 (heuristic by Driebeck and Tomlin [default]) or 5 (hybrid pseudocost heuristic)");
532 error (
"__glpk__: PARAM.btrack must be 1 (depth first search) or 2 (breadth first search) or 3 (best local bound) or 4 (best projection heuristic [default]");
538 error (
"__glpk__: PARAM.presol must be 0 (do NOT use LP presolver) or 1 (use LP presolver [default])");
543 if (lpsolver < 1 || lpsolver > 2)
544 error (
"__glpk__: PARAM.lpsolver must be 1 (simplex method) or 2 (interior point method)");
550 error (
"__glpk__: PARAM.rtest must be 17 (standard ratio test) or 34 (Harris' two-pass ratio test [default])");
561 save_pb = save_pb != 0;
612 extra.
assign (
"lambda", lambda);
613 extra.assign (
"redcosts", redcosts);
616 extra.assign (
"time", time);
617 extra.assign (
"status", status);
619 return ovl (xmin, fmin, errnum, extra);
623 octave_unused_parameter (
args);
octave_idx_type cols(void) const
octave_idx_type rows(void) const
OCTAVE_EXPORT octave_value_list isa nd deftypefn *return ovl(args(0).is_integer_type())
OCTINTERP_API void print_usage(void)
octave_idx_type numel(void) const
Number of elements in the array.
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, int scale, const control_params *par, double *xmin, double *fmin, int *status, double *lambda, double *redcosts, double *time)
void error(const char *fmt,...)
octave_idx_type * cidx(void)
octave_idx_type rows(void) const
octave_idx_type nnz(void) const
Actual number of nonzero terms.
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
F77_RET_T const F77_INT F77_CMPLX const F77_INT F77_CMPLX * B
#define OCTAVE_GLPK_GET_INT_PARAM(NAME, VAL)
nd deftypefn *octave_map m
#define OCTAVE_GLPK_GET_REAL_PARAM(NAME, VAL)
void resize(const dim_vector &dv, const T &rfv)
the sparsity preserving column transformation such that that defines the pivoting threshold can be given in which case it defines the c
octave_int< T > xmin(const octave_int< T > &x, const octave_int< T > &y)
charNDArray max(char d, const charNDArray &m)
octave_idx_type * ridx(void)
=val(i)}if ode{val(i)}occurs in table i
void assign(const std::string &k, const octave_value &val)
void scale(Matrix &m, double x, double y, double z)
#define DEFUN_DLD(name, args_name, nargout_name, doc)
const T * fortran_vec(void) const
Vector representing the dimensions (size) of an Array.
void err_disabled_feature(const std::string &fcn, const std::string &feature, const std::string &pkg)
void resize(octave_idx_type n, const double &rfv=0)
F77_RET_T const F77_INT F77_CMPLX * A