We consider a chemical reaction network governed by mass action kinetics and composed of N different species which can reversibly form heterodimers. A fast iterative algorithm is introduced to compute the equilibrium concentrations of such networks. We show that the convergence is guaranteed by the Banach fixed point theorem. As a practical example, of relevance for a quantitative analysis of microarray data, we consider a reaction network formed by N~10^6 mutually hybridizing different mRNA sequences. We show that, despite the large number of species involved, the convergence to equilibrium is very rapid for most species. The origin of slow convergence for some specific subnetworks is discussed. This provides some insights for improving the performance of the algorithm.