$title Warm-starting quad-precision MINOS (DQQ,SEQ=414)

$onText
Illustrate the DQQ procedure. This procedure is used for multiscale
LP models (i.e. models with coefficients and/or variables that vary
over an extreme range). It takes advantage of the quad-precision
solver QUADMINOS to handle the extreme range of values encountered.

In the reference above, several instances are solved by DQQ: the data
for some of these instances can be processed here by specifying the
GDX input when running, e.g.

  > gams dqq --GDXNAME=pilot

We include data for the following instances (in order of increasing difficulty):
  pilot4, pilot, pilot87, tma_me, glcaerwt, glcalift
These instances were obtained from MPS files via the mps2gms utility.

Ma, D., Yang, L., Fleming, R.M.T., Thiele, I., Palsson, B.O.,
Saunders, M.A.  Reliable and efficient solution of genome-scale models
of Metabolism and macromolecular Expression, Scientific Reports 7,
Article number 40863 (2017).
doi:10.1038/srep40863. http://rdcu.be/oCpn

Keywords: linear programming, warm starting, multiscale model, GAMS language features,
          quadruple-precision floating-point arithmetic, QUADMINOS
$offText

$if not set GDXNAME $set GDXNAME pilot4

$ifThen not set LPNAME
$  set LPNAME MINOS
$  set OPTNUM 11
$else
$  set OPTNUM 0
$endIf

Set
   i     'all rows in MPS order'
   ig(i) 'greater-than-or equal rows'
   il(i) 'less-than-or equal rows'
   ie(i) 'equality rows'
   ir(i) 'ranged rows';

Set
   j     'all columns in MPS order'
   jc(j) 'continuous columns';

Parameter
   c(j)     'objective coefs'
   cobj     'objective constant'
   b(i)     'right hand sides'
   ac(i,jc) 'matrix coefs: continuous variables';

Set mps2gms;

Parameter mps2gmsstats(mps2gms);

Equation
   eobj  'objective function'
   eg(i) 'greater-than-or equal equs'
   el(i) 'less-than-or equal equs'
   ee(i) 'equality equs'
   er(i) 'ranged equs';

Variable obj 'objective variable';

Positive Variable
   xc(j) 'continuous variables'
   r(i)  'ranged row variables';

$gdxIn %GDXNAME%
$load i j mps2gms mps2gmsstats
$load ig il ie ir
$load jc
$load cobj c b
$load ac
$load xc r
$gdxIn

eobj..   obj =e= sum(jc, c(jc)*xc(jc)) + cobj;

eg(ig).. sum(jc, ac(ig,jc)*xc(jc)) =g= b(ig);

el(il).. sum(jc, ac(il,jc)*xc(jc)) =l= b(il);

ee(ie).. sum(jc, ac(ie,jc)*xc(jc)) =e= b(ie);

er(ir).. sum(jc, ac(ir,jc)*xc(jc)) =e= r(ir);

Model m / all /;

option limCol = 0, limRow = 0, solPrint = off;

$onText
Implement DQQ procedure:
 1. Solve with Double-precision MINOS, scaling on,  optFile = 11
 2. Solve with Quad-precision   MINOS, scaling on,  optFile = 12
 3. Solve with Quad-precision   MINOS, scaling off, optFile = 13
$offText

$onEcho > minos.o11
* step 1 (D) of DQQ
scale option             2
feasibility tolerance    1e-7
optimality tolerance     1e-7
expand frequency         100000
LU factor tolerance      1.9
LU update tolerance      1.9
$offEcho

$onEcho > quadminos.o12
* step 2 (Q1) of DQQ
scale option             2
feasibility tolerance    1e-15
optimality tolerance     1e-15
expand frequency         100000
LU factor tolerance      10.0
LU update tolerance      10.0
$offEcho

$onEcho > quadminos.o13
* step 3 (Q2) of DQQ
scale option             0
feasibility tolerance    1e-15
optimality tolerance     1e-15
expand frequency         100000
LU factor tolerance      5.0
LU update tolerance      5.0
$offEcho

Parameter xbar1(jc), xbar2(jc), xbar3(jc);

Scalar obj1, itr1, obj2, itr2, dx2 / NA /, obj3, itr3, dx3 / NA /;

* step 1: solve with normal (Double-precision) MINOS
m.optFile = %OPTNUM%;
option lp = %LPNAME%;
solve m using lp minimizing obj;
obj1 = m.objVal;
itr1 = m.iterUsd;
xbar1(jc) = xc.L(jc);

* step 2: solve with quad-precision MINOS
m.optFile = 12;
option lp = quadminos;
solve m using lp minimizing obj;
obj2 = m.objVal;
itr2 = m.iterUsd;
xbar2(jc) = xc.L(jc);

* step 3: solve with quad-precision MINOS, no scaling
m.optFile = 13;
option lp = quadminos;
solve m using lp minimizing obj;
obj3 = m.objVal;
itr3 = m.iterUsd;
xbar3(jc) = xc.L(jc);

dx2 = smax{jc, abs(xbar1(jc) - xbar2(jc))};
dx3 = smax{jc, abs(xbar2(jc) - xbar3(jc))};

File log /''/;
putClose log ' '
   / 'Instance name: %GDXNAME%'
   / 'Step 1:  itr =' itr1:10:0  '   obj =' obj1:18:8
   / 'Step 2:  itr =' itr2:10:0  '   obj =' obj2:18:8 '   dx =' dx2:14:8
   / 'Step 3:  itr =' itr3:10:0  '   obj =' obj3:18:8 '   dx =' dx3:14:8
   /;