  A preconditioning strategy is proposed for the iterative solve of large
numbers of linear systems with variable matrix and right-hand side which arise
during the computation of solution statistics of stochastic elliptic partial
differential equations with random variable coefficients sampled by Monte
Carlo. Building on the assumption that a truncated Karhunen-Lo\`{e}ve expansion
of a known transform of the random variable coefficient is known, we introduce
a compact representation of the random coefficient in the form of a Voronoi
quantizer. The number of Voronoi cells, each of which is represented by a
centroidal variable coefficient, is set to the prescribed number $P$ of
preconditioners. Upon sampling the random variable coefficient, the linear
system assembled with a given realization of the coefficient is solved with the
preconditioner whose centroidal variable coefficient is the closest to the
realization. We consider different ways to define and obtain the centroidal
variable coefficients, and we investigate the properties of the induced
preconditioning strategies in terms of average number of solver iterations for
sequential simulations, and of load balancing for parallel simulations. Another
approach, which is based on deterministic grids on the system of stochastic
coordinates of the truncated representation of the random variable coefficient,
is proposed with a stochastic dimension which increases with the number $P$ of
preconditioners. This approach allows to bypass the need for preliminary
computations in order to determine the optimal stochastic dimension of the
truncated approximation of the random variable coefficient for a given number
of preconditioners.
