An efficient method for computing thermodynamic equilibrium states at the micromagnetic length scale is introduced, based on the Markov chain Monte Carlo method.Trial moves include not only rotations of spins, but also a change in their magnetization length.The method is parameterized using the longitudinal susceptibility, reproduces the same Maxwell-Boltzmann distribution as the stochastic Landau-Lifshitz-Bloch equation, and is applicable both below and above the Curie temperature. The algorithm is fully parallel, can be executed on graphical processing units, and efficiently includes the long range dipolar interaction. This method is generally useful for computing finite-temperature relaxation states both for uniform and non-uniform temperature profiles, and can be considered as complementary to zero-temperature micromagnetic energy minimization solvers, with comparable computation time. However, unlike quasi-zero temperature approaches which do not take into account the magnetization length distribution and stochasticity, the method is better suited for structures with unbroken symmetry around the applied field axis, granular films, and at higher temperatures and fields. In particular, applications to finite-temperature hysteresis loop modelling, chiral magnetic thin films, granular magnetic media, and artificial spin ices are discussed.