Bilinear models allow to characterize more precisely the behavior of physical processes in electric power systems and thereby improve the quality of their monitoring and control. However, the use of even medium size bilinear models usually requires a lot of computations. In order to improve the efficiency of bilinear analysis of power systems, this paper proposes an accelerated algorithm for calculating the Gramians of bilinear dynamical systems, which also takes into account the effects of sparse matrices. The conventional algorithms iterates the whole matrix. In contrast, in the proposed algorithm we continue iterate only large enough components of the solution matrix and exclude from iterations those components of the Gramian for which corrections become small. The algorithm was implemented and tested for four bilinear models.