The calculation of overlaps between many-electron wave functions at different nuclear geometries during nonadiabatic dynamics simulations requires the evaluation of a large number of determinants of matrices that differ only in a few rows/columns. While this calculation is fast for small systems, its cost grows faster than the alternative electronic structure calculation used to obtain the wave functions. For wave functions that can be written as a CIS expansion, all determinants can be computed using the set of level-2 minors of the reference matrix. However, this is still a costly computation for large systems. In this paper, we provide an algorithm for efficiently calculating all level-2 minors of a matrix by reutilizing and updating the LU factorization for the determinants of the minors. This approach results in a parallel version of the algorithm that is up to an order of magnitude faster then the current best parallel implementation. The algorithm thus allows the computation of exact wave function overlaps for relatively large systems, with a high density of states, at virtually no cost compared with the electronic structure calculations. Furthermore, the new algorithm opens the path to further investigations in efficient computing of the exact wave function overlaps for complex wave functions such as MR-CIS and MR-CISD.