The particular properties of systems of linear equations arising in the context of the Stochastic Finite Element Method (SFEM) motivate the customization of existing iterative solution algorithms. The implementation described in this paper has aimed at optimizing data management, MAT-VEC operations and preconditioning strategies. It turns out that SFEM-systems can be solved with much less effort than their size suggests. The main idea is based on the fact that the full system matrix consists of few, relatively small submatrices with identical dimensions and sparsity pattern. This makes it very efficient to perform matrix–vector multiplications at the submatrix level and to avoid the assembly of the full coefficient matrix.