I have contacted a friend who suggested some algorithms to try out for solving the linear system:
algebraic brute force (does not address the nonlinearity, reported to work on moderately sized problems): gmres with ILU(0) preconditioning.
If ILU(0) does not work:
We are solving K*x = b where K has a block structure: K = | A B | | B^T -C | (C can be positive semi-definite or zero - our case)
Instead, an "augmented Lagrangian" technique would lead to Qx = b with Q = | A + alphaB*B^t B | | B^T -C | alpha>=0, small
Either the augmented system can be solved, or Q could be used in ILU(0) preconditioner.