Row reduction macro

Question

This is more a suggestion/request than a question:

Someone should write a LaTeX macro that automatically row reduces a matrix until it's in (reduced) echelon form and typesets all the steps. (As far as I can tell, none such exists.)

I'm thinking of something like the gauss package, except that the row reductions themselves are carried out automatically, like in the Linear Algebra Toolkit. This would be similar to the \polylongdiv command in the polynom package, where all one needs to do is enter the polynomials to be divided and the macro carries out the algorithm and typesets the steps.

Of course you might be wondering why I don't just do it myself. Well I guess my answer to that is: some (linear) combination of laziness, busyness, not being the right person for the job, etc.

Thanks, regards, respect, and even a little love. :*

Answer 1

This answer has some macros picked up from https://tex.stackexchange.com/a/143035/4686. I am not too happy with some internal data structure, but I decided to live it standing.

The https://tex.stackexchange.com/a/143035/4686 computes determinants, inverses, etc..., either exactly or with float operations.

Here I focus on exact computations. The matrix entries may be integers, fractions, decimal numbers, or in scientific notation, but they are handled exactly. Hence, there is no question of numerical instability here. Regarding input, the format is with semi-colon separated rows of comma separated coefficients.

The last edit improves some internal aspects, has a better example for A=PLUQ, and redoes the initial example of Row Reduction to use for display truncated, not rounded, decimal expansions as they are followed with dots.

PLUQ ANSWER

The code typesets with TeX and also outputs to a file in Maple matrix notation the final result, for example.

A:=Matrix([[3, 1, -7, 5, 0, 9, -9, 7, -5], [-9, -4, 22, -14, 9, 2, 7, -6, -8], [-6, -3, 15, -9, 9, 11, -2, 1, -13], [-5, 8, 2, -18, 7, -1, 8, -7, 0], [4, 6, -14, 2, 1, -5, 6, 5, -3], [-11, 5, 17, -27, 16, 10, 6, -6, -13]]);
P:=Matrix([[1, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0], [0, 0, 0, 0, 1, 0], [0, 0, 1, 0, 0, 0], [0, 0, 0, 1, 0, 0], [0, 0, 0, 0, 0, 1]]);
L:=Matrix([[1, 0, 0, 0, 0, 0], [-3, 1, 0, 0, 0, 0], [-5/3, -29/3, 1, 0, 0, 0], [4/3, -14/3, 43/94, 1, 0, 0], [-2, 1, 0, 0, 1, 0], [-11/3, -26/3, 1, 0, 0, 1]]);
U:=Matrix([[3, 1, 0, 9, -7, 5, -9, 7, -5], [0, -1, 9, 29, 1, 1, -20, 15, -23], [0, 0, 94, 883/3, 0, 0, -601/3, 449/3, -692/3], [0, 0, 0, -1533/94, 0, 0, 1533/94, -263/94, 87/47], [0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0]]);
Q:=Matrix([[1, 0, 0, 0, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 1, 0, 0, 0, 0], [0, 0, 0, 0, 0, 1, 0, 0, 0], [0, 0, 1, 0, 0, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 1, 0, 0], [0, 0, 0, 0, 0, 0, 0, 1, 0], [0, 0, 0, 0, 0, 0, 0, 0, 1]]);

Now we can copy paste into Maple and check that indeed A=PLUQ:

> with(LinearAlgebra):
> MatrixAdd(A,-P.L.U.Q);
                  [0    0    0    0    0    0    0    0    0]
                  [                                         ]
                  [0    0    0    0    0    0    0    0    0]
                  [                                         ]
                  [0    0    0    0    0    0    0    0    0]
                  [                                         ]
                  [0    0    0    0    0    0    0    0    0]
                  [                                         ]
                  [0    0    0    0    0    0    0    0    0]
                  [                                         ]
                  [0    0    0    0    0    0    0    0    0]

Notice that in a PLUQ decomposition, a P and a Q will appear with my code only if necessary.

\documentclass[a4paper]{article}
\usepackage[hscale=0.85, vscale=0.85]{geometry}
\usepackage{xintfrac}
\usepackage{xinttools}
\usepackage{array} 
% \usepackage {siunitx}
% \usepackage {numprint}

