In micromagnetic simulations, the demagnetization field is by far the computationally most expensive field component and often a limiting factor in large multilayer systems. We present an exact method to calculate the demagnetization field of magnetic layers with arbitrary thicknesses. In this approach we combine the widely used fast-Fourier-transform based circular convolution method with an explicit convolution using a generalized form of the Newell formulas. We implement the method both for central processors and graphics processors and find that significant speedups for irregular multilayer geometries can be achieved. Using this method we optimize the geometry of a magnetic random-access memory cell by varying a single specific layer thickness and simulate a hysteresis curve to determine the resulting switching field.