# Solving large sparse system

I am working on a problem with very large sparse matrices. I'd like to compute $$A^{-1} B$$, that is a crucial part of converting DAE to ODE (and there is no workaround). Here size of $$A$$ is 2E+5 x 2E+5 with 0.7% density (especially all diagonal entries are non-zeros). And $$B$$ is 2E+5 x 1E+6 whose density is only 1E-6.

I tried to solve this problem in python by scipy.sparse.spsolve, which eventually gave me a memory issue. The same issue happened while using an iterative solver. Unless I save each column to the disk and free the memory before computing the next column in the solution, the memory issue will persist.

So does anybody know a memory-efficient sparse solver or suggest a better way to solve this problem?

• Original question on StackOverflow – Anton Menshov Jul 9 at 14:32
• I think the question as asked is not really solvable. You have a nice sparse $\mathbf B$ that you are clobbering into a dense mess $\mathbf C = \mathbf A^{-1} \mathbf B$. You must pursue methods that do not require this .. something that can work with $\mathbf C$ in its present/operator form, as a cascade of sparse matvec and sparse backsolve. You might be able to gather ideas by posting a new question that presents the larger context. – rchilton1980 Jul 9 at 14:44
• Can you explain why you need to compute this full matrix? I'm speculating you have an ODE/DAE of the form $A y' = f(y) = B y$. Many ODE/DAE solvers are specifically designed to handle mass matrices and maintain sparsity. Explicit methods will require a linear solve with $A$ per stage, which can be done with an iterative method or a precomputed sparse decomposition. Implicit methods often involve linear solves with matrices of the form $A - h \gamma J$, where $J$ is the Jacobian of $f$. As long as $J$ isn't dense, this should avoid these memory issues. – Steven Roberts Jul 9 at 15:17
• Hi first thanks for all your replies. @StevenRoberts I would select some "coupled dofs" by investigating the non-zero pattern in this $A^{-1}B$ matrix. So that's why I need to know how resultant matrix looks like. – Lambda Z Jul 13 at 9:05

$$A^{-1}B$$ is a $$2\cdot 10^5 \times 1\cdot 10^6$$ matrix that is likely to be full, so the size of your output is 1.6 TB. I don't think you have other alternatives than writing it to the disk.
• Except, of course, to just never compute $A^{-1}B$ explicitly as a matrix :-) – Wolfgang Bangerth Jul 10 at 15:21