Wednesday, June 18, 2008

Invert: An External Matrix Inversion program for GAMS

Matrix Inversion


It is possible to perform a matrix inversion using a standard LP solver. See http://www.amsterdamoptimization.com/pdf/lineq.pdf for a detailed description and examples of how this can be achieved. Sometimes it is quicker however just to call an external program to calculate the inverse of a square matrix. The program uses LAPACK's DGESV subroutine.



Download


Attach:invert.exe



Usage


The command line parameters are explained when running invert without parameters:




C:\>invert
----------------------------------------
INVERT: matrix inversion
Erwin Kalvelagen, Amsterdam Optimization
----------------------------------------
Usage
> invert gdxin i a gdxout inva
where
gdxin : name of gdxfile with matrix
i : name of set used in matrix
a : name of 2 dimensional parameter inside gdxin
gdxout : name of gdxfile for results (inverse matrix)
inva : name of 2 dimensional parameter inside gdxout

Calculates the inverse of a non-singular matrix a(i,j) where
i and j are aliased sets. inva will contain the inverse
matrix.


C:\>



Example 1


The input looks like:




$ontext

Finds the inverse of a matrix through an external program

Erwin Kalvelagen, march 2005

Reference: model gauss.gms from the model library
http://www.gams.com/modlib/libhtml/gauss.htm

$offtext

set i /i1*i3 /;
alias (i,j);

table a(i,j) 'original matrix'
i1 i2 i3
i1 1 2 3
i2 1 3 4
i3 1 4 3
;

parameter inva(i,j) 'inverse of a';

execute_unload 'a.gdx',i,a;
execute '=invert.exe a.gdx i a b.gdx inva';
execute_load 'b.gdx',inva;

display a,inva;


The program has an updated GDXIO interface that both works with old 22.5 GAMS and earlier and with new GAMS 22.6 and later. The GAMS people broke the GDX interface so these versions are really incompatible. However, some of my clients were using a combination of older and newer GAMS systems, so I had to spend lots of time to write a GDX interface that both work with all these GAMS systems. Here is a log from a run with 22.6:




--- Job invert1.gms Start 02/22/08 13:29:46
GAMS Rev 149 Copyright (C) 1987-2007 GAMS Development. All rights reserved
Licensee: Erwin Kalvelagen G070509/0001CE-WIN
GAMS Development Corporation DC4572
--- Starting compilation
--- invert1.gms(28) 3 Mb
--- Starting execution: elapsed 0:00:00.016
--- invert1.gms(25) 4 Mb
----------------------------------------
INVERT: matrix inversion
Erwin Kalvelagen, Amsterdam Optimization
----------------------------------------
DLL:GDX Library Dec 24, 2007 WIN.FR.NA 22.6 239.000.000.vis P3PC
GDXin:a.gdx
Input set:i
Input parameter:a
LAPACK DGESV, n= 3
GDXout:b.gdx
Output parameter:inva
Done
--- invert1.gms(26) 4 Mb
--- GDXin=c:\tmp\New Folder\b.gdx
--- invert1.gms(28) 4 Mb
*** Status: Normal completion
--- Job invert1.gms Stop 02/22/08 13:29:46 elapsed 0:00:00.063




Example 2

$ontext

Finds the inverse of a matrix through an external program

Erwin Kalvelagen, march 2005

Reference: model gauss.gms from the model library
http://www.gams.com/modlib/libhtml/gauss.htm

$offtext

set i /i1*i5/;
alias (i,j);

parameter a(i,j) 'Hilbert matrix';
a(i,j) = 1/(ord(i)+ord(j)-1);

parameter inva(i,j) 'inverse of a';

execute_unload 'a.gdx',i,a;
execute '=invert.exe a.gdx i a b.gdx inva';
execute_load 'b.gdx',inva;

display a,inva;



The output with GAMS 22.5 is:




--- Job invert2.gms Start 02/22/08 13:32:23
GAMS Rev 148 Copyright (C) 1987-2007 GAMS Development. All rights reserved
Licensee: Erwin Kalvelagen G070509/0001CE-WIN
GAMS Development Corporation DC4572
--- Starting compilation
--- invert2.gms(24) 3 Mb
--- Starting execution
--- invert2.gms(21) 4 Mb
----------------------------------------
INVERT: matrix inversion
Erwin Kalvelagen, Amsterdam Optimization
----------------------------------------
DLL:_GAMS_GDX_237_2007-01-09
GDXin:a.gdx
Input set:i
Input parameter:a
LAPACK DGESV, n= 5
GDXout:b.gdx
Output parameter:inva
Done
--- invert2.gms(22) 4 Mb
--- GDXin=c:\tmp\New Folder\b.gdx
--- invert2.gms(24) 4 Mb
*** Status: Normal completion
--- Job invert2.gms Stop 02/22/08 13:32:23 elapsed 0:00:00.078



Notes
  • For Lapack info see: http://www.netlib.org/lapack/
  • The program can efficiently invert matrices up to (1000 x 1000). For bigger matrices, the dense inversion routine becomes slower and memory requirements become large.
  • Source: Main program Attach:invert.f90



No comments:

Post a Comment