I wanted to know which of the two methods gives more accurate results when solving linear equations problems, Gauss elimination or LU factorization? And what happens if we use partial pivoting while applying the first method, is this more accurate than the LU factorization?
[Math] Gauss elimination vs LU factorization
linear algebranumerical methods
Related Solutions
Gaussian Elimination helps to put a matrix in row echelon form, while Gauss-Jordan Elimination puts a matrix in reduced row echelon form. For small systems (or by hand), it is usually more convenient to use Gauss-Jordan elimination and explicitly solve for each variable represented in the matrix system. However, Gaussian elimination in itself is occasionally computationally more efficient for computers. Also, Gaussian elimination is all you need to determine the rank of a matrix (an important property of each matrix) while going through the trouble to put a matrix in reduced row echelon form is not worth it to only solve for the matrix's rank.
EDIT: Here are some abbreviations to start off with: REF = "Row Echelon Form". RREF = "Reduced Row Echelon Form."
In your question, you say you reduce a matrix A to a diagonal matrix where every nonzero value equals 1. For this to happen, you must perform row operations to "pivot" along each entry along the diagonal. Such row operations usually involve multiplying/dividing by nonzero scalar multiples of the row, or adding/subtracting nonzero scalar multiples of one row from another row. My interpretation of REF is just doing row operations in such a way to avoid dividing rows by their pivot values (to make the pivot become 1). If you go through each pivot (the numbers along the diagonal) and divide those rows by their leading coefficient, then you will end up in RREF. See these Khan Academy videos for worked examples.
In a system $Ax=B$, $x$ can only be solved for if $A$ is invertible. Invertible matrices have several important properties. The most useful property for your question is that their RREF is the identity matrix (a matrix with only 1's down the diagonal and 0's everywhere else). If you row-reduce a matrix and it does not become an identity matrix in RREF, then that matrix was non-invertible. Non-invertible matrices (also known as singular matrices) are not as helpful when trying to solve a system exactly.
Rule of Thumb/TL;DR: When doing calculations using floating point numbers (such as double, single, and float data types in many common programming languages), use partial pivoting unless you know you're safe without it and complete pivoting only when you know that you need it.
Longer Explanation: A square matrix $A$ has an $LU$ factorization (without pivoting) if, and only if, no zero is encountered in the pivot position when computing an $LU$ factorization of $A$. However, when does computations using floating point numbers a pivot that is nearly zero can lead to dramatic rounding errors. The simple workaround is to always permute the rows of the matrix such that the largest nonzero entry in a column is chosen as the pivot entry. This ensures that a nearly zero is never chosen. Complete pivoting goes even further by using row and column permutations to select the largest entry in the entire matrix as the pivot entry.
The above paragraph is a lose, intuitive picture of why pivoting is necessary. One can also prove sharp error bounds by carefully tracking the propagation of errors throughout the entire $LU$ factorization calculation. One way of structuring this error bound is a so-called backward error estimate. In a backwards error estimate for solving a linear system of equations $Ax = b$, one bounds the perturbation $E$ necessary to make the computed solution $\hat{x}$ produced by doing Gaussian elimination followed by back substitution an exact solution to a nearby system of linear equations $(A+E)\hat{x} = b$. A backwards error estimate can be very revealing if, for example, the matrix $A$ is determined by measurements of an engineering system with some error tolerance. If the entries of $A$ are only known $\pm 1\%$ and the backwards error is less than $0.1\%$, then the numerical errors made during our computations are smaller than the measurement errors and we've done a good job. TL;DR the quantity $E$ is a reasonable quantity to measure the error in an $LU$ factorization and resulting linear solve.
For Gaussian elimination without pivoting, the backwards error can be arbitrarily bad. Fortunately, for partial pivoting, the backwards error $E$ can be bounded as $\|E\|_\infty \le 6n^3 \rho \|A\|_\infty u + \mbox{higher order terms in $u$}$ ${}^\dagger$. Here, $\|\cdot\|_\infty$ is the operator matrix $\infty$-norm and $u$ is the unit roundoff which quantifies the accuracy of floating point computations ($u \approx 10^{-16}$ for IEE double precision arithmetic). The quantity $\rho$ is known as the growth factor for partial pivoting. While it is possible for $\rho$ to be as large as $2^{n-1}$, in practice $\rho$ usually grows very modestly with $n$.${}^\dagger$ Concerning the fact that $\rho$ grows very modestly in most applications, legendary numerical analyst Kahan wrote "Intolerable pivot-growth [with partial pivoting] is a phenomenon that happens only to numerical analysts who are looking for that phenomenon."${}^\$$
Nonetheless, one can write down matrices for which partial pivoting will fail to give an accurate answer due to an exponential growth factor $\rho = 2^{n-1}$. Wilkinson showed that the growth factor for complete pivoting is much smaller in the worst case $\rho \le n^{1/2} (2\cdot 3^{1/2}\cdots n^{1/n-1})^{1/2}$. In practice, $\rho$ for complete pivoting is almost always less than $100$.${}^\dagger$
After having delved into the details, you can see that this is a subtle business, and even in the classical field of numerical linear algebra there is somewhat of a gap between theory and experiment. In general, Gaussian elimination with partial pivoting is very reliable. Unless you know you can get away without pivoting (symmetric positive definite and diagonally dominant matrices are notable examples), one should use partial pivoting to get an accurate result. (Or compensate with something clever. For example, SuperLU uses "static pivoting" where one does "best guess" pivoting before starting the $LU$ factorization and then does no pivoting during the factorization. The loss of accuracy of this approach is compensated by using a few steps of iterative refinement.)
If partial pivoting isn't accurate enough, one can move to using complete pivoting instead for its lower growth factor. As user3417 points out, there are other ways of solving $Ax = b$ other than using $LU$ factorization-based approaches and these may be faster and more accurate than Gaussian elimination with complete pivoting. For example, $QR$ factorization runs in the $O(n^3)$ operations and has no growth factor. There may be special cases when one truly does want to use an $LU$ factorization: for example, a Gaussian elimination-based approach can be used to construct structure-preserving factorizations of a Cauchy-like matrix. In this case, complete pivoting (or its close cousin rook pivoting) may be the best approach.
${}^\dagger$ Reference Golub and Van Loan's Matrix Computations Fourth Edition Chapter 3.4
${}^\$$ Quoted in Higham's Accuracy and Stability of Numerical Algorithms Second Edition Chapter 9
Best Answer
I think the book Numerical Linear Algebra by Trefethen and Bau has a pretty good explanation of this topic. I've summarized some of the content from chapters 20 and 21.
Gaussian Elimination is unstable
Consider the matrix $$ A = \begin{bmatrix} 10^{-20} & 1 \\ 1 & 1\end{bmatrix} $$
This has exact LU decomposition, $$ L = \begin{bmatrix} 1 & 0 \\ 10^{20} & 1\end{bmatrix},~~~~ U = \begin{bmatrix} 10^{-20} & 1 \\ 0 & 1-10^{20}\end{bmatrix} $$
However, suppose we make small (relative) rounding errors and end up representing $1-10^{20}$ as $-10^{20}$. Then, $$ \begin{bmatrix} 1 & 0 \\ 10^{20} & 1\end{bmatrix}\begin{bmatrix} 10^{-20} & 1 \\ 0 & -10^{20}\end{bmatrix} = \begin{bmatrix} 10^{-20} & 1 \\ 1 & 0 \end{bmatrix} $$ which is far from $A$.
LU is Gaussian Elimination without Pivoting
To expand on my comment, when you do Gaussian Elimination without pivoting, each step can be written in terms of matrix multiplication on the left by a lower triangular matrix. The end result is an upper triangular matrix so written in matrix form Gaussian elimination looks something like $$ L_{m-1}\cdots L_2L_1 A = U $$
The product of lower triangular matrices is still lower triangular, as is the inverse, so setting $L^{-1} = L_{m-1}\cdots L_2L_1$ gives the LU decomposition.
As a concrete example, let's consider the matrix $$ A = \begin{bmatrix} 2 & 1 & 1\\ 4 & 3 & 2\\ 8 & 7 & 9 \end{bmatrix} $$
The first step of Gaussian elimination would be to subtract twice the first row from the second, and 4 times the first from the third. We can write this as, $$ L_1 = \begin{bmatrix} 1 \\ -2 & 1\\ -4 & & 1 \end{bmatrix} $$ You can verify that, $$ L_1A = \begin{bmatrix} 2 & 1 & 1 \\ 0 & 1 & 0 \\ 0 & 3 & 5 \end{bmatrix} $$
The next step of Gaussian elimination is to subtract 3 times the second row from the third row. Again, we can write this as, $$ L_2 = \begin{bmatrix} 1 \\ & 1 \\ & -3 & 1 \end{bmatrix} $$ and verify that, $$ L_2(L_1A) = \begin{bmatrix} 2 & 1 & 1 \\ 0 & 1 & 0 \\ 0 & 0 & 5 \end{bmatrix} $$
It turns out that the invese of $L_i$ is quite simple, you just negate the entries below the main diagonal. Similarly, the product of such matrices is also simple, you just add the entries below the diagonal. You can again verify that, $$ L = (L_2L_1)^{-1} = L_1^{-1}L_2^{-1} = \begin{bmatrix} 1 \\ 2 & 1\\ 4 & & 1 \end{bmatrix} \begin{bmatrix} 1 \\ & 1 \\ & 3 & 1 \end{bmatrix} = \begin{bmatrix} 1 \\ 2 & 1 \\ 4 & 3 & 1 \end{bmatrix} $$
PLU is Gaussian Elimination with Pivoting
Now, if you pivot at each step, you can view this as swapping the rows of the working matrix. This is the same as left multiplication by a permutation matrix. So doing gaussian elimination with pivoting will give a set of operations like $$ L_{m-1}P_{m-1}\cdots L_2P_2L_1P_1 A = U $$
If we define, $$ L_k' = P_{m-1}\cdots P_{k+1}L_kP_{k+1}^{-1}\cdots P_{m-1}^{-1} $$ you will see we can write the above series of operations as $$ (L_{m-1}'\cdots L_2'L_1')(P_{m-1}^{-1}\cdots P_1^{-1}) A = U $$
Which you can rearrange to get a PLU decomposition.