The safety assessment of offshore infrastructures, such as submarine cables and monopile foundations, requires a quantitative understanding of submarine landslides and their solid–fluid interaction mechanisms. To address this, an improved semi-implicit solid–fluid coupled double-point Material Point Method (DP-MPM) is proposed to capture collapse and fluidization during the transition from quasi-static deformation to a flowing soil–water mixture. Water-saturated soil is modeled within the framework of porous media theory, with water treated as an incompressible fluid. Incompressibility is enforced using a fractional-step projection method applied to the momentum equations, implemented in the time discretization scheme. This treatment effectively mitigates pressure oscillations associated with weak incompressibility constraints and allows for a substantial increase in the stable time step. To enhance numerical accuracy, second-order B-spline basis functions and the Affine Particle-In-Cell (APIC) mapping scheme are employed for both phases. In this study, an MPM algorithm is developed with a primary focus on numerical stability and computational efficiency. For that end, we implement an implicit drag term through an analytical solution of the resulting quadratic equation, thereby avoiding the use of a Newton–Raphson iteration. In addition, a Finite Element Method (FEM)–based Geometric Multigrid (MG) preconditioner is implemented for a parallel, matrix-free Preconditioned Conjugate Gradient (PCG) solver applied to the pressure Poisson equation. This approach eliminates the substantial memory overhead associated with explicit coefficient matrix assembly. Furthermore, to ensure a well-conditioned Poisson system, a particle-based signed-distance function is employed to capture the free surface. A ghost fluid method is then used to impose Dirichlet pressure boundary conditions consistently across the fluid–air interface. The proposed DP-MPM framework is implemented using the JAX framework to enable parallel computation on graphics processing unit (GPU) platforms. The proposed implementation is verified against analytical solutions as well as numerical simulations of submerged column collapse and submerged slope collapse.