We study numerical methods and algorithms for time-dependent fractional-in-space diffusion problems. The considered anomalous diffusion is modelled by the fractional Laplacian (−Δ)α, 0<α<1, following the integral definition. Fractional diffusion is non-local, and the finite element method (FEM) discretization in space leads to a dense stiffness matrix. It is well known that numerically solving such non-local boundary value problems is expensive. Difficulties increase significantly when the problem is time-dependent. The aim of the article is to develop computationally efficient methods and algorithms. There are two main features of our approach. Hierarchical semi-separable (HSS) compression is applied for an approximate solution of the arising linear systems. For time discretization, we use an adaptive forward–backward Euler scheme. The properties of the composite algorithm thus obtained are investigated. In particular, the block representation of HSS compression allowed us to upgrade the HSS solver to efficiently handle varying diagonally perturbed transition matrices corresponding to changing time steps. The contribution of the paper is threefold. The methods are completely constructive, which allows for a clearly structured description of the algorithms. A theoretical estimate of the computational complexity is presented. It shows the advantages of the adaptive time stepping in combination with the HSS solver. Theoretical results are supported by representative numerical experiments. Both sequential and parallel scalability and efficiency are analyzed. The presented results provide convincing proof of the concept of the proposed methods and algorithms.