I have been working on some problems on Gaussian processes and came across the formula $$ \int \mathcal{N}\left(a ; Bc, D\right)\mathcal{N}\left(c; e, F\right) dc = \mathcal{N} \left(a; Be, D+BFB^T\right), $$
where the densities are multivariate normal and vectors/matrices are of the appropriate size. It's not that clear to me (without getting into messy algebra) where this is coming from. Does anyone have a reference for this result?