A search for Higgs boson pair production in the bb¯ bb¯ final state is carried out with up to 36.1 fb−1 of LHC proton-proton collision data collected at s=13 TeV with the ATLAS detector in 2015 and 2016. Three benchmark signals are studied: a spin-2 graviton decaying into a Higgs boson pair, a scalar resonance decaying into a Higgs boson pair, and Standard Model non-resonant Higgs boson pair production. Two analyses are carried out, each implementing a particular technique for the event reconstruction that targets Higgs bosons reconstructed as pairs of jets or single boosted jets. The resonance mass range covered is 260–3000 GeV. The analyses are statistically combined and upper limits on the production cross section of Higgs boson pairs times branching ratio to bb¯ bb¯ are set in each model. No significant excess is observed; the largest deviation of data over prediction is found at a mass of 280 GeV, corresponding to 2.3 standard deviations globally. The observed 95% confidence level upper limit on the non-resonant production is 13 times the Standard Model prediction.[Figure not available: see fulltext.].