In this paper, we present an efficient \GQR\ algorithm for solving the linear response eigenvalue problem $\scrH\bx ={\lambda}\bx $, where $\scrH$ is $\bPi^-$-symmetric with respect to $\Gamma_0 = \diag(I_n,-I_n)$. Based on newly introduced $\Gamma$-orthogonal transformations, the \GQR\ algorithm preserves the $\bPi^-$-symmetric structure of $\scrH$ throughout the whole process, which guarantees the computed eigenvalues to appear pairwise $(\lambda,-\lambda)$ as they should.
With the help of a newly established implicit $\Gamma$-orthogonality theorem, we incorporate the implicit multi-shift technique to accelerate the convergence of the \GQR\ algorithm. Numerical experiments are given to show the effectiveness of
the algorithm.