Gauss Elimination with Pivoting

Concept of Pivoting

While performing the Gauss Elimination method for a linear system, we might encounter cases in which the pivot element is equal to zero or has a magnitude quite small to other coefficients. If the pivot element is zero then it introduces infinities in the system and if it has a small magnitude then it introduces floating-point precision error. Both cases are undesirable and need to be fixed.

In such cases, first, we search for a new pivot element as the absolute maximum of the coefficients in the current column, below the current row. Then we exchange the current row with the row containing the new pivot element and continue as previously. For example:

$$\begin{bmatrix} 0 & 1 & 1 \\ -4 & 3 & 5 \\ 3 & -2 & 7 \\ \end{bmatrix} \begin{bmatrix} x_1\\ x_2\\ x_3 \end{bmatrix} \begin{bmatrix} 2\\ 7\\ 9\\ \end{bmatrix}$$

In the above system, we start with the pivot element $a_{11}$. But since $a_{11}=0$, we need to exchange the pivot. For that the coefficients in the $1^{st}$ column below $1^{st}$ row are $-4,\ 3$. We take $-4$ as the new pivot element as $|-4|=4$ is the maximum. Thus we exchange $R_1 \leftrightarrow R_2$ and the new system becomes.

$$\begin{bmatrix} -4 & 3 & 5 \\ 0 & 1 & 1 \\ 3 & -2 & 7 \\ \end{bmatrix} \begin{bmatrix} x_1\\ x_2\\ x_3 \end{bmatrix} \begin{bmatrix} 7\\ 2\\ 9\\ \end{bmatrix}$$

From here we can continue as basic Gauss Elimination method. If we encounter zero pivot element in lower rows then we should again perform the pivoting.

Types of Pivoting

Partial Pivoting

In partial pivoting, whenever we encounter zero pivot element we exchange the row to maximize the magnitude of the pivot element. It involves only the exchange of rows and thus is the preferred method. We can also always maximize the magnitude of the pivot element without checking whether it is zero or not. Doing so will reduce the floating-point precision error but at the cost of performance.

Partial Pivoting Algorithm

The following implements Gauss Elimination with partial pivoting in which we always maximize the magnitude of the pivot element by row exchanges. Steps 6 and 7 are the additional steps.

1. Start
2. Declare and read the size of the system $N$
3. Declare variable for system, $S$ as a 2D array of size $N\times (N+1)$. Such that $S_{i,j}$ represents the element at $i^{th}$ row and $j^{th}$ column.
4. For each elements of $S$ input a value from user
5. Start from $k=1$
6. Find $m$ such that $S_{m,k}$ is maximum among $S_{k,k},\ S_{k+1,k},\ ...,\ S_{N,k}$
7. Exchange rows $R_{k}\leftrightarrow R_{m}$
8. Take $R_{k}$ as pivot row and then eliminate $x_{k}$ from the rows $R_{i}$ for $i=(k+1),...,N$ by applying $R_{i}\leftarrow R_{i}- \frac{S_{i,k}}{S_{k,k}}R_{k}$ modifying $S$
9. Increase $k$ by $1$. Goto step 6 if $k \text{<} N$
10. Find $x_{j}=\frac{S_{j,N+1}-\sum_{i=j+1}^{N}x_{i}S_{j,i}}{S_{j,j}}$ for $j=N,N-1,...,1$
11. Stop

Partial Pivoting Pseudo code

Read N
Declare S = Array(N, N+1)

-- Read the input
for i=1:N
  for j=1:N+1
    Read S[i][j]
  endfor
endfor

-- Perform elimination
for k=1:N-1
  -- Find the maximum row
  maxMag = 0
  maxRow = k
  for m=k:N
    if abs(S[m][k])>maxMag:
      maxMag = abs(S[m][k])
      maxRow = m
    end
  end

  -- Exchange the rows
  S[maxRow], S[k] = S[k], S[maxRow]

  for i=k+1:N
    factor = S[i][k]/S[k][k]
    for j=k+1:N+1
      S[i][j] = S[i][j] - factor * S[k][j]
    end
  end
