Question - Large Matrices - Limitations

[saisandeepdammati]

Howdy!

I am trying to solve a higher-order elliptic equation using LU decomposition on a 3D uniform grid. The grid is N^3 in size with 'N' being the number of cells along on direction and the grid is uniformly distributed across several MPI processors. On decomposing the elliptic equation using finite difference, it results in a sparse matrix A of size (N^3 x N3), b matrix of size [N^3 x 1] and solution matrix, x of size [N^3 x 1]. I am want to solve this system using SuperLU_dist, however, I see an issue when N becomes large say N^3 = 2048^3. The sparse matrix 'A' in SuperLU, has member variables (say fst_row) of type int and the max value that one can store in a 4 byte int is 2 billion. In a sparse matrix of size N^3 = 2048^3, the number of rows in A can be easily more than 2 billion and lead to integer overflow. I have the following questions:

1. Is there a work around for this? In particular, can 'fst_row' variables in 'A' be changed to say long int?

2. Can SuperLU_dist handle sparse matrix sizes in the range of (2048^3 x 2048^3) to (8192^3 x 8192^3)?

Thank you very much.

[xiaoyeli]

For such a big problem, you can not use SuperLU (or any explicit matrix solver). Even if SuperLU allows 64-bit int for N, the MPI communication layer cannot communicate a vector of size larger than 2 billion.

You should think about developing a matrix-free solver, mostly likely an iterative solver, that only needs matrix vector multiplication, but without forming A matrix explicit. You can also use domain decomposition, and use SuperLU (or other explicit solver) to solve the smaller subdomain problems as a preconditioner.

[saisandeepdammati]

Howdy!

Thank you very much for your insightful response. In my case, the sparse matrix 'A' does not change in my CFD simulation and only the RHS vector 'b' changes, so therefore, in principle I can perform LU decomposition on 'A' only once at the beginning of the simulation and store the 'L' and 'U' matrices and then use them to perform back substitution on 'b' when needed to obtain my solution vector 'x'.

This way, I need to perform LU decomposition of 'A' only once at the beginning. Even if this LU decomposition is slow, since it is performed only once, it should be fine. So, my only question is if SuperLU_dist can handle a sparse matrix of size (1 billion x 1 billion) distributed equally across say 100,000 cores?

[xiaoyeli]

Theoretically possible, but for such a large matrix, the amount of fill in the L & U factors will be huge, making it impractical.
For a typical 3D kxkxk discretized PDE, using nested dissection ordering, the number of fill-ins is O(k^4), the flop count is O(k^6).