# HG changeset patch # User jwe # Date 1111508306 0 # Node ID 9b776f5a33ebafcac64d1e6bfe226f424c5bf5ea # Parent f19646530e62462beaf202ba7109c4000a3c9436 [project @ 2005-03-22 16:16:30 by jwe] diff -r f19646530e62 -r 9b776f5a33eb scripts/ChangeLog --- a/scripts/ChangeLog Tue Mar 22 15:13:59 2005 +0000 +++ b/scripts/ChangeLog Tue Mar 22 16:18:26 2005 +0000 @@ -1,3 +1,11 @@ +2005-03-22 John W. Eaton + + * 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 * strings/split.m: Quick return for empty second arg. diff -r f19646530e62 -r 9b776f5a33eb scripts/Makefile.in --- a/scripts/Makefile.in Tue Mar 22 15:13:59 2005 +0000 +++ b/scripts/Makefile.in Tue Mar 22 16:18:26 2005 +0000 @@ -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) diff -r f19646530e62 -r 9b776f5a33eb scripts/optimization/Makefile.in --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/optimization/Makefile.in Tue Mar 22 16:18:26 2005 +0000 @@ -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 diff -r f19646530e62 -r 9b776f5a33eb scripts/optimization/glpk.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/optimization/glpk.m Tue Mar 22 16:18:26 2005 +0000 @@ -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 diff -r f19646530e62 -r 9b776f5a33eb scripts/optimization/glpkparams.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/optimization/glpkparams.m Tue Mar 22 16:18:26 2005 +0000 @@ -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) +% + diff -r f19646530e62 -r 9b776f5a33eb scripts/optimization/glpktest1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/optimization/glpktest1 Tue Mar 22 16:18:26 2005 +0000 @@ -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) diff -r f19646530e62 -r 9b776f5a33eb scripts/optimization/glpktest2 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/optimization/glpktest2 Tue Mar 22 16:18:26 2005 +0000 @@ -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) diff -r f19646530e62 -r 9b776f5a33eb src/ChangeLog --- a/src/ChangeLog Tue Mar 22 15:13:59 2005 +0000 +++ b/src/ChangeLog Tue Mar 22 16:18:26 2005 +0000 @@ -1,3 +1,8 @@ +2005-03-22 John W. Eaton + + * DLD-FUNCTIONS/__glpk__.cc: New file. + * Makefile.in (DLD_XSRC): Add it to the list. + 2005-03-21 John W. Eaton * octave.cc (maximum_braindamage): diff -r f19646530e62 -r 9b776f5a33eb src/DLD-FUNCTIONS/__glpk__.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/DLD-FUNCTIONS/__glpk__.cc Tue Mar 22 16:18:26 2005 +0000 @@ -0,0 +1,653 @@ +/*- -------------------------------------------------------------------- + -- + -- Copyright (C) 2005, Nicolo' Giorgetti, All rights reserved. + -- E-mail: . + -- + -- 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 +#include +#include + +//-- Octave headers +#include +#include +#include +#include + +//-- 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"< 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 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;imem_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 rn (mrowsA*mrowsc+1); + Array cn (mrowsA*mrowsc+1); + ColumnVector a (mrowsA*mrowsc+1, 0.0); + + int nz=0; + for(int i=0;i0) 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 freeLB (mrowsc); + for(int i=0;i freeUB (mrowsc); + for(int i=0;i 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; +} + diff -r f19646530e62 -r 9b776f5a33eb src/Makefile.in --- a/src/Makefile.in Tue Mar 22 15:13:59 2005 +0000 +++ b/src/Makefile.in Tue Mar 22 16:18:26 2005 +0000 @@ -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))