更新-2

Question 1

这个答案有一些宏摘自https://tex.stackexchange.com/a/143035/4686。我对一些内部数据结构不太满意，但我决定保留它。

这https://tex.stackexchange.com/a/143035/4686计算行列式、逆等等...可以精确计算或者通过浮点运算计算。

在这里我重点关注精确的计算。矩阵条目可能是整数、分数、小数或科学计数法，但它们的处理方式准确。因此，这里不存在数值不稳定的问题。关于输入，格式是用分号分隔的行，用逗号分隔系数。

最后的编辑改进了一些内部方面，提供了一个更好的例子A=PLUQ，并重做了行减少的初始例子，用于显示截断的、而不是四舍五入的十进制扩展，因为它们后面带有点。

PLUQ 答案

例如，代码使用 TeX 排版，并将最终结果输出到 Maple 矩阵符号的文件中。

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]]);

现在我们可以复制粘贴到 Maple 中并检查确实如此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]

请注意，在 PLUQ 分解中，仅在必要时，P 和 Q 才会与我的代码一起出现。

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

行减少（初始）答案

我通过编辑稍微改进了代码的一些内部方面。

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

条目可能是十进制数，例如37.156。

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

示例（当我添加点时，我使用截断而不是四舍五入）：

\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

这个答案有一些宏摘自https://tex.stackexchange.com/a/143035/4686。我对一些内部数据结构不太满意，但我决定保留它。

这https://tex.stackexchange.com/a/143035/4686计算行列式、逆等等...可以精确计算或者通过浮点运算计算。

在这里我重点关注精确的计算。矩阵条目可能是整数、分数、小数或科学计数法，但它们的处理方式准确。因此，这里不存在数值不稳定的问题。关于输入，格式是用分号分隔的行，用逗号分隔系数。

最后的编辑改进了一些内部方面，提供了一个更好的例子A=PLUQ，并重做了行减少的初始例子，用于显示截断的、而不是四舍五入的十进制扩展，因为它们后面带有点。

PLUQ 答案

例如，代码使用 TeX 排版，并将最终结果输出到 Maple 矩阵符号的文件中。

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]]);

现在我们可以复制粘贴到 Maple 中并检查确实如此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]

请注意，在 PLUQ 分解中，仅在必要时，P 和 Q 才会与我的代码一起出现。

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

行减少（初始）答案

我通过编辑稍微改进了代码的一些内部方面。

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

条目可能是十进制数，例如37.156。

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

示例（当我添加点时，我使用截断而不是四舍五入）：

\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

Question 2

更新-2

我听到有人说吉文斯轮换。

% 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

对于 10×10 的随机矩阵，范数A-二维码约为 4e-4。范数QQᵀ-我约为 2e-4。

更新-1：新答案

我实现了三种分解：

A=鲁
A=附加值（即部分旋转）
A=普鲁士（即完全旋转）

如果A是米经过n，然后磷，大号是米经过米；乌是相同的A; 和问是n经过n。

优点

访问矩阵元素的复杂度是哦（1）（\csname假设哦（1））。因此分解的复杂度为哦（米²n）。
输入利用\pgfplotstableread来自pgfplotstable。因此它接受内联表、文件、已加载的表，甚至是创建的表\pgfplotstablenew。您还可以向其传递选项。（例如过滤）
输出利用\pgfplotstabletypeset了同一个包。或者您可以将矩阵转换回表格并执行任何您想做的事情。
计算由完成\pgfmathparse。我假设 FPU 已打开。但可以重新实现。
有一个调试宏可以输出矩阵的原始数据。您可以将这些数据复制并粘贴到任何现代矩阵计算器中。
根据维基百科，即使部分旋转在实践中也是数值稳定的。我测试了一个 10 x 10 的随机矩阵，并检查A-普鲁士在 Sage 中；标准约为 1.1e-6。（这与 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}

调试模式

旧答案

我想试一试

\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

.......

由于篇幅限制，其余内容被删除。

Answer

更新-2

我听到有人说吉文斯轮换。

% 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

对于 10×10 的随机矩阵，范数A-二维码约为 4e-4。范数QQᵀ-我约为 2e-4。

更新-1：新答案

我实现了三种分解：

A=鲁
A=附加值（即部分旋转）
A=普鲁士（即完全旋转）

如果A是米经过n，然后磷，大号是米经过米；乌是相同的A; 和问是n经过n。

优点

访问矩阵元素的复杂度是哦（1）（\csname假设哦（1））。因此分解的复杂度为哦（米²n）。
输入利用\pgfplotstableread来自pgfplotstable。因此它接受内联表、文件、已加载的表，甚至是创建的表\pgfplotstablenew。您还可以向其传递选项。（例如过滤）
输出利用\pgfplotstabletypeset了同一个包。或者您可以将矩阵转换回表格并执行任何您想做的事情。
计算由完成\pgfmathparse。我假设 FPU 已打开。但可以重新实现。
有一个调试宏可以输出矩阵的原始数据。您可以将这些数据复制并粘贴到任何现代矩阵计算器中。
根据维基百科，即使部分旋转在实践中也是数值稳定的。我测试了一个 10 x 10 的随机矩阵，并检查A-普鲁士在 Sage 中；标准约为 1.1e-6。（这与 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}

调试模式

旧答案

我想试一试

\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

.......

由于篇幅限制，其余内容被删除。

更新-2

答案1

PLUQ 答案

行减少（初始）答案

答案2

更新-2

更新-1：新答案

优点

旧答案

相关内容