! Computing the rank of a matrix.  (RankFind.lng);
! Keywords: Matrix rank, Rank of Matrix;
SETS:
 dimr: listr;  ! The row dimension;
 dimc: listc;  ! The col dimension;
 rxc( dimr, dimc): A, WORK;
ENDSETS
DATA:
  tolpiv = 1e-14; ! Pivot tolerance;

! Case 1;
!Case01  dimr = 1..4;
!Case01  dimc = 1..4;
!Case01  A=
  5  7  3 -1
  1 -2  0  4
  1  2  4  5
  2 -4  0  8;

! Case 2;
!Case02  dimr = 1..4;
!Case02  dimc = 1..4;
!Case02  A=
  4  8  4 -2
  4 -8  0 16
 -2 -4 -2  1
  2 -4  0  8;

! Case 3;
!Case03  dimr = 1..2;
!Case03  dimc = 1..4;
!Case03 A=
  4  8  4 -2
  4 -8  0 16;

! Case 4;
!Case04;  dimr = 1..4!Case04;  dimc = 1..3!Case04; A=
1	3	4	
2	2	4
3	4	7
4	5	9;
ENDDATA

procedure comprank:
! Compute the rank of a dense matrix.
! Inputs:
     A = a matrix,
     tolpiv = a pivot tolerance, e.g., 10^(-13),
  Output:
    rank = rank of the matrix;

  nremr = @size( dimr); ! row dimension;
  nremc = @size( dimc); ! col dimension;
! Copy A into the WORK matrix;
  @FOR( rxc( i, j): WORK( i, j) = A( i, j));

! Load all the rows and cols onto the unpivoted lists;
  @FOR( dimr( i): listr( i) = i;);
  @FOR( dimc( j): listc( j) = j;);

  rank = 0; ! Count number of successful pivots;
 ! Keep pivoting while we have acceptable pivot;
  more = 1;
  @FOR( dimr( k) | more:
 ! Find the largest unpivoted element;
     pivMAX = 0;
     @FOR( dimr( kr) | kr #le# nremr:
        icand = listr( kr);
       @FOR( dimc( kc) | kc #le# nremc:
          jcand = listc( kc);
          @IFC( @ABS( work( icand, jcand)) #gt# pivmax:
            pivmax = @abs( work( icand, jcand)); !Found new big;
            ipos = kr; ! Remember list posn of row;
            jpos = kc; ! Remember list posn of col;
              );
            );
         );
! Did we find another acceptable pivot element? ;
      @IFC( pivmax #ge# tolpiv:
! Yes, get ready to do the pivot;
      rank = rank + 1;
      ipiv = listr( ipos); ! Get pivot row;
      jpiv = listc( jpos); ! Get pivot col;
! Remove pivot row and column from unpivoted lists;
      listr( ipos) = listr( nremr);
      nremr = nremr - 1;
      listc( jpos) = listc( nremc);
      nremc = nremc - 1;
! Now actually do the pivot, updating the unpivoted rows and cols;
      divisor = WORK( ipiv, jpiv);
      @FOR( dimr( kr) | kr #le# nremr:
        ii = listr( kr);
        mult = WORK( ii, jpiv)/ divisor;
        @FOR( dimc( kc) | kc #le# nremc:
          jj = listc( kc);
          WORK( ii, jj) = WORK( ii, jj) - mult*WORK( ipiv, jj);
            );
          );
    @else
       more = 0; ! Could not find a pivot. Get out;
          )
      );

  endprocedure

CALC:
  @SET( 'TERSEO',2);    ! Output level (0:verb, 1:terse, 2:only errors, 3:none);
  @SET( 'OROUTE',1);    ! 1: Route output immediately to the window line by line;
! Compute rank of matrix A;
  comprank;
ENDCALC
DATA:
  @TEXT() = 'Given matrix:';
  @TEXT() = @table( A);
  @TEXT() = ' The rank= ', rank;
ENDDATA