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.