|
5234
|
1 /* |
|
|
2 |
|
|
3 Copyright (C) 2005 Nicolo' Giorgetti |
|
|
4 |
|
|
5 This file is part of Octave. |
|
|
6 |
|
|
7 Octave is free software; you can redistribute it and/or modify it |
|
|
8 under the terms of the GNU General Public License as published by the |
|
|
9 Free Software Foundation; either version 2, or (at your option) any |
|
|
10 later version. |
|
|
11 |
|
|
12 Octave is distributed in the hope that it will be useful, but WITHOUT |
|
|
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
|
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
|
15 for more details. |
|
|
16 |
|
|
17 You should have received a copy of the GNU General Public License |
|
|
18 along with Octave; see the file COPYING. If not, write to the Free |
|
5307
|
19 Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA |
|
|
20 02110-1301, USA. |
|
5234
|
21 |
|
|
22 */ |
|
|
23 |
|
|
24 #ifdef HAVE_CONFIG_H |
|
|
25 #include <config.h> |
|
|
26 #endif |
|
5232
|
27 |
|
|
28 #include <cfloat> |
|
|
29 #include <csetjmp> |
|
|
30 #include <ctime> |
|
|
31 |
|
5234
|
32 #include "defun-dld.h" |
|
|
33 #include "error.h" |
|
5235
|
34 #include "gripes.h" |
|
|
35 #include "oct-map.h" |
|
5234
|
36 #include "oct-obj.h" |
|
5235
|
37 #include "pager.h" |
|
|
38 |
|
|
39 #if defined (HAVE_GLPK) |
|
5232
|
40 |
|
6333
|
41 extern "C" |
|
|
42 { |
|
6804
|
43 #if defined (HAVE_GLPK_GLPK_H) |
|
|
44 #include <glpk/glpk.h> |
|
|
45 #else |
|
5234
|
46 #include <glpk.h> |
|
6804
|
47 #endif |
|
6333
|
48 |
|
6381
|
49 #ifdef GLPK_PRE_4_14 |
|
6333
|
50 |
|
6287
|
51 #ifndef _GLPLIB_H |
|
|
52 #include <glplib.h> |
|
|
53 #endif |
|
|
54 #ifndef lib_set_fault_hook |
|
|
55 #define lib_set_fault_hook lib_fault_hook |
|
|
56 #endif |
|
|
57 #ifndef lib_set_print_hook |
|
|
58 #define lib_set_print_hook lib_print_hook |
|
|
59 #endif |
|
5232
|
60 |
|
6333
|
61 #else |
|
|
62 |
|
|
63 void _glp_lib_print_hook (int (*func)(void *info, char *buf), void *info); |
|
|
64 void _glp_lib_fault_hook (int (*func)(void *info, char *buf), void *info); |
|
|
65 |
|
|
66 #endif |
|
6472
|
67 } |
|
6333
|
68 |
|
5232
|
69 #define NIntP 17 |
|
|
70 #define NRealP 10 |
|
|
71 |
|
5234
|
72 int lpxIntParam[NIntP] = { |
|
5240
|
73 0, |
|
5234
|
74 1, |
|
|
75 0, |
|
|
76 1, |
|
|
77 0, |
|
|
78 -1, |
|
|
79 0, |
|
|
80 200, |
|
|
81 1, |
|
|
82 2, |
|
|
83 0, |
|
|
84 1, |
|
|
85 0, |
|
|
86 0, |
|
|
87 2, |
|
|
88 2, |
|
|
89 1 |
|
5232
|
90 }; |
|
|
91 |
|
5234
|
92 int IParam[NIntP] = { |
|
|
93 LPX_K_MSGLEV, |
|
|
94 LPX_K_SCALE, |
|
|
95 LPX_K_DUAL, |
|
|
96 LPX_K_PRICE, |
|
|
97 LPX_K_ROUND, |
|
|
98 LPX_K_ITLIM, |
|
|
99 LPX_K_ITCNT, |
|
|
100 LPX_K_OUTFRQ, |
|
|
101 LPX_K_MPSINFO, |
|
|
102 LPX_K_MPSOBJ, |
|
|
103 LPX_K_MPSORIG, |
|
|
104 LPX_K_MPSWIDE, |
|
|
105 LPX_K_MPSFREE, |
|
|
106 LPX_K_MPSSKIP, |
|
|
107 LPX_K_BRANCH, |
|
|
108 LPX_K_BTRACK, |
|
|
109 LPX_K_PRESOL |
|
5232
|
110 }; |
|
|
111 |
|
|
112 |
|
5234
|
113 double lpxRealParam[NRealP] = { |
|
|
114 0.07, |
|
|
115 1e-7, |
|
|
116 1e-7, |
|
|
117 1e-9, |
|
|
118 -DBL_MAX, |
|
|
119 DBL_MAX, |
|
|
120 -1.0, |
|
|
121 0.0, |
|
|
122 1e-6, |
|
|
123 1e-7 |
|
5232
|
124 }; |
|
|
125 |
|
5234
|
126 int RParam[NRealP] = { |
|
|
127 LPX_K_RELAX, |
|
|
128 LPX_K_TOLBND, |
|
|
129 LPX_K_TOLDJ, |
|
|
130 LPX_K_TOLPIV, |
|
|
131 LPX_K_OBJLL, |
|
|
132 LPX_K_OBJUL, |
|
|
133 LPX_K_TMLIM, |
|
|
134 LPX_K_OUTDLY, |
|
|
135 LPX_K_TOLINT, |
|
|
136 LPX_K_TOLOBJ |
|
5232
|
137 }; |
|
|
138 |
|
5241
|
139 static jmp_buf mark; //-- Address for long jump to jump to |
|
5232
|
140 |
|
5234
|
141 int |
|
5235
|
142 glpk_fault_hook (void * /* info */, char *msg) |
|
5234
|
143 { |
|
5240
|
144 error ("CRITICAL ERROR in GLPK: %s", msg); |
|
5234
|
145 longjmp (mark, -1); |
|
5232
|
146 } |
|
|
147 |
|
5234
|
148 int |
|
5235
|
149 glpk_print_hook (void * /* info */, char *msg) |
|
5232
|
150 { |
|
5241
|
151 message (0, "%s", msg); |
|
5234
|
152 return 1; |
|
5232
|
153 } |
|
|
154 |
|
5234
|
155 int |
|
|
156 glpk (int sense, int n, int m, double *c, int nz, int *rn, int *cn, |
|
|
157 double *a, double *b, char *ctype, int *freeLB, double *lb, |
|
|
158 int *freeUB, double *ub, int *vartype, int isMIP, int lpsolver, |
|
|
159 int save_pb, double *xmin, double *fmin, double *status, |
|
|
160 double *lambda, double *redcosts, double *time, double *mem) |
|
5232
|
161 { |
|
5240
|
162 int errnum; |
|
5234
|
163 int typx = 0; |
|
|
164 int method; |
|
|
165 |
|
|
166 clock_t t_start = clock(); |
|
|
167 |
|
6381
|
168 #ifdef GLPK_PRE_4_14 |
|
6333
|
169 lib_set_fault_hook (0, glpk_fault_hook); |
|
|
170 #else |
|
|
171 _glp_lib_fault_hook (glpk_fault_hook, 0); |
|
|
172 #endif |
|
5232
|
173 |
|
5234
|
174 if (lpxIntParam[0] > 1) |
|
6381
|
175 #ifdef GLPK_PRE_4_14 |
|
6333
|
176 lib_set_print_hook (0, glpk_print_hook); |
|
|
177 #else |
|
|
178 _glp_lib_print_hook (glpk_print_hook, 0); |
|
|
179 #endif |
|
5234
|
180 |
|
|
181 LPX *lp = lpx_create_prob (); |
|
|
182 |
|
5232
|
183 |
|
5234
|
184 //-- Set the sense of optimization |
|
|
185 if (sense == 1) |
|
|
186 lpx_set_obj_dir (lp, LPX_MIN); |
|
|
187 else |
|
|
188 lpx_set_obj_dir (lp, LPX_MAX); |
|
|
189 |
|
|
190 //-- If the problem has integer structural variables switch to MIP |
|
|
191 if (isMIP) |
|
|
192 lpx_set_class (lp, LPX_MIP); |
|
5232
|
193 |
|
5234
|
194 lpx_add_cols (lp, n); |
|
|
195 for (int i = 0; i < n; i++) |
|
|
196 { |
|
5232
|
197 //-- Define type of the structural variables |
|
5234
|
198 if (! freeLB[i] && ! freeUB[i]) |
|
|
199 lpx_set_col_bnds (lp, i+1, LPX_DB, lb[i], ub[i]); |
|
|
200 else |
|
|
201 { |
|
|
202 if (! freeLB[i] && freeUB[i]) |
|
|
203 lpx_set_col_bnds (lp, i+1, LPX_LO, lb[i], ub[i]); |
|
|
204 else |
|
|
205 { |
|
|
206 if (freeLB[i] && ! freeUB[i]) |
|
|
207 lpx_set_col_bnds (lp, i+1, LPX_UP, lb[i], ub[i]); |
|
|
208 else |
|
|
209 lpx_set_col_bnds (lp, i+1, LPX_FR, lb[i], ub[i]); |
|
|
210 } |
|
|
211 } |
|
|
212 |
|
|
213 // -- Set the objective coefficient of the corresponding |
|
5232
|
214 // -- structural variable. No constant term is assumed. |
|
5234
|
215 lpx_set_obj_coef(lp,i+1,c[i]); |
|
5232
|
216 |
|
5234
|
217 if (isMIP) |
|
|
218 lpx_set_col_kind (lp, i+1, vartype[i]); |
|
|
219 } |
|
|
220 |
|
|
221 lpx_add_rows (lp, m); |
|
|
222 |
|
|
223 for (int i = 0; i < m; i++) |
|
|
224 { |
|
5232
|
225 /* If the i-th row has no lower bound (types F,U), the |
|
|
226 corrispondent parameter will be ignored. |
|
|
227 If the i-th row has no upper bound (types F,L), the corrispondent |
|
|
228 parameter will be ignored. |
|
|
229 If the i-th row is of S type, the i-th LB is used, but |
|
|
230 the i-th UB is ignored. |
|
|
231 */ |
|
5234
|
232 |
|
|
233 switch (ctype[i]) |
|
|
234 { |
|
|
235 case 'F': |
|
|
236 typx = LPX_FR; |
|
|
237 break; |
|
|
238 |
|
|
239 case 'U': |
|
|
240 typx = LPX_UP; |
|
|
241 break; |
|
|
242 |
|
|
243 case 'L': |
|
|
244 typx = LPX_LO; |
|
|
245 break; |
|
5232
|
246 |
|
5234
|
247 case 'S': |
|
|
248 typx = LPX_FX; |
|
|
249 break; |
|
5232
|
250 |
|
5234
|
251 case 'D': |
|
|
252 typx = LPX_DB; |
|
|
253 break; |
|
|
254 } |
|
|
255 |
|
|
256 lpx_set_row_bnds (lp, i+1, typx, b[i], b[i]); |
|
|
257 |
|
|
258 } |
|
|
259 |
|
|
260 lpx_load_matrix (lp, nz, rn, cn, a); |
|
5232
|
261 |
|
5234
|
262 if (save_pb) |
|
|
263 { |
|
6484
|
264 static char tmp[] = "outpb.lp"; |
|
|
265 if (lpx_write_cpxlp (lp, tmp) != 0) |
|
5234
|
266 { |
|
5240
|
267 error ("__glpk__: unable to write problem"); |
|
5234
|
268 longjmp (mark, -1); |
|
|
269 } |
|
|
270 } |
|
|
271 |
|
|
272 //-- scale the problem data (if required) |
|
|
273 //-- if (scale && (!presol || method == 1)) lpx_scale_prob(lp); |
|
|
274 //-- LPX_K_SCALE=IParam[1] LPX_K_PRESOL=IParam[16] |
|
|
275 if (lpxIntParam[1] && (! lpxIntParam[16] || lpsolver != 1)) |
|
|
276 lpx_scale_prob (lp); |
|
|
277 |
|
|
278 //-- build advanced initial basis (if required) |
|
|
279 if (lpsolver == 1 && ! lpxIntParam[16]) |
|
|
280 lpx_adv_basis (lp); |
|
|
281 |
|
|
282 for(int i = 0; i < NIntP; i++) |
|
|
283 lpx_set_int_parm (lp, IParam[i], lpxIntParam[i]); |
|
5232
|
284 |
|
5234
|
285 for (int i = 0; i < NRealP; i++) |
|
|
286 lpx_set_real_parm (lp, RParam[i], lpxRealParam[i]); |
|
|
287 |
|
|
288 if (lpsolver == 1) |
|
|
289 method = 'S'; |
|
|
290 else |
|
|
291 method = 'T'; |
|
5232
|
292 |
|
5234
|
293 switch (method) |
|
|
294 { |
|
|
295 case 'S': |
|
|
296 { |
|
|
297 if (isMIP) |
|
|
298 { |
|
|
299 method = 'I'; |
|
5240
|
300 errnum = lpx_simplex (lp); |
|
|
301 errnum = lpx_integer (lp); |
|
5234
|
302 } |
|
|
303 else |
|
5240
|
304 errnum = lpx_simplex(lp); |
|
5234
|
305 } |
|
5232
|
306 break; |
|
5234
|
307 |
|
|
308 case 'T': |
|
5240
|
309 errnum = lpx_interior(lp); |
|
5234
|
310 break; |
|
|
311 |
|
|
312 default: |
|
6381
|
313 #ifdef GLPK_PRE_4_14 |
|
5234
|
314 insist (method != method); |
|
6333
|
315 #else |
|
6484
|
316 static char tmp[] = "method != method"; |
|
|
317 glpk_fault_hook (0, tmp); |
|
6333
|
318 #endif |
|
5234
|
319 } |
|
|
320 |
|
5240
|
321 /* errnum assumes the following results: |
|
|
322 errnum = 0 <=> No errors |
|
|
323 errnum = 1 <=> Iteration limit exceeded. |
|
|
324 errnum = 2 <=> Numerical problems with basis matrix. |
|
5234
|
325 */ |
|
5240
|
326 if (errnum == LPX_E_OK) |
|
5234
|
327 { |
|
|
328 if (isMIP) |
|
|
329 { |
|
5241
|
330 *status = lpx_mip_status (lp); |
|
5234
|
331 *fmin = lpx_mip_obj_val (lp); |
|
|
332 } |
|
|
333 else |
|
|
334 { |
|
|
335 if (lpsolver == 1) |
|
|
336 { |
|
5241
|
337 *status = lpx_get_status (lp); |
|
5234
|
338 *fmin = lpx_get_obj_val (lp); |
|
|
339 } |
|
|
340 else |
|
|
341 { |
|
5241
|
342 *status = lpx_ipt_status (lp); |
|
5234
|
343 *fmin = lpx_ipt_obj_val (lp); |
|
|
344 } |
|
|
345 } |
|
5232
|
346 |
|
5234
|
347 if (isMIP) |
|
|
348 { |
|
|
349 for (int i = 0; i < n; i++) |
|
|
350 xmin[i] = lpx_mip_col_val (lp, i+1); |
|
|
351 } |
|
|
352 else |
|
|
353 { |
|
|
354 /* Primal values */ |
|
|
355 for (int i = 0; i < n; i++) |
|
|
356 { |
|
|
357 if (lpsolver == 1) |
|
|
358 xmin[i] = lpx_get_col_prim (lp, i+1); |
|
|
359 else |
|
|
360 xmin[i] = lpx_ipt_col_prim (lp, i+1); |
|
|
361 } |
|
5232
|
362 |
|
5234
|
363 /* Dual values */ |
|
|
364 for (int i = 0; i < m; i++) |
|
|
365 { |
|
|
366 if (lpsolver == 1) |
|
|
367 lambda[i] = lpx_get_row_dual (lp, i+1); |
|
|
368 else |
|
|
369 lambda[i] = lpx_ipt_row_dual (lp, i+1); |
|
|
370 } |
|
|
371 |
|
|
372 /* Reduced costs */ |
|
|
373 for (int i = 0; i < lpx_get_num_cols (lp); i++) |
|
|
374 { |
|
|
375 if (lpsolver == 1) |
|
|
376 redcosts[i] = lpx_get_col_dual (lp, i+1); |
|
|
377 else |
|
|
378 redcosts[i] = lpx_ipt_col_dual (lp, i+1); |
|
|
379 } |
|
|
380 } |
|
|
381 |
|
5241
|
382 *time = (clock () - t_start) / CLOCKS_PER_SEC; |
|
6333
|
383 |
|
6381
|
384 #ifdef GLPK_PRE_4_14 |
|
5241
|
385 *mem = (lib_env_ptr () -> mem_tpeak); |
|
6333
|
386 #else |
|
|
387 *mem = 0; |
|
|
388 #endif |
|
5234
|
389 |
|
|
390 lpx_delete_prob (lp); |
|
|
391 return 0; |
|
|
392 } |
|
|
393 |
|
|
394 lpx_delete_prob (lp); |
|
|
395 |
|
5241
|
396 *status = errnum; |
|
5234
|
397 |
|
5240
|
398 return errnum; |
|
5232
|
399 } |
|
|
400 |
|
5235
|
401 #endif |
|
|
402 |
|
5240
|
403 #define OCTAVE_GLPK_GET_REAL_PARAM(NAME, IDX) \ |
|
|
404 do \ |
|
|
405 { \ |
|
|
406 if (PARAM.contains (NAME)) \ |
|
|
407 { \ |
|
|
408 Cell tmp = PARAM.contents (NAME); \ |
|
|
409 \ |
|
|
410 if (! tmp.is_empty ()) \ |
|
|
411 { \ |
|
|
412 lpxRealParam[IDX] = tmp(0).scalar_value (); \ |
|
|
413 \ |
|
|
414 if (error_state) \ |
|
|
415 { \ |
|
|
416 error ("glpk: invalid value in param." NAME); \ |
|
|
417 return retval; \ |
|
|
418 } \ |
|
|
419 } \ |
|
|
420 else \ |
|
|
421 { \ |
|
|
422 error ("glpk: invalid value in param." NAME); \ |
|
|
423 return retval; \ |
|
|
424 } \ |
|
|
425 } \ |
|
|
426 } \ |
|
|
427 while (0) |
|
|
428 |
|
|
429 #define OCTAVE_GLPK_GET_INT_PARAM(NAME, VAL) \ |
|
|
430 do \ |
|
|
431 { \ |
|
|
432 if (PARAM.contains (NAME)) \ |
|
|
433 { \ |
|
|
434 Cell tmp = PARAM.contents (NAME); \ |
|
|
435 \ |
|
|
436 if (! tmp.is_empty ()) \ |
|
|
437 { \ |
|
|
438 VAL = tmp(0).int_value (); \ |
|
|
439 \ |
|
|
440 if (error_state) \ |
|
|
441 { \ |
|
|
442 error ("glpk: invalid value in param." NAME); \ |
|
|
443 return retval; \ |
|
|
444 } \ |
|
|
445 } \ |
|
|
446 else \ |
|
|
447 { \ |
|
|
448 error ("glpk: invalid value in param." NAME); \ |
|
|
449 return retval; \ |
|
|
450 } \ |
|
|
451 } \ |
|
|
452 } \ |
|
|
453 while (0) |
|
|
454 |
|
5235
|
455 DEFUN_DLD (__glpk__, args, , |
|
5245
|
456 "-*- texinfo -*-\n\ |
|
|
457 @deftypefn {Loadable Function} {[@var{values}] =} __glpk__ (@var{args})\n\ |
|
|
458 Internal interface for the GNU GLPK library.\n\ |
|
|
459 You should be using using the @code{glpk} function instead.\n\ |
|
|
460 @end deftypefn") |
|
5232
|
461 { |
|
5234
|
462 // The list of values to return. See the declaration in oct-obj.h |
|
|
463 octave_value_list retval; |
|
5232
|
464 |
|
5235
|
465 #if defined (HAVE_GLPK) |
|
|
466 |
|
5234
|
467 int nrhs = args.length (); |
|
|
468 |
|
5240
|
469 if (nrhs != 9) |
|
5234
|
470 { |
|
5823
|
471 print_usage (); |
|
5232
|
472 return retval; |
|
5234
|
473 } |
|
|
474 |
|
5237
|
475 //-- 1nd Input. A column array containing the objective function |
|
|
476 //-- coefficients. |
|
5880
|
477 volatile int mrowsc = args(0).rows(); |
|
5234
|
478 |
|
5237
|
479 Matrix C (args(0).matrix_value ()); |
|
5240
|
480 |
|
|
481 if (error_state) |
|
|
482 { |
|
|
483 error ("__glpk__: invalid value of C"); |
|
|
484 return retval; |
|
|
485 } |
|
|
486 |
|
5234
|
487 double *c = C.fortran_vec (); |
|
5455
|
488 Array<int> rn; |
|
|
489 Array<int> cn; |
|
|
490 ColumnVector a; |
|
5515
|
491 volatile int mrowsA; |
|
5455
|
492 volatile int nz = 0; |
|
5232
|
493 |
|
5237
|
494 //-- 2nd Input. A matrix containing the constraints coefficients. |
|
5234
|
495 // If matrix A is NOT a sparse matrix |
|
5631
|
496 if (args(1).is_sparse_type ()) |
|
5455
|
497 { |
|
5603
|
498 SparseMatrix A = args(1).sparse_matrix_value (); // get the sparse matrix |
|
5455
|
499 |
|
|
500 if (error_state) |
|
|
501 { |
|
|
502 error ("__glpk__: invalid value of A"); |
|
|
503 return retval; |
|
|
504 } |
|
5234
|
505 |
|
5455
|
506 mrowsA = A.rows (); |
|
|
507 octave_idx_type Anc = A.cols (); |
|
5604
|
508 octave_idx_type Anz = A.nzmax (); |
|
5455
|
509 rn.resize (Anz+1); |
|
|
510 cn.resize (Anz+1); |
|
|
511 a.resize (Anz+1, 0.0); |
|
|
512 |
|
|
513 if (Anc != mrowsc) |
|
|
514 { |
|
|
515 error ("__glpk__: invalid value of A"); |
|
|
516 return retval; |
|
|
517 } |
|
|
518 |
|
|
519 for (octave_idx_type j = 0; j < Anc; j++) |
|
|
520 for (octave_idx_type i = A.cidx(j); i < A.cidx(j+1); i++) |
|
|
521 { |
|
|
522 nz++; |
|
|
523 rn(nz) = A.ridx(i) + 1; |
|
|
524 cn(nz) = j + 1; |
|
5603
|
525 a(nz) = A.data(i); |
|
5455
|
526 } |
|
|
527 } |
|
5631
|
528 else |
|
|
529 { |
|
|
530 Matrix A (args(1).matrix_value ()); // get the matrix |
|
|
531 |
|
|
532 if (error_state) |
|
|
533 { |
|
|
534 error ("__glpk__: invalid value of A"); |
|
|
535 return retval; |
|
|
536 } |
|
|
537 |
|
|
538 mrowsA = A.rows (); |
|
|
539 rn.resize (mrowsA*mrowsc+1); |
|
|
540 cn.resize (mrowsA*mrowsc+1); |
|
|
541 a.resize (mrowsA*mrowsc+1, 0.0); |
|
|
542 |
|
|
543 for (int i = 0; i < mrowsA; i++) |
|
|
544 { |
|
|
545 for (int j = 0; j < mrowsc; j++) |
|
|
546 { |
|
|
547 if (A(i,j) != 0) |
|
|
548 { |
|
|
549 nz++; |
|
|
550 rn(nz) = i + 1; |
|
|
551 cn(nz) = j + 1; |
|
|
552 a(nz) = A(i,j); |
|
|
553 } |
|
|
554 } |
|
|
555 } |
|
|
556 |
|
|
557 } |
|
5232
|
558 |
|
5237
|
559 //-- 3rd Input. A column array containing the right-hand side value |
|
5234
|
560 // for each constraint in the constraint matrix. |
|
5237
|
561 Matrix B (args(2).matrix_value ()); |
|
5240
|
562 |
|
|
563 if (error_state) |
|
|
564 { |
|
|
565 error ("__glpk__: invalid value of b"); |
|
|
566 return retval; |
|
|
567 } |
|
|
568 |
|
5234
|
569 double *b = B.fortran_vec (); |
|
|
570 |
|
5237
|
571 //-- 4th Input. An array of length mrowsc containing the lower |
|
5234
|
572 //-- bound on each of the variables. |
|
5237
|
573 Matrix LB (args(3).matrix_value ()); |
|
5240
|
574 |
|
|
575 if (error_state) |
|
|
576 { |
|
|
577 error ("__glpk__: invalid value of lb"); |
|
|
578 return retval; |
|
|
579 } |
|
|
580 |
|
5234
|
581 double *lb = LB.fortran_vec (); |
|
|
582 |
|
|
583 //-- LB argument, default: Free |
|
|
584 Array<int> freeLB (mrowsc); |
|
|
585 for (int i = 0; i < mrowsc; i++) |
|
|
586 { |
|
5455
|
587 if (xisinf (lb[i])) |
|
5234
|
588 { |
|
|
589 freeLB(i) = 1; |
|
|
590 lb[i] = -octave_Inf; |
|
|
591 } |
|
|
592 else |
|
|
593 freeLB(i) = 0; |
|
|
594 } |
|
|
595 |
|
5237
|
596 //-- 5th Input. An array of at least length numcols containing the upper |
|
5234
|
597 //-- bound on each of the variables. |
|
5237
|
598 Matrix UB (args(4).matrix_value ()); |
|
5234
|
599 |
|
5240
|
600 if (error_state) |
|
|
601 { |
|
|
602 error ("__glpk__: invalid value of ub"); |
|
|
603 return retval; |
|
|
604 } |
|
|
605 |
|
5234
|
606 double *ub = UB.fortran_vec (); |
|
5232
|
607 |
|
5234
|
608 Array<int> freeUB (mrowsc); |
|
|
609 for (int i = 0; i < mrowsc; i++) |
|
|
610 { |
|
5455
|
611 if (xisinf (ub[i])) |
|
5234
|
612 { |
|
|
613 freeUB(i) = 1; |
|
|
614 ub[i] = octave_Inf; |
|
|
615 } |
|
|
616 else |
|
|
617 freeUB(i) = 0; |
|
|
618 } |
|
|
619 |
|
5237
|
620 //-- 6th Input. A column array containing the sense of each constraint |
|
|
621 //-- in the constraint matrix. |
|
|
622 charMatrix CTYPE (args(5).char_matrix_value ()); |
|
5240
|
623 |
|
|
624 if (error_state) |
|
|
625 { |
|
|
626 error ("__glpk__: invalid value of ctype"); |
|
|
627 return retval; |
|
|
628 } |
|
|
629 |
|
5237
|
630 char *ctype = CTYPE.fortran_vec (); |
|
|
631 |
|
|
632 //-- 7th Input. A column array containing the types of the variables. |
|
|
633 charMatrix VTYPE (args(6).char_matrix_value ()); |
|
5232
|
634 |
|
5240
|
635 if (error_state) |
|
|
636 { |
|
|
637 error ("__glpk__: invalid value of vtype"); |
|
|
638 return retval; |
|
|
639 } |
|
|
640 |
|
5234
|
641 Array<int> vartype (mrowsc); |
|
5235
|
642 volatile int isMIP = 0; |
|
5234
|
643 for (int i = 0; i < mrowsc ; i++) |
|
|
644 { |
|
|
645 if (VTYPE(i,0) == 'I') |
|
|
646 { |
|
|
647 isMIP = 1; |
|
|
648 vartype(i) = LPX_IV; |
|
|
649 } |
|
|
650 else |
|
5235
|
651 vartype(i) = LPX_CV; |
|
5234
|
652 } |
|
|
653 |
|
5237
|
654 //-- 8th Input. Sense of optimization. |
|
|
655 volatile int sense; |
|
|
656 double SENSE = args(7).scalar_value (); |
|
5240
|
657 |
|
|
658 if (error_state) |
|
|
659 { |
|
|
660 error ("__glpk__: invalid value of sense"); |
|
|
661 return retval; |
|
|
662 } |
|
|
663 |
|
5237
|
664 if (SENSE >= 0) |
|
|
665 sense = 1; |
|
|
666 else |
|
|
667 sense = -1; |
|
|
668 |
|
5234
|
669 //-- 9th Input. A structure containing the control parameters. |
|
|
670 Octave_map PARAM = args(8).map_value (); |
|
|
671 |
|
5240
|
672 if (error_state) |
|
|
673 { |
|
|
674 error ("__glpk__: invalid value of param"); |
|
|
675 return retval; |
|
|
676 } |
|
|
677 |
|
5234
|
678 //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
|
679 //-- Integer parameters |
|
|
680 //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
|
681 |
|
|
682 //-- Level of messages output by the solver |
|
5240
|
683 OCTAVE_GLPK_GET_INT_PARAM ("msglev", lpxIntParam[0]); |
|
|
684 if (lpxIntParam[0] < 0 || lpxIntParam[0] > 3) |
|
5234
|
685 { |
|
5240
|
686 error ("__glpk__: param.msglev must be 0 (no output [default]) or 1 (error messages only) or 2 (normal output) or 3 (full output)"); |
|
|
687 return retval; |
|
|
688 } |
|
5232
|
689 |
|
5234
|
690 //-- scaling option |
|
5240
|
691 OCTAVE_GLPK_GET_INT_PARAM ("scale", lpxIntParam[1]); |
|
|
692 if (lpxIntParam[1] < 0 || lpxIntParam[1] > 2) |
|
5234
|
693 { |
|
5240
|
694 error ("__glpk__: param.scale must be 0 (no scaling) or 1 (equilibration scaling [default]) or 2 (geometric mean scaling)"); |
|
|
695 return retval; |
|
5234
|
696 } |
|
5232
|
697 |
|
5234
|
698 //-- Dual dimplex option |
|
5240
|
699 OCTAVE_GLPK_GET_INT_PARAM ("dual", lpxIntParam[2]); |
|
|
700 if (lpxIntParam[2] < 0 || lpxIntParam[2] > 1) |
|
5234
|
701 { |
|
5240
|
702 error ("__glpk__: param.dual must be 0 (do NOT use dual simplex [default]) or 1 (use dual simplex)"); |
|
|
703 return retval; |
|
5234
|
704 } |
|
5232
|
705 |
|
5234
|
706 //-- Pricing option |
|
5240
|
707 OCTAVE_GLPK_GET_INT_PARAM ("price", lpxIntParam[3]); |
|
|
708 if (lpxIntParam[3] < 0 || lpxIntParam[3] > 1) |
|
5234
|
709 { |
|
5240
|
710 error ("__glpk__: param.price must be 0 (textbook pricing) or 1 (steepest edge pricing [default])"); |
|
|
711 return retval; |
|
|
712 } |
|
5232
|
713 |
|
5234
|
714 //-- Solution rounding option |
|
5240
|
715 OCTAVE_GLPK_GET_INT_PARAM ("round", lpxIntParam[4]); |
|
|
716 if (lpxIntParam[4] < 0 || lpxIntParam[4] > 1) |
|
5234
|
717 { |
|
5240
|
718 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)"); |
|
|
719 return retval; |
|
5234
|
720 } |
|
|
721 |
|
|
722 //-- Simplex iterations limit |
|
5240
|
723 OCTAVE_GLPK_GET_INT_PARAM ("itlim", lpxIntParam[5]); |
|
5232
|
724 |
|
5234
|
725 //-- Simplex iterations count |
|
5240
|
726 OCTAVE_GLPK_GET_INT_PARAM ("itcnt", lpxIntParam[6]); |
|
5234
|
727 |
|
|
728 //-- Output frequency, in iterations |
|
5240
|
729 OCTAVE_GLPK_GET_INT_PARAM ("outfrq", lpxIntParam[7]); |
|
5234
|
730 |
|
|
731 //-- Branching heuristic option |
|
5240
|
732 OCTAVE_GLPK_GET_INT_PARAM ("branch", lpxIntParam[14]); |
|
|
733 if (lpxIntParam[14] < 0 || lpxIntParam[14] > 2) |
|
5234
|
734 { |
|
5240
|
735 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]"); |
|
|
736 return retval; |
|
|
737 } |
|
5232
|
738 |
|
5234
|
739 //-- Backtracking heuristic option |
|
5240
|
740 OCTAVE_GLPK_GET_INT_PARAM ("btrack", lpxIntParam[15]); |
|
|
741 if (lpxIntParam[15] < 0 || lpxIntParam[15] > 2) |
|
5234
|
742 { |
|
5240
|
743 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]"); |
|
|
744 return retval; |
|
|
745 } |
|
5232
|
746 |
|
5234
|
747 //-- Presolver option |
|
5240
|
748 OCTAVE_GLPK_GET_INT_PARAM ("presol", lpxIntParam[16]); |
|
|
749 if (lpxIntParam[16] < 0 || lpxIntParam[16] > 1) |
|
5234
|
750 { |
|
5240
|
751 error ("__glpk__: param.presol must be 0 (do NOT use LP presolver) or 1 (use LP presolver [default])"); |
|
|
752 return retval; |
|
5234
|
753 } |
|
|
754 |
|
5237
|
755 //-- LPsolver option |
|
|
756 volatile int lpsolver = 1; |
|
5240
|
757 OCTAVE_GLPK_GET_INT_PARAM ("lpsolver", lpsolver); |
|
|
758 if (lpsolver < 1 || lpsolver > 2) |
|
5237
|
759 { |
|
5240
|
760 error ("__glpk__: param.lpsolver must be 1 (simplex method) or 2 (interior point method)"); |
|
|
761 return retval; |
|
5237
|
762 } |
|
|
763 |
|
|
764 //-- Save option |
|
|
765 volatile int save_pb = 0; |
|
5240
|
766 OCTAVE_GLPK_GET_INT_PARAM ("save", save_pb); |
|
|
767 save_pb = save_pb != 0; |
|
5237
|
768 |
|
5234
|
769 //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
|
770 //-- Real parameters |
|
|
771 //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
|
772 |
|
|
773 //-- Ratio test option |
|
5240
|
774 OCTAVE_GLPK_GET_REAL_PARAM ("relax", 0); |
|
5232
|
775 |
|
5234
|
776 //-- Relative tolerance used to check if the current basic solution |
|
|
777 //-- is primal feasible |
|
5240
|
778 OCTAVE_GLPK_GET_REAL_PARAM ("tolbnd", 1); |
|
5234
|
779 |
|
|
780 //-- Absolute tolerance used to check if the current basic solution |
|
|
781 //-- is dual feasible |
|
5240
|
782 OCTAVE_GLPK_GET_REAL_PARAM ("toldj", 2); |
|
5232
|
783 |
|
5234
|
784 //-- Relative tolerance used to choose eligible pivotal elements of |
|
|
785 //-- the simplex table in the ratio test |
|
5240
|
786 OCTAVE_GLPK_GET_REAL_PARAM ("tolpiv", 3); |
|
5234
|
787 |
|
5240
|
788 OCTAVE_GLPK_GET_REAL_PARAM ("objll", 4); |
|
5234
|
789 |
|
5240
|
790 OCTAVE_GLPK_GET_REAL_PARAM ("objul", 5); |
|
5232
|
791 |
|
5240
|
792 OCTAVE_GLPK_GET_REAL_PARAM ("tmlim", 6); |
|
5234
|
793 |
|
5240
|
794 OCTAVE_GLPK_GET_REAL_PARAM ("outdly", 7); |
|
5234
|
795 |
|
5240
|
796 OCTAVE_GLPK_GET_REAL_PARAM ("tolint", 8); |
|
5234
|
797 |
|
5240
|
798 OCTAVE_GLPK_GET_REAL_PARAM ("tolobj", 9); |
|
5234
|
799 |
|
|
800 //-- Assign pointers to the output parameters |
|
|
801 ColumnVector xmin (mrowsc); |
|
|
802 ColumnVector fmin (1); |
|
|
803 ColumnVector status (1); |
|
|
804 ColumnVector lambda (mrowsA); |
|
|
805 ColumnVector redcosts (mrowsc); |
|
|
806 ColumnVector time (1); |
|
|
807 ColumnVector mem (1); |
|
5232
|
808 |
|
5234
|
809 int jmpret = setjmp (mark); |
|
5235
|
810 |
|
5234
|
811 if (jmpret == 0) |
|
5235
|
812 glpk (sense, mrowsc, mrowsA, c, nz, rn.fortran_vec (), |
|
|
813 cn.fortran_vec (), a.fortran_vec (), b, ctype, |
|
|
814 freeLB.fortran_vec (), lb, freeUB.fortran_vec (), |
|
|
815 ub, vartype.fortran_vec (), isMIP, lpsolver, |
|
|
816 save_pb, xmin.fortran_vec (), fmin.fortran_vec (), |
|
|
817 status.fortran_vec (), lambda.fortran_vec (), |
|
|
818 redcosts.fortran_vec (), time.fortran_vec (), |
|
|
819 mem.fortran_vec ()); |
|
5232
|
820 |
|
5234
|
821 Octave_map extra; |
|
|
822 |
|
5238
|
823 if (! isMIP) |
|
|
824 { |
|
|
825 extra.assign ("lambda", octave_value (lambda)); |
|
|
826 extra.assign ("redcosts", octave_value (redcosts)); |
|
|
827 } |
|
|
828 |
|
5234
|
829 extra.assign ("time", octave_value (time)); |
|
|
830 extra.assign ("mem", octave_value (mem)); |
|
5232
|
831 |
|
5234
|
832 retval(3) = extra; |
|
|
833 retval(2) = octave_value(status); |
|
|
834 retval(1) = octave_value(fmin); |
|
|
835 retval(0) = octave_value(xmin); |
|
|
836 |
|
5235
|
837 #else |
|
|
838 |
|
|
839 gripe_not_supported ("glpk"); |
|
|
840 |
|
|
841 #endif |
|
|
842 |
|
5234
|
843 return retval; |
|
5232
|
844 } |
