\label{methods}
This section will provide detailed insights into the implementation of the solvers. It will also introduce the experimental data, discuss the generation of Helmholtz matrices, and outline the specifications of the machines used for conducting these experiments.

\begin{figure}[h!]
    \centering
    \includegraphics[width=1.8\linewidth]{images/Figure_3d_sol.png}
    \caption{3D Visualization of the Helmholtz Equation Solution}
    \label{fig:3dsol}
\end{figure}

\subsection{Accelerating Mathematical Operations in the \cg{} method}
Our primary aim is to develop a CG subsolver in Python using the PyOpenCL library, with the intention of surpassing other solvers implemented using NumPy and Scipy libraries. As previously mentioned, this approach enhances code readability, simplifying maintenance and potential future enhancements. In our application, CG serves as an inner solver to address subdomains generated from the \raas{} method, with FGMRES being our outer Krylov subspace solver. The core of the \cg{} method involves several fundamental operations such as vector additions, scalar-vector multiplications, dot products, and matrix-vector multiplications. When these operations are performed on large data sets or high-dimensional spaces, they can become computationally intensive. Therefore, our implementation consists of five different kernels tailored to make the most of the GPU's computing power for various mathematical operations. As the Helmholtz matrices can be either real or complex, we have 2 different sets of kernels to support operations involving both real and complex values. These kernels enhance performance by enabling faster convergence and facilitating the handling of larger problems that are typical in scientific computing. It is also worth mentioning that the kernels have support for multiple RHSs where the actual performance improvement can be observed. Below, we list these kernels and explain where they are used in the \cg{} method:
\begin{itemize}
    \item \textbf{Axpy} - This kernel is used to calculate $x + \alpha d$ and $r - \alpha q$ in the algorithm \ref{alg:cg}.
    \item \textbf{Aypx} - This kernel is used to calculate $r + \beta d$ in the algorithm \ref{alg:cg}.
    \item \textbf{Spmv} - This kernel is used to compute matrix-vector multiplication. In the CG algorithm \ref{alg:cg}, this kernel calculates $Ax$ and $Ad$.
    \item \textbf{Sub} - This kernel calculates vector subtraction. In the algorithm \ref{alg:cg}, $b-Ax$ uses sub kernel where we already know the result of $Ax$.
    \item \textbf{VDot} - This kernel calculates the vector-vector multiplication. It is used to calculate $delta$ and $dq$ dot operations in the algorithm \ref{alg:cg}.
\end{itemize}
The mentioned kernels are written in a language based on the C programming language. However, OpenCL C, the language used for writing kernels, has some differences and restrictions compared to standard C. These kernels can be utilized within different tools to leverage the computation power of GPU.
Moreover, you can also see the kernel code of the Axpy operation for real numbers in the listing \ref{axpykernel}.

\begin{lstlisting}[language=C, caption=Kernel code of Axpy operation for real numbers, label=axpykernel]
__kernel void axpy(__global float *x,
                   __global float *y,
                   __constant float *a,
                   const int aSign,
                   const int size) {

    const int __idx = get_global_id(0);
    if (__idx < size) {
        if (aSign)
            for (int r = 0; r < N_RHS; r++)
                y[__idx + r * size] += a[r] * x[__idx + r * size];
        else
            for (int r = 0; r < N_RHS; r++)
                y[__idx + r * size] -= a[r] * x[__idx + r * size];
    }
}
\end{lstlisting}
\subsection{Host-device communication}
Host-device communication involves transferring data between the host (typically the CPU) and the device (e.g., GPU) before and after kernel execution. PyOpenCL provides APIs for managing data transfer between the host and the device, as well as for queuing kernel execution commands. APIs provided by PyOpenCL allow us to create a context with a specified GPU. The chosen GPU is responsible for the execution of the kernel. I have separated the initialization of the context and loading of kernels from the actual code in order to avoid the initialization time of context and loading of kernels for each iteration of the FGMRES method. 

\subsection{GPGPU PyOpenCL solver for \bcg{} method}
Our implementation presents 2 solutions for the \bcg{} algorithm with support for complex-valued inputs and multiple right-hand-side. One of them solves the linear equation of $Ax = b$ with a single right-hand side. Here, the $x$ is a single vector that satisfies the equation. When multiple right-hand-sides (RHS) are involved, the problem extends to solving \( AX = B \), where \( B \) is now a matrix of RHS vectors, and \( X \) is a matrix of solution vectors corresponding to each RHS vector. While solving the linear equation with multiple RHSs, computations are done by sending the whole RHS vector of the original equation with the submatrix, and the results are mapped into the whole solution vector. Our setup allows us to use multiple RHS solutions having the exact same subdomains. This allows us to get rid of the iterations through the subdomains.

\subsection{Helmholtz matrices}\label{helmholtz}
Helmholtz matrices arise when discretizing the Helmholtz equation using numerical methods such as finite difference methods, finite element methods (FEM), or boundary element methods (BEM) \cite{bootland2021comparison}. The matrix represents the linear system that approximates the Helmholtz equation. Helmholtz matrices are an essential aspect of numerical methods related to the Helmholtz equation, a fundamental partial differential equation (PDE). In mathematical physics, these matrices represent wave propagation phenomena, such as sound and electromagnetic waves. 

