Automated Gaussian elimination with partial pivoting/back-solving

First, we introduce a few Mathematica commands that we will use in the partial pivoting algorithm below.

mylist = {1.2, 3.4, 5.1, 0.3}

Max[mylist]

Position[mylist, 5.1]

Position[mylist, Max[mylist]]

a = {{1, 5, 2, 3}, {3, 2, 3, 5}, {-4, 3, 2, 
            1}};
columnlist = Abs[Table[a[[j, 1]], {j, 1, 3}]]

p = Position[columnlist, Max[columnlist]][[1, 1]]

The following implements Gaussian elimination with partial pivoting.  Print statements are used to display the sequence of matrices that result from the algorithm.

a = {{1, 15, 2, 3}, {3, 2, 3, 5}, {-4, 3, 2, 1}}; 
Print[a//MatrixForm];
n = Length[a]; 
at=Abs[Transpose[a]];
(* Gaussian elimination *)
Do[
	(* Find pivot index p using Max *)
	columnlist=Abs[Table[a[[j,i]],{j,i,n}]];
	p=Position[columnlist,Max[columnlist]][[1,1]]+(i-1);
	If[a[[p,i]]==0,Return[Print["unique solution not guaranteed"]]];
	(* Pivot if needed *)
	If[p != i, 
	   	temprow = a[[p]]; 
	    a[[p]] = a[[i]]; 
        a[[i]] = temprow;
        Print[a//MatrixForm]];
    (* After any pivot, do elimination *) 
    Do[m = a[[j,i]]/a[[i,i]]; 
    a[[j]] = a[[j]] - m*a[[i]], {j, i + 1, n}];
    Print[a//MatrixForm],
  {i, 1, n - 1}];
(* Back solving *)
If[a[[n,n]] == 0, Return["unique solution not guaranteed"]]; 
x = {}; 
x = Append[x, a[[n,n + 1]]/a[[n,n]]]; 
Do[x = Prepend[x, (a[[i,n + 1]] - Sum[a[[i,j]]*x[[j-n-1]], 
       {j, i + 1, n}])/a[[i,i]]], {i, n - 1, 1, -1}];
x

As a check, we can use Mathematica's RowReduce command.

RowReduce[a]