Upcoming and existing large-scale surveys of galaxies require accurate theoretical predictions of the dark matter clustering statistics for thousands of mock galaxy catalogs . We demonstrate that this goal can be achieve with our new Parallel Particle-Mesh ( PM ) N -body code ( PPM-GLAM ) at a very low computational cost . We run about 15,000 simulations with \sim 2 billion particles that provide \sim 1 % accuracy of the dark matter power spectra P ( k ) for wave-numbers up to k \sim 1 h { Mpc } ^ { -1 } . Using this large data-set we study the power spectrum covariance matrix , the stepping stone for producing mock catalogs . In contrast to many previous analytical and numerical results , we find that the covariance matrix normalised to the power spectrum C ( k,k ^ { \prime } ) / P ( k ) P ( k ^ { \prime } ) has a complex structure of non-diagonal components . It has an upturn at small k , followed by a minimum at k \approx 0.1 - 0.2 h { Mpc } ^ { -1 } . It also has a maximum at k \approx 0.5 - 0.6 h { Mpc } ^ { -1 } . The normalised covariance matrix strongly evolves with redshift : C ( k,k ^ { \prime } ) \propto \delta ^ { \alpha } ( t ) P ( k ) P ( k ^ { \prime } ) , where \delta is the linear growth factor and \alpha \approx 1 - 1.25 , which indicates that the covariance matrix depends on cosmological parameters . We also show that waves longer than 1 \mbox { $h ^ { -1 } $Gpc } have very little impact on the power spectrum and covariance matrix . This significantly reduces the computational costs and complexity of theoretical predictions : relatively small volume \sim ( 1 \mbox { $h ^ { -1 } $Gpc } ) ^ { 3 } simulations capture the necessary properties of dark matter clustering statistics . All the power spectra obtained from many thousands of our simulations are publicly available .