Data Modelling & AI Data Structure & Algorithm

Gaussian Elimination to Solve Linear Equations

21 June 2025

0

The article focuses on using an algorithm for solving a system of linear equations. We will deal with the matrix of coefficients. Gaussian Elimination does not work on singular matrices (they lead to division by zero).

Input: For N unknowns, input is an augmented 
       matrix of size N x (N+1). One extra 
       column is for Right Hand Side (RHS)
  mat[N][N+1] = {{3.0, 2.0,-4.0, 3.0},
                {2.0, 3.0, 3.0, 15.0},
                {5.0, -3, 1.0, 14.0}
               };
Output: Solution to equations is:
        3.000000
        1.000000
        2.000000
Explanation:
Given matrix represents following equations
3.0X₁ + 2.0X₂ - 4.0X₃ =  3.0
2.0X₁ + 3.0X₂ + 3.0X₃ = 15.0
5.0X₁ - 3.0X₂ +    X₃ = 14.0
There is a unique solution for given equations, 
solutions is, X₁ = 3.0, X₂ = 1.0, X₃ = 2.0,

Row echelon form: Matrix is said to be in r.e.f. if the following conditions hold:

The first non-zero element in each row, called the leading coefficient, is 1.
Each leading coefficient is in a column to the right of the previous row leading coefficient.
Rows with all zeros are below rows with at least one non-zero element.

matrix1

Reduced row echelon form: Matrix is said to be in r.r.e.f. if the following conditions hold –

All the conditions for r.e.f.
The leading coefficient in each row is the only non-zero entry in its column.

The algorithm is majorly about performing a sequence of operations on the rows of the matrix. What we would like to keep in mind while performing these operations is that we want to convert the matrix into an upper triangular matrix in row echelon form. The operations can be:

Swapping two rows
Multiplying a row by a non-zero scalar
Adding to one row a multiple of another

The process:

Forward elimination: reduction to row echelon form. Using it one can tell whether there are no solutions, or unique solution, or infinitely many solutions.
Back substitution: further reduction to reduced row echelon form.

Algorithm:

Partial pivoting: Find the kth pivot by swapping rows, to move the entry with the largest absolute value to the pivot position. This imparts computational stability to the algorithm.
For each row below the pivot, calculate the factor f which makes the kth entry zero, and for every element in the row subtract the fth multiple of the corresponding element in the kth row.
Repeat above steps for each unknown. We will be left with a partial r.e.f. matrix.

Below is the implementation of the above algorithm.

C++

// C++ program to demonstrate working of Gaussian Elimination
// method
#include<bits/stdc++.h>
using namespace std;
 
#define N 3        // Number of unknowns
 
// function to reduce matrix to r.e.f.  Returns a value to 
// indicate whether matrix is singular or not
int forwardElim(double mat[N][N+1]);
 
// function to calculate the values of the unknowns
void backSub(double mat[N][N+1]);
 
// function to get matrix content
void gaussianElimination(double mat[N][N+1])
{
    /* reduction into r.e.f. */
    int singular_flag = forwardElim(mat);
 
    /* if matrix is singular */
    if (singular_flag != -1)
    {
        printf("Singular Matrix.\n");
 
        /* if the RHS of equation corresponding to
           zero row  is 0, * system has infinitely
           many solutions, else inconsistent*/
        if (mat[singular_flag][N])
            printf("Inconsistent System.");
        else
            printf("May have infinitely many "
                   "solutions.");
 
        return;
    }
 
    /* get solution to system and print it using
       backward substitution */
    backSub(mat);
}
 
// function for elementary operation of swapping two rows
void swap_row(double mat[N][N+1], int i, int j)
{
    //printf("Swapped rows %d and %d\n", i, j);
 
    for (int k=0; k<=N; k++)
    {
        double temp = mat[i][k];
        mat[i][k] = mat[j][k];
        mat[j][k] = temp;
    }
}
 
// function to print matrix content at any stage
void print(double mat[N][N+1])
{
    for (int i=0; i<N; i++, printf("\n"))
        for (int j=0; j<=N; j++)
            printf("%lf ", mat[i][j]);
 
    printf("\n");
}
 
