Rotation matrices are still ubiquitous even in systems which use quaternions: if you want to apply a dense linear transformation to a vector (be it a rotation or otherwise), pretty much nothing will beat (in terms of speed) a matrix-vector multiplication. For instance even if you work out the rotation matrices using quaternions, an engine will end up shipping a rotation matrix to a graphics card to actually transform the bulk of vertices in a scene.