We present an algorithm to reduce the computational complexity for plane-wave codes used in electronic structure calculations. The proposed algorithm avoids the diagonalization of large Hermitian matrices arising in such problems. The computational time for the diagonalization procedure typically grows as the cube of the number of atoms, or the number of eigenvalues required. To reduce this computational demand, we approximate directly the occupation operator corresponding to the eigenvectors associated with the occupied states in a certain subspace without actually computing these eigenvectors. A smoothed Chebyshev-Jackson expansion of the Heaviside function of the Hamiltonian matrix is used to represent the occupation operator. This procedure requires only matrix-vector products and is intrinsically parallelizable.