// function to reduce matrix to r.e.f.
int forwardElim(double mat[N][N+1])
{
    for (int k=0; k<N; k++)
    {
        // Initialize maximum value and index for pivot
        int i_max = k;
        int v_max = mat[i_max][k];
 
        /* find greater amplitude for pivot if any */
        for (int i = k+1; i < N; i++)
            if (abs(mat[i][k]) > v_max)
                v_max = mat[i][k], i_max = i;
 
        /* if a principal diagonal element  is zero,
         * it denotes that matrix is singular, and
         * will lead to a division-by-zero later. */
        if (!mat[k][i_max])
            return k; // Matrix is singular
 
        /* Swap the greatest value row with current row */
        if (i_max != k)
            swap_row(mat, k, i_max);
 
 
        for (int i=k+1; i<N; i++)
        {
            /* factor f to set current row kth element to 0,
             * and subsequently remaining kth column to 0 */
            double f = mat[i][k]/mat[k][k];
 
            /* subtract fth multiple of corresponding kth
               row element*/
            for (int j=k+1; j<=N; j++)
                mat[i][j] -= mat[k][j]*f;
 
            /* filling lower triangular matrix with zeros*/
            mat[i][k] = 0;
        }
 
        //print(mat);        //for matrix state
    }
    //print(mat);            //for matrix state
    return -1;
}
 
// function to calculate the values of the unknowns
void backSub(double mat[N][N+1])
{
    double x[N];  // An array to store solution
 
    /* Start calculating from last equation up to the
       first */
    for (int i = N-1; i >= 0; i--)
    {
        /* start with the RHS of the equation */
        x[i] = mat[i][N];
 
        /* Initialize j to i+1 since matrix is upper
           triangular*/
        for (int j=i+1; j<N; j++)
        {
            /* subtract all the lhs values
             * except the coefficient of the variable
             * whose value is being calculated */
            x[i] -= mat[i][j]*x[j];
        }
 
        /* divide the RHS by the coefficient of the
           unknown being calculated */
        x[i] = x[i]/mat[i][i];
    }
 
    printf("\nSolution for the system:\n");
    for (int i=0; i<N; i++)
        printf("%lf\n", x[i]);
}
 
// Driver program
int main()
{
    /* input matrix */
    double mat[N][N+1] = {{3.0, 2.0,-4.0, 3.0},
                          {2.0, 3.0, 3.0, 15.0},
                          {5.0, -3, 1.0, 14.0}
                         };
 
    gaussianElimination(mat);
 
    return 0;
}

Java

// Java program to demonstrate working of Gaussian Elimination
// method
import java.io.*;
class GFG
{
 
  public static int N = 3; // Number of unknowns
 
  // function to get matrix content
  static void gaussianElimination(double mat[][])
  {
 
    /* reduction into r.e.f. */
    int singular_flag = forwardElim(mat);
 
    /* if matrix is singular */
    if (singular_flag != -1) 
    {
      System.out.println("Singular Matrix.");
 
      /* if the RHS of equation corresponding to
               zero row  is 0, * system has infinitely
               many solutions, else inconsistent*/
      if (mat[singular_flag][N] != 0)
        System.out.print("Inconsistent System.");
      else
        System.out.print(
        "May have infinitely many solutions.");
 
      return;
    }
 
    /* get solution to system and print it using
           backward substitution */
    backSub(mat);
  }
 
  // function for elementary operation of swapping two
  // rows
  static void swap_row(double mat[][], int i, int j)
  {
    // printf("Swapped rows %d and %d\n", i, j);
 
    for (int k = 0; k <= N; k++)
    {
      double temp = mat[i][k];
      mat[i][k] = mat[j][k];
      mat[j][k] = temp;
    }
  }
 
  // function to print matrix content at any stage
  static void print(double mat[][])
  {
    for (int i = 0; i < N; i++, System.out.println())
      for (int j = 0; j <= N; j++)
        System.out.print(mat[i][j]);
    System.out.println();
  }
 
  // function to reduce matrix to r.e.f.
  static int forwardElim(double mat[][])
  {
    for (int k = 0; k < N; k++) 
    {
 
      // Initialize maximum value and index for pivot
      int i_max = k;
      int v_max = (int)mat[i_max][k];
 
      /* find greater amplitude for pivot if any */
      for (int i = k + 1; i < N; i++)
        if (Math.abs(mat[i][k]) > v_max) 
        {
          v_max = (int)mat[i][k];
          i_max = i;
        }
 
      /* if a principal diagonal element  is zero,
             * it denotes that matrix is singular, and
             * will lead to a division-by-zero later. */
      if (mat[k][i_max] == 0)
        return k; // Matrix is singular
 
      /* Swap the greatest value row with current row
             */
      if (i_max != k)
        swap_row(mat, k, i_max);
 
      for (int i = k + 1; i < N; i++)
      {
 
        /* factor f to set current row kth element
                 * to 0, and subsequently remaining kth
                 * column to 0 */
        double f = mat[i][k] / mat[k][k];
 
        /* subtract fth multiple of corresponding
                   kth row element*/
        for (int j = k + 1; j <= N; j++)
          mat[i][j] -= mat[k][j] * f;
 
        /* filling lower triangular matrix with
                 * zeros*/
        mat[i][k] = 0;
      }
 
      // print(mat);        //for matrix state
    }
 
    // print(mat);            //for matrix state
    return -1;
  }
 
