Excited states exhibiting double-excitation character are notoriously difficult to model using conventional single-reference methods, such as adiabatic time-dependent density functional theory (TD-DFT) or equation-of-motion coupled cluster (EOM-CC). In addition, these states are typical experimentally "dark", making their detection in photoabsorption spectra very challenging. Nonetheless, they play a key role in the faithful description of many physical, chemical, and biological processes. In the present work, we provide accurate reference excitation energies for transitions involving a substantial amount of double excitation using a series of increasingly large diffuse-containing atomic basis sets. Our set gathers 20 vertical transitions from 14 small- and medium-size molecules (acrolein, benzene, beryllium atom, butadiene, carbon dimer and trimer, ethylene, formaldehyde, glyoxal, hexatriene, nitrosomethane, nitroxyl, pyrazine, and tetrazine). Depending on the size of the molecule, selected configuration interaction (sCI) and/or multiconfigurational (CASSCF, CASPT2, (X)MS-CASPT2, and NEVPT2) calculations are performed in order to obtain reliable estimates of the vertical transition energies. In addition, coupled cluster approaches including at least contributions from iterative triples (such as CC3, CCSDT, CCSDTQ, and CCSDTQP) are assessed. Our results clearly evidence that the error in CC methods is intimately related to the amount of double-excitation character of the transition. For "pure" double excitations (i.e., for transitions which do not mix with single excitations), the error in CC3 can easily reach 1 eV, while it goes down to a few tenths of an electronvolt for more common transitions (such as in trans-butadiene) involving a significant amount of singles. As expected, CC approaches including quadruples yield highly accurate results for any type of transition. The quality of the excitation energies obtained with multiconfigurational methods is harder to predict. We have found that the overall accuracy of these methods is highly dependent on both the system and the selected active space. The inclusion of the σ and σ* orbitals in the active space, even for transitions involving mostly π and π* orbitals, is mandatory in order to reach high accuracy. A theoretical best estimate (TBE) is reported for each transition. We believe that these reference data will be valuable for future methodological developments aiming at accurately describing double excitations.