  Fully implicit Runge-Kutta (IRK) methods have many desirable properties as
time integration schemes in terms of accuracy and stability, but high-order IRK
methods are not commonly used in practice with numerical PDEs due to the
difficulty of solving the stage equations. This paper introduces a theoretical
and algorithmic preconditioning framework for solving the systems of equations
that arise from IRK methods applied to linear numerical PDEs (without algebraic
constraints). This framework also naturally applies to discontinuous Galerkin
discretizations in time. Under quite general assumptions on the spatial
discretization that yield stable time integration, the preconditioned operator
is proven to have condition number bounded by a small, order-one constant,
independent of the spatial mesh and time-step size, and with only weak
dependence on number of stages/polynomial order; for example, the
preconditioned operator for 10th-order Gauss IRK has condition number less than
two, independent of the spatial discretization and time step. The new method
can be used with arbitrary existing preconditioners for backward Euler-type
time stepping schemes, and is amenable to the use of three-term recursion
Krylov methods when the underlying spatial discretization is symmetric. The new
method is demonstrated to be effective on various high-order finite-difference
and finite-element discretizations of linear parabolic and hyperbolic problems,
demonstrating fast, scalable solution of up to 10th order accuracy. The new
method consistently outperforms existing block preconditioning approaches, and
in several cases, the new method can achieve 4th-order accuracy using Gauss
integration with roughly half the number of preconditioner applications and
wallclock time as required using standard diagonally implicit RK methods.
