We present an efficient program for the exact diagonalization solution of the pairing Hamiltonian in spherical systems with rotational invariance based on the SU(2) quasi-spin algebra. The basis vectors with quasi-spin symmetry considered are generated by using an iterative algorithm. Then the Hamiltonian matrix constructed on this basis is diagonalized with the Lanczos algorithm. All non-zero matrix elements of the Hamiltonian matrix are evaluated "on the fly" by the scattering operator and hash search acting on the basis vectors. The OpenMP parallel program thus developed, PairDiagSph, can efficiently calculate the ground-state eigenvalue and eigenvector of general spherical pairing Hamiltonians. Systems with dimension up to 10(8) can be calculated in few hours on standard desktop computers. Program summary Program Title: PairDiagSph. CPC Library link to program files: https://doi.org/10.17632/t5txsg3s5g.1 Code Ocean capsule: https://codeocean.com/capsule/9159628 Licensing provisions: CC by NC 3.0. Programming language: Fortran 95. Nature of problem: The exact diagonalization of spherical pairing Hamiltonian can be achieved in the quasi-spin space. Solution method: The program generates the basis vectors via the adjacency excitation algorithm, and diagonalizes the spherical pairing Hamiltonian by the Lanczos + QR algorithm. Additional comments including restrictions and unusual features: The total number of spherical must be less than 64; The maximal dimension that can be handled is restricted by the local RAM capacity.