Thank you for the suggestion! OK so I’m not very fresh on topological spaces (compactness and all that), so I’ll try to avoid discussing these abstract concepts, but don’t hesitate to tell me if I should think more about it and to correct me if I’m wrong. Again, I’ll lay out a detailed take on what I’m doing or trying to do, so apologies for the length.

I’ll start from the equation for the Jacobian that I already had above:

J(x) = \begin{bmatrix}
(\partial_1 G_1)(x) & (\partial_2 G_1)(x) & \cdots & (\partial_N G_1)(x)\\
(\partial_1 G_2)(x) & (\partial_2 G_2)(x) & \cdots & (\partial_N G_2)(x)\\
\vdots & & \ddots & \vdots \\
(\partial_1 G_N)(x) & (\partial_2 G_N)(x) & \cdots & (\partial_N G_N)(x)\\
\end{bmatrix}
- \begin{bmatrix}
T_1 & & & \\
& T_2 & & \\
& & \ddots & \\
& & & T_N\\
\end{bmatrix}

where the G_i(x_1,\ldots,x_n) are functions of the vectors x_i and the \partial_iG_j are the corresponding (sub-block) Jacobians, which are all diagonals. The T_i are matrices that, as you say, represent differential operators (but discretized onto a finite 3D grid and vectorized (as in, rearranged from 3D to 1D)).

A sketch of what it might look like in ASCII art :

```
⠑⢄ ⠑⢄ ⠳⣄
J(x) = ⠑⢄ ⠑⢄ ⠑⢄ + ⢿⣷
⠑⢄ ⠑⢄ ⠑⢄
= ∇G(x) + T
```

I’m assuming there was a bit of confusion with the names, but I guess by G you meant the first term in that equation above (that I write as `∇G(x)`

). Assuming that’s correct, then I think T might fit your criteria, although my T contains at least one T_i block on the diagonal that is divergence-free, such that `T_i * ones(n) .== 0`

(where `n = size(T_i, 1)`

). This means the null space of T_i — and thus that of T — is not trivial, which I think might be a problem for doing LU (or iLU too?) on T. (It is thanks to the diagonal terms that arise from the \partial_iG_i that the null space of the whole Jacobian J(x) is actually trivial and allows me to “invert” it.)

However, I think I can get away with a preconditioner by rearranging things a bit. Indeed, I can remove a constant piece from the \partial_i G_j and move them into T such that it is invertible and amenable to iLU preconditioning. That’s because the blocks for which the T_i are non-invertible also come with a corresponding G_i that has a linear term. So I can write them as G_i(x_1,\ldots,x_N) = \tilde{G}_i(x_1, \ldots, x_N) + \lambda_i x_i, where the \lambda_i x_i term has Jacobian \lambda_i I. Sticking this \lambda_i I into a new \tilde{T}_i = T_i - \lambda_i I allows me to rewrite J(x) as

J(x) = \begin{bmatrix}
(\partial_1 \tilde{G}_1)(x) & (\partial_2 G_1)(x) & \cdots & (\partial_N G_1)(x)\\
(\partial_1 G_2)(x) & (\partial_2 \tilde{G}_2)(x) & \cdots & (\partial_N G_2)(x)\\
\vdots & & \ddots & \vdots \\
(\partial_1 G_N)(x) & (\partial_2 G_N)(x) & \cdots & (\partial_N \tilde{G}_N)(x)\\
\end{bmatrix}
- \begin{bmatrix}
\tilde{T}_1 & & & \\
& \tilde{T}_2 & & \\
& & \ddots & \\
& & & \tilde{T}_N\\
\end{bmatrix}

where the \tilde{T}_i are all invertible so iLU should work (I think). But I’m not sure:

- Are the criteria that you mention met? (Do they strictly need to be?)
- Is \nabla\tilde{G}(x) compact?
- Is the resolvent of T compact? (I don’t know what the resolvent is actually so I’d love a small explanation if possible, otherwise, I have been trying to go to bed for 1hr+ now so I’ll leave it at that and come back to this later )

Another question comes to mind though. When LU + backsolve work (i.e., when the matrices are small or sparse enough and I have enough memory), then shouldn’t I stick with that instead of going for Krylov? I thought GMRES would be slower in this case. (FWIW, GMRES and co. are part of composability that I am looking for in the future, but I don’t think I need it right now.)