  The Hartree-Fock-Bogoliubov (HFB) theory is the starting point for treating
superconducting systems. However, the computational cost for solving large
scale HFB equations can be much larger than that of the Hartree-Fock equations,
particularly when the Hamiltonian matrix is sparse, and the number of electrons
$N$ is relatively small compared to the matrix size $N_{b}$. We first provide a
concise and relatively self-contained review of the HFB theory for general
finite sized quantum systems, with special focus on the treatment of spin
symmetries from a linear algebra perspective. We then demonstrate that the pole
expansion and selected inversion (PEXSI) method can be particularly well suited
for solving large scale HFB equations. For a Hubbard-type Hamiltonian, the cost
of PEXSI is at most $\Or(N_b^2)$ for both gapped and gapless systems, which can
be significantly faster than the standard cubic scaling diagonalization
methods. We show that PEXSI can solve a two-dimensional Hubbard-Hofstadter
model with $N_b$ up to $2.88\times 10^6$, and the wall clock time is less than
$100$ s using $17280$ CPU cores. This enables the simulation of physical
systems under experimentally realizable magnetic fields, which cannot be
otherwise simulated with smaller systems.
