Fractional diffusion equations model phenomena exhibiting anomalous diffusion that cannot be modeled accurately by second-order diffusion equations. Because of the non-local property of fractional differential operators, numerical methods for space-fractional diffusion equations generate dense or even full coefficient matrices with complicated structures. Traditionally, these methods were solved with Gaussian elimination, which requires computational work of O(N3) per time step and O(N2) of memory to store where N is the number of spatial grid points in the discretization. The significant computational work and memory requirement of these methods makes a numerical simulation of three-dimensional space-fractional diffusion equations computationally prohibitively expensive.In this paper we develop an efficient and faithful solution method for the implicit finite difference discretization of time-dependent space-fractional diffusion equations in three space dimensions, by carefully analyzing the structure of the coefficient matrix of the finite difference method and delicately decomposing the coefficient matrix into a combination of sparse and structured dense matrices. The fast method has a computational work count of O(NlogN) per iteration and a memory requirement of O(N), while retaining the same accuracy as the underlying finite difference method solved with Gaussian elimination. Numerical experiments of a three-dimensional space-fractional diffusion equation show the utility of the fast method.