Если в системе линейных уравнений матрица является симметричной и к тому же положительно определённой ( xAx > 0, если x != 0 ), то в этом случае можно использовать для решения специальные методы.

Одним из таких методов является разложение Холецкого. Он представляет исходную матрицу в виде произведения: A = LL^T, где L - нижняя треугольная матрица. Для того, чтобы найти решение системы Ax = b, последовательно решаются системы Ly = b, L^Tx = y. Здесь разложение Холецкого реализовано в виде класса SM_Chol. Конструктор осуществляет разложение исходной матрицы и хранит результат в одномерном массиве. Для диагональных элементов матрицы L записаны обратные значения. Метод determinant вычисляет определитель матрицы, а метод solve решает систему уравнений. Заметим, что в качестве аргументов для метода solve могут быть указатели на один и тот же массив. В этом случае исходные данные заменятся на вектор решения. Для того, чтобы найти обратную матрицу нужно n раз применить метод solve к соответствующему столбцу единичной матрицы, при этом будет получаться столбец обратной матрицы.

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

Другим методом является LDL^T разложение, где L - нижняя треугольная матрица с единичной диагональю, а D - диагональная матрица. В отличии от разложения Холецкого этот метод не использует операцию извлечения квадратного корня. Класс SM_LDLt имеет интерфейс аналогичный интерфейсу класса SM_Chol.

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

Оба класса были написаны на основе соответствующих программ из книги "Справочник алгоритмов на языке Алгол" ( Уилкинсон, Райнш ).

Если матрица сильно разреженна, то для неё применяют специальный вариант метода LDL^T:

bool slu_LDLt ( nat n, const nat * m, double * a, const double * b, double * x );
bool slu_LDLt ( nat n, const Suite<SortItem<nat, double> > * data, const double * b, double * x );

Функция slu_LDLtO ( O - оптимизирующая ) вначале переставляет строки и столбцы матрицы для того, чтобы увеличить количество лидирующих нулей и тем самым ускоряет вычисления:

bool slu_LDLtO ( nat n, const Suite<SortItem<nat, double> > * data, const double * b, double * x );

Также для разреженных матриц применяют метод сопряжённых градиентов:

void slu_cg ( nat n, const Suite<SortItem<nat, double> > * data, const double * b, double * x );

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