Наиболее распространённым методом решения систем линейных уравнений является метод Гаусса с частичным выбором ведущего элемента по столбцам. Первый вариант этого метода здесь реализован в виде класса SLU_Gauss, написанного на основе двух фортран-программ DECOMP и SOLVE из книги "Машинные методы математических вычислений" (Дж.Форсайт, М.Малькольм, К.Моулер). Этот класс позволяет вычислить определитель матрицы, оценить её число обусловленности, решить систему уравнений. Для того, чтобы найти обратную матрицу нужно n раз применить метод solve к соответствующему столбцу единичной матрицы, при этом будет получаться столбец обратной матрицы.

typedef nat int nat;

class SLU_Gauss
{
    const nat n;
    DynArray<nat> ipvt;
    DynArray2<double> a;
    double cond;
// Запрет конструктора копии и оператора присваивания:
    SLU_Gauss ( SLU_Gauss & );
    void operator = ( SLU_Gauss & );
public:
    SLU_Gauss ( nat k, const double * const * a );
    bool solve ( const double * b, double * x ) const; // b[n], x[n]
    bool solve ( const double * const * a, const double * b, double * x ) const;
    double condition () const { return cond; }
    double determinant () const;
};

Конструктор класса приводит исходную матрицу к треугольному виду и приблизительно вычисляет число обусловленности, которое здесь определено, как произведение 1-нормы матрицы A ( вычисляется явно ) на 1-норму матрицы A^-1 ( оценивается приближённо ). Если число обусловленности очень большое, т.е. cond + 1 равен cond, то решать такую систему смысла нет и в этом случае методы solve возвращают значение false. Заметим, что в качестве аргументов для методов solve могут быть указатели на один и тот же массив ( b и x ). Метод solve с двумя параметрами просто решает систему уравнений, а метод solve с тремя параметрами ( первый из которых - это указатель на исходную матрицу ) вначале вызывает solve с двумя параметрами, а затем уточняет результат уменьшая квадратичную норму невязки. В моих экспериментах заметный эффект от уточнения наблюдался в редких случаях.

Второй вариант реализации метода Гаусса в отличии от первого не использует операторы new-delete:

bool slu_gauss ( ArrRef2<double> & data ); // несколько правых столбцов
bool slu_gauss ( ArrRef2<double> & data, double * x ); // один правый столбец

Здесь квадратная матрица и столбцы свободных членов записаны в один двухмерный массив data, который в процессе вычислений изменяется, и на месте правых столбцов будут находится столбцы с решениями.

Следущий вариант метода Гаусса выбирает ведущие элементы по всей указанной подматрице:

//*************************** 24.03.2011 ******************************//
//
//      Решение систем линейных уравнений методом Гаусса
//      Общий выбор ведущего элемента
//      nc - к-во столбцов в левой части системы уравнений
//      col - массив индексов столбцов
//
//*************************** 24.03.2011 ******************************//

bool slu_gauss ( ArrRef2<double> & data, const nat nc, ArrRef<nat> & col );

Исходники находятся в файле mathem.cpp.

Следущий вариант метода Гаусса выбирает ведущие элементы по строкам:

//*************************** 19.02.2017 ******************************//
//
//      Метод исключений Гаусса.
//      Выбор ведущего элемента по строкам.
//      nRow, nCol - размеры матрицы.
//      index - массив для индексов выбранных столбцов (nCol).
//      mCol - размер подматрицы, где выбираются ведущие элементы.
//
//*************************** 02.06.2019 ******************************//

bool sluGaussRow ( IMatrix<T> & data, nat nRow, nat nCol, unsigned * index, nat mCol );
bool sluGaussRow ( ArrRef2<T> & data, nat nRow, nat nCol, unsigned * index, nat mCol );

Исходники находятся в файле LinAlg.h.