I've been struggling with some simulations that involve Lagrange multipliers, specifically the Neumann type boundary condition for homogenization (3, 6, 9 Lagrange multipliers)

So, these problems generally look like:

[K G] [u] = [f]
[G' C] [y]   [g]

where g and C are often zero (but we need not concern ourselves with this)
Here, y, is a few Lagrange multipliers.

The main problem is of course that even a single Lagrange multiplier messes up the system, and forces sub-optimal solvers.

K * u + G * y = f   ==>  u = K^(-1) * [ f - G*y ]

G' * u + C * y = g  ==>   - [G' * K^(-1) * G + C] * y = g - G' * K^(-1) * f 

So, we need to solve:

K^(-1) * f
K^(-1) * G

so the cost depends directly on the number of columns of G. So, for my current case, with just 3 Lagrange multipliers, it's just 3 extra KSP-steps (as preconditioner can be reused).

I'm leaning towards adding an input to NRSolver to make use of block-matrix structure to explicitly solve the linear system in this manner. Alternatively, a sparse linear solver type that does this as a "wrapper".

Some concerns with the NRSolver approach (which might still be the easiest):
If we assemble the full system, then take out K and G from this system, we basically double the memory consumption for storing these temporary sub-matrices. However, our infrastructure for giving the stiffness matrix pointer to NRSolver, building the internal structure, etc. doesn't really allow us to change this easily.

--------------

An alternative is to just set up Schur-complement with PETSc. But I struggle with the details here.

I have tried for many hours of getting something useful from the Schur-complement methods in PETSc, but with no success. I'm testing it out on a relatively small problem, ~5k unknowns with 3 (global) Lagrange multipliers.
I have made use of the following flags:

-pc_type fieldsplit -pc_fieldsplit_detect_saddle_point -pc_fieldsplit_schur

(which are essentially required), and then there are lots and lots of other flags to test

Many things I just can't get to work at all, even though they are suggested in the manual. For example the option

-fieldsplit_1_pc_type lsc

is just not useable at all, as it generates an error that MatMatMult is not supported for the matrix type (which is called by PCSetUp_LSC()).

I have managed to obtain answers, but they are 10-100 times slower than just the simple;

-pc_type none -ksp_type minres -mat_type sbaij

But there are dozens of flags to use, and I can barely get any of them to work. But, there is no guarantee that there are good ways to solve this at all, since the Schur-complement stuff in PETSc assumes that there are many Lagrange multipliers.