  // function to calculate the values of the unknowns
  static void backSub(double mat[][])
  {
    double x[]
      = new double[N]; // An array to store solution
 
    /* Start calculating from last equation up to the
           first */
    for (int i = N - 1; i >= 0; i--)
    {
 
      /* start with the RHS of the equation */
      x[i] = mat[i][N];
 
      /* Initialize j to i+1 since matrix is upper
               triangular*/
      for (int j = i + 1; j < N; j++) 
      {
 
        /* subtract all the lhs values
                 * except the coefficient of the variable
                 * whose value is being calculated */
        x[i] -= mat[i][j] * x[j];
      }
 
      /* divide the RHS by the coefficient of the
               unknown being calculated */
      x[i] = x[i] / mat[i][i];
    }
 
    System.out.println();
    System.out.println("Solution for the system:");
    for (int i = 0; i < N; i++)
    {
      System.out.format("%.6f", x[i]);
      System.out.println();
    }
  }
 
  // Driver program
  public static void main(String[] args)
  {
 
    /* input matrix */
    double mat[][] = { { 3.0, 2.0, -4.0, 3.0 },
                      { 2.0, 3.0, 3.0, 15.0 },
                      { 5.0, -3, 1.0, 14.0 } };
 
    gaussianElimination(mat);
  }
}
 
// This code is contributed by Dharanendra L V.

Python3

# Python3 program to demonstrate working of
# Gaussian Elimination method
N = 3
 
# function to get matrix content
def gaussianElimination(mat):
 
    # reduction into r.e.f.
    singular_flag = forwardElim(mat)
 
    # if matrix is singular
    if (singular_flag != -1):
 
        print("Singular Matrix.")
 
        # if the RHS of equation corresponding to
        #   zero row  is 0, * system has infinitely
        #   many solutions, else inconsistent*/
        if (mat[singular_flag][N]):
            print("Inconsistent System.")
        else:
            print("May have infinitely many solutions.")
 
        return
 
    # get solution to system and print it using
    #   backward substitution
    backSub(mat)
 
# function for elementary operation of swapping two rows
def swap_row(mat, i, j):
 
    for k in range(N + 1):
 
        temp = mat[i][k]
        mat[i][k] = mat[j][k]
        mat[j][k] = temp
 
# function to reduce matrix to r.e.f.
def forwardElim(mat):
    for k in range(N):
       
        # Initialize maximum value and index for pivot
        i_max = k
        v_max = mat[i_max][k]
 
        # find greater amplitude for pivot if any
        for i in range(k + 1, N):
            if (abs(mat[i][k]) > v_max):
                v_max = mat[i][k]
                i_max = i
 
        # if a principal diagonal element  is zero,
        # it denotes that matrix is singular, and
        # will lead to a division-by-zero later.
        if not mat[k][i_max]:
            return k    # Matrix is singular
 
        # Swap the greatest value row with current row
        if (i_max != k):
            swap_row(mat, k, i_max)
 
        for i in range(k + 1, N):
 
            # factor f to set current row kth element to 0,
            # and subsequently remaining kth column to 0 */
            f = mat[i][k]/mat[k][k]
 
            # subtract fth multiple of corresponding kth
            # row element*/
            for j in range(k + 1, N + 1):
                mat[i][j] -= mat[k][j]*f
 
            # filling lower triangular matrix with zeros*/
            mat[i][k] = 0
 
        # print(mat);        //for matrix state
 
    # print(mat);            //for matrix state
    return -1
 
# function to calculate the values of the unknowns
def backSub(mat):
 
    x = [None for _ in range(N)]    # An array to store solution
 
    # Start calculating from last equation up to the
    #  first */
    for i in range(N-1, -1, -1):
 
        # start with the RHS of the equation */
        x[i] = mat[i][N]
 
        # Initialize j to i+1 since matrix is upper
        #  triangular*/
        for j in range(i + 1, N):
           
            # subtract all the lhs values
            # except the coefficient of the variable
            # whose value is being calculated */
            x[i] -= mat[i][j]*x[j]
 
        # divide the RHS by the coefficient of the
        #  unknown being calculated
        x[i] = (x[i]/mat[i][i])
 
    print("\nSolution for the system:")
    for i in range(N):
        print("{:.8f}".format(x[i]))
 
# Driver program
 
# input matrix
mat = [[3.0, 2.0, -4.0, 3.0], [2.0, 3.0, 3.0, 15.0], [5.0, -3, 1.0, 14.0]]
gaussianElimination(mat)
 
# This code is contributed by phasing17

C#

// C# program to demonstrate working 
// of Gaussian Elimination method
using System;
 
