Vesicles are locally-inextensible fluid membranes that can sustain bending. We consider the dynamics of flows of vesicles suspended in Stokesian fluids. We use a boundary integral formulation for the fluid that results in a set of nonlinear integro-differential equations for the vesicle dynamics. The motion of the vesicles is determined by balancing the nonlocal hydrodynamic forces with the elastic forces due to bending and tension. Numerical simulations of such vesicle motions are quite challenging. On one hand, explicit time-stepping schemes suffer from a severe stability constraint due to the stiffness related to high-order spatial derivatives and a milder constraint due to a transport-like stability condition. On the other hand, an implicit scheme can be expensive because it requires the solution of a set of nonlinear equations at each time step. We present a set of numerical techniques for efficient simulation of vesicle flows. The distinctive features of these numerical methods include (1) using boundary integral method accelerated with the fast multipole method (2) spectral (spherical harmonic) discretization of deforming surfaces in space (3) an algorithm for surface reparameterization ensuring stability of the time-stepping scheme and spectral filtering accuracy while minimizing computational costs. By introducing these algorithmic components, we obtain a time-stepping scheme that experimentally is unconditionally stable and has a low cost per time step. We present numerical results to analyze the cost and convergence rates of the scheme.