I've always found the field split options in PETSc to be quite a headache, and haven't had much success with it myself.
The link you give was actually me asking that question, and, well, I never ended up doing this, but opted for 2 other strategies instead:

1. Try -ksp_type minres and hope for the best. You may need to scale the equations a bit (pick units so that the matrix isn't to badly scaled). No guarantee that this will work at all, but I have had luck occasionally.

2. Problem specific, but, if there are relatively few lagrange multipliers, it's doable;

K*x_a + C*x_b = r_a
C'*x_a = r_b

==> 
x_a  = K^-1 * (r_a - C*x_b)
==> 
C' * K^-1 * (r_a - C*x_b) = r_b
==>
C' * K^-1 * C*x_b = C' * K^-1 * r_a - r_b

Solve:
D = K^-1 * C

Use the fact that
D' = C' * K^-1
and then you can solve for x_b
D' * C *x_b = D' * r_a - r_b
and then for x_a
x_a  = K^-1 * (r_a - C*x_b)

This method is used in OOFEM when it comes to solving the tangent sensitivity problem in the homogenization B.C.s (in the case of the Neumann type homogenization b.c. you get at most 6 RHS's to solve for, and 6 lagrange multipliers).

If you have hundreds+ of lagrange multipliers, this is not an option, because solving D becomes prohibitively expensive.

3. The issues (no support for this operation for this object type!) i suspect is related to the matrix type.
The zero-diagonals defined is another issue; many of the solvers in PETSc expect the 0 alog the diagonal to be allocated, even though it's not used. You may have some luck with -mat_type and maybe -mat_shift

Though, I think, for any of this to be successfull, OOFEM needs to do something a lot more involved, perhaps the Schur-complement code needs to be implemented in OOFEM in some more or less clever way. Similarly, geometric multigrid. But these are not simple things to implement, and put a lot of requirement in the input as well.