В предыдущей статье было показано, что при решении СЛАУ с симметричной разреженной матрицей наличие лидирующих нулей приводит к уменьшению количества вычислений. В этой статье будет представлен алгоритм, предназначенный для увеличения количества лидирующих нулей данной матрицы. Если переставить i-ую и j-ую строки, а также i-ый и j-ый столбцы, то матрица останется симметричной. Такие перестановки называют симметричными. Они могут менять количество лидирующих нулей и, если их правильно применять, то количество лидирующих нулей можно увеличить. Другими словами, нам надо сделать так, чтобы все ненулевые члены по возможности находились возле главной диагонали. В частности, если известно, что матрица - ленточная, то делать ничего не надо.
Предлагается следующий алгоритм.

Вначале выбираем столбец ( или строку, что неважно, так как матрица симметричная ) с минимальным числом ненулевых элементов. Если таких столбцов несколько, то выбирается какой-то из них. При помощи симметричной перестановки делаем этот столбец первым.
Таким образом количество нулей в этом столбце будет максимальным. Далее строки в которых были не нули игнорируем. Находим столбец с минимальным числом ненулевых элементов без учёта этих строк и делаем его следующим. И так далее пока не пройдём всю матрицу.

Ниже помимо краткого текстового описания программы приводится много кода на С++, который сам по себе является точным описанием алгоритма.

Вначале введём несколько вспомогательных классов и функций.

Шаблон структур SortItem предназначен для сравнения объектов по выделенной части. В поле head записывается та часть, которая сравнивается, а в поле tail та, которая не сравнивается.

template <class H, class T=H> struct SortItem
{
    H head;
    T tail;

    SortItem () {}
    SortItem ( const H & h ) : head(h) {} 
    SortItem ( const H & h, const T & t ) : head(h), tail(t) {} 
};

template <class H, class T> inline
bool operator < ( const SortItem<H,T> & a, const SortItem<H,T> & b )
{
    return a.head < b.head;
}

template <class H, class T> inline
bool operator > ( const SortItem<H,T> & a, const SortItem<H,T> & b )
{
    return a.head > b.head;
}

template <class H, class T> inline
bool operator <= ( const SortItem<H,T> & a, const SortItem<H,T> & b )
{
    return a.head <= b.head;
}

template <class H, class T> inline
bool operator >= ( const SortItem<H,T> & a, const SortItem<H,T> & b )
{
    return a.head >= b.head;
}

template <class H, class T> inline
bool operator == ( const SortItem<H,T> & a, const SortItem<H,T> & b )
{
    return a.head == b.head;
}

template <class H, class T> inline
bool operator != ( const SortItem<H,T> & a, const SortItem<H,T> & b )
{
    return a.head != b.head;
}

Шаблон _swap делает обмен данными у двух объектов:

template<class T> inline void _swap ( T & i1, T & i2 )
{
    const T i ( i1 );
    i1 = i2;
    i2 = i;
}

Ещё нам понадобится следующая специализация шаблона _swap:

typedef unsigned int nat; // натуральное число

inline void _swap ( SortItem<nat, std::pair<nat, nat> *> & p1, 
                    SortItem<nat, std::pair<nat, nat> *> & p2 )
{
    const SortItem<nat, std::pair<nat, nat> *> p ( p1 );
    p1 = p2;
    p2 = p;
    _swap ( p1.tail->second, p2.tail->second );
}

Шаблон классов MinHeap представляет собой очередь с минимальным приоритетом, реализованную посредством конструкции "пирамида" ( binary heap ). В ней можно элементы добавлять, удалять, менять, переставлять.

template <class T>
class BaseHeap // Очередь с приоритетом
{
protected:
    nat _size, max_size;
    T * heap;
// Запрет конструктора копии и оператора присваивания
    BaseHeap ( const BaseHeap & );
    void operator = ( const BaseHeap & );
public:
    explicit BaseHeap ( nat n ) : _size(0), max_size(n), heap ( new T[n+1] ) {}
    ~BaseHeap () { delete[] heap; }
// Указатель на i-тый элемент ( 0, если i выходит за предел )
    const T * operator[] ( nat i ) const { return _size > i ? heap + i + 1 : 0; }
          T * operator[] ( nat i )       { return _size > i ? heap + i + 1 : 0; }
// Количество элементов в очереди
    nat size() const { return _size; }
// Очистить очередь
    void clear() { _size = 0; }
};

template <class T> // Очередь с минимальным приоритетом
class MinHeap : public BaseHeap<T>
{
// Запрет конструктора копии и оператора присваивания
    MinHeap ( const MinHeap & );
    void operator = ( const MinHeap & );
public:
    explicit MinHeap ( nat n ) : BaseHeap<T> ( n ) {}
// Поднять i-тый элемент
    void raise ( nat i )
    {
        if ( i >= _size ) return;
        for ( ++i; i > 1; )
        {
            const nat j = i >> 1;
            if ( heap[j] <= heap[i] ) break;
            _swap ( heap[i], heap[j] );
            i = j;
        }
    }
// Опустить i-тый элемент
    void down ( nat i )
    {
        for (++i;;)
        {
            const nat i1 = i + i;
            if ( i1 > _size ) break;
            const nat i2 = i1 + 1;
            const nat j = i2 > _size ? i1 : heap[i1] <= heap[i2] ? i1 : i2;
            if ( heap[i] <= heap[j] ) break;
            _swap ( heap[i], heap[j] );
            i = j;
        }
    }
// Удалить i-тый элемент
    bool del ( nat i )
    {
        if ( i >= _size ) return false;
        const nat i1 = i + 1;
        if ( i1 == _size ) { _size--; return true; }
        _swap ( heap[i1], heap[_size] );
        if ( heap[i1] > heap[_size--] )
            down ( i );
        else
            raise ( i );
        return true;
    }
// Добавить элемент в очередь
    bool operator << ( const T & o )
    {
        if ( _size == max_size ) return false;
        heap[++_size] = o;
        raise ( _size - 1 );
        return true;
    }
// Убрать минимальный элемент из очереди
    bool operator >> ( T & o )
    {
        if ( _size == 0 ) return false;
        o = heap[1];
        _swap ( heap[1], heap[_size--] );
        down ( 0 );
        return true;
    }
};

