Gauss Jordan Elimination Algorithm
The Gauss-Jordan Elimination Algorithm is a fundamental linear algebra technique used to solve systems of linear equations and find the inverse of a matrix. It is an extension of Gaussian elimination, which is based on the principle of row operations to create an upper triangular matrix. The Gauss-Jordan method goes a step further by transforming the given matrix into a Reduced Row Echelon Form (RREF), which has a diagonal of ones and zeros below and above the diagonal. This simplified form allows for a more straightforward determination of the solutions to the system of equations or the inverse of a matrix.
The algorithm begins by performing a sequence of three row operations: row swapping, row scaling, and row addition. These operations are applied iteratively to the augmented matrix, which is formed by appending the constant terms of each equation to the original matrix. The objective is to obtain a matrix with a diagonal of ones and zeros elsewhere in each column, except for the last column of the augmented matrix. Row swapping reorders the rows to avoid a zero pivot element, while row scaling ensures the pivot element is one, and row addition eliminates other non-zero elements in the pivot column. Upon reaching RREF, the solutions to the system of linear equations can be easily read off, or if the matrix is invertible, its inverse can be obtained from the augmented matrix. The Gauss-Jordan elimination algorithm is widely used in various fields, including engineering, computer science, and economics, for solving problems requiring linear systems and matrix inversions.
using System;
namespace Algorithms.Numeric
{
/// <summary>
/// TODO.
/// </summary>
public class GaussJordanElimination
{
private int RowCount { get; set; }
/// <summary>
/// Method to find a linear equation system using gaussian elimination.
/// </summary>
/// <param name="matrix">The key matrix to solve via algorithm.</param>
/// <returns>whether the input matrix has a unique solution or not.
/// and solves on the given matrix. </returns>
public bool Solve(double[,] matrix)
{
RowCount = matrix.GetUpperBound(0) + 1;
if (!CanMatrixBeUsed(matrix))
{
throw new ArgumentException("Please use a n*(n+1) matrix with Length > 0.");
}
var pivot = PivotMatrix(ref matrix);
if (!pivot)
{
return false;
}
Elimination(ref matrix);
return ElementaryReduction(ref matrix);
}
/// <summary>
/// To make simple validation of the matrix to be used.
/// </summary>
/// <param name="matrix">Multidimensional array matrix.</param>
/// <returns>True: if algorithm can be use for given matrix;
/// False: Otherwise. </returns>
private bool CanMatrixBeUsed(double[,] matrix) => matrix?.Length == RowCount * (RowCount + 1) && RowCount > 1;
/// <summary>
/// To prepare given matrix by pivoting rows.
/// </summary>
/// <param name="matrix">Input matrix.</param>
/// <returns>Matrix.</returns>
private bool PivotMatrix(ref double[,] matrix)
{
for (var col = 0; col + 1 < RowCount; col++)
{
if (matrix[col, col] == 0)
{
// To find a non-zero coefficient
var rowToSwap = FindNonZeroCoefficient(ref matrix, col);
if (matrix[rowToSwap, col] != 0)
{
var tmp = new double[RowCount + 1];
for (var i = 0; i < RowCount + 1; i++)
{
// To make the swap with the element above.
tmp[i] = matrix[rowToSwap, i];
matrix[rowToSwap, i] = matrix[col, i];
matrix[col, i] = tmp[i];
}
}
else
{
// To return that the matrix doesn't have a unique solution.
return false;
}
}
}
return true;
}
private int FindNonZeroCoefficient(ref double[,] matrix, int col)
{
var rowToSwap = col + 1;
// To find a non-zero coefficient
for (; rowToSwap < RowCount; rowToSwap++)
{
if (matrix[rowToSwap, col] != 0)
{
return rowToSwap;
}
}
return col + 1;
}
/// <summary>
/// Applies REF.
/// </summary>
/// <param name="matrix">Input matrix.</param>
private void Elimination(ref double[,] matrix)
{
for (var srcRow = 0; srcRow + 1 < RowCount; srcRow++)
{
for (var destRow = srcRow + 1; destRow < RowCount; destRow++)
{
var df = matrix[srcRow, srcRow];
var sf = matrix[destRow, srcRow];
for (var i = 0; i < RowCount + 1; i++)
{
matrix[destRow, i] = matrix[destRow, i] * df - matrix[srcRow, i] * sf;
}
}
}
}
/// <summary>
/// To continue reducing the matrix using RREF.
/// </summary>
/// <param name="matrix">Input matrix.</param>
/// <returns>True if it has a unique solution; false otherwise.</returns>
private bool ElementaryReduction(ref double[,] matrix)
{
for (var row = RowCount - 1; row >= 0; row--)
{
var element = matrix[row, row];
if (element == 0)
{
return false;
}
for (var i = 0; i < RowCount + 1; i++)
{
matrix[row, i] /= element;
}
for (var destRow = 0; destRow < row; destRow++)
{
matrix[destRow, RowCount] -= matrix[destRow, row] * matrix[row, RowCount];
matrix[destRow, row] = 0;
}
}
return true;
}
}
}
