Diagrammatic Monte Carlo methods provide robust routines for accurate computations of correlated electronic systems in the thermodynamical limit. Recently, its versatility was extended to SU(N) Hubbard model, where the core is a novel dynamical programming approach to the summation of all connected Feynman diagrams. We present several generalizations of it with more interaction vertices and symmetry broken terms. The framework treats SU(N) symmetry breaking both from nonuniform, flavor-dependent chemical potentials and from spontaneously broken phases induced by shift parameters. We also provide an end-to-end GPU acceleration path with dedicated CUDA C++ optimizations independent of a previous CUDA acceleration approach, where the parallelization strategy is different. We performed detailed numerical study of new algorithms involved in this article and exposed numerical instabilities of connected determinant formalism, which we solved in multiple ways. Together, these advances establish a scalable, high-performance DiagMC toolbox for multi-flavor correlated systems with and without symmetry breaking.