To address issues such as insufficient accuracy of traditional sound source localization algorithms in complex environments, high complexity, and slow convergence of sparse Bayesian solvers, a parallel subarray block sparse Bayesian learning (PS-BSBL) algorithm for sound source localization is proposed. An offline subarray partitioning (OSP) algorithm is designed, which divides microphone arrays and spatial grids into matched subarrays and subgrid spaces through covariance matrix adaptive evolution strategy (CMA-ES)-optimized array weight design and Leiden graph clustering algorithm. Each subarray solves its corresponding spatial problem in parallel, reducing the global solution scale. The PS-BSBL framework is constructed, employing inverse Gaussian priors and Gamma priors to model sparse hyperparameters and noise precision, respectively, with variational Bayesian methods for posterior approximation. For scenarios where Bessel functions cannot be eliminated, a multi-step inertial Newton iteration (MSI) is proposed to accelerate hyperparameter updates, enhancing pruning efficiency and convergence speed. Experimental results show that the proposed OSP algorithm reduces the average running time of fast variational block sparse Bayesian learning (F-BSBL), expectation-maximization BSBL (BSBL-EM), and weighted BSBL (W-BSBL) by 97.61%, 87.63%, and 87.82%, respectively， compared to block diagonal approximation algorithms, while normalized mean squared error (NMSE) decreases by an average of 86.53%, 16.64%, and 24.13%. The PS-BSBL algorithm improves convergence speed by 99.80% and 99.85% over group sparse fast marginalization BSBL (g-FMLM) and fast marginalization BSBL (BSBL-FM), respectively, with NMSE reductions of 51.30% and 8.92%, respectively. Under signal-to-noise ratios from -10 to 20 dB, PS-BSBL achieves lower angle root mean squared error (RMSE) than L1 norm and singular value decomposition-based sparse reconstruction (L1SVD), off-grid sparse Bayesian learning (OGSBL), convolutional neural network-based (CNN), and residual network-based (ResNet) methods by averages of 65.11%, 20.86%, 81.59%, and 78.45%, respectively, demonstrating superior robustness in low SNR, underdetermined, and large-scale array scenarios.