Recent graphic processing units (GPUs) have remarkable raw computing power, which can be used for very computationally challenging problems. Like in micromagnetic simulations, where the magnetostatic field computation to analyze the magnetic behavior at very small time and space scale demands a huge computation time. This paper presents a multidimensional FFT-based parallel implementation of a magnetostatic field computation on GPUs. We have developed a specialized 3D FFT library for magnetostatic field calculation on GPUs. This made it possible to fully exploit the symmetries inherent in the field calculation and other optimizations specific to the GPUs architecture. We have compared our results with the widely used CPU-based parallel OOMMF program and with an equivalent serial implementation on CPU. The results have shown a speedup of up to 95x and 8.7x for single and 66x and 4.6x for double precision floating point accuracy against equivalent serial implementation and OOMMF, respectively.