Prestack depth migration for seismic reflection data is commonly used tool for imaging complex geological structures such as salt domes, faults, thrust belts, and stratigraphic structures. Phase shift plus interpolation (PSPI) algorithm is a useful tool to directly solve a wave equation and the results have natural properties of the wave equation. Amplitude and phase characteristics, in particular, are better preserved. The PSPI algorithm is widely used in hydrocarbon exploration because of its simplicity, efficiency, and reduced efforts for computation. However, meaningful depth image of 3D subsurface requires parallel computing to handle heavy computing time and great amount of input data. We implemented a parallelized version of 3D PSPI for prestack depth migration using Open-Multi-Processing (Open MP) library. We verified its performance through applications to 3D SEG/EAGE salt model with a small scale Linux cluster. Phase-shift was performed in the vertical and horizontal directions, respectively, and then interpolated at each node. This gave a single image gather according to shot gather. After summation of each single image gather, we got a 3D stacked image in the depth domain. The numerical model example shows good agreement with the original geological model.