As a foundation for numerical solvers in computational electromagnetics, particularly for multiphysics and electromagnetic compatibility applications, Jefimenko’s equations offer a generalized solution to Maxwell’s equations, enabling the direct computation of electromagnetic fields from time-dependent source distributions without the boundary-condition artifacts inherent to grid-based methods. However, the numerical integration of these equations is computationally intensive, typically scaling as $\mathcal{O}\left(N_{\mathrm{s}}{N}_{\mathrm{o}}\right)$ for $N_{\mathrm{s}}$ source points and $N_{\mathrm{o}}$ observation points. In this paper, we present JefiFast, a highly optimized graphics processing unit (GPU) implementation that significantly outperforms the state-of-the-art JefiGPU algorithm. We identify that previous implementations are strictly memory-bound due to inefficient global memory transactions and a lack of data reuse. JefiFast addresses these bottlenecks through four key optimizations: (i) a packed memory layout (PML) using an array-of-structures approach to ensure coalesced memory access for source densities and their derivatives; (ii) geometry-aware shared memory tiling strategies that maximize L2 (level-2) cache hit rates and on-chip data reuse; (iii) pre-computation of time derivatives to minimize redundant arithmetic operations; and (iv) a robust observation domain decomposition strategy that enables linear scaling across multiple GPUs. Benchmarks demonstrate that JefiFast achieves speedups ranging from $4.08$ times (for ${30}^3$ grids on a single NVIDIA V100 graphic processor) to $84.51$ times (for ${50}^3$ grids on 4 NVIDIA V100 processors) compared to the baseline. Notably, for a ${50}^3$ grid on a single GPU, JefiFast reduces execution time from about 51 min to just about 2.6 min ($19.54$ times speedup). These performance advances make high-resolution relativistic heavy-ion collision simulations feasible in near real-time.