Extrapolation cascadic multigrid (EXCMG) method with the conjugate gradient smoother is shown to be an efficient solver for large sparse symmetric positive definite systems resulting from linear finite element discretization of second-order elliptic boundary value problems [Pan et al. J. Comput. Phys. 344 (2017) 499–515]. In this paper, we generalize the EXCMG method to solve a class of spatial fractional diffusion equations (SFDEs) with variable coefficients. Both steady-state and time-dependent problems are considered. First of all, space-fractional derivatives defined in Riemann–Liouville sense are discretized by using the weighted average of shifted Grünwald formula, which results in a fully dense and nonsymmetric linear system for the steady-state problem, or a semi-discretized ODE system for the time-dependent problem. Then, to solve the former problem, we propose the EXCMG method with the biconjugate gradient stabilized smoother to deal with the dense and nonsymmetric linear system. Next, such technique is extended to solve the latter problem since it becomes fully discrete when the Crank-Nicolson scheme is introduced to handle the temporal derivative. Finally, several numerical examples are reported to show that the EXCMG method is an efficient solver for both steady-state and time-dependent SFDEs, and performs much better than the V-cycle multigrid method with banded-splitting smoother for time-dependent SFDEs [Lin et al. J. Comput. Phys. 336 (2017) 69–86].