Studying the optoelectronic structure of materials can require the computation of several thousands of the smallest positive eigenpairs of a pseudo-hermitian Hamiltonian. Iterative eigensolvers may be preferred over direct methods for this task since their complexity is a function of the desired fraction of the spectrum. In addition, they generally rely on highly optimized and scalable kernels such as matrix-vector multiplications that leverage the massive parallelism and the computational power of modern exascale systems. The Chebyshev Accelerated Subspace iteration Eigensolver (ChASE) is able to compute several thousands of the most extreme eigenpairs of dense hermitian matrices with proven scalability over massive parallel accelerated clusters. This work presents an extension of ChASE to solve for a portion of the smallest positive eigenpairs of pseudo-hermitian Hamiltonians as they appear in the treatment of excitonic materials. By exploiting the numerical structure and spectral properties of the Hamiltonian matrix, we preserve the characteristic positive-negative symmetry in the treatment of the eigenvectors and propose an oblique variant of Rayleigh-Ritz projection that features quadratic convergence of the Ritz values with no explicit construction of the dual basis. Additionally, we introduce a parallel implementation of the recursive matrix-product operation appearing in the Chebyshev filter with limited amount of global communications. Our development is supported by a full numerical analysis and experimental tests.