The elastic contact of rough surfaces and the subsurface stresses caused by the contact have been analyzed by means of a numerical model based on fast Fourier transforms (FFT) and minimization of complementary energy. The elastic contact has been modeled mathematically as a linear complementarity problem and solved by a robust algorithm, Conjugate Gradient Method, while the force-displacement relation is determined through a FFT approach. After solving for the pressure distribution, the subsurface stress field is obtained by calculating the stresses due to the application of a point force, and then integrating over the contact region. In comparison with the matrix based method published in recent years, the numerical approach presented in this study is more efficient, more stable and requires less memory. It has great potential in application to problems with general contact geometry and three-dimensional surface roughness. The results show that high frequency roughness could lead to very sharp impulses and significant oscillations in pressure distribution. The amplitude and location of the maximum shear stress in the subsurface region are constantly changing when the contacting rough surfaces are in relative motion.Presented at the 53rd Annual Meeting in Detroit, Michigan May 17–21, 1998