This talk presents (parallel) preconditioners and multigrid solvers for solving linear systems of equations arising from stochastic polynomial chaos (PC) formulations of the diffusion equation with random coefficients. These preconditioners/solvers are an extension of the preconditioner developed in [2] for strongly coupled systems of elliptic partial equations (PDEs) that are norm equivalent to systems that can be factored into an algebraic coupling component and a diagonal differential component, as is the case for the systems produced by PC formulations of the diffusion equation. The first preconditioner, which is applied to the norm equivalent system, is obtained by sparsifying the inverse of the algebraic coupling component. This sparsification leads to an efficient method for solving these PC systems on the large scale, even for problems with large statistical variations in the random coefficients. An extension of this preconditioner leads to stand-alone multigrid methods that are applied directly to the actual system of PDEs rather than to the norm equivalent system. These multigrid methods exploit the algebraic/differential factorization of the norm equivalent systems to produce variable decoupled systems on the coarse levels, and lead to simplified software implementation through re-use of existing high-performance software such as the Hypre library package ([1]). Two-grid matrix bounds are established and numerical results are given.