If GAMS supported a subroutine call, it might take the form:
call sub(a,b,c,d)
where a, b, c and d are subroutine arguments, and sub.gms is a separate GAMS program. We can accomplish the same functionality using a $batinclude or $libinclude routine call:
$batinclude sub a b c d
The problem with a batinclude program is that it shares the calling program namespace.
While a user can write code which generates her own GDX data file, executes gams sub.gms and retrieves subroutine output, this sort of programming is tedious and error prone. The CALLBACK utility is designed to automate the process and provide improved reliability through type checking and the provision of diagnostic error messages. CALLBACK is tool which makes it easier to use GAMS programs as subroutines from a GAMS main program. Data is automatically transferred using GDX files between the main program and the subroutine argument list is verified for type, dimensionality and initialization.
Required arguments to a $batinclude routine using the CALLBACK approach are defined by position, and the optional arguments follow. The optional argument list consists of items which are uniquely defined by their type and dimension. In the current implementation, you therefore cannot have two optional arguments which are one-dimensional parameters. You may, however, have optional arguments which include one scalar, one vector and one two-dimensional matrix parameter, and you may have a similar sets of variables, equations and sets. (I hope that in the future this syntax can be improved. It works, but it is difficult to explain. I would much prefer to let the optional items be named, but we don't yet have this functionality for $batinclude and $libinclude calls in GAMS.)
For an example, see the interface to LAPACK's DGEEV routine: http://mpsge.org/callback/dgeev.zip
Here is a second example: an interface to Michael Saunders' LUSOL package.
$title Illustrate Use of the DGEEV Eigenvalue Interface
set i Dimension of the test matrix /i1*i100/;
* Row and column indices must be the same:
alias (i,j);
parameter a(i,j) Random matrix (row and column domain must be the same);
* Generate a reasonably sparse random matrix:
loop((i,j), a(i,j)$(max(0,uniform(0,1)-0.8)) = uniform(0,1););
* If required, make the matrix symmetric:
$if set symmetric a(i,j) = (a(i,j)+a(j,i))/2;
set ev Order set of labels for eigen values /ev1*ev15/,
evt Eigenvalue data types /
r Real,
cr Complex -- real component
ci Complex -- imaginary component /;
parameter eval(ev,evt) Eigen values,
evec(i,evt,ev) Eigen vectors;
$batinclude dgeev a eval evec ev
display eval;
option evec:3:1:2;
display evec;
$title DGEEV: A Wrapper for the LAPACK Eigenvalue Routine DGEEV
* Lib/batinclude Arguments:
* $batinclude dgeev a eval [evec] [ev]
* Required arguments:
* arg1: a(*,*) Matrix
* arg2: eval(ev,evt) Eigen values
* Optional arguments (in any order, distinguished by type and dimension)
* opt: evec(i,evt,ec) Eigen vectors -- 3 dimensional parameter
* opt: ev Eigen value labels -- 1 dimensional set
* ---------------------------------------------------------------------------------------------
* 1. The five dimensional tuple "interface" defines inputs and outputs:
* ---------------------------------------------------------------------------------------------
$ifthen.interface set interface
set argopt Required argument position or optional /arg1*arg7,opt/,
type Argument types /"set","par","var","equ"/,
dim Argument dimensions ("*"=variable) /"*",0*10/,
inout Argument direction /in,out,inout/;
set interface(argopt,*,type,dim,inout) Inputs and outputs /
"arg1"."A". "par"."2"."in",
"arg2"."eval". "par"."2"."out",
"opt"."evec". "par"."3"."out",
"opt"."ev". "set"."1"."in"/;
$exit
$endif.interface
* ---------------------------------------------------------------------------------------------
* 2. Boilerplate code which ultimately produces input GDX data, executes "gams dgeev" and
* retrieves output GDX data.
* ---------------------------------------------------------------------------------------------
$ifthen.callback not "%system.incparent%"==""
$setargs args
$batinclude "%system.fp%callback.gms" %args%
$exit
$endif.callback
* ---------------------------------------------------------------------------------------------
* 3. Standalone libinclude utility dgeev.gms begins here.
* ---------------------------------------------------------------------------------------------
set u(*) Row and column labels for a;
$if not exist 'dgeev.gdx' $abort 'Cannot find data file: dgeev.gdx'
$gdxin 'dgeev.gdx';
$loaddc u=*
*.$call 'copy /y dgeev.gdx dgeev_input.gdx';
$eval maxdim card(u)
set seq /1*%maxdim%/;
alias (u,u1,u2);
set i(u1) Rows in matrix a,
j(u2) Columns in matrix a;
parameter a(u1,u2) Matrix a;
$loaddc a
* Retrieve the matrix domain:
option i<a, j<a;
* Check that there are an equal number of rows and columns:
abort$(card(i)<>card(j)) "Error in DGEEV. Matrix a is not square:",i,j;
$ifthen.evdef set ev
set ev(u) Identifiers for the eigen values;
$loaddc ev
$else.evdef
set ev /1*%maxdim%/;
$endif.evdef
* Use the scratch directory for temporary files:
file kdat /%gams.scrdir%dgeev_input.scr/; put kdat;
kdat.nd=0; kdat.nw=0;
put card(i)/;
put card(a)/;
loop((i,j)$a(i,j),
kdat.nd=0; kdat.nw=0; kdat.nr=1;
put i.pos, ' ',j.pos, ' ';
kdat.nr=2; kdat.nd=8;
put a(i,j)/;
);
putclose;
execute '=dgeev "%gams.scrdir%dgeev_input.scr" "%gams.scrdir%dgeev_output.scr"';
abort$errorlevel 'Abort: call to dgeev.exe has failed. See %gams.wdir%dgeev.lst.';
execute '=gams "%gams.scrdir%dgeev_output.scr" gdx="%gams.scrdir%dgeev_gdx.scr" o="%gams.scrdir%dgeev_list.scr"';
abort$errorlevel 'Abort: call to gams has failed. See "%gams.scrdir%dgeev_list.scr"';
parameter wr(seq) Real component of eigenvalues,
wi(seq) Imaginary eigen values,
vr(seq,seq) Right eigen vectors;
execute_load '%gams.scrdir%dgeev_gdx.scr',wr,wi,vr;
set real(seq) Real eigenvalues,
complex(seq) Complex eigenvalues;
parameter nreal Number of real eigenvalues,
ncomplex Number of complex eigenvalues
ndim Number of dimensions;
real(seq)$(wr(seq)<>0 and wi(seq)=0) = yes;
* Store only one of each complex conjugate pair:
complex(seq) = yes$(wi(seq)<>0);
loop(seq$complex(seq), complex(seq+1) = no;);
nreal = card(real);
ncomplex = card(complex);
ndim = card(i);
display nreal, ncomplex, ndim;
set evt Complex eigen data type /
r Real,
cr Complex (real component)
ci Complex (imaginary component) /;
parameter rank(seq) Rank order of eigenvalues
modulus(seq) Rank order of complex eigenvalues;
* First rank order the real eigenvalues:
rank(real) = wr(real);
execute_unload '%gams.scrdir%rank.scr',rank;
execute '=gdxrank "%gams.scrdir%rank.scr" "%gams.scrdir%rank_output.scr"';
abort$errorlevel 'Abort: call to gdxrank.exe has failed.';
execute_load '%gams.scrdir%rank_output.scr', rank;
* Then rank complex eigenvalues using the modulus:
modulus(complex) = sign(wr(complex)) * sqrt(sqr(wr(complex)) + sqr(wi(complex)));
execute_unload '%gams.scrdir%rank.scr',modulus;
execute '=gdxrank "%gams.scrdir%rank.scr" "%gams.scrdir%rank_output.scr"';
abort$errorlevel 'Abort: call to gdxrank.exe has failed.';
execute_load '%gams.scrdir%rank_output.scr', modulus;
rank(complex) = modulus(complex)
set in(seq) Numeric indices for eigenvectors;
in(seq) = yes$(seq.val <= card(i));
parameter eval(ev,evt) Eigen values;
$ifthen.evecdef set evec
parameter evec(*,evt,ev) Eigen vectors;
loop((real(seq),ev)$(rank(seq) = ev.pos),
eval(ev,"r") = wr(seq);
loop((i,in)$(i.pos=in.val),
evec(i,"r",ev)$(not wi(seq)) = vr(in,seq);)
);
loop((complex(seq),ev)$(rank(seq)+nreal = ev.pos),
eval(ev,"cr") = wr(seq);
eval(ev,"ci") = wi(seq);
loop((i,in)$(i.pos=in.val),
evec(i,"cr",ev) = vr(in,seq);
evec(i,"ci",ev) = vr(in,seq+1);)
);
execute_unload 'dgeev.gdx',eval,evec;
$else.evecdef
loop((real(seq),ev)$(rank(seq) = ev.pos),
eval(ev,"r") = wr(seq););
loop((complex(seq),ev)$(rank(seq)+nreal = ev.pos),
eval(ev,"cr") = wr(seq);
eval(ev,"ci") = wi(seq););
execute_unload 'dgeev.gdx',eval;
$endif.evecdef
$title Program which Manages Callbacks to GAMS "Subroutines"
$ifthen.setup not "%system.incparent%"==""
$setnames "%system.incparent%" fp fn fa
$set gluecode "%fp%%fn%_gluecode.gms"
$ifthen.gluecode not exist "%gluecode%"
$call 'gams "%system.incname%" --program="%fn%" --input="%system.incparent%" --output="%gluecode%"'
$if errorlevel 1 $abort 'Error running %system.incname%'
$endif.gluecode
$setargs args
$batinclude %gluecode% %args%
$exit
$endif.setup
$title Glue Code Generator for GAMS Subroutine Callback
$log input : %input%
$log program : %program%
$log output : %output%
$if not set input $abort GLUECODEGEN error: input has not been defined.
$if not set program $abort GLUECODEGEN error: program has not been defined.
$if not set output $abort GLUECODEGEN error: output file has not been defined.
* An argument list can be provided to customize the call:
$if not set arglist $set arglist
* Retrieve the set defining the subroutine interface:
set argopt Argument assignments (and optional) /arg1*arg7,opt/,
arg(argopt) Argument list /arg1*arg7/,
type Argument types /"set","par","var","equ"/,
dim Argument dimensions /"*",0*10/,
inout Argument direction /in,out,inout/;
alias (u,*);
$set interface "%gams.scrdir%interface_gdx.scr"
$call '=gams "%input%" --interface="%interface%" gdx="%interface%"'
$if errorlevel 1 $abort 'Error generated by gams "%input%" --interface="%interface%". See listing file.'
set interface(argopt,u,type,dim,inout) Interface definition
$gdxin "%interface%"
$loaddc interface
set id(u) List of argument identifiers;
option id<interface;
* Verify logical consistency of the interface -- several consistency checks
parameter nreq(arg) Maximum one item assigned for each argument;
nreq(arg) = sum(interface(arg,id,type,dim,inout),1);
abort$(smax(arg,nreq(arg))>1) "More than one identifier assigned to argument list.";
set reqarg(arg) Required arguments;
reqarg(arg) = yes$(ord(arg)<=card(nreq));
set optarg(arg) Arguments which should not be assigned;
optarg(arg) = yes$(ord(arg)>card(reqarg) and ord(arg)<=card(interface));
abort$sum(optarg$nreq(optarg),1) "Error: required arguments are not assigned in sequence.", nreq,optarg;
set multipledef(u) Optional argument which with multiple definitions;
multipledef(id) = yes$(sum(interface(argopt,id,type,dim,inout),1) > 1);
abort$card(multipledef) multipledef;
parameter nopt(type,dim) Number of optional items of a particular type and dimension (maximum=1!);
nopt(type,dim) = sum(interface("opt",id,type,dim,inout),1);
set notunique(u) Optional arguments which are not uniquely defined by type and dimension;
loop(interface("opt",id,type,dim,inout)$(nopt(type,dim)>1), notunique(id) = yes;);
option interface:0:0:1;
abort$card(notunique) interface,notunique;
* Generate the controller code:
file kgms /%output%/; put kgms; kgms.lw=0; kgms.nd=0; kgms.nw=0;
put '$setargs '; loop(arg, put arg.tl,' ';); put //;
put '$setlocal args '; loop(arg, put '%',arg.tl,'% ';); put /;
put /'* Verify no extra aruments are provided:'//;
put '$if not "%arg',(card(reqarg)+card(optarg)+1),'%"=="" $abort Error in invocation -- too many arguments: %args%'/;
put /'* Verify consistency of required arguments:'/;
loop(interface(reqarg(arg),id,type,dim,inout),
put /'$ifthen"%',arg.tl,'%"==""'/;
put '$abort Error in %program% call: missing argument ',arg.tl,' = ',id.tl,'.'/;
put '$endif'/;
put '$if not ',type.tl,'type %',arg.tl,'% $abort Error in %program%: ',id.tl,'=%',arg.tl,'% is not ',dim.tl,'-dimensional ',type.tl/;
if (not sameas(dim,"*"),
put '$if not dimension ',dim.tl,' %',arg.tl,'% $abort Error in %program%: ',id.tl,'=%',arg.tl,'% is not ',dim.tl,'-dimensional ',type.tl/; ));;
put /'* Find ids to which optional arguments have been assigned:'/;
loop(optarg, put '$setlocal ',optarg.tl,'id'/;);
loop(optarg,
put /'$ifthen.',optarg.tl,' not "%',optarg.tl,'%"==""'/;
loop(interface("opt",id,type,dim,inout),
put /'$ifthen.',id.tl,'type',optarg.tl,' ',type.tl,'type %',optarg.tl,'%'/;
if (sameas(dim,"*"),
put '$setlocal ',optarg.tl,'id ',id.tl/;
else
put '$ifthen.',id.tl,'dim', optarg.tl,' dimension ',dim.tl,' %',optarg.tl,'%'/;
put '$setlocal ',optarg.tl,'id ',id.tl/;
put '$endif.',id.tl,'dim', optarg.tl/;
);
put '$endif.',id.tl,'type',optarg.tl/;);
put /'$ifthen.empty',optarg.tl,' "%',optarg.tl,'%"==""'/;
put '$abort " Error in %program%: argument ',optarg.tl,'=%',optarg.tl,'% does not match optional arguments: ';
loop(interface("opt",id,type,dim,inout),
put id.tl,' (',dim.tl,' dim ',type.tl,'), ';); put @(file.cc-2),' '/;
put '$endif.empty',optarg.tl/;
put /'$endif.',optarg.tl/;)
put /'* Verify that all inputs have been defined:'/;
set ininout(inout)/in,inout/;
loop(interface(reqarg,id,type,dim,ininout),
put '$if not defined %',reqarg.tl,'% $abort "Input item not defined: %',reqarg.tl,'% (= ',id.tl,')"'/;);
loop(interface("opt",id,type,dim,ininout),
loop(optarg,
put '$ifthen.defined_',optarg.tl,' "%',optarg.tl,'id%"=="',id.tl,'"'/;
put '$if not defined %',optarg.tl,'% $abort "Input item not defined: %',optarg.tl,'% (= ',id.tl,')"'/;
put '$endif.defined_',optarg.tl/; ));
put /'* Write a GDX file with the program inputs:'/;
put '$setlocal inputs ';
loop(interface(reqarg,id,type,dim,ininout), put '%',reqarg.tl,'%=',id.tl,', ';); put @(file.cc-2),' '/;
loop(interface("opt",id,type,dim,ininout),
loop(optarg,
put '$if "%',optarg.tl,'id%"=="',id.tl,'" $setlocal inputs %inputs%, %',optarg.tl,'%=',id.tl/;));
put "execute_unload '%program%.gdx', %",'inputs%;'/;
put /'$log Entering %program% with %inputs%'/;
put /'* Generate the command line argument list:'/;
put '$setglobal ide "ide=%','gams.ide% lo=%','gams.lo% errorlog=%','gams.errorlog% errmsg=1"'/;
put /'* Provide flags for optional arguments:'/;
loop(interface("opt",id,type,dim,ininout),
loop(optarg,
put '$if not "%',optarg.tl,'id%"=="" $set ide %','ide% --%',optarg.tl,'id%=',optarg.tl/;));
put /'* Execute the program:'/;
put '$log --- Executing ',"'",'=gams "%program%" %','ide%',"'"/;
put 'execute ',"'",'=gams "%program%" %','ide%',"';"/;
put 'abort$errorlevel ',"'Abort: call to %program% has failed. See %","gams.wdir%","%program%.lst.'"/;
put /'* Clear the items to be returned:'/;
set outinout(inout)/out,inout/;
loop(interface(reqarg,u,type,dim,outinout), put 'option clear=%',reqarg.tl,'%'/;)
loop(interface("opt",id,type,dim,outinout),
loop(optarg,
put '$if "%',optarg.tl,'id%"=="',id.tl,'" option clear=%',optarg.tl,'%'/;));
put /'* Read the program outputs:'/;
put '$setlocal outputs ';
loop(interface(reqarg,id,type,dim,outinout), put '%',reqarg.tl,'%=',id.tl,', ';); put @(file.cc-2),' '/;
loop(interface("opt",id,type,dim,outinout),
loop(optarg,
put '$if "%',optarg.tl,'id%"=="',id.tl,'" $setlocal outputs %outputs%, %',optarg.tl,'%=',id.tl/;));
put /'$log Exiting %program% with %outputs%'/;
put "execute_load '%program%.gdx', %",'outputs%;'/;
$setargs arg1 arg2 arg3 arg4 arg5 arg6 arg7
$setlocal args %arg1% %arg2% %arg3% %arg4% %arg5% %arg6% %arg7%
* Verify no extra aruments are provided:
$if not "%arg5%"=="" $abort Error in invocation -- too many arguments: %args%
* Verify consistency of required arguments:
$ifthen"%arg1%"==""
$abort Error in dgeev call: missing argument arg1 = A.
$endif
$if not partype %arg1% $abort Error in dgeev: A=%arg1% is not 2-dimensional par
$if not dimension 2 %arg1% $abort Error in dgeev: A=%arg1% is not 2-dimensional par
$ifthen"%arg2%"==""
$abort Error in dgeev call: missing argument arg2 = eval.
$endif
$if not partype %arg2% $abort Error in dgeev: eval=%arg2% is not 2-dimensional par
$if not dimension 2 %arg2% $abort Error in dgeev: eval=%arg2% is not 2-dimensional par
* Find ids to which optional arguments have been assigned:
$setlocal arg3id
$setlocal arg4id
$ifthen.arg3 not "%arg3%"==""
$ifthen.evectypearg3 partype %arg3%
$ifthen.evecdimarg3 dimension 3 %arg3%
$setlocal arg3id evec
$endif.evecdimarg3
$endif.evectypearg3
$ifthen.evtypearg3 settype %arg3%
$ifthen.evdimarg3 dimension 1 %arg3%
$setlocal arg3id ev
$endif.evdimarg3
$endif.evtypearg3
$ifthen.emptyarg3 "%arg3%"==""
$abort " Error in dgeev: argument arg3=%arg3% does not match optional arguments: evec (3 dim par), ev (1 dim set)
$endif.emptyarg3
$endif.arg3
$ifthen.arg4 not "%arg4%"==""
$ifthen.evectypearg4 partype %arg4%
$ifthen.evecdimarg4 dimension 3 %arg4%
$setlocal arg4id evec
$endif.evecdimarg4
$endif.evectypearg4
$ifthen.evtypearg4 settype %arg4%
$ifthen.evdimarg4 dimension 1 %arg4%
$setlocal arg4id ev
$endif.evdimarg4
$endif.evtypearg4
$ifthen.emptyarg4 "%arg4%"==""
$abort " Error in dgeev: argument arg4=%arg4% does not match optional arguments: evec (3 dim par), ev (1 dim set)
$endif.emptyarg4
$endif.arg4
* Verify that all inputs have been defined:
$if not defined %arg1% $abort "Input item not defined: %arg1% (= A)"
$ifthen.defined_arg3 "%arg3id%"=="ev"
$if not defined %arg3% $abort "Input item not defined: %arg3% (= ev)"
$endif.defined_arg3
$ifthen.defined_arg4 "%arg4id%"=="ev"
$if not defined %arg4% $abort "Input item not defined: %arg4% (= ev)"
$endif.defined_arg4
* Write a GDX file with the program inputs:
$setlocal inputs %arg1%=A
$if "%arg3id%"=="ev" $setlocal inputs %inputs%, %arg3%=ev
$if "%arg4id%"=="ev" $setlocal inputs %inputs%, %arg4%=ev
execute_unload 'dgeev.gdx', %inputs%;
$log Entering dgeev with %inputs%
* Generate the command line argument list:
$setglobal ide "ide=%gams.ide% lo=%gams.lo% errorlog=%gams.errorlog% errmsg=1"
* Provide flags for optional arguments:
$if not "%arg3id%"=="" $set ide %ide% --%arg3id%=arg3
$if not "%arg4id%"=="" $set ide %ide% --%arg4id%=arg4
* Execute the program:
$log --- Executing '=gams "dgeev" %ide%'
execute '=gams "dgeev" %ide%';
abort$errorlevel 'Abort: call to dgeev has failed. See %gams.wdir%dgeev.lst.'
* Clear the items to be returned:
option clear=%arg2%
$if "%arg3id%"=="evec" option clear=%arg3%
$if "%arg4id%"=="evec" option clear=%arg4%
* Read the program outputs:
$setlocal outputs %arg2%=eval
$if "%arg3id%"=="evec" $setlocal outputs %outputs%, %arg3%=evec
$if "%arg4id%"=="evec" $setlocal outputs %outputs%, %arg4%=evec
$log Exiting dgeev with %outputs%
execute_load 'dgeev.gdx', %outputs%;