Within the NRQCD factorization approach the production of heavy quarkonia c
an be calculated in perturbative QCD in terms of a few production-process-i
ndependent NRQCD matrix elements. We test the universality of these long-di
stance matrix elements by comparing their values as determined from charmon
ium production at p (p) over bar, ep, and e(+)e(-) colliders, and in Z(0),
Y, and B decays.