Сравнение пар по первому элементу:

inline bool compare ( const std::pair<nat, double> & p1, 
                      const std::pair<nat, double> & p2 )
{
    return p1.first <= p2.first;
}

На вход функции slu_LDLtO матрица поступает в виде массива строк data. Каждая строка - это массив пар. Первое число пары - номер столбца, второе - значение элемента матрицы. Далее на входе следуют массивы свободных членов и неизвестных.

bool slu_LDLtO ( nat n, const std::vector<std::pair<nat, double>> * data, 
                 const double * b, double * x )
{
    nat i, j, k, l;
    std::vector<nat> index(n), index2(n), fi(n), m(n);
    // Переставим столбцы и строки для увеличения лидирующих нулей
    {
        MinHeap< SortItem<nat, std::pair<nat, nat> *> > heap ( n );
        std::vector<std::pair<nat, nat> > elem ( n );
        for ( i = 0; i < n; ++i )
        {
            std::pair<nat, nat> & e = elem[i];
            e.first = i;
            e.second = heap.size();
            heap << SortItem<nat, std::pair<nat, nat> *> ( data[i].size(), &elem[i] );
        }
        std::vector<bool> flag ( n, 0 );
        SortItem<nat, std::pair<nat, nat> *> t;
        for ( l = 0; l < n; ++l )
        {
            heap >> t;
            t.tail->second = n;
            index[l] = t.tail->first;
            const std::vector<std::pair<nat, double>> & a = data[t.tail->first];
            for ( i = 0; i < a.size(); ++i ) flag[a[i].first] = true;
            for ( i = 0; i < a.size(); ++i )
            {
                const std::pair<nat, nat> & e = elem[a[i].first];
                if ( e.second == n ) continue;
                const std::vector<std::pair<nat, double>> & b = data[e.first];
                nat sum = 0;
                for ( j = 0; j < b.size(); ++j ) if ( ! flag[b[j].first] ) ++sum;
                heap[e.second]->head = sum;
                heap.raise ( e.second );
            }
        }
    }
// Формирование параметров и вызов функции slu_LDLt
    for ( i = 0; i < n; ++i ) index2[index[i]] = i;
    for ( i = 0; i < n; ++i )
    {
        const std::vector<std::pair<nat, double> > & arr = data[i];
        nat min = index2[arr[0].first];
        for ( nat j = 1; j < arr.size(); ++j ) if ( min > index2[arr[j].first] ) min = index2[arr[j].first];
        fi[i] = min;
    }
    nat nbuf = n * ( n + 1 ) / 2;
    nat max = 0;
    for ( i = 0; i < n; ++i )
    {
        nbuf -= m[i] = fi[index[i]];
        if ( max < data[i].size() ) max = data[i].size();
    }
    std::vector<double> b2(n), buf(nbuf);
    for ( i = 0; i < n; ++i ) b2[index2[i]] = b[i];
    std::vector<std::pair<nat, double>> arr2 ( max );
    for ( l = i = 0; i < n; ++i )
    {
        const std::vector<std::pair<nat, double> > & arr = data[index[i]];
        const nat na = arr.size();
        arr2.resize ( na );
        for ( j = 0; j < na; ++j )
        {
            arr2[j].first = index2[arr[j].first];
            arr2[j].second = arr[j].second;
        }
        std::sort ( begin(arr2), end(arr2), compare );
        for ( j = m[i], k = 0; j <= i; ++j )
        {
            buf[l++] = k < na && j == arr2[k].first ? arr2[k++].second: 0;
        }
    }
    std::vector<double> y(n);
    if ( ! slu_LDLt ( n, &m[0], &buf[0], &b2[0], &y[0]) ) return false;
    for ( i = 0; i < n; ++i ) x[i] = y[index2[i]];
    return true;
}

Вначале заполняется очередь heap и создаётся массив flag для пометки строк у которых появились ненулевые элементы. Затем в цикле из очереди достаётся столбец с минимальным числом ненулевых элементов. В массиве flag отмечаются те строки, которые в этом столбце содержат ненулевые элементы. Также пересчитываются те столбцы в очереди у которых отмечены строки с ненулевыми элементами. Когда мы пройдём все столбцы, то получим перестановки записанные в массив index. Также создадим массив обратных перестановок index2. Для того, чтобы решить систему уравнений, нужно подготовить параметры для функции slu_LDLt. Этому посвящена оставшаяся часть программы.
В зависимости от структуры входной матрицы функция slu_LDLtO может быть быстрее, чем slu_LDLt, а может быть медленнее. Поэтому лучше вначале испытать данный тип матрицы, а потом по результату выбрать функцию.