  The field of network synchronization has seen tremendous growth following the
introduction of the master stability function (MSF) formalism, which enables
the efficient stability analysis of synchronization in large oscillator
networks. However, to make further progress we must overcome the limitations of
this celebrated formalism, which focuses on global synchronization and requires
both the oscillators and their interaction functions to be identical, while
many systems of interest are inherently heterogeneous and exhibit complex
synchronization patterns. Here, we establish a generalization of the MSF
formalism that can characterize the stability of any cluster synchronization
pattern, even when the oscillators and/or their interaction functions are
nonidentical. The new framework is based on finding the finest simultaneous
block diagonalization of matrices in the variational equation and does not rely
on information about network symmetry. This leads to an algorithm that is
error-tolerant and orders of magnitude faster than existing symmetry-based
algorithms. As an application, we rigorously characterize the stability of
chimera states in networks with multiple types of interactions.
