We consider an adaptive algorithm for finite element methods for the isogeometric analysis (IGAFEM) of elliptic (possibly non-symmetric) second-order partial differential equations in arbitrary space dimension d ≥ 2. We employ hierarchical B-splines of arbitrary degree and different order of smoothness. We propose a refinement strategy to generate a sequence of locally refined meshes and corresponding discrete solutions. Adaptivity is driven by some weighted residual a posteriori error estimator. We prove linear convergence of the error estimator (respectively, the sum of energy error plus data oscillations) with optimal algebraic rates. Numerical experiments underpin the theoretical findings.