  The real-space density-functional perturbation theory (DFPT) for the
computations of the response properties with respect to the atomic displacement
and homogeneous electric field perturbation has been recently developed and
implemented into the all-electron, numeric atom-centered orbitals electronic
structure package FHI-aims. It is found that the bottleneck for large scale
applications is the computation of the response density matrix, which scales as
$O(N^3)$. Here for the response properties with respect to the homogeneous
electric field, we present an efficient parallel linear scaling algorithm for
the response density matrix calculation. Our scheme is based on the
second-order trace-correcting purification and the parallel sparse
matrix-matrix multiplication algorithms. The new scheme reduces the formal
scaling from $O(N^3)$ to $O(N)$, and shows good parallel scalability over tens
of thousands of cores. As demonstrated by extensive validation, we achieve a
rapid computation of accurate polarizabilities using DFPT. Finally, the
computational efficiency of this scheme has been illustrated by making the
scaling tests and scalability tests on massively parallel computer systems.
