# Find Elementary Matrices that produce RREF

I have a matrix, M, which I am reducing to RREF using RowReduce[M]. I would like to know the elementary matrices which perform the row reducing operations on M during RowReduce[M] to reduce it to RREF, or at least the product of them all.

Updated to Solve Ax=b and display all transformations to reduced Echelon form, and simplified the API. This now also prints the matrix inverse as by product, since it now uses an full augmented matrix.

The idea is to use LUDecomposition first to obtain the permutation rows used for pivoting. Once the order of the rows used is known, then forward elimination is used to generate Echelon form, then backward elimination is used to produce the final reduced echelon form and the solution vector.

Full augmented matrix is used so that the RHS of the augmented matrix will contain the matrix inverse at the end. Verified same inverse is produced as Mathematica Inverse. This only works on matrices that have non-zero determinant.

Here are some usage examples

## Example 1

(*solve Ax=b *)
mat = {{2, 7, 3}, {1, 3, 2}, {3, 7, 9}};
b = {11, 2, -12};
displayRREF[mat, b] ## Example 2

mat={{17,42,-36},{13,45,-34},{12,47,-35}};
b={213,226,197};
displayRREF[mat,b] ## Example 3

mat={{5,2,18,4},{0,1,2,5},{4,1,12,6},{2,3,8,9}};
b={1,2,3,4};
displayRREF[mat,b] ## Function

displayRREF[mat_?(MatrixQ[#] &), b_?(VectorQ[#] &)] :=
Module[{i, j, multiplier, pivot, augmented, m = Length@mat, lu, p, c, tmp, inverse},
(*version 3/10/2017*)

(*check if matrix is singular.Per Daniel Lichtblau post seen Wolfram site*)
(*this is better method than using Det*)
If[MatrixQ[mat] &&
MatrixRank[mat] == Length[mat] == Length[mat[]] === False,
Return["Sorry, but matrix is singular!"]
];

If[Length@b != Length@mat,
Return[
"Size of b vector not the same as number of rows in A matrix"]
];

{lu, p, c} = LUDecomposition[mat];
tmp = lu SparseArray[{i_, j_} /; j >= i -> 1, {Length@mat,
Length@mat}];
(*Print["Mathematica says Echelon form is ", MatrixForm@tmp];*)
tmp = mat[[p, All]];
(*Print["Mathematica says inverse Matrix is", MatrixForm@Inverse@
mat];*)
augmented = Join[mat, Transpose[{b}], 2];
augmented = ArrayFlatten[{{augmented, IdentityMatrix[Length@tmp]}}];

Print[">>>>>>Starting forward Gaussian elimination phase using ",
augmented[[All, 1 ;; m + 1]] //
matWithDiv[m + 1, Background -> LightOrange] ,
MatrixForm[augmented[[All, m + 2 ;;]]]];
Do[
Print["pivot now is (", pivot, ",", pivot, ")" ];
Do[
multiplier = augmented[[j, pivot]]/augmented[[pivot, pivot]];
Print["will now zero out element (", j, ",", pivot,
") by subtracting ", multiplier, " times row ", pivot,
" from row ", j];
augmented[[j, pivot ;;]] =
augmented[[j, pivot ;;]] -
multiplier*augmented[[pivot, pivot ;;]];
Print[
augmented[[All, 1 ;; m + 1]] //
matWithDiv[m + 1, Background -> LightOrange] ,
MatrixForm[augmented[[All, m + 2 ;;]]]]
, {j, pivot + 1, m}
]
, {pivot, 1, m}
];

Print[">>>>>>Starting backward elimination phase"];

Do[
Do[
multiplier = augmented[[j, pivot]]/augmented[[pivot, pivot]];
Print["will now zero out element (", j, ",", pivot,
") by subtracting ", multiplier, " times row ", pivot,
" from row ", j];
augmented[[j, pivot ;;]] =
augmented[[j, pivot ;;]] -
multiplier*augmented[[pivot, pivot ;;]];
Print[
augmented[[All, 1 ;; m + 1]] //
matWithDiv[m + 1, Background -> LightOrange] ,
MatrixForm[augmented[[All, m + 2 ;;]]]]
, {j, 1, pivot - 1}
]
, {pivot, 2, m}
];

Print[">>>>>>Starting Final phase, convert reduced echelon to     identity matrix"];
Do[
augmented[[j, ;;]] = augmented[[j, ;;]]/augmented[[j, j]]
, {j, 1, m}
];

Print[augmented[[All, 1 ;; m + 1]] //
matWithDiv[m + 1, Background -> LightOrange] ,
MatrixForm[augmented[[All, m + 2 ;;]]]];
(*flip at inverse Matrix now using p from LUDecomposition,
but using column wise*)
Print["Inverse Matrix is ",
MatrixForm[ augmented[[All, m + 2 ;;]] ]];
Print["Solution  vector is ", MatrixForm[augmented[[All, m + 1]]]]
];