We present a GPU-implementation of linear response theory in the Kohn-Sham DFT approximation with an adoption of localized Gaussian atomic orbitals-it encompasses the generalized eigenvalue equation as well as the complex polarization propagator. In the iterative linear response equation solvers, we adopt the use of symmetric and anti-symmetric trial vectors as introduced by Kauczor et al. for the reduced space algorithms. The key compute intensive operation is the construction of (auxiliary) Fock matrices from density matrices that are either symmetric or anti-symmetric depending on the symmetry of the trial vectors. This Fock matrix construction, in turn, has two separate parts of significance to the computational cost namely the contraction or the ERI tensor with said density matrices and the numerical DFT kernel integration, both of which are carried out with use of GPUs in the present work.