Accurate calculations of shortwave reflectances in clear-sky aerosol-laden atmospheres are necessary for various applications in atmospheric sciences. However computational cost becomes increasingly important for some applications such as data assimilation of top-of-atmosphere reflectances in models of atmospheric composition. This study aims to provide a benchmark that can help in assessing these two requirements in combination. We describe a protocol and input data for 44,080 cases involving various solar and viewing geometries, four different surfaces (one oceanic bidirectional reflectance function and three albedo values for a Lambertian surface), eight aerosol optical depths, five wavelengths, and four aerosol types. We first consider two models relying on the discrete ordinate method: VLIDORT (in vector and scalar configurations) and DISORT (scalar configuration only). We use VLIDORT in its vector configuration as a reference model and quantify the loss of accuracy due to (i) neglecting the effect of polarisation in the DISORT and VLIDORT (scalar) models and (ii) decreasing the number of streams in DISORT. We further test two other models: the 6SV2 model relying on the successive orders of scattering method and FLOTSAM, a new model under development by two of the authors. Typical mean fractional errors of 2.8 and 2.4 % for 6SV2 and FLOTSAM are found, respectively. Computational cost depends on the input parameters but also on the code implementation and application as some models solve the radiative transfer equations for a range of geometries while others do not. All necessary input and output data are provided as Supplementary Material as a potential resource for interested developers and users of radiative transfer models.