We develop an inverse matrix method to solve for resonance masses from a dispersion relation obeyed by a correlation function. Given the operator product expansion (OPE) of a correlation function in the deep Euclidean region, we obtain the nonperturbative spectral density, which exhibits resonance structures naturally. The value of the gluon condensate in the OPE is fixed by producing the ρ meson mass in the formalism, and then input into the dispersion relations for the scalar, pseudoscalar and tensor glueballs. It is shown that the low-energy limit of the correlation function for the scalar glueball, derived from the spectral density, discriminates the lattice estimate for the triple-gluon condensate from the single-instanton estimate. The spectral densities for the scalar and pseudoscalar glueballs reveal a double-peak structure: the peak located at lower mass implies that the f0(500) and f0(980) (η ad η ) mesons contain small amount of gluonium components, and should be included into scalar (pseudoscalar) mixing frameworks. Another peak determines the scalar (pseudoscalar) glueball mass around 1.50 (1.75) GeV with a broad width about 200 MeV, suggesting that the f0(1370), f0(1500) and f0(1710) (η( 1760)) mesons are the glue-rich states. We also predict the topological susceptability χ 1/4 t = 75-78 MeV, deduced from the correlation function for the pseudoscalar glueball at zero momentum. Our analysis gives no resonance solution for the tensor glueball, which may be attributed to the insufficient nonperturbative condensate information in the currently available OPE.