class GFG{
     
// Number of unknowns
public static int N = 3; 
 
// Function to get matrix content
static void gaussianElimination(double [,]mat)
{
 
    /* reduction into r.e.f. */
    int singular_flag = forwardElim(mat);
     
    /* if matrix is singular */
    if (singular_flag != -1) 
    {
        Console.WriteLine("Singular Matrix.");
         
        /* if the RHS of equation corresponding to
           zero row  is 0, * system has infinitely
           many solutions, else inconsistent*/
        if (mat[singular_flag,N] != 0)
            Console.Write("Inconsistent System.");
        else
            Console.Write("May have infinitely " + 
                          "many solutions.");
         
        return;
    }
     
    /* get solution to system and print it using
    backward substitution */
    backSub(mat);
}
 
// Function for elementary operation of swapping two
// rows
static void swap_row(double [,]mat, int i, int j)
{
     
    // printf("Swapped rows %d and %d\n", i, j);
     
    for(int k = 0; k <= N; k++)
    {
        double temp = mat[i, k];
        mat[i, k] = mat[j, k];
        mat[j, k] = temp;
    }
}
     
// Function to print matrix content at any stage
static void print(double [,]mat)
{
    for(int i = 0; i < N; i++, Console.WriteLine())
        for(int j = 0; j <= N; j++)
            Console.Write(mat[i, j]);
             
    Console.WriteLine();
}
 
// Function to reduce matrix to r.e.f.
static int forwardElim(double [,]mat)
{
    for(int k = 0; k < N; k++) 
    {
         
        // Initialize maximum value and index for pivot
        int i_max = k;
        int v_max = (int)mat[i_max, k];
         
        /* find greater amplitude for pivot if any */
        for(int i = k + 1; i < N; i++)
        {
            if (Math.Abs(mat[i, k]) > v_max) 
            {
                v_max = (int)mat[i, k];
                i_max = i;
            }
         
            /* If a principal diagonal element  is zero,
            *  it denotes that matrix is singular, and
            *  will lead to a division-by-zero later. */
            if (mat[k, i_max] == 0)
                return k; // Matrix is singular
             
            /* Swap the greatest value row with 
               current row
            */
            if (i_max != k)
                swap_row(mat, k, i_max);
             
            for(int a = k + 1; a < N; a++)
            {
             
                /* factor f to set current row kth element
                *  to 0, and subsequently remaining kth
                *  column to 0 */
                double f = mat[i, k] / mat[k, k];
                 
                /* subtract fth multiple of corresponding
                   kth row element*/
                for(int j = k + 1; j <= N; j++)
                    mat[i, j] -= mat[k, j] * f;
                 
                /* filling lower triangular matrix with
                *  zeros*/
                mat[i, k] = 0;
            }
        }
        // print(mat);        //for matrix state
    }
     
    // print(mat);            //for matrix state
    return -1;
}
 
// Function to calculate the values of the unknowns
static void backSub(double [,]mat)
{
     
    // An array to store solution
    double []x = new double[N]; 
     
    /* Start calculating from last equation up to the
       first */
    for(int i = N - 1; i >= 0; i--)
    {
     
        /* start with the RHS of the equation */
        x[i] = mat[i,N];
         
        /* Initialize j to i+1 since matrix is upper
        triangular*/
        for(int j = i + 1; j < N; j++) 
        {
         
            /* subtract all the lhs values
             * except the coefficient of the variable
             * whose value is being calculated */
            x[i] -= mat[i,j] * x[j];
        }
         
        /* divide the RHS by the coefficient of the
        unknown being calculated */
        x[i] = x[i] / mat[i,i];
    }
     
    Console.WriteLine();
    Console.WriteLine("Solution for the system:");
    for(int i = 0; i < N; i++)
    {
        Console.Write("{0:F6}", x[i]);
        Console.WriteLine();
    }
}
 
// Driver code
public static void Main(String[] args)
{
     
    /* input matrix */
    double [,]mat = { { 3.0, 2.0, -4.0, 3.0 },
                      { 2.0, 3.0, 3.0, 15.0 },
                      { 5.0, -3, 1.0, 14.0 } };
     
    gaussianElimination(mat);
}
}
 
// This code is contributed by shikhasingrajput

Javascript

// JavaScript program to demonstrate working of Gaussian Elimination
// method
 
let N = 3;
 
// function to get matrix content
function gaussianElimination(mat)
{
    /* reduction into r.e.f. */
    let singular_flag = forwardElim(mat);
 
    /* if matrix is singular */
    if (singular_flag != -1)
    {
        console.log("Singular Matrix.");
 
        /* if the RHS of equation corresponding to
           zero row  is 0, * system has infinitely
           many solutions, else inconsistent*/
        if (mat[singular_flag][N])
            console.log("Inconsistent System.");
        else
            console.log("May have infinitely many solutions.");
 
        return;
    }
 
    /* get solution to system and print it using
       backward substitution */
    backSub(mat);
}
 
