В этой статье мы будем рассматривать решения СЛАУ вида Ax = b, где A - симметричная разреженная матрица. Такие матрицы появляются, например, при решении задач методом наименьших квадратов. Для симметричных СЛАУ разработаны специальные методы, такие, как метод Холецкого и LDL^T разложение ( см. [1][2] ). Так как первый из них применим к более узкому классу матриц, чем второй, поэтому далее будем рассматривать только LDL^T разложение, хотя выводы этой статьи применимы к обоим методам.

При LDL^T разложении матрица A представляется в виде произведения LDL^T, где L - нижняя треугольная матрица с единицами на главной диагонали, D - диагональная матрица,
L^T- транспонированная матрица L. Затем последовательно решаются уравнения Ly = b и DL^Tx = y. При этом затрачивается около n³/6 умножений ( где n - размерность матрицы А ),
что в 2 раза меньше, чем в методе Гаусса. Это неудивительно, так как при большой размерности в симметричной матрице содержится приблизительно в 2 раза меньше информации, чем в матрице общего вида.
Далее, чтобы уменьшить количество вычислений, воспользуемся разреженностью матрицы А. Элементы матриц D и L вычисляются по формулам:

Из этих формул следует, что если в какой-то строке матрицы А первые k элементов были нулями, то и в матрице L в той же строке первые k элементов тоже будут нулями. А это значит, что эти элементы вычислять не надо, и, если их много, то получается большая экономия в вычислениях. Это основная идея этой статьи. Назовём эти нули - "лидирующими нулями".

Ниже приведена программа на языке С++, которая решает СЛАУ с учётом лидирующих нулей. Это реальная программа с проверками на переполнение и исчезновение порядка
для чисел с плавающей точкой. Матрица А поступает на вход программы в компактном виде. Во-первых, опускаются лидирующие нули, а во-вторых, указываются элементы только до диагонального, так как матрица симметричная. Новые матрицы L и D записываются на место исходной матрицы А, причём для матрицы D записываются обратные диагональные элементы.
Входные параметры функции slu_LDLt:
n - размерность матрицы А,
m - массив с количеством лидирующих нулей по строкам,
a - матрица А записанная в виде массива по строкам фрагментами от первого ненулевого до диагонального элемента,
b - массив свободных членов,
x - массив неизвестных.

typedef unsigned int nat;

template<class T> inline const T & _max ( const T & i1, const T & i2 )
{
    return i1 > i2 ? i1 : i2;
}

bool slu_LDLt ( nat n, const nat * m, double * a, const double * b, double * x )
{
    nat i, j, k, l, nn = 0;
// Заполнение индексов диагональных элементов
    std::vector<nat> di ( n );
    for ( i = 0; i < n; ++i )
    {
        if ( m[i] > i )
            return false;
        di[i] = nn += i - m[i];
        ++nn;
    }
// Построение матриц L и D
    for ( i = l = 0; i < n; ++i )
    {
        const nat mi = m[i];
        const nat ni = di[i] - i;
        const double * ai = a + ni;
        for ( j = mi; j <= i; ++j )
        {
            double * aj = a + di[j] - j;
            double x = ai[j];
            if ( i == j )
            {
                for ( k = mi; k < j; ++k )
                {
                    double & ajk = aj[k];
                    const double y = ajk;
                    ajk *= a[di[k]];
                    x -= y * ajk;
                }
                const double ax = fabs ( x );
                if ( ax < 1e-100 )
                {
                    return false;
                }
                a[l] = ax > 1e100 ? 0 : 1 / x;
            }
            else
            {
                for ( k = _max ( mi, m[j] ); k < j; ++k ) x -= ai[k] * aj[k];
                const double ax = fabs ( x );
                if ( ax > 1e100 )
                {
                    return false;
                }
                a[l] = ax < 1e-100 ? 0 : x;
            }
            ++l;
        }
    }
// Решение системы Ly = b
    for ( i = j = 0; i < n; ++i, ++j )
    {
        double t = b[i];
        for ( k = m[i]; k < i; ++k ) t -= a[j++] * x[k];
        x[i] = t;
    }
// Решение системы DLtx = y;
    for ( i = 0; i < n; ++i ) x[i] *= a[di[i]];
    for ( l = nn; --i > 0; )
    {
        const double t = x[i];
        const double * p = a + ( l -= i + 1 );
        for ( l += k = m[i]; k < i; ++k ) x[k] -= p[k] * t;
    }
    return true;
}

В случае отсутствия лидирующих нулей - это будет обычная программа для LDL^T разложения. А если их будет много, то время работы программы резко уменьшится по сравнению с обычной. В следующей статье будет показано, как можно увеличить количество лидирующих нулей, переставляя столбцы и строки исходной матрицы А.

Литература.

1. Уилкинсон, Райнш. "Справочник алгоритмов на языке АЛГОЛ. Линейная алгебра".

2. Голуб, Ван Лоун. "Матричные вычисления".