Mercurial > hg > octave-lyh
changeset 5232:9b776f5a33eb
[project @ 2005-03-22 16:16:30 by jwe]
| author | jwe |
|---|---|
| date | Tue, 22 Mar 2005 16:18:26 +0000 |
| parents | f19646530e62 |
| children | bdf892d3b024 |
| files | scripts/ChangeLog scripts/Makefile.in scripts/optimization/Makefile.in scripts/optimization/glpk.m scripts/optimization/glpkparams.m scripts/optimization/glpktest1 scripts/optimization/glpktest2 src/ChangeLog src/DLD-FUNCTIONS/__glpk__.cc src/Makefile.in |
| diffstat | 10 files changed, 1245 insertions(+), 3 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,11 @@ +2005-03-22 John W. Eaton <jwe@octave.org> + + * optimization: New directory. + * Makefile.in (SUBDIRS): Add it to the list. + * optimization/Makefile.in: New file. + * optimization/glpk.m, optimization/glpkparams.m, + optimization/glpktest1, optimization/glpktest2: New files. + 2005-03-16 Soren Hauberg <soren@hauberg.org> * strings/split.m: Quick return for empty second arg.
--- a/scripts/Makefile.in +++ b/scripts/Makefile.in @@ -30,8 +30,8 @@ skip-autoheader DOCSTRINGS SUBDIRS = audio control deprecated elfun finance general image io \ - linear-algebra miscellaneous plot polynomial quaternion \ - set signal sparse specfun special-matrix startup \ + linear-algebra miscellaneous optimization plot polynomial \ + quaternion set signal sparse specfun special-matrix startup \ statistics strings time DISTSUBDIRS = $(SUBDIRS)
new file mode 100644 --- /dev/null +++ b/scripts/optimization/Makefile.in @@ -0,0 +1,65 @@ +# +# Makefile for octave's scripts/optimization directory +# +# John W. Eaton +# jwe@che.utexas.edu +# Department of Chemical Engineering +# The University of Texas at Austin + +TOPDIR = ../.. + +script_sub_dir = optimization + +srcdir = @srcdir@ +top_srcdir = @top_srcdir@ +VPATH = @srcdir@ + +include $(TOPDIR)/Makeconf + +INSTALL = @INSTALL@ +INSTALL_PROGRAM = @INSTALL_PROGRAM@ +INSTALL_DATA = @INSTALL_DATA@ + +SOURCES = *.m + +EXTRAS = glpktest1 glpktest2 + +DISTFILES = Makefile.in $(SOURCES) $(EXTRAS) + +FCN_FILES = $(wildcard $(srcdir)/*.m) +FCN_FILES_NO_DIR = $(notdir $(FCN_FILES)) + +all: +.PHONY: all + +install install-strip: + $(do-script-install) +.PHONY: install install-strip + +uninstall: + $(do-script-uninstall) +.PHONY: uninstall + +clean: +.PHONY: clean + +tags: $(SOURCES) + ctags $(SOURCES) + +TAGS: $(SOURCES) + etags $(SOURCES) + +mostlyclean: clean +.PHONY: mostlyclean + +distclean: clean + rm -f Makefile +.PHONY: distclean + +maintainer-clean: distclean + rm -f tags TAGS +.PHONY: maintainer-clean + +dist: + ln $(DISTFILES) ../../`cat ../../.fname`/scripts/optimization +.PHONY: dist
new file mode 100644 --- /dev/null +++ b/scripts/optimization/glpk.m @@ -0,0 +1,309 @@ +function varargout=glpk(varargin) +## GLPK - An Octave Interface for the GNU GLPK library +## +## This routine calls the glpk library to solve an LP/MIP problem. A typical +## LP problem has following structure: +## +## [min|max] C'x +## s.t. +## Ax ["="|"<="|">="] b +## {x <= UB} +## {x >= LB} +## +## The calling syntax is: +## [XOPT,FOPT,STATUS,EXTRA]=glpkmex(SENSE,C,... +## A,B,CTYPE,LB,UB,... +## VARTYPE,PARAM,LPSOLVER,SAVE) +## +## For a quick reference to the syntax just type glpk at command prompt. +## +## The minimum number of input arguments is 4 (SENSE,C,A,B). In this case we +## assume all the constraints are '<=' and all the variables are continuous. +## +## --- INPUTS --- +## +## SENSE: indicates whether the problem is a minimization +## or maximization problem. +## SENSE = 1 minimize +## SENSE = -1 maximize. +## +## C: A column array containing the objective function +## coefficients. +## +## A: A matrix containing the constraints coefficients. +## +## B: A column array containing the right-hand side value for +## each constraint in the constraint matrix. +## +## CTYPE: A column array containing the sense of each constraint +## in the constraint matrix. +## CTYPE(i) = 'F' Free (unbounded) variable +## CTYPE(i) = 'U' "<=" Variable with upper bound +## CTYPE(i) = 'S' "=" Fixed Variable +## CTYPE(i) = 'L' ">=" Variable with lower bound +## CTYPE(i) = 'D' Double-bounded variable +## (This is case sensitive). +## +## LB: An array of at least length numcols containing the lower +## bound on each of the variables. +## +## UB: An array of at least length numcols containing the upper +## bound on each of the variables. +## +## VARTYPE: A column array containing the types of the variables. +## VARTYPE(i) = 'C' continuous variable +## VARTYPE(i) = 'I' Integer variable +## (This is case sensitive). +## +## PARAM: A structure containing some parameters used to define +## the behavior of solver. For more details type +## HELP GLPKPARAMS. +## +## LPSOLVER: Selects which solver using to solve LP problems. +## LPSOLVER=1 Revised Simplex Method +## LPSOLVER=2 Interior Point Method +## If the problem is a MIP problem this flag will be ignored. +## +## SAVE: Saves a copy of the problem if SAVE<>0. +## The file name can not be specified and defaults to "outpb.lp". +## The output file is CPLEX LP format. +## +## --- OUTPUTS --- +## +## XOPT: The optimizer. +## +## FOPT: The optimum. +## +## STATUS: Status of the optimization. +## +## - Simplex Method - +## Value Code +## 180 LPX_OPT solution is optimal +## 181 LPX_FEAS solution is feasible +## 182 LPX_INFEAS solution is infeasible +## 183 LPX_NOFEAS problem has no feasible solution +## 184 LPX_UNBND problem has no unbounded solution +## 185 LPX_UNDEF solution status is undefined +## +## - Interior Point Method - +## Value Code +## 150 LPX_T_UNDEF the interior point method is undefined +## 151 LPX_T_OPT the interior point method is optimal +## * Note that additional status codes may appear in +## the future versions of this routine * +## +## - Mixed Integer Method - +## Value Code +## 170 LPX_I_UNDEF the status is undefined +## 171 LPX_I_OPT the solution is integer optimal +## 172 LPX_I_FEAS solution integer feasible but +## its optimality has not been proven +## 173 LPX_I_NOFEAS no integer feasible solution +## +## EXTRA: A data structure containing the following fields: +## LAMBDA Dual variables +## REDCOSTS Reduced Costs +## TIME Time (in seconds) used for solving LP/MIP problem in seconds. +## MEM Memory (in bytes) used for solving LP/MIP problem. +## +## +## In case of error the glpkmex returns one of the +## following codes (these codes are in STATUS). For more informations on +## the causes of these codes refer to the GLPK reference manual. +## +## Value Code +## 204 LPX_E_FAULT unable to start the search +## 205 LPX_E_OBJLL objective function lower limit reached +## 206 LPX_E_OBJUL objective function upper limit reached +## 207 LPX_E_ITLIM iterations limit exhausted +## 208 LPX_E_TMLIM time limit exhausted +## 209 LPX_E_NOFEAS no feasible solution +## 210 LPX_E_INSTAB numerical instability +## 211 LPX_E_SING problems with basis matrix +## 212 LPX_E_NOCONV no convergence (interior) +## 213 LPX_E_NOPFS no primal feas. sol. (LP presolver) +## 214 LPX_E_NODFS no dual feas. sol. (LP presolver) +## + +% --- CODE --- +% If there is no input output the version and syntax +if nargin==0 + printf("glpk: An Octave interface to the GLPK library\n"); + printf("Version: 1.0\n"); + printf("\nSyntax: [xopt,fopt,status,extra]=glpk(sense,c,a,b,ctype,lb,ub,vartype,param,lpsolver,save)\n"); + return; +endif + +% If there are less than 4 arguments output an error message +if nargin<4 + error("At least 4 inputs required (sense,c,a,b)\n"); + return; +endif + +% At least 2 outputs required +if nargout<2 + error("2 outputs required\n"); + return; +endif + +% 1) Sense of optimization +sense=varargin{1}; + +% 2) Cost vector +c=varargin{2}; + +if (all(size(c)>1) | iscomplex(c) | ischar(c)) + error("C must be a real vector\n"); + return; +endif +nx=length(c); +if size(c,1) ~= nx + c=c'; +endif + +% 3) Matrix constraint +a=varargin{3}; +if isempty(a) + error("A cannot be an empty matrix\n"); + return; +endif +[nc, nxa]=size(a); +if (ischar(a) | iscomplex(a) | (nxa ~= nx)) + error("A must be a real valued %d by %d matrix\n",nc,nx); + return; +endif + +% 4) RHS +b=varargin{4}; +if isempty(b) + error("B cannot be an empty vector\n"); + return; +endif +if (ischar(b) | iscomplex(b) | (length(b) ~= nc)) + error("B must be a real valued %d by 1 vector\n",nc); + return; +endif + +% 5) Sense of each constraint +ctype=[]; +if nargin>4 + ctype=varargin{5}; + if isempty(ctype) + ctype=char('U'*ones(nc,1)); + elseif (isnumeric(ctype) | all(size(ctype)>1) | (length(ctype) ~= nc)) + error("CTYPE must be a char valued %d by 1 column vector\n",nc); + return; + else + for i=1:nc + if (ctype(i)!='F' & ctype(i)!='U' & ctype(i)!='S' & ctype(i)!='L' & ctype(i)!='D') + error("CTYPE must contain only F,U,S,L and D\n"); + return; + endif + endfor + endif +else + ctype=char('U'*ones(nc,1)); +end + +% 6) Vector with the lower bound of each variable +lb=[]; +if nargin>5 + lb=varargin{6}; + if isempty(lb) + lb=-Inf*ones(nx,1); + elseif (ischar(lb) | iscomplex(lb) | all(size(lb)>1) | (length(lb)~=nx)) + error("LB must be a real valued %d by 1 column vector\n",nx); + return; + endif +else + lb=-Inf*ones(nx,1); +end + +% 7) Vector with the upper bound of each variable +ub=[]; +if nargin>6 + ub=varargin{7}; + if isempty(ub) + ub=Inf*ones(nx,1); + elseif (ischar(ub) | iscomplex(ub) | all(size(ub)>1) | (length(ub)~=nx)) + error("UB must be a real valued %d by 1 column vector\n",nx); + return; + endif +else + ub=Inf*ones(nx,1); +end + +% 8) Vector with the type of variables +vartype=[]; +if nargin>7 + vartype=varargin{8}; + if isempty(vartype) + vartype=char('C'*ones(nx,1)); + elseif (isnumeric(vartype) | all(size(vartype)>1) | (length(vartype)~=nx)) + error("VARTYPE must be a char valued %d by 1 column vector\n",nx); + return; + else + for i=1:nx + if (vartype(i)!='C' & vartype(i)!='I') + error("VARTYPE must contain only C or I\n"); + return; + endif + endfor + endif +else + vartype=char('C'*ones(nx,1)); % As default we consider continuous vars +endif + +% 9) Parameters vector +param=[]; +if nargin>8 + param=varargin{9}; + if !isstruct(param) + error("PARAM must be a structure\n"); + return; + endif +else + param=struct; +endif + +% 10) Select solver method: simplex or interior point +lpsolver=[]; +if nargin>9 + lpsolver=varargin{10}; + if (!isnumeric(lpsolver) | all(size(lpsolver)>1)) + error("LPSOLVER must be a real scalar value\n"); + return; + endif +else + lpsolver=1; +endif + +% 11) Save the problem +savepb=[]; +if nargin>10 + savepb=varargin{11}; + if (!isnumeric(savepb) | all(size(savepb)>1)) + error("LPSOLVER must be a real scalar value\n"); + return; + endif +else + savepb=0; +end + +if nargin>11 + warning("Extra parameters ignored\n"); +endif + +try + [xopt, fmin, status, extra]=glpkoct(sense,c,a,b,ctype,lb,ub,vartype,param,lpsolver,savepb); +catch + error("Problems with glpkoct"); +end_try_catch + +varargout{1}=xopt; +varargout{2}=fmin; +varargout{3}=status; +varargout{4}=extra; + + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/optimization/glpkparams.m @@ -0,0 +1,159 @@ +% GLPKMEX Parameters list +% +% This document describes all control parameters currently implemented +% in the GLPK, an Octave interface for the GLPK library. Symbolic names +% of control parameters and corresponding codes of GLPK are given on the +% left. Types, default values, and descriptions are given on the right. +% +% ----------------------- +% 1 Integer parameters +% ----------------------- +% +% msglev type: integer, default: 1 (LPX_K_MSGLEV) +% Level of messages output by solver routines: +% 0 no output +% 1 error messages only +% 2 normal output +% 3 full output (includes informational messages) +% +% scale type: integer, default: 1 (LPX_K_SCALE) +% Scaling option: +% 0 no scaling +% 1 equilibration scaling +% 2 geometric mean scaling, then equilibration scaling +% +% dual type: integer, default: 0 (LPX_K_DUAL) +% Dual simplex option: +% 0 do not use the dual simplex +% 1 if initial basic solution is dual feasible, use the +% dual simplex +% +% price type: integer, default: 1 (LPX_K_PRICE) +% Pricing option (for both primal and dual simplex): +% 0 textbook pricing +% 1 steepest edge pricing +% +% round type: integer, default: 0 (LPX_K_ROUND) +% Solution rounding option: +% 0 report all primal and dual values "as is" +% 1 replace tiny primal and dual values by exact zero +% +% itlim type: integer, default: -1 (LPX_K_ITLIM) +% Simplex iterations limit. If this value is positive, it is +% decreased by one each time when one simplex iteration has +% been performed, and reaching zero value signals the solver +% to stop the search. Negative value means no iterations limit. +% +% itcnt type: integer, initial: 0 (LPX_K_ITCNT) +% Simplex iterations count.This count is increased by one +% each time when one simplex iteration has beenperformed. +% +% outfrq type: integer, default: 200 (LPX_K_OUTFRQ) +% Output frequency, in iterations. This parameter specifies +% how frequently the solver sends information about the solution +% to the standard output. +% +% branch type: integer, default: 2 (LPX_K_BRANCH) +% Branching heuristic option (for MIP only): +% 0 branch on the first variable +% 1 branch on the last variable +% 2 branch using a heuristic by Driebeck and Tomlin +% +% btrack type: integer, default: 2 (LPX_K_BTRACK) +% Backtracking heuristic option (for MIP only): +% 0 depth first search +% 1 breadth first search +% 2 backtrack using the best projection heuristic +% +% presol type: int, default: 1 (LPX_K_PRESOL) +% If this flag is set, the routine lpx_simplex solves the +% problem using the built-in LP presolver. Otherwise the LP +% presolver is not used. +% +% +% ----------------------- +% 2 Real parameters +% ----------------------- +% +% relax type: real, default: 0.07 (LPX_K_RELAX) +% Relaxation parameter used in the ratio test. If it is zero, the +% textbook ratio test is used. If it is non-zero (should be +% positive), Harris' two-pass ratio test is used. In the latter +% case on the first pass of the ratio test basic variables (in +% the case of primal simplex) or reduced costs of non-basic +% variables (in the case of dual simplex) are allowed to slightly +% violate their bounds, but not more than (RELAX ˇ TOLBND) or +% (RELAX ˇTOLDJ) (thus, RELAX is a percentage of TOLBND or TOLDJ). +% +% tolbnd type: real, default: 10e-7 (LPX_K_TOLBND) +% Relative tolerance used to check ifthe current basic solution +% is primal feasible (Do not change this parameter without detailed +% understanding its purpose). +% +% toldj type: real, default: 10e-7 (LPX_K_TOLDJ) +% Absolute tolerance used to check if the current basic solution +% is dual feasible (Do not change this parameter without detailed +% understanding its purpose). +% +% tolpiv type: real, default: 10e-9 (LPX_K_TOLPIV) +% Relative tolerance used to choose eligible pivotal elements of +% the simplex table (Do not change this parameter without detailed +% understanding its purpose). +% +% objll type: real, default: -DBL_MAX (LPX_K_OBJLL) +% Lower limit of the objective function.If on the phase II the +% objective function reaches this limit and continues decreasing, +% the solver stops the search.(Used in the dual simplex only) +% +% objul type: real, default: +DBL_MAX (LPX_K_OBJUL) +% Upper limit of the objective function. If on the phase II the +% objective function reaches this limit and continues increasing, +% the solver stops the search.(Used in the dual simplex only.) +% +% tmlim type: real, default: -1.0 (LPX_K_TMLIM) +% Searching time limit, in seconds. If this value is positive, +% it is decreased each time when one simplex iteration has been +% performed by the amount of time spent for the iteration, and +% reaching zero value signals the solver to stop the search. +% Negative value means no time limit. +% +% outdly type: real, default: 0.0 (LPX_K_OUTDLY) +% Output delay, in seconds. This parameter specifies how long +% the solver should delay sending information about the solution +% to the standard output. Non-positive value means no delay. +% +% tolint type: real, default: 10e-5 (LPX_K_TOLINT) +% Relative tolerance used to check ifthe current basic solution is +% integer feasible.(Do not change this parameter without detailed +% understanding its purpose). +% +% tolobj type: real, default: 10e-7 (LPX_K_TOLOBJ) +% Relative tolerance used to check if the value of the objective +% function is not better than in the best known integer feasible +% solution. (Do not change this parameter without detailed +% understanding its purpose) +% +% +% ----------------------- +% 3 Octave Example +% ----------------------- +% +% % Problem data +% s=-1; +% c=[10,6,4]'; +% a=[1,1,1;... +% 10,4,5;... +% 2,2,6]; +% b=[100,600,300]'; +% ctype=['U','U','U']'; +% lb=[0,0,0]'; +% ub=[]; +% vartype=['C','C','C']'; +% +% % Setting parameters +% param.msglev=1; % error messages only +% param.itlim=100; % Simplex iterations limit = 100 +% +% [xmin,fmin,status,lambda,extra]=glpkmex(s,c,a,b,ctype,lb,ub,vartype,param) +% +
new file mode 100644 --- /dev/null +++ b/scripts/optimization/glpktest1 @@ -0,0 +1,17 @@ +clear; + +disp('1st problem'); +s=-1; +c=[10,6,4]'; +a=[1,1,1;... + 10,4,5;... + 2,2,6]; +b=[100,600,300]'; +ctype=['U','U','U']'; +lb=[0,0,0]'; +ub=[]'; +vartype=['C','C','C']'; +param.msglev=3; +lpsolver=1; +save_pb=1; +[xmin,fmin,status,extra]=glpk(s,c,a,b,ctype,lb,ub,vartype,param,lpsolver,save_pb)
new file mode 100644 --- /dev/null +++ b/scripts/optimization/glpktest2 @@ -0,0 +1,25 @@ +clear; + +disp('2nd problem'); +s=1; +c=[-1,-1]'; +a=[-2,5;2,-2]; +b=[5;1]; +ctype=['U','U']'; +lb=[0;0]; ub=[]; +vartype=['I';'I']; +param.msglev=1; +[xmin,fmin,status,extra]=glpk(s,c,a,b,ctype,lb,ub,vartype,param) +pause; + +disp('3rd problem'); +s=1; +c=[0 0 0 -1 -1]'; +a=[-2 0 0 1 0;... + 0 1 0 0 2;... + 0 0 1 3 2]; +b=[4 12 18]'; +ctype=['S','S','S']'; +lb=[0,0,0,0,0]'; ub=[]; +vartype=['C','C','C','C','C']'; +[xmin,fmin,status,extra]=glpk(s,c,a,b,ctype,lb,ub,vartype)
--- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,8 @@ +2005-03-22 John W. Eaton <jwe@octave.org> + + * DLD-FUNCTIONS/__glpk__.cc: New file. + * Makefile.in (DLD_XSRC): Add it to the list. + 2005-03-21 John W. Eaton <jwe@octave.org> * octave.cc (maximum_braindamage):
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/__glpk__.cc @@ -0,0 +1,653 @@ +/*- -------------------------------------------------------------------- + -- + -- 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. + -- + -- ---------------------------------------------------------------------*/ + +#include <cfloat> +#include <csetjmp> +#include <ctime> + +//-- Octave headers +#include <oct.h> +#include <octave/ov-struct.h> +#include <octave/config.h> +#include <octave/error.h> + +//-- GLPK C header +extern "C"{ +#include "glpk.h" +} + +#define OCTOUT octave_stdout +#define OCTERR octave_stdout +#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 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 +}; + +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 glpkmex_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) +{ + 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) +{ + + LPX *lp; + int i,j; + int error; + clock_t t_start; + int typx=0; + int method; + + t_start = clock(); + + lib_set_fault_hook(NULL,glpkmex_fault_hook); + + + 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); + + //-- 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++){ + + //-- 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 + // -- structural variable. No constant term is assumed. + 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 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 + parameter will be ignored. + 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]); + + } + lpx_load_matrix(lp,nz,rn,cn,a); + + 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(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(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); + } + break; + 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=(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; + + lpx_delete_prob(lp); + return(0); + } + lpx_delete_prob(lp); + *status=(double)error; + return(error); +} + + + + + +DEFUN_DLD(glpkoct, args, nlhs, "glpkoct: OCT interface for the GLPK library. Don't use glpkoct, use glpk.m instead") +{ + // The list of values to return. See the declaration in oct-obj.h + octave_value_list retval; + + 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); + + 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; +// int *jc,*ir; +// double *pr; +// int nelc,count,row; +// +// /* NOTE: nnz is the actual number of nonzeros and is stored as the +// last element of the jc array where the size of the jc array is the +// number of columns + 1 */ +// nz = *(mxGetJc(A_IN) + mrowsc); +// jc = mxGetJc(A_IN); +// ir = mxGetIr(A_IN); +// pr = mxGetPr(A_IN); +// +// rn=(int *)calloc(nz+1,sizeof(int)); +// cn=(int *)calloc(nz+1,sizeof(int)); +// a=(double *)calloc(nz+1,sizeof(double)); +// +// count=0; row=0; +// for(i=1;i<=mrowsc;i++){ +// nelc=jc[i]-jc[i-1]; +// for(j=0;j<nelc;j++){ +// count++; +// rn[count]=ir[row]+1; +// cn[count]=i; +// a[count]=pr[row]; +// row++; +// } +// } +// } + + //-- 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; + } + + + //-- 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> 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]=(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; + } + + //-- 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; + } + + //-- 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; + } + + //-- Simplex iterations limit + if(PARAM.contains("itlim")){ + octave_value tmp=PARAM.contents(PARAM.seek("itlim"))(0); + lpxIntParam[5]=(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; + } + + //-- 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; + } + + //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + //-- Real parameters + //-- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + //-- 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 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(); } + + + //-- 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); + + 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 + + retval(3)=extra; + + return retval; +} +
--- a/src/Makefile.in +++ b/src/Makefile.in @@ -49,7 +49,8 @@ getpwent.cc getrusage.cc givens.cc hess.cc inv.cc kron.cc \ lpsolve.cc lsode.cc lu.cc minmax.cc pinv.cc qr.cc \ quad.cc qz.cc rand.cc schur.cc sort.cc sparse.cc spdet.cc \ - splu.cc spparms.cc sqrtm.cc svd.cc syl.cc time.cc gplot.l + splu.cc spparms.cc sqrtm.cc svd.cc syl.cc time.cc gplot.l \ + __glpk__.cc DLD_SRC := $(addprefix DLD-FUNCTIONS/, $(DLD_XSRC))