  In this paper, a two-sided variable-coefficient space-fractional diffusion
equation with fractional Neumann boundary condition is considered. To conquer
the weak singularity caused by the nonlocal space-fractional differential
operators, by introducing an auxiliary fractional flux variable and using
piecewise linear interpolations, a fractional block-centered finite difference
(BCFD) method on general nonuniform grids is proposed. However, like other
numerical methods, the proposed method still produces linear algebraic systems
with unstructured dense coefficient matrices under the general nonuniform
grids.Consequently, traditional direct solvers such as Gaussian elimination
method shall require $\mathcal{O}(M^2)$ memory and $\mathcal{O}(M^3)$
computational work per time level, where $M$ is the number of spatial unknowns
in the numerical discretization. To address this issue, we combine the
well-known sum-of-exponentials (SOE) approximation technique with the
fractional BCFD method to propose a fast version fractional BCFD algorithm.
Based upon the Krylov subspace iterative methods, fast matrix-vector
multiplications of the resulting coefficient matrices with any vector are
developed, in which they can be implemented in only $\mathcal{O}(MN_{exp})$
operations per iteration, where $N_{exp}\ll M$ is the number of exponentials in
the SOE approximation. Moreover, the coefficient matrices do not necessarily
need to be generated explicitly, while they can be stored in
$\mathcal{O}(MN_{exp})$ memory by only storing some coefficient vectors.
Numerical experiments are provided to demonstrate the efficiency and accuracy
of the method.
