With tens of petaflops supercomputers already in operation and exaflops machines expected to appear within the next 10 years, efficient parallel computational methods are required to take advantage of such extreme-scale machines. In this paper, we present a three-dimensional domain decomposition scheme for enabling large-scale electronic calculations based on density functional theory (DFT) on massively parallel computers. It is composed of two methods: (i) atom decomposition method and (ii) grid decomposition method. In the former, we develop a modified recursive bisection method based on inertia tensor moment to reorder the atoms along a principal axis so that atoms that are close in real space are also close on the axis to ensure data locality. The atoms are then divided into sub-domains depending on their projections onto the principal axis in a balanced way among the processes. In the latter, we define four data structures for the partitioning of grids that are carefully constructed to make data locality consistent with that of the clustered atoms for minimizing data communications between the processes. We also propose a decomposition method for solving the Poisson equation using three-dimensional FFT in Hartree potential calculation, which is shown to be better than a previously proposed parallelization method based on a two-dimensional decomposition in terms of communication efficiency. For evaluation, we perform benchmark calculations with our open-source DFT code, OpenMX, paying particular attention to the O(N) Krylov subspace method. The results show that our scheme exhibits good strong and weak scaling properties, with the parallel efficiency at 131,072 cores being 67.7% compared to the baseline of 16,384 cores with 131,072 diamond atoms on the K computer.