Consider a system of the following form:
$$(A+K)x=b$$
where $A$ is symmetric, positive definite and block diagonal (in fact, a block diagonal matrix made of stiffness matrices arising from FEM discretizations at multiple timesteps, grouped together). Multiplying our linear system by $P^t=K^tA^{-1}+I$, we get $$(K+K^t + K^tA^{-1}K+A)x = P^tb$$
We assume that $K$ is non-symmetric, has full rank and that $K+K^t$ is positive definite, so that this is a symmetric and positive definite system of equations. The only inconvenient here is the presence of $A^{-1}$, and the fact that one needs to solve a (block diagonal) system to obtain $P^tb$.
Yet another way is to further rewrite everything recognizing a Schur complement:
$$ \left\lbrack \begin{array}{cc} A& -K^t\\ -K^t& -K-K^t-A \end{array} \right\rbrack \left\lbrack \begin{array}{c} p\\ x\end{array} \right\rbrack= \left\lbrack \begin{array}{c}-b\\ -b\end{array}\right\rbrack $$
Which one of the three systems would you rather work with, and with which solver would you work? I am interested in wall-clock time and parallelizability
We assume our blocks to be large, e.g. of size $1e8-1e10$.
To connect the first two formulations, an idea might be to solve formulation one with pre-conditioned GMRES, and use $P^t$ as preconditioner. Applying $P^t$ can be done with the (preconditioned) conjugate gradient, and this is highly parallelizable since $A$ has a diagonal block structure.
But in this way, since we are solving essentially formulation two, we'd still need to come up with a good preconditioner for it.
Edit. Also because formulation two seems much worse conditioned than formulation one, as numerical experiments show.
