We investigate a one dimensional chain of 2N harmonic oscillators in which neighboring sites have their energies redistributed randomly. The sites −N and N are in contact with thermal reservoirs at different temperature τ − and τ + . Kipnis, Marchioro, and Presutti [18] proved that this model satisfies Fourier's law and that in the hydrodynamical scaling limit, when N → ∞, the stationary state has a linear energy density profileθ(u), u ∈ [−1, 1]. We derive the large deviation function S(θ(u)) for the probability of finding, in the stationary state, a profile θ(u) different fromθ(u). The function S(θ) has striking similarities to, but also large differences from, the corresponding one of the symmetric exclusion process. Like the latter it is nonlocal and satisfies a variational equation. Unlike the latter it is not convex and the Gaussian normal fluctuations are enhanced rather than suppressed compared to the local equilibrium state. We also briefly discuss more general model and find the features common in these two and other models whose S(θ) is known.