The linear response eigenvalue problem (LREP), which arises from many scientific and engineering fields, is quite challenging to solve numerically, especially when it has zero eigenvalues and the discrete sparse/dense system is large-scale. Based on a direct sum decomposition of biorthogonal invariant subspaces, using the structure of generalized nullspace, we propose a Bi-Orthogonal Structure-Preserving subspace iterative solver, which is stable, efficient, and of excellent parallel scalability. The biorthogonality is of essential importance and preserved well by a modified Gram-Schmidt biorthogonalization (MGS-Biorth) algorithm on the discrete level. We naturally deflate out the already converged eigenvectors by searching the next eigenpairs in the biorthogonal complementary subspace without introducing artificial parameters. When the number of requested eigenpairs, denoted as nev, is very large, we propose a moving mechanism to compute the eigenpairs batch by batch so that the projection matrix size is small and independent of nev. The Bogoliubov-de Gennes equations (BdG) is a special LREP that arises from the perturbation of the Bose-Einstein Condensates (BEC) around stationary states. Our algorithm is applied to the BdG of rotating BEC and of spin-1 BEC. Ample numerical examples are provided to demonstrate the stability, efficiency, and parallel scalability.