Recently, fractional calculations have gained extensive attention within computational mathematics and applied science. Their numerical treatment is an important research area as such equations pose substantial challenges to existing algorithms. In the first part of this thesis, we review some basic concepts of the fractional derivatives. The computational challenge of the fractional partial differential equations is also discussed. The main research topic of this thesis is the fast solvers for fractional partial differential equations. In the second part, we construct an approximate inverse preconditioner for the discretized linear systems arising from the space-fractional diffusion equations with piecewise continuous coefficients, where the interpolation formula is applied to approximate the eigenvalues of circulant matrices. Therefore, the discontinuity of the diffusion coefficients does not influence the efficiency of the preconditioner. Theoretically, the spectra of the resulting preconditioned matrices are shown to be clustered around one. Numerical examples are provided to demonstrate the effectiveness of our method. The third part concerns circulant preconditioners for discretized matrices arising from finite difference schemes for a kind of fractional diffusion equations. The fractional differential operator is comprised of left-sided and right-sided derivatives with order in (½, 1). The resulting discretized matrices preserve Toeplitz-like structure and hence their matrix-vector multiplications can be computed efficiently by the fast Fourier transform, Using the circulant preconditioning technique, we show that the spectra of the circulant preconditioned matrices are clustered around one under some conditions. Numerical experiments are presented to demonstrate that the preconditioning technique is very efficient. In the fourth part, we propose a fast algorithm for the differentiation of the Caputo fractional derivative and apply it to numerically solve the fractional sub-diffusion equations. The storage requirements and computational workload are of Ο(m log²n)and Ο(mn log²n), respectively, where m is the total number of the spatial grid and n is the size of the temporal grid. The resulting finite difference scheme is theoretically proved to be stable. Numerical examples are presented to show the efficiency of our fast algorithm.