CARP-CG: a computational study

Forward modelling of the wave equation is a key ingredient in seismic full waveform inversion (FWI). Simulation in the time domain and solution of the wave equation in the frequency domain are two competing approaches to modelling. Frequency domain approaches can further be categorized as using either direct or iterative solvers. For 3D FWI, iterative solvers in the frequency domain are attractive, partly because they require less memory than the other methods. This is due to the fact that there is no need to store the wavefield at each time step, or compute a factorization of the Helmholtz operator that will not be as sparse as the original matrix. One iterative solver that has been applied to the Helmholtz equation is CARP-CG. CARP-CG uses Kaczmarz row projections for each block of a domain decomposition scheme to precondition the Helmholtz system into being symmetric and positive semidefinite. The method of conjugate gradients is then used to solve the preconditioned system. We present a comparison of the performance of CARP-CG implemented in several languages (MATLAB, C, FORTRAN, julia) and in two different hardware environments: the LIMA HPC cluster hosted at UBC, and the Checkers cluster which is part of the Westgrid consortium. Parallelization of the algorithm via domain decomposition implemented with MPI (distributed memory) and OMP (shared memory) is also examined.