// function for elementary operation of swapping two rows
function swap_row(mat, i, j)
{
    //printf("Swapped rows %d and %d\n", i, j);
 
    for (var k=0; k<=N; k++)
    {
        let temp = mat[i][k];
        mat[i][k] = mat[j][k];
        mat[j][k] = temp;
    }
}
 
// function to print matrix content at any stage
function print(mat)
{
    for (var i=0; i<N; i++, console.log(""))
        for (var j=0; j<=N; j++)
            process.stdout.write("" + mat[i][j]);
     
    console.log("");
}
 
// function to reduce matrix to r.e.f.
function forwardElim(mat)
{
    for (var k=0; k<N; k++)
    {
        // Initialize maximum value and index for pivot
        var i_max = k;
        var v_max = mat[i_max][k];
 
        /* find greater amplitude for pivot if any */
        for (var i = k+1; i < N; i++)
            if (Math.abs(mat[i][k]) > v_max)
                v_max = mat[i][k], i_max = i;
 
        /* if a principal diagonal element  is zero,
         * it denotes that matrix is singular, and
         * will lead to a division-by-zero later. */
        if (!mat[k][i_max])
            return k; // Matrix is singular
 
        /* Swap the greatest value row with current row */
        if (i_max != k)
            swap_row(mat, k, i_max);
 
 
        for (var i=k+1; i<N; i++)
        {
            /* factor f to set current row kth element to 0,
             * and subsequently remaining kth column to 0 */
            let f = mat[i][k]/mat[k][k];
 
            /* subtract fth multiple of corresponding kth
               row element*/
            for (var j=k+1; j<=N; j++)
                mat[i][j] -= mat[k][j]*f;
 
            /* filling lower triangular matrix with zeros*/
            mat[i][k] = 0;
        }
 
        //print(mat);        //for matrix state
    }
    //print(mat);            //for matrix state
    return -1;
}
 
// function to calculate the values of the unknowns
function backSub(mat)
{
    let x = new Array(N);  // An array to store solution
 
    /* Start calculating from last equation up to the
       first */
    for (var i = N-1; i >= 0; i--)
    {
        /* start with the RHS of the equation */
        x[i] = mat[i][N];
 
        /* Initialize j to i+1 since matrix is upper
           triangular*/
        for (var j=i+1; j<N; j++)
        {
            /* subtract all the lhs values
             * except the coefficient of the variable
             * whose value is being calculated */
            x[i] -= mat[i][j]*x[j];
        }
 
        /* divide the RHS by the coefficient of the
           unknown being calculated */
        x[i] = Math.round(x[i]/mat[i][i]);
    }
 
    console.log("\nSolution for the system:");
    for (var i=0; i<N; i++)
        console.log(x[i].toFixed(8));
}
 
// Driver program
 
/* input matrix */
let mat =   [[3.0, 2.0,-4.0, 3.0],
             [2.0, 3.0, 3.0, 15.0],
             [5.0, -3, 1.0, 14.0]];
 
gaussianElimination(mat);
 
 
// This code is contributed by phasing17

PHP

<?php
// PHP program to demonstrate working 
// of Gaussian Elimination method
 
$N = 3; // Number of unknowns
 
 
// function to get matrix content
function gaussianElimination($mat)
{
    global $N;
     
    /* reduction into r.e.f. */
    $singular_flag = forwardElim($mat);
 
    /* if matrix is singular */
    if ($singular_flag != -1)
    {
        print("Singular Matrix.\n");
 
        /* if the RHS of equation corresponding to
        zero row is 0, * system has infinitely
        many solutions, else inconsistent*/
        if ($mat[$singular_flag][$N])
            print("Inconsistent System.");
        else
            print("May have infinitely many solutions.");
 
        return;
    }
 
    /* get solution to system and print it using
    backward substitution */
    backSub($mat);
}
 
// function for elementary operation
// of swapping two rows
function swap_row(&$mat, $i, $j)
{
    global $N;
    //printf("Swapped rows %d and %d\n", i, j);
 
    for ($k = 0; $k <= $N; $k++)
    {
        $temp = $mat[$i][$k];
        $mat[$i][$k] = $mat[$j][$k];
        $mat[$j][$k] = $temp;
    }
}
 
// function to print matrix content at any stage
function print1($mat)
{
    global $N;
    for ($i=0; $i<$N; $i++, print("\n"))
        for ($j=0; $j<=$N; $j++)
            print($mat[$i][$j]);
 
    print("\n");
}
 
