The use of linear multivariate mixed models in animal genetics, leads to very large, sparse linear systems of equations. The sparse, symmetric coefficient matrix is too large to be constructed explicitly. We describe a parallel, iterative linear equation solver for large sparse systems, developed by DIAS and UNI-C. The solver takes advantage of the structure of the multivariate mixed model equations, and is based on Gauss-Seidel and second order Jacobi iteration. It is parallelized for distributed memory architectures.