Modern industrial computational fluid dynamics tools employing either implicit time steppers or linear harmonic methods have a strong requirement for the efficient iterative solution of large and sparse linear systems of equations, often relying on Krylov subspace methods, necessitating the best possible preconditioning approach to minimise the computational burden. In the context of aircraft aeroelasticity, the aerodynamic response to structural vibrations is commonly computed in the linearised frequency domain (LFD) [1]. Each structural mode (typically hundreds) and each frequency (tens) contribute to the generation of generalised aerodynamic forces useful for flutter analysis, ultimately requiring the solution of thousands of linear systems. The LFD incarnation of the flow solver TAU, developed by the German Aerospace Centre and used herein, usually relies on a solver stack that employs a Krylov method, such as generalized minimal residual, combined with incomplete lower-upper (ILU) factorisation of a suitably approximated coefficient matrix, with distributed memory parallelisation enabled through message passing interface [1]. As part of this work, we have implemented and assessed Jacobi-type iterations for preconditioning, that iteratively apply an ILU factorisation as the generalised diagonal. To obtain performance gains in the linear solution, we investigate a wide range of functionality, including parallel partitioning approaches, bandwidth optimisation (representative matrix sparsity patterns for natural and reverse Cuthill-McKee (RCM) orderings are shown in Figure 1), and choice of matrix approximation. The main test case is the NASA Common Research Model [2] with more than 100 million degrees of freedom. Scalability studies are undertaken across different levels of MPI parallelisation. For the cases investigated, we show improved performance and robustness when using the Krylov-Jacobi-ILU solver stack.