Ab initio electronic structure calculations within Kohn-Sham density functional theory requires a solution for the Kohn-Sham equation. However, the traditional self-consistent field (SCF) approach of solving the equation using iterative diagonalization exhibits an inherent cubic scaling behavior and becomes prohibitive for large systems. The Chebyshev-filtered subspace iteration (CheFSI) method holds considerable promise for large-system calculations by substantially accelerating the SCF procedure. Here, we employed a combination of the real space finite-difference formulation and CheFSI to solve the Kohn-Sham equation, and implemented this approach in ab initio Real-space Electronic Structure (ARES) software in a multi-processor, parallel environment. An improved scheme was proposed to generate the initial subspace of Chebyshev filtering in ARES efficiently, making it suitable for large-scale simulations. The accuracy, stability, and efficiency of the ARES software were illustrated by simulations of large-scale crystalline systems containing thousands of atoms.