[Tex/LaTex] 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

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.