Here is the model:
$ontext
 
   Invert the Pei test matrix
 
   Erwin Kalvelagen, july 2008
 
   Command line arguments:
      --n=integer   order of the matrix
      --alpha=real  scalar defining the Pei matrix
      --method=1    invert using equations A AINV = I
      --method=2    invert using equations A AINV = I using advanced basis
      --method=3    invert using external solver
 
 
   Reference:
 
    Erwin Kalvelagen
    SOLVING SYSTEMS OF LINEAR EQUATIONS WITH GAMS
    http://www.amsterdamoptimization.com/pdf/lineq.pdf
  
    ML Pei,
    A test matrix for inversion procedures,
    Communications of the ACM,
    Volume 5, 1962, page 508.
 
$offtext
 
*
* defaults
*
$if not set n $set n 100
$if not set method $set method 1
$if not set alpha $set alpha 1
 
scalar alpha /%alpha%/;
 
set i /i1*i%n%/;
alias (i,j,k);
 
parameter a(i,j);
a(i,j)=1;
a(i,i)=1+alpha;
 
parameter ident(i,j);
ident(i,i)=1;
 
variables
   ainv(i,j)
   dummy
;
 
equations
   e(i,j)
   edummy
;
 
edummy.. dummy =e= 0;
e(i,j)..  sum(k, a(i,k)*ainv(k,j)) =e= ident(i,j);
 
model m /all/;
 
option iterlim=1000000;
 
 
 
$if "%method%"==1 solve m minimizing dummy using lp;
 
$ifthen "%method%"==2
 
edummy.m=1;
e.m(i,j)=1;
dummy.m=0;
ainv.m(i,j)=0;
solve m minimizing dummy using lp;
 
$endif
 
$ifthen "%method%"==3
 
execute_unload 'a.gdx',i,a;
execute '=invert.exe a.gdx i a b.gdx pinva';
parameter pinva(i,j);
execute_load 'b.gdx',pinva;
 
$endif
We have implemented three methods of calculating the inverse.
- Method=1 uses the LP solver to solve A * invA = I.
- Method=2 is as method=1 but we use an advanced basis to speed up the solver.
- Method=3 is calling the external solver INVERT.EXE
For more information see: http://www.amsterdamoptimization.com/pdf/lineq.pdf.
Some timings are listed below (all times in seconds):
             method=1   method=2   method=3
   n=50       0.637       0.433      0.027
   n=100      8.267       4.036      0.055
   n=200    313.236      53.395      0.118
The timings are taken from the last line in the log file. E.g.:
--- Job Untitled_36.gms Start 07/19/08 15:46:36
GAMS Rev 227  Copyright (C) 1987-2008 GAMS Development. All rights reserved
Licensee: Erwin Kalvelagen                               G070509/0001CE-WIN
          GAMS Development Corporation                               DC4572
--- Starting compilation
--- Untitled_36.gms(81) 3 Mb
--- Starting execution: elapsed 0:00:00.004
--- Untitled_36.gms(69) 5 Mb
 ----------------------------------------
 INVERT: matrix inversion
 Erwin Kalvelagen, Amsterdam Optimization
 ----------------------------------------
 DLL:GDX Library      May  1, 2008 22.7.2 WIN 4701.4799 VIS x86/MS Windows
 GDXin:a.gdx
 Input set:i
 Input parameter:a
 LAPACK DGESV, n=         200
 GDXout:b.gdx
 Output parameter:pinva
 Done
--- Untitled_36.gms(71) 5 Mb
--- GDXin=C:\projects\test\b.gdx
--- Untitled_36.gms(71) 6 Mb
*** Status: Normal completion
--- Job Untitled_36.gms Stop 07/19/08 15:46:36 elapsed 0:00:00.118
 
 
No comments:
Post a Comment