In this talk we discuss the analysis and fast implementation of numerical algorithms for the solution of fractional linear elliptic problems on a bounded domain with a homogeneous Dirichlet boundary condition. We focus on the spectral fractional Laplacian and the corresponding Caffarelli-Silvestre extension. We make use of tensor product finite element spaces with hp finite elements in the extended direction and piecewise linear or hp finite elements in the bounded domain and prove convergence estimates. We describe fast methods for solving the resulting linear system resulting in algorithms of almost optimal complexity. Numerical results will be shown that confirm the theoretical predictions.