  In this paper we propose an algorithm for the formation of matrices of
isogeometric Galerkin methods. The algorithm is based on three ideas. The first
is that we perform the external loop over the rows of the matrix. The second is
that we calculate the row entries by weighted quadrature. The third is that we
exploit the (local) tensor product structure of the basis functions. While all
ingredients have a fundamental role for computational efficiency, the major
conceptual change of paradigm with respect to the standard implementation is
the idea of using weighted quadrature: the test function is incorporated in the
integration weight while the trial function, the geometry parametrization and
the PDEs coefficients form the integrand function. This approach is very
effective in reducing the computational cost, while maintaining the optimal
order of approximation of the method. Analysis of the cost is confirmed by
numerical testing, where we show that, for $p$ large enough, the time required
by the floating point operations is less than the time spent in unavoidable
memory operations (the sparse matrix allocation and memory write). The proposed
algorithm allows significant time saving when assembling isogeometric Galerkin
matrices for all the degrees of the test spline space and paves the way for a
use of high-degree $k$-refinement in isogeometric analysis.
