gauss.gms : Matrix Inversion with Full Pivoting

Description

This example demonstrates the use of Loops and Dynamic definition
of sets in elementary transformations using Gaussian Elimination with
full pivot selection.


Reference

  • GAMS Development Corporation, Formulation and Language Example.

Small Model of Type : GAMS


Category : GAMS Model library


Main file : gauss.gms

$Title Matrix Inversion with full Pivoting (GAUSS,SEQ=71)

$Ontext

   This example demonstrates the use of Loops and Dynamic definition
   of sets in elementary transformations using Gaussian Elimination with
   full pivot selection.


GAMS Development Corporation, Formulation and Language Example.

$Offtext

Sets  i  row labels    / row-1 * row-3 /
      j  column labels / col-1 * col-3 / ;  alias (i,ip), (j,jp)

Table a(i,j) original matrix

       col-1  col-2  col-3
row-1    1      2      3
row-2    1      3      4
row-3    1      4      3

Parameters  b(j,i)    inverse of a
            bp(i,j)   permuted and transposed inverse of a
            pair(i,j) pivoting sequence and permutation
            rank      rank of matrix a
            adet      absolute value of determinant of matrix a

            piv, big, tol ;

Sets  r(i)  pivot row candidates   , npr(i)  non pivot rows
      s(j)  pivot column candidates, nps(j)  non pivot columns ;

r(i) = yes; s(j) = yes; bp(i,j) = a(i,j); rank = 0;
adet = 1; tol = 1e-5;

Loop(j, big = smax((r,s), abs(bp(r,s))); big$(big lt tol) = 0;
        npr(i) = yes; nps(jp) = yes;

        Loop((r,s)$(big and big eq abs(bp(r,s))),

                rank = rank+1; pair(r,s) = rank; piv = 1/bp(r,s);
                big = 0; adet = adet/piv; npr(r) = no; nps(s) = no;

                bp(  r,nps)  =  bp(r,nps)*piv;
                bp(npr,nps)  =  bp(npr,nps) - bp(r,nps)*bp(npr,s);
                bp(npr,  s)  = -bp(npr,s)*piv;
                bp(  r,  s)  =  piv             );

        r(r) = npr(r); s(s) = nps(s)                 );

b(j,i) = sum((ip,jp)$(pair(i,jp) and pair(ip,j)), bp(ip,jp));
adet   = abs(adet);

Display a, rank, adet, pair, bp, b;

Parameters ir(i,ip) a times b,  ic(j,jp) b times a;

   ir(i,ip) = sum(j, a(i,j)*b(j,ip));
   ic(j,jp) = sum(i, b(j,i)*a(i,jp));

Display ir, ic;