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.