Solution (source code)

= Solution

Solving the $m$ unrelated shifted systems by separate <Gaussian elimination> costs $\boxed{\tfrac23m^4+O(m^3)}$ <floating-point operations>.

Instead compute $H=Q^TAQ$ once in $O(m^3)$ work. For each shift, solve
$$
(H-\lambda_jI)z_j=Q^Tb_j,\qquad x_j=Qz_j.
$$
The shift changes only the diagonal, so $H-\lambda_jI$ remains an <upper Hessenberg matrix>. In each elimination column only the entry immediately below the pivot needs removal. With partial pivoting the only possible nonzero candidate rows are the current row and the next, and each update touches $O(m)$ entries. Factorization therefore costs $O(m^2)$ per shift; the triangular solve and both orthogonal-vector transformations also cost $O(m^2)$ each.

Thus $m$ shifted solves plus the initial <Hessenberg reduction> cost $\boxed{O(m^3)}$ in total. Pivoting is still needed in general: nonsingularity excludes a singular system but does not by itself guarantee that every unpivoted elimination pivot is nonzero.