I am interested in knowing some efficient techniques for solving the following extended Saddle point problem. \begin{align} \begin{bmatrix} A & B^T & C^T \\ B & 0 & 0 \\ C & 0 & 0 \\ \end{bmatrix} \begin{Bmatrix} x \\ y \\ z \\ \end{Bmatrix} = \begin{Bmatrix} F_x \\ F_y \\ F_z \\ \end{Bmatrix} \end{align}

x is velocity, y is pressure, and z is Lagrange multipliers for imposing interface conditions on immersed solids. Such a matrix system is obtained using the fictitious domain/distributed Lagrange multiplier method for incompressible flow problems.

I am looking at two different schemes. In both the schemes, matrices A and B remain the same at each time step but for matrix C the values, as well as the locations of non-zero entries, change from one time step to the other. The number of rows of matrix C is significantly small compared to those of A and B.

Scheme 1: Matrix A is diagonal (and invertible). Since A is diagonal, I am employing a Schur complement based solver using SuperLU solver from the Eigen library. I am wondering if the efficiency can be improved further since only matrix C, which is significantly smaller in size compared to B, changes at each time step.

Scheme 2: Matrix A is sparse and symmetric (and invertible). At the moment I use the sparse direct solver PARDISO. This is certainly not an efficient approach since matrices A and B remain the same. I do not know what type of Schur complement approach would be efficient.

I am looking for some alternative strategies that can be efficient. I very much appreciate any inputs in this regard.

Note: I tried Uzawa type schemes in the past. They are neither efficient nor robust for generic problems; I do not like to fiddle with the relaxation parameters in the Uzawa scheme and its variants.

  • 1
    $\begingroup$ What is the fluid flow problem you are solving? Darcy? Stokes? Can you give more information about the continuous problem? $\endgroup$ Commented Mar 31, 2021 at 23:35
  • $\begingroup$ Hi @AbdullahAliSivas. I am solving incompressible Navier-Stokes. For scheme 1, viscous and convection terms are treated explicitly. For scheme2, convection is treated explicitly and the viscous term is treated implicitly. $\endgroup$
    – Chenna K
    Commented Apr 1, 2021 at 10:00
  • $\begingroup$ You are open to use iterative solvers, right? $\endgroup$ Commented Apr 1, 2021 at 14:30
  • $\begingroup$ @AbdullahAliSivas. Yes. Absolutely. $\endgroup$
    – Chenna K
    Commented Apr 1, 2021 at 14:58

2 Answers 2


For both schemes: I am going to lump $B$ and $C$ into a single matrix $B$ for convenience of communicating the idea. Likewise, $y$ and $z$ are now $y$, and $F_y$ and $F_z$ are $F_y$.

$Ax+B^Ty=F_x$ or $x = A^{-1}(F_x-B^Ty)$

$Bx=F_y$ or $BA^{-1}(F_x-B^Ty)=F_y$


Note that $BA^{-1}B^T$ is symmetric and invertible, and can be factorized using a cholesky decomposition. To construct $BA^{-1}B^T$ (assuming the $A$ is amenable to a cholesky decomposition):

$$BA^{-1}B^T = BL^{-T}L^{-1}B^T=(L^{-1}B^T)^T(L^{-1}B^T)$$ which is amenable to a symmetric Rank-K update, which is fast, efficient, and present in most linear algebra libraries. In BLAS, it is usually called something like syrk; in Eigen, it is called by invoking selfAdjointView<Lower>().rankUpdate(). If $A$ is diagonal, this is much simpler.

EDIT: Regarding your comment that $B$ is actually quite large:

$$ \begin{pmatrix} B \\ C \\ \end{pmatrix}A^{-1} \begin{pmatrix} B^T C^T \\ \end{pmatrix}= \begin{pmatrix} BA^{-1}B^T & BA^{-1}C^T \\ CA^{-1}B^T & CA^{-1}C^T \\ \end{pmatrix} $$

The top left block of this matrix will be constant. The top right is redundant (in a lower triangular storage scheme). If you use syrk to update the top left and bottom right corners, they will store and compute only the lower triangular portion, which saves a lot of time and space.

