Consider a system of the following form:


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.

  • $\begingroup$ Could you elaborate and write out the structure of Pt? $\endgroup$ Commented Mar 16, 2023 at 13:53
  • $\begingroup$ @rchilton1980of course, I edited the question $\endgroup$
    – Lilla
    Commented Mar 16, 2023 at 14:05
  • $\begingroup$ Guessing K is not symmetric? $\endgroup$ Commented Mar 16, 2023 at 14:07
  • $\begingroup$ No, $K$ has a block structure like the time matrix from the implicit Euler method. And the blocks are not constant, and not symmetric. @rchilton1980 $\endgroup$
    – Lilla
    Commented Mar 16, 2023 at 14:52


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.