A scalable matrix solver was developed for the moving particle hydrodynamics for incompressible flows (MPH-I) method. Since the MPH-I method can calculate both incompressible and highly viscous flows while ensuring stability through physical consistency, a wide range of industrial applications is expected. However, in its implicit calculation, both the pressure and velocity must be solved simultaneously via a linear equation with a nondefinite symmetric coefficient matrix. In this study, this nondefinite linear system was converted into a symmetric positive definite (SPD) system where only the velocity is unknown. This conversion enabled us to solve the system with well-known solvers such as the conjugated gradient (CG) and conjugated residual (CR) methods. For scalability, bucket-based multigrid preconditioned CG and CR solvers were developed for the SPD system. To handle multidimensionality during preconditioning, an extended Jacobi smoother that is even applicable in a nondiagonally dominant matrix system was proposed. The numerical efficiency was confirmed via a simple high-viscosity incompressible dam break calculation, and the scalability within the presented case was confirmed. In addition, the performance under shared memory parallel computations was studied.