GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
__glpk__.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2005-2013 Nicolo' Giorgetti
4 Copyright (C) 2013 Sébastien Villemot <sebastien@debian.org>
5 
6 This file is part of Octave.
7 
8 Octave is free software; you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by the
10 Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12 
13 Octave is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 for more details.
17 
18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, see
20 <http://www.gnu.org/licenses/>.
21 
22 */
23 
24 #ifdef HAVE_CONFIG_H
25 #include <config.h>
26 #endif
27 
28 #include <cfloat>
29 #include <csetjmp>
30 #include <ctime>
31 
32 #include "lo-ieee.h"
33 
34 #include "defun-dld.h"
35 #include "error.h"
36 #include "gripes.h"
37 #include "oct-map.h"
38 #include "oct-obj.h"
39 #include "pager.h"
40 
41 #if defined (HAVE_GLPK)
42 
43 extern "C"
44 {
45 #if defined (HAVE_GLPK_GLPK_H)
46 #include <glpk/glpk.h>
47 #else
48 #include <glpk.h>
49 #endif
50 }
51 
53 {
54  int msglev;
55  int dual;
56  int price;
57  int itlim;
58  int outfrq;
59  int branch;
60  int btrack;
61  int presol;
62  int rtest;
63  int tmlim;
64  int outdly;
65  double tolbnd;
66  double toldj;
67  double tolpiv;
68  double objll;
69  double objul;
70  double tolint;
71  double tolobj;
72 };
73 
74 static jmp_buf mark; //-- Address for long jump to jump to
75 
76 int
77 glpk (int sense, int n, int m, double *c, int nz, int *rn, int *cn,
78  double *a, double *b, char *ctype, int *freeLB, double *lb,
79  int *freeUB, double *ub, int *vartype, int isMIP, int lpsolver,
80  int save_pb, int scale, const control_params *par,
81  double *xmin, double *fmin, int *status,
82  double *lambda, double *redcosts, double *time)
83 {
84  int typx = 0;
85  int errnum = 0;
86 
87  clock_t t_start = clock ();
88 
89  glp_prob *lp = glp_create_prob ();
90 
91  //-- Set the sense of optimization
92  if (sense == 1)
93  glp_set_obj_dir (lp, GLP_MIN);
94  else
95  glp_set_obj_dir (lp, GLP_MAX);
96 
97  glp_add_cols (lp, n);
98  for (int i = 0; i < n; i++)
99  {
100  //-- Define type of the structural variables
101  if (! freeLB[i] && ! freeUB[i])
102  {
103  if (lb[i] != ub[i])
104  glp_set_col_bnds (lp, i+1, GLP_DB, lb[i], ub[i]);
105  else
106  glp_set_col_bnds (lp, i+1, GLP_FX, lb[i], ub[i]);
107  }
108  else
109  {
110  if (! freeLB[i] && freeUB[i])
111  glp_set_col_bnds (lp, i+1, GLP_LO, lb[i], ub[i]);
112  else
113  {
114  if (freeLB[i] && ! freeUB[i])
115  glp_set_col_bnds (lp, i+1, GLP_UP, lb[i], ub[i]);
116  else
117  glp_set_col_bnds (lp, i+1, GLP_FR, lb[i], ub[i]);
118  }
119  }
120 
121  // -- Set the objective coefficient of the corresponding
122  // -- structural variable. No constant term is assumed.
123  glp_set_obj_coef(lp,i+1,c[i]);
124 
125  if (isMIP)
126  glp_set_col_kind (lp, i+1, vartype[i]);
127  }
128 
129  glp_add_rows (lp, m);
130 
131  for (int i = 0; i < m; i++)
132  {
133  /* If the i-th row has no lower bound (types F,U), the
134  corrispondent parameter will be ignored.
135  If the i-th row has no upper bound (types F,L), the corrispondent
136  parameter will be ignored.
137  If the i-th row is of S type, the i-th LB is used, but
138  the i-th UB is ignored.
139  */
140 
141  switch (ctype[i])
142  {
143  case 'F':
144  typx = GLP_FR;
145  break;
146 
147  case 'U':
148  typx = GLP_UP;
149  break;
150 
151  case 'L':
152  typx = GLP_LO;
153  break;
154 
155  case 'S':
156  typx = GLP_FX;
157  break;
158 
159  case 'D':
160  typx = GLP_DB;
161  break;
162  }
163 
164  glp_set_row_bnds (lp, i+1, typx, b[i], b[i]);
165 
166  }
167 
168  glp_load_matrix (lp, nz, rn, cn, a);
169 
170  if (save_pb)
171  {
172  static char tmp[] = "outpb.lp";
173  if (glp_write_lp (lp, NULL, tmp) != 0)
174  {
175  error ("__glpk__: unable to write problem");
176  longjmp (mark, -1);
177  }
178  }
179 
180  //-- scale the problem data
181  if (!par->presol || lpsolver != 1)
182  glp_scale_prob (lp, scale);
183 
184  //-- build advanced initial basis (if required)
185  if (lpsolver == 1 && !par->presol)
186  glp_adv_basis (lp, 0);
187 
188  /* For MIP problems without a presolver, a first pass with glp_simplex
189  is required */
190  if ((!isMIP && lpsolver == 1)
191  || (isMIP && !par->presol))
192  {
193  glp_smcp smcp;
194  glp_init_smcp (&smcp);
195  smcp.msg_lev = par->msglev;
196  smcp.meth = par->dual;
197  smcp.pricing = par->price;
198  smcp.r_test = par->rtest;
199  smcp.tol_bnd = par->tolbnd;
200  smcp.tol_dj = par->toldj;
201  smcp.tol_piv = par->tolpiv;
202  smcp.obj_ll = par->objll;
203  smcp.obj_ul = par->objul;
204  smcp.it_lim = par->itlim;
205  smcp.tm_lim = par->tmlim;
206  smcp.out_frq = par->outfrq;
207  smcp.out_dly = par->outdly;
208  smcp.presolve = par->presol;
209  errnum = glp_simplex (lp, &smcp);
210  }
211 
212  if (isMIP)
213  {
214  glp_iocp iocp;
215  glp_init_iocp (&iocp);
216  iocp.msg_lev = par->msglev;
217  iocp.br_tech = par->branch;
218  iocp.bt_tech = par->btrack;
219  iocp.tol_int = par->tolint;
220  iocp.tol_obj = par->tolobj;
221  iocp.tm_lim = par->tmlim;
222  iocp.out_frq = par->outfrq;
223  iocp.out_dly = par->outdly;
224  iocp.presolve = par->presol;
225  errnum = glp_intopt (lp, &iocp);
226  }
227 
228  if (!isMIP && lpsolver == 2)
229  {
230  glp_iptcp iptcp;
231  glp_init_iptcp (&iptcp);
232  iptcp.msg_lev = par->msglev;
233  errnum = glp_interior (lp, &iptcp);
234  }
235 
236  if (errnum == 0)
237  {
238  if (isMIP)
239  {
240  *status = glp_mip_status (lp);
241  *fmin = glp_mip_obj_val (lp);
242  }
243  else
244  {
245  if (lpsolver == 1)
246  {
247  *status = glp_get_status (lp);
248  *fmin = glp_get_obj_val (lp);
249  }
250  else
251  {
252  *status = glp_ipt_status (lp);
253  *fmin = glp_ipt_obj_val (lp);
254  }
255  }
256 
257  if (isMIP)
258  {
259  for (int i = 0; i < n; i++)
260  xmin[i] = glp_mip_col_val (lp, i+1);
261  }
262  else
263  {
264  /* Primal values */
265  for (int i = 0; i < n; i++)
266  {
267  if (lpsolver == 1)
268  xmin[i] = glp_get_col_prim (lp, i+1);
269  else
270  xmin[i] = glp_ipt_col_prim (lp, i+1);
271  }
272 
273  /* Dual values */
274  for (int i = 0; i < m; i++)
275  {
276  if (lpsolver == 1)
277  lambda[i] = glp_get_row_dual (lp, i+1);
278  else
279  lambda[i] = glp_ipt_row_dual (lp, i+1);
280  }
281 
282  /* Reduced costs */
283  for (int i = 0; i < glp_get_num_cols (lp); i++)
284  {
285  if (lpsolver == 1)
286  redcosts[i] = glp_get_col_dual (lp, i+1);
287  else
288  redcosts[i] = glp_ipt_col_dual (lp, i+1);
289  }
290  }
291 
292  *time = (clock () - t_start) / CLOCKS_PER_SEC;
293  }
294 
295  glp_delete_prob (lp);
296 
297  return errnum;
298 }
299 
300 #endif
301 
302 #define OCTAVE_GLPK_GET_REAL_PARAM(NAME, VAL) \
303  do \
304  { \
305  octave_value tmp = PARAM.getfield (NAME); \
306  \
307  if (tmp.is_defined ()) \
308  { \
309  if (! tmp.is_empty ()) \
310  { \
311  VAL = tmp.scalar_value (); \
312  \
313  if (error_state) \
314  { \
315  error ("glpk: invalid value in PARAM." NAME); \
316  return retval; \
317  } \
318  } \
319  else \
320  { \
321  error ("glpk: invalid value in PARAM." NAME); \
322  return retval; \
323  } \
324  } \
325  } \
326  while (0)
327 
328 #define OCTAVE_GLPK_GET_INT_PARAM(NAME, VAL) \
329  do \
330  { \
331  octave_value tmp = PARAM.getfield (NAME); \
332  \
333  if (tmp.is_defined ()) \
334  { \
335  if (! tmp.is_empty ()) \
336  { \
337  VAL = tmp.int_value (); \
338  \
339  if (error_state) \
340  { \
341  error ("glpk: invalid value in PARAM." NAME); \
342  return retval; \
343  } \
344  } \
345  else \
346  { \
347  error ("glpk: invalid value in PARAM." NAME); \
348  return retval; \
349  } \
350  } \
351  } \
352  while (0)
353 
354 DEFUN_DLD (__glpk__, args, ,
355  "-*- texinfo -*-\n\
356 @deftypefn {Loadable Function} {[@var{values}] =} __glpk__ (@var{args})\n\
357 Undocumented internal function.\n\
358 @end deftypefn")
359 {
360  // The list of values to return. See the declaration in oct-obj.h
361  octave_value_list retval;
362 
363 #if defined (HAVE_GLPK)
364 
365  int nrhs = args.length ();
366 
367  if (nrhs != 9)
368  {
369  print_usage ();
370  return retval;
371  }
372 
373  //-- 1nd Input. A column array containing the objective function
374  //-- coefficients.
375  volatile int mrowsc = args(0).rows ();
376 
377  Matrix C (args(0).matrix_value ());
378 
379  if (error_state)
380  {
381  error ("__glpk__: invalid value of C");
382  return retval;
383  }
384 
385  double *c = C.fortran_vec ();
386  Array<int> rn;
387  Array<int> cn;
388  ColumnVector a;
389  volatile int mrowsA;
390  volatile int nz = 0;
391 
392  //-- 2nd Input. A matrix containing the constraints coefficients.
393  // If matrix A is NOT a sparse matrix
394  if (args(1).is_sparse_type ())
395  {
396  SparseMatrix A = args(1).sparse_matrix_value (); // get the sparse matrix
397 
398  if (error_state)
399  {
400  error ("__glpk__: invalid value of A");
401  return retval;
402  }
403 
404  mrowsA = A.rows ();
405  octave_idx_type Anc = A.cols ();
406  octave_idx_type Anz = A.nnz ();
407  rn.resize (dim_vector (Anz+1, 1));
408  cn.resize (dim_vector (Anz+1, 1));
409  a.resize (Anz+1, 0.0);
410 
411  if (Anc != mrowsc)
412  {
413  error ("__glpk__: invalid value of A");
414  return retval;
415  }
416 
417  for (octave_idx_type j = 0; j < Anc; j++)
418  for (octave_idx_type i = A.cidx (j); i < A.cidx (j+1); i++)
419  {
420  nz++;
421  rn(nz) = A.ridx (i) + 1;
422  cn(nz) = j + 1;
423  a(nz) = A.data(i);
424  }
425  }
426  else
427  {
428  Matrix A (args(1).matrix_value ()); // get the matrix
429 
430  if (error_state)
431  {
432  error ("__glpk__: invalid value of A");
433  return retval;
434  }
435 
436  mrowsA = A.rows ();
437  rn.resize (dim_vector (mrowsA*mrowsc+1, 1));
438  cn.resize (dim_vector (mrowsA*mrowsc+1, 1));
439  a.resize (mrowsA*mrowsc+1, 0.0);
440 
441  for (int i = 0; i < mrowsA; i++)
442  {
443  for (int j = 0; j < mrowsc; j++)
444  {
445  if (A(i,j) != 0)
446  {
447  nz++;
448  rn(nz) = i + 1;
449  cn(nz) = j + 1;
450  a(nz) = A(i,j);
451  }
452  }
453  }
454 
455  }
456 
457  //-- 3rd Input. A column array containing the right-hand side value
458  // for each constraint in the constraint matrix.
459  Matrix B (args(2).matrix_value ());
460 
461  if (error_state)
462  {
463  error ("__glpk__: invalid value of B");
464  return retval;
465  }
466 
467  double *b = B.fortran_vec ();
468 
469  //-- 4th Input. An array of length mrowsc containing the lower
470  //-- bound on each of the variables.
471  Matrix LB (args(3).matrix_value ());
472 
473  if (error_state || LB.length () < mrowsc)
474  {
475  error ("__glpk__: invalid value of LB");
476  return retval;
477  }
478 
479  double *lb = LB.fortran_vec ();
480 
481  //-- LB argument, default: Free
482  Array<int> freeLB (dim_vector (mrowsc, 1));
483  for (int i = 0; i < mrowsc; i++)
484  {
485  if (xisinf (lb[i]))
486  {
487  freeLB(i) = 1;
488  lb[i] = -octave_Inf;
489  }
490  else
491  freeLB(i) = 0;
492  }
493 
494  //-- 5th Input. An array of at least length numcols containing the upper
495  //-- bound on each of the variables.
496  Matrix UB (args(4).matrix_value ());
497 
498  if (error_state || UB.length () < mrowsc)
499  {
500  error ("__glpk__: invalid value of UB");
501  return retval;
502  }
503 
504  double *ub = UB.fortran_vec ();
505 
506  Array<int> freeUB (dim_vector (mrowsc, 1));
507  for (int i = 0; i < mrowsc; i++)
508  {
509  if (xisinf (ub[i]))
510  {
511  freeUB(i) = 1;
512  ub[i] = octave_Inf;
513  }
514  else
515  freeUB(i) = 0;
516  }
517 
518  //-- 6th Input. A column array containing the sense of each constraint
519  //-- in the constraint matrix.
520  charMatrix CTYPE (args(5).char_matrix_value ());
521 
522  if (error_state)
523  {
524  error ("__glpk__: invalid value of CTYPE");
525  return retval;
526  }
527 
528  char *ctype = CTYPE.fortran_vec ();
529 
530  //-- 7th Input. A column array containing the types of the variables.
531  charMatrix VTYPE (args(6).char_matrix_value ());
532 
533  if (error_state)
534  {
535  error ("__glpk__: invalid value of VARTYPE");
536  return retval;
537  }
538 
539  Array<int> vartype (dim_vector (mrowsc, 1));
540  volatile int isMIP = 0;
541  for (int i = 0; i < mrowsc ; i++)
542  {
543  if (VTYPE(i,0) == 'I')
544  {
545  isMIP = 1;
546  vartype(i) = GLP_IV;
547  }
548  else
549  vartype(i) = GLP_CV;
550  }
551 
552  //-- 8th Input. Sense of optimization.
553  volatile int sense;
554  double SENSE = args(7).scalar_value ();
555 
556  if (error_state)
557  {
558  error ("__glpk__: invalid value of SENSE");
559  return retval;
560  }
561 
562  if (SENSE >= 0)
563  sense = 1;
564  else
565  sense = -1;
566 
567  //-- 9th Input. A structure containing the control parameters.
568  octave_scalar_map PARAM = args(8).scalar_map_value ();
569 
570  if (error_state)
571  {
572  error ("__glpk__: invalid value of PARAM");
573  return retval;
574  }
575 
576  control_params par;
577 
578  //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
579  //-- Integer parameters
580  //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
581 
582  //-- Level of messages output by the solver
583  par.msglev = 1;
584  OCTAVE_GLPK_GET_INT_PARAM ("msglev", par.msglev);
585  if (par.msglev < 0 || par.msglev > 3)
586  {
587  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)");
588  return retval;
589  }
590 
591  //-- scaling option
592  volatile int scale = 16;
593  OCTAVE_GLPK_GET_INT_PARAM ("scale", scale);
594  if (scale < 0 || scale > 128)
595  {
596  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");
597  return retval;
598  }
599 
600  //-- Dual simplex option
601  par.dual = 1;
602  OCTAVE_GLPK_GET_INT_PARAM ("dual", par.dual);
603  if (par.dual < 1 || par.dual > 3)
604  {
605  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)");
606  return retval;
607  }
608 
609  //-- Pricing option
610  par.price = 34;
611  OCTAVE_GLPK_GET_INT_PARAM ("price", par.price);
612  if (par.price != 17 && par.price != 34)
613  {
614  error ("__glpk__: PARAM.price must be 17 (textbook pricing) or 34 (steepest edge pricing [default])");
615  return retval;
616  }
617 
618  //-- Simplex iterations limit
620  OCTAVE_GLPK_GET_INT_PARAM ("itlim", par.itlim);
621 
622  //-- Output frequency, in iterations
623  par.outfrq = 200;
624  OCTAVE_GLPK_GET_INT_PARAM ("outfrq", par.outfrq);
625 
626  //-- Branching heuristic option
627  par.branch = 4;
628  OCTAVE_GLPK_GET_INT_PARAM ("branch", par.branch);
629  if (par.branch < 1 || par.branch > 5)
630  {
631  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)");
632  return retval;
633  }
634 
635  //-- Backtracking heuristic option
636  par.btrack = 4;
637  OCTAVE_GLPK_GET_INT_PARAM ("btrack", par.btrack);
638  if (par.btrack < 1 || par.btrack > 4)
639  {
640  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]");
641  return retval;
642  }
643 
644  //-- Presolver option
645  par.presol = 1;
646  OCTAVE_GLPK_GET_INT_PARAM ("presol", par.presol);
647  if (par.presol < 0 || par.presol > 1)
648  {
649  error ("__glpk__: PARAM.presol must be 0 (do NOT use LP presolver) or 1 (use LP presolver [default])");
650  return retval;
651  }
652 
653  //-- LPsolver option
654  volatile int lpsolver = 1;
655  OCTAVE_GLPK_GET_INT_PARAM ("lpsolver", lpsolver);
656  if (lpsolver < 1 || lpsolver > 2)
657  {
658  error ("__glpk__: PARAM.lpsolver must be 1 (simplex method) or 2 (interior point method)");
659  return retval;
660  }
661 
662  //-- Ratio test option
663  par.rtest = 34;
664  OCTAVE_GLPK_GET_INT_PARAM ("rtest", par.rtest);
665  if (par.rtest != 17 && par.rtest != 34)
666  {
667  error ("__glpk__: PARAM.rtest must be 17 (standard ratio test) or 34 (Harris' two-pass ratio test [default])");
668  return retval;
669  }
670 
672  OCTAVE_GLPK_GET_INT_PARAM ("tmlim", par.tmlim);
673 
674  par.outdly = 0;
675  OCTAVE_GLPK_GET_INT_PARAM ("outdly", par.outdly);
676 
677  //-- Save option
678  volatile int save_pb = 0;
679  OCTAVE_GLPK_GET_INT_PARAM ("save", save_pb);
680  save_pb = save_pb != 0;
681 
682  //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
683  //-- Real parameters
684  //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
685 
686  //-- Relative tolerance used to check if the current basic solution
687  //-- is primal feasible
688  par.tolbnd = 1e-7;
689  OCTAVE_GLPK_GET_REAL_PARAM ("tolbnd", par.tolbnd);
690 
691  //-- Absolute tolerance used to check if the current basic solution
692  //-- is dual feasible
693  par.toldj = 1e-7;
694  OCTAVE_GLPK_GET_REAL_PARAM ("toldj", par.toldj);
695 
696  //-- Relative tolerance used to choose eligible pivotal elements of
697  //-- the simplex table in the ratio test
698  par.tolpiv = 1e-10;
699  OCTAVE_GLPK_GET_REAL_PARAM ("tolpiv", par.tolpiv);
700 
702  OCTAVE_GLPK_GET_REAL_PARAM ("objll", par.objll);
703 
705  OCTAVE_GLPK_GET_REAL_PARAM ("objul", par.objul);
706 
707  par.tolint = 1e-5;
708  OCTAVE_GLPK_GET_REAL_PARAM ("tolint", par.tolint);
709 
710  par.tolobj = 1e-7;
711  OCTAVE_GLPK_GET_REAL_PARAM ("tolobj", par.tolobj);
712 
713  //-- Assign pointers to the output parameters
714  ColumnVector xmin (mrowsc, octave_NA);
715  double fmin = octave_NA;
716  ColumnVector lambda (mrowsA, octave_NA);
717  ColumnVector redcosts (mrowsc, octave_NA);
718  double time;
719  int status, errnum = 0;
720 
721  int jmpret = setjmp (mark);
722 
723  if (jmpret == 0)
724  errnum = glpk (sense, mrowsc, mrowsA, c, nz, rn.fortran_vec (),
725  cn.fortran_vec (), a.fortran_vec (), b, ctype,
726  freeLB.fortran_vec (), lb, freeUB.fortran_vec (), ub,
727  vartype.fortran_vec (), isMIP, lpsolver, save_pb, scale,
728  &par, xmin.fortran_vec (), &fmin, &status,
729  lambda.fortran_vec (), redcosts.fortran_vec (), &time);
730 
731  octave_scalar_map extra;
732 
733  if (! isMIP)
734  {
735  extra.assign ("lambda", lambda);
736  extra.assign ("redcosts", redcosts);
737  }
738 
739  extra.assign ("time", time);
740  extra.assign ("status", status);
741 
742  retval(3) = extra;
743  retval(2) = errnum;
744  retval(1) = fmin;
745  retval(0) = xmin;
746 
747 #else
748 
749  gripe_not_supported ("glpk");
750 
751 #endif
752 
753  return retval;
754 }
755 
756 /*
757 ## No test needed for internal helper function.
758 %!assert (1)
759 */