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