This paper focuses on designing and implementing parallel Moving Least Squares (MLS) interpolation algorithms by exploiting the Graphics Processing Unit (GPU) for the usage in Meshfree methods. The MLS method is an approach for scattered points' approximation / interpolation, which is commonly employed as the shape functions in various Meshfree methods. To improve the computational efficiency in building stiffness matrices in Meshfree methods, we are specifically interested in parallelizing the MLS interpolation on the GPU. The accelerated Meshfree methods can be employed to numerically analyze large deformations of soil and rock masses such as the landslides. In this paper, we first introduce our previously proposed method for finding the k nearest neighboring points located within the local region of the interest point in MLS. We then develop one sequential and three parallel implementations for each of three variations of the MLS interpolation, including the serial implementation, the parallel implementation on the multi-core CPU, the parallel implementation on a single GPU, and the parallel implementation on multi-GPUs.To evaluate the computational performance of our GPU implementations, five groups of benchmark tests are conducted in both two-dimensions and three-dimensions. We observe that our GPU implementations can achieve satisfied speedups over the sequential CPU implementation for varied sizes of testing data. To benefit the community, all source code and testing data related to the presented parallel MLS interpolation are publicly available.