\catcode`_ 11
\makeatletter

\newwrite\MATout
\immediate\openout\MATout=\jobname.pluqout\relax

% (the typeout format is for input in Maple for example)
\def\MATtypeout {\MATtypeoutwith {\MATtypeoutone}}%
\def\MATtypeoutone #1{\xintPRaw{\xintRawWithZeros{#1}}}% (lacking an \xintPRawWithZeros)
\def\MATtypeoutwith #1#2#3{%
     \edef\I{\xintSeq {1}{#3[I]}}% indices for rows
     \edef\J{\xintSeq {1}{#3[J]}}% indices for columns
     \immediate\write\MATout{#2:=Matrix([[%
       \xintListWithSep {], [}{\xintApply { \MAT_typeout_row {#1}#3}{\I}}%
    ]]);}%
}%
\def\MAT_typeout_row #1#2#3{%
    \xintListWithSep {, }{\xintApply { \MAT_typeout_one {#1}#2{#3}}{\J}}%
}%
\def\MAT_typeout_one #1#2#3#4{#1{#2[#3,#4]}}%

% we don't need all of them
\newcount\MAT_cnta
\newcount\MAT_cntb
\newcount\MAT_cntc
\newcount\MAT_cntd
\newcount\MAT_cnte

% Usage: \MATset\myMatrix{semi-colon separated rows of comma separated values}
% example.
% \MATset\MatrixA { 1/3 , 1/4, 1/5 ; 
%                   1/6 , 1/7 , 1/8 ; 
%                   1/9 , 1/10 , 1/11 ; }
% The final semi-colon is optional.

% We indeed focus here on manipulating matrices with rational entries, the
% code at https://tex.stackexchange.com/a/143035/4686 has the set-up for
% floating point numbers too (in an arbitrary, user decided precision).

\def\MATset {\def\MAT_xintin {\xintRaw}\MATset_ }%

\def\MATset_ #1#2{%
    \def\MATset_name{#1}%
    \edef\MAT_tmpa {#2}%
    \MAT_cnta \xint_c_ % sets \MAT_cnta to zero
    \expandafter\MATset_a 
    \romannumeral0\expandafter\xintzapspaces\expandafter{\MAT_tmpa};!;%
}%
\def\MATset_a {\futurelet\XINT_token\MATset_b }%
\def\MATset_b #1;{\def\MAT_tmpa{#1}%
                  \ifx\XINT_token;\expandafter\MATset_w
                  \else
                  \ifx\XINT_token!%
                          \expandafter\expandafter\expandafter\MATset_x
                     \else
                          \expandafter\expandafter\expandafter\MATset_c
                  \fi\fi }%
\def\MATset_w !;{\MATset_x }%
\def\MATset_x {\expandafter\def
  \csname MAT@\expandafter\string\MATset_name {I}\expandafter\endcsname
  \expandafter {\the\MAT_cnta }%
               \expandafter\def
  \csname MAT@\expandafter\string\MATset_name {J}\expandafter\endcsname 
  \expandafter {\the\MAT_cntb }%
  \expandafter\edef \MATset_name [##1]%
     {\noexpand\csname MAT@\expandafter\string\MATset_name 
               \noexpand\MAT_in ##1,\noexpand\xint_bye,\endcsname }%
}%
% a bit convoluted, no comments.
\def\MAT_in #1,#2,{\xint_bye #2\xint_gobble_iv\xint_bye
                   {\the\numexpr #1}{\the\numexpr #2}\xint_gobble_iii 
                   {\xintZapSpaces{#1}}}%
\def\MATset_c {\advance\MAT_cnta \xint_c_i % row count ++
               \MAT_cntb \xint_c_ % column count initially zero
               \expandafter\MATset_d\romannumeral0\expandafter
               \xintzapspaces\expandafter {\MAT_tmpa},!,}%
\def\MATset_d {\futurelet\XINT_token\MATset_e }%
\def\MATset_e #1,{\ifx\XINT_token!\expandafter\MATset_a
  \else
      \advance\MAT_cntb \xint_c_i
      \expandafter\def
  \csname MAT@\expandafter\string\MATset_name 
            {\the\MAT_cnta}{\the\MAT_cntb}\expandafter\endcsname
  \expandafter{\romannumeral-`0\MAT_xintin{\xintZapSpacesB{#1}}}%
  \expandafter\MATset_d\fi
}%

% removed \toks2 et \toks4 usage from https://tex.stackexchange.com/a/143035/4686
\def\MATlet #1#2{%
    \edef\MAT@seqI{\xintSeq {1}{#2[I]}}%
    \edef\MAT@seqJ{\xintSeq {1}{#2[J]}}%
    \xintFor* ##1 in {\MAT@seqI}
    \do{\xintFor* ##2 in {\MAT@seqJ}
        \do{\expandafter\let
     \csname MAT@\string#1{##1}{##2}\expandafter\endcsname
     \csname MAT@\string#2{##1}{##2}\endcsname
     }}%
  \expandafter\edef\csname MAT@\string#1{I}\endcsname {#2[I]}%
  \expandafter\edef\csname MAT@\string#1{J}\endcsname {#2[J]}%
  \edef #1[##1]%
     {\noexpand\csname 
      MAT@\string#1\noexpand\MAT_in ##1,\noexpand\xint_bye,\endcsname }%
}%

% We need identity matrices.
% again copied as is from https://tex.stackexchange.com/a/143035/4686
% IDENTITY MATRIX
% usage \MATid\foo{37} defines a 37 times 37 identity matrix.
\def\MATid {\def\MAT_tmpf{/1}\MAT_id }%
%\def\MATfloatid {\def\MAT_tmpf{}\MAT_id }%
% This identity matrix insists on coefficients written internally
% 0[0] or 1[0], this is a remnant of
% https://tex.stackexchange.com/a/143035/4686 whose aim is is minuscule
% optimization when these numbers are involved in computations done by
% the xintfrac macros.
\def\MAT_id #1#2{%
    \MAT_cntc #2\relax
    \MAT_cnta \xint_c_i % 1
    \xintloop
      {\expandafter\def\expandafter\MAT_tmpa \expandafter{\the\MAT_cnta}%
       \MAT_cntb \xint_c_i % 1
       \xintloop
         \expandafter\edef
         \csname MAT@\string#1{\MAT_tmpa}{\the\MAT_cntb}\endcsname 
           {\ifnum\MAT_cntb=\MAT_cnta 1\else 0\fi \MAT_tmpf[0]}%
       \ifnum\MAT_cntb<\MAT_cntc
         \advance\MAT_cntb \xint_c_i
       \repeat
     \ifnum\MAT_cnta<\MAT_cntc
        \advance\MAT_cnta \xint_c_i
    }\repeat
    \expandafter\def\csname MAT@\string#1{I}\expandafter\endcsname
            \expandafter {\the\MAT_cntc}%
    \expandafter\def\csname MAT@\string#1{J}\expandafter\endcsname 
            \expandafter {\the\MAT_cntc}%
    \edef #1[##1]%
       {\noexpand\csname 
         MAT@\string#1\noexpand\MAT_in ##1,\noexpand\xint_bye,\endcsname }%
}%

% EXCHANGING ROWS OR COLUMNS OF A GIVEN MATRIX
\def\MATexchangecol #1#2#3{%
    \MAT_cnta=#3[I]\relax
    \MAT_cntb=\xint_c_i % 1
    \xintloop
       \expandafter\let\expandafter\MAT@tmp
       \csname MAT@\string#3{\the\MAT_cntb}{#1}\endcsname
       \expandafter\let
       \csname MAT@\string#3{\the\MAT_cntb}{#1}\expandafter\endcsname
       \csname MAT@\string#3{\the\MAT_cntb}{#2}\endcsname
       \expandafter\let
       \csname MAT@\string#3{\the\MAT_cntb}{#2}\endcsname
       \MAT@tmp
    \ifnum\MAT_cntb<\MAT_cnta
       \advance\MAT_cntb\xint_c_i
    \repeat
}%
% perhaps only columns "to the right" actually need exchange in usage of this
\def\MATexchangerow #1#2#3{%
    \MAT_cnta=#3[J]\relax
    \MAT_cntb=\xint_c_i % 1
    \xintloop
       \expandafter\let\expandafter\MAT@tmp
       \csname MAT@\string#3{#1}{\the\MAT_cntb}\endcsname
       \expandafter\let
       \csname MAT@\string#3{#1}{\the\MAT_cntb}\expandafter\endcsname
       \csname MAT@\string#3{#2}{\the\MAT_cntb}\endcsname
       \expandafter\let
       \csname MAT@\string#3{#2}{\the\MAT_cntb}\endcsname
       \MAT@tmp
    \ifnum\MAT_cntb<\MAT_cnta
       \advance\MAT_cntb\xint_c_i
    \repeat
}%
\def\MATexchangerowspecial #1#2#3{%#1>#2, only columns <#2 need update
    \MAT_cnta=#2\relax
    \MAT_cntb=\xint_c_ % 0
    \xintloop
    \advance\MAT_cntb\xint_c_i
    \ifnum\MAT_cntb<\MAT_cnta
       \expandafter\let\expandafter\MAT@tmp
       \csname MAT@\string#3{#1}{\the\MAT_cntb}\endcsname
       \expandafter\let
       \csname MAT@\string#3{#1}{\the\MAT_cntb}\expandafter\endcsname
       \csname MAT@\string#3{#2}{\the\MAT_cntb}\endcsname
       \expandafter\let
       \csname MAT@\string#3{#2}{\the\MAT_cntb}\endcsname
       \MAT@tmp
    \repeat
}%


% Usage:
% \MATpluq\A (\A previously defined by \MATset)
% Effect: sets \P, \L, \U, \Q, to matrices in the sense of \MATset,
% so that "A=PLUQ" and it writes all matrices out
% to some file. See initial answer about row reduction for typesetting
% in document.
% The code is a simple adaptation of this initial answer. Now I use \MATpluq
% prefix.
\def\MATpluq #1{%
%  \begingroup
    \MATlet\@U#1%
    \edef\MATpluq@rows{\@U[I]}% nb of rows
    \edef\MATpluq@cols{\@U[J]}% nb of columns.
    \MATid\@P\MATpluq@rows
    \MATid\@L\MATpluq@rows
    \MATid\@Q\MATpluq@cols
    \def\MATpluq@pivrow {0}%
    \def\MATpluq@pivcol {0}%
    %\edef\MATpluq@name {\string#1}%
    \let\MATpluq@ifcontinue\iftrue
%  Starting the reduction.
    \MATtypeout{^^JA}#1%
\[A = \MATdisplay\@U\]
    \xintloop
% Nota Bene: in the PLUQ reduction, the pivots are anyhow organized
% along the main diagonal so pivrow and pivcol will be kept in sync over
% the execution of the algorithm but we use two variables nevertheless.
       \edef\MATpluq@pivrow{\the\numexpr\MATpluq@pivrow+\xint_c_i}%
       \edef\MATpluq@pivcol{\the\numexpr\MATpluq@pivcol+\xint_c_i}%
       \MATpluq@dopiv
    \MATpluq@ifcontinue
    \repeat
%  Done. The rank of the matrix is \the\numexpr\MATpluq@pivrow-\xint_c_i.\par
%  \endgroup
    \MATtypeout{P}\@P
    \MATtypeout{L}\@L
    \MATtypeout{U}\@U
    \MATtypeout{Q}\@Q
\[ P = \MATdisplay\@P\]
\[ L = \MATdisplay\@L\qquad U = \MATdisplay\@U\]
\[ Q = \MATdisplay\@Q\]
}

\def\MATpluq@done {\let\MATpluq@ifcontinue\iffalse}

% Remark on algorithm: I hesitated about doing column permutations first,
% rather than row permutations with the idea to recognize faster an entirely
% vanishing row, so that we can put it at the end and ignore it entirely, in
% effect reducing the number of rows by one, and possibly making algorithm
% faster. But for simplicity I just keep algorithm close to the one as in my
% initial answer. We only have to keep track in \P, \L, \Q of the needed
% operations.

\def\MATpluq@dopiv{%
    \let\MATpluq@row\MATpluq@pivrow
    \let\MATpluq@col\MATpluq@pivcol
    \ifnum\MATpluq@row>\MATpluq@rows\relax
        \MATpluq@done
    \else
        \ifnum\MATpluq@col>\MATpluq@cols\relax
            \MATpluq@done
        \else
            \expandafter\expandafter\expandafter\MATpluq@dopiv@i
        \fi
    \fi
}

\def\MATpluq@dopiv@i{%
    \edef\MATpluq@piv@value{\@U[\MATpluq@row,\MATpluq@col]}%
    \xintifZero{\MATpluq@piv@value}
        \MATpluq@dopiv@steprow
        \MATpluq@dopiv@ii
}

\def\MATpluq@dopiv@steprow{%
     \ifnum\MATpluq@row=\MATpluq@rows\relax
         \par No pivot found in column \MATpluq@col.\par
         \let\MATpluq@row\MATpluq@pivrow
         \expandafter\MATpluq@dopiv@stepcol
     \else
         \edef\MATpluq@row{\the\numexpr\MATpluq@row+\xint_c_i}%
         \expandafter\MATpluq@dopiv@i
     \fi
}

\def\MATpluq@dopiv@stepcol{%
     \ifnum\MATpluq@col=\MATpluq@cols\relax
         \MATpluq@done
     \else
         \edef\MATpluq@col{\the\numexpr\MATpluq@col+\xint_c_i}%
         \expandafter\MATpluq@dopiv@i
     \fi
}

% found a pivot
\def\MATpluq@dopiv@ii{%
Pivot \MATpluqprintonevalue{\MATpluq@piv@value} at \MATpluq@row, \MATpluq@col.\par
    \ifnum\MATpluq@col>\MATpluq@pivcol\relax
Exchange of columns \MATpluq@pivcol\space and \MATpluq@col.\par
      \MATexchangerow{\MATpluq@col}{\MATpluq@pivcol}\@Q
      \MATexchangecol{\MATpluq@col}{\MATpluq@pivcol}\@U
\[U = \MATdisplay\@U\qquad Q = \MATdisplay\@Q\]
    \fi
    \ifnum\MATpluq@pivrow=\MATpluq@rows\relax
      \edef\MATpluq@pivrow{\the\numexpr\MATpluq@pivrow+\xint_c_i}%
      \MATpluq@done
    \else
      \expandafter\MATpluq@dopiv@iii
    \fi
}

\def\MATpluq@dopiv@iii{%
    \ifnum\MATpluq@row>\MATpluq@pivrow\relax
Exchange of rows \MATpluq@pivrow\space and \MATpluq@row.\par
      \MATexchangecol{\MATpluq@row}{\MATpluq@pivrow}\@P
      \MATexchangerow{\MATpluq@row}{\MATpluq@pivrow}\@U
      \MATexchangerowspecial{\MATpluq@row}{\MATpluq@pivrow}\@L
\[L = \MATdisplay\@L\qquad U = \MATdisplay\@U\]
\[P = \MATdisplay\@P\]
    \fi
    \MAT_cntc\MATpluq@pivrow\relax% we are guaranteed < nb of rows
    \xintloop
      \advance\MAT_cntc\xint_c_i 
      \edef\MATpluq@entry{\@U[\MAT_cntc,\MATpluq@pivcol]}%
      \xintifZero\MATpluq@entry
     {% nothing to do, the L coeff is already set to zero
     }%
     {\edef\MATpluq@ratio
         {\xintIrr{\xintDiv{\MATpluq@entry}{\MATpluq@piv@value}}[0]}%
      \expandafter\let
      \csname MAT@\string\@L{\the\MAT_cntc}{\MATpluq@pivcol}\endcsname
      \MATpluq@ratio
Subtract from row \the\MAT_cntc\space the pivot row multiplied by
\MATpluqprintonevalue{\MATpluq@ratio}.\par
      \@namedef{MAT@\string\@U{\the\MAT_cntc}{\MATpluq@pivcol}}{0[0]}%
      \MAT_cntd\MATpluq@pivcol\relax
      \xintloop
      \advance\MAT_cntd\xint_c_i
      \unless\ifnum\MATpluq@cols<\MAT_cntd
        \expandafter\edef
        \csname MAT@\string\@U{\the\MAT_cntc}{\the\MAT_cntd}\endcsname
          {\xintIrr{%
           \xintSub{\@U[\MAT_cntc,\MAT_cntd]}
                   {\xintMul{\MATpluq@ratio}{\@U[\MATpluq@pivrow,\MAT_cntd]}}%
           }[0]}%
      \repeat
      }%
    \unless\ifnum\MATpluq@rows=\MAT_cntc
    \repeat
\[L = \MATdisplay\@L\qquad U = \MATdisplay\@U\]
}

\def\MATpluqprintonevalue{\xintPRaw}
%\def\MATpluqdisplay#1{\[\MATdisplay#1\]}%

%% MATH MODE MATRIX DISPLAY

\makeatother

\newcommand\MATdisplay [1][1.25]{\MATdisplaywith [#1]{\MATdisplayone}}
\def\MATdisplayone {\xintSignedFrac}
\newcolumntype\MATdisplaycoltype {c}
\newcolumntype\MATdisplaypreamble [1]{@{}*{#1[J]}\MATdisplaycoltype@{}}
\newcommand\MATdisplaywith [3][1.25]
   {\left(\def\arraystretch{#1}%
    \begin{array}{\MATdisplaypreamble {#3}}
         \xintListWithSep {\\}
       {\xintApply { \MAT_display_row {#2}#3}{\xintSeq {1}{#3[I]}}}
    \end{array}\right)%
}%
\def\MAT_display_row #1#2#3{%
    \xintListWithSep {&}
     {\xintApply{ \MAT_display_one {#1}#2{#3}}{\xintSeq {1}{#2[J]}}}%
}%
\def\MAT_display_one #1#2#3#4{#1{#2[#3,#4]}}%

\catcode`_ 8   

\begin{document}\pagestyle{empty}
\MATset\MatrixA { 1/3 , 1/4, 1/5 ; 
                  1/6 , 1/7 , 1/8 ; 
                  1/9 , 1/10 , 1/11 ; }

\MATpluq\MatrixA

See pluqout file.\clearpage

\MATset\A {
3, -7, 5, 0, 1, 0, 1;
-9, -8, -2, 9, -1, 9, -4;
4, 6, 0, -1, -2, -1, -3;
-5, 2, -6, 7, 8, 7, 8;
-1, -2, -1, -3, 4, 6, 0;
7, 8, 7, 8, -5, 2, -6;
}

\MATpluq\A

See pluqout file.\clearpage

\MATset\A {
2, 0, 3, 0;
1, 0, 0, 0;
0, 0, 4, 0;
0, 2, 0, 1;
}

\MATpluq\A

See pluqout file.\clearpage

\MATset\MatrixB {
 3, 1, -7, 5, 0, 9, -9, 7, -5;
 -9, -4, -8, -2, 9, 2, 8, -6, -8;
 4, -3, 6, 0, -1, 5, -4, -3, 4;
 -5, 8, 2, -6, 7, -1, 1, -7, 0;
 3, 6, -2, -1, 8, -2, -6, 7, -7;
 4, 6, 3, -9, 1, -5, 0, 5, -3;
}

\MATpluq\MatrixB

See pluqout file.\clearpage

\MATset\MatrixC {
  3,  1,  -7,   5,  0,  9, -9,  7,  -5;
 -9, -4,  22, -14,  9,  2,  7, -6,  -8;
 -6, -3,  15,  -9,  9, 11, -2,  1, -13;
 -5,  8,   2, -18,  7, -1,  8, -7,   0;
  4,  6, -14,   2,  1, -5,  6,  5,  -3;
-11,  5,  17, -27, 16, 10,  6, -6, -13;
}

\MATpluq\MatrixC

See pluqout file for the matrices in Maple format.\clearpage

\immediate\closeout\MATout
\end{document}

ROW REDUCTION (INITIAL) ANSWER

I have improved a bit some internal aspects of the code in an edit.

\documentclass{article}
\usepackage{xintfrac}
\usepackage{xinttools}
\usepackage{array} 

\catcode`_ 11
\makeatletter

\newcount\MAT_cnta
\newcount\MAT_cntb
\newcount\MAT_cntc
\newcount\MAT_cntd
\newcount\MAT_cnte

% Usage: \MATset\myMatrix{semi-colon separated rows of comma separated values}
% example.
% \MATset\MatrixA { 1/3 , 1/4, 1/5 ; 
%                   1/6 , 1/7 , 1/8 ; 
%                   1/9 , 1/10 , 1/11 ; }

\def\MATset {\def\MAT_xintin {\xintRaw}\MATset_ }%

\def\MATset_ #1#2{%
    \def\MATset_name{#1}%
    \edef\MAT_tmpa {#2}%
    \MAT_cnta \xint_c_ % sets \MAT_cnta to zero
    \expandafter\MATset_a 
    \romannumeral0\expandafter\xintzapspaces\expandafter{\MAT_tmpa};!;%
}%
\def\MATset_a {\futurelet\XINT_token\MATset_b }%
\def\MATset_b #1;{\def\MAT_tmpa{#1}%
                  \ifx\XINT_token;\expandafter\MATset_w
                  \else
                  \ifx\XINT_token!%
                          \expandafter\expandafter\expandafter\MATset_x
                     \else
                          \expandafter\expandafter\expandafter\MATset_c
                  \fi\fi }%
\def\MATset_w !;{\MATset_x }%
\def\MATset_x {\expandafter\def
  \csname MAT@\expandafter\string\MATset_name {I}\expandafter\endcsname
  \expandafter {\the\MAT_cnta }%
               \expandafter\def
  \csname MAT@\expandafter\string\MATset_name {J}\expandafter\endcsname 
  \expandafter {\the\MAT_cntb }%
  \expandafter\edef \MATset_name [##1]%
     {\noexpand\csname MAT@\expandafter\string\MATset_name 
               \noexpand\MAT_in ##1,\noexpand\xint_bye,\endcsname }%
}%
%
\def\MAT_in #1,#2,{\xint_bye #2\xint_gobble_iv\xint_bye
                   {\the\numexpr #1}{\the\numexpr #2}\xint_gobble_iii 
                   {\xintZapSpaces{#1}}}%
%
\def\MATset_c {\advance\MAT_cnta \xint_c_i % row count ++
               \MAT_cntb \xint_c_ % column count initially zero
               \expandafter\MATset_d\romannumeral0\expandafter
               \xintzapspaces\expandafter {\MAT_tmpa},!,}%
\def\MATset_d {\futurelet\XINT_token\MATset_e }%
\def\MATset_e #1,{\ifx\XINT_token!\expandafter\MATset_a
  \else
      \advance\MAT_cntb \xint_c_i
      \expandafter\def
  \csname MAT@\expandafter\string\MATset_name 
            {\the\MAT_cnta}{\the\MAT_cntb}\expandafter\endcsname
  \expandafter{\romannumeral-`0\MAT_xintin{\xintZapSpacesB{#1}}}%
  \expandafter\MATset_d\fi
}%

\def\MATlet #1#2{%
    \edef\MAT@seqI{\xintSeq {1}{#2[I]}}%
    \edef\MAT@seqJ{\xintSeq {1}{#2[J]}}%
    \xintFor* ##1 in {\MAT@seqI}
    \do{\xintFor* ##2 in {\MAT@seqJ}
        \do{\expandafter\let
     \csname MAT@\string#1{##1}{##2}\expandafter\endcsname
     \csname MAT@\string#2{##1}{##2}\endcsname
     }}%
  \expandafter\edef\csname MAT@\string#1{I}\endcsname {#2[I]}%
  \expandafter\edef\csname MAT@\string#1{J}\endcsname {#2[J]}%
  \edef #1[##1]%
     {\noexpand\csname 
      MAT@\string#1\noexpand\MAT_in ##1,\noexpand\xint_bye,\endcsname }%
}%

\def\MATrowreduce #1{%
  \begingroup
    \edef\MATrr@rows{#1[I]}%
    \edef\MATrr@cols{#1[J]}%
    \def\MATrr@pivrow {0}%
    \def\MATrr@pivcol {0}%
    \MATlet\@U #1%
    \let\MATrr@ifcontinue\iftrue
Starting the reduction.
\MATrrdisplaymatrix\@U
    \xintloop
       \edef\MATrr@pivrow{\the\numexpr\MATrr@pivrow+\xint_c_i}%
       \edef\MATrr@pivcol{\the\numexpr\MATrr@pivcol+\xint_c_i}%
       \MATrr@dopiv
    \MATrr@ifcontinue
    \repeat
Done. The rank of the matrix is \the\numexpr\MATrr@pivrow-\xint_c_i.\par
  \endgroup
}

\def\MATrr@done {\let\MATrr@ifcontinue\iffalse}

\def\MATrr@dopiv{%
    \let\MATrr@row\MATrr@pivrow
    \let\MATrr@col\MATrr@pivcol
    \ifnum\MATrr@row>\MATrr@rows\relax
        \MATrr@done
    \else
        \ifnum\MATrr@col>\MATrr@cols\relax
            \MATrr@done
        \else
            \expandafter\expandafter\expandafter\MATrr@dopiv@i
        \fi
    \fi
}

\def\MATrr@dopiv@i{%
    \edef\MATrr@piv@value{\@U[\MATrr@row,\MATrr@pivcol]}%
    \xintifZero{\MATrr@piv@value}
        \MATrr@dopiv@steprow
        \MATrr@dopiv@ii
}

\def\MATrr@dopiv@steprow{%
     \ifnum\MATrr@row=\MATrr@rows\relax
         \let\MATrr@row\MATrr@pivrow
\par No pivot found in column \MATrr@pivcol.\par
         \expandafter\MATrr@dopiv@stepcol
     \else
         \edef\MATrr@row{\the\numexpr\MATrr@row+\xint_c_i}%
         \expandafter\MATrr@dopiv@i
     \fi
}

\def\MATrr@dopiv@stepcol{%
     \ifnum\MATrr@pivcol=\MATrr@cols\relax
         \MATrr@done
     \else
         \edef\MATrr@pivcol{\the\numexpr\MATrr@pivcol+\xint_c_i}%
         \expandafter\MATrr@dopiv@i
     \fi
}

\def\MATrr@dopiv@ii{%
    \ifnum\MATrr@pivrow=\MATrr@rows\relax
      \edef\MATrr@pivrow{\the\numexpr\MATrr@pivrow+\xint_c_i}\MATrr@done
    \else
      \expandafter\MATrr@dopiv@iii
    \fi
}

\def\MATrr@dopiv@iii{%
Now using the pivot with value \MATrrprintonevalue{\MATrr@piv@value}
at row \MATrr@row\space and column \MATrr@pivcol.\par
    \ifnum\MATrr@row>\MATrr@pivrow\relax
Exchange of row \MATrr@row\space with row \MATrr@pivrow.\par
      \MAT_cntb=\MATrr@pivcol\relax
      \xintloop
         \expandafter\let\expandafter\MAT@tmp
         \csname MAT@\string\@U{\MATrr@row}{\the\MAT_cntb}\endcsname
         \expandafter\let
         \csname MAT@\string\@U{\MATrr@row}{\the\MAT_cntb}\expandafter\endcsname
         \csname MAT@\string\@U{\MATrr@pivrow}{\the\MAT_cntb}\endcsname
         \expandafter\let
         \csname MAT@\string\@U{\MATrr@pivrow}{\the\MAT_cntb}\endcsname
         \MAT@tmp
      \ifnum\MATrr@cols>\MAT_cntb
         \advance\MAT_cntb\xint_c_i
      \repeat
\MATrrdisplaymatrix\@U\par
    \fi
    \MAT_cntc\MATrr@pivrow
    \xintloop
      \advance\MAT_cntc\xint_c_i
      \edef\MATrr@entry{\@U[\MAT_cntc,\MATrr@pivcol]}%
      \xintifZero\MATrr@entry
     {}%
     {\edef\MATrr@ratio{\xintIrr{\xintDiv{\MATrr@entry}{\MATrr@piv@value}}[0]}%
Subtract from row \the\MAT_cntc\space the pivot row multiplied by
\MATrrprintonevalue{\MATrr@ratio}.\par
      \@namedef{MAT@\string\@U{\the\MAT_cntc}{\MATrr@pivcol}}{0[0]}%
      \MAT_cntd\MATrr@pivcol\relax
      \xintloop
      \advance\MAT_cntd\xint_c_i
      \unless\ifnum\MATrr@cols<\MAT_cntd
        \expandafter\edef
        \csname MAT@\string\@U{\the\MAT_cntc}{\the\MAT_cntd}\endcsname
          {\xintIrr{%
           \xintSub{\@U[\MAT_cntc,\MAT_cntd]}
                   {\xintMul{\MATrr@ratio}{\@U[\MATrr@pivrow,\MAT_cntd]}}%
           }[0]}%
      \repeat
      }%
    \unless\ifnum\MATrr@rows=\MAT_cntc
    \repeat
\MATrrdisplaymatrix\@U
}

\def\MATrrprintonevalue{\xintPRaw}
\def\MATrrdisplaymatrix #1{\[\MATdisplay#1\]}%

%% MATH MODE MATRIX DISPLAY

\makeatother

\newcommand\MATdisplay [1][1.25]{\MATdisplaywith [#1]{\MATdisplayone}}
\def\MATdisplayone {\xintSignedFrac}

\newcolumntype\MATdisplaycoltype {c}
\newcolumntype\MATdisplaypreamble [1]{@{}*{#1[J]}\MATdisplaycoltype@{}}

\newcommand\MATdisplaywith [3][1.25]
   {\left(\def\arraystretch{#1}%
    \begin{array}{\MATdisplaypreamble {#3}}
         \xintListWithSep {\\}
       {\xintApply { \MAT_display_row {#2}#3}{\xintSeq {1}{#3[I]}}}
    \end{array}\right)%
}%

\def\MAT_display_row #1#2#3{%
    \xintListWithSep {&}
     {\xintApply{ \MAT_display_one {#1}#2{#3}}{\xintSeq {1}{#2[J]}}}%
}%

\def\MAT_display_one #1#2#3#4{#1{#2[#3,#4]}}%

\catcode`_ 8   

\begin{document}
\MATset\MatrixC {
 3, 1, -7, 5, 0, 9, -9, 7, -5;
 -9, -4, -8, -2, 9, 2, 8, -6, -8;
 -6, -3, -15, 3, 9, 11, -1, 1, -13;
 -5, 8, 2, -6, 7, -1, 1, -7, 0;
 4, 6, 3, -9, 1, -5, 0, 5, -3;
 -11, 5, -13, -3, 16, 10, 0, -6, -13;
}

\MATrowreduce\MatrixC

\end{document}

Entries may be decimal numbers like 37.156.

\def\MATrrprintonevalue{\xintRound{2}}
\def\MATrrdisplaymatrix #1{\[\MATdisplaywith{\xintRound{2}}#1\]}%

Example (as I add dots, I use truncating rather than rounding):

\def\MATrrprintonevalue#1{\xintTrunc{3}{#1}\dots\ (=\xintPRaw{#1})}
\def\MATrrdisplaymatrix #1{\[\MATdisplay#1=\MATdisplaywith{\TruncWithDots{3}}#1\]}%
\def\TruncWithDots #1#2{\xintTrunc{#1}{#2}...}
\MATset\MatrixA { 1/3 , 1/4, 1/5 ; 
                  1/6 , 1/7 , 1/8 ; 
                  1/9 , 1/10 , 0.09 ; }

\MATrowreduce\MatrixA

Answer 2

Update-2

I heard someone said Givens rotations.

% Givens rotation
% I assume #1 < #2
% does not use theta because it is unstable
% #4 is cosine and #5 is sine
\def\pgflabgivensrotaterow #1 and row #2 by #3 and #4 in #5{
    \pgfkeys{/lab/#5/w/.get=\pgflabw}
    \pgfplotsforeachungrouped\g@j in{1,...,\pgflabw}{
        \pgfkeys{/lab/#5/#1/\g@j/.get=\pgflabtempentrya}
        \pgfkeys{/lab/#5/#2/\g@j/.get=\pgflabtempentryb}
        \pgfmathparse{#3*\pgflabtempentrya-#4*\pgflabtempentryb}
        \pgfkeys{/lab/#5/#1/\g@j/.let=\pgfmathresult}
        \pgfmathparse{#4*\pgflabtempentrya+#3*\pgflabtempentryb}
        \pgfkeys{/lab/#5/#2/\g@j/.let=\pgfmathresult}
    }
}

% I assume #1 < #2
% does not use theta because it is unstable
% #4 is cosine and #5 is sine
\def\pgflabgivensrotatecol #1 and col #2 by #3 and #4 in #5{
    \pgfkeys{/lab/#5/h/.get=\pgflabh}
    \pgfplotsforeachungrouped\g@i in{1,...,\pgflabh}{
        \pgfkeys{/lab/#5/\g@i/#1/.get=\pgflabtempentrya}
        \pgfkeys{/lab/#5/\g@i/#2/.get=\pgflabtempentryb}
        \pgfmathparse{#3*\pgflabtempentrya-#4*\pgflabtempentryb}
        \pgfkeys{/lab/#5/\g@i/#1/.let=\pgfmathresult}
        \pgfmathparse{#4*\pgflabtempentrya+#3*\pgflabtempentryb}
        \pgfkeys{/lab/#5/\g@i/#2/.let=\pgfmathresult}
    }
}

% A = QR decomposition
\def\pgflabQRdecompose #1 as #2 times #3{
    \pgfkeys{/lab/#1/w/.get=\pgflabW}
    \pgfkeys{/lab/#1/h/.get=\pgflabH}
    % decide the loop boundary
    \edef\pgflab@H-1{\the\numexpr\pgflabH-1}
    \ifnum\pgflab@H-1>\pgflabW
        \edef\pgflab@H-1{\pgflabW}
    \fi
    % set Q as identity
    % set #2 as identity
    \pgflabneweyeof {\pgflabH} by {\pgflabH} as {#2}
    % copy A to R
    % copy #1 to #3
    \pgflabcopymatrix {#1} to {#3}
    % forget A, do job at Q and R
    % forget #1, do job at #2 and #3
    \pgfplotsforeachungrouped\d@i in{1,...,\pgflab@H-1}{
        \edef\d@@i+1{\the\numexpr\d@i+1}
        \pgfplotsforeachungrouped\d@j in{\d@@i+1,...,\pgflabh}{
            \pgfkeys{/lab/#3/\d@i/\d@i/.get=\pgflabtempentrya}
            \pgfkeys{/lab/#3/\d@j/\d@i/.get=\pgflabtempentryb}
            \pgfmathsetmacro\pgflabtempradius{sqrt(\pgflabtempentrya*\pgflabtempentrya+\pgflabtempentryb*\pgflabtempentryb)}
            \pgfmathsetmacro\pgflabtempcos{ \pgflabtempentrya/\pgflabtempradius} % cosine
            \pgfmathsetmacro\pgflabtempsin{-\pgflabtempentryb/\pgflabtempradius} % sine
            \pgflabgivensrotaterow {\d@i} and row {\d@j} by {\pgflabtempcos} and {\pgflabtempsin} in {#3}
            \pgflabgivensrotatecol {\d@i} and col {\d@j} by {\pgflabtempcos} and {\pgflabtempsin} in {#2}
            eliminate one entry. check Q and R\par
            $Q=\pgflabtypeset{#2};$
            $R=\pgflabtypeset{#3};$
        }
    }
}
\pgflabread{A}{
    0   0   0   1
    1   0   0   0
    0   1   0   0
    0   0   1   0
}
\pgflabQRdecompose A as Q times R

For a 10 by 10 random matrix, the norm of A - QR is about 4e-4. The norm of QQᵀ - I is about 2e-4.

Update-1: New Answer

I implement three decompositions:

• A = LU
• A = PLU (i.e. partial pivoting)
• A = PLUQ (i.e. complete pivoting)

If A is m by n, then P, L are m by m; U is the same as A; and Q is n by n.

Advantages

• The complexity of accessing a matrix entry is O(1). (Assuming \csname is O(1)). So the complexity of decompositions is O(m²n).

• The input utilize \pgfplotstableread from pgfplotstable. So it accepts inline-table, file, loaded table, and even the table created by \pgfplotstablenew. You can also pass options to it. (such as filtering)

• The output utilize \pgfplotstabletypeset from the same package. Or you can convert the matrix back to a table and do whatever you want.

• The calculation is done by \pgfmathparse. I assume FPU is on. But one can reimplement that.

• There is a debug macro that output the raw data of matrices. You can copy and paste those data into whatever modern matrix calculator.

• According to Wikipeida, even partial pivoting is numerically stable in practice. I test a 10 by 10 random matrix and check A - PLUQ in sage; the norm is about 1.1e-6. (This is about the precision of FPU)

---

\documentclass{article}
\usepackage[a3paper,landscape,margin=1cm]{geometry}
\usepackage{pgfplotstable,mathtools}
\pgfplotsset{compat=newest}
\pgfkeys{/pgf/fpu,/pgf/number format/fixed}
\begin{document}


\makeatletter
% \pgfmatrix... is used
% we use \pgflab...

% call pgfplotstable to read the data
% put options in [] if desired
% the options go to \pgfplotstableread
\def\pgflabread{
    \pgfutil@ifnextchar[
        {\pgflabread@opt}
        {\pgflabread@opt[]}
}

% #1: optional option
% #2: a name of the matrix... usually A
\def\pgflabread@opt[#1]#2{
    \edef\pgflabname{#2}
    \pgfplotstableread[header=false,#1]
}

% we did not provide a macro to pgfplotstable to store the table
% we give it a temporary one called \pgflabtemptable
% and then copy it to our data structure
\long\def\pgfplotstableread@impl@collectfirstarg#1#2{
    \pgfplotstableread@impl@{#1}{#2}\pgflabtemptable
    \pgflabconverttable\pgflabtemptable to matrix{\pgflabname}
}

% this helps us to deal with pgfleys
\pgfkeys{/handlers/.let/.code=\pgfkeyslet{\pgfkeyscurrentpath}{#1}}

% copy pgfplotstable table to our data structure in pgfkeys
% #1: the macro that pgfplotstable used to store the table
% #2: a name of the matrix
\def\pgflabconverttable#1to matrix#2{
    % extract height and width
    \pgfplotstablegetrowsof#1\xdef\pgflabh{\pgfplotsretval}\pgfkeys{/lab/#2/h/.let=\pgflabh}
    %%%height = \pgflabh \par
    \pgfplotstablegetcolsof#1\xdef\pgflabw{\pgfplotsretval}\pgfkeys{/lab/#2/w/.let=\pgflabw}
    %%%width = \pgflabw \par
    % extract entries
    % \c@i and \c@j cannot be used outside
    \pgfplotsforeachungrouped\c@i in{1,...,\pgflabh}{
        \pgfplotsforeachungrouped\c@j in{1,...,\pgflabw}{
            % since fpu is on, this is easier way to do 9-1
            \pgfplotstablegetelem{\the\numexpr\c@i-1}{\the\numexpr\c@j-1}\of\pgflabtemptable
            \pgfkeys{/lab/#2/\c@i/\c@j/.let=\pgfplotsretval}
            %%%\pgfplotsretval,
        }
        %%%; \par
    }
}
\pgflabread{A}{
    3   1   -7  5   0
    -9  -4  -8  -2  9
    4   -3  6   0   -1
    -5  8   2   -6  7
}

% the opposite of the previous one
% #1: the name of the matrix
% #2: a macro for pgfplotstable to store the table
\def\pgflabconvertmatrix #1 to table #2{
    % makeup meta data
    \expandafter\def\csname\string#2@@table@name\endcsname{<inline_table>}
    % build a new list of columns
    \pgfkeys{/lab/#1/h/.get=\pgflabh}
    \pgfkeys{/lab/#1/w/.get=\pgflabw}
    \pgfplotslistnew#2{0,...,\the\numexpr\pgflabw-1}
    % fill in columns
    \pgfplotsforeachungrouped\c@j in{1,...,\pgflabw}{
        \pgfplotslistnewempty\pgflabtempcolumn
        \pgfplotsforeachungrouped\c@i in{1,...,\pgflabh}{
            \pgfkeys{/lab/#1/\c@i/\c@j/.get=\pgflabtempentry}
            \expandafter\pgfplotslistpushback\pgflabtempentry\to\pgflabtempcolumn
        }
        \edef\c@k{\the\numexpr\c@j-1}
        \expandafter\let\csname\string#2@\c@k\endcsname\pgflabtempcolumn
    }
}

% typeset the matrix by \pgfplotstabletypeset
\def\pgflabtypeset{
    \pgfutil@ifnextchar[
        {\pgflabtypeset@opt}
        {\pgflabtypeset@opt[]}
}

% #1: optional option
% #2: the name of the matrix
\def\pgflabtypeset@opt[#1]#2{
    \pgflabconvertmatrix #2 to table \pgflabtemptable
    \pgfplotstabletypeset[every head row/.style={output empty row}]\pgflabtemptable
}
Matrix A is
$A=\pgflabtypeset{A}$

% define row operation: switch
% does not check boundary
\def\pgflabswitchrow #1 and row #2 in #3{
    \pgfkeys{/lab/#3/w/.get=\pgflabw}
    \pgfplotsforeachungrouped\s@j in{1,...,\pgflabw}{
        \pgfkeys{/lab/#3/#1/\s@j/.get=\pgflabtempentrya}
        \pgfkeys{/lab/#3/#2/\s@j/.get=\pgflabtempentryb}
        \pgfkeys{/lab/#3/#1/\s@j/.let=\pgflabtempentryb}
        \pgfkeys{/lab/#3/#2/\s@j/.let=\pgflabtempentrya}
    }
}
\bigskip
\pgflabswitchrow 1 and row 3 in A
switch row 1 and row 3;
$A=\pgflabtypeset{A}$

% define column operation: switch
% does not check boundary
\def\pgflabswitchcol #1 and col #2 in #3{
    \pgfkeys{/lab/#3/h/.get=\pgflabh}
    \pgfplotsforeachungrouped\s@i in{1,...,\pgflabh}{
        \pgfkeys{/lab/#3/\s@i/#1/.get=\pgflabtempentrya}
        \pgfkeys{/lab/#3/\s@i/#2/.get=\pgflabtempentryb}
        \pgfkeys{/lab/#3/\s@i/#1/.let=\pgflabtempentryb}
        \pgfkeys{/lab/#3/\s@i/#2/.let=\pgflabtempentrya}
    }
}
\bigskip
\pgflabswitchcol 2 and col 3 in A
switch col 2 and col 3;
$A=\pgflabtypeset{A}$

% define row operation: multiplication
% does not check boundary
\def\pgflabmultiplyrow #1 by #2 in #3{
    \pgfkeys{/lab/#3/w/.get=\pgflabw}
    \pgfplotsforeachungrouped\m@j in{1,...,\pgflabw}{
        \pgfkeys{/lab/#3/#1/\m@j/.get=\pgflabtempentry}
        \pgfmathparse{\pgflabtempentry*#2}
        \pgfkeys{/lab/#3/#1/\m@j/.let=\pgfmathresult}
    }
}
\bigskip
\pgflabmultiplyrow 3 by -1 in A
multiply row 3 by -1;
$A=\pgflabtypeset{A}$

% define row operation: addition
% does not check boundary
\def\pgflabaddrow #1 by row #2 times #3 in #4{
    \pgfkeys{/lab/#4/w/.get=\pgflabw}
    \pgfplotsforeachungrouped\a@j in{1,...,\pgflabw}{
        \pgfkeys{/lab/#4/#1/\a@j/.get=\pgflabtempentrya}
        \pgfkeys{/lab/#4/#2/\a@j/.get=\pgflabtempentryb}
        \pgfmathparse{\pgflabtempentrya+\pgflabtempentryb*#3}
        \pgfkeys{/lab/#4/#1/\a@j/.let=\pgfmathresult}
    }
}
\bigskip
\pgflabaddrow 2 by row 3 times 2 in A
add row 2 by row 3 times 2;
$A=\pgflabtypeset{A}$

% define column operation: addition
% does not check boundary
\def\pgflabaddcol #1 by col #2 times #3 in #4{
    \pgfkeys{/lab/#4/h/.get=\pgflabh}
    \pgfplotsforeachungrouped\a@i in{1,...,\pgflabh}{
        \pgfkeys{/lab/#4/\a@i/#1/.get=\pgflabtempentrya}
        \pgfkeys{/lab/#4/\a@i/#2/.get=\pgflabtempentryb}
        \pgfmathparse{\pgflabtempentrya+\pgflabtempentryb*#3}
        \pgfkeys{/lab/#4/\a@i/#1/.let=\pgfmathresult}
    }
}
\bigskip
\pgflabaddcol 5 by col 4 times -1 in A
add col 5 by row 4 times -1;
$A=\pgflabtypeset{A}$

% new identity matrix
\def\pgflabneweyeof #1 by #2 as #3{
    \def\pgflabh{#1}\pgfkeys{/lab/#3/h/.let=\pgflabh}
    \def\pgflabw{#2}\pgfkeys{/lab/#3/w/.let=\pgflabw}
    \pgfplotsforeachungrouped\n@i in{1,...,\pgflabh}{
        \pgfplotsforeachungrouped\n@j in{1,...,\pgflabw}{
            \ifnum\n@i=\n@j
                \pgfkeys{/lab/#3/\n@i/\n@i/.initial=1}
            \else
                \pgfkeys{/lab/#3/\n@i/\n@j/.initial=0}
            \fi
        }
    }
}
\bigskip
\pgflabneweyeof 4 by 4 as I
identity matrix;
$A=\pgflabtypeset{I}$

\bigskip
\pgflabneweyeof 3 by 5 as B
rectangular identity matrix;
$B=\pgflabtypeset{B}$

% copy matrix
\def\pgflabcopymatrix #1 to #2{
    \pgfkeys{/lab/#1/h/.get=\pgflabh}\pgfkeys{/lab/#2/h/.let=\pgflabh}
    \pgfkeys{/lab/#1/w/.get=\pgflabw}\pgfkeys{/lab/#2/w/.let=\pgflabw}
    \pgfplotsforeachungrouped\n@i in{1,...,\pgflabh}{
        \pgfplotsforeachungrouped\n@j in{1,...,\pgflabw}{
            \pgfkeys{/lab/#1/\n@i/\n@j/.get=\pgflabtempentry}
            \pgfkeys{/lab/#2/\n@i/\n@j/.let=\pgflabtempentry}
        }
    }
}
\bigskip
\pgflabcopymatrix A to B
copy matrix A to B;
$B=\pgflabtypeset{B}$

% LU decomposition
% if encounter 0, probably will result in inf or nan
\def\pgflabLUdecompose #1 as #2 times #3{
    \pgfkeys{/lab/#1/h/.get=\pgflabh@u}
    \pgfkeys{/lab/#1/w/.get=\pgflabw@u}
    % decide the loop boundary
    \edef\pgflabh@v{\the\numexpr\pgflabh@u-1}
    \ifnum\pgflabh@v>\pgflabw@u
        \edef\pgflabh@v{\pgflabw@u}
    \fi
    % set L as identity
    % set #2 as identity
    \pgflabneweyeof {\pgflabh@u} by {\pgflabh@u} as #2
    % copy A to U
    % copy #1 to #3
    \pgflabcopymatrix #1 to #3
    % forget A, do job at L and U
    % forget #1, do job at #2 and #3
    \pgfplotsforeachungrouped\d@i in{1,...,\pgflabh@v}{
        \edef\d@@i+1{\the\numexpr\d@i+1}
        \pgfplotsforeachungrouped\d@j in{\d@@i+1,...,\pgflabh@u}{
            % use (\d@i,\d@i) to eliminate (\d@j,\d@i)
            \pgfkeys{/lab/#3/\d@i/\d@i/.get=\pgflabtempentrya}
            \pgfkeys{/lab/#3/\d@j/\d@i/.get=\pgflabtempentryb}
            \pgfmathsetmacro\pgflabtempratio{\pgflabtempentryb/\pgflabtempentrya}
            \pgflabaddcol {\d@i} by col {\d@j} times {\pgflabtempratio} in {#2}
            \pgflabaddrow {\d@j} by row {\d@i} times {-\pgflabtempratio} in {#3}

            \medskip
            eliminate one entry. check L and U \par
            $L=\pgflabtypeset{#2};$
            $U=\pgflabtypeset{#3};$
        }
    }
}
\clearpage
$A=\pgflabtypeset{A}$
\pgflabLUdecompose A as L times U

% find pivot in the specific column
% find pivot in the range (#1,#1) to (#1,end)
% does not check boundary
\def\pgflabfindpivotatcol #1 in #2{
    \pgfkeys{/lab/#2/h/.get=\pgflabh}
    \def\pgflabtempmax{-inf}
    \def\pgflabtempindex{0}
    \pgfplotsforeachungrouped\f@i in{#1,...,\pgflabh}{
        \pgfkeys{/lab/#2/\f@i/#1/.get=\pgflabtempentry}
        % compare the abs value
        \pgfmathsetmacro\pgflabtempentry{abs(\pgflabtempentry)}
        \pgfmathparse{\pgflabtempmax<\pgflabtempentry}
        % update if necessary
        \ifpgfmathfloatcomparison
            \let\pgflabtempmax\pgflabtempentry
            \let\pgflabtempindex\f@i
        \fi
    }
}
\clearpage
$A=\pgflabtypeset{A}$
find pivot at specific column: \par
\pgflabfindpivotatcol 1 in A
at col 1 it is \pgflabtempmax at row \pgflabtempindex \par
\pgflabfindpivotatcol 2 in A
at col 2 it is \pgflabtempmax at row \pgflabtempindex \par
\pgflabfindpivotatcol 3 in A
at col 2 it is \pgflabtempmax at row \pgflabtempindex

% A = PLU decomposition
% partial pivoting
\def\pgflabPLUdecompose #1 as #2 times #3 times #4{
    \pgfkeys{/lab/#1/h/.get=\pgflabH}
    \pgfkeys{/lab/#1/w/.get=\pgflabW}
    % decide the loop boundary
    \edef\pgflab@H-1{\the\numexpr\pgflabH-1}
    \ifnum\pgflab@H-1>\pgflabW
        \edef\pgflab@H-1{\pgflabW}
    \fi
    % set P as identity
    % set #2 as identity
    \pgflabneweyeof {\pgflabH} by {\pgflabH} as {#2}
    % set L as identity
    % set #3 as identity
    \pgflabneweyeof {\pgflabH} by {\pgflabH} as {#3}
    % copy A to U
    % copy #1 to #4
    \pgflabcopymatrix {#1} to {#4}
    % forget A, do job at P and L and U
    % forget #1, do job at #2 and #3 and #4
    \pgfplotsforeachungrouped\d@i in{1,...,\pgflab@H-1}{
        \pgflabfindpivotatcol {\d@i} in {#4}
        \pgflabswitchrow {\d@i} and row {\pgflabtempindex} in {#4}
        \pgflabswitchcol {\d@i} and col {\pgflabtempindex} in {#3}
        \pgflabswitchrow {\d@i} and row {\pgflabtempindex} in {#3}
        \pgflabswitchcol {\d@i} and col {\pgflabtempindex} in {#2}
        \par\medskip
        switch \d@i{} and \pgflabtempindex\par
        $P=\pgflabtypeset{#2};$
        $L=\pgflabtypeset{#3};$
        $U=\pgflabtypeset{#4};$
        \edef\d@@i+1{\the\numexpr\d@i+1}
        \pgfplotsforeachungrouped\d@j in{\d@@i+1,...,\pgflabH}{
            % use (\d@i,\d@i) to eliminate (\d@j,\d@i)
            \pgfkeys{/lab/#4/\d@i/\d@i/.get=\pgflabtempentrya}
            \pgfkeys{/lab/#4/\d@j/\d@i/.get=\pgflabtempentryb}
            \pgfmathsetmacro\pgflabtempratio{\pgflabtempentryb/\pgflabtempentrya}
            \pgflabaddcol {\d@i} by col {\d@j} times {\pgflabtempratio} in {#3}
            \pgflabaddrow {\d@j} by row {\d@i} times {-\pgflabtempratio} in {#4}
        }
        \par\medskip
        eliminate one column. check P and L and U \par
        $P=\pgflabtypeset{#2};$
        $L=\pgflabtypeset{#3};$
        $U=\pgflabtypeset{#4};$
    }
}
\pgflabread{A}{
    3   1   -7  5   0
    -9  -4  -8  -2  9
    4   -3  6   0   -1
    -5  8   2   -6  7
}
\bigskip
\pgflabPLUdecompose A as P times L times U

% find pivot in the specific column and row
% find pivot in the range (#1,#1) to (end,end)
% does not check boundary
\def\pgflabfindpivotafter#1 in #2{
    \pgfkeys{/lab/#2/h/.get=\pgflabh}
    \pgfkeys{/lab/#2/w/.get=\pgflabw}
    \def\pgflabtempmax{-inf}
    \def\pgflabtempindex{0}
    \def\pgflabtempjndex{0}
    \pgfplotsforeachungrouped\f@i in{#1,...,\pgflabh}{
        \pgfplotsforeachungrouped\f@j in{#1,...,\pgflabw}{
            \pgfkeys{/lab/#2/\f@i/\f@j/.get=\pgflabtempentry}
            % compare the abs value
            \pgfmathsetmacro\pgflabtempentry{abs(\pgflabtempentry)}
            \pgfmathparse{\pgflabtempmax<\pgflabtempentry}
            % update if necessary
            \ifpgfmathfloatcomparison
                \let\pgflabtempmax\pgflabtempentry
                \let\pgflabtempindex\f@i
                \let\pgflabtempjndex\f@j
            \fi
        }
    }
}

% A = PLUQ decomposition
% partial pivoting
\def\pgflabPLUQdecompose #1 as #2 times #3 times #4 times #5{
    \pgfkeys{/lab/#1/h/.get=\pgflabH}
    \pgfkeys{/lab/#1/w/.get=\pgflabW}
    % decide the loop boundary
    \edef\pgflab@H-1{\the\numexpr\pgflabH-1}
    \ifnum\pgflab@H-1>\pgflabW
        \edef\pgflab@H-1{\pgflabW}
    \fi
    % set P as identity
    % set #2 as identity
    \pgflabneweyeof {\pgflabH} by {\pgflabH} as {#2}
    % set L as identity
    % set #3 as identity
    \pgflabneweyeof {\pgflabH} by {\pgflabH} as {#3}
    % copy A to U
    % copy #1 to #4
    \pgflabcopymatrix {#1} to {#4}
    % set Q as identity
    % set #5 as identity
    \pgflabneweyeof {\pgflabW} by {\pgflabW} as {#5}
    % forget A, do job at P and L and U
    % forget #1, do job at #2 and #3 and #4
    \pgfplotsforeachungrouped\d@i in{1,...,\pgflab@H-1}{
        \pgflabfindpivotafter {\d@i} in #4

        \pgflabswitchrow {\d@i} and row {\pgflabtempindex} in {#4}
        \pgflabswitchcol {\d@i} and col {\pgflabtempindex} in {#3}
        \pgflabswitchrow {\d@i} and row {\pgflabtempindex} in {#3}
        \pgflabswitchcol {\d@i} and col {\pgflabtempindex} in {#2}
        {}
        \pgflabswitchcol {\d@i} and col {\pgflabtempjndex} in {#4}
        \pgflabswitchrow {\d@i} and row {\pgflabtempjndex} in {#5}

        switch (\d@i{},\d@i{}) and (\pgflabtempindex,\pgflabtempjndex) \par
        $P=\pgflabtypeset{#2};$
        $L=\pgflabtypeset{#3};$
        $U=\pgflabtypeset{#4};$
        $Q=\pgflabtypeset{#5};$
        \edef\d@@i+1{\the\numexpr\d@i+1}
        \pgfplotsforeachungrouped\d@j in{\d@@i+1,...,\pgflabH}{
            % use (\d@i,\d@i) to eliminate (\d@j,\d@i)
            \pgfkeys{/lab/#4/\d@i/\d@i/.get=\pgflabtempentrya}
            \pgfkeys{/lab/#4/\d@j/\d@i/.get=\pgflabtempentryb}
            \pgfmathsetmacro\pgflabtempratio{\pgflabtempentryb/\pgflabtempentrya}
            \pgflabaddcol {\d@i} by col {\d@j} times {\pgflabtempratio} in {#3}
            \pgflabaddrow {\d@j} by row {\d@i} times {-\pgflabtempratio} in {#4}
        }

        eliminate one column. check P and L and U and Q\par
        $P=\pgflabtypeset{#2};$
        $L=\pgflabtypeset{#3};$
        $U=\pgflabtypeset{#4};$
        $Q=\pgflabtypeset{#5};$
    }
}
\pgflabread{A}{
    3   -7  5   0   1   0   1
    -9  -8  -2  9   -1  9   -4
    4   6   0   -1  -2  -1  -3
    -5  2   -6  7   8   7   8
    -1  -2  -1  -3  4   6   0
    7   8   7   8   -5  2   -6
}
\bigskip
$A=\pgflabtypeset{A}$
\pgflabPLUQdecompose A as P times L times U times Q




















% new matrix with desired entry
% entry can contain \n@i and \n@j
\def\pgflabnewmatrixof #1 by #2 with #3 as #4{
    \def\pgflabh{#1}\pgfkeys{/lab/#4/h/.let=\pgflabh}
    \def\pgflabw{#2}\pgfkeys{/lab/#4/w/.let=\pgflabw}
    \pgfplotsforeachungrouped\n@i in{1,...,\pgflabh}{
        \pgfplotsforeachungrouped\n@j in{1,...,\pgflabw}{
            \pgfmathparse{#3}
            \pgfkeys{/lab/#4/\n@i/\n@j/.let=\pgfmathresult}
        }
    }
}
\clearpage
\pgflabnewmatrixof 10 by 10 with rand as C
\pgflabPLUQdecompose C as P times L times U times Q

% debug macro
% we can pass it to sage
% but we need to replace negative sign by ascii's -
\def\pgflabrawoutput#1{%
    \pgfkeys{/lab/#1/h/.get=\pgflabh}%
    \pgfkeys{/lab/#1/w/.get=\pgflabw}%
    matrix([%
    \pgfplotsforeachungrouped\t@i in{1,...,\pgflabh}{%
        [%
        \pgfplotsforeachungrouped\t@j in{1,...,\pgflabw}{%
            \pgfkeys{/lab/#1/\t@i/\t@j/.get=\pgflabtempentry}%
            \pgfmathparse{\pgflabtempentry}%
            \pgfmathfloattosci{\pgfmathresult}%
            \mbox{\pgfmathresult}%
            \ifnum\t@j<\pgflabw,\hskip1ptplus3pt\allowbreak\fi
        }%
        ]%
        \ifnum\t@i<\pgflabh,\hskip1ptplus3pt\allowbreak\fi
    }%
    ])%
}
\clearpage
C=\pgflabrawoutput{C};\par
P=\pgflabrawoutput{P};\par
L=\pgflabrawoutput{L};\par
U=\pgflabrawoutput{U};\par
Q=\pgflabrawoutput{Q};\par
(C-P*L*U*Q).norm()

\end{document}

debug mode

Old Answer

I would like to try

\documentclass{article}
\usepackage{pgfplotstable,mathtools}
\pgfplotsset{compat=newest}
\pgfkeys{/pgf/fpu}
\begin{document}



% we are lazy
% let pgfplotstable read the matrix
\pgfplotstableread[header=false]{
    8   1   6   8
    3   5   7   5
    4   9   2   7
}\matrixA

% we will store data by pgfleys
% create a handy handler
\pgfkeys{/handlers/.let/.code=\pgfkeyslet{\pgfkeyscurrentpath}{#1}}
% PS: /.initial is more like \def, but we want \xdef or \edef or \let

% but we also need some fast macros
\pgfplotstablegetrowsof\matrixA \xdef\matrixheight{\pgfplotsretval}\pgfkeys{/matrix/A/height/.let=\matrixheight}
\pgfplotstablegetcolsof\matrixA \xdef\matrixwidth{\pgfplotsretval} \pgfkeys{/matrix/A/width/.let=\matrixwidth}

% check data
Matrix $A$ is \pgfkeys{/matrix/A/height} by \pgfkeys{/matrix/A/width}.
In other words: \par Matrix $A$ is \matrixheight{} by \matrixwidth{}.



% store the entries into pgfkeys
\pgfplotsforeachungrouped\i in{1,...,\matrixheight}{
    \pgfplotsforeachungrouped\j in{1,...,\matrixwidth}{
        % since fpu is on, this is easier way to do 9+1
        \pgfplotstablegetelem{\the\numexpr\i-1}{\the\numexpr\j-1}\of\matrixA
        \pgfkeys{/matrix/A/\i/\j/.let=\pgfplotsretval}
    }
}

% check data
\bigskip The matrix entries are: \par
\pgfplotsforeachungrouped\i in{1,...,\matrixheight}{
    \pgfplotsforeachungrouped\j in{1,...,\matrixwidth}{
        \pgfkeys{/matrix/A/\i/\j}, 
    }
    ; \par
}



% define row operation: switch
\def\rowoperationswitch#1and#2 {
    \pgfplotsforeachungrouped\j in{1,...,\matrixwidth}{
        \pgfkeys{/matrix/A/#1/\j/.get=\tempmatrixentryA}
        \pgfkeys{/matrix/A/#2/\j/.get=\tempmatrixentryB}
        \pgfkeys{/matrix/A/#1/\j/.let=\tempmatrixentryB}
        \pgfkeys{/matrix/A/#2/\j/.let=\tempmatrixentryA}
    }
}

% try and check
\rowoperationswitch3and2
\bigskip After switching row3 and row2, the matrix entries are: \par
\pgfplotsforeachungrouped\i in{1,...,\matrixheight}{
    \pgfplotsforeachungrouped\j in{1,...,\matrixwidth}{
        \pgfkeys{/matrix/A/\i/\j}, 
    }
    ; \par
}



% define row operation: multiplication
\def\rowoperationmultiply#1by#2 {
    \pgfplotsforeachungrouped\j in{1,...,\matrixwidth}{
        \pgfkeys{/matrix/A/#1/\j/.get=\tempmatrixentry}
        \pgfmathparse{\tempmatrixentry*#2}
        \pgfkeys{/matrix/A/#1/\j/.let=\pgfmathresult}
    }
}

% try and check
\rowoperationmultiply3by9
\bigskip After multiplying row3 by 9, the matrix entries are: \par
\pgfplotsforeachungrouped\i in{1,...,\matrixheight}{
    \pgfplotsforeachungrouped\j in{1,...,\matrixwidth}{
        \pgfkeys{/matrix/A/\i/\j}, 
    }
    ; \par
}
remember: fpu is on! \par
Human readable version \par
\def\pgfmathprintmatrix{
    \pgfplotsforeachungrouped\i in{1,...,\matrixheight}{
        \indent
        \pgfplotsforeachungrouped\j in{1,...,\matrixwidth}{
            \pgfkeys{/matrix/A/\i/\j/.get=\tempmatrixentry}
            \pgfmathparse{\tempmatrixentry}
            \clap{\pgfmathprintnumber{\pgfmathresult}}\hskip20pt
        }
        \par
    }
}
\pgfmathprintmatrix



% define row operation: addition
\def\rowoperationadd#1by#2times#3 {
    \pgfplotsforeachungrouped\j in{1,...,\matrixwidth}{
        \pgfkeys{/matrix/A/#1/\j/.get=\tempmatrixentryA}
        \pgfkeys{/matrix/A/#2/\j/.get=\tempmatrixentryB}
        \pgfmathparse{\tempmatrixentryA+\tempmatrixentryB*#3}
        \pgfkeys{/matrix/A/#1/\j/.let=\pgfmathresult}
    }
}

% try and check
\rowoperationadd1by2times-1
\bigskip After adding row2 by row1 times -1, the matrix entries are: \par
\pgfmathprintmatrix



% We do RREF by hand
\pgfkeys{/pgf/number format/fixed}
\bigskip We do RREF by hand \par
add 2 by 1 times -1: \par
\rowoperationadd2by1times-1
\pgfmathprintmatrix

\medskip add 3 by 1 times -6.75: \par
\rowoperationadd3by1times-6.75
\pgfmathprintmatrix

\medskip add 3 by 2 times -5.8235: \par
\rowoperationadd3by2times-5.8235
\pgfmathprintmatrix



% renew A
\pgfplotsforeachungrouped\i in{1,...,\matrixheight}{
    \pgfplotsforeachungrouped\j in{1,...,\matrixwidth}{
        \pgfplotstablegetelem{\the\numexpr\i-1}{\the\numexpr\j-1}\of\matrixA % lazy~~
        \pgfkeys{/matrix/A/\i/\j/.let=\pgfplotsretval}
    }
}
\clearpage Restart with $A$ \par
\pgfmathprintmatrix



% Automatic RREF without row switching
% \I is different form \i
\xdef\matrixheightminusone{\the\numexpr\matrixheight-1}
\pgfplotsforeachungrouped\I in{1,...,\matrixheightminusone}{
    \pgfplotsforeachungrouped\J in{\I,...,\matrixheightminusone}{
        \xdef\J{\the\numexpr\J+1}
        \pgfkeys{/matrix/A/\I/\I/.get=\tempmatrixentryA}
        \pgfkeys{/matrix/A/\J/\I/.get=\tempmatrixentryB}
        \bigskip
        entry [\I][\I] is \tempmatrixentryA \par
        entry [\J][\I] is \tempmatrixentryB \par
        \pgfmathparse{-\tempmatrixentryB/\tempmatrixentryA}
        \xdef\temprowscaler{\pgfmathresult}
        \message{^^J^^J\I,\J,\temprowscaler^^J^^J}
        add row\J{} by row\I{} times \pgfmathprintnumber{\temprowscaler} \par
        \rowoperationadd\J by\I times{\temprowscaler}
        \pgfmathprintmatrix
    }
}



% renew A
\pgfplotsforeachungrouped\i in{1,...,\matrixheight}{
    \pgfplotsforeachungrouped\j in{1,...,\matrixwidth}{
        \pgfplotstablegetelem{\the\numexpr\i-1}{\the\numexpr\j-1}\of\matrixA % lazy~~
        \pgfkeys{/matrix/A/\i/\j/.let=\pgfplotsretval}
    }
}
\clearpage Restart with $A$ \par
\pgfmathprintmatrix



% maybe we need pivoting
\def\rowoperationfindpivot{
    % find the maximal element in this column
    \def\maxofthiscolumn{-inf}
    \def\maxofthiscolumnindex{0}
    \pgfplotsforeachungrouped\K in{\I,...,\matrixheight}{
        \pgfkeys{/matrix/A/\K/\I/.get=\tempmatrixentry}
        % compare
        \pgfmathparse{abs(\tempmatrixentry)}
        \let\tempmatrixabsentry\pgfmathresult
        \pgfmathparse{\maxofthiscolumn<\tempmatrixabsentry}
        % update if necessary
        \ifpgfmathfloatcomparison

            \let\maxofthiscolumn\tempmatrixabsentry
            \let\maxofthiscolumnindex\K
        \fi
    }
}
\xdef\I{1}
\rowoperationfindpivot
For column \I, the maximum is \pgfmathprintnumber{\maxofthiscolumn} at row \maxofthiscolumnindex



% Automatic RREF with partial pivot
\def\RREFwithpivoting{
    \pgfplotsforeachungrouped\I in{1,...,\matrixheightminusone}{
        \rowoperationfindpivot
        \rowoperationswitch\I and{\maxofthiscolumnindex}
        \bigskip
        For column \I, the maximum is \pgfmathprintnumber{\maxofthiscolumn} at row \maxofthiscolumnindex \par
        so we switch row\I{} and row\maxofthiscolumnindex, the matrix entries are: \par
        \pgfmathprintmatrix
        \pgfplotsforeachungrouped\J in{\I,...,\matrixheightminusone}{
            \xdef\J{\the\numexpr\J+1}
            \pgfkeys{/matrix/A/\I/\I/.get=\tempmatrixentryA}
            \pgfkeys{/matrix/A/\J/\I/.get=\tempmatrixentryB}
            \bigskip
            \pgfmathparse{-\tempmatrixentryB/\tempmatrixentryA}
            \xdef\temprowscaler{\pgfmathresult}
            add row\J{} by row\I{} times \pgfmathprintnumber{\temprowscaler} \par
            \rowoperationadd\J by\I times{\temprowscaler}
            \pgfmathprintmatrix
        }
    }
}
\RREFwithpivoting

.......

The rest is deleted because of the length limitation.