Magnetohydrodynamic turbulence is important in many high-energy astrophysical systems, where instabilities can amplify the local magnetic field over very short timescales(1,2). Specifically, the magnetorotational instability and dynamo action(3-6) have been suggested as a mechanism for the growth of magnetar-strength magnetic fields (of 10(15) gauss and above) and for powering the explosion(7-10) of a rotating massive star(11,12). Such stars are candidate progenitors of type Ic-bl hypernovae(13,14), which make up all supernovae that are connected to long gamma-ray bursts(15,16). The magnetorotational instability has been studied with local high-resolution shearing-box simulations in three dimensions(17-19), and with global two-dimensional simulations(20), but it is not known whether turbulence driven by this instability can result in the creation of a large-scale, ordered and dynamically relevant field. Here we report results from global, three-dimensional, general-relativistic magnetohydrodynamic turbulence simulations. We show that hydromagnetic turbulence in rapidly rotating protoneutron stars produces an inverse cascade of energy. We find a large-scale, ordered toroidal field that is consistent with the formation of bipolar magnetorotationally driven outflows. Our results demonstrate that rapidly rotating massive stars are plausible progenitors for both type Ic-bl supernovae(13,21,22) and long gamma-ray bursts, and provide a viable mechanism for the formation of magnetars(23,24). Moreover, our findings suggest that rapidly rotating massive stars might lie behind potentially magnetar-powered superluminous supernovae(25,26).