FAQ - часто задаваемые вопросы

#include <vector>
#include <algorithm>

//Назначение: решение СЛАУ вида Matrix * X = Fx
//Возвращаемое значение: true - если существует единственное решение,
//false - в остальных случаях
bool Gauss(
    const std::vector< std::vector<double> > &Matrix,
    const std::vector<double> &Fx,
    /*out*/ std::vector<double> &X)
{
    if (Fx.size() != Matrix.size()) return false;

    std::vector< std::vector<double> > matrix(Matrix);
    std::vector<double> fx(Fx);

    X.clear();
    X.resize(matrix.size());

    //прямой ход
    for (unsigned i = 0; i < matrix.size(); ++i)
    {
        for (unsigned Row = i + 1; Row < matrix.size(); ++Row)
        {
            //поменяем строки, если элемент на главной диагонали равен нулю
            if (0 == matrix)
            {
                unsigned rowNotZero;
                for (rowNotZero = i; rowNotZero < matrix.size(); ++rowNotZero)
                {
                    if (matrix[rowNotZero] != 0)    break;
                }

                if (rowNotZero == matrix.size())    return false;
                std::swap(matrix, matrix[rowNotZero]);
                std::swap(fx, fx[rowNotZero]);
            }

            double Coef = matrix[Row] / matrix;
            for (unsigned Col = i; Col < matrix[Row].size(); ++Col)
            {
                matrix[Row][Col] -= Coef * matrix[Col];
            }

            fx[Row] -= Coef * fx;
        }
    }

    //обратный ход
    for (int i = matrix.size() - 1; i >= 0; --i)
    {
        double sum = fx;
        for (unsigned j = i+1; j < matrix.size(); ++j)
        {
            sum -= matrix[j] * X[j];
        }

        X = sum / matrix;
    }

    return true;
}

#include <vector>
#include <algorithm>

//Назначение: решение СЛАУ вида Matrix * X = Fx
//Возвращаемое значение: true - если существует единственное решение,
//false - в остальных случаях
bool JordanGauss(
    const std::vector< std::vector<double> > &Matrix,
    const std::vector<double> &Fx,
    /*out*/ std::vector<double> &X)
{
    if (Fx.size() != Matrix.size()) return false;

    std::vector< std::vector<double> > matrix(Matrix);
    std::vector<double> fx(Fx);

    X.clear();
    X.resize(matrix.size());

    //приведем к треугольной матрице
    for (unsigned i = 0; i < matrix.size(); ++i)
    {
        for (unsigned Row = i + 1; Row < matrix.size(); ++Row)
        {
            //поменяем строки, если элемент на главной диагонали равен нулю
            if (0 == matrix)
            {
                unsigned rowNotZero;
                for (rowNotZero = i; rowNotZero < matrix.size(); ++rowNotZero)
                {
                    if (matrix[rowNotZero] != 0)    break;
                }

                if (rowNotZero == matrix.size())    return false;
                std::swap(matrix, matrix[rowNotZero]);
                std::swap(fx, fx[rowNotZero]);
            }

            double Coef = matrix[Row] / matrix;
            for (unsigned Col = i; Col < matrix[Row].size(); ++Col)
            {
                matrix[Row][Col] -= Coef * matrix[Col];
            }

            fx[Row] -= Coef * fx;
        }
    }

    //получим единицы на главной диагонали
    for (int Row = 0; Row < static_cast<int>(matrix.size()); ++Row)
    {
        fx[Row] /= matrix[Row][Row];
        for (int Col = matrix[Row].size() - 1; Col >= Row; --Col)
        {
            matrix[Row][Col] /= matrix[Row][Row];
        }
    }

    //приведем к единичной матрице слева
    for (int i = matrix.size() - 1; i > 0; --i)
    {
        for (int Row = i - 1; Row >= 0; --Row)
        {
            double Coef = matrix[Row] / matrix;
            for (int Col = matrix[Row].size() - 1; Col >= i; --Col)
            {
                matrix[Row][Col] -= Coef * matrix[Col];
            }

            fx[Row] -= Coef * fx;
        }
    }

    X = fx;

    return true;
}