In this paper, we explore the problem of learning the parameters of Gaussian mixture models using the expectation-maximization (EM) algorithm. We propose a new parallel formulation of the EM algorithm that utilizes a static decomposition and distributes the learning set and the matrix storing posterior probabilities among MPI processes. During both the E-step and M-step of an EM iteration, each MPI process spawns a team of OpenMP threads that process data in blocks sized to fit the last-level cache. The calculations necessary for obtaining weighted Gaussian densities in the E-step and mixture parameters in the M-step are conducted using optimized level 3 BLAS operations. Hierarchical all-reduce operations are employed to compute the sums required in the M-step, initially among the OpenMP threads within each MPI process and subsequently among all the MPI processes in the parallel application. In our computational experiments, we compared this proposed approach with the traditional method that employs level 2 BLAS operations and assessed its strong scaling on 64 nodes of a compute cluster. The results demonstrate that the proposed method is 1.6 to 4.7 times faster than the conventional approach utilizing level 2 BLAS. The parallel efficiency of our approach on 64 nodes ranges from 60% to 73%.