Microscopic stochastic Langevin-type spin dynamics equations provide a convenient and tractable model describing the relaxation of spin and spin-lattice ensembles. We develop a robust and numerically stable algorithm for integrating the Langevin spin dynamics equations, and explore, both numerically and analytically, a range of applications of the method. We show that the algorithm conserves the magnitude of the spin vector irrespectively of the amplitude of the thermal noise. Using the Furutsu-Novikov theorem, we derive a system of deterministic differential equations for the ensemble-average moments of solutions of the Langevin spin dynamics equations, and explore the dynamics of relaxation of a spin ensemble toward the equilibrium Gibbs distribution. Analytical solutions of the moments equations make it possible to estimate the time scales of spin thermalization and spin-spin self-correlation, which we investigate as functions of the damping parameter and temperature.