We extend the recently proposed heat-bath configuration interaction (HCI) method [Holmes, Tubman, Umrigar, J. Chem. Theory Comput. 2016, 12, 3674], by introducing a semistochastic algorithm for performing multireference Epstein-Nesbet perturbation theory, in order to completely eliminate the severe memory bottleneck of the original method. The proposed algorithm has several attractive features. First, there is no sign problem that plagues several quantum Monte Carlo methods. Second, instead of using Metropolis-Hastings sampling, we use the Alias method to directly sample determinants from the reference wave function, thus avoiding correlations between consecutive samples. Third, in addition to removing the memory bottleneck, semistochastic HCI (SHCI) is faster than the deterministic variant for many systems if a stochastic error of 0.1 mHa is acceptable. Fourth, within the SHCI algorithm one can trade memory for a modest increase in computer time. Fifth, the perturbative calculation is embarrassingly parallel. The SHCI algorithm extends the range of applicability of the original algorithm, allowing us to calculate the correlation energy of very large active spaces. We demonstrate this by performing calculations on several first row dimers including F2 with an active space of (14e, 108o), Mn-Salen cluster with an active space of (28e, 22o), and Cr2 dimer with up to a quadruple-ζ basis set with an active space of (12e, 190o). For these systems we were able to obtain better than 1 mHa accuracy with a wall time of merely 55 s, 37 s, and 56 min on 1, 1, and 4 nodes, respectively.