  Taking partial traces for computing reduced density matrices, or related
functions, is a ubiquitous procedure in the quantum mechanics of composite
systems. In this article, we present a thorough description of this function
and analyze the number of elementary operations (ops) needed, under some
possible alternative implementations, to compute it on a classical computer. As
we notice, it is worthwhile doing some analytical developments in order to
avoid making null multiplications and sums, what can considerably reduce the
ops. For instance, for a bipartite system
$\mathcal{H}_{a}\otimes\mathcal{H}_{b}$ with dimensions
$d_{a}=\dim\mathcal{H}_{a}$ and $d_{b}=\dim\mathcal{H}_{b}$ and for
$d_{a},d_{b}\gg1$, while a direct use of partial trace definition applied to
$\mathcal{H}_{b}$ requires $\mathcal{O}(d_{a}^{6}d_{b}^{6})$ ops, its optimized
implementation entails $\mathcal{O}(d_{a}^{2}d_{b})$ ops. In the sequence, we
regard the computation of partial traces for general multipartite systems and
describe Fortran code provided to implement it numerically. We also consider
the calculation of reduced density matrices via Bloch's parametrization with
generalized Gell Mann's matrices.
