Frontier scientific and AI workloads now reach 1019−1025 fused multiply–add (FMA) operations per run (on the order of 2×1019−2×1025 FLOPs). At today's ∼10 pJ per FMA, this corresponds to approximately 108−1014 joules of arithmetic energy. At this scale, energy becomes the limiting resource for continued growth in computational workloads, motivating a re-evaluation of long-standing algorithmic assumptions. It is often assumed that reversible computing only matters near the Landauer limit. Building on prior physical arguments that full energy recovery is only possible when computation preserves information, we demonstrate that this same requirement governs floating-point numerical kernels: overwriting state enforces a non-zero energy floor, even under ideal recovery. Thus, eliminating this wall in practice requires that the numerical algorithm itself be injective. We therefore present the first reversible floating-point realizations of core dense numerical kernels—matrix multiplication, LU factorization, and conjugate-gradient iteration—that retain rounding information rather than discarding it. Implemented directly in IEEE arithmetic, they achieve machine-precision forward–reverse agreement on well- and ill-conditioned problems with minimal auxiliary state. A toggle-based model with measured switching costs and realistic recovery factors predicts 103−104× reductions in arithmetic energy. These results establish injective floating-point kernels as a foundation for energy-recovering numerical computation, and indicate that realizing this potential will require sustained co-design across applied mathematics, computer science, and hardware engineering.

We present a matrix-free parallel scalable multilevel deflation preconditioned method for heterogeneous time-harmonic wave problems. Building on the higher-order deflation preconditioning proposed by Dwarka and Vuik (SIAM J. Sci. Comput. 42(2):A901-A928, 2020; J. Comput. Phys. 469:111327, 2022) for highly indefinite time-harmonic waves, we adapt these techniques for parallel implementation in the context of solving large-scale heterogeneous problems with minimal pollution error. Our proposed method integrates the Complex Shifted Laplacian preconditioner with deflation approaches. We employ higher-order deflation vectors and re-discretization schemes derived from the Galerkin coarsening approach for a matrix-free parallel implementation. We suggest a robust and efficient configuration of the matrix-free multilevel deflation method, which yields a close to wavenumber-independent convergence and good time efficiency. Numerical experiments demonstrate the effectiveness of our approach for increasingly complex model problems. The matrix-free implementation of the preconditioned Krylov subspace methods reduces memory consumption, and the parallel framework exhibits satisfactory parallel performance and weak parallel scalability. This work represents a significant step towards developing efficient, scalable, and parallel multilevel deflation preconditioning methods for large-scale real-world applications in wave propagation.

The Helmholtz equation is related to seismic exploration, sonar, antennas, and medical imaging applications. It is one of the most challenging problems to solve in terms of accuracy and convergence due to the scalability issues of the numerical solvers. For 3D large-scale applications, high-performance parallel solvers are also needed. In this paper, a matrix-free parallel iterative solver is presented for the three-dimensional (3D) heterogeneous Helmholtz equation. We consider the preconditioned Krylov subspace methods for solving the linear system obtained from finite-difference discretization. The Complex Shifted Laplace Preconditioner (CSLP) is employed since it results in a linear increase in the number of iterations as a function of the wavenumber. The preconditioner is approximately inverted using one parallel 3D multigrid cycle. For parallel computing, the global domain is partitioned blockwise. The matrix-vector multiplication and preconditioning operator are implemented in a matrix-free way instead of constructing large, memory-consuming coefficient matrices. Numerical experiments of 3D model problems demonstrate the robustness and outstanding strong scaling of our matrix-free parallel solution method. Moreover, the weak parallel scalability indicates our approach is suitable for realistic 3D heterogeneous Helmholtz problems with minimized pollution error.

We propose a matrix-free parallel two-level deflation method combined with the Complex Shifted Laplacian Preconditioner (CSLP) for two-dimensional heterogeneous Helmholtz problems encountered in seismic exploration, antennas, and medical imaging. These problems pose challenges in terms of accuracy and convergence due to scalability issues with numerical solvers. Motivated by the limitations imposed by excessive computational time and memory constraints when employing a sequential solver with constructed matrices, we parallelize the two-level deflation method without constructing any matrices. Our approach utilizes preconditioned Krylov subspace methods and approximates the CSLP preconditioner with a parallel geometric multigrid V-cycle. For the two-level deflation, standard inter-grid deflation vectors and further high-order deflation vectors are considered. As another main contribution, the matrix-free Galerkin coarsening approach and a novel re-discretization scheme as well as high-order finite-difference schemes on the coarse grid are studied to obtain wavenumber-independent convergence. The optimal settings for an efficient coarse-grid problem solver are investigated. Numerical experiments of model problems show that the wavenumber independence has been obtained for medium wavenumbers. The matrix-free parallel framework shows satisfactory weak and strong parallel scalability.

