nRows[mat_] := Dimensions[mat][[1]]

nColumns[mat_] := Dimensions[mat][[2]]

addColumn[mat_?MatrixQ, col_List] /; nRows[mat] == Length[col] :=
    Transpose[ Append[Transpose[mat], col] ]

eliminateColumn[ mat_?MatrixQ, i_Integer ] /; 1 <= i <= nColumns[mat] :=
    Module[{res = mat, m = nRows[mat], j, c},
        Do[ c = -res[[j, i]]/res[[i, i]];
            res[[j]] = res[[j]] + c res[[i]]
          , {j, i+1, m} ];
        res
    ]

GaussianElimination[ mat_?MatrixQ ] :=
    Module[{res = mat, m = Min[nRows[mat], nColumns[mat]], i},
        Do[ res = eliminateColumn[res, i], {i, 1, m-1} ];
        res
    ]

backSubstitution[mat_?MatrixQ] /; nColumns[mat] == 1 + nRows[mat] :=
    Module[{m = nRows[mat], i, k, x},
        x = Table[0, {m}];
        Do[ x[[i]] = (mat[[i, -1]] -
                      Sum[ x[[k]] mat[[i, k]], {k, i+1, m} ])/
                     mat[[i, i]];
          , {i, m, 1, -1} ];
        x
    ]

linearSystem[ a_?MatrixQ, b_?VectorQ ] /; nRows[a] === nColumns[a] === Length[b] :=
    backSubstitution[ GaussianElimination[ addColumn[a, b] ] ]