You can either accept the resulting block matrix and solve it, or you can use a block inversion scheme to determine the action of that inverse using the cholesky factorizations of the top left and bottom right corners. This may be useful since you only have to factorize the top left corner once.

  • $\begingroup$ Hi @Charlie S, thank you for your comment. Dense matrices are not an option since the B matrix is of considerable size; can go up to or even beyond 100,000 in 3D. $\endgroup$
    – Chenna K
    Commented Apr 1, 2021 at 10:07
  • $\begingroup$ The approach you described is exactly what I use for scheme 1. It works decently with sparse matrices in Eigen. The only issue is that I need to recompute the Schur complement at every time step since matrix C changes at every time step. I am curious to know if I can improve the performance as matrix B remains the same. $\endgroup$
    – Chenna K
    Commented Apr 1, 2021 at 10:55
  • $\begingroup$ That is a very impressive $B$! I was hoping this was not the case, but indeed you can exploit the constant nature of $B$ to save a lot of time, at the expense of a giant headache in coding the solution! I have edited my answer to reflect this. $\endgroup$
    – Charlie S
    Commented Apr 1, 2021 at 13:56
  • $\begingroup$ Hi @Charlie S, I had derived the approach you described in the edited version but didn't know how to implement it efficiently. Your comment pointed me in the right direction of rank update. I will explore this further. Thank you very much for your help! $\endgroup$
    – Chenna K
    Commented Apr 2, 2021 at 8:14
  • 1
    $\begingroup$ @ChennaK In Eigen, I don't think the PardisoSupport module exposes the the ability to apply $L^{-1}$, which is unfortunate. The 5th argument in internal::pardiso_run_selector<StorageIndex>::run is the phase argument, which controls the solve operation (33 is solve, 331 is forward substitution, 333 is backward). You may have to extend the PardisoSupport module or write your own wrappers to use it effectively. $\endgroup$
    – Charlie S
    Commented Apr 2, 2021 at 14:14

Disclaimer: I wrote this answer in a rush, may not be up to the standards. Also, I don't have enough information about your problem to give more detailed advice. However, the information presented here should be enough to get you going.

I will suggest you to use MINRES with an appropriate preconditioner.

For a moment, let's assume that this is not a twofold saddle point problem, and you have the following problem to solve

$\begin{equation} \begin{bmatrix} A & B^T \\ B & 0 \end{bmatrix} \begin{bmatrix} U \\ P \end{bmatrix} = \begin{bmatrix} F \\ G \end{bmatrix}. \end{equation}$

I am assuming that you are using an inf-sup stable discretization of the Navier-Stokes equations, (for scheme 2) which would imply that $BA^{-1}B^T$ and $Q_p$, the mass matrix defined on the pressure space, are spectrally equivalent. Moreover, you can prove that $Q_p$ and $diag(Q_p)$ are spectrally equivalent, and spectral equivalence is a transitive relation. Also, with scheme 2, $A$ is just a vector Laplacian, hence you could use (algebraic) multigrid to very efficiently invert it. Such a preconditioner can be described through the action of its inverse as

$P^{-1} = \begin{bmatrix} AMG(A) & 0 \\ B & diag(Q_p)^{-1} \end{bmatrix}.$

Note that I am using $P^{-1}$ to mean as the action of the inverse of

$P= \begin{bmatrix} A & 0 \\ B & diag(Q_p) \end{bmatrix}.$

If $diag(Q_p)$ is too slow in terms of iterations, you can use $AMG(Q_p)$ to obtain a better preconditioner (in terms of iterations) at slightly higher computational cost.

Now, we need to deal with the second saddle point problem. I am not familiar with immersed solid literature, but I think there must be an inf-sup condition for that part too. Then by using Theorem 3.1 of (Howell JS, Walkington NJ (2011) Inf-sup conditions for twofold saddle point problems.), you can find another spectrally equivalent matrix -potentially a mass matrix depending on problem parameters-, let's call it $Z$. Finally, the following preconditioner probably would be optimal for your problem

$P= \begin{bmatrix} A & 0 & 0\\ B & diag(Q_p) & 0\\ C & 0 & Z \end{bmatrix}.$

Depending on the form of $Z$, you may be able to approximate it with its diagonal part or using an algebraic multigrid method (I suggest BoomerAMG). Or maybe if the nonzero structure of $Z$ is nice, for example block-diagonal with small size diagonal blocks, you can use a direct solver at low computational cost.

  • 1
    $\begingroup$ Hi @Abdullah Ali Sivas, appreciate your inputs. The only thing that stops me from using this approach is the lack of a mass-type matrix for the Lagrange multiplier field. I will think about if I can come up with something similar. $\endgroup$
    – Chenna K
    Commented Apr 3, 2021 at 13:00
  • $\begingroup$ Depending on the space of the Lagrange multipliers, there may be one. At least chose $Z=hI$ and give it a go. $\endgroup$ Commented Apr 4, 2021 at 1:32

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.