We present a matrix-free parallel iterative solver for the Helmholtz equation related to applications in seismic problems and study its parallel performance. We apply Krylov subspace methods, GMRES, Bi-CGSTAB and IDR(s), to solve the linear system obtained from a second-order finite difference discretization. The Complex Shifted Laplace Preconditioner (CSLP) is employed to improve the convergence of Krylov solvers. The preconditioner is approximately inverted by multigrid iterations. For parallel computing, the global domain is partitioned blockwise. The standard MPI library is employed for data communication. The matrix-vector multiplication and preconditioning operator are implemented in a matrix-free way instead of constructing large, memory-consuming coefficient matrices. These adjustments lead to direct improvements in terms of memory consumption. Numerical experiments of model problems show that the matrix-free parallel solution method has satisfactory parallel performance and weak scalability. It allows us to solve larger problems in parallel to obtain more accurate numerical solutions.

The bottleneck in designing iterative solvers for the Helmholtz equation lies in balancing the trade-off between accuracy and scalability. Both the accuracy of the numerical solution and the number of iterations to reach convergence deteriorate in higher dimensions and increase with the wavenumber. To address these issues in this dissertation, we formulated three research pillars: accuracy, wavenumber independent convergence and linear complexity. Below, we summarize the core findings of this dissertation: WAVENUMBER INDEPENDENT CONVERGENCE We develop the first preconditioning technique which leads to close to wavenumber independent convergence for very large wavenumbers in 1D, 2D and 3D. Building on a two-level deflation projection method, we incorporated Quadratic Rational Bezier curves to construct the deflation space and vectors (Chapter 7). As a result, the near-zero eigenvalues of the coarse grid operator remain aligned with the fine-grid operator, keeping the spectrum of the preconditioned system clustered, leading to superior convergence properties compared to previous methods. LINEAR COMPLEXITY For over 30 years, applied mathematicians have tried to make convergent (standard) multigrid solvers for the Helmholtz equation. Multigrid solvers use sequences of smaller problem sizes and are computationally cheap and easy to implement. Unfortunately, multigrid methods diverge for Helmholtz and solving this issue remained an open problem. Using standard smoothing techniques, combined with similar higher-order coarse spaces, we constructed a fully convergent V- and W-cycle algorithm (Chapter 9). The key features of the algorithm are the use of higher-order transfer operators (instead of deflation vectors in the previous application) and a complex shift in the smoothing operator. While the method converges and the preliminary results have been proven, much research can still be conducted in this area, as this could support a paradigm shift in solving the complexity issue for very large wavenumbers in 2D and 3D. In light of this, we extended the two-level deflation solver to a multi-level deflation solver to address both the issue of wavenumber and problem size dependence (Chapter 8). In this part, we show better convergence properties and provide numerical experiments on challenging 2D and 3D test problems to corroborate the theoretical results. ACCURACY Finally, we developed an unprecedented way to study the accuracy of the numerical solutions by studying the eigenvalues of systems where the analytical solution is known (Chapter 5). Expressing the pollution error in terms of these eigenmodes, enabled theoretical accuracy studies and dispersion corrections in higher dimensions, irrespective of the wave propagation angles. Something which was previously impossible. We also studied the application of Isogeometric Analysis (IgA) to improve the accuracy and reduce the pollution error (Chapter 6). Our results showed that the use of IgA was able to significantly suppress the pollution error compared to Finite Elements Discretizations of the same order.

In recent years, domain decomposition based preconditioners have become popular tools to solve the Helmholtz equation. Notorious for causing a variety of convergence issues, the Helmholtz equation remains a challenging PDE to solve numerically. Even for simple model problems, the resulting linear system after discretisation becomes indefinite and tailored iterative solvers are required to obtain the numerical solution efficiently. At the same time, the mesh must be kept fine enough in order to prevent numerical dispersion ‘polluting’ the solution [4]. This leads to very large linear systems, further amplifying the need to develop economical solver methodologies.

Recent research efforts aimed at iteratively solving time-harmonic waves have focused on a broad range of techniques to accelerate convergence. In particular, for the famous Helmholtz equation, deflation techniques have been studied to accelerate the convergence of Krylov subspace methods. In this work, we extend the two-level deflation method to a multilevel deflation method for (heterogeneous) Helmholtz and elastic wave problems. By using higher-order deflation vectors, we show that up to the level where the coarse-grid linear systems remain indefinite, the near-zero eigenvalues of the these coarse-grid operators remain aligned with the near-zero eigenvalues of the fine-grid operator, keeping the spectrum of the preconditioned system away from the origin. Combining this with the well-known CSLP-preconditioner, we obtain a scalable solver for the highly indefinite linear systems. This can be attributed to a close to wave number independent convergence and an optimal use of the CSLP-preconditioner on the indefinite levels. There, we approximate the CSLP-preconditioner, while allowing the complex shift to be small, by using inner Bi-CGSTAB iterations instead of a multigrid F-cycle. The proposed method shows very promising results for the more challenging two- and three-dimensional heterogeneous time-harmonic wave problems.

