Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It only takes a minute to sign up.
Sign up to join this community
Does Lapack have a routine that, given symmetric $A=A^T$ and $B$, computes the symmetric matrix $B^{-1}AB^{-T}$ (while preserving symmetry exactly)?
It would be enough to have this routine for triangular $B$, because of course the first step here will be computing a factorization of $B$ (LU or LDL, according to whether $B$ is symmetric, too).
Motivation: I think it would be useful in this question, especially if one factorizes the expression as suggested in the second part of my answer.
Unfortunately I don't think there's a great way to do this, at least not without some effort.
In the event that $\mathbf A$ is nonsingular, it might be useful to point out that the desired matrix $\mathbf C = \mathbf B^{-1} \mathbf A \mathbf B^{-T}$ is the inverse of the symmetric Schur complement $\mathbf S = \mathbf B^T \mathbf A^{-1} \mathbf B$. I find that BLAS/LAPACK has much better support for forming a matrix like $\mathbf S$ than one like $\mathbf C$, so maybe it's profitable to compute $\mathbf S$ first, then apply a symmetry-preserving inversion routine? The utility of this idea depends on the conditioning of $\mathbf A$ and $\mathbf B$.
If you're lucky enough that $\mathbf A$ is positive definite in addition to symmetric, you can use [potrf] to Cholesky factor/overwrite $\mathbf A = \mathbf L \mathbf L^T$. Then $\mathbf S = \mathbf B^T \left(\mathbf L \mathbf L^T \right)^{-1} \mathbf B = \left( \mathbf L^{-1} \mathbf B \right) ^T \left( \mathbf L^{-1} \mathbf B \right)$. Using [trsm], you can overwrite $\tilde{\mathbf B} = \mathbf L^{-1} \mathbf B$, then use [syrk] to compute $\mathbf S = \tilde{\mathbf B}^T \tilde{\mathbf B}$ into a temporary (or just overwrite $\mathbf A$). If $\mathbf B$ has full column rank, then $\mathbf S$ is also positive definite, and the desired output $\mathbf C = \mathbf S^{-1}$ can be computed using [potrf] followed by [potri]. The symmetry of $\mathbf A$, $\mathbf S$ and $\mathbf C$ is enforced explicitly at the API level, because all of [syrk], [potrf] and [potri] operate only upon a single triangle (with the opposing triangle assumed to match).
When $\mathbf A$ is invertible but indefinite, you can follow the same basic idea but computing $\mathbf S$ is more involved due to pivoting considerations. If you'd like to see that procedure too, just leave a comment and I will add it as an edit.
