We develop and analyze an efficient numerical method for computing the response of the solution of an elliptic problem with randomly perturbed coefficients. We use a variational analysis based on the adjoint operator to deal with the perturbations in data. To deal with perturbations in the diffusion coefficient, we construct a piecewise constant approximation to the random perturbation then use domain decomposition to decompose the problem into sub-problems on which the diffusion coefficient is constant. To compute local solutions of the sub-problems, we use the infinite series for the inverse of a perturbation of an invertible matrix to devise a fast way to compute the effects of variation in the parameter. Finally, we derive a posteriori error estimates that take into account all the sources of error and derive a new adaptive algorithm that provides a quantitative way to distribute computational resources between all of the sources.