Helmholtz matrices are typically large and sparse because discretization of the domain involves many grid points or elements but has interaction with only those of its local neighbors. These matrices can be indefinite or even nearly singular, especially at higher frequencies, which makes them challenging for standard linear solvers \cite{ernst2011difficult}.

We have a function implemented in Python using NumPy that discretizes Helmholtz Equation in a square domain $\Omega = (0,1)^2$ using the Finite Element Method (FEM) with Robin boundary conditions on all sides before the start of the solution process. The Helmholtz equation is defined as the following \cite{graham2017recent}:
\begin{equation}  
-\Delta u - (k^2 + i\epsilon) u = f \quad \text{on} \, \Omega,
\end{equation}
where $\Delta$ is the Laplacian, $i$ is the imaginary unit, $\epsilon$ is a damping parameter, and $k$ wave number. 
On the boundary $\partial \Omega$, the condition is given by:
\begin{equation}  
\partial_n u - ik u = 0 \quad \text{on}  \, \partial \Omega,
\end{equation}  
where $\partial_n u$ denotes the normal derivative of $u$.

The implemented function has the following inputs that can be modified to explore more about the Helmholtz equations and their convergence behavior in different setups:
\begin{itemize}
  \item $N$: Number of grid points in each direction on the square domain $\Omega$.
  \item $k$: Wave number, $k = \frac{\omega}{c}$ with $\omega$ representing the frequency of a wave and $c$ the variable wave speed.
  \item $\epsilon$: The damping parameter quantifies a wave's energy loss due to absorption, scattering, or other effects in the medium.
\end{itemize}
For simplicity, we are using uniform wave speed in our experiments. However, having a uniform wave speed is also quite common in both theoretical and applied physics. Some examples of Helmholtz problems where the wave speed is almost uniform are Underground Seismic Waves, Medical Ultrasound, and Astronomical Observations. How these changes in these parameters affect Helmholtz matrices' generation will be described later.
The output of the function is a list of exact same-size sparse matrices. The reason for having exact same-size matrices is to ease the work of the domain decomposition method and to be able to solve the problem with multiple right-hand-sides.

Figure \ref{fig:3dsol} illustrates the computed real part of the solution to the Helmholtz equation. The visualization highlights the dynamic behavior of wave propagation within the medium, showing variations in amplitude and phase that are critical to understanding the physical phenomena being modeled. Each peak and trough is represented with a color gradient, which corresponds to the real value intensity, ranging from -4000 to 4000, as detailed in the accompanying color bar.

\subsection{Experimental Environments}
Experiments and benchmarks are conducted at the High-Performance Computing (HPC) Center of Tartu University\cite{https://doi.org/10.23673/ph6n-0144}. For these tests, I utilized the Rocket cluster at the HPC. I specifically requested the Falcon GPU nodes for computational tasks, excluding the Falcon 2 node due to the errors while running OpenCL code. Also, each of these nodes is equipped with 24x NVIDIA Tesla V100 GPUs. The Tesla V100 GPUs come with 5120 CUDA cores, which are fundamental to its ability to handle parallel processing tasks efficiently. Below, you will find the SBATCH script used to submit the job:
\begin{verbatim}
#!/bin/bash
#SBATCH -Jp_helm
#SBATCH --partition=gpu
#SBATCH --exclude=falcon2
#SBATCH --gres=gpu:tesla:4
##SBATCH -N 2
##SBATCH --ntasks-per-node=1
#SBATCH -t 10:00:00
##SBATCH --exclusive
##SBATCH --mem=47000
####SBATCH --nodelist=falcon

S=4
N=128
NIT=512

module load py-numpy/1.22.0
module load py-mpi4py/3.1.4
module load cuda/11.1.0
module load py-scipy/1.8.1
module load py-pyopencl/2020.2.2
module load cmake/3.23.1

cd /gpfs/space/home/ziya/distributed-systems
srun hostname
./p_h-PY_C-CL.py ${S} ${N} 2 ${NIT}
\end{verbatim}
Additionally, you can see the modules that have been loaded to establish an environment suitable for executing the code. The py-pyopencl library is an essential module for executing OpenCL (Open Computing Language) programs in Python. 

\subsection{Evaluation} \label{evaluation}
Four different subsolvers, namely Exact Subsolver (SciPy), PyOpenCL-CG for a single RHS, PyOpenCL-CG for multiple RHSs, and NumPy-CG, are compared in different terms. The time complexity of different implementations was the first metric, varying with the number of subdomains and their sizes. 
The influence of extended overlaps obtained from applying the \raas{} preconditioner was the second metric to observe the convergence of the solution. Additionally, we performed tests with different wavenumbers to see how well the Krylov subspace methods, with the inner iterative solver CG and the outer iterative solver FGMRES, perform under varying medium properties. 