LU.cs

using System;

namespace Algorithms.Numeric.Decomposition;

/// <summary>
///     LU-decomposition factors the "source" matrix as the product of lower triangular matrix
///     and upper triangular matrix.
/// </summary>
public static class Lu
{
    /// <summary>
    ///     Performs LU-decomposition on "source" matrix.
    ///     Lower and upper matrices have same shapes as source matrix.
    ///     Note: Decomposition can be applied only to square matrices.
    /// </summary>
    /// <param name="source">Square matrix to decompose.</param>
    /// <returns>Tuple of lower and upper matrix.</returns>
    /// <exception cref="ArgumentException">Source matrix is not square shaped.</exception>
    public static (double[,] L, double[,] U) Decompose(double[,] source)
    {
        if (source.GetLength(0) != source.GetLength(1))
        {
            throw new ArgumentException("Source matrix is not square shaped.");
        }

        var pivot = source.GetLength(0);
        var lower = new double[pivot, pivot];
        var upper = new double[pivot, pivot];

        for (var i = 0; i < pivot; i++)
        {
            for (var k = i; k < pivot; k++)
            {
                double sum = 0;

                for (var j = 0; j < i; j++)
                {
                    sum += lower[i, j] * upper[j, k];
                }

                upper[i, k] = source[i, k] - sum;
            }

            for (var k = i; k < pivot; k++)
            {
                if (i == k)
                {
                    lower[i, i] = 1;
                }
                else
                {
                    double sum = 0;

                    for (var j = 0; j < i; j++)
                    {
                        sum += lower[k, j] * upper[j, i];
                    }

                    lower[k, i] = (source[k, i] - sum) / upper[i, i];
                }
            }
        }

        return (L: lower, U: upper);
    }

    /// <summary>
    ///     Eliminates linear equations system represented as A*x=b, using LU-decomposition,
    ///     where A - matrix of equation coefficients, b - vector of absolute terms of equations.
    /// </summary>
    /// <param name="matrix">Matrix of equation coefficients.</param>
    /// <param name="coefficients">Vector of absolute terms of equations.</param>
    /// <returns>Vector-solution for linear equations system.</returns>
    /// <exception cref="ArgumentException">Matrix of equation coefficients is not square shaped.</exception>
    public static double[] Eliminate(double[,] matrix, double[] coefficients)
    {
        if (matrix.GetLength(0) != matrix.GetLength(1))
        {
            throw new ArgumentException("Matrix of equation coefficients is not square shaped.");
        }

        var pivot = matrix.GetLength(0);
        var upperTransform = new double[pivot, 1]; // U * upperTransform = coefficients
        var solution = new double[pivot]; // L * solution = upperTransform
        (double[,] l, double[,] u) = Decompose(matrix);

        for (var i = 0; i < pivot; i++)
        {
            double pivotPointSum = 0;

            for (var j = 0; j < i; j++)
            {
                pivotPointSum += upperTransform[j, 0] * l[i, j];
            }

            upperTransform[i, 0] = (coefficients[i] - pivotPointSum) / l[i, i];
        }

        for (var i = pivot - 1; i >= 0; i--)
        {
            double pivotPointSum = 0;

            for (var j = i; j < pivot; j++)
            {
                pivotPointSum += solution[j] * u[i, j];
            }

            solution[i] = (upperTransform[i, 0] - pivotPointSum) / u[i, i];
        }

        return solution;
    }
}