In cases when the nature of the system and electronic transitions are such that a largescale QM region is called for, then one, in addition, typically faces the challenge of high densities of states, even if the gap energy itself is sizeable. This is computationally prohibitive for the bottom-up approaches used to solve the generalized eigenvalue equation in timedependent DFT. One way to circumvent this issue is to adopt the complex polarization propagator (CPP) formulation of response theory where the account made for the finite lifetimes of the excited states makes response functions physically sound also in resonance regions of the spectrum with real and imaginary parts being related by Kramers-Kronig relations. The underlying excited states are not resolved in the calculation which instead focuses on determining the response functions pertinent to the spectroscopy at hand for a discretized frequency window enclosing the sought for spectral bands. With a solution of the associated set of response equations made in parallel, the convergence can be highly accelerated. The CPP theory has been implemented in several quantum chemical programs and thereby for most standard electronic theory methods. In this work, we will employ the DFT implementation of real and complex linear response theory made in the VeloxChem software and introduce a GPU-accelerated construction of the regular Fock (or Kohn-Sham) matrix needed for the self-consistent field (SCF) optimization of the reference state as well as the auxiliary Fock matrices needed in the response part of the calculation.