The Gaussian elimination method is one of the most important and ubiquitous algorithms that can help deduce important information about the given system of linear equations based on the augmented matrix:
 If there is a solution, the method will find it.
 If there is more than one solution, the method will find a general formula for all of them.
 If there are no solutions, the method will let us know.
This method has been historically around and used by ancient Babylonian and Chinese mathematicians since 179 CE. It was first described in the Nine Chapters of Mathematical Art. The author is known as Fang Cheng Shu, See the historical notes.
Many people have contributed to Gaussian elimination, including Isaac Newton. However, it was named by the American mathematician George Forsythe (19171972) in honor of the German mathematician and physicist Carl Friedrich Gauss (17771855). The full history of Gaussian algorithm can be found within the Joseph F. Grcar article.
Row Echelone Forms
The Gaussian elimination method is basically a series of operations carried out on a given matrix, in order to mathematically simplify it to its echelon form. When it is applied to solve a linear system A x = b, it consists of two steps: forward elimination (also frequently called Gaussian elimination procedure) to reduce the matrix to upper triangular form, and back substitution. Therefore, solving a linear system of algebraic equations using this elimination procedure can also be called forward elimination and back substitution and abbreviated as FEBS.0 
p_{1}


0 
p_{2}


0 
p_{3}


0 
p_{4}


0  
0 
....

p_{r}

 All nonzero rows are above any rows of all zeroes.
 Each leading entry, called the pivot, of a row is in a column to the right of the leading entry of the row above it.
 All entries in a column below a leading entry are zeroes.
 Any present zero rows are placed all the way at the bottom of the matrix.
 The very first value of the matrix (termed as the leading entry or pivot) is a nonzero term (some texts prefer it be ‘1’, but it is not necessarily).
Any nonzero matrix may be row reduced (that is, transformed by elementary row operations) into more than one matrix in echelon form, using different sequences of row operations. However, the leading entries are always in the same positions in any echelon form obtained from a given matrix.
Theorem: A linear system has a solution
if and only if the rightmost column of the associated augmented matrix is
not a pivot columnthat is, if and only if its echelon form does not
contain a row of the form
[ 0 0 ... 0 ⊚ ] with ⊚
nonzero.
If a linear system is consistent, then the solution set contains either a
unique solution (with no free variables) or infinitely many solutions, when
there is at least one free variable.
Row echelon form states that the Gaussian elimination method has been specifically applied to the rows of the matrix. It uses only those operations that preserve the solution set of the system, known as elementary row operations:
 Addition of a multiple of one equation to another.
Symbolically: (equation j) \( \mapsto \) (equation j) + k (equation i).  Multiplication of an equation by a nonzero constant k.
Symbolically: (equation j) \( \mapsto \) k (equation j).  Interchange of two equations>
Symbolically: (equation j) \( \Longleftrightarrow \) (equation i).
Row echelon and Reduced row echelon forms are the resulting matrices of the Gaussian elimination method. By applying Gaussian elimination to any matrix, one can easily deduce the following information about the given matrix:
 the rank of the matrix;
 the determinant of the matrix;
 the inverse (invertible square matrices only);
 the kernel vector (also known as ‘null vector’).
