|
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 |
|
|
19 Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. |
|
|
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 |
|
5234
|
40 extern "C" { |
|
|
41 #include <glpk.h> |
|
5232
|
42 } |
|
|
43 |
|
|
44 #define OCTOUT octave_stdout |
|
|
45 #define OCTERR octave_stdout |
|
|
46 #define NIntP 17 |
|
|
47 #define NRealP 10 |
|
|
48 |
|
5234
|
49 int lpxIntParam[NIntP] = { |
|
|
50 1, |
|
|
51 1, |
|
|
52 0, |
|
|
53 1, |
|
|
54 0, |
|
|
55 -1, |
|
|
56 0, |
|
|
57 200, |
|
|
58 1, |
|
|
59 2, |
|
|
60 0, |
|
|
61 1, |
|
|
62 0, |
|
|
63 0, |
|
|
64 2, |
|
|
65 2, |
|
|
66 1 |
|
5232
|
67 }; |
|
|
68 |
|
5234
|
69 int IParam[NIntP] = { |
|
|
70 LPX_K_MSGLEV, |
|
|
71 LPX_K_SCALE, |
|
|
72 LPX_K_DUAL, |
|
|
73 LPX_K_PRICE, |
|
|
74 LPX_K_ROUND, |
|
|
75 LPX_K_ITLIM, |
|
|
76 LPX_K_ITCNT, |
|
|
77 LPX_K_OUTFRQ, |
|
|
78 LPX_K_MPSINFO, |
|
|
79 LPX_K_MPSOBJ, |
|
|
80 LPX_K_MPSORIG, |
|
|
81 LPX_K_MPSWIDE, |
|
|
82 LPX_K_MPSFREE, |
|
|
83 LPX_K_MPSSKIP, |
|
|
84 LPX_K_BRANCH, |
|
|
85 LPX_K_BTRACK, |
|
|
86 LPX_K_PRESOL |
|
5232
|
87 }; |
|
|
88 |
|
|
89 |
|
5234
|
90 double lpxRealParam[NRealP] = { |
|
|
91 0.07, |
|
|
92 1e-7, |
|
|
93 1e-7, |
|
|
94 1e-9, |
|
|
95 -DBL_MAX, |
|
|
96 DBL_MAX, |
|
|
97 -1.0, |
|
|
98 0.0, |
|
|
99 1e-6, |
|
|
100 1e-7 |
|
5232
|
101 }; |
|
|
102 |
|
5234
|
103 int RParam[NRealP] = { |
|
|
104 LPX_K_RELAX, |
|
|
105 LPX_K_TOLBND, |
|
|
106 LPX_K_TOLDJ, |
|
|
107 LPX_K_TOLPIV, |
|
|
108 LPX_K_OBJLL, |
|
|
109 LPX_K_OBJUL, |
|
|
110 LPX_K_TMLIM, |
|
|
111 LPX_K_OUTDLY, |
|
|
112 LPX_K_TOLINT, |
|
|
113 LPX_K_TOLOBJ |
|
5232
|
114 }; |
|
|
115 |
|
|
116 jmp_buf mark; //-- Address for long jump to jump to |
|
5234
|
117 int fperr; //-- Global error number |
|
5232
|
118 |
|
|
119 |
|
5234
|
120 int |
|
5235
|
121 glpk_fault_hook (void * /* info */, char *msg) |
|
5234
|
122 { |
|
|
123 OCTERR << "*** SEVERE CRITICAL ERROR *** from GLPK !\n\n"<<msg<<" %s\n"; |
|
|
124 longjmp (mark, -1); |
|
5232
|
125 } |
|
|
126 |
|
5234
|
127 int |
|
5235
|
128 glpk_print_hook (void * /* info */, char *msg) |
|
5232
|
129 { |
|
5234
|
130 OCTERR << msg << "\n"; |
|
|
131 return 1; |
|
5232
|
132 } |
|
|
133 |
|
|
134 |
|
5234
|
135 int |
|
|
136 glpk (int sense, int n, int m, double *c, int nz, int *rn, int *cn, |
|
|
137 double *a, double *b, char *ctype, int *freeLB, double *lb, |
|
|
138 int *freeUB, double *ub, int *vartype, int isMIP, int lpsolver, |
|
|
139 int save_pb, double *xmin, double *fmin, double *status, |
|
|
140 double *lambda, double *redcosts, double *time, double *mem) |
|
5232
|
141 { |
|
5234
|
142 int error; |
|
|
143 int typx = 0; |
|
|
144 int method; |
|
|
145 |
|
|
146 clock_t t_start = clock(); |
|
|
147 |
|
|
148 lib_set_fault_hook (NULL, glpk_fault_hook); |
|
5232
|
149 |
|
5234
|
150 if (lpxIntParam[0] > 1) |
|
|
151 lib_set_print_hook (NULL, glpk_print_hook); |
|
|
152 |
|
|
153 LPX *lp = lpx_create_prob (); |
|
|
154 |
|
5232
|
155 |
|
5234
|
156 //-- Set the sense of optimization |
|
|
157 if (sense == 1) |
|
|
158 lpx_set_obj_dir (lp, LPX_MIN); |
|
|
159 else |
|
|
160 lpx_set_obj_dir (lp, LPX_MAX); |
|
|
161 |
|
|
162 //-- If the problem has integer structural variables switch to MIP |
|
|
163 if (isMIP) |
|
|
164 lpx_set_class (lp, LPX_MIP); |
|
5232
|
165 |
|
5234
|
166 lpx_add_cols (lp, n); |
|
|
167 for (int i = 0; i < n; i++) |
|
|
168 { |
|
5232
|
169 //-- Define type of the structural variables |
|
5234
|
170 if (! freeLB[i] && ! freeUB[i]) |
|
|
171 lpx_set_col_bnds (lp, i+1, LPX_DB, lb[i], ub[i]); |
|
|
172 else |
|
|
173 { |
|
|
174 if (! freeLB[i] && freeUB[i]) |
|
|
175 lpx_set_col_bnds (lp, i+1, LPX_LO, lb[i], ub[i]); |
|
|
176 else |
|
|
177 { |
|
|
178 if (freeLB[i] && ! freeUB[i]) |
|
|
179 lpx_set_col_bnds (lp, i+1, LPX_UP, lb[i], ub[i]); |
|
|
180 else |
|
|
181 lpx_set_col_bnds (lp, i+1, LPX_FR, lb[i], ub[i]); |
|
|
182 } |
|
|
183 } |
|
|
184 |
|
|
185 // -- Set the objective coefficient of the corresponding |
|
5232
|
186 // -- structural variable. No constant term is assumed. |
|
5234
|
187 lpx_set_obj_coef(lp,i+1,c[i]); |
|
5232
|
188 |
|
5234
|
189 if (isMIP) |
|
|
190 lpx_set_col_kind (lp, i+1, vartype[i]); |
|
|
191 } |
|
|
192 |
|
|
193 lpx_add_rows (lp, m); |
|
|
194 |
|
|
195 for (int i = 0; i < m; i++) |
|
|
196 { |
|
5232
|
197 /* If the i-th row has no lower bound (types F,U), the |
|
|
198 corrispondent parameter will be ignored. |
|
|
199 If the i-th row has no upper bound (types F,L), the corrispondent |
|
|
200 parameter will be ignored. |
|
|
201 If the i-th row is of S type, the i-th LB is used, but |
|
|
202 the i-th UB is ignored. |
|
|
203 */ |
|
5234
|
204 |
|
|
205 switch (ctype[i]) |
|
|
206 { |
|
|
207 case 'F': |
|
|
208 typx = LPX_FR; |
|
|
209 break; |
|
|
210 |
|
|
211 case 'U': |
|
|
212 typx = LPX_UP; |
|
|
213 break; |
|
|
214 |
|
|
215 case 'L': |
|
|
216 typx = LPX_LO; |
|
|
217 break; |
|
5232
|
218 |
|
5234
|
219 case 'S': |
|
|
220 typx = LPX_FX; |
|
|
221 break; |
|
5232
|
222 |
|
5234
|
223 case 'D': |
|
|
224 typx = LPX_DB; |
|
|
225 break; |
|
|
226 } |
|
|
227 |
|
|
228 lpx_set_row_bnds (lp, i+1, typx, b[i], b[i]); |
|
|
229 |
|
|
230 } |
|
|
231 |
|
|
232 lpx_load_matrix (lp, nz, rn, cn, a); |
|
5232
|
233 |
|
5234
|
234 if (save_pb) |
|
|
235 { |
|
|
236 if (lpx_write_cpxlp (lp, "outpb.lp") != 0) |
|
|
237 { |
|
|
238 OCTERR << "Unable to write problem\n"; |
|
|
239 longjmp (mark, -1); |
|
|
240 } |
|
|
241 } |
|
|
242 |
|
|
243 //-- scale the problem data (if required) |
|
|
244 //-- if (scale && (!presol || method == 1)) lpx_scale_prob(lp); |
|
|
245 //-- LPX_K_SCALE=IParam[1] LPX_K_PRESOL=IParam[16] |
|
|
246 if (lpxIntParam[1] && (! lpxIntParam[16] || lpsolver != 1)) |
|
|
247 lpx_scale_prob (lp); |
|
|
248 |
|
|
249 //-- build advanced initial basis (if required) |
|
|
250 if (lpsolver == 1 && ! lpxIntParam[16]) |
|
|
251 lpx_adv_basis (lp); |
|
|
252 |
|
|
253 for(int i = 0; i < NIntP; i++) |
|
|
254 lpx_set_int_parm (lp, IParam[i], lpxIntParam[i]); |
|
5232
|
255 |
|
5234
|
256 for (int i = 0; i < NRealP; i++) |
|
|
257 lpx_set_real_parm (lp, RParam[i], lpxRealParam[i]); |
|
|
258 |
|
|
259 if (lpsolver == 1) |
|
|
260 method = 'S'; |
|
|
261 else |
|
|
262 method = 'T'; |
|
5232
|
263 |
|
5234
|
264 switch (method) |
|
|
265 { |
|
|
266 case 'S': |
|
|
267 { |
|
|
268 if (isMIP) |
|
|
269 { |
|
|
270 method = 'I'; |
|
|
271 error = lpx_simplex (lp); |
|
|
272 error = lpx_integer (lp); |
|
|
273 } |
|
|
274 else |
|
5235
|
275 error = lpx_simplex(lp); |
|
5234
|
276 } |
|
5232
|
277 break; |
|
5234
|
278 |
|
|
279 case 'T': |
|
|
280 error = lpx_interior(lp); |
|
|
281 break; |
|
|
282 |
|
|
283 default: |
|
|
284 insist (method != method); |
|
|
285 } |
|
|
286 |
|
|
287 /* error assumes the following results: |
|
|
288 error=0 <=> No errors |
|
|
289 error=1 <=> Iteration limit exceeded. |
|
|
290 error=2 <=> Numerical problems with basis matrix. |
|
|
291 */ |
|
|
292 if (error == LPX_E_OK) |
|
|
293 { |
|
|
294 if (isMIP) |
|
|
295 { |
|
|
296 *status = static_cast<double> (lpx_mip_status (lp)); |
|
|
297 *fmin = lpx_mip_obj_val (lp); |
|
|
298 } |
|
|
299 else |
|
|
300 { |
|
|
301 if (lpsolver == 1) |
|
|
302 { |
|
|
303 *status = static_cast<double> (lpx_get_status (lp)); |
|
|
304 *fmin = lpx_get_obj_val (lp); |
|
|
305 } |
|
|
306 else |
|
|
307 { |
|
|
308 *status = static_cast<double> (lpx_ipt_status (lp)); |
|
|
309 *fmin = lpx_ipt_obj_val (lp); |
|
|
310 } |
|
|
311 } |
|
5232
|
312 |
|
5234
|
313 if (isMIP) |
|
|
314 { |
|
|
315 for (int i = 0; i < n; i++) |
|
|
316 xmin[i] = lpx_mip_col_val (lp, i+1); |
|
|
317 } |
|
|
318 else |
|
|
319 { |
|
|
320 /* Primal values */ |
|
|
321 for (int i = 0; i < n; i++) |
|
|
322 { |
|
|
323 if (lpsolver == 1) |
|
|
324 xmin[i] = lpx_get_col_prim (lp, i+1); |
|
|
325 else |
|
|
326 xmin[i] = lpx_ipt_col_prim (lp, i+1); |
|
|
327 } |
|
5232
|
328 |
|
5234
|
329 /* Dual values */ |
|
|
330 for (int i = 0; i < m; i++) |
|
|
331 { |
|
|
332 if (lpsolver == 1) |
|
|
333 lambda[i] = lpx_get_row_dual (lp, i+1); |
|
|
334 else |
|
|
335 lambda[i] = lpx_ipt_row_dual (lp, i+1); |
|
|
336 } |
|
|
337 |
|
|
338 /* Reduced costs */ |
|
|
339 for (int i = 0; i < lpx_get_num_cols (lp); i++) |
|
|
340 { |
|
|
341 if (lpsolver == 1) |
|
|
342 redcosts[i] = lpx_get_col_dual (lp, i+1); |
|
|
343 else |
|
|
344 redcosts[i] = lpx_ipt_col_dual (lp, i+1); |
|
|
345 } |
|
|
346 } |
|
|
347 |
|
|
348 *time = static_cast<double> (clock () - t_start) / CLOCKS_PER_SEC; |
|
|
349 *mem = static_cast<double> (lib_env_ptr () -> mem_tpeak); |
|
|
350 |
|
|
351 lpx_delete_prob (lp); |
|
|
352 return 0; |
|
|
353 } |
|
|
354 |
|
|
355 lpx_delete_prob (lp); |
|
|
356 |
|
|
357 *status= static_cast<double> (error); |
|
|
358 |
|
|
359 return error; |
|
5232
|
360 } |
|
|
361 |
|
5235
|
362 #endif |
|
|
363 |
|
|
364 DEFUN_DLD (__glpk__, args, , |
|
5234
|
365 "__glpk__: internal interface for the GLPK library.\n\ |
|
|
366 You should be using using glpk instead") |
|
5232
|
367 { |
|
5234
|
368 // The list of values to return. See the declaration in oct-obj.h |
|
|
369 octave_value_list retval; |
|
5232
|
370 |
|
5235
|
371 #if defined (HAVE_GLPK) |
|
|
372 |
|
5234
|
373 int nrhs = args.length (); |
|
|
374 |
|
|
375 if (nrhs < 1) |
|
|
376 { |
|
5232
|
377 OCTERR<<"Use the script glpk for the optimization\n"; |
|
|
378 return retval; |
|
5234
|
379 } |
|
|
380 |
|
|
381 //-- 1st Input. Sense of optimization. |
|
5235
|
382 volatile int sense; |
|
5234
|
383 double SENSE = args(0).scalar_value (); |
|
5235
|
384 if (SENSE >= 0) |
|
5234
|
385 sense = 1; |
|
|
386 else |
|
|
387 sense = -1; |
|
|
388 |
|
|
389 //-- 2nd Input. A column array containing the objective function |
|
|
390 //-- coefficients. |
|
|
391 int mrowsc = args(1).rows(); |
|
|
392 |
|
|
393 Matrix C (args(1).matrix_value ()); |
|
|
394 double *c = C.fortran_vec (); |
|
5232
|
395 |
|
5234
|
396 //-- 3rd Input. A matrix containing the constraints coefficients. |
|
|
397 // If matrix A is NOT a sparse matrix |
|
|
398 // if(!mxIsSparse(A_IN)){ |
|
|
399 int mrowsA = args(2).rows(); |
|
|
400 Matrix A (args(2).matrix_value ()); // get the matrix |
|
|
401 Array<int> rn (mrowsA*mrowsc+1); |
|
|
402 Array<int> cn (mrowsA*mrowsc+1); |
|
|
403 ColumnVector a (mrowsA*mrowsc+1, 0.0); |
|
|
404 |
|
5235
|
405 volatile int nz = 0; |
|
5234
|
406 for (int i = 0; i < mrowsA; i++) |
|
|
407 { |
|
|
408 for (int j = 0; j < mrowsc; j++) |
|
|
409 { |
|
|
410 if (A(i,j) != 0) |
|
|
411 { |
|
|
412 nz++; |
|
|
413 rn(nz) = i + 1; |
|
|
414 cn(nz) = j + 1; |
|
|
415 a(nz) = A(i,j); |
|
|
416 } |
|
|
417 } |
|
|
418 } |
|
|
419 |
|
5232
|
420 // DON'T DELETE THIS PART... REPRESENTS THE SPARSE MATRICES MANIPULATION |
|
|
421 // }else{ |
|
|
422 // int i,j; |
|
|
423 // int *jc,*ir; |
|
|
424 // double *pr; |
|
|
425 // int nelc,count,row; |
|
|
426 // |
|
|
427 // /* NOTE: nnz is the actual number of nonzeros and is stored as the |
|
|
428 // last element of the jc array where the size of the jc array is the |
|
|
429 // number of columns + 1 */ |
|
|
430 // nz = *(mxGetJc(A_IN) + mrowsc); |
|
|
431 // jc = mxGetJc(A_IN); |
|
|
432 // ir = mxGetIr(A_IN); |
|
|
433 // pr = mxGetPr(A_IN); |
|
|
434 // |
|
|
435 // rn=(int *)calloc(nz+1,sizeof(int)); |
|
|
436 // cn=(int *)calloc(nz+1,sizeof(int)); |
|
|
437 // a=(double *)calloc(nz+1,sizeof(double)); |
|
|
438 // |
|
|
439 // count=0; row=0; |
|
|
440 // for(i=1;i<=mrowsc;i++){ |
|
|
441 // nelc=jc[i]-jc[i-1]; |
|
|
442 // for(j=0;j<nelc;j++){ |
|
|
443 // count++; |
|
|
444 // rn[count]=ir[row]+1; |
|
|
445 // cn[count]=i; |
|
|
446 // a[count]=pr[row]; |
|
|
447 // row++; |
|
|
448 // } |
|
|
449 // } |
|
|
450 // } |
|
|
451 |
|
5234
|
452 //-- 4th Input. A column array containing the right-hand side value |
|
|
453 // for each constraint in the constraint matrix. |
|
|
454 Matrix B (args(3).matrix_value ()); |
|
|
455 double *b = B.fortran_vec (); |
|
|
456 |
|
|
457 //-- 5th Input. A column array containing the sense of each constraint |
|
|
458 //-- in the constraint matrix. |
|
|
459 charMatrix CTYPE (args(4).char_matrix_value ()); |
|
|
460 char *ctype = CTYPE.fortran_vec (); |
|
|
461 |
|
|
462 //-- 6th Input. An array of length mrowsc containing the lower |
|
|
463 //-- bound on each of the variables. |
|
|
464 Matrix LB (args(5).matrix_value ()); |
|
|
465 double *lb = LB.fortran_vec (); |
|
|
466 |
|
|
467 //-- LB argument, default: Free |
|
|
468 Array<int> freeLB (mrowsc); |
|
|
469 for (int i = 0; i < mrowsc; i++) |
|
|
470 { |
|
|
471 if (isinf (lb[i])) |
|
|
472 { |
|
|
473 freeLB(i) = 1; |
|
|
474 lb[i] = -octave_Inf; |
|
|
475 } |
|
|
476 else |
|
|
477 freeLB(i) = 0; |
|
|
478 } |
|
|
479 |
|
|
480 //-- 7th Input. An array of at least length numcols containing the upper |
|
|
481 //-- bound on each of the variables. |
|
|
482 Matrix UB (args(6).matrix_value ()); |
|
|
483 |
|
|
484 double *ub = UB.fortran_vec (); |
|
5232
|
485 |
|
5234
|
486 Array<int> freeUB (mrowsc); |
|
|
487 for (int i = 0; i < mrowsc; i++) |
|
|
488 { |
|
|
489 if (isinf (ub[i])) |
|
|
490 { |
|
|
491 freeUB(i) = 1; |
|
|
492 ub[i] = octave_Inf; |
|
|
493 } |
|
|
494 else |
|
|
495 freeUB(i) = 0; |
|
|
496 } |
|
|
497 |
|
|
498 //-- 8th Input. A column array containing the types of the variables. |
|
|
499 charMatrix VTYPE (args(7).char_matrix_value ()); |
|
5232
|
500 |
|
5234
|
501 Array<int> vartype (mrowsc); |
|
5235
|
502 volatile int isMIP = 0; |
|
5234
|
503 for (int i = 0; i < mrowsc ; i++) |
|
|
504 { |
|
|
505 if (VTYPE(i,0) == 'I') |
|
|
506 { |
|
|
507 isMIP = 1; |
|
|
508 vartype(i) = LPX_IV; |
|
|
509 } |
|
|
510 else |
|
5235
|
511 vartype(i) = LPX_CV; |
|
5234
|
512 } |
|
|
513 |
|
|
514 //-- 9th Input. A structure containing the control parameters. |
|
|
515 Octave_map PARAM = args(8).map_value (); |
|
|
516 |
|
|
517 //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
|
518 //-- Integer parameters |
|
|
519 //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
|
520 |
|
|
521 //-- Level of messages output by the solver |
|
|
522 if (PARAM.contains ("msglev")) |
|
|
523 { |
|
|
524 octave_value tmp = PARAM.contents (PARAM.seek ("msglev"))(0); |
|
|
525 |
|
|
526 double numtmp = tmp.scalar_value (); |
|
|
527 if (numtmp != 0 && numtmp != 1 && numtmp != 2 && numtmp != 3) |
|
|
528 { |
|
|
529 OCTOUT << "'msglev' parameter must be only:\n\t0 - no output,\n\t1 - error messages only),\n\t2 - normal output,\n\t3 - full output [default]\n"; |
|
|
530 return retval; |
|
|
531 } |
|
|
532 |
|
|
533 lpxIntParam[0] = static_cast<int> (numtmp); |
|
5232
|
534 } |
|
|
535 |
|
5234
|
536 //-- scaling option |
|
|
537 if (PARAM.contains ("scale")) |
|
|
538 { |
|
|
539 octave_value tmp = PARAM.contents (PARAM.seek ("scale"))(0); |
|
|
540 double numtmp = tmp.scalar_value (); |
|
|
541 if (numtmp != 0 && numtmp != 1 && numtmp != 2) |
|
|
542 { |
|
|
543 OCTOUT << "'scale' parameter must be only:\n\t0 - no scaling,\n\t1 - equilibration scaling,\n\t2 - geometric mean scaling\n"; |
|
|
544 return retval; |
|
|
545 } |
|
|
546 lpxIntParam[1] = static_cast<int> (numtmp); |
|
|
547 } |
|
5232
|
548 |
|
5234
|
549 //-- Dual dimplex option |
|
|
550 if (PARAM.contains ("dual")) |
|
|
551 { |
|
|
552 octave_value tmp = PARAM.contents (PARAM.seek ("dual"))(0); |
|
|
553 double numtmp = tmp.scalar_value (); |
|
|
554 if (numtmp != 0 && numtmp != 1) |
|
|
555 { |
|
|
556 OCTOUT<<"'dual' parameter must be only:\n\t0 - do not use the dual simplex [default],\n\t1 - use dual simplex\n"; |
|
|
557 return retval; |
|
|
558 } |
|
|
559 lpxIntParam[2] = static_cast<int> (numtmp); |
|
|
560 } |
|
5232
|
561 |
|
5234
|
562 //-- Pricing option |
|
|
563 if (PARAM.contains ("price")) |
|
|
564 { |
|
|
565 octave_value tmp = PARAM.contents (PARAM.seek ("price"))(0); |
|
|
566 double numtmp = tmp.scalar_value(); |
|
|
567 if (numtmp != 0 && numtmp != 1) |
|
|
568 { |
|
|
569 OCTOUT << "'price' parameter must be only:\n\t0 - textbook pricing,\n\t1 - steepest edge pricing [default]\n"; |
|
|
570 return retval; |
|
|
571 } |
|
|
572 lpxIntParam[3] = static_cast<int> (numtmp); |
|
5232
|
573 } |
|
|
574 |
|
5234
|
575 //-- Solution rounding option |
|
|
576 if (PARAM.contains ("round")) |
|
|
577 { |
|
|
578 octave_value tmp = PARAM.contents (PARAM.seek ("round"))(0); |
|
|
579 double numtmp = tmp.scalar_value (); |
|
|
580 if (numtmp != 0 && numtmp != 1) |
|
|
581 { |
|
|
582 OCTOUT << "'round' parameter must be only:\n\t0 - report all primal and dual values [default],\n\t1 - replace tiny primal and dual values by exact zero\n"; |
|
|
583 return retval; |
|
|
584 } |
|
|
585 lpxIntParam[4] = static_cast<int> (numtmp); |
|
|
586 } |
|
|
587 |
|
|
588 //-- Simplex iterations limit |
|
|
589 if (PARAM.contains ("itlim")) |
|
|
590 { |
|
|
591 octave_value tmp = PARAM.contents (PARAM.seek ("itlim"))(0); |
|
|
592 lpxIntParam[5] = static_cast<int> (tmp.scalar_value ()); |
|
|
593 } |
|
5232
|
594 |
|
5234
|
595 //-- Simplex iterations count |
|
|
596 if (PARAM.contains ("itcnt")) |
|
|
597 { |
|
|
598 octave_value tmp = PARAM.contents (PARAM.seek ("itcnt"))(0); |
|
|
599 lpxIntParam[6] = static_cast<int> (tmp.scalar_value ()); |
|
|
600 } |
|
|
601 |
|
|
602 //-- Output frequency, in iterations |
|
|
603 if (PARAM.contains ("outfrq")) |
|
|
604 { |
|
|
605 octave_value tmp = PARAM.contents (PARAM.seek ("outfrq"))(0); |
|
|
606 lpxIntParam[7] = static_cast<int> (tmp.scalar_value ()); |
|
|
607 } |
|
|
608 |
|
|
609 //-- Branching heuristic option |
|
|
610 if (PARAM.contains("branch")) |
|
|
611 { |
|
|
612 octave_value tmp = PARAM.contents (PARAM.seek ("branch"))(0); |
|
|
613 double numtmp = tmp.scalar_value (); |
|
|
614 if (numtmp != 0 && numtmp != 1 && numtmp != 2) |
|
|
615 { |
|
|
616 OCTOUT << "'branch' parameter must be only (for MIP only):\n\t0 - branch on the first variable,\n\t1 - branch on the last variable,\n\t2 - branch using a heuristic by Driebeck and Tomlin [default]\n"; |
|
|
617 return retval; |
|
|
618 } |
|
|
619 lpxIntParam[14] = static_cast<int> (numtmp); |
|
5232
|
620 } |
|
|
621 |
|
5234
|
622 //-- Backtracking heuristic option |
|
|
623 if (PARAM.contains ("btrack")) |
|
|
624 { |
|
|
625 octave_value tmp = PARAM.contents (PARAM.seek ("btrack"))(0); |
|
|
626 double numtmp = tmp.scalar_value (); |
|
|
627 if (numtmp != 0 && numtmp != 1 && numtmp != 2) |
|
|
628 { |
|
|
629 OCTOUT << "'btrack' parameter must be only (for MIP only):\n\t0 - depth first search,\n\t1 - breadth first search,\n\t2 - backtrack using the best projection heuristic\n"; |
|
|
630 return retval; |
|
|
631 } |
|
|
632 lpxIntParam[15] = static_cast<int> (numtmp); |
|
5232
|
633 } |
|
|
634 |
|
5234
|
635 //-- Presolver option |
|
|
636 if (PARAM.contains ("presol")) |
|
|
637 { |
|
|
638 octave_value tmp = PARAM.contents (PARAM.seek ("presol"))(0); |
|
|
639 double numtmp = tmp.scalar_value (); |
|
|
640 if (numtmp != 0 && numtmp != 1) |
|
|
641 { |
|
|
642 OCTOUT << "'presol' parameter must be only:\n\t0 - LP presolver is ***NOT*** used,\n\t1 - LP presol is used\n"; |
|
|
643 return retval; |
|
|
644 } |
|
|
645 lpxIntParam[16] = static_cast<int> (numtmp); |
|
|
646 } |
|
|
647 |
|
|
648 //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
|
649 //-- Real parameters |
|
|
650 //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
|
651 |
|
|
652 //-- Ratio test option |
|
|
653 if (PARAM.contains ("relax")) |
|
|
654 { |
|
|
655 octave_value tmp = PARAM.contents (PARAM.seek ("relax"))(0); |
|
|
656 lpxRealParam[0] = tmp.scalar_value (); |
|
|
657 } |
|
5232
|
658 |
|
5234
|
659 //-- Relative tolerance used to check if the current basic solution |
|
|
660 //-- is primal feasible |
|
|
661 if (PARAM.contains ("tolbnd")) |
|
|
662 { |
|
|
663 octave_value tmp = PARAM.contents (PARAM.seek ("tolbn"))(0); |
|
|
664 lpxRealParam[1] = tmp.scalar_value (); |
|
|
665 } |
|
|
666 |
|
|
667 //-- Absolute tolerance used to check if the current basic solution |
|
|
668 //-- is dual feasible |
|
|
669 if (PARAM.contains ("toldj")) |
|
|
670 { |
|
|
671 octave_value tmp = PARAM.contents (PARAM.seek ("toldj"))(0); |
|
|
672 lpxRealParam[2] = tmp.scalar_value(); |
|
|
673 } |
|
5232
|
674 |
|
5234
|
675 //-- Relative tolerance used to choose eligible pivotal elements of |
|
|
676 //-- the simplex table in the ratio test |
|
|
677 if (PARAM.contains ("tolpiv")) |
|
|
678 { |
|
|
679 octave_value tmp = PARAM.contents (PARAM.seek ("tolpiv"))(0); |
|
|
680 lpxRealParam[3] = tmp.scalar_value (); |
|
|
681 } |
|
|
682 |
|
|
683 if (PARAM.contains ("objll")) |
|
|
684 { |
|
|
685 octave_value tmp = PARAM.contents (PARAM.seek ("objll"))(0); |
|
|
686 lpxRealParam[4] = tmp.scalar_value (); |
|
|
687 } |
|
|
688 |
|
|
689 if (PARAM.contains ("objul")) |
|
|
690 { |
|
|
691 octave_value tmp = PARAM.contents (PARAM.seek ("objul"))(0); |
|
|
692 lpxRealParam[5] = tmp.scalar_value (); |
|
|
693 } |
|
5232
|
694 |
|
5234
|
695 if (PARAM.contains ("tmlim")) |
|
|
696 { |
|
|
697 octave_value tmp = PARAM.contents (PARAM.seek ("tmlim"))(0); |
|
|
698 lpxRealParam[6] = tmp.scalar_value (); |
|
|
699 } |
|
|
700 |
|
|
701 if (PARAM.contains ("outdly")) |
|
|
702 { |
|
|
703 octave_value tmp = PARAM.contents (PARAM.seek ("outdly"))(0); |
|
|
704 lpxRealParam[7] = tmp.scalar_value (); |
|
|
705 } |
|
|
706 |
|
|
707 if (PARAM.contains ("tolint")) |
|
|
708 { |
|
|
709 octave_value tmp = PARAM.contents (PARAM.seek ("tolint"))(0); |
|
|
710 lpxRealParam[8] = tmp.scalar_value (); |
|
|
711 } |
|
|
712 |
|
|
713 if (PARAM.contains ("tolobj")) |
|
|
714 { |
|
|
715 octave_value tmp = PARAM.contents (PARAM.seek ("tolobj"))(0); |
|
|
716 lpxRealParam[9] = tmp.scalar_value (); |
|
|
717 } |
|
|
718 |
|
|
719 //-- 10th Input. If the problem is a LP problem you may select which solver |
|
|
720 //-- use: RSM (Revised Simplex Method) or IPM (Interior Point Method). |
|
|
721 //-- If the problem is a MIP problem this field will be ignored. |
|
|
722 octave_value tmp = args(9).scalar_value (); |
|
|
723 int lpsolver = static_cast<int> (tmp.scalar_value ()); |
|
|
724 |
|
|
725 //-- 11th Input. Saves a copy of the problem if SAVE<>0. |
|
|
726 tmp = args(10).scalar_value(); |
|
|
727 int save_pb = (tmp.scalar_value() != 0); |
|
5232
|
728 |
|
5234
|
729 //-- Assign pointers to the output parameters |
|
|
730 ColumnVector xmin (mrowsc); |
|
|
731 ColumnVector fmin (1); |
|
|
732 ColumnVector status (1); |
|
|
733 ColumnVector lambda (mrowsA); |
|
|
734 ColumnVector redcosts (mrowsc); |
|
|
735 ColumnVector time (1); |
|
|
736 ColumnVector mem (1); |
|
5232
|
737 |
|
5234
|
738 int jmpret = setjmp (mark); |
|
5235
|
739 |
|
5234
|
740 if (jmpret == 0) |
|
5235
|
741 glpk (sense, mrowsc, mrowsA, c, nz, rn.fortran_vec (), |
|
|
742 cn.fortran_vec (), a.fortran_vec (), b, ctype, |
|
|
743 freeLB.fortran_vec (), lb, freeUB.fortran_vec (), |
|
|
744 ub, vartype.fortran_vec (), isMIP, lpsolver, |
|
|
745 save_pb, xmin.fortran_vec (), fmin.fortran_vec (), |
|
|
746 status.fortran_vec (), lambda.fortran_vec (), |
|
|
747 redcosts.fortran_vec (), time.fortran_vec (), |
|
|
748 mem.fortran_vec ()); |
|
5232
|
749 |
|
5234
|
750 Octave_map extra; |
|
|
751 |
|
|
752 extra.assign ("lambda", octave_value (lambda)); |
|
|
753 extra.assign ("redcosts", octave_value (redcosts)); |
|
|
754 extra.assign ("time", octave_value (time)); |
|
|
755 extra.assign ("mem", octave_value (mem)); |
|
5232
|
756 |
|
5234
|
757 retval(3) = extra; |
|
|
758 retval(2) = octave_value(status); |
|
|
759 retval(1) = octave_value(fmin); |
|
|
760 retval(0) = octave_value(xmin); |
|
|
761 |
|
5235
|
762 #else |
|
|
763 |
|
|
764 gripe_not_supported ("glpk"); |
|
|
765 |
|
|
766 #endif |
|
|
767 |
|
5234
|
768 return retval; |
|
5232
|
769 } |
