  We develop a robust matrix-free, communication avoiding parallel, high-degree
polynomial preconditioner for the Conjugate Gradient method for large and
sparse symmetric positive definite linear systems. We discuss the selection of
a scaling parameter aimed at avoiding unwanted clustering of eigenvalues of the
preconditioned matrices at the extrema of the spectrum. We use this
preconditioned framework to solve a $3 \times 3$ block system arising in the
simulation of fluid flow in large-size discrete fractured networks. We apply
our polynomial preconditioner to a suitable Schur complement related with this
system, which can not be explicitly computed because of its size and density.
Numerical results confirm the excellent properties of the proposed
preconditioner up to very high polynomial degrees. The parallel implementation
achieves satisfactory scalability by taking advantage from the reduced number
of scalar products and hence of global communications.
