# Row reduction macro

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. :*

• rref is numerically an ill conditioned algorithm. It becomes even worse with TeX precision. The industry standard is LU factorization. Commented Mar 11, 2017 at 16:12
• @percusse Good point. But I just want something that can handle the simple matrices the students are required to row reduce in the class I teach. Know what I mean? Commented Mar 11, 2017 at 16:18
• Doesn't matter, you'll get wrong results in even 3x3 matrices. multiply first row with 2/3 and add to the second row even that would screw up all integer valued entries. Commented Mar 11, 2017 at 16:29
• The sagetex package relies on the computer algebra system SAGE to calculate the row echelon form of a matrix, even if it has variables and gives you access to a lot of other commands on matrices. However, it doesn't give the steps to get the answer.
– DJP
Commented Mar 11, 2017 at 21:44
• Are you sure you need everything to be done inside TeX itself? Why not rely on an external program for doing all the computation (and communicate with the program from TeX maybe)? Commented Mar 22, 2017 at 17:07

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.

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):
[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
\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
\repeat
\ifnum\MAT_cnta<\MAT_cntc
}\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
\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
\repeat
}%
\def\MATexchangerowspecial #1#2#3{%#1>#2, only columns <#2 need update
\MAT_cnta=#2\relax
\MAT_cntb=\xint_c_ % 0
\xintloop
\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
\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
\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}


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
\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
\repeat
\MATrrdisplaymatrix\@U\par
\fi
\MAT_cntc\MATrr@pivrow
\xintloop
\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
\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


• It is important to stress relative to PLUQ (and even to Row Reduction) that I am doing exact computations, thus there is no need for the algorithm to look for the largest possible pivot, with the attendant fact that entries of L would be smaller than one in absolute value. Thus I do not follow here the algorithm fashionable in numerical analysis (by the way, it may in exceptional cases have stability problems, the U matrix containing exponentially big entries respective to the size of the matrices). I just pick up the first available pivot irrespective of its size.
– user4686
Commented Mar 25, 2017 at 14:49
• don't you want Bruhat decomposition of invertible matrices ? that will be a challenge to numerical analysis, but no problem for exact computations.
– user4686
Commented Mar 25, 2017 at 14:57
• by the way Maple's LUDecomposition operates like this: For Matrices with... at least one floating-point entry, pivots ... according to absolute magnitude. For Matrices with ... rational... entries, pivots ... the first nonzero element in the current column.. and indeed this is what I did here. One can check that with p, l, u := LUDecomposition(A); one does get p=P, l = L, and u = U.Q, where P, L, U, Q are computed as above. (Maple16)
– user4686
Commented Mar 25, 2017 at 18:01
• At first glance, this looks great. I'll need to give it a more thorough inspection before I reward you with the 500 reputation and risqué pics though. ;) Commented Mar 25, 2017 at 22:30

# 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}
\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};$
}
}
}
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.

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.

• 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 . 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}
\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
\pgfutil@ifnextchar[
}

% #1: optional option
% #2: a name of the matrix... usually A
\edef\pgflabname{#2}
}

% 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
\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
}
}
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
}
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}$

% 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}$

% 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};$
}
}
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};$
}
}
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

I would like to try

\documentclass{article}
\usepackage{pgfplotstable,mathtools}
\pgfkeys{/pgf/fpu}
\begin{document}

% we are lazy
% let pgfplotstable read the matrix
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
\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

\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
\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
\pgfmathprintmatrix

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

\medskip add 3 by 2 times -5.8235: \par
\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
\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
\pgfmathprintmatrix
}
}
}
\RREFwithpivoting

.......


The rest is deleted because of the length limitation.

• By the way, the product of pivots is -862566.82. And by Wolfram|Alpha the determinant of the first six columns is 862575. So the error is .0009% Commented Mar 22, 2017 at 21:09
• Hi, I would like to learn how to write matrix-related code like yours and the one of the other answer. Which of the following is more convenient to learn to do this? Programing in Tex or in LuaLaTex? Commented Dec 11, 2021 at 16:02
• If you can use LuaLaTeX use it. Commented Dec 11, 2021 at 23:04