GaussJordanElimination.cs

namespace Algorithms.Numeric;

/// <summary>
///     Algorithm used to find the inverse of any matrix that can be inverted.
/// </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;
    }
}