A search for heavy resonances decaying into pairs of standard model Higgs bosons is performed using data from proton-proton collisions at a centre-of-mass energy of 13 TeV collected by the CMS experiment in 2016, and corresponding to an integrated luminosity of 35.9 fb. The final state consists of both Higgs bosons decaying to b quark-antiquark pairs. For resonance masses above 1 TeV the Higgs bosons are Lorentz-boosted and each is reconstructed as one hadronic jet. The signal is characterized as a peak over the invariant mass spectrum of dijet events from standard model multijet processes. The signal strengths for different assumptions of resonance masses are estimated by a combined likelihood fit of the background and the signal shapes to the data. The results are consistent with the standard model expectations, and are interpreted as upper limits on the -channel production cross sections of narrow bulk gravitons and scalar radions in warped extra-dimensional models for resonance masses between 800 and 3000 GeV.