Mercurial > hg > octave-nkf
changeset 5234:a791b8b777e4
[project @ 2005-03-22 18:37:06 by jwe]
author | jwe |
---|---|
date | Tue, 22 Mar 2005 18:37:06 +0000 |
parents | bdf892d3b024 |
children | 5f0ad69b5c8c |
files | src/ChangeLog src/DLD-FUNCTIONS/__glpk__.cc |
diffstat | 2 files changed, 650 insertions(+), 540 deletions(-) [+] |
line wrap: on
line diff
--- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,6 +1,12 @@ 2005-03-22 John W. Eaton <jwe@octave.org> + + * DLD-FUNCTIONS/__glpk__.cc (F__glpk__, glpk): + Adapt to Octave coding style. + (glpk): Move decls closer to first use. + (F__glpk__): Eliminate unnecessary loop seting inf values. * DLD-FUNCTIONS/__glpk__.cc: New file. + * Makefile.in (DLD_XSRC): Add it to the list. 2005-03-21 John W. Eaton <jwe@octave.org>
--- a/src/DLD-FUNCTIONS/__glpk__.cc +++ b/src/DLD-FUNCTIONS/__glpk__.cc @@ -1,40 +1,39 @@ -/*- -------------------------------------------------------------------- - -- - -- Copyright (C) 2005, Nicolo' Giorgetti, All rights reserved. - -- E-mail: <giorgetti@dii.unisi.it>. - -- - -- This file is part of GLPKOCT an Octave interface to GLPK. - -- - -- GLPK is free software; you can redistribute it and/or modify it - -- under the terms of the GNU General Public License as published by - -- the Free Software Foundation; either version 2, or (at your option) - -- any later version. - -- - -- GLPK is distributed in the hope that it will be useful, but WITHOUT - -- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY - -- or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public - -- License for more details. - -- - -- You should have received a copy of the GNU General Public License - -- along with GLPK; see the file COPYING. If not, write to the Free - -- Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA - -- 02111-1307, USA. - -- - -- ---------------------------------------------------------------------*/ +/* + +Copyright (C) 2005 Nicolo' Giorgetti + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif #include <cfloat> #include <csetjmp> #include <ctime> -//-- Octave headers -#include <oct.h> -#include <octave/ov-struct.h> -#include <octave/config.h> -#include <octave/error.h> +#include "defun-dld.h" +#include "error.h" +#include "oct-obj.h" -//-- GLPK C header -extern "C"{ -#include "glpk.h" +extern "C" { +#include <glpk.h> } #define OCTOUT octave_stdout @@ -42,152 +41,154 @@ #define NIntP 17 #define NRealP 10 -int lpxIntParam[NIntP]= { - 1, - 1, - 0, - 1, - 0, - -1, - 0, - 200, - 1, - 2, - 0, - 1, - 0, - 0, - 2, - 2, - 1 +int lpxIntParam[NIntP] = { + 1, + 1, + 0, + 1, + 0, + -1, + 0, + 200, + 1, + 2, + 0, + 1, + 0, + 0, + 2, + 2, + 1 }; -int IParam[NIntP]={ - LPX_K_MSGLEV, - LPX_K_SCALE, - LPX_K_DUAL, - LPX_K_PRICE, - LPX_K_ROUND, - LPX_K_ITLIM, - LPX_K_ITCNT, - LPX_K_OUTFRQ, - LPX_K_MPSINFO, - LPX_K_MPSOBJ, - LPX_K_MPSORIG, - LPX_K_MPSWIDE, - LPX_K_MPSFREE, - LPX_K_MPSSKIP, - LPX_K_BRANCH, - LPX_K_BTRACK, - LPX_K_PRESOL +int IParam[NIntP] = { + LPX_K_MSGLEV, + LPX_K_SCALE, + LPX_K_DUAL, + LPX_K_PRICE, + LPX_K_ROUND, + LPX_K_ITLIM, + LPX_K_ITCNT, + LPX_K_OUTFRQ, + LPX_K_MPSINFO, + LPX_K_MPSOBJ, + LPX_K_MPSORIG, + LPX_K_MPSWIDE, + LPX_K_MPSFREE, + LPX_K_MPSSKIP, + LPX_K_BRANCH, + LPX_K_BTRACK, + LPX_K_PRESOL }; -double lpxRealParam[NRealP]={ - 0.07, - 1e-7, - 1e-7, - 1e-9, - -DBL_MAX, - DBL_MAX, - -1.0, - 0.0, - 1e-6, - 1e-7 +double lpxRealParam[NRealP] = { + 0.07, + 1e-7, + 1e-7, + 1e-9, + -DBL_MAX, + DBL_MAX, + -1.0, + 0.0, + 1e-6, + 1e-7 }; -int RParam[NRealP]={ - LPX_K_RELAX, - LPX_K_TOLBND, - LPX_K_TOLDJ, - LPX_K_TOLPIV, - LPX_K_OBJLL, - LPX_K_OBJUL, - LPX_K_TMLIM, - LPX_K_OUTDLY, - LPX_K_TOLINT, - LPX_K_TOLOBJ +int RParam[NRealP] = { + LPX_K_RELAX, + LPX_K_TOLBND, + LPX_K_TOLDJ, + LPX_K_TOLPIV, + LPX_K_OBJLL, + LPX_K_OBJUL, + LPX_K_TMLIM, + LPX_K_OUTDLY, + LPX_K_TOLINT, + LPX_K_TOLOBJ }; jmp_buf mark; //-- Address for long jump to jump to -int fperr; //-- Global error number +int fperr; //-- Global error number -int glpkmex_fault_hook(void *info, char *msg) -{ - OCTERR<<"*** SEVERE CRITICAL ERROR *** from GLPK !\n\n"<<msg<<" %s\n"; - longjmp( mark, -1 ); +int +glpk_fault_hook (void *info, char *msg) +{ + OCTERR << "*** SEVERE CRITICAL ERROR *** from GLPK !\n\n"<<msg<<" %s\n"; + longjmp (mark, -1); } -int glpkmex_print_hook(void *info, char *msg) +int +glpk_print_hook (void *info, char *msg) { - OCTERR<<msg<<"\n"; - return 1; + OCTERR << msg << "\n"; + return 1; } -int glpk(int sense,int n, int m, double *c,int nz,int *rn,int *cn, - double *a,double *b, char *ctype,int *freeLB, double *lb, - int *freeUB, double *ub, int *vartype, int isMIP, int lpsolver, - int save_pb, double *xmin, double *fmin, double *status, - double *lambda, double *redcosts, double *time, double *mem) +int +glpk (int sense, int n, int m, double *c, int nz, int *rn, int *cn, + double *a, double *b, char *ctype, int *freeLB, double *lb, + int *freeUB, double *ub, int *vartype, int isMIP, int lpsolver, + int save_pb, double *xmin, double *fmin, double *status, + double *lambda, double *redcosts, double *time, double *mem) { - - LPX *lp; - int i,j; - int error; - clock_t t_start; - int typx=0; - int method; + int error; + int typx = 0; + int method; + + clock_t t_start = clock(); + + lib_set_fault_hook (NULL, glpk_fault_hook); - t_start = clock(); - - lib_set_fault_hook(NULL,glpkmex_fault_hook); - + if (lpxIntParam[0] > 1) + lib_set_print_hook (NULL, glpk_print_hook); + + LPX *lp = lpx_create_prob (); + - if (lpxIntParam[0] > 1){ - lib_set_print_hook(NULL,glpkmex_print_hook); - } - - lp=lpx_create_prob(); - - - //-- Set the sense of optimization - if (sense==1) lpx_set_obj_dir(lp,LPX_MIN); - else lpx_set_obj_dir(lp,LPX_MAX); + //-- Set the sense of optimization + if (sense == 1) + lpx_set_obj_dir (lp, LPX_MIN); + else + lpx_set_obj_dir (lp, LPX_MAX); + + //-- If the problem has integer structural variables switch to MIP + if (isMIP) + lpx_set_class (lp, LPX_MIP); - //-- If the problem has integer structural variables switch to MIP - if(isMIP) lpx_set_class(lp,LPX_MIP); - - lpx_add_cols(lp,n); - for(i=0;i<n;i++){ - + lpx_add_cols (lp, n); + for (int i = 0; i < n; i++) + { //-- Define type of the structural variables - if (!freeLB[i] && !freeUB[i]){ - lpx_set_col_bnds(lp,i+1,LPX_DB,lb[i],ub[i]); - }else{ - if (!freeLB[i] && freeUB[i]){ - lpx_set_col_bnds(lp,i+1,LPX_LO,lb[i],ub[i]); - }else{ - if (freeLB[i] && !freeUB[i]){ - lpx_set_col_bnds(lp,i+1,LPX_UP,lb[i],ub[i]); - }else{ - lpx_set_col_bnds(lp,i+1,LPX_FR,lb[i],ub[i]); - } - } - } - // -- Set the objective coefficient of the corresponding + if (! freeLB[i] && ! freeUB[i]) + lpx_set_col_bnds (lp, i+1, LPX_DB, lb[i], ub[i]); + else + { + if (! freeLB[i] && freeUB[i]) + lpx_set_col_bnds (lp, i+1, LPX_LO, lb[i], ub[i]); + else + { + if (freeLB[i] && ! freeUB[i]) + lpx_set_col_bnds (lp, i+1, LPX_UP, lb[i], ub[i]); + else + lpx_set_col_bnds (lp, i+1, LPX_FR, lb[i], ub[i]); + } + } + + // -- Set the objective coefficient of the corresponding // -- structural variable. No constant term is assumed. - lpx_set_obj_coef(lp,i+1,c[i]); + lpx_set_obj_coef(lp,i+1,c[i]); - if(isMIP){ - lpx_set_col_kind(lp,i+1,vartype[i]); - } - } - - lpx_add_rows(lp,m); - for(i=0;i<m;i++){ - + if (isMIP) + lpx_set_col_kind (lp, i+1, vartype[i]); + } + + lpx_add_rows (lp, m); + + for (int i = 0; i < m; i++) + { /* If the i-th row has no lower bound (types F,U), the corrispondent parameter will be ignored. If the i-th row has no upper bound (types F,L), the corrispondent @@ -195,161 +196,218 @@ If the i-th row is of S type, the i-th LB is used, but the i-th UB is ignored. */ - switch(ctype[i]){ - case 'F': typx=LPX_FR; break; - case 'U': typx=LPX_UP; break; - case 'L': typx=LPX_LO; break; - case 'S': typx=LPX_FX; break; - case 'D': typx=LPX_DB; - } - lpx_set_row_bnds(lp,i+1,typx,b[i],b[i]); + + switch (ctype[i]) + { + case 'F': + typx = LPX_FR; + break; + + case 'U': + typx = LPX_UP; + break; + + case 'L': + typx = LPX_LO; + break; - } - lpx_load_matrix(lp,nz,rn,cn,a); + case 'S': + typx = LPX_FX; + break; - if (save_pb){ - if(lpx_write_cpxlp(lp, "outpb.lp") != 0){ - OCTERR<<"Unable to write problem\n"; - longjmp( mark, -1 ); - } - } + case 'D': + typx = LPX_DB; + break; + } + + lpx_set_row_bnds (lp, i+1, typx, b[i], b[i]); + + } + + lpx_load_matrix (lp, nz, rn, cn, a); - //-- scale the problem data (if required) - //-- if (scale && (!presol || method == 1)) lpx_scale_prob(lp); - //-- LPX_K_SCALE=IParam[1] LPX_K_PRESOL=IParam[16] - if (lpxIntParam[1] && (!lpxIntParam[16] || lpsolver!=1)){ - lpx_scale_prob(lp); - } - //-- build advanced initial basis (if required) - if (lpsolver == 1 && !lpxIntParam[16]){ - lpx_adv_basis(lp); - } - - for(i=0;i<NIntP;i++){ - lpx_set_int_parm(lp,IParam[i],lpxIntParam[i]); - } - for(i=0;i<NRealP;i++){ - lpx_set_real_parm(lp,RParam[i],lpxRealParam[i]); - } + if (save_pb) + { + if (lpx_write_cpxlp (lp, "outpb.lp") != 0) + { + OCTERR << "Unable to write problem\n"; + longjmp (mark, -1); + } + } + + //-- scale the problem data (if required) + //-- if (scale && (!presol || method == 1)) lpx_scale_prob(lp); + //-- LPX_K_SCALE=IParam[1] LPX_K_PRESOL=IParam[16] + if (lpxIntParam[1] && (! lpxIntParam[16] || lpsolver != 1)) + lpx_scale_prob (lp); + + //-- build advanced initial basis (if required) + if (lpsolver == 1 && ! lpxIntParam[16]) + lpx_adv_basis (lp); + + for(int i = 0; i < NIntP; i++) + lpx_set_int_parm (lp, IParam[i], lpxIntParam[i]); - if(lpsolver==1) method='S'; - else method='T'; + for (int i = 0; i < NRealP; i++) + lpx_set_real_parm (lp, RParam[i], lpxRealParam[i]); + + if (lpsolver == 1) + method = 'S'; + else + method = 'T'; - switch(method){ - case 'S': - if(isMIP){ - method='I'; - error=lpx_simplex(lp); - error=lpx_integer(lp); - }else{ - error=lpx_simplex(lp); - } + switch (method) + { + case 'S': + { + if (isMIP) + { + method = 'I'; + error = lpx_simplex (lp); + error = lpx_integer (lp); + } + else + error=lpx_simplex(lp); + } break; - case 'T': - error=lpx_interior(lp); - break; - default: - insist(method != method); - } + + case 'T': + error = lpx_interior(lp); + break; + + default: + insist (method != method); + } + + /* error assumes the following results: + error=0 <=> No errors + error=1 <=> Iteration limit exceeded. + error=2 <=> Numerical problems with basis matrix. + */ + if (error == LPX_E_OK) + { + if (isMIP) + { + *status = static_cast<double> (lpx_mip_status (lp)); + *fmin = lpx_mip_obj_val (lp); + } + else + { + if (lpsolver == 1) + { + *status = static_cast<double> (lpx_get_status (lp)); + *fmin = lpx_get_obj_val (lp); + } + else + { + *status = static_cast<double> (lpx_ipt_status (lp)); + *fmin = lpx_ipt_obj_val (lp); + } + } - /* - error assumes the following results: - error=0 <=> No errors - error=1 <=> Iteration limit exceeded. - error=2 <=> Numerical problems with basis matrix. - */ - if(error==LPX_E_OK){ - if(isMIP){ - *status=(double)lpx_mip_status(lp); - *fmin=lpx_mip_obj_val(lp); - }else{ - if(lpsolver==1){ - *status=(double)lpx_get_status(lp); - *fmin=lpx_get_obj_val(lp); - }else{ - *status=(double)lpx_ipt_status(lp); - *fmin=lpx_ipt_obj_val(lp); - } - } - if(isMIP){ - for(i=0;i<n;i++) xmin[i]=lpx_mip_col_val(lp,i+1); - }else{ - /* Primal values */ - for(i=0;i<n;i++){ - if(lpsolver==1) xmin[i]=lpx_get_col_prim(lp,i+1); - else xmin[i]=lpx_ipt_col_prim(lp,i+1); - } - /* Dual values */ - for(i=0; i<m; i++){ - if(lpsolver==1) lambda[i]=lpx_get_row_dual(lp,i+1); - else lambda[i]=lpx_ipt_row_dual(lp,i+1); - } - /* Reduced costs */ - for(i=0; i<lpx_get_num_cols(lp); i++){ - if(lpsolver==1) redcosts[i]=lpx_get_col_dual(lp,i+1); - else redcosts[i]=lpx_ipt_col_dual(lp,i+1); - } - } - *time=((double)(clock() - t_start))/CLOCKS_PER_SEC; - *mem=(double)lib_env_ptr()->mem_tpeak; + if (isMIP) + { + for (int i = 0; i < n; i++) + xmin[i] = lpx_mip_col_val (lp, i+1); + } + else + { + /* Primal values */ + for (int i = 0; i < n; i++) + { + if (lpsolver == 1) + xmin[i] = lpx_get_col_prim (lp, i+1); + else + xmin[i] = lpx_ipt_col_prim (lp, i+1); + } - lpx_delete_prob(lp); - return(0); - } - lpx_delete_prob(lp); - *status=(double)error; - return(error); + /* Dual values */ + for (int i = 0; i < m; i++) + { + if (lpsolver == 1) + lambda[i] = lpx_get_row_dual (lp, i+1); + else + lambda[i] = lpx_ipt_row_dual (lp, i+1); + } + + /* Reduced costs */ + for (int i = 0; i < lpx_get_num_cols (lp); i++) + { + if (lpsolver == 1) + redcosts[i] = lpx_get_col_dual (lp, i+1); + else + redcosts[i] = lpx_ipt_col_dual (lp, i+1); + } + } + + *time = static_cast<double> (clock () - t_start) / CLOCKS_PER_SEC; + *mem = static_cast<double> (lib_env_ptr () -> mem_tpeak); + + lpx_delete_prob (lp); + return 0; + } + + lpx_delete_prob (lp); + + *status= static_cast<double> (error); + + return error; } - - - - -DEFUN_DLD(glpkoct, args, nlhs, "glpkoct: OCT interface for the GLPK library. Don't use glpkoct, use glpk.m instead") +DEFUN_DLD (__glpk__, args, nlhs, + "__glpk__: internal interface for the GLPK library.\n\ +You should be using using glpk instead") { - // The list of values to return. See the declaration in oct-obj.h - octave_value_list retval; + // The list of values to return. See the declaration in oct-obj.h + octave_value_list retval; - int nrhs=args.length(); - - if(nrhs<1){ + int nrhs = args.length (); + + if (nrhs < 1) + { OCTERR<<"Use the script glpk for the optimization\n"; return retval; - } - - //-- 1st Input. Sense of optimization. - int sense; - double SENSE = args(0).scalar_value (); - if (SENSE>=0) sense=1; - else sense =-1; - - //-- 2nd Input. A column array containing the objective function - //-- coefficients. - int mrowsc=args(1).rows(); - - Matrix C(args(1).matrix_value()); - double *c=C.fortran_vec(); - - //-- 3rd Input. A matrix containing the constraints coefficients. - // If matrix A is NOT a sparse matrix - // if(!mxIsSparse(A_IN)){ - int mrowsA=args(2).rows(); - Matrix A(args(2).matrix_value()); // get the matrix - Array<int> rn (mrowsA*mrowsc+1); - Array<int> cn (mrowsA*mrowsc+1); - ColumnVector a (mrowsA*mrowsc+1, 0.0); + } + + //-- 1st Input. Sense of optimization. + int sense; + double SENSE = args(0).scalar_value (); + if (SENSE> = 0) + sense = 1; + else + sense = -1; + + //-- 2nd Input. A column array containing the objective function + //-- coefficients. + int mrowsc = args(1).rows(); + + Matrix C (args(1).matrix_value ()); + double *c = C.fortran_vec (); - int nz=0; - for(int i=0;i<mrowsA;i++){ - for(int j=0;j<mrowsc;j++){ - if(A(i,j) != 0){ - nz++; - rn(nz)=i+1; - cn(nz)=j+1; - a(nz)=A(i,j); - } - } - } + //-- 3rd Input. A matrix containing the constraints coefficients. + // If matrix A is NOT a sparse matrix + // if(!mxIsSparse(A_IN)){ + int mrowsA = args(2).rows(); + Matrix A (args(2).matrix_value ()); // get the matrix + Array<int> rn (mrowsA*mrowsc+1); + Array<int> cn (mrowsA*mrowsc+1); + ColumnVector a (mrowsA*mrowsc+1, 0.0); + + int nz = 0; + for (int i = 0; i < mrowsA; i++) + { + for (int j = 0; j < mrowsc; j++) + { + if (A(i,j) != 0) + { + nz++; + rn(nz) = i + 1; + cn(nz) = j + 1; + a(nz) = A(i,j); + } + } + } + // DON'T DELETE THIS PART... REPRESENTS THE SPARSE MATRICES MANIPULATION // }else{ // int i,j; @@ -382,272 +440,318 @@ // } // } - //-- 4th Input. A column array containing the right-hand side value - // for each constraint in the constraint matrix. - Matrix B(args(3).matrix_value()); - double *b=B.fortran_vec(); - - for(int i=0; i< mrowsA; i++){ - if (isinf(b[i]) && b[i]<0) b[i]=-octave_Inf; - if (isinf(b[i]) && b[i]>0) b[i]=octave_Inf; - } + //-- 4th Input. A column array containing the right-hand side value + // for each constraint in the constraint matrix. + Matrix B (args(3).matrix_value ()); + double *b = B.fortran_vec (); + + //-- 5th Input. A column array containing the sense of each constraint + //-- in the constraint matrix. + charMatrix CTYPE (args(4).char_matrix_value ()); + char *ctype = CTYPE.fortran_vec (); + + //-- 6th Input. An array of length mrowsc containing the lower + //-- bound on each of the variables. + Matrix LB (args(5).matrix_value ()); + double *lb = LB.fortran_vec (); + + //-- LB argument, default: Free + Array<int> freeLB (mrowsc); + for (int i = 0; i < mrowsc; i++) + { + if (isinf (lb[i])) + { + freeLB(i) = 1; + lb[i] = -octave_Inf; + } + else + freeLB(i) = 0; + } + + //-- 7th Input. An array of at least length numcols containing the upper + //-- bound on each of the variables. + Matrix UB (args(6).matrix_value ()); + + double *ub = UB.fortran_vec (); - - //-- 5th Input. A column array containing the sense of each constraint - //-- in the constraint matrix. - charMatrix CTYPE(args(4).char_matrix_value()); - - char *ctype=CTYPE.fortran_vec(); - - //-- 6th Input. An array of length mrowsc containing the lower - //-- bound on each of the variables. - Matrix LB(args(5).matrix_value()); - - double *lb=LB.fortran_vec(); - - //-- LB argument, default: Free - Array<int> freeLB (mrowsc); - for(int i=0;i<mrowsc;i++){ - if(isinf(lb[i])){ - freeLB(i)=1; - lb[i]=-octave_Inf; - }else freeLB(i)=0; - } - - //-- 7th Input. An array of at least length numcols containing the upper - //-- bound on each of the variables. - Matrix UB(args(6).matrix_value()); - - double *ub=UB.fortran_vec(); - - Array<int> freeUB (mrowsc); - for(int i=0;i<mrowsc;i++){ - if(isinf(ub[i])){ - freeUB(i)=1; - ub[i]=octave_Inf; - }else freeUB(i)=0; - } - - //-- 8th Input. A column array containing the types of the variables. - charMatrix VTYPE(args(7).char_matrix_value()); + Array<int> freeUB (mrowsc); + for (int i = 0; i < mrowsc; i++) + { + if (isinf (ub[i])) + { + freeUB(i) = 1; + ub[i] = octave_Inf; + } + else + freeUB(i) = 0; + } + + //-- 8th Input. A column array containing the types of the variables. + charMatrix VTYPE (args(7).char_matrix_value ()); - Array<int> vartype (mrowsc); - int isMIP; - for (int i = 0; i < mrowsc ; i++){ - if(VTYPE(i,0)=='I'){ - isMIP=1; - vartype(i)=LPX_IV; - }else{ - vartype(i)=LPX_CV; - } + Array<int> vartype (mrowsc); + int isMIP; + for (int i = 0; i < mrowsc ; i++) + { + if (VTYPE(i,0) == 'I') + { + isMIP = 1; + vartype(i) = LPX_IV; + } + else + { + vartype(i) = LPX_CV; + } + } + + //-- 9th Input. A structure containing the control parameters. + Octave_map PARAM = args(8).map_value (); + + //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + //-- Integer parameters + //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + //-- Level of messages output by the solver + if (PARAM.contains ("msglev")) + { + octave_value tmp = PARAM.contents (PARAM.seek ("msglev"))(0); + + double numtmp = tmp.scalar_value (); + if (numtmp != 0 && numtmp != 1 && numtmp != 2 && numtmp != 3) + { + 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"; + return retval; + } + + lpxIntParam[0] = static_cast<int> (numtmp); } - //-- 9th Input. A structure containing the control parameters. - Octave_map PARAM = args(8).map_value(); - - //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - //-- Integer parameters - //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - - //-- Level of messages output by the solver - if(PARAM.contains("msglev")){ - octave_value tmp=PARAM.contents(PARAM.seek("msglev"))(0); - - double numtmp=tmp.scalar_value(); - if((numtmp != 0) && (numtmp != 1) && (numtmp != 2) && (numtmp != 3)){ - 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"; - return retval; - } - lpxIntParam[0]=(int) numtmp; - } - - //-- scaling option - if(PARAM.contains("scale")){ - octave_value tmp=PARAM.contents(PARAM.seek("scale"))(0); - double numtmp=tmp.scalar_value(); - if((numtmp != 0) && (numtmp != 1) && (numtmp != 2)){ - OCTOUT<<"'scale' parameter must be only:\n\t0 - no scaling,\n\t1 - equilibration scaling,\n\t2 - geometric mean scaling\n"; - return retval; - } - lpxIntParam[1]=(int) numtmp; - } + //-- scaling option + if (PARAM.contains ("scale")) + { + octave_value tmp = PARAM.contents (PARAM.seek ("scale"))(0); + double numtmp = tmp.scalar_value (); + if (numtmp != 0 && numtmp != 1 && numtmp != 2) + { + OCTOUT << "'scale' parameter must be only:\n\t0 - no scaling,\n\t1 - equilibration scaling,\n\t2 - geometric mean scaling\n"; + return retval; + } + lpxIntParam[1] = static_cast<int> (numtmp); + } - //-- Dual dimplex option - if(PARAM.contains("dual")){ - octave_value tmp=PARAM.contents(PARAM.seek("dual"))(0); - double numtmp=tmp.scalar_value(); - if((numtmp != 0) && (numtmp != 1)){ - OCTOUT<<"'dual' parameter must be only:\n\t0 - do not use the dual simplex [default],\n\t1 - use dual simplex\n"; - return retval; - } - lpxIntParam[2]=(int) numtmp; - } + //-- Dual dimplex option + if (PARAM.contains ("dual")) + { + octave_value tmp = PARAM.contents (PARAM.seek ("dual"))(0); + double numtmp = tmp.scalar_value (); + if (numtmp != 0 && numtmp != 1) + { + OCTOUT<<"'dual' parameter must be only:\n\t0 - do not use the dual simplex [default],\n\t1 - use dual simplex\n"; + return retval; + } + lpxIntParam[2] = static_cast<int> (numtmp); + } - //-- Pricing option - if(PARAM.contains("price")){ - octave_value tmp=PARAM.contents(PARAM.seek("price"))(0); - double numtmp=tmp.scalar_value(); - if((numtmp != 0) && (numtmp != 1)){ - OCTOUT<<"'price' parameter must be only:\n\t0 - textbook pricing,\n\t1 - steepest edge pricing [default]\n"; - return retval; - } - lpxIntParam[3]=(int) numtmp; - } - - //-- Solution rounding option - if(PARAM.contains("round")){ - octave_value tmp=PARAM.contents(PARAM.seek("round"))(0); - double numtmp=tmp.scalar_value(); - if((numtmp != 0) && (numtmp != 1)){ - 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"; - return retval; - } - lpxIntParam[4]=(int) numtmp; + //-- Pricing option + if (PARAM.contains ("price")) + { + octave_value tmp = PARAM.contents (PARAM.seek ("price"))(0); + double numtmp = tmp.scalar_value(); + if (numtmp != 0 && numtmp != 1) + { + OCTOUT << "'price' parameter must be only:\n\t0 - textbook pricing,\n\t1 - steepest edge pricing [default]\n"; + return retval; + } + lpxIntParam[3] = static_cast<int> (numtmp); } - //-- Simplex iterations limit - if(PARAM.contains("itlim")){ - octave_value tmp=PARAM.contents(PARAM.seek("itlim"))(0); - lpxIntParam[5]=(int) tmp.scalar_value(); } + //-- Solution rounding option + if (PARAM.contains ("round")) + { + octave_value tmp = PARAM.contents (PARAM.seek ("round"))(0); + double numtmp = tmp.scalar_value (); + if (numtmp != 0 && numtmp != 1) + { + 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"; + return retval; + } + lpxIntParam[4] = static_cast<int> (numtmp); + } + + //-- Simplex iterations limit + if (PARAM.contains ("itlim")) + { + octave_value tmp = PARAM.contents (PARAM.seek ("itlim"))(0); + lpxIntParam[5] = static_cast<int> (tmp.scalar_value ()); + } - //-- Simplex iterations count - if(PARAM.contains("itcnt")){ - octave_value tmp=PARAM.contents(PARAM.seek("itcnt"))(0); - lpxIntParam[6]=(int) tmp.scalar_value(); } - - //-- Output frequency, in iterations - if(PARAM.contains("outfrq")){ - octave_value tmp=PARAM.contents(PARAM.seek("outfrq"))(0); - lpxIntParam[7]=(int) tmp.scalar_value(); } - - //-- Branching heuristic option - if(PARAM.contains("branch")){ - octave_value tmp=PARAM.contents(PARAM.seek("branch"))(0); - double numtmp=tmp.scalar_value(); - if((numtmp != 0) && (numtmp != 1) && (numtmp != 2)){ - 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"; - return retval; - } - lpxIntParam[14]=(int) numtmp; + //-- Simplex iterations count + if (PARAM.contains ("itcnt")) + { + octave_value tmp = PARAM.contents (PARAM.seek ("itcnt"))(0); + lpxIntParam[6] = static_cast<int> (tmp.scalar_value ()); + } + + //-- Output frequency, in iterations + if (PARAM.contains ("outfrq")) + { + octave_value tmp = PARAM.contents (PARAM.seek ("outfrq"))(0); + lpxIntParam[7] = static_cast<int> (tmp.scalar_value ()); + } + + //-- Branching heuristic option + if (PARAM.contains("branch")) + { + octave_value tmp = PARAM.contents (PARAM.seek ("branch"))(0); + double numtmp = tmp.scalar_value (); + if (numtmp != 0 && numtmp != 1 && numtmp != 2) + { + 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"; + return retval; + } + lpxIntParam[14] = static_cast<int> (numtmp); } - //-- Backtracking heuristic option - if(PARAM.contains("btrack")){ - octave_value tmp=PARAM.contents(PARAM.seek("btrack"))(0); - double numtmp=tmp.scalar_value(); - if((numtmp != 0) && (numtmp != 1) && (numtmp != 2)){ - 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"; - return retval; - } - lpxIntParam[15]=(int) numtmp; - } - - //-- Presolver option - if(PARAM.contains("presol")){ - octave_value tmp=PARAM.contents(PARAM.seek("presol"))(0); - double numtmp=tmp.scalar_value(); - if((numtmp != 0) && (numtmp != 1)){ - OCTOUT<<"'presol' parameter must be only:\n\t0 - LP presolver is ***NOT*** used,\n\t1 - LP presol is used\n"; - return retval; - } - lpxIntParam[16]=(int) numtmp; + //-- Backtracking heuristic option + if (PARAM.contains ("btrack")) + { + octave_value tmp = PARAM.contents (PARAM.seek ("btrack"))(0); + double numtmp = tmp.scalar_value (); + if (numtmp != 0 && numtmp != 1 && numtmp != 2) + { + 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"; + return retval; + } + lpxIntParam[15] = static_cast<int> (numtmp); } - //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - //-- Real parameters - //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + //-- Presolver option + if (PARAM.contains ("presol")) + { + octave_value tmp = PARAM.contents (PARAM.seek ("presol"))(0); + double numtmp = tmp.scalar_value (); + if (numtmp != 0 && numtmp != 1) + { + OCTOUT << "'presol' parameter must be only:\n\t0 - LP presolver is ***NOT*** used,\n\t1 - LP presol is used\n"; + return retval; + } + lpxIntParam[16] = static_cast<int> (numtmp); + } + + //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + //-- Real parameters + //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + //-- Ratio test option + if (PARAM.contains ("relax")) + { + octave_value tmp = PARAM.contents (PARAM.seek ("relax"))(0); + lpxRealParam[0] = tmp.scalar_value (); + } - //-- Ratio test option - if(PARAM.contains("relax")){ - octave_value tmp=PARAM.contents(PARAM.seek("relax"))(0); - lpxRealParam[0]=tmp.scalar_value(); } + //-- Relative tolerance used to check if the current basic solution + //-- is primal feasible + if (PARAM.contains ("tolbnd")) + { + octave_value tmp = PARAM.contents (PARAM.seek ("tolbn"))(0); + lpxRealParam[1] = tmp.scalar_value (); + } + + //-- Absolute tolerance used to check if the current basic solution + //-- is dual feasible + if (PARAM.contains ("toldj")) + { + octave_value tmp = PARAM.contents (PARAM.seek ("toldj"))(0); + lpxRealParam[2] = tmp.scalar_value(); + } - //-- Relative tolerance used to check if the current basic solution - //-- is primal feasible - if(PARAM.contains("tolbnd")){ - octave_value tmp=PARAM.contents(PARAM.seek("tolbn"))(0); - lpxRealParam[1]=tmp.scalar_value(); } + //-- Relative tolerance used to choose eligible pivotal elements of + //-- the simplex table in the ratio test + if (PARAM.contains ("tolpiv")) + { + octave_value tmp = PARAM.contents (PARAM.seek ("tolpiv"))(0); + lpxRealParam[3] = tmp.scalar_value (); + } + + if (PARAM.contains ("objll")) + { + octave_value tmp = PARAM.contents (PARAM.seek ("objll"))(0); + lpxRealParam[4] = tmp.scalar_value (); + } + + if (PARAM.contains ("objul")) + { + octave_value tmp = PARAM.contents (PARAM.seek ("objul"))(0); + lpxRealParam[5] = tmp.scalar_value (); + } - //-- Absolute tolerance used to check if the current basic solution - //-- is dual feasible - if(PARAM.contains("toldj")){ - octave_value tmp=PARAM.contents(PARAM.seek("toldj"))(0); - lpxRealParam[2]=tmp.scalar_value(); } - - //-- Relative tolerance used to choose eligible pivotal elements of - //-- the simplex table in the ratio test - if(PARAM.contains("tolpiv")){ - octave_value tmp=PARAM.contents(PARAM.seek("tolpiv"))(0); - lpxRealParam[3]=tmp.scalar_value(); } - - if(PARAM.contains("objll")){ - octave_value tmp=PARAM.contents(PARAM.seek("objll"))(0); - lpxRealParam[4]=tmp.scalar_value(); } - - if(PARAM.contains("objul")){ - octave_value tmp=PARAM.contents(PARAM.seek("objul"))(0); - lpxRealParam[5]=tmp.scalar_value(); } - - if(PARAM.contains("tmlim")){ - octave_value tmp=PARAM.contents(PARAM.seek("tmlim"))(0); - lpxRealParam[6]=tmp.scalar_value(); } - - if(PARAM.contains("outdly")){ - octave_value tmp=PARAM.contents(PARAM.seek("outdly"))(0); - lpxRealParam[7]=tmp.scalar_value(); } - - if(PARAM.contains("tolint")){ - octave_value tmp=PARAM.contents(PARAM.seek("tolint"))(0); - lpxRealParam[8]=tmp.scalar_value(); } - - if(PARAM.contains("tolobj")){ - octave_value tmp=PARAM.contents(PARAM.seek("tolobj"))(0); - lpxRealParam[9]=tmp.scalar_value(); } - + if (PARAM.contains ("tmlim")) + { + octave_value tmp = PARAM.contents (PARAM.seek ("tmlim"))(0); + lpxRealParam[6] = tmp.scalar_value (); + } + + if (PARAM.contains ("outdly")) + { + octave_value tmp = PARAM.contents (PARAM.seek ("outdly"))(0); + lpxRealParam[7] = tmp.scalar_value (); + } + + if (PARAM.contains ("tolint")) + { + octave_value tmp = PARAM.contents (PARAM.seek ("tolint"))(0); + lpxRealParam[8] = tmp.scalar_value (); + } + + if (PARAM.contains ("tolobj")) + { + octave_value tmp = PARAM.contents (PARAM.seek ("tolobj"))(0); + lpxRealParam[9] = tmp.scalar_value (); + } + + //-- 10th Input. If the problem is a LP problem you may select which solver + //-- use: RSM (Revised Simplex Method) or IPM (Interior Point Method). + //-- If the problem is a MIP problem this field will be ignored. + octave_value tmp = args(9).scalar_value (); + int lpsolver = static_cast<int> (tmp.scalar_value ()); + + //-- 11th Input. Saves a copy of the problem if SAVE<>0. + tmp = args(10).scalar_value(); + int save_pb = (tmp.scalar_value() != 0); - //-- 10th Input. If the problem is a LP problem you may select which solver - //-- use: RSM (Revised Simplex Method) or IPM (Interior Point Method). - //-- If the problem is a MIP problem this field will be ignored. - octave_value tmp=args(9).scalar_value(); - int lpsolver = (int) tmp.scalar_value(); - - //-- 11th Input. Saves a copy of the problem if SAVE<>0. - tmp=args(10).scalar_value(); - int save_pb = (tmp.scalar_value() != 0); + //-- Assign pointers to the output parameters + ColumnVector xmin (mrowsc); + ColumnVector fmin (1); + ColumnVector status (1); + ColumnVector lambda (mrowsA); + ColumnVector redcosts (mrowsc); + ColumnVector time (1); + ColumnVector mem (1); - //-- Assign pointers to the output parameters - ColumnVector xmin (mrowsc); - ColumnVector fmin (1); - ColumnVector status (1); - ColumnVector lambda (mrowsA); - ColumnVector redcosts (mrowsc); - ColumnVector time (1); - ColumnVector mem (1); - - int jmpret = setjmp( mark ); - if (jmpret==0){ - int error=glpk(sense,mrowsc,mrowsA,c,nz, - rn.fortran_vec(),cn.fortran_vec(),a.fortran_vec(),b,ctype, - freeLB.fortran_vec(),lb,freeUB.fortran_vec(),ub, - vartype.fortran_vec(),isMIP,lpsolver,save_pb, - xmin.fortran_vec(),fmin.fortran_vec(),status.fortran_vec(), - lambda.fortran_vec(),redcosts.fortran_vec(), - time.fortran_vec(),mem.fortran_vec()); - } + int jmpret = setjmp (mark); + if (jmpret == 0) + { + int error = glpk (sense, mrowsc, mrowsA, c, nz, rn.fortran_vec (), + cn.fortran_vec (), a.fortran_vec (), b, ctype, + freeLB.fortran_vec (), lb, freeUB.fortran_vec (), + ub, vartype.fortran_vec (), isMIP, lpsolver, + save_pb, xmin.fortran_vec (), fmin.fortran_vec (), + status.fortran_vec (), lambda.fortran_vec (), + redcosts.fortran_vec (), time.fortran_vec (), + mem.fortran_vec ()); + } - // Set the output parameters - retval(0)=octave_value(xmin); // returns xmin - retval(1)=octave_value(fmin); // returns fmin - retval(2)=octave_value(status); // returns status - - // Extra informations - Octave_map extra("lambda",octave_value(lambda)); //returns lambda - extra.assign("redcosts",octave_value(redcosts)); //returns the reduced costs - extra.assign("time",octave_value(time)); //time to solve the optimization pb - extra.assign("mem",octave_value(mem)); //memory used + Octave_map extra; + + extra.assign ("lambda", octave_value (lambda)); + extra.assign ("redcosts", octave_value (redcosts)); + extra.assign ("time", octave_value (time)); + extra.assign ("mem", octave_value (mem)); - retval(3)=extra; - - return retval; + retval(3) = extra; + retval(2) = octave_value(status); + retval(1) = octave_value(fmin); + retval(0) = octave_value(xmin); + + return retval; } -