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