We present a boundary integral equation/conjugate gradient formulation to study the propagation of seismic waves through complex geological structures. The method is aimed at extending the range of applications of boundary integral equations or boundary element methods to geological models of relatively large size or complexity. We show that the system of equations that expresses the boundary conditions at the medium interfaces and that is inherent to the biem or bem approach can be drastically reduced in size and that only 10 to 20% of the terms of this system contribute significantly to the solution. The boundary conditions may thus be expressed in the form of a very sparse linear system that can be inverted iteratively by the conjugate gradient method. We use this approach to investigate the effect of a sedimentary layer of varying thickness on local and regional seismic wave propagation. We show that the amplitude of ground motion and the duration of shaking are drastically increased when the sediment/bedrock interface is rough and irregular relatively to the case where it is flat.