Kohn-Sham density functional theory (KSDFT) is the most widely used electronic structure theory for molecules and condensed matter systems. For a system with N electrons, the standard method for solving KSDFT requires solving N eigenvectors for an O(N) * O(N) Kohn-Sham Hamiltonian matrix. The computational cost for such procedure is expensive and scales as O(N^3). We have developed pole expansion plus selected inversion (PEXSI) method, in which KSDFT is solved by evaluating the selected elements of the inverse of a series of sparse symmetric matrices, and the overall algorithm scales at most O(N^2) for all materials including insulators, semiconductors and metals. The PEXSI method can be used with orthogonal or nonorthogonal basis set, and the physical quantities including electron density, energy, atomic force, density of states, and local density of states are calculated accurately without using the eigenvalues and eigenvectors. The recently developed massively parallel PEXSI method has been implemented in SIESTA, one of the most popular electronic structure software packages using atomic orbital basis sets. The resulting method can allow accurate treatment of electronic structure in a unprecedented scale. We demonstrate the application of the method for solving graphene-like structures with more than 30,000 atoms, and the method can be efficiently parallelized 10,000 - 100,000 processors on Department of Energy (DOE) high performance machines.