Here is an updated version which shows step by step the forward Gaussian Elimination steps for any size matrix.
This does not do the reduced echelon part (reverse direction). This can be easily added if needed.
Example 1
mat = {{1, 2, 1}, {-2, -3, 1}, {3, 5, 0}};
displayRREF[mat]

Example 2
mat = {{1, 2, 2, 4}, {1, 3, 3, 5}, {2, 6, 5, 6}};
displayRREF[mat]

Example 3
mat = {{-7, -6, -12, -33}, {5, 5, 7, 24}, {1, 0, 4, 5}};
displayRREF[mat]

Example 4
mat = {{1, -1, 2, 1}, {2, 1, 1, 8}, {1, 1, 0, 5}};
displayRREF[mat]

Example 5
mat = {{2, 1, 7, -7, 2}, {-3, 4, -5, -6, 3}, {1, 1, 4, -5, 2}};
displayRREF[mat]

Example 6
mat = {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}};
displayRREF[mat]

Code
Note, code was written in a code cell. Not an input cell. So formating can get messed up when using this back in an input cell. I've also put the notebook in case hard to copy this here
displayRREF[matIn_?(MatrixQ[#] &)] :=
Module[{mat = matIn, nRows,nCols,p, newp, tmp,scale},
(*version Nov 21, 2021*)
{nRow,nCol} = Dimensions[mat];
Print[">>>>>>Starting forward Gaussian elimination phase using ", MatrixForm[mat]];
Do[
If[p<nRow,Print["pivot now is (", p, ",", p, ")" ]];
If[ mat[[p,p]] === 0 && p<nRow,
newp = FirstPosition[mat[[p;;nRow, p]], _?(# != 0 &)];
If[ newp===Missing["NotFound"],
Print["Unable to continue. Can not find non-zero pivot"];
Return[Module];
];
newp = p+First@newp-1;
tmp = mat[[p,All]];
mat[[p,All]] = mat[[newp,All]];
mat[[newp,All]] = tmp;
Print["Since pivot is zero, then we exchange row ",newp," with row ",p,". ===> ",MatrixForm[mat]];
];
If[ mat[[p,p]] != 1,
newp = FirstPosition[mat[[ p;;nRow, p]], _?(Abs[#] == 1 &)];
If[ newp=!=Missing["NotFound"],
newp = p+First@newp-1;
tmp = mat[[p,All]];
mat[[p,All]] = mat[[newp,All]];
mat[[newp,All]] = tmp;
Print["Swapping row ",p," with row ",newp," ===> ",MatrixForm[mat]];
If[mat[[p,p]]==-1,
mat[[p,All]] =- mat[[p,All]];
Print["Scaling row ",p," by -1 ===>",MatrixForm[mat]]
]
,
If[mat[[p,p]]=!=0,
mat[[p,All]] = mat[[p,All]]/mat[[p,p]];
Print["Scaling pivot row so that the pivot element is one. ===> ",MatrixForm[mat]]
]
]
];
(*now elimination is done to zero all rows below the pivot row*)
Do[
scale = mat[[j,p]]*mat[[p,p]];
mat[[j,All]] = mat[[j,All]] - scale*mat[[p,All]];
Print["R(",j,") -> R(",j,") - (", scale,") * R(",p,") ===> ", MatrixForm[mat]]
,
{j, p+1, nRow}
]
,
{p, 1, nRow}
]
]
LUDdcomposition
@Nasser. I will note that LUD extends to rectangular matrices, e.g.In[9]:= LUDecomposition[{{1, 2, 3}, {4, 3, 1}}] Out[9]= {{{1, 2, 3}, {4, -5, -11}}, {1, 2}, 0}
. I do not know offhand whether the code at that link relies on the matrix being square though. $\endgroup$