Efficiently solving large-scale linear systems is a critical challenge in electromagnetic simulations, particularly when using the Crank-Nicolson Finite-Difference Time-Domain method. Existing iterative solvers are commonly employed to handle the resulting sparse systems but suffer from slow convergence due to the ill-conditioned nature of the double-curl operator. Approximate preconditioners, like SOR and Incomplete LU decomposition (ILU) provide insufficient convergence, while direct solvers are impractical due to excessive memory requirements. To address this, we propose FlashMP, a novel preconditioning system that designs a subdomain exact solver based on discrete transforms. FlashMP provides an efficient GPU implementation that achieves multi-GPU scalability through domain decomposition. Evaluations on AMD MI60 GPU clusters (up to$\mathbf{1 0 0 0 ~ G P U s}$) show that FlashMP reduces iteration counts by up to$16 \times$and achieves speedups of$2.5 \times$to$4.9 \times$compared to baseline implementations in state-of-the-art libraries Hypre. Weak scalability tests show parallel efficiencies up to 84.1 %.

Accurately calculating the electronic structure of strongly correlated chemical systems necessitates a detailed description of both static and dynamical electron correlations, posing a significant challenge in ab initio quantum chemistry. Although the high memory and computational demands generally limit these calculations to relatively modest systems, the advanced computational capabilities of modern GPUs provide new avenues to expand these limits. However, complex control flows inherent to computation notably impair performance on GPUs. Furthermore, the significant disparity in computational load across different branches leads to load imbalance, challenging the large-scale simulations.