$TITLE: M5-5.GMS: invert a matrix using mcp;

$ontext
Interestingly, inverting a matrix can be converted into an MCP problem.
An the method is extremely sparse in coding and very fast: a 500x500
matrix can be inverted in a couple of second.

But first, we show a simple 2x2 problem written out in full to show
the key to the majic.  Then we illustrate the efficient version on
a 3x3.  Thanks to Tom Rutherford and Edward Balistreri for the latter.
$offtext

SETS R /R1*R2/
     C /C1*C2/;


TABLE M(R,C)
       C1    C2
R1      2     1
R2      1     2;

VARIABLES
 MINV11  element 11 of the inverse of M
 MINV21  element 21 of the inverse of M
 MINV12  element 12 of the inverse of M
 MINV22  element 22 of the inverse of M;

EQUATIONS
 EL11
 EL21
 EL12
 EL22;

* note: first two equations solve for two unknowns: MINV11, MINV21

EL11.. M("R1","C1")*MINV11 + M("R1","C2")*MINV21 =E= 1;
EL21.. M("R2","C1")*MINV11 + M("R2","C2")*MINV21 =E= 0;

* note: second two equations for for two unknowns: MINV12, MINV22

EL12.. M("R1","C1")*MINV12 + M("R1","C2")*MINV22 =E= 0;
EL22.. M("R2","C1")*MINV12 + M("R2","C2")*MINV22 =E= 1;

MODEL INVERSE /EL11.MINV11, EL21.MINV21, EL12.MINV12, EL22.MINV22/;
SOLVE INVERSE USING MCP;



* this specific example is useful in understanding the following general
* method: solve nxn sub-problems for each column of the inverse matrix


SETS I      row index       /1*3/
     N(I)   active row;

ALIAS (I,J,K);

TABLE A(I,J) matrix to be inverted
        1       2        3
1       4       1       -1
2       0       3        2
3       3       0        7;

PARAMETERS
 IM(I,J) identity matrix
 B(I,J)  inverse of A;

IM(I,I)= 1;

VARIABLE
 X(I) current solution column of B;

EQUATION
 INV(I) definition of inverse on column I of B;

INV(I).. SUM(K, A(I,K)*X(K)) - 1$N(I) =E= 0;

MODEL INVERT /INV.X/;

LOOP(J,

N(I) = YES$(ORD(I) eq ORD(J));

SOLVE INVERT USING MCP;

B(I,J) = X.L(I);

);


* check that we have the inverse
* also shows how to do matrix multiplication

PARAMETER
 VERIFY(I,J)  A times B: should be a matrix of zeros;

*VERIFY(I,J) = SUM(K, A(I,K)*B(K,J)) - IM(I,J);
VERIFY(I,J) = SUM(K, A(I,K)*B(K,J));

DISPLAY A, B, VERIFY;