# Viable algorithms for efficiently solving a block matrix system with non-uniform sparsity structure

I am trying to use Newton's method to get a stationary solution for a system of equations of the following form:

$$\begin{Bmatrix} \frac{\partial x}{\partial t} \\ 0 \end{Bmatrix} = \begin{Bmatrix} f(x, y) \\ g(x, y) \end{Bmatrix}$$

For vectors $$x$$ and $$y$$.

Due to the nature of the equations I'm trying to solve, the Jacobian can be expressed in the form:

$$J = \begin{bmatrix} \frac{\partial f}{\partial x} & \frac{\partial f}{\partial y} \\ \frac{\partial g}{\partial x} & \frac{\partial g}{\partial y} \end{bmatrix} = \begin{bmatrix} A & B \\ C & D \end{bmatrix}$$

Such that $$A$$ and $$C$$ are dense matrices, while $$B$$ and $$D$$ are sparse. It should also be noted that the elements in the diagonals of $$A$$ and $$C$$ are much larger than those in $$B$$ and $$D$$ (by ~100 times).

Newton's method works perfectly when I solve the linear problem using Gaussian elimination, and the eigenvalues of the Jacobian are such that no severe numerical problems are faced with that algorithm.

I need to use an iterative method for larger problems, however, and the linear problem isn't appropriately solved with GMRES, IDR(s) or biconjugate gradient stabilized method when the matrix is taken as a whole, so I suppose another technique or precaution could be necessary.

I have tried preconditioning with an incomplete LU decomposition and with a diagonal preconditioner, and failed with both. When individually solving for blocks $$A$$ and $$D$$, however, one notices that the linear system $$A\Delta x = - f(x, y)$$ exhibits the same behaviour as the complete Jacobian, while $$D \Delta y = - g(x, y)$$ is easily solved by any of the Krylov subspace algorithms I cited above.

Is there any algorithm that could help solve the linear system?

• Have you tried using Schur’s complement, the version involving the inverse of D? Jun 15 at 19:33
• @AmitHochman, in order to use Schur's complement, I'd have to compute the product $BD^{-1}C$, for which I'd have to invert D, right? I'd like to avoid that if possible, since I can't guarantee that D is tridiagonal (or of any other structure that allows for O(n) inversion) for all cases Jun 15 at 20:07
• You don’t have to invert D, but you would have to solve DX=B, or D^TX=C, whichever is easier. Jun 15 at 20:24
• There are methods, such as Krylov subspace recycling, which you could use to lower the cost of solving the same system with different right hand sides. Jun 15 at 20:27
• The form of your system is called in literature differential algebraic equation DAE. This might help in your literature search.
– Bort
Jun 16 at 12:25