Fractional diffusion equations model phenomena exhibiting anomalous diffusion that cannot be modeled accurately by second-order diffusion equations. Because of the nonlocal property of fractional differential operators, numerical methods for space-fractional diffusion equations generate complicated dense or full coefficient matrices. Consequently, these numerical methods were traditionally solved by Gaussian elimination, which requires computational work of O(N³) per time step and O(N²) of memory, where N is the number of spatial grid points in the discretization. The significant computational work and memory requirement of the numerical methods impose a serious challenge for the numerical simulation of two- and especially three-dimensional space-fractional diffusion equations. We develop a fast and yet accurate solution method for the implicit finite difference discretization of space-fractional diffusion equations in two 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(N log N) 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 show that the fast method has a significant reduction of CPU time, from two months and eight days as consumed by the traditional method to less than 40 minutes, with less than one ten-thousandth of the memory required by the traditional method, in the context of a two-dimensional space-fractional diffusion equation with 512 x 512 = 262,144 spatial nodes and 512 time steps. This demonstrates the utility of the method. [ABSTRACT FROM AUTHOR]