  We present our preliminary results for the computation of the
non-perturbative running of renormalized quark masses in $N_f = 3$ QCD, between
the electroweak and hadronic scales, using standard finite-size scaling
techniques. The computation is carried out to very high precision, using
massless $\mathcal{O}(a)$-improved Wilson quarks. Following the strategy
adopted by the ALPHA Collaboration for the running coupling, different schemes
are used above and below a scale $\mu_0 \sim m_b$, which differ by using either
the Schr\"odinger Functional or Gradient Flow renormalized coupling. We discuss
our results for the running in both regions, and the procedure to match the two
schemes.
