How to find the steps in finding the reduced row echelon form of a matrix programmatically?

Question

I would like to figure out a way to see the row operations RowReduce uses to arrive at the reduced row echelon form of a matrix. I am aware of the step-by-step solutions in Wolfram Alpha and Wolfram Alpha Style Notebooks available from the File tab. I want to compute a programmatic sequence of row operations. I think there might be a way to use Reap and Sow to arrive at my goal. My grasp of Reap and Sow is weak, unfortunately.

A similar question covers row reduced function for square matrices. I would like to take any matrix and find the row echelon form.

Answer 1

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