c++ - siguientes - sistema de ecuaciones 2x2 por matriz inversa

Resolución de grandes sistemas lineales con matrices dispersas de bloques.

Question

Quiero resolver Ax = b donde A es una matriz de bloques simétrica definida positiva muy grande cuadrada y x y b son vectores. Cuando digo grande me refiero a una matriz nxn una n tan grande como 300,000 .

Aquí hay un ejemplo de una matriz mucho más pequeña pero representativa que quiero resolver.

Y aquí está la misma matriz ampliada para mostrar que está compuesta de bloques de matrices densas.

Anteriormente (vea here , here , y tan atrás como here ) usé el solucionador Cholesky de Eigen que funcionó bien para n<10000 pero con n=300000 el solucionador Cholesky es demasiado lento. Sin embargo, no aproveché el hecho de que tengo una matriz de bloques. Aparentemente, existen algoritmos para resolver matrices de bloques dispersos (por ejemplo , factorización de bloques y bloques ).

Me gustaría saber específicamente si Eigen ha optimizado los algoritmos, utilizando métodos de factorización o iterativos, para matrices de bloques densos dispersos que puedo emplear?

¿También puede sugerir otros algoritmos que podrían ser ideales para resolver mi matriz? Me refiero a que, por lo que entiendo, al menos para la factorización, encontrar una matriz de permutación es difícil en NP, por lo que existen muchos métodos heurísticos diferentes y, por lo que puedo decir, las personas desarrollan una intuición de diferentes estructuras matriciales (por ejemplo, matriz en bandas ) y qué algoritmos funcionan mejor. con ellos. Todavía no tengo esta intuición.

Actualmente estoy estudiando el uso de un método de gradiente conjugado . Lo he implementado yo mismo en C ++ pero todavía no aprovecho el hecho de que la matriz es una matriz de bloques.

//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;
}

Answer 1

Creo que ARPACK fue diseñado para ser eficiente para tareas como encontrar la solución de un sistema de ecuaciones Ax = b cuando A es escasa, como en su caso. Puede encontrar el código fuente de ARPACK here .

Answer 2

echa un vistazo a SuiteSparse, http://faculty.cse.tamu.edu/davis/suitesparse.html El autor, Tim Davis, es famoso en el campo de las matrices dispersas. También recibió un premio de Google por la calidad del código. El rendimiento de su código también es sobresaliente.

Aclamaciones