gistfile1.txt

/// <summary>
/// Решает СЛАУ методом Гаусса.
/// Возвращает true если решение найдено, false если однозначного решения нет
/// (такое может произойти, если существуют строки, целиком содержащие нули или такие строки
/// появляются в результате элементарных преобразований).
/// Если метод вернул true, то в последнем столбце матрицы будут записаны найденные значения коэффициентов,
/// а матрица будет приведена к канонической форме.
/// </summary>
/// <param name="m">Расширенная матрица M+1 x M - последний столбец содержит вектор свободных членов</param>
/// <param name="epsilon">Значение для оценки близких к нулю чисел (0.00001, к примеру)</param>
/// <param name="debug">Будут выведены диагностические сообщения, если debug задан в true</param>
public static bool SolveUsingGauss(double[,] m, double epsilon = DEFAULT_EPSILON, bool debug = false) {
    if (m == null)
        throw new ArgumentNullException("m");
    // в СЛАУ должно быть как минимум 2 переменных
    if (m.Length < 6)
        throw new ArgumentException("Source matrix is too small.");
    int X = m.GetLength(0); // размер матрицы по X (включая столбец свободных членов)
    int Y = m.GetLength(1); // размер матрицы по Y
    if (X != Y + 1)
        throw new ArgumentException("Only square matrixes are supported.");

    if (debug) {
        Console.WriteLine("Source martix :");
        DumpMatrix(m);
        Console.WriteLine("Elimination :");
    }

    // для каждой строки выполняем прямой ход метода Гаусса
    // условие y < Y казалось бы можно оптимизировать до y < Y - 1, однако в этом случае
    // последняя строка не будет проверена на наличие ненулевого первого коэффициента
    for (int y = 0; y < Y; y++) {
        // если текущий элемент главной диагонали m[y, y] = 0, то ищем внизу матрицы строку
        // с ненулевым элементом в этой позиции и меняем их местами
        // если такой строки нет, то это означает, что однозначного решения система не имеет,
        // и функция возвращает false
        if (Math.Abs(m[y, y] - .0d) < epsilon) {
            bool swap_found = false;

            for (int row = y + 1; row < Y; row++) {
                if (Math.Abs(m[y, row] - .0d) > epsilon) {
                    for (int x = 0; x < X; x++) {
                        double tmp = m[x, row];
                        m[x, row] = m[x, y];
                        m[x, y] = tmp;
                    }
                    swap_found = true;
                    break;
                }
            }

            if (!swap_found) {
                if (debug) {
                    Console.WriteLine("There are no unique solution, function will return FALSE");
                }
                return false;
            }
        }

        // работаем со строками, находящимися ниже текущей - приводим их первые коэффициенты к нулю
        for (int k = y + 1; k < Y; k++) {
            double df = m[y, y];
            double sf = m[y, k];
            
            if (debug) {
                Console.WriteLine("sf/df = {0, 14}/{1, 14}", sf, df);
            }

            m[y, k] = .0d;
            for (int i = y + 1; i < X; i++) {
                m[i, k] = (m[i, k]*df - m[i, y]*sf)/df;
            }

            if (debug) {
                DumpMatrix(m);
            }
        }
    }

    if (debug) {
        Console.WriteLine("Back insert :");
    }

    // выполняем обратный ход Гаусса
    for (int y = Y - 1; y >= 0; y--) {
        // делим текущую строку на первый коэффициент (чтобы сделать его единицей)
        double d = m[y, y];
        for (int x = y; x < X; x++)
            m[x, y] /= d;
        // вычисляем значение текущей переменной и записываем его в столбец свободных членов
        for (int x = y + 1; x < X - 1; x++) {
            // вычитаем из самого правого столбца все остальные (кроме первого)
            m[X - 1, y] -= m[x, y];
            // и обнуляем использованные слагаемые чтобы привести матрицу к стандартному виду
            m[x, y] = .0d;
        }
        // применяем вычисленное значение переменной к строкам выше текущей
        double variable = m[X - 1, y];
        if (Math.Abs(variable - 0) < epsilon) {
            if (debug) {
                Console.WriteLine("There are no unique solution, function will return FALSE");
            }
            return false;
        }
        for (int i = 0; i < y; i++) {
            m[y, i] *= variable;
        }

        if (debug) {
            DumpMatrix(m);
        }
    }

    return true;
}