LUSOL maintains LU factors of a square or rectangular sparse matrix. The website for LUSOL is: http://stanford.edu/group/SOL/software/lusol/ The code is distributed under the terms of the [Common Public License][CPL].
@echo off
gfortran -fPIC -O3 -c lusol_precision.f90
gfortran -fPIC -O3 -c lusol.f90
gfortran -fPIC -O3 -c main.f90
gfortran main.o lusol_precision.o lusol.o -o lusol.exe
Source files and Windows executable:
29-Oct-2014 8:03:22p 227,520 lusol.f90 29-Oct-2014 8:03:22p 1,286 lusol_precision.f90 27-Nov-2014 9:48:34p 13,801 main.f90 27-Nov-2014 11:17:24p 237,645 lusol.exe
$title Demonstrate use of LUSOL interface
Set i Row indices /1*4/,
j Column indices /a*d/;
* NB: Sets i and j must have the same dimension,
* but they need not have the same set of elements.
* Define some options here:
$set rp_eps 2.2E-016
$eval small %rp_eps%**0.8
$eval utol1 %rp_eps%**0.67
$eval utol2 %rp_eps%**0.67
$eolcom !
parameter luopt(*) /
lprint 10 ! Print level
* Print level:
* < 0 suppresses output.
* = 0 gives error messages.
* >= 10 gives statistics about the LU factors.
* >= 50 gives debug output from lu1fac
* (the pivot row and column and the
* no. of rows and columns involved at
* each elimination step).
pivotscheme 0 ! Pivoting scheme index (0=TPP)
* 0 TPP Threshold Partial Pivoting.
* 1 TRP Threshold Rook Pivoting.
* 2 TCP Threshold Complete Pivoting.
* 3 TSP Threshold Symmetric Pivoting.
* 4 TDP Threshold Diagonal Pivoting.
maxcol 5 ! Maximum number of columns searched
Ltol1 10 ! Max Lij allowed during factor
Ltol2 10 ! Max Lij allowed during updates
small %small% ! Absolute tolerance for treating reals as zero
Utol1 %Utol1% ! Absolute tol for flagging small diagonals of U
Utol2 %Utol2% ! Relative tol for flagging small diagonals of U
Uspace 3 ! Factor limiting waste space in U
dens1 0.3 ! Density at which the Markowitz should reach maxcol columns
dens2 0.1 ! Density at which the Markowitz strategy search only 1 column
/;
$gdxout 'luoptions.gdx'
$unload luopt
$gdxout
parameter a(i,j) Random matrix to be factorized;
a(i,j) = uniform(0,1);
alias (i,ii);
parameter ident(i,ii) Identity matrix (RHS for extracing inverse),
ainv(j,i) Inverse of a;
ident(i,i) = 1;
* First, do a factor-solve:
$batinclude lusolve a ident ainv
parameter chk Check of LUSOL solution;
chk(i,ii) = sum(j, a(i,j)*ainv(j,ii));
display "Factor/solve:", ainv, chk;
* Next, do a factorization with options file written above
* and save the factorization in lusol.bin. No rhs vector
* is used in this invocation:
$set cmdline --luopt=luoptions.gdx
$batinclude lufactor a lusol.bin
* Third, do a solve with the existing factors. Even though we use
* factors from lusol.bin, lusol.gms requires that the matrix a
* be provided so that the solution can be verified.
* N.B. The --luopt only affects factorization. It is ignored when
* solving from an existing LU factorization:
* After saving the temporary file, the LUSOL.EXE program can
* be invoked directly as:
* lusol "<dir>\lusol_solve.scr"
* where <dir> is the scratch directory, typically named 225b
$set cmdline keep=1
$batinclude lusolve a ident ainv lusol.bin
parameter chk Check of LUSOL solution;
chk(i,ii) = sum(j, a(i,j)*ainv(j,ii));
display "Factor and then solve:", ainv, chk;
$title LUSOLVE: LUSOL Wrapper to Solve a Linear System
$setargs a b x lu
$if not set cmdline $set cmdline
$if not "%lu%"=="" $set cmdline %cmdline% --lu=%lu%
$batinclude lusol %a% %b%
execute_load 'lusol.gdx',%x%=x;
$title LUFACTOR: LUSOL Wrapper to Factor a Matrix (generates LUSOL.BIN)
$setargs a lu
$if not set cmdline $set cmdline
$if not "%lu%"=="" $set cmdline %cmdline% --lu="%lu%"
$batinclude lusol %a%
$title GAMS Interface to LUSOL
* Usage:
* $set arglist lu=<lu>
* $batinclude lusol a [b]
* Required arguments:
* arg1: a(*,*) Matrix
* Optional arguments (in any order, distinguished by type and dimension)
* opt: b(i[,k]) RHS vector or array (1 or more dimensional parameter)
* NB: - If matrix b is omitted, matrix a must be provided
* and no LU factors may be specified.
* - If LU factors are provided, matrix b must be provided
* and the existing factorization is employed.
* Contents of lusol.gdx following a factorization (no b):
* i Row indices of a and b
* j Column indices of a
* a(i,j) Original matrix
* fstatus(info) LUSOL status for factor
* Contents of lusol.gdx following a factor/solve:
* k Column domain of b
* b(i[,k]) RHS vector(s)
* x(j[,k]) Solution vector
* sstatus(info[,k]) LUSOL status for solve
* The following is an identity satisfied by the LUSOL output data:
* sum(j, a(i,j)*x(j,k)) = b(i,k);
* ---------------------------------------------------------------------------------------------
* 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",
"opt" ."b". "par"."*"."in" /;
$exit
$endif.interface
* ---------------------------------------------------------------------------------------------
* 2. Boilerplate code which ultimately produces input GDX data, executes "gams lusol" 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 lusol.gms begins here.
* ---------------------------------------------------------------------------------------------
$ifthen.lu set lu
$if not exist %lu% $abort "LUSOL error: LU factor file (%lu%) does not exist."
$endif.lu
alias (u,u1,u2,*);
set i(u1) Rows in matrix a,
ib(u) Rows in b,
noti(u) Rows in b which are not in a,
j(u2) Columns in matrix a;
parameter aa(u1,u2) Matrix a;
set luoptions LU options /
lprint "Print level",
maxcol "Maximum number of columns searched",
keepLU "Means the numerical factors will be computed",
Ltol1 "Max Lij allowed during factor",
Ltol2 "Max Lij allowed during updates",
small "Absolute tolerance for treating reals as zero",
Utol1 "Absolute tol for flagging small diagonals of U",
Utol2 "Relative tol for flagging small diagonals of U",
Uspace "Factor limiting waste space in U",
dens1 "Density at which the Markowitz should reach maxcol columns",
dens2 "Density at which the Markowitz strategy search only 1 column",
pivotscheme "Pivoting scheme index (0=TPP)" /;
set info Information returned by the factorization routine /
inform "Return code from last call to any LU routine.",
nsing "No. of singularities marked in the output array w(*).",
jsing "Column index of last singularity.",
minlen "Minimum recommended value for lena.",
maxlen "?",
lena "Allocation for LENA at factorization",
nelem "Number of nonzeros at factorization",
nupdat "No. of updates performed by the lu8 routines.",
nrank "No. of nonempty rows of U.",
ndens1 "No. of columns remaining when the density of the factorized reached dens1.",
ndens2 "No. of columns remaining when the density of the matrix being factorized reached dens2."
jumin "The column index associated with DUmin.",
numL0 "No. of columns in initial L.",
lenL0 "Size of initial L (no. of nonzeros).",
lenU0 "Size of initial U.",
lenL "Size of current L.",
lenU "Size of current U.",
lrow "Length of row file.",
ncp "No. of compressions of LU data structures.",
mersum "lu1fac: sum of Markowitz merit counts.",
nUtri "lu1fac: triangular rows in U.",
nLtri "lu1fac: triangular rows in L.",
Amax "Maximum element in A.",
Lmax "Maximum multiplier in current L."
Umax "Maximum element in current U."
DUmax "Maximum diagonal in U."
DUmin "Minimum diagonal in U."
Akmax "Maximum element generated at any stage during TCP factorization."
growth "TPP: Umax/Amax TRP, TCP, TSP: Akmax/Amax",
resid "Residual after solve with U or U" /;
parameter fstatus(info) Information returned by LUSOL factorization;
* ------------------------------------------------------------------------------------
* Read the input file:
$call 'gdxdump "lusol.gdx" NoData NoHeader >"%gams.scrdir%lusol_input.scr"'
$include "%gams.scrdir%lusol_input.scr"
* ------------------------------------------------------------------------------------
aa(u1,u2)$a(u1,u2) = a(u1,u2);
option i<aa, j<aa;
* Check that there are an equal number of rows and columns:
abort$(card(i)<>card(j)) "Error in LUSOL. Matrix a is not square:",i,j;
* Generate a set of numerical indices which can be used
* to define numerical indices for rows, columns and rhs vectors:
$if not set lenascale $set lenascale 10
$if not set lena $eval lena max( %lenascale%*card(a), 10000 )
set ia Indices in the factorization /1*%lena%/,
in(ia) Numeric indices for rows and columns;
in(ia)$(ia.val <= card(i)) = yes;
alias (in,jn);
set imap(*,ia) Mapping from row labels to numeric indices ,
jmap(*,ia) Mapping from col labels to numeric indices ;
imap(i,in)$(i.pos = in.val) = yes;
jmap(j,jn)$(j.pos = jn.val) = yes;
option imap:0:0:1, jmap:0:0:1;
display imap, jmap;
$ifthen.b defined b
$set bdim 0
$if dimension 1 b $set bdim 1
$if dimension 2 b $set bdim 2
$if dimension 3 b $set bdim 3
$if dimension 4 b $set bdim 4
$if dimension 5 b $set bdim 5
$if dimension 6 b $set bdim 6
$if dimension 7 b $set bdim 7
$if dimension 8 b $set bdim 8
$if dimension 9 b $set bdim 9
$if %bdim%==0 $abort "LUSOL dimension error: parameter b dimension not between 1 and 9?"
$eval cdim %bdim%-1
$ifthen.cdim not %cdim%==0
alias (cd_1,cd_2,cd_3,cd_4,cd_5,cd_6,cd_7,*);
$if %cdim%==1 $set cd cd_1
$if %cdim%==2 $set cd cd_1,cd_2
$if %cdim%==3 $set cd cd_1,cd_2,cd_3
$if %cdim%==4 $set cd cd_1,cd_2,cd_3,cd_4
$if %cdim%==5 $set cd cd_1,cd_2,cd_3,cd_4,cd_5
$if %cdim%==6 $set cd cd_1,cd_2,cd_3,cd_4,cd_5,cd_6
$if %cdim%==7 $set cd cd_1,cd_2,cd_3,cd_4,cd_5,cd_6,cd_7
set k(%cd%) Dimensions of RHS vector b;
parameter bb(u,%cd%) Parameter used to extract domain of b;
bb(u,%cd%)$b(u,%cd%) = b(u,%cd%);
option k<bb, ib<bb;
noti(ib) = yes$(not i(ib));
abort$card(noti) "LUSOL error: rows in b are not found in a.",i,ib,noti;
set kn(ia) Numeric indices of the RHS;
kn(ia)$(ia.val <= card(k)) = yes;
set kmap(%cd%,ia) Mapping from rhs vectors to numbers
parameter kcount Counter for rhs vectors /0/;
loop(k, kcount = kcount + 1;
kmap(k,kn) = yes$(kn.val = kcount););
option kmap:0:0:1;
display kmap;
$else.cdim
* b is a column vector. Check that all nonzeros are
* from rows in matrix a:
parameter bb(u) Parameter used to extract domain of b;
bb(u)$b(u) = b(u);
option ib<bb;
noti(ib) = yes$(not i(ib));
abort$card(noti) "LUSOL error: rows in b are not found in a.",i,ib,noti;
$endif.cdim
$endif.b
* ------------------------------------------------------------------------------------
parameter luopt(luoptions) LUSOL Options;
$ifthen.luopt set luopt
$gdxin '%luopt%'
$loaddc luopt
$else.luopt
luopt(luoptions)=0;
$endif.luopt
* Assign default LU parameters if we are factorizing:
* Values set for ip and rp correspond to those set in lusol_precision.f90.
* Must have ip=4.
$set rp_huge4 3.40282347E+38
$set rp_huge8 1.79769313486231571E+308
$set rp_huge16 1.18973149535723176508575932662800702E+4932
$set rp_eps4 1.19209290E-07
$set rp_eps8 2.22044604925031308E-016
$set rp_eps16 1.92592994438723585305597794258492732E-0034
$offdigit
$set ip_huge 2147483647
$set rp_eps %rp_eps8%
$set rp_huge %rp_huge8%
luopt("lprint")$(not luopt("lprint")) = 0;
* Print level:
* < 0 suppresses output.
* = 0 gives error messages.
* >= 10 gives statistics about the LU factors.
* >= 50 gives debug output from lu1fac
* (the pivot row and column and the
* no. of rows and columns involved at
* each elimination step).
luopt("maxcol")$(not luopt("maxcol")) = 5;
* Maximum number of columns searched:
* Applies for a Markowitz-type
* search for the next pivot element.
* For some of the factorization, the
* number of rows searched is
* maxrow = maxcol - 1.
$if not set keeplu $set keeplu 1
luopt("keepLU")= %keeplu%;
* keepLU=1 means the numerical factors will be computed if possible.
* keepLU=0 means L and U will be discarded but other information
* such as the row and column permutations will be returned.
* keepLU=0 requires less storage.
luopt("Ltol1")$(not luopt("Ltol1")) = 10;
* Max Lij allowed during Factor:
* TPP 10.0 or 100.0
* TRP 4.0 or 10.0
* TCP 5.0 or 10.0
* TSP 4.0 or 10.0
* With TRP and TCP (Rook and Complete Pivoting), values less than 25.0
* may be expensive on badly scaled data. However, values less than 10.0
* may be needed to obtain a reliable rank-revealing factorization.
luopt("Ltol2")$(not luopt("Ltol2")) = 10;
* Max Lij allowed during Updates.
luopt("small")$(not luopt("small")) = %rp_eps%**0.8;
* Absolute tolerance for treating reals as zero.
luopt("Utol1")$(not luopt("Utol1")) = %rp_eps%**0.67;
* Absolute tol for flagging small diagonals of U.
luopt("Utol2")$(not luopt("Utol2")) = %rp_eps%**0.67;
* Relative tol for flagging small diagonals of U.
luopt("Uspace")$(not luopt("Uspace")) = 3;
* Factor limiting waste space in U.
* In lu1fac, the row or column lists are compressed if their length
* exceeds Uspace times the length of either file after the last
* compression.
luopt("dens1")$(not luopt("dens1")) = 0.3;
* The density at which the Markowitz should reach maxcol columns
* Use 0.3 unless you are experimenting with the pivot strategy.)
luopt("dens2")$(not luopt("dens2")) = 0.1;
* the density at which the Markowitz strategy should search only 1
* column, or (if storage is available) the density at which all
* remaining rows and columns will be processed by a dense LU code. For
* example, if dens2 = 0.1 and lena is large enough, a dense LU will be
* used once more than 10 per cent of the remaining matrix is nonzero.
luopt("PivotScheme")$(not luopt("PivotScheme")) = 0;
* 0 TPP Threshold Partial Pivoting.
* 1 TRP Threshold Rook Pivoting.
* 2 TCP Threshold Complete Pivoting.
* 3 TSP Threshold Symmetric Pivoting.
* 4 TDP Threshold Diagonal Pivoting.
* Note: TDP not yet implemented.
* TRP and TCP are more expensive than TPP but more stable and better at
* revealing rank. Take care with setting parmlu(1), especially with TCP.
* NOTE: TSP and TDP are for symmetric matrices that are either definite
* or quasi-definite. TSP is effectively TRP for symmetric matrices.
* TDP is effectively TCP for symmetric matrices.
* Max Lij allowed during factor -- depends on pivoting scheme.
display luopt;
set parmset /1*30/;
parameter luparm(parmset), parmlu(parmset);
luparm("2")$luopt("lprint") = luopt("lprint");
luparm("8")$luopt("keepLU") = luopt("keepLU");
luparm("3")$luopt("maxcol") = luopt("maxcol");
luparm("6")$luopt("PivotScheme") = luopt("PivotScheme");
parmlu("1")$luopt("Ltol1") = luopt("Ltol1");
parmlu("2")$luopt("Ltol2") = luopt("Ltol2");
parmlu("3")$luopt("small") = luopt("small");
parmlu("4")$luopt("Utol1") = luopt("Utol1");
parmlu("5")$luopt("Utol2") = luopt("Utol2");
parmlu("6")$luopt("Uspace") = luopt("Uspace");
parmlu("7")$luopt("dens1") = luopt("dens1");
parmlu("8")$luopt("dens2") = luopt("dens2");
* ------------------------------------------------------------------------------------
* Only a is defined -- perform a factor:
$if not defined b $goto factor
* Solve with existing factors:
$if set lu $goto solve
* Both a and b are in the input file. Perform a factor / solve:
$goto factorsolve
* ------------------------------------------------------------------------------------
* Factorize the linear system (no rhs vector provided):
$label factor
file kdat /"%gams.scrdir%lusol_factor.scr"/;
put kdat;
kdat.lw=0; kdat.nz = 1.0e-13;
kdat.nw=0; kdat.nd=0;
put 'factor'/;
put '%gams.scrdir%lusol_factor_out.scr'/;
put card(i)/;
put card(a)/;
put '%lena%'/;
kdat.nr=2; kdat.nw=20; kdat.nd=12;
loop((i,j)$a(i,j), loop((imap(i,in),jmap(j,jn)),
put in.tl,' ',jn.tl,a(i,j)/;));
loop(parmset,
kdat.nr=0; kdat.nw=0; kdat.nd=0; put luparm(parmset);
kdat.nr=2; kdat.nw=20; kdat.nd=12; put parmlu(parmset)/;);
putclose;
* Generate the factorization:
execute 'lusol "%gams.scrdir%lusol_factor.scr"';
execute 'gams "%gams.scrdir%lusol_factor_out.scr" gdx="%gams.scrdir%lusol_factor_gdx.scr" o="%gams.scrdir%lusol_factor_lst.scr"';
execute_load '%gams.scrdir%lusol_factor_gdx.scr',fstatus;
execute_unload 'lusol.gdx',a,i,j,fstatus;
$exit
* ------------------------------------------------------------------------------------
* Factorize the linear system and solve:
$label factorsolve
file kdat /"%gams.scrdir%lusol_factorsolve.scr"/;
put kdat;
kdat.lw=0; kdat.nz = 1.0e-13;
kdat.nw=0; kdat.nd=0;
put 'factorsolve'/;
put '%gams.scrdir%lusol_factorsolve_out.scr'/;
put card(i)/;
put card(a)/;
put '%lena%'/;
$if defined k put card(kn)/;
$if not defined k put '1'/;
put card(b)/;
kdat.nr=2; kdat.nw=20; kdat.nd=12;
loop((i,j)$a(i,j), loop((imap(i,in),jmap(j,jn)),
put in.tl,' ',jn.tl,a(i,j)/;));
loop(parmset,
kdat.nr=0; kdat.nw=0; kdat.nd=0; put luparm(parmset);
kdat.nr=2; kdat.nw=20; kdat.nd=12; put parmlu(parmset)/;);
$ifthen defined k
loop((imap(i,in),kmap(k,kn))$b(i,k), put in.tl,' ',kn.tl, b(i,k)/;);
$else
loop(imap(i,in)$b(i), put in.tl,' 1', b(i)/;);
$endif
putclose;
execute 'lusol "%gams.scrdir%lusol_factorsolve.scr"';
execute 'gams "%gams.scrdir%lusol_factorsolve_out.scr" gdx="%gams.scrdir%lusol_factorsolve_gdx.scr" o="%gams.scrdir%lusol_factorsolve_lst.scr"';
parameter xn(ia,ia) Solution values returned by LUSOL,
sstatusn(info,ia) Solution status returned by LUSOL;
execute_load '%gams.scrdir%lusol_factorsolve_gdx.scr',fstatus,sstatusn,xn;
$ifthen defined k
parameter x(*,%cd%) Indexed solution value
sstatus(info,%cd%) Indexed solution status;
loop((jmap(j,jn),kmap(k,kn)),
x(j,k)= xn(jn,kn);
sstatus(info,k) = sstatusn(info,kn););
execute_unload 'lusol.gdx',a,b,i,j,k,fstatus,sstatus,x;
$else
parameter x(j) Indexed solution value
sstatus(info) Indexed solution status;
loop(jmap(j,jn),
x(j) = xn(jn,"1");
sstatus(info) = sstatusn(info,"1"););
execute_unload 'lusol.gdx',a,b,i,j,fstatus,sstatus,x;
$endif
$exit
* ------------------------------------------------------------------------------------
$label solve
file kdat /"%gams.scrdir%lusol_solve.scr"/;
put kdat;
kdat.lw=0; kdat.nz = 1.0e-13;
kdat.nw=0; kdat.nd=0;
put 'solve'/;
put '%lu%'/;
put '%gams.scrdir%lusol_solve_out.scr'/;
put card(i)/;
put card(a)/;
put '%lena%'/;
$if defined k put card(kn)/;
$if not defined k put '1'/;
put card(b)/;
kdat.nr=2; kdat.nw=20; kdat.nd=12;
$ifthen defined k
loop((imap(i,in),kmap(k,kn))$b(i,k), put in.tl,' ',kn.tl, b(i,k)/;);
$else
loop(imap(i,in)$b(i), put in.tl,' 1', b(i)/;);
$endif
putclose;
execute 'lusol "%gams.scrdir%lusol_solve.scr"';
execute 'gams "%gams.scrdir%lusol_solve_out.scr" gdx="%gams.scrdir%lusol_solve_gdx.scr" o="%gams.scrdir%lusol_solve_lst.scr"';
parameter xn(ia,ia) Solution values returned by LUSOL,
sstatusn(info,ia) Solution status returned by LUSOL;
execute_load '%gams.scrdir%lusol_solve_gdx.scr',sstatusn,xn;
$ifthen defined k
parameter x(*,%cd%) Indexed solution value
sstatus(info,%cd%) Indexed solution status;
loop((jmap(j,jn),kmap(k,kn)),
x(j,k)= xn(jn,kn);
sstatus(info,k) = sstatusn(info,kn););
execute_unload 'lusol.gdx',a,b,i,j,k,sstatus,x;
$else
parameter x(j) Indexed solution value
sstatus(info) Indexed solution status;
loop(jmap(j,jn),
x(j) = xn(jn,"1");
sstatus(info) = sstatusn(info,"1"););
execute_unload 'lusol.gdx',a,b,i,j,sstatus,x;
$endif
program main
use lusol_precision, only : ip, rp
use lusol
implicit none
integer(ip) :: i, j, k, n, nelem, nb, nbnz, lena, inform, nindc, nindr, nalu
integer(ip) :: luparm(30), mode
character (len=256) :: inpfile,outfile,lufile,operation
real(rp) :: parmlu(30), small
integer(ip), allocatable :: indc(:) , indr(:), &
p(:) , q(:), &
lenc(:) , lenr(:) , &
iploc(:) , iqloc(:) , &
ipinv(:) , iqinv(:) , &
locc(:) , locr(:), luparmval(:,:)
real(rp), allocatable :: a(:), w(:), x(:,:), v(:), b(:,:), parmluval(:,:)
integer(ip) :: kinp =10, kgms=11, kbin=12
if (command_argument_count().eq.0) stop 'No input file.'
call getarg(1,inpfile)
kinp = 10
kbin = 11
kgms = 12
open(kinp, file=inpfile)
read(kinp,'(a)') operation
if (operation.eq.'factorsolve') call factorsolve
if (operation.eq.'factor') call factor
if (operation.eq.'solve') call solve
contains
subroutine factorsolve
read(kinp,'(a)') outfile
read(kinp,*) n
read(kinp,*) nelem
read(kinp,*) lena
read(kinp,*) nb
read(kinp,*) nbnz
call reada
call readparm
call readb
close(kinp)
call allocatelu
open(kgms,file=outfile)
write(kgms,'(a)') '$ontext'
write(kgms,'(a)') 'LU1FAC LOG:'
luparm(1) = kgms
call lu1fac(n, n, nelem, lena, luparm, parmlu, &
a , indc , indr , p , q , &
lenc , lenr , locc , locr , &
iploc, iqloc, ipinv, iqinv, w , inform )
write(*,'(4(a,i0))') 'LUSOL: n=',n,', nelem=',nelem,', lena=',lena,' -- INFORM (lu1fac) = ',inform
write(kgms,'(a)') '$offtext'
call writefstatus
call writelu
write(kgms,'(//,a)') '$ontext'
write(kgms,'(a)') 'LU6SOL LOG:'
allocate (x(1:n,1:nb))
allocate (parmluval(1:30,1:nb))
allocate (luparmval(1:30,1:nb))
do k=1,nb
v(1:n) = b(1:n,k)
w(1:n) = 0.0
mode = 5
call lu6sol( mode, n, n, v, w, &
lena, luparm, parmlu, &
a, indc, indr, p, q, &
lenc, lenr, locc, locr, &
inform )
write(*,'(a,i0,a,i0,a,i0)') 'LUSOL: ',k,' of ',nb,' -- INFORM (lu6sol) = ',inform
x(1:n,k) = w(1:n)
luparmval(1:30,k) = luparm(1:30)
parmluval(1:30,k) = parmlu(1:30)
end do
write(kgms,'(a//)') '$offtext'
call writesol
close(kgms)
call deallocate_workspace
return
end subroutine factorsolve
subroutine factor
read(kinp,'(a)') outfile
read(kinp,*) n
read(kinp,*) nelem
read(kinp,*) lena
call reada
call readparm
close(kinp)
call allocatelu
open(kgms,file=outfile)
write(kgms,'(a)') '$ontext'
write(kgms,'(a)') 'LU1FAC LOG:'
luparm(1) = kgms
call lu1fac(n, n, nelem, lena, luparm, parmlu, &
a , indc , indr , p , q , &
lenc , lenr , locc , locr , &
iploc, iqloc, ipinv, iqinv, w , inform )
write(*,'(4(a,i0))') 'LUSOL: n=',n,', nelem=',nelem,', lena=',lena,' -- INFORM (lu1fac) = ',inform
write(kgms,'(a)') '$offtext'
call writefstatus
close(kgms)
call writelu
call deallocate_workspace
return
end subroutine factor
subroutine solve
read(kinp,'(a)') lufile
read(kinp,'(a)') outfile
read(kinp,*) n
read(kinp,*) nelem
read(kinp,*) lena
read(kinp,*) nb
read(kinp,*) nbnz
call readb
close(kinp)
call allocatelu
call readlu
open(kgms,file=outfile)
write(kgms,'(//,a)') '$ontext'
write(kgms,'(a)') 'LU6SOL LOG:'
luparm(1) = kgms
! Workspace for solution values and solution status:
allocate (parmluval(1:30,1:nb))
allocate (luparmval(1:30,1:nb))
allocate (x(1:n,1:nb))
do k=1,nb
v(1:n) = b(1:n,k)
w(1:n) = 0.0d0
mode = 5
call lu6sol( mode, n, n, v, w, &
lena, luparm, parmlu, &
a, indc, indr, p, q, &
lenc, lenr, locc, locr, &
inform )
Write(*,'(A,i0,a,i0,a,i0)') 'LUSOL: ',k,' of ',nb,' -- INFORM (lu6sol) = ',inform
x(1:n,k) = w(1:n)
luparmval(1:30,k) = luparm(1:30)
parmluval(1:30,k) = parmlu(1:30)
end do
write(kgms,'(a//)') '$offtext'
call writesol
close(kgms)
call deallocate_workspace
return
end subroutine solve
subroutine reada
allocate (a(1:lena))
allocate (indc(1:lena))
allocate (indr(1:lena))
do k=1,nelem
read(kinp,*) indc(k), indr(k), a(k)
end do
return
end
subroutine readparm
do k=1,30
read(kinp,*) luparm(k),parmlu(k)
end do
return
end
subroutine readb
allocate (b(1:n,1:nb))
do k=1,nbnz
read(kinp,*) i,j,b(i,j)
enddo
return
end
subroutine allocatelu
allocate (p(1:n))
allocate (q(1:n))
allocate (lenc(1:n))
allocate (lenr(1:n))
allocate (iploc(1:n))
allocate (iqloc(1:n))
allocate (ipinv(1:n))
allocate (iqinv(1:n))
allocate (locc(1:n))
allocate (locr(1:n))
allocate (w(1:n))
allocate (v(1:n))
return
end
subroutine writefstatus
write(kgms,'(a)') 'parameter fstatus(*) /'
write(kgms,'(a,i0)') 'inform ',luparm(10) ! Return code from last call to any LU routine.
write(kgms,'(a,i0)') 'nsing ',luparm(11) ! No. of singularities marked in the output array w(*).
write(kgms,'(a,i0)') 'jsing ',luparm(12) ! Column index of last singularity.
write(kgms,'(a,i0)') 'minlen ',luparm(13) ! Minimum recommended value for lena.
write(kgms,'(a,i0)') 'maxlen ',luparm(14) ! ?
write(kgms,'(a,i0)') 'lena ',lena ! Allocated length
write(kgms,'(a,i0)') 'nelem ',nelem ! Number of elements
write(kgms,'(a,i0)') 'nupdat ',luparm(15) ! No. of updates performed by the lu8 routines.
write(kgms,'(a,i0)') 'nrank ',luparm(16) ! No. of nonempty rows of U.
write(kgms,'(a,i0)') 'ndens1 ',luparm(17) ! No. of columns remaining when the density of
! the matrix being factorized reached dens1.
write(kgms,'(a,i0)') 'ndens2 ',luparm(18) ! No. of columns remaining when the density of
! the matrix being factorized reached dens2.
write(kgms,'(a,i0)') 'jumin ',luparm(19) ! The column index associated with DUmin.
write(kgms,'(a,i0)') 'numL0 ',luparm(20) ! No. of columns in initial L.
write(kgms,'(a,i0)') 'lenL0 ',luparm(21) ! Size of initial L (no. of nonzeros).
write(kgms,'(a,i0)') 'lenU0 ',luparm(22) ! Size of initial U.
write(kgms,'(a,i0)') 'lenL ',luparm(23) ! Size of current L.
write(kgms,'(a,i0)') 'lenU ',luparm(24) ! Size of current U.
write(kgms,'(a,i0)') 'lrow ',luparm(25) ! Length of row file.
write(kgms,'(a,i0)') 'ncp ',luparm(26) ! No. of compressions of LU data structures.
write(kgms,'(a,i0)') 'mersum ',luparm(27) ! lu1fac: sum of Markowitz merit counts.
write(kgms,'(a,i0)') 'nUtri ',luparm(28) ! lu1fac: triangular rows in U.
write(kgms,'(a,i0)') 'nLtri ',luparm(29) ! lu1fac: triangular rows in L.
write(kgms,'(a,1pe20.12)') 'Amax ',parmlu(10) ! Maximum element in A.
write(kgms,'(a,1pe20.12)') 'Lmax ',parmlu(11) ! Maximum multiplier in current L.
write(kgms,'(a,1pe20.12)') 'Umax ',parmlu(12) ! Maximum element in current U.
write(kgms,'(a,1pe20.12)') 'DUmax ',parmlu(13) ! Maximum diagonal in U.
write(kgms,'(a,1pe20.12)') 'DUmin ',parmlu(14) ! Minimum diagonal in U.
write(kgms,'(a,1pe20.12)') 'Akmax ',parmlu(15) ! Maximum element generated at any stage
! during TCP factorization.
write(kgms,'(a,1pe20.12)') 'growth ',parmlu(16) ! TPP: Umax/Amax TRP, TCP, TSP: Akmax/Amax
write(kgms,'(a,1pe20.12)') 'resid ',parmlu(20) ! lu6sol: residual after solve with U or U
write(kgms,'(a)') '/;'
return
end
subroutine readlu
integer(ip) :: n_lu, nelem_lu, lena_lu
open(kbin,file=lufile,form='unformatted',access='sequential')
read(kbin) n_lu
read(kbin) nelem_lu
read(kbin) lena_lu
if (n.ne.n_lu) stop "Error: LU factors are inconsisten: n is different"
if (nelem.ne.nelem_lu) stop "Error: LU factors are inconsistent: nelem is different."
if (lena.ne.lena_lu) stop "Error: LU factors are inconsistent: lena is different."
read(kbin) luparm
read(kbin) parmlu
read(kbin) p
read(kbin) q
read(kbin) lenc
read(kbin) lenr
read(kbin) iploc
read(kbin) iqloc
read(kbin) ipinv
read(kbin) iqinv
read(kbin) locc
read(kbin) locr
! write(*,'(a,1pe20.12)') 'readlu -- parmlu(3) = ',parmlu(3)
allocate (a(1:lena))
allocate (indc(1:lena))
allocate (indr(1:lena))
a(1:lena) = 0.0
indc(1:lena) = 0
indr(1:lena) = 0
read(kbin) nindc,nindr,nalu
do k=1,nindc
read(kbin) i,indc(i)
end do
do k=1,nindr
read(kbin) i, indr(i)
end do
do k=1,nalu
read(kbin) i, a(i)
end do
close(kbin)
return
end subroutine readlu
subroutine writelu
open(kbin,file='lusol.bin',form='unformatted',access='sequential')
write(kbin) n
write(kbin) nelem
write(kbin) lena
write(kbin) luparm
write(kbin) parmlu
write(kbin) p
write(kbin) q
write(kbin) lenc
write(kbin) lenr
write(kbin) iploc
write(kbin) iqloc
write(kbin) ipinv
write(kbin) iqinv
write(kbin) locc
write(kbin) locr
nindc=0
nindr=0
nalu=0
!write(*,'(a,1pe20.12)') 'writelu -- parmlu(3) = ',parmlu(3)
small = parmlu(3) ! Drop tolerance
do i=1,lena
if (indc(i).ne.0) nindc = nindc + 1
if (indr(i).ne.0) nindr = nindr + 1
if (abs(a(i))>=small) nalu = nalu + 1
end do
write(kbin) nindc,nindr,nalu
do i=1,lena
if (indc(i).ne.0) write(kbin) i,indc(i)
end do
do i=1,lena
if (indr(i).ne.0) write(kbin) i,indr(i)
end do
do i=1,lena
if (abs(a(i))>=small) write(kbin) i,a(i)
end do
close(kbin)
return
end subroutine
subroutine writesol
small = parmlu(3) ! Drop tolerance
write(kgms,'(a)') 'parameter xn(*,*) /'
do i=1,n
do k=1,nb
if (abs(x(i,k))>=small) write(kgms,'(i0,a,i0,1x,1pe22.12)') i,'.',k,x(i,k)
end do
end do
write(kgms,'(a)') '/;'
write(kgms,'(a)') 'parameter sstatusn(*,*) /'
do k=1,nb
write(kgms,'(a,i0,1x,i0)') 'inform.',k,luparmval(10,k) ! Return code from last call to any LU routine.
write(kgms,'(a,i0,1x,i0)') 'nsing.',k,luparmval(11,k) ! No. of singularities marked in the output array w(*).
write(kgms,'(a,i0,1x,i0)') 'jsing.',k,luparmval(12,k) ! Column index of last singularity.
write(kgms,'(a,i0,1x,i0)') 'minlen.',k,luparmval(13,k) ! Minimum recommended value for lena.
write(kgms,'(a,i0,1x,i0)') 'maxlen.',k,luparmval(14,k) ! ?
write(kgms,'(a,i0,1x,i0)') 'lena.',k,lena ! Allocated length
write(kgms,'(a,i0,1x,i0)') 'nupdat.',k,luparmval(15,k) ! No. of updates performed by the lu8 routines.
write(kgms,'(a,i0,1x,i0)') 'nrank.',k,luparmval(16,k) ! No. of nonempty rows of U.
write(kgms,'(a,i0,1x,i0)') 'ndens1.',k,luparmval(17,k) ! No. of columns remaining when the density of
write(kgms,'(a,i0,1x,i0)') 'ndens2.',k,luparmval(18,k) ! No. of columns remaining when the density of
write(kgms,'(a,i0,1x,i0)') 'jumin.',k,luparmval(19,k) ! The column index associated with DUmin.
write(kgms,'(a,i0,1x,i0)') 'numL0.',k,luparmval(20,k) ! No. of columns in initial L.
write(kgms,'(a,i0,1x,i0)') 'lenL0.',k,luparmval(21,k) ! Size of initial L (no. of nonzeros).
write(kgms,'(a,i0,1x,i0)') 'lenU0.',k,luparmval(22,k) ! Size of initial U.
write(kgms,'(a,i0,1x,i0)') 'lenL.',k,luparmval(23,k) ! Size of current L.
write(kgms,'(a,i0,1x,i0)') 'lenU.',k,luparmval(24,k) ! Size of current U.
write(kgms,'(a,i0,1x,i0)') 'lrow.',k,luparmval(25,k) ! Length of row file.
write(kgms,'(a,i0,1x,i0)') 'ncp.',k,luparmval(26,k) ! No. of compressions of LU data structures.
write(kgms,'(a,i0,1x,i0)') 'mersum.',k,luparmval(27,k) ! lu1fac: sum of Markowitz merit counts.
write(kgms,'(a,i0,1x,i0)') 'nUtri.',k,luparmval(28,k) ! lu1fac: triangular rows in U.
write(kgms,'(a,i0,1x,i0)') 'nLtri.',k,luparmval(29,k) ! lu1fac: triangular rows in L.
write(kgms,'(a,i0,1x,1pe20.12)') 'Amax.',k,parmluval(10,k) ! Maximum element in A.
write(kgms,'(a,i0,1x,1pe20.12)') 'Lmax.',k,parmluval(11,k) ! Maximum multiplier in current L.
write(kgms,'(a,i0,1x,1pe20.12)') 'Umax.',k,parmluval(12,k) ! Maximum element in current U.
write(kgms,'(a,i0,1x,1pe20.12)') 'DUmax.',k,parmluval(13,k) ! Maximum diagonal in U.
write(kgms,'(a,i0,1x,1pe20.12)') 'DUmin.',k,parmluval(14,k) ! Minimum diagonal in U.
write(kgms,'(a,i0,1x,1pe20.12)') 'Akmax.',k,parmluval(15,k) ! Maximum element generated at any stage
write(kgms,'(a,i0,1x,1pe20.12)') 'growth.',k,parmluval(16,k) ! TPP: Umax/Amax TRP, TCP, TSP: Akmax/Amax
write(kgms,'(a,i0,1x,1pe20.12)') 'resid.',k,parmluval(20,k) ! lu6sol: residual after solve with U or U
end do
write(kgms,'(a)') '/;'
return
end
subroutine deallocate_workspace
deallocate (p)
deallocate (q)
deallocate (lenc)
deallocate (lenr)
deallocate (iploc)
deallocate (iqloc)
deallocate (ipinv)
deallocate (iqinv)
deallocate (locc)
deallocate (locr)
deallocate (w)
deallocate (v)
deallocate (a)
deallocate (indc)
deallocate (indr)
if (allocated(parmluval)) deallocate (parmluval)
if (allocated(luparmval)) deallocate (luparmval)
if (allocated(b)) deallocate (b)
if (allocated(x)) deallocate (x)
return
end
end program main
$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.
* 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 /'* Pass along user provided arguments:'/;
put '$if not set cmdline $set cmdline'/;
put /'* Generate the command line argument list:'/;
put '$setlocal ide ide=%','gams.ide% lo=%','gams.lo% errorlog=%','gams.errorlog% errmsg=1 %cmdline%'/;
put /'* Provide flags for optional arguments:'/;
loop(interface("opt",id,type,dim,ininout),
loop(optarg,
put '$if not "%',optarg.tl,'id%"=="" $setlocal ide %','ide% --%',optarg.tl,'id%=%',optarg.tl,'%'/;));
put /'* Execute the program:'/;
put '$log --- Executing ',"'",'=gams "%input%" %','ide%',"'"/;
put 'execute ',"'",'=gams "%input%" %','ide%',"';"/;
put 'abort$errorlevel ',"'Abort: call to %input% 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,'%'/;));
if (sum(interface(argopt,u,type,dim,outinout),1),
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 "%arg3%"=="" $abort Error in invocation -- too many arguments: %args%
* Verify consistency of required arguments:
$ifthen"%arg1%"==""
$abort Error in lusol call: missing argument arg1 = A.
$endif
$if not partype %arg1% $abort Error in lusol: A=%arg1% is not 2-dimensional par
$if not dimension 2 %arg1% $abort Error in lusol: A=%arg1% is not 2-dimensional par
* Find ids to which optional arguments have been assigned:
$setlocal arg2id
$ifthen.arg2 not "%arg2%"==""
$ifthen.btypearg2 partype %arg2%
$setlocal arg2id b
$endif.btypearg2
$ifthen.emptyarg2 "%arg2%"==""
$abort " Error in lusol: argument arg2=%arg2% does not match optional arguments: b (* dim par)
$endif.emptyarg2
$endif.arg2
* Verify that all inputs have been defined:
$if not defined %arg1% $abort "Input item not defined: %arg1% (= A)"
$ifthen.defined_arg2 "%arg2id%"=="b"
$if not defined %arg2% $abort "Input item not defined: %arg2% (= b)"
$endif.defined_arg2
* Write a GDX file with the program inputs:
$setlocal inputs %arg1%=A
$if "%arg2id%"=="b" $setlocal inputs %inputs%, %arg2%=b
execute_unload 'lusol.gdx', %inputs%;
$log Entering lusol with %inputs%
* Pass along user provided arguments:
$if not set cmdline $set cmdline
* Generate the command line argument list:
$setlocal ide ide=%gams.ide% lo=%gams.lo% errorlog=%gams.errorlog% errmsg=1 %cmdline%
* Provide flags for optional arguments:
$if not "%arg2id%"=="" $setlocal ide %ide% --%arg2id%=%arg2%
* Execute the program:
$log --- Executing '=gams "d:\svn.mpsge.org\lsa\lusol\distrib\batinclude\lusol.gms" %ide%'
execute '=gams "d:\svn.mpsge.org\lsa\lusol\distrib\batinclude\lusol.gms" %ide%';
abort$errorlevel 'Abort: call to d:\svn.mpsge.org\lsa\lusol\distrib\batinclude\lusol.gms has failed. See %gams.wdir%lusol.lst.';