The modified Poisson-Boltzmann (MPB) equations are often used to describe the equilibrium particle distribution of ionic systems. In this work, we propose a fast algorithm to solve the MPB equations with the self Green's function as the self-energy in two/three dimensions, where the solution of the self Green's function poses a computational bottleneck due to the requirement of solving a high-dimensional partial differential equation. Our algorithm combines the selected inversion with hierarchical interpolative factorization for the self Green's function. By strategically exploiting the locality and low-rank characteristics of the corresponding operators, algorithms with linear complexity in two dimensions and quasi-linear complexity in three dimensions are achieved. Extensive numerical results are conducted to demonstrate the accuracy and efficiency of the proposed algorithm for problems in two/three dimensions.
