The survey is devoted to numerical solution of the equation A(alpha)u = f, 0 < alpha < 1, where A is a symmetric positive definite operator corresponding to a second order elliptic boundary value problem in a bounded domain Omega in R-d. The fractional power A(alpha) is a non-local operator and is defined though the spectrum of A. Due to growing interest and demand in applications of sub-diffusion models to physics and engineering, in the last decade, several numerical approaches have been proposed, studied, and tested. We consider discretizations of the elliptic operator A by using an N-dimensional finite element space V-h or finite differences over a uniform mesh with N points. In the case of finite element approximation we get a symmetric and positive definite operator A(h) : V-h -> V-h, which results in an operator equation A(h)(alpha)u(h) = f(h) for u(h) is an element of V-h. The numerical solution of this equation is based on the following three equivalent representations of the solution: (1) Dunford-Taylor integral formula (or its equivalent Balakrishnan formula, (2.5)), (2) extension of the a second order elliptic problem in Omega x (0, infinity) subset of Rd+1 [17, 55] (with a local operator) or as a pseudo-parabolic equation in the cylinder (x, t) is an element of Omega x (0, 1), [70, 29], (3) spectral representation (2.6) and the best uniform rational approximation (BURA) of z(alpha) on [0,1], [37, 40]. Though substantially different in origin and their analysis, these methods can be interpreted as some rational approximation of A(h)(-alpha). In this paper we present the main ideas of these methods and the corresponding algorithms, discuss their accuracy, computational complexity and compare their efficiency and robustness.