In this talk, we construct a fully discrete scheme of the linear FE method in both temporal and spatial directions, derive many characterizations on the coefficient matrix and numerically verify that the fully FE approximation possesses the saturation error order under L2 norm.  We present an estimation like $1+\mathcal{O}(\tau^\alpha h^{-2\beta})$ on the condition number of the coefficient matrix. Finally, we develop and analyze an adaptive algebraic multigrid (AMG) method with low algorithmic complexity. We reveal a reference formula to measure the strength-of-connection tolerance which severely affect the robustness of AMG methods in handling fractional diffusion equations, and illustrate the well robustness and high efficiency of the proposed algorithm compared with the classical AMG, conjugate gradient and Jacobi iterative methods.