// function to reduce matrix to r.e.f.
function forwardElim(&$mat)
{
    global $N;
    for ($k=0; $k<$N; $k++)
    {
        // Initialize maximum value and index for pivot
        $i_max = $k;
        $v_max = $mat[$i_max][$k];
 
        /* find greater amplitude for pivot if any */
        for ($i = $k+1; $i < $N; $i++)
            if (abs($mat[$i][$k]) > $v_max)
                {
                    $v_max = $mat[$i][$k];
                    $i_max = $i;
                }
 
        /* if a principal diagonal element is zero,
        * it denotes that matrix is singular, and
        * will lead to a division-by-zero later. */
        if (!$mat[$k][$i_max])
            return $k; // Matrix is singular
 
        /* Swap the greatest value row with current row */
        if ($i_max != $k)
            swap_row($mat, $k, $i_max);
 
 
        for ($i = $k + 1; $i < $N; $i++)
        {
            /* factor f to set current row kth element to 0,
            * and subsequently remaining kth column to 0 */
            $f = $mat[$i][$k]/$mat[$k][$k];
 
            /* subtract fth multiple of corresponding kth
            row element*/
            for ($j = $k + 1; $j <= $N; $j++)
                $mat[$i][$j] -= $mat[$k][$j] * $f;
 
            /* filling lower triangular matrix with zeros*/
            $mat[$i][$k] = 0;
        }
 
        //print(mat); //for matrix state
    }
    //print(mat);     //for matrix state
    return -1;
}
 
// function to calculate the values of the unknowns
function backSub(&$mat)
{
    global $N;
    $x = array_fill(0, $N, 0); // An array to store solution
 
    /* Start calculating from last equation up to the
    first */
    for ($i = $N - 1; $i >= 0; $i--)
    {
        /* start with the RHS of the equation */
        $x[$i] = $mat[$i][$N];
 
        /* Initialize j to i+1 since matrix is upper
        triangular*/
        for ($j = $i + 1; $j < $N; $j++)
        {
            /* subtract all the lhs values
            * except the coefficient of the variable
            * whose value is being calculated */
            $x[$i] -= $mat[$i][$j] * $x[$j];
        }
 
        /* divide the RHS by the coefficient of the
        unknown being calculated */
        $x[$i] = $x[$i] / $mat[$i][$i];
    }
 
    print("\nSolution for the system:\n");
    for ($i = 0; $i < $N; $i++)
        print(number_format(strval($x[$i]), 6)."\n");
}
 
// Driver program
 
    /* input matrix */
    $mat = array(array(3.0, 2.0,-4.0, 3.0),
                        array(2.0, 3.0, 3.0, 15.0),
                        array(5.0, -3, 1.0, 14.0));
 
    gaussianElimination($mat);
 
// This code is contributed by mits
?>

Output

Solution for the system:
3.000000
1.000000
2.000000

Illustration:

process

Time Complexity: Since for each pivot we traverse the part to its right for each row below it, O(n)*(O(n)*O(n)) = O(n³).
We can also apply Gaussian Elimination for calculating:

Rank of a matrix
Determinant of a matrix
Inverse of an invertible square matrix

Approach:

Another approach for implementing Gaussian Elimination is using Partial Pivoting. In partial pivoting, at each iteration, the algorithm searches for the row with the largest pivot element (absolute value) and swaps that row with the current row. This ensures that the largest absolute value is used as the pivot element, reducing the rounding errors that can occur when dividing by a small pivot element. The steps involved in Partial Pivoting Gaussian Elimination are:

Start with the given augmented matrix.

Find the row with the largest absolute value in the first column and swap that row with the first row.

Divide the first row by the pivot element to make it equal to 1.

Use the first row to eliminate the first column in all other rows by subtracting the appropriate multiple of the first row from each row.

Repeat steps 2-4 for the remaining columns, using the row with the largest absolute value as the pivot row for each column.

Back-substitute to obtain the values of the unknowns.

implementation of above approach:

C++

#include <iostream>
#include <cmath>
using namespace std;
 
const int MAXN = 100;
 
void partial_pivot(double A[MAXN][MAXN+1], int n) {
    for (int i = 0; i < n; i++) {
        int pivot_row = i;
        for (int j = i+1; j < n; j++) {
            if (abs(A[j][i]) > abs(A[pivot_row][i])) {
                pivot_row = j;
            }
        }
        if (pivot_row != i) {
            for (int j = i; j <= n; j++) {
                swap(A[i][j], A[pivot_row][j]);
            }
        }
        for (int j = i+1; j < n; j++) {
            double factor = A[j][i] / A[i][i];
            for (int k = i; k <= n; k++) {
                A[j][k] -= factor * A[i][k];
            }
        }
    }
}
 
void back_substitute(double A[MAXN][MAXN+1], int n, double x[MAXN]) {
    for (int i = n-1; i >= 0; i--) {
        double sum = 0;
        for (int j = i+1; j < n; j++) {
            sum += A[i][j] * x[j];
        }
        x[i] = (A[i][n] - sum) / A[i][i];
    }
}
 