end

-- Perform back substitution
for k=N:1
  sum = 0
  for i=k+1:N
    sum += S[k][i]
  end
  Print x[k] = (S[k][N+1]-sum)/S[k][k]
end

Partial Pivoting Program

C

/*
 * Gauss Elimination Method with Partial Pivoting
*/

#include <stdio.h>
#include<stdlib.h>
#include<math.h>
#define MAX_SYSTEM_SIZE 20

int main()
{
  // The size of linear system
  int N = 3;

  // Preallocate enough memory for the system
  double S[MS][MS+1];
  double X[MS];

  // Get the size of the system
  printf("Enter the size of linear system: ");
  scanf("%d", &N);

  // Limit the max system size
  if (N >= MS)
  {
    printf("Please enter size below %d", MS);
    return -1;
  }

  // Ask for the system
  printf("\nEnter the matrix in augmented matrix form:\n");
  for (int i = 0; i < N; i++)
  {
    for (int j = 0; j < N + 1; j++)
    {
      printf("Enter S[%d][%d] = ", i + 1, j + 1);
      scanf("%lf", &S[i][j]);
    }
  }

  // Perform Elimination
  for (int k = 0; k < N - 1; k++)
  {
    // Find out the maximum magnitude pivot
    double maxMag = 0;
    int maxRow = k;
    for(int m=k;m<N;m++){
      if(fabs(S[m][k])>maxMag){
        maxMag = fabs(S[m][k]);
        maxRow = m;
      }
    }

    // Exchange the rows
    for(int j=0;j<N+1;j++){
      double temp = S[k][j];
      S[k][j] = S[maxRow][j];
      S[maxRow][j] = temp;
    }

    for (int i = k + 1; i < N; i++)
    {
      // Precalculate factor as S[i][k] will 
      // change in the following loop
      double factor = S[i][k] / S[k][k];
      for (int j = k; j < N + 1; j++)
      {
        S[i][j] = S[i][j] - factor * S[k][j];
      }
    }
  }

  // Perform Back substitution
  for (int k = N - 1; k >= 0; k--)
  {
    double sum = 0;
    for (int i = k + 1; i < N; i++)
    {
      sum += S[k][i] * X[i];
    }
    X[k] = (S[k][N] - sum) / S[k][k];
  }

  printf("\nSolution for the system is: \n");
  for (int i = 0; i < N; i++)
  {
    printf("x%d=%lf\n", i + 1, X[i]);
  }

  return 0;
}

Python

#
# Gauss Elimination Method With Partial Pivoting
#

N = int(input("Enter the size of the system "))
S = [[0]*(N+1) for _ in range(N)]
X = [0]*N

# Ask for input
print("Enter the system in augmented matrix form:")
for i in range(N):
  for j in range(N+1):
    S[i][j] = float(input(f"Enter S[{i+1}][{j+1}] = "))

# Perform elimination
for k in range(N-1):
  # Find the max magnitude pivot element
  maxMag = 0
  maxRow = k
  for m in range(k, N):
    if abs(S[m][k]) > maxMag:
      maxMag = abs(S[m][k])
      maxRow = m

  # Exchange the rows
  S[k], S[maxRow] = S[maxRow], S[k]

  for i in range(k+1, N):
    # Precalculate factor as S[i][k] will change in the following loop
    factor = S[i][k]/S[k][k]
    for j in range(k, N+1):
      S[i][j] = S[i][j] - factor*S[k][j]

# Perform back substitution
for k in reversed(range(0, N)):
  sum = 0
  for i in range(k+1, N):
    sum += S[k][i]*X[i]
  
  X[k] = (S[k][N]-sum)/S[k][k]

# Print the solution
for i in range(N):
  print(f"x{i+1} = {X[i]}")

Complete Pivoting

In complete pivoting, we maximize the pivot element along the row as well as the column. This involves the exchange of both rows and columns. When we exchange the column, the order of unknowns also needs to change. Complete pivoting results in the lowest precision error. However, it is rarely used due to the higher complexity of implementation and slower execution.