We present a method for scalable and fully 3D magnetic field simultaneouslocalisation and mapping (SLAM) using local anomalies in the magnetic field asa source of position information. These anomalies are due to the presence offerromagnetic material in the structure of buildings and in objects such asfurniture. We represent the magnetic field map using a Gaussian process modeland take well-known physical properties of the magnetic field into account. Webuild local magnetic field maps using three-dimensional hexagonal block tiling.To make our approach computationally tractable we use reduced-rank Gaussianprocess regression in combination with a Rao--Blackwellised particle filter. Weshow that it is possible to obtain accurate position and orientation estimatesusing measurements from a smartphone, and that our approach provides a scalablemagnetic SLAM algorithm in terms of both computational complexity and mapstorage.