Fast finite-difference convolution for 3D problems in layered media

Abstract: We developed fast direct solver for 3D Helmholtz and Maxwell equations in layered medium. The algorithm is based on the ideas of cyclic reduction for separable matrices. For the grids with major uniform part (within the survey domain in the problems of geophysical prospecting, for example) and small non-uniform part (PML and coarsening to approximate problems in infinite domain) the computational cost of our approach is $O(N_xN_ylog(N_xN_y)N_z)$. For general non-uniform grids the cost is $O(N^{3/2}_xN^{3/2}_yN_z)$. The first asymptotics coincide with the cost of FFT-based methods, which can be applied for uniform gridding (in x and y) only. Our approach is significantly more efficient compared to the algorithms based on discrete Fourier transform which cost is $O(N^2_xN^2_yN_z)$. The algorithm can be easily extended for solving the elasticity problems as well.