Understanding the interaction of vortices with inclusions in type-II superconductors is a major outstanding challenge both for fundamental science and energy applications. At application-relevant scales, the long-range interactions between a dense configuration of vortices and the dependence of their behavior on external parameters, such as temperature and an applied magnetic field, are all important to the net response of the superconductor. Capturing these features, in general, precludes analytical description of vortex dynamics and has also made numerical simulation prohibitively expensive. Here we report on a highly optimized iterative implicit solver for the time-dependent Ginzburg-Landau equations suitable for investigations of type-II superconductors on massively parallel architectures. Its main purpose is to study vortex dynamics in disordered or geometrically confined mesoscopic systems. In this work, we present the discretization and time integration scheme in detail for two types of boundary conditions. We describe the necessary conditions for a stable and physically accurate integration of the equations of motion. Using an inclusion pattern generator, we can simulate complex pinning landscapes and the effect of geometric confinement. We show that our algorithm, implemented on a GPU, can provide static and dynamic solutions of the Ginzburg-Landau equations for mesoscopically large systems over thousands of time steps in a matter of hours. Using our formulation, studying scientifically-relevant problems is a computationally reasonable task.