Ising machines are domain-specific computers that solve combinatorial optimization problems (COPs). They utilize an Ising model to represent a COP and search for the optimal spin configuration of the Ising model to solve the COP. Most Ising machines are based on simulated annealing (SA) and update the spin configuration according to single-spin-flip Markov Chain Monte Carlo methods. However, supporting a multi-spin flip is important to enhance the search performance of SA-based Ising machines when constrained COPs are solved. In this paper, we extend the merge method, which was introduced to provide a multi-spin-flip capability for SA-based Ising machines, and formulate it as a series of matrix multiplications that can be executed on a graphics processing unit (GPU) efficiently. We then construct a GPU-based Ising machine that implements the GPU-oriented merge method together with an extended SA algorithm. We finally demonstrate its superiority over the state-of-the-art GPU-based Ising machine for quadratic knapsack problems.INDEX TERMS constrained combinatorial optimization, graphics processing unit (GPU), Ising machine, Ising model, multi-spin flip, simulated annealing (SA)