Theorem: If a homogeneous linear system has n unknowns, and if the row echelon form of its augmented matrix has r nonzero rows, then the system has nr free variables.
Theorem: Let \( {\bf B} = \left[ {\bf A} \, \vert \, {\bf b} \right] \) be the augmented matrix corresponding to the linear equation Ax=b, and suppose B is row equivalent (using a sequence of elementary row operations) to the new augmented matrix \( {\bf C} = \left[ {\bf U} \, \vert \, {\bf c} \right] , \) which corresponds to the linear system Ux=c. Then the two linear systems have precisely the same solution set.
Naive Gaussian elimination algorithm for Ax=b with a square matrix A (pseudocode)
for i=1 to n1
for j=i+1 to n
m=a(j,i)/a(i,i)
for k=i+1 to n
a(j,k) = a(j,k) m*a(i,k)
endfor
b(j) = b(j)  m*b(i)
endfor
endfor
The outmost loop (the i loop) ranges over the columns of the matrix;
the last column is skipped because we do not need to perform any eliminations
there because there are no elements below the diagonal. (if we were doing
elimination on a nonsquare matrix with more rows than columns, we would have
to include the last column in this loop.)
The middle loop (the j loop) ranges down the ith column, below the diagonal (hence j ranges only from i+1 to nthe dimension of the matrix A). We first compute the multiplier m, for each row. This is the constant by which we multiply the ith row in order to eliminate the a_{ji} element. Note that we overwrite the previous values with the new ones, and we do not actually carry out computation that makes a_{ji} zero. Also this loop is where the rightside vector is modified to reflect the elimination step.
The innermost loop (the k loop) ranges across the jth row, starting after the ith column, modifying each element appropriately to reflect the elimination of a_{ji}.
Finally, we must be aware that the algorithm does not actually create the zeroes in the lower triangular half of B; this would be wasteful of computer time since we don't need to have these zeroes in place for the algorithm to work. The algorithm works because we utilize only the upper triangular part of the matrix from this point forward, so the lower triangular elements need never be referenced.
This algorithm requires \( \frac{n}{3} \left( n^2 1 \right) \) ﬂoating point multiplications and divisions for operations on the coefﬁcient matrix and \( \frac{n}{2} \left( n 1 \right) \) multiplications for operations on the righthand terms, where after the triangular set has to be solved with \( \frac{n}{2} \left( n +1 \right) \) operations.
Backward solution algorithm for A (pseudocode)
x(n) = b(n)/a*n,n)
for i=n1 to 1
sum = 0
for j=i+1 to n
sum = sum + a(i,j)*x(j)
endfor
x(i) = (b(i)  sum)/a(i,i)
endfor
This algorithm simply matches backward up the diagonal, computing each
x_{i} in turn. Finally, we are computing
Some equations have to be interchanged if the corner elements \( A_{11}, \ A_{22}^{(1)} \) are not all zeroes to allow Gauss elimination to work. In the following, \( A_{ij}^{(n)} \) is the element after the nth iteration. One method is: if \( A_{kk}^{(k1)} = 0, \) than search for an element \( A_{pk}^{(k1)} \) with p > k that is not zero and interchange the pth and the nth equation. This strategy fails only if the set is singular and has no solution at all.
For example, consider the system of equations
Obviously one of the problems with the process of row reducing a matrix is the potential for human arithmetic errors. It is clear that computers can do this job much better than any human beeing. We can deepen our understanding of how the row reduction process works, and simultaneously eliminate arithmetic mistakes, by using a computer algebra system in a stepbystep fashion.
From computational point of view, the pivot should be the largest element. If we look at the first column, the largest entry is in the second row. To move it into the first row and choose "4" as a pivot, we need to swap first two rows. It is convenient, but not necessarily, to have 1 in this position, which we call pivot.
Mathematica stores vectors and matrices as lists and lists of lists, respectively. We enter the augmented matrix, say M, rowwise in Mathematica with the command
MatrixForm
or TraditionalForm
.
We can access row 1 of the matrix A with the syntax M[[1]], and similarly work with any other row of the matrix. This allows us to execute elementary row operations.
To work on the row reduction in the desired example, we ﬁrst want to swap rows 1 and 3; this is accomplished by using a temporary variable “temp” in the following way:
temp = A1[[1]];
A1[[1]] = A1[[3]];
A1[[3]] = temp; MatrixForm[A1]
A2[[2]] = 4*A1[[1]] + A1[[2]];
A2[[3]] = 2*A1[[1]] + A1[[3]];
MatrixForm[A2]
Another option is to choose "2" in the first row as pivot because it is not zero. We want to start with a 1 in the first position of the first column. To achieve it, we can either divide each entry in the first row by 2 or switch the first and the third rows because the last row has "1" in the first position already. If we stick with the latter decision, we may want to begin by scaling the first row by 1/2.
The first argument is the index of the row to be replaced, the second argument is the row to form a multiple of, and the final argument is the scale factor. Thus, M.add_multiple_of_row(n,m,A) would replace row n with (row n)+A*(row m).
Since the last entry in the first column is already zero we can move on to the second column. Our first step is to get a 1 in the second position of the second column. Normally we would do this by scaling the row but in this case we can swap the second and third row
Now we want to eliminate the 5 below the 1. This is done by multiplying the second row by −5 and adding it to the third row.
To get a 1 as the last entry in the third column we can scale the third row by −1/2
At this point, the matrix is in echelon form (well, having the 1's down the diagonal of the matrix is not required, but it does make our work easier). All that remains to find the solution is to put the matrix in reduced echelon form which requires that we replace all entries in the first three columns that are not on the main diagonal (where the 1's are) with zeros. We will start with the third column and work our way up and to the left. Remember that this is an augmented matrix and we are going to ignore the rightmost column; it just "goes along for the ride."
Therefore, our solution is complete. We see that the solution to our linear system is
\( x_1 =2, \ x_2 = 2, \ \mbox{and} \ x_3 = 3. \) ■
The following subroutine was prepared by Nasser M. Abbasi (see his website for possible updates: https://12000.org/my_notes/rref/index.htm).
Det[A]
M = Join[A, b, 2]
displayRREF
pivot now is (1,1)
R(2) > R(2) − (4) * R(1) \( \displaystyle \begin{pmatrix} 1 & 2& 3&0 \\ 0&3&6&3 \\ 7&8&9&6 \end{pmatrix} \)
R(3) > R(3) − (7) * R(1) \( \displaystyle \begin{pmatrix} 1 & 2& 3&0 \\ 0&3&6&3 \\ 0&6&12&6 \end{pmatrix} \)
pivot now is (2,2)
Scaling pivot row so that the pivot element is one. ===> \( \displaystyle \begin{pmatrix} 1 & 2& 3&0 \\ 0&1&2&1 \\ 0&6&12&6 \end{pmatrix} \)
R(3) > R(3) − (−6) * R(2) \( \displaystyle \begin{pmatrix} 1 & 2& 3&0 \\ 0&1&2&1 \\ 0&0&0&0 \end{pmatrix} \)
LinearSolve[A, c]

For each set of matrices, determine which elementary operations can transform
matrix A to matrix B.

\[ {\bf A} = \begin{bmatrix} 1 & 5 & 3 & 2\\ 2 & 2 & 3 & 5\\ 0 & 2 & 4 & 3\\ \end{bmatrix} \qquad {\bf B} = \begin{bmatrix} 1 & 5 & 3 & 2\\ 0 & 8 & 9 & 1\\ 0 & 2 & 4 & 3\\ \end{bmatrix} ; \]

\[ {\bf A} = \begin{bmatrix} 2 & 3 & 7 & 1\\ 0 & 4 & 3 & 5\\ 1 & 2 & 7 & 4\\ \end{bmatrix} \qquad {\bf B} = \begin{bmatrix} 2 & 3 & 7 & 1\\ 1 & 2 & 7 & 4\\ 0 & 4 & 3 & 5\\ \end{bmatrix} ; \]

\[ {\bf A} = \begin{bmatrix} 3 & 2 & 1 & 1\\ 4 & 6 & 5 & 2\\ 1 & 7 & 2 & 3\\ \end{bmatrix} \qquad {\bf B} = \begin{bmatrix} 0 & 19 & 7 & 10\\ 4 & 6 & 5 & 2\\ 1 & 7 & 2 & 3\\ \end{bmatrix} . \]


For the following problems, use row operations to solve the system of equations.

\begin{equation*} x_{1}+2\,x_{2}=3 \end{equation*} \begin{equation*} 2\,x_{1}+3\,x_{2}=4 \end{equation*}

\begin{equation*} 3\,x_{1}+x_{2}=6 \end{equation*} \begin{equation*} 3\,x_{1}+3\,x_{2}=2 \end{equation*}


In the following problems, solve the systems of linear equations given
their augmented matrices.

\begin{equation*} \begin{bmatrix} 1 & 1 & 2 & 1 &  & 5\\ 0 & 5 & 1 & 3 &  & 4\\ 0 & 0 & 3 & 1 &  & 5\\ 0 & 0 & 0 & 2 &  & 8\\ \end{bmatrix} \end{equation*}

\begin{equation*} \begin{bmatrix} 1 & 2 & 2 &  & 4\\ 0 & 3 & 2 &  & 1\\ 0 & 0 & 0 &  & 3\\ \end{bmatrix} \end{equation*}


For the following problems, solve the system of equations.

\begin{equation*} 2\,x_{2}+x_{3}=3 \end{equation*} \begin{equation*} 2\,x_{1}+x_{2} \,2x_{3}=3 \end{equation*} \begin{equation*} 4\,x_{1}+4\,x_{2}3\,x_{3}=3 \end{equation*}

\begin{eqnarray*} 3\,x_{1}2\,x_{2}1\,x_{3} &=1 , \\ 2\,x_{1}1\,x_{2}+2\,x_{3} &=4 , \\ x_{2}2\,x_{3} &=6. \end{eqnarray*}


For the following problems, find the value or values of f, if any,
which make the matrix an augmented matrix of a consistent linear system.

\[ \begin{bmatrix} 2 & f & 3\\ 4 & 8 & 6 \end{bmatrix} ; \]

\[ \begin{bmatrix} 1 & f & 3\\ 2 & 4 & 8 \end{bmatrix} ; \]

\[ \begin{bmatrix} 2 & 1 & f\\ 4 & 2 & 6\\ \end{bmatrix} . \]


For the following problems, use row operations to solve the system of equations.
 \begin{equation*} x_{1}+x_{2}x_{3}=9 \end{equation*} \begin{equation*} x_{2}+3x_{3}=3 \end{equation*} \begin{equation*} x_{1}2x_{3}=2 \end{equation*}
 \begin{equation*} 4x_{1}+8x_{2}4x_{3}=4 \end{equation*} \begin{equation*} 3x_{1}+8x_{2}+5_{3}=11 \end{equation*} \begin{equation*} 2x_{1}+x_{2}+12x_{3}=17 \end{equation*}
 \begin{equation*} x_{1}+2x_{2}+3x_{3}=4 \end{equation*} \begin{equation*} 2x_{1}+3x_{2}+3x_{3}=1 \end{equation*} \begin{equation*} 4x_{1}+2x_{2}+3x_{3}=1 \end{equation*}

Determine whether the following augmented matrices represent consistent
systems of equations and how many solutions exist.
 A= \begin{bmatrix} 1 & 4 & 2 &  & 8\\ 5 & 7 & 5 &  & 6\\ 3 & 2 & 6 &  & 6\\ \end{bmatrix}
 A= \begin{bmatrix} 1 & 1 & 1 &  & 3\\ 2 & 0 & 1 &  & 1\\ 0 & 1 & 1 &  & 0\\ \end{bmatrix}
 A= \begin{bmatrix} 1 & 1 & 1 &  & 1\\ 1 & 2 & 2 &  & 4\\ 3 & 2 & 1 &  & 12\\ \end{bmatrix}
 A= $\begin{bmatrix} 3 & 1 & 3 &  & 1\\ 1 & 2 & 2 &  & 1\\ 2 & 3 & 1 &  & 1\\ \end{bmatrix}$

Determine whether the system with complex coefficients is
consistent. If it is, find the general solution.
 \( \displaystyle \begin{split} z_1  {\bf j}\,z_2 &= 2  {\bf j} , \\ 3\,z_1 + \left( 1  2{\bf j} \right) z_2 &= 5  2 {\bf j} ; \end{split} \)
 \( \displaystyle \begin{split} \left( 1  2{\bf j} \right) z_1  {\bf j}\, z_2 &= 1 + 2{\bf j} , \\ \left( 1 + 3{\bf j} \right) z_1  \left( 1 + {\bf j} \right) z_2 &= {\bf j} 3 ; \end{split} \)
 \( \displaystyle \begin{split} \left( 1 + {\bf j} \right) z_1 + 2\,z_2 &= 1 + {\bf j} , \\ \left( 1  {\bf j} \right) z_1 + {\bf j} \, z_2 &= 1  {\bf j} ; \end{split} \)
 \( \displaystyle \begin{split} z_1 + 2\, z_2 + \left( 3 = 2{\bf j} \right) z_3 &= 5 + {\bf j} , \\ 2{\bf j}\, z_1 + 7{\bf j}\, z+2 + \left( 4 + 7 {\bf j} \right) z_3 &= 2 + 21 {\bf j} . \end{split} \)
 Anton, Howard (2005), Elementary Linear Algebra (Applications Version) (9th ed.), Wiley International
 Beezer, R., A First Course in Linear Algebra, 2015.
 Beezer, R., A Second Course in Linear Algebra, 2013.
 Casady, Chelsea, Gaussian Elimination and its History, Liberty University
 Davis, E., Linear Algebra and Probability for Computer Science Applications, CRC Press, 2012.
 Grcar, Joseph F., Mathematicians of Gaussian Elimination, Notes of American Mathematical Society, 58, No 6, 2011, 782792.
 Gaussian elimination for dummies
 Examples of Gaussian Elimination, Dartmouth College.
 Higham, Nicholas, Gaussian Elimination, Manchester Institute for Mathematical Sciences, School of Mathematics, The University of Manchester, 2011.
 Systems of Linear Equations
 Elementary Row Operations
 Elementary Row Operations for Matrices
 Automatic Row Reduction Input Form, A program that performs Gaussian elimination similarly to a human working on paper.
 Reduce Row Echelon Form on line by Kardi Teknomo.
 An online tool for performing elementary operations with rows and columns,
 GaussJordan Elimination Calculator.