Geometrically thick tori with constant specific angular momentum have been widely used in the last decades to construct numerical models of accretion flows onto black holes. Such discs are prone to a global non-axisymmetric hydrodynamical instability, known as Papaloizou-Pringle instability (PPI), which can redistribute angular momentum and also lead to an emission of gravitational waves. It is, however, not clear yet how the development of the PPI is affected by the presence of a magnetic field and by the concurrent development of the magnetorotational instability (MRI). We present the first non-linear analysis using three-dimensional GRMHD simulations of the interplay between the PPI and the MRI and compare results between simulations of non-magnetized and magnetized tori. In the purely hydrodynamic case, the PPI selects the $m=1$ azimuthal mode as the fastest growing and non-linearly dominant mode. We show that, for tori with a weak toroidal magnetic field, the development of the MRI leads to the suppression of large-scale modes and redistributes power across smaller scales. For strong enough (but still subthermal) magnetic fields the $m=1$ mode is no longer dominant, and the PPI appears to be completely inhibited.