We present a complete linear scaling method for hybrid Kohn-Sham density functional theory electronic structure calculations and demonstrate its performance. Particular attention is given to the linear scaling computation of the Kohn-Sham exchange-correlation matrix directly in sparse form within the generalized gradient approximation. The described method makes efficient use of sparse data structures at all times and scales linearly with respect to both computational time and memory usage. Benchmark calculations at the BHandHLYP/3-21G level of theory are presented for polypeptide helix molecules with up to 53 250 atoms. Threshold values for computational approximations were chosen on the basis of their impact on the occupied subspace so that the different parts of the calculations were carried out at balanced levels of accuracy. The largest calculation used 307 204 Gaussian basis functions on a single computer with 72 GB of memory. Benchmarks for three-dimensional water clusters are also included, as well as results using the 6-31G** basis set.