Scaling through abstractions – high-performance vectorial wave simulations for seismic inversion with Devito
TABLE I: Per-grid-point FLOPs of the finite-difference TTI wave-equation stencil with different spatial discretization orders.
['spatial order', 'w/o optimizations', 'w/ optimizations']
Owing to the very high number of floating-point operations (FLOPs) needed per grid point for the weighted rotated Laplacian, this anisotropic wave-equation is extremely challenging to implement. The version without FLOP-reducing optimizations is a direct translation of the discretized operators into stencil expressions (see tti-so8-unoptimized.c). The version with optimizations employs transformations such as common sub-expressions elimination, factorization, and cross-iteration redundancy elimination – the latter being key in removing redundancies introduced by mixed derivatives. Implementing all of these techniques manually is inherently difficult and laborious. Further, to obtain the desired performance improvements it is necessary to orchestrate them with aggressive loop fusion (for data locality), tiling (for data locality and tensor temporaries), and potentially ad-hoc vectorization strategies (if rotating registers are used). what is emphasized here is that users can easily take full advantage of these optimizations without needed to concern themselves with the details.