\label{discussion}
In this section, we compare our implementations using both exact subsolves and the Numpy CG solver with respect to time complexity. We also conducted experiments to examine the performance of these solvers under various subdomain configurations, as well as different characteristics of preconditioners and Helmholtz matrices. We will discuss the experiments that met our initial expectations and those that did not align with them. Additionally, we will provide insights into potential future research directions. 
\subsection{Performance of PyOpenCL implementations in terms of time complexity}
As outlined in previous sections, our goal was to accelerate computational processes by harnessing the capabilities of GPUs. Our experiments demonstrated that the PyOpenCL implementation of the Conjugate Gradient method used as an inner solver for both single and multiple right-hand sides (RHSs), significantly outperformed traditional approaches. These results confirm that the computational power of GPUs was effectively utilized, resulting in substantial performance improvements, especially in bigger subdomain setups.

In those experiments, we have set different numbers for CG iterations considering the domain size. However, it is obvious that the number of iterations is not always precise. Setting a higher number of iterations for an algorithm might increase the likelihood of convergence, but this approach can lead to resource wastage if the number is unnecessarily high. Conversely, setting too low a number may prevent the algorithm from converging at all, as it might not perform enough computations to reach an accurate enough solution. This balance requires careful tuning of the iteration count to optimize both computational efficiency and the reliability of results. One solution to that would be to have a stopping criteria for the algorithm. Each iteration of the CG algorithm typically involves a check to determine if the solution has converged sufficiently close to the true solution. This check usually requires calculating the norm of the residual (the difference between the left and right sides of the equation being solved), which is computationally intensive. Therefore, stopping criteria should be analyzed in detail to avoid performance issues.

\subsection{Unexpected behavior of the solvers in higher wavenumbers}
Algorithms were evaluated across various wavenumbers to assess their performance in different media. Surprisingly, higher wavenumbers resulted in fewer iterations and shorter execution times. This phenomenon can be explained as follows: Physically, higher wavenumbers correspond to shorter wavelengths, which may interact with the medium in ways that simplify the solution's structure, thereby reducing the complexity of the calculations required for convergence. Additionally, the matrix representation of the problem at higher wavenumbers might become better conditioned, facilitating faster convergence by reducing the number of iterations needed.

\subsection{Future work}
There are multiple directions in which further improvement could be done for the implemented solver. Currently, we are using even subdomains where the subdomain size is the same. However, an investigation can be done on how to adjust the solvers to work with uneven subdomains by adding a restriction
and an interpolation operator, which translates between different sizes of subdomains into the setup used here for accelerated performance with multiple right-hand-sides solved in one go. For example, areas with higher complexity or finer details might require smaller subdomains to achieve more accurate results, whereas less complex areas can be represented with larger subdomains.

There can be enhancements to the code, such as using multiple GPUs per computation node. Each subdomain being computed independently and in parallel in different GPUs might result in better performance of the code. 

Also, the implementation of the GPGPU CG solver can be enhanced by introducing stopping criteria based on the value of $\delta$ in the Algorithm \ref{alg:cg}. Within CG iterations, the value of $\delta$ (although non-monotonically) decreases. The value can be monitored for each right-hand side vector until all have reduced beneath a given threshold that guarantees the best performance for the overall (outer) iterative process.

The FGMRES solver can use different types of preconditioners at different iterations, greatly enhancing the adaptability and effectiveness of linear system solutions. In fact, it allows one to choose between different preconditioners suitably configured to the characteristics of different subdomains in the solution of a computational problem. In turn, using AI techniques, it is possible to conduct a dynamic choice or even design a preconditioner, which is grounded on real-time analysis of the solver's performance and on the changing conditions of the domain. AI will help automatically tune the strategy of preconditioning so that it optimizes the rates of convergence and computational efficiency. By selecting the most effective preconditioner for each iteration, AI can help minimize the number of iterations needed to reach convergence, thereby speeding up the solution process.
