(* by Nasser M. Abbasi *)
(****** Version Oct 15, 2023. Fixed display again and improved row exchange*)
(****** Version Oct 14, 2023. Add displayREF which does forward elimination only*)
(****** Version Oct 10, 2023 of displayRREF. Add flag so pivot not normalized to one and kept as is*)

(*****Version July 17, 2023 of displayRREF. Add flag to turn on/off display of matrices*)
(*bug reports are welcome*)

SetDirectory[NotebookDirectory[]];
CurrentValue[$FrontEnd,{"PrintAction"}]={"PrintToNotebook"};
CurrentValue[$FrontEnd,{"MessageOptions","KernelMessageAction"}]={"PrintToNotebook","Beep"};

(*updated august 20, 2023 to fix a minor bug*)

displayREF[Ain_?(MatrixQ[#]&),displayMat_:True,normalizePivot_:True]:=Module[
{multiplier,j,i,pivotRow,pivotCol,nRows,nCols,p,tmp,startAgain,A=Ain,n,m,pivotsFound={},keepEliminating,nIter,entry,debug=False},
    Print[MatrixForm[A]];
    {nRows,nCols} = Dimensions[A];
    
    keepEliminating=True;
    
    n=1; m=1;
    nIter = 0;
    While[keepEliminating,
        nIter++;
        If[nIter>100, (*safe guard*)
             Return["Internal error. Something went wrong. Or very large system?",Module]
        ];
        
        If[debug,Print["Row = ",n," column = ",m]];
        
        If[m==nCols,
             keepEliminating=False 
        ,
             (*If[displayMat,Print@makeNiceMatrix[A,{n,m}]];*)   
             If[A[[n,m]] =!= 0,
                  Print["Pivot is A(",n,",",m,")"];
                  If[displayMat,Print@makeNiceMatrix[A,{n,m}]];
                  If[normalizePivot,
                        If[A[[n,m]] =!= 1,
                            A[[n,All]] = A[[n,All]]/A[[n,m]];
                            A=Simplify[A];
                            Print["Making the pivot 1 using using row(",n,")= row(",n,")/A(",n,",",m,")"];
                            If[displayMat,Print@makeNiceMatrix[A,{n,m}]];
                        ]
                 ];
                 If[n<nRows,
                       Do[
                            If[A[[j,m]] =!= 0,
                                  multiplier = A[[j,m]]/A[[n,m]];
                                  Print["Zeroing out element A(",j,",",m,") using row(",j,")=",multiplier,"*row(",n,")-row(",j,")"];
                                  A[[j,m;;]] = A[[j,m;;]] - multiplier*A[[n,m;;]];
                                  A=Simplify[A];
                                  If[displayMat,Print@makeNiceMatrix[A,{n,m}]]
                            ]
                       ,{j,n+1,nRows}
                       ];
                  ];
                  pivotsFound = AppendTo[pivotsFound,{n,m}];
                  If[n==nRows,
                      keepEliminating=False
                  ,
                      n++;
                      If[m<nCols,m++]
                  ]
            ,
                  (*pivot is zero. If we can find non-zero pivot row below, then exchange rows*)
                  Print["Entry {",n,",",m,") is zero. Looking to find non-zero pivot below in order to exchange rows"];
                  If[n==nRows&&m==nCols, 
                        If[debug,Print["keepEliminating is set to false since reached end of matrix"]];
                        keepEliminating=False
                  ,
                       (*If we can find non-zero pivot row below, then exchange rows*)
                       If[n<nRows,
                             If[debug,Print["Calling FirstPosition[A[[n+1;;,m]],_?(# =!=0)&]"]];
                             If[debug,Print["Where A=",A," and n+1=",n+1," and m=",m]];
                             p=SparseArray[A[[n+1;;,m]]]["NonzeroPositions"];
                             If[debug,Print["p=",p]];
                             If[Length[p]==0,
                                  If[m<nCols,
                                      If[debug,Print["p===Missing. So m<nCols, making column now from m=",m," to ",m+1]];
                                      m++
                                  ,
                                      keepEliminating=False
                                  ]
                             ,
                                  (*found non zero pivot below. Exchange rows*)
                                  p=First@p[[1]];
                                  
                                  Print["Found non zero pivot"];
                                  tmp = A[[n,All]];
                                  A[[n,All]] = A[[p+n,All]];
                                  A[[p+n,All]] = tmp;
                                  A=Simplify[A];
                                  Print["Exchanging row(",n,") and row(",p+n,")"];
                                  If[displayMat,Print@makeNiceMatrix[A,{n,m}]]
                            ]
                      ,
                            If[m<nCols,
                                  m++
                            ,
                                  keepEliminating=False
                            ]
                      ]
                  ]
             ]
        ]
    ];
    
    {A,pivotsFound}
];

(*----------------------------------------------------------*)
displayRREF[Ain_?(MatrixQ[#]&),displayMat_:True,normalizePivot_:True]:=Module[
 {multiplier,j,i,pivotRow,pivotCol,nRows,nCols,p,tmp,startAgain,A=Ain,n,m,pivotsFound,keepEliminating,nIter,entry,debug=False},
     {A,pivotsFound} = displayREF[Ain,displayMat,normalizePivot];
        
    If[debug,Print["Completed forward"]];
    If[debug,Print["pivotsFound",pivotsFound]];
    (*If[displayMat,Print@makeNiceMatrix[A,{n,m}]];*)
    Print[">>>>>>Starting backward elimination phase. The pivots are ",pivotsFound];
    Do[ 
        pivotRow=First@entry;
        pivotCol=Last@entry;
        If[debug,Print["pivotRow= ",pivotRow]];
        If[debug,Print[" pivotCol= ", pivotCol]];
       
        If[pivotRow>1,
            Do[
                 If[ A[[i,pivotCol]] =!= 0,
                      Print["Zeroing out element A(",i,",",pivotCol,") using row(",i,")=row(",i,")-A(",i,",",pivotCol,")*row(",pivotRow,")"];
                      A[[i,;;]] = A[[i,;;]] - A[[i,pivotCol]]*A[[pivotRow,;;]];
                      A=Simplify[A];
                      If[displayMat,Print@makeNiceMatrix[A,{pivotRow,pivotCol}]]
                 ]
            ,
            {i,pivotRow-1,1,-1} 
           ]
        ]
        ,
        {entry,pivotsFound}    
   ];
   
   {A,pivotsFound}
];

(*----------------------------------------------------------*)
makeSolutionSpecialCase[A_?(MatrixQ[#]&),b_?(VectorQ[#]&),pivotCols_List]:=Module[
   {nRows,nCols,nLeadingVariables,nFreeVariables,n,m,k,variables={},eq,freeVariables,sol={}},
   
    Print["Pivot columns are ",MatrixForm[pivotCols]];
    
    ClearAll[x,t]; (*did not make them local, to prevent $ from showing in print*)
    {nRows,nCols} = Dimensions[A];
    nLeadingVariables = Length[pivotCols];
    nFreeVariables =  nCols-nLeadingVariables;    
    Print["There are ",nLeadingVariables," leading variables and ",nFreeVariables," free variables. These are "];
    
    Array[t, nFreeVariables];
    Array[x, nCols];
    
    m=0;k=0;
    Do[
        If[Not[MemberQ[pivotCols,n]],
            m++;
            Print[x[n]," is a free variable. Let ",x[n],"=",t[m]];
            AppendTo[variables,t[m]];
            AppendTo[sol,0]
        ,
            Print[x[n]," is a leading variable"];
            AppendTo[variables,x[n]];
            AppendTo[sol,b[[++k]]]
        ]
    ,{n,1,nCols}
    ];
    
    freeVariables=(t[#]&/@Range[nFreeVariables]);
    Print["Hence the system after RREF is the following>>>>>>"];
    Print[MatrixForm[A.variables],"=",MatrixForm[b]];
    Print["There is different solution for different value of the free variables."];
    Print["Setting free variable ", freeVariables," to zero gives"];
    variables=variables/.((t[#]->0)&/@Range[nFreeVariables]);
    Print[MatrixForm[A.variables],"=",MatrixForm[b]];
    Print["Therefore the final solution is "];
    Print[MatrixForm[x[#]&/@Range[nCols]],"=",MatrixForm[sol]]
   
]

makeNiceMatrix[mat_?MatrixQ,pivot_List]:=Module[{g,nRow,nCol},
    {nRow,nCol} = Dimensions[mat];
     g = Grid[mat,Frame->{None,None,{pivot->True}}];
   
    MatrixForm[{{g}}]
]
(*thanks to http://mathematica.stackexchange.com/questions/60613/how-to-add-a-vertical-line-to-a-matrix*)
(*makes a dash line inside Matrix*)
Format[matWithDiv[n_,opts:OptionsPattern[Grid]][m_?MatrixQ]]:=MatrixForm[{{Grid[m,opts,Dividers->{n->{Red,Dashed}}]}}];

Example 1  REF only with no normalization

M={{1,2,3,0},{4,5,6,3},{7,8,9,6}};
displaymat=True;normalizePivot=False;
{result,pivots}=displayREF[M,displaymat,normalizePivot]

Example 2  REF only with no normalization

M={{0,2,3,0},{4,5,6,3},{7,8,9,6}};
displaymat=True;normalizePivot=False;
{result,pivots}=displayREF[M,displaymat,normalizePivot]

Example 3  REF only with normalization. Row exchange example

M={{0,2,3,0},{4,5,6,3},{7,8,9,6}};
displaymat=True;normalizePivot=True;
{result,pivots}=displayREF[M,displaymat,normalizePivot]

Example 4  displayRREF with normalization

mat={{2,3,2,1},{0,1,2,-3},{3,3,3,2},{-1,1,1,-4}};
displaymat=True;
normalizePivot=True;
{result,pivots}=displayRREF[mat,displaymat,normalizePivot]

Example 5 displayRREF with no normalization

mat={{2,3,2,1},{0,1,2,-3},{3,3,3,2},{-1,1,1,-4}};
displaymat=True;
normalizePivot=False;
{result,pivots}=displayREF[mat,displaymat,normalizePivot]

Example 6 displayRREF with normalization

mat={{1,2,2,4},{1,3,3,5},{2,6,5,6}};
displaymat=True; normalizePivot=True;
{result,pivots}=displayRREF[mat]

Example 7 displayRREF with normalization

mat={{-7,-6,-12,-33},{5,5,7,24},{1,0,4,5}};
displaymat=True; normalizePivot=True;
{result,pivots}=displayRREF[mat,displaymat,normalizePivot]

Example 8 displayRREF on symbolic matrix

ClearAll[s];
mat={{s,Sqrt[s],3,10,s^2},{1,s,2*s,10,1},{0,2,3,10,s^9},{1/s,5,3,8,2+s}};
displaymat=True; normalizePivot=True;
{result,pivots}=displayRREF[mat,displaymat,normalizePivot]

Example 9 displayRREF with normalization, row exchange is needed

mat={{0,-6,-12,-33},{5,5,7,24},{1,0,4,5}};
displaymat=True; normalizePivot=True;
{result,pivots}=displayRREF[mat,displaymat,normalizePivot]

RowReduce[mat]