Наиболее рапространённым методом решения систем линейных уравнений является метод Гаусса с частичным выбором ведущего элемента по столбцам. Первый вариант этого метода здесь реализован в виде класса 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, который в процессе вычислений изменяется, и на месте правых столбцов будут находится столбцы с решениями.

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

Следущий вариант метода Гаусса реализован в виде шаблонов:

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

template <class Matrix>
bool _doGaussRow ( Matrix & data, nat nRow, nat nCol, nat index[], nat mRow, nat mCol )
{
    if ( nRow == 0 || nCol <= nRow ) return false;
    if ( mRow > nRow ) mRow = nRow;
    if ( mCol > nCol ) mCol = nCol;
    if ( mRow > mCol ) return false;
// Прямой ход
    nat i, j, k;
    for ( k = 0; k < mRow; ++k )
    {
// Поиск максимального по модулю члена в k-ой строке
        nat m = 0;
        double * rk = & data[k][0];
        double max = fabs ( rk[0] );
        for ( i = 1; i < mCol; ++i )
        {
            const double t = fabs ( rk[i] );
            if ( max < t ) max = t, m = i;
        }
        if ( max == 0 ) return false;
        index[k] = m;
// Нормализация строки
        const double p = 1. / rk[m];
        for ( i = 0; i < nCol; ++i ) rk[i] *= p;
        rk[m] = 1;
// Вычитание строк
        for ( j = k+1; j < nRow; ++j )
        {
            double * rj = & data[j][0];
            const double t = rj[m];
            for ( i = 0; i < nCol; ++i )
            {
                if ( fabs ( rj[i] -= rk[i] * t ) < 1e-290 ) rj[i] = 0;
            }
        }
    }
// Обратная подстановка
    for ( j = mRow; --j > 0; )
    {
        const double * rj = & data[j][0];
        for ( i = 0; i < j; ++i )
        {
            double * ri = & data[i][0];
            const double t = ri[index[j]];
            for ( k = 0; k < nCol; ++k ) ri[k] -= rj[k] * t;
        }
    }
    return true;
}

template <class T> class IMatrix;

template <class T> inline
bool doGaussRow ( IMatrix<T> & data, nat nRow, nat nCol, nat * index, nat mRow, nat mCol )
{
    return _doGaussRow ( data, nRow, nCol, index, mRow, mCol );
}

template <class T> class ArrRef2;

template <class T> inline
bool doGaussRow ( ArrRef2<T> & data, nat nRow, nat nCol, nat * index, nat mRow, nat mCol )
{
    return _doGaussRow ( data, nRow, nCol, index, mRow, mCol );
}

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