Many physical problems can be modeled by partial differential equations with stochastic coefficients. Different coefficients in the same problem may be perturbed by different random functions. In this talk, we discuss  multilevel Schwarz type domain decomposition preconditioned Krylov subspace methods for the numerical solution of second order elliptic equations with multiple stochastic coefficients. We show that the Karhunen-Loeve expansion and double orthogonal polynomials can be used to reformulate the stochastic problem into a large number of uncoupled deterministic equations. We also discuss a PETSc based parallel implementation of a recycling Krylov subspace method and show numerical results obtained on machines with a large number of processors.