We consider the stiffness matrices arising from the Galerkin B-spline isogeometric analysis discretization of classical elliptic problems. By exploiting their specific spectral properties, compactly described by a symbol, we design an efficient multigrid method for the fast solution of the related linear systems. The convergence rate of general-purpose multigrid methods, based on classical stationary smoothers, is optimal (i.e., bounded independently of the matrix size), but it also worsens exponentially with respect to the spline degree. The symbol allows us to give a detailed theoretical explanation of this exponential worsening in the case of the two-grid scheme. In addition, thanks to a specific factorization of the symbol, we are able to design an ad hoc multigrid method with an effective preconditioned CG or GMRES smoother at the finest level, in the spirit of the multi-iterative idea. The convergence rate of this multi-iterative multigrid method is not only optimal but also robust (i.e., bounded su...