Recently, the escalating computational capability provided by advanced GPU technology for numerical simulations is well-suited for tackling large-scale engineering challenges. The lattice Boltzmann flux solver (LBFS), as a relatively new fluid-solving method, combines the parallelization characteristics of the lattice Boltzmann method (LBM) with the ability to handle non-uniform grids. Building upon these advantages, this study aligns its flux computational steps with the demands of the GPU parallel computing environment, thus developing an efficient flux-reconstructed lattice Boltzmann flux solver (FRLBFS), the improved algorithm not only ensures precision but also achieves a significant enhancement in efficiency, reaching up to a remarkable 500-fold improvement. Additionally, this work combines the immersed boundary method, significantly improving the efficiency of addressing hydrodynamic problems with multiple structures. Through numerical validation, under identical GPU hardware conditions, the proposed method in this paper outperform LBM in terms of computational efficiency. Lastly, simulations of flow past an array of eight cylinders at different Reynolds numbers are conducted, instantaneous contour plots at various dimensionless time points and time-dependent curves of drag and lift coefficients for different cylinders are provided, to demonstrate the complex vortex shedding surrounding multiple cylinders and its extraordinary computational efficiency.