int main() {
    int n = 3;
    double A[MAXN][MAXN+1] = {{3.0, 2.0,-4.0, 3.0},
                              {2.0, 3.0, 3.0, 15.0},
                              {5.0, -3, 1.0, 14.0}
                             };
    double x[MAXN];
     
    partial_pivot(A, n);
    back_substitute(A, n, x);
     
    cout << "Solution for the system:\n";
    for (int i = 0; i < n; i++) {
       
        cout << x[i] << endl;
    }
     
    return 0;
}

Java

public class Main {
     
    private static final int MAXN = 100;
     
    public static void partialPivot(double[][] A, int n) {
        for (int i = 0; i < n; i++) {
            int pivotRow = i;
            for (int j = i+1; j < n; j++) {
                if (Math.abs(A[j][i]) > Math.abs(A[pivotRow][i])) {
                    pivotRow = j;
                }
            }
            if (pivotRow != i) {
                for (int j = i; j <= n; j++) {
                    double temp = A[i][j];
                    A[i][j] = A[pivotRow][j];
                    A[pivotRow][j] = temp;
                }
            }
            for (int j = i+1; j < n; j++) {
                double factor = A[j][i] / A[i][i];
                for (int k = i; k <= n; k++) {
                    A[j][k] -= factor * A[i][k];
                }
            }
        }
    }
     
    public static void backSubstitute(double[][] A, int n, double[] x) {
        for (int i = n-1; i >= 0; i--) {
            double sum = 0;
            for (int j = i+1; j < n; j++) {
                sum += A[i][j] * x[j];
            }
            x[i] = (A[i][n] - sum) / A[i][i];
        }
    }
     
    public static void main(String[] args) {
        int n = 3;
        double[][] A = {{3.0, 2.0,-4.0, 3.0},
                        {2.0, 3.0, 3.0, 15.0},
                        {5.0, -3, 1.0, 14.0}
                       };
        double[] x = new double[MAXN];
         
        partialPivot(A, n);
        backSubstitute(A, n, x);
         
        System.out.println("Solution for the system:");
        for (int i = 0; i < n; i++) {
            System.out.println(x[i]);
        }
    }
}

Python3

import numpy as np
MAXN = 100
 
# Function to perform partial pivot for Gaussian elimination
def partial_pivot(A, n):
    # Iterate through each row in the matrix
    for i in range(n):
        pivot_row = i
        # Find the row with the maximum absolute value in the current column
        for j in range(i + 1, n):
            if abs(A[j][i]) > abs(A[pivot_row][i]):
                pivot_row = j
        # Swap the current row with the row having the maximum absolute value
        if pivot_row != i:
            A[[i, pivot_row]] = A[[pivot_row, i]]
        # Perform Gaussian elimination on the matrix
        for j in range(i + 1, n):
            factor = A[j][i] / A[i][i]
            A[j] -= factor * A[i]
 
# Function to perform back substitution to solve the system of equations
def back_substitute(A, n):
    x = np.zeros(n)
    # Iterate through each row in reverse order
    for i in range(n - 1, -1, -1):
        sum_val = sum(A[i][j] * x[j] for j in range(i + 1, n))
        # Solve for x[i] using the previously calculated values of x
        x[i] = (A[i][n] - sum_val) / A[i][i]
    return x
# Driver code
if __name__ == "__main__":
    n = 3
    # Coefficient matrix augmented with the constant terms
    A = np.array([[3.0, 2.0, -4.0, 3.0],
                  [2.0, 3.0, 3.0, 15.0],
                  [5.0, -3, 1.0, 14.0]])
     
    # Perform Gaussian elimination with partial pivot
    partial_pivot(A, n)
    x = back_substitute(A, n)
     
    print("Solution for the system:")
    for i in range(n):
        print(x[i])
#This code is contributed by Vikram_Shirsat

C#

using System;
 
class Program {
  const int MAXN = 100;
 
  // Function to find the partial pivot
  static void partial_pivot(double[, ] A, int n)
  {
    for (int i = 0; i < n; i++) {
      int pivot_row = i;
      for (int j = i + 1; j < n; j++) {
        if (Math.Abs(A[j, i])
            > Math.Abs(A[pivot_row, i])) {
          pivot_row = j;
        }
      }
      if (pivot_row != i) {
        for (int j = i; j <= n; j++) {
          double temp = A[i, j];
          A[i, j] = A[pivot_row, j];
          A[pivot_row, j] = temp;
        }
      }
      for (int j = i + 1; j < n; j++) {
        double factor = A[j, i] / A[i, i];
        for (int k = i; k <= n; k++) {
          A[j, k] -= factor * A[i, k];
        }
      }
    }
  }
 
