  Generalized-ensemble Monte Carlo simulations such as the multicanonical
method and similar techniques are among the most efficient approaches for
simulations of systems undergoing discontinuous phase transitions or with
rugged free- energy landscapes. As Markov chain methods, they are inherently
serial computationally. It was demonstrated recently, however, that a
combination of independent simulations that communicate weight updates at
variable intervals allows for the efficient utilization of parallel
computational resources for multicanonical simulations. Implementing this
approach for the many-thread architecture provided by current generations of
graphics processing units (GPUs), we show how it can be efficiently employed
with of the order of $10^4$ parallel walkers and beyond, thus constituting a
versatile tool for Monte Carlo simulations in the era of massively parallel
computing. We provide the fully documented source code for the approach applied
to the paradigmatic example of the two-dimensional Ising model as starting
point and reference for practitioners in the field.
