@ioannisPApapadopoulos @KarsKnook
I've been playing with QL instead of reverse Cholesky for non-SPD operators. It suffers from fill in in the top row of blocks so the complexity would be $O(np+n^3)$.
BUT if we did strong form (so non-symmetric matrices) I believe everything apart from the first block column will be diagonal. Then QL will be optimal complexity.
Note what we need is ContinuousPolynomial{-1} which would be P^(1,1) instead of Legendre but also with delta functions between the elements.