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