$title 'Test invert and numpy matrix_rank on rank-deficient inputs' (INVERT02,SEQ=392)

$onText

Test the invert utility and some embedded code Python with the numpy 
algorithm matrix_rank on rank-deficient inputs.

Contributor: Erwin Kalvelagen and Steve Dirkse, July 2008. 
Adjusted by Michael Bussieck, January 2024.

$offText

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

parameter
  A(i,j)
  rankDeficient(i,j)
  inv(i,j)           'inverse matrix'
  chk(i,j)           'check the product'
  rank               'matrix rank'
  ;

A(i,i) = 1;


executeTool.checkErrorLevel 'linalg.invert i A inv';
$onImplicitAssign

chk(i,j) = sum{k, A(i,k)*inv(k,j)};
chk(i,j) = round(chk(i,j),14);
display A,inv,chk;
chk(i,i) = chk(i,i) - 1;
abort$[card(chk)] 'A * inv <> identity';

loop {r,
* create a rank-r matrix from A, and check that we get the right
* return code from invert
  rankDeficient(i,j) = A(i,j)$[ord(j) <= ord(r)];
  embeddedCode Python:
  import gams.transfer as gt
  import numpy as np

  m = gt.Container(gams.db, system_directory=r"%gams.sysdir% ".strip())
  rankDeficient = m.data["rankDeficient"].toDense()
  gams.set('rank', [float(np.linalg.matrix_rank(rankDeficient))])
  endEmbeddedCode rank
  abort$(ord(r) <> rank) 'Bad rank returned from numpy matrix_rank',
    rank, rankDeficient, r;
};