The kernel-free boundary integral (KFBI) method has successfully solved partial differential equations (PDEs) on irregular domains. Diverging from traditional boundary integral methods, the computation of boundary integrals in KFBI is executed through the resolution of equivalent simple interface problems on Cartesian grids, utilizing fast algorithms. While existing implementations of KFBI methods predominantly utilize CPU platforms, GPU architecture’s superior computational capabilities and extensive memory bandwidth offer an efficient resolution to computational bottlenecks. This paper delineates the algorithms adapted for both single-GPU and multiple-GPU applications. On a single GPU, assigning individual threads can control correction, interpolation, and jump calculations. The algorithm is expanded to multiple GPUs to enhance the processing of larger-scale problems. The arrowhead decomposition method is employed in multiple-GPU settings, ensuring optimal computational efficiency and load balancing. Numerical examples show that the proposed algorithm is second-order accurate and efficient. Single-GPU solver runs 50-200 times faster than traditional CPU, and the parallel efficiency of a multiple-GPU distributed solver within the same NUMA node reaches up to 80%.