  // Function to perform the back
  // substitution
  static void back_substitute(double[, ] A, int n,
                              double[] x)
  {
    for (int i = n - 1; i >= 0; i--) {
      double sum = 0;
      for (int j = i + 1; j < n; j++) {
        sum += A[i, j] * x[j];
      }
      x[i] = (A[i, n] - sum) / A[i, i];
    }
  }
 
  // Driver Code
  static void Main(string[] args)
  {
    int n = 3;
    double[, ] A
      = new double[, ] { { 3.0, 2.0, -4.0, 3.0 },
                        { 2.0, 3.0, 3.0, 15.0 },
                        { 5.0, -3, 1.0, 14.0 } };
    double[] x = new double[n];
 
    partial_pivot(A, n);
    back_substitute(A, n, x);
 
    Console.WriteLine("Solution for the system:");
    for (int i = 0; i < n; i++) {
      Console.WriteLine(Math.Round(x[i]));
    }
  }
}

Javascript

// This program solves a system of linear equations using Gaussian Elimination with partial pivoting
const MAXN = 100;
 
// Function to perform partial pivoting on the matrix A
function partial_pivot(A, n) {
    for (let i = 0; i < n; i++) {
        let pivot_row = i;
        // Find the row with the largest absolute value in the ith column
        for (let j = i + 1; j < n; j++) {
            if (Math.abs(A[j][i]) > Math.abs(A[pivot_row][i])) {
                pivot_row = j;
            }
        }
        // Swap rows if necessary
        if (pivot_row != i) {
            for (let j = i; j <= n; j++) {
                [A[i][j], A[pivot_row][j]] = [A[pivot_row][j], A[i][j]];
            }
        }
        // Perform Gaussian Elimination on the matrix
        for (let j = i + 1; j < n; j++) {
            let factor = A[j][i] / A[i][i];
            for (let k = i; k <= n; k++) {
                A[j][k] -= factor * A[i][k];
            }
        }
    }
}
 
// Function to perform back substitution on the matrix A to find the solution vector x
function back_substitute(A, n, x) {
    for (let i = n - 1; i >= 0; i--) {
        let sum = 0;
        for (let j = i + 1; j < n; j++) {
            sum += A[i][j] * x[j];
        }
        x[i] = (A[i][n] - sum) / A[i][i];
    }
}
 
// Driver code
let n = 3;
let A = [
    [3.0, 2.0, -4.0, 3.0],
    [2.0, 3.0, 3.0, 15.0],
    [5.0, -3, 1.0, 14.0]
];
let x = Array(n);
 
partial_pivot(A, n);
back_substitute(A, n, x);
 
console.log("Solution for the system:");
for (let i = 0; i < n; i++) {
    console.log(Math.round(x[i]));
}

Output

Solution for the system:
3
1
2

complexity

The time complexity of Partial Pivoting Gaussian Elimination is O(N^3) for an N x N matrix.

This article is contributed by Yash Varyani. Please write comments if you find anything incorrect, or you want to share more information about the topic discussed above.

Recommended

Solve DSA problems on GfG Practice.

Solve Problems

Feeling lost in the world of random DSA topics, wasting time without progress? It’s time for a change! Join our DSA course, where we’ll guide you on an exciting journey to master DSA efficiently and on schedule.
Ready to dive in? Explore our Free Demo Content and join our DSA course, trusted by over 100,000 neveropen!

Gaussian Elimination to Solve Linear Equations

C++

Java

Python3

C#

Javascript

PHP

Approach:

implementation of above approach:

C++

Java

Python3

C#

Javascript

complexity

Smarter Retrieval for RAG: Late Chunking with Jina Embeddings v2 and Milvus

From Word2Vec to LLM2Vec: How to Choose the Right Embedding Model for RAG

How to Debug Slow Search Requests in Milvus

LEAVE A REPLY Cancel reply

Most Popular

Android 16 QPR2 Beta 3 lands with a flurry of bug fixes

Google is working on dedicated ‘Bills’ and ‘Travel’ folders for Gmail

Mint Mobile’s big bet on 5G home internet might change everything

Interviewed With Kyle Smith – Founder and CEO of Escalated by Shauli Zacks

EDITOR PICKS

Android 16 QPR2 Beta 3 lands with a flurry of bug fixes

Google is working on dedicated ‘Bills’ and ‘Travel’ folders for Gmail

Mint Mobile’s big bet on 5G home internet might change everything

POPULAR POSTS

Android 16 QPR2 Beta 3 lands with a flurry of bug fixes

Google is working on dedicated ‘Bills’ and ‘Travel’ folders for Gmail

Mint Mobile’s big bet on 5G home internet might change everything

POPULAR CATEGORY

ABOUT US

FOLLOW US