Finding fast yet accurate numerical solutions to the Helmholtz equation remains a challenging task. The pollution error (i.e. the discrepancy between the numerical and analytical wave number k) requires the mesh resolution to be kept fine enough to obtain accurate solutions. A recent study showed that the use of Isogeometric Analysis (IgA) for the spatial discretization significantly reduces the pollution error. However, solving the resulting linear systems by means of a direct solver remains computationally expensive when large wave numbers or multiple dimensions are considered. An alternative lies in the use of (preconditioned) Krylov subspace methods. Recently, the use of the exact Complex Shifted Laplacian Preconditioner (CSLP) with a small complex shift has shown to lead to wave number independent convergence while obtaining more accurate numerical solutions using IgA. In this paper, we propose the use of deflation techniques combined with an approximated inverse of the CSLP using a geometric multigrid method. Numerical results obtained for one- and two-dimensional model problems, including constant and non-constant wave numbers, show scalable convergence with respect to the wave number and approximation order p of the spatial discretization. Furthermore, when kh is kept constant, the proposed approach leads to a significant reduction of the computational time compared to the use of the multigrid-approximated or exact inverse of the CSLP with a small shift, in particular for three-dimensional model problems.

In researching the Helmholtz equation, the focus has either been on the accuracy of the numerical solution (pollution) or the acceleration of the convergence of a preconditioned Krylov-based solver (scalability). While it is widely recognized that the convergence properties can be investigated by studying the eigenvalues, information from the eigenvalues is not used in studying the numerical dispersion which drives the pollution error. Our aim is to bring the topics of accuracy and scalability together for the first time; instead of approaching the pollution error in the conventional sense of being the result of a discrepancy between the exact and numerical wavenumber, we show that the dispersion which drives the pollution error can also be decomposed in terms of the eigenvectors and eigenvalues. Using these novel insights, we construct sharper upper bounds for the total error independent of the grid resolution. While the pollution error can be minimized in one-dimension by introducing a dispersion correction, the latter is not possible in higher dimensions, even for very simple model problems. For our model problem, a correction on the eigenvalues enables us to remove the pollution error and study it in full detail, both in one- and two-dimensions.

Recent research efforts aimed at iteratively solving the Helmholtz equation have focused on incorporating deation techniques for accelerating the convergence of Krylov subpsace methods. The requisite for these efforts lies in the fact that the widely used and well-acknowledged complex shifted Laplacian preconditioner (CSLP) shifts the eigenvalues of the preconditioned system towards the origin as the wave number increases. The two-level-deation preconditioner combined with CSLP showed encouraging results in moderating the rate at which the eigenvalues approach the origin. However, for large wave numbers the initial problem resurfaces and the near-zero eigenvalues reappear. Our findings reveal that the reappearance of these near-zero eigenvalues occurs if the near-singular eigenmodes of the fine-grid operator and the coarse-grid operator are not properly aligned. This misalignment is caused by accumulating approximation errors during the inter-grid transfer operations. We propose the use of higher-order approximation schemes to construct the deation vectors. The results from rigorous Fourier analysis and numerical experiments confirm that our newly proposed scheme outperforms any other deation-based preconditioner for the Helmholtz problem. In particular, the spectrum of the adjusted preconditioned operator stays fixed near one. These results can be generalized to general shifted indefinite systems with random right-hand sides. For the first time, the convergence properties for very large wave numbers (k = 106 in one dimension and k = 103 in two dimensions) have been studied, and the convergence is close to wave number independence. Wave number independence for three dimensions has been obtained for wave numbers up to k = 75. The new scheme additionally shows very promising results for the more challenging Marmousi problem. Irrespective of the strongly varying wave number, we obtain a constant number of iterations and a reduction in computational time as the results remain robust without the use of the CSLP preconditioner.

Recent research efforts aimed at iteratively solving the Helmholtz equation have focused on incorporating deflation techniques for accelerating the convergence of Krylov subspace methods. In this work, we extend the two-level deflation method in [6] to a multilevel deflation method. By using higher-order deflation vectors, we show that up to the level where the coarse-grid linear systems remain indefinite, the near-zero eigenvalues of the these coarse-grid operators remain aligned with the fine-grid operator keeping the spectrum of the preconditioned system fixed away from the origin. Combining this with the well-known CSLP-preconditioner, we obtain a scalable solver with theoretical linear complexity for the highly indefinite Helmholtz equation. This can be attributed to a fixed number of iterations independent of the wave number and an optimal use of the CSLP-preconditioner. We approximate the CSLP-preconditioner, while allowing the complex shift to be small. The proposed configuration additionally shows very promising results for the more challenging Marmousi problem.