The radiative transfer equations have three spatial dimensions, two radiational dimensions, one frequency variable, and possibly time. In this work, without loss of generality, they are reduced to a nonlinear integro-differential system and solved by an iterative scheme which is shown to be monotone and convergent. At each step, two large integrals of convolution type need to be computed. The kernel matrices are assembled in compressed form using Hierarchical matrices, with the help of the H-Matrix high performance parallel library Htool. The whole method is implemented in FreeFEM [3], a free open-source software designed for multi-physics simulations. The resulting method is of complexity n log(n) where n is the number of vertices of the physical domain. As an example, the method is applied to the temperature and radiation in the valley of Chamonix (France), the challenge being the large variation of density, altitude and light reflection due to snow and clouds. Our method is able to handle 240K physical points, all directions of radiation and 683 frequencies in less than 35 minutes on an Apple M1 laptop [1]. The method has also been extended to handle reflective boundary conditions [2].