  We develop a fast solver for the spectral element method (SEM) applied to the
two-sided fractional diffusion equation on uniform, geometric and graded
meshes. By approximating the singular kernel with a degenerate kernel, we
construct a hierarchical matrix (H-matrix) to represent the stiffness matrix of
the SEM and provide error estimates verified numerically. We can solve
efficiently the H-matrix approximation problem using a hierarchical LU
decomposition method, which reduces the computational cost to $O(R^2 N_d
\log^2N) +O(R^3 N_d \log N)$, where $R$ it is the rank of submatrices of the
H-matrix approximation, $N_d$ is the total number of degrees of freedom and $N$
is the number of elements. However, we lose the high accuracy of the SEM. Thus,
we solve the corresponding preconditioned system by using the H-matrix
approximation problem as a preconditioner, recovering the high order accuracy
of the SEM. The condition number of the preconditioned system is independent of
the polynomial degree $P$ and grows with the number of elements, but at modest
values of the rank $R$ is below order 10 in our experiments, which represents a
reduction of more than 11 orders of magnitude from the unpreconditioned system;
this reduction is higher in the two-sided fractional derivative compared to
one-sided fractional derivative. The corresponding cost is $O(R^2 N_d \log^2
N)+O(R^3 N_d \log N)+O(N_d^2)$. Moreover, by using a structured mesh (uniform
or geometric mesh), we can further reduce the computational cost to $O(R^2
N_d\log^2 N) +O(R^3 N_d \log N)+ O(P^2 N\log N)$ for the preconditioned system.
We present several numerical tests to illustrate the proposed algorithm using
$h$ and $p$ refinements.
