Решение больших линейных систем с блочными разреженными матрицами

Я хочу решить Ax = b, где A - очень большая квадратная положительно определенная симметричная блочная матрица, а x и b - векторы. Когда я говорю "большой", я имею в виду для матрицы nxn a n размером 300,000.

Вот пример гораздо меньшей, но репрезентативной матрицы, которую я хочу решить.

И вот та же матрица, увеличенная, показывая, что она состоит из блоков плотных матриц.

Я ранее (см. здесь, здесь, а еще здесь) использовал Eigen Cholesky solver, который отлично работал для n<10000, но с помощью n=300000 холеский решатель был слишком медленным. Однако я не воспользовался тем, что у меня есть блок-матрица. По-видимому, существуют алгоритмы для решения разреженных блочных матриц (например, блокировать холескую факторизацию).

Я хотел бы знать конкретно, если Eigen оптимизировал алгоритмы, используя факторизацию или итеративные методы, для разреженных плотных матриц блоков, которые я могу использовать?

Также вы можете предложить другие алгоритмы, которые могут быть идеальными для решения моей матрицы? Я имею в виду, насколько я понимаю, по крайней мере, для факторизации найти матрицу перестановок NP трудно, так много различных эвристических методов существуют, и насколько я могу сказать, люди создают интуицию различных матричных структур (например, banded matrix) и какие алгоритмы работают лучше всего с ними. У меня еще нет этой интуиции.

В настоящее время я изучаю метод метода сопряжения градиента. Я сам реализовал это на С++, но все еще не воспользовался тем, что матрица является блочной матрицей.

//solve A*rk = xk
//Eigen::SparseMatrix<double> A;
//Eigen::VectorXd rk(n);
Eigen::VectorXd pk = rk;

double rsold = rk.dot(rk);
int maxiter = rk.size();
for (int k = 0; k < maxiter; k++) {
    Eigen::VectorXd Ap = A*pk;  
    double ak = rsold /pk.dot(Ap);
    xk += ak*pk, rk += -ak*Ap;      
    double rsnew = rk.dot(rk);
    double xlength = sqrt(xk.dot(xk));
    //relaxing tolerance when x is large
    if (sqrt(rsnew) < std::max(std::min(tolerance * 10000, tolerance * 10 * xlength), tolerance)) {
        rsold = rsnew;
        break;
    }
    double bk = rsnew / rsold;
    pk *= bk;
    pk += rk;
    rsold = rsnew;
}