  This paper proposes a combination of a hybrid CPU--GPU and a pure GPU
software implementation of a direct algorithm for solving shifted linear
systems $(A - \sigma I)X = B$ with large number of complex shifts $\sigma$ and
multiple right-hand sides. Such problems often appear e.g. in control theory
when evaluating the transfer function, or as a part of an algorithm performing
interpolatory model reduction, as well as when computing pseudospectra and
structured pseudospectra, or solving large linear systems of ordinary
differential equations. The proposed algorithm first jointly reduces the
general full $n\times n$ matrix $A$ and the $n\times m$ full right-hand side
matrix $B$ to the controller Hessenberg canonical form that facilitates
efficient solution: $A$ is transformed to a so-called $m$-Hessenberg form and
$B$ is made upper-triangular. This is implemented as blocked highly parallel
CPU--GPU hybrid algorithm; individual blocks are reduced by the CPU, and the
necessary updates of the rest of the matrix are split among the cores of the
CPU and the GPU. To enhance parallelization, the reduction and the updates are
overlapped. In the next phase, the reduced $m$-Hessenberg--triangular systems
are solved entirely on the GPU, with shifts divided into batches. The benefits
of such load distribution are demonstrated by numerical experiments. In
particular, we show that our proposed implementation provides an excellent
basis for efficient implementations of computational methods in systems and
control theory, from evaluation of transfer function to the interpolatory model
reduction.
