\label{background}
This section will give an outlook on the application of GPGPU programming using OpenCL to accelerate the Krylov subspace methods for solving Helmholtz equations and provide background information about the separate components of the implementation. 

\subsection{General-purpose computing on graphics processing units (GPGPU) computations}
Recently, the Graphical Processing Unit (GPU) has shown great potential for the future of computation. Initially designed for fast image rendering in graphical applications, GPUs are now widely used for non-graphical applications, a field known as GPU computing. The computational capability of GPUs is evolving into powerful parallel computing units. Hence, General-purpose computing on graphics processing units (GPGPU) is a technique that involves utilizing the processing power of graphics processing units (GPUs) for general-purpose computing tasks that were previously performed by the central processing unit (CPU).

The architecture of GPUs is designed to handle a large number of computations simultaneously. It makes GPUs more powerful than CPUs in certain types of computations. The main difference between CPU and GPU architecture is the number of cores. A CPU typically has a few cores optimized for sequential processing, while a GPU has thousands of smaller cores optimized for parallel processing. This makes GPUs more suitable for tasks that require big parallelisms, such as image and video processing, scientific simulations, and machine learning. 

Furthermore, GPUs have much faster memory and memory bandwidth than CPUs. That allows GPUs to move data between the processor and memory quickly. This is critical for performance in data-intensive applications. Additionally, modern GPUs are designed to be highly programmable and can be programmed using a variety of programming languages, including CUDA, OpenCL, and others.

CUDA (Compute Unified Device Architecture) is a parallel computing platform and programming model created by NVIDIA. It allows developers to use GPUs for general-purpose computing tasks in addition to their traditional use in graphics processing. The CUDA is a parallel computing platform developed by NVIDIA and optimized only for NVIDIA GPUs. However, depending on the underlying hardware, it might bring some performance overhead and limitations. On the other hand, OpenCL is an open standard for parallel programming across different hardware architectures, including GPUs, CPUs, and FPGAs. This means that OpenCL can be used with various hardware vendors. It makes OpenCL more adaptable and practical compared to CUDA \cite{du2012cuda,fang2011comprehensive}.\\

\begin{figure}[h!]
    \centering
    \includegraphics[scale=0.5]{images/memory-model-opencl.png}
    \caption{OpenCL memory model. LDS is an abbreviation for Local data storage}
    \label{fig:opencl-memory-model}
\end{figure}
\subsection{Execution of a program in OpenCL} 
In today's era of big data and complex computational tasks, the demand for high-performance computing solutions has never been greater. Traditional approaches to parallel computing often face challenges in achieving portability across different hardware platforms. OpenCL addresses these challenges by providing a standardized framework for parallel programming that enables developers to write efficient, portable code across a wide range of devices. It is used in a wide range of applications and use cases, including scientific simulations, image and video processing, machine learning, and computational finance. Real-world examples demonstrate the performance benefits and scalability achieved by using OpenCL in these applications.


OpenCL follows a host-device model, where the host CPU coordinates the execution of parallel kernels on compute devices such as GPUs and other accelerators. The architecture of OpenCL consists of compute devices, compute units, work-items, work-groups, and different memories. The memory hierarchy includes global, local, and private memory spaces; see Figure \ref{fig:opencl-memory-model}. This architecture enables fine-grained parallelism and efficient data movement between the host and compute devices. The programming model of OpenCL involves writing kernels, which are parallel functions written in the OpenCL C language. Kernels are executed on compute devices in parallel by a large number of work-items, with each work-item representing a single execution thread. Data parallelism is achieved by partitioning the computation across multiple work-items, while task parallelism is achieved by launching multiple kernels concurrently \cite{kaeli2015heterogeneous, stone2010opencl}.

To start writing an OpenCL program, we must ensure that OpenCL SDK is correctly installed for our platform. It includes the necessary libraries, header files, and tools for writing and building OpenCL programs. The term platform in OpenCL is the hardware and software environment in which OpenCL runs. There can be multiple platforms on a single computer. Each platform has its own set of devices. A compute device is the physical hardware on which OpenCL code runs. Examples of devices include CPUs, GPUs, and FPGAs. GPUs consist of a large number of compute units (CUs), and each CU has a set of cores that can execute instructions in parallel. These cores are also called processing elements. In addition to the processing cores, GPU architecture has several memory types, including global memory, shared memory, and registers. Global memory is accessible to all cores and is used to store data that needs to be accessed by all threads in a program. Shared memory is a fast, low-latency memory shared between threads within a thread block. Registers are fast, on-chip memory that stores data and intermediate values during computations. Additionally, each processing element has its own private memory, accessible only by one work item or thread \cite{komatsu2010evaluating}.

In OpenCL programs, the host device communicates with the compute device by sending commands through a command queue. These commands can be for kernel execution or data transfers. Each command is executed in the same order as in the queue and asynchronously, meaning the host can add commands to the queue while the device executes previous commands. The compute device receives the command from the command queue to execute the kernel on a workgroup. A workgroup is a group of work items that execute, communicate, and synchronize using shared memory and barrier synchronization on the same compute unit. A workgroup is defined by its size, specified by the developer when the OpenCL kernel is launched. On the other hand, a wavefront is a group of work items that execute the same instruction in a lock-step fashion, and its size depends on the hardware architecture of the GPU, which is usually 32 or 64. Once all the workgroups have been executed, the compute device sends the results back to the host device \cite{tompson2012introduction}.

There are several ways to develop OpenCL programs. The most renowned is C language-based development. While C provides a robust API for writing efficient parallel code across different computing devices, its complexity can cause some challenges. Hence, alternative options are explored and discussed in the following section.

\subsection{GPU programming with PyOpenCL} 
Nowadays, we are witnessing the different complex algorithms that solve real-life applications. Therefore, an increase in the time complexity of algorithms has led to a need to accelerate them using GPUs.
PyOpenCL is a Python library used to program heterogeneous systems with OpenCL. It provides a way for Python programmers to harness the power of OpenCL from within their Python scripts, allowing them to write high-performance code that runs on various hardware architectures. It also allows the use of multiple GPUs for a single Python program. There are many advantages to using the PyOpenCL Python library over the C programming language. I will list some of them below:
\begin{itemize}
    \item \textbf{Simplicity} - Python is a high-level programming language with a simpler syntax than C, which is also easy to debug. This can make development faster and more efficient.
    \item \textbf{Easy-to-use API} - PyOpenCL provides a simpler and more intuitive API than the low-level C API, making it easier to write and understand code.
    \item \textbf{Easy integration with different libraries} - PyOpenCL integrates well with the Python ecosystem, which includes libraries like NumPy, SciPy, and Matplotlib. This makes it a popular choice for scientific computing applications.
    \item \textbf{Cross-platform support} - PyOpenCL is designed to work on various operating systems and hardware platforms, including Windows, macOS, Linux, and others. 
\end{itemize}

The PyOpenCL platform has different sections for both host code and device code. The host code encompasses typical Python operations like variable definition and regular code execution. In contrast, the device code is responsible for actual implementation. This device code operates within kernels. PyOpenCL views GPUs and accelerators as assemblies of multiple compute devices, each subdivided into units containing numerous processing elements. This organizational structure facilitates the execution of functions.
Development in PyOpencCl starts with defining the context that binds all the objects inside the code. It can be created by providing the list of devices as \mintinline{python}{cl.Context(devices=[devices])}. The context contains all the information about devices, kernels, memory, and variables in the created context.
The command queue defines the order in which the commands are executed. Queue is defined as \mintinline{python}{cl.CommandQueue(ctx)}, and used to enqueue commands for execution on a specific device associated with the context. We can pass data to kernels by buffers. It is defined as \mintinline{python}{cl.Buffer(ctx, mf.READ_WRITE, size=256, hostbuf=a_values)}. In the mentioned command, we define the read-write permission of the buffer and its size. Moreover, we pass the value \mintinline{python}{a_values} that is a numpy array. The most important part of PyOpenCL is defining the kernels. The kernel code is written as a string inside \mintinline{python}{cl.Program(ctx, kernel)} \cite{pyopencl-docs}.\\
\newpage
\begin{lstlisting}[language=Python, caption=PyOpenCL code for matrix multiplication, label=pyopenclmatmul]
program = cl.Program(ctx, """
__kernel void mult_matrices(const unsigned int size, __global float * first, __global float * second, __global float * result) {
    int i = get_global_id(1);
    int j = get_global_id(0);
    result[i + size * j] = 0;
    for (int k = 0; k < size; k++) {
        result[j + size * i] += first[k + size * i] * second[j + size * k];
    }
} """).build()

program.multmatrices(queue, first.shape, None, np.int32(len(first)), first_buf, second_buf, result_buf)

cl.enqueue_copy(queue, result, result_buf)
\end{lstlisting}
\hfill \break
The provided implementation [\ref{pyopenclmatmul}] demonstrates a Matrix Multiplication program utilizing PyOpenCL. It initializes two random matrices using NumPy arrays, passed to the kernel via buffers. The resulting matrix is stored in a destination matrix labeled 'result'. PyOpenCL uses the `enqueue\_copy` function to copy data between different memory objects. Specifically, it copies data from one buffer to another. This function allows you to efficiently transfer data between different memory regions, such as between the host and device memory or between different regions of device memory. In our code, \mintinline{python}{cl.enqueue_copy(queue, result, result_buf)} line copies the data from buffer `result\_buf` to the NumPy array.\\

\begin{lstlisting}[language=Python, caption=Python code for matrix multiplication, label=pythonmatmul]
for i in range(first[0].size):
    for j in range(second[0].size):
        for k in range(second[0].size):
            result[i][j] += first[i][k] * second[k][j]
\end{lstlisting}
\hfill \break
The data presented in Table \ref{table:matrixmulcomparison} shows that the PyOpenCL implementation considerably outperforms the regular Python implementation [\ref{pythonmatmul}]. For instance, while the PyOpenCL implementation executes the multiplication in merely 0.183 seconds for a matrix size of 1024, the regular Python code requires a significantly longer time of 1379 seconds to accomplish the same task. 
Observations indicate that PyOpenCL demonstrates performance that increases exponentially as the size of the matrices increases. This phenomenon arises from the efficient method of passing an entire vector from the matrix as a single work item to the kernel. It indicates that using PyOpenCL for matrix multiplication offers considerable performance improvements over regular Python code.

\begin{table}[htbp]
    \centering
    \caption{Execution Time Comparison: PyOpenCL vs. Python for Matrix Multiplication}
    \begin{tabularx}{\textwidth}{@{}XXXX@{}}
    \toprule
    \textbf{Size of Matrix} & \textbf{PyOpenCL (sec)} & \textbf{Normal Python (sec)} \\ \midrule
    128                      & 0.146                    & 2.613                        \\
    256                      & 0.148                    & 20.627                       \\
    512                      & 0.151                    & 170.378                      \\
    1024                     & 0.183                    & 1378.951                     \\ \bottomrule
    \end{tabularx}
    \label{table:matrixmulcomparison}
\end{table}

\subsection{Krylov subspace methods}
Solving large linear systems is a fundamental task in scientific computing, engineering, and other fields. Direct methods, while accurate, face scalability challenges for large-scale problems due to their computational complexity and memory requirements. Iterative solvers, particularly Krylov subspace methods, offer a promising alternative for efficiently solving large sparse linear systems. Krylov subspace methods are a family of iterative solvers used to solve large, sparse linear systems of equations. These methods are particularly effective when direct methods, such as LU decomposition or Cholesky factorization, become computationally unusable due to the size and sparsity of the coefficient matrix.

The fundamental idea behind Krylov subspace methods is to iteratively construct a sequence of approximations to the solution vector by generating a sequence of subspaces known as Krylov subspaces. These subspaces are formed by repeatedly applying the coefficient matrix of the linear system to an initial vector, see \cite{saad2003iterative}.

Given a matrix $A$ and a starting vector $r_0$, the Krylov subspace of dimension $n$, denoted as $\mathcal{K}_n(A, r_0)$, is the linear span of the vectors formed by repeatedly applying the matrix $A$ to $r_0$ up to $n$ times. Mathematically, it is defined as:
\begin{equation}
 \mathcal{K}_n(A, r_0) = \text{span}\{r_0, Ar_0, A^2r_0, \ldots, A^{n-1}r_0\}. 
\end{equation}
Krylov subspaces have several key properties:

\begin{itemize}
    \item \textbf{Dimension}: The dimension of the Krylov subspace is at most $n$, where $n$ is the number of iterations or the desired dimension of approximation.
    
    \item \textbf{Basis}: The vectors forming the basis of the Krylov subspace are obtained by repeatedly applying the matrix $A$ to the initial vector $r_0$.
    
    \item \textbf{Approximation Space}: Krylov subspaces provide a way to approximate the action of $A$ on any vector. By using a low-dimensional Krylov subspace, we can obtain an approximation to the action of $A$ that can be computationally more efficient than directly computing $A$ times a vector.
    
    \item \textbf{Convergence}: Iterative methods like the Conjugate Gradient method utilize Krylov subspaces to iteratively approximate the solution of linear systems. The convergence properties of these methods are related to the properties of the Krylov subspaces, such as their dimension and the properties of the matrix $A$.
\end{itemize}

Let's delve into a bit more detail on how the Krylov subspace is utilized in the \cg{} method. The vectors generated within the Krylov subspace, $K_n(A, r_0)$, are not only linear combinations of powers of the initial residual $r_0$ but also inherently orthogonal to each other concerning the matrix $A$. This orthogonality property is leveraged in the \cg{} method to ensure that the search directions at each iteration are conjugate to each other. In the context of the \cg{} method, when we say that vectors are conjugate to each other, we mean that they are orthogonal with respect to the matrix $A$ being used in the linear system. Mathematically, two vectors $u$ and $v$ are considered conjugate with respect to a symmetric positive definite matrix $A$ if their inner product (also known as the dot product) satisfies the following condition:
\begin{equation}
     u^T A v = 0. 
\end{equation}

This condition indicates that the vectors $u$ and $v$ are perpendicular to each other in the space defined by the matrix $A$. In other words, they are orthogonal with respect to the inner product induced by $A$.

Therefore, if $p_k$ and $p_{k+1}$ are two consecutive search directions in the CG algorithm, then $p_k^T A p_{k+1} = 0$. Conjugacy ensures that the CG method moves efficiently through the search space, avoiding redundant directions and converging to the solution in fewer iterations.

The vectors generated in the Krylov subspace form a basis for constructing approximate solutions to the linear system. Krylov subspace methods iteratively refine these approximate solutions to converge towards the true solution.
The iterative refinement process involves updating the current approximation to the solution based on the residual error of the linear system. The solution gradually converges towards the true solution by improving the approximation using information from the Krylov subspace.

The choice of initial vector and the size of the Krylov subspace can influence the convergence and efficiency of Krylov subspace methods. Additionally, the selection of the orthogonalization method and other algorithmic parameters can impact the performance and robustness of these iterative solvers.

\begin{figure}[!h]
    \centering
    \includegraphics[scale=0.6]{images/quadratic-form-graph.png}
    \caption{The graph of quadratic form $f(x)$. The solution of the $Ax = b$ gives the minimum point to the surface}
    \label{fig:quadraticform}
\end{figure}

\subsubsection{The Conjugate Gradient method}
Iterative methods play a crucial role in efficiently solving large systems of linear equations. One such iterative method is the \cg{} (CG) algorithm, which is widely used due to its desirable convergence properties and low memory requirements. It was first introduced by Hestenes and Stiefel in 1952. It has since become a widely used method in numerical linear algebra due to its efficiency and effectiveness on large and sparse matrices. The \cg{} method iteratively finds the solution by minimizing the residual error through a sequence of conjugate directions. However, solving large-scale systems using \cg{} can be computationally demanding and time-consuming.


The method works by minimizing a quadratic form (Figure \ref{fig:quadraticform}) defined by the linear equation system.
Instead of explicitly forming and manipulating the matrix $A$, the CG method operates within the Krylov subspace to iteratively refine the solution. It starts with an initial guess $x_0$ and iteratively updates the solution $x$ by moving in the conjugate direction $p$ within the Krylov subspace. At each iteration, the process computes a search direction conjugate to the previous search direction and uses this direction to update the solution. The method converges in at most n steps, where n is the dimension of the matrix if no rounding errors occur. 

The convex optimization problem in quadratic form can be written as:
\begin{equation}
f(x) = 1/2x^{T}Ax - b^{T}x + c.
\end{equation}

Here $A$ is matrix $x$, and $b$ are vectors, and $c$ is scalar. If we consider that A is symmetric and positive, define $f(x)$ can be minimized by the solution (Figure \ref{fig:quadraticform}):
\begin{equation}
Ax=b.  
\end{equation}

In the algorithm \ref{alg:cg}, $A$ is the coefficient matrix, $b$ is the right-hand side vector, $x$ is the $n$-th approximate solution, $r$ is the $n$-th residual, $d$ is the $n$-th search direction considering we have iteration $n$, and $alpha$ and $beta$ are scalar coefficients. The algorithm iteratively updates the approximate solution $x$ by finding the optimal search direction $d$ and the optimal step size $\alpha$ and then uses them to compute the new approximation $x$. This process continues until the residual $r$ becomes sufficiently small, indicating that a good approximation has been found.\cite{fletcher1964function,lanczos1952solution}.

\begin{algorithm}
\caption{ Conjugate Gradient algorithm}\label{alg:cg}
\begin{algorithmic}
\Ensure $A$ is a symmetric positive-define matrix
\State $r \gets b - Ax$
\State $d \gets r$
\State $\delta \gets r^{T}r$
\While{not converged}
\State $q \gets Ad$
\State $\alpha \gets \dfrac {\delta}{d^{T}q}$
\State $x \gets x + \alpha d$
\State $r \gets r - \alpha q$
\State $\delta_{old} \gets \delta$
\State $\delta \gets r^{T}r$
\State $\beta \gets \dfrac{\delta}{\delta_{old}}$
\State $d \gets r + \beta d$
\EndWhile
\State \Return x
\end{algorithmic}
\end{algorithm}

The Conjugate Gradient (CG) method theoretically converges to the exact solution of the linear system in at most $n$ iterations for symmetric positive definite matrices. However, in practice, the number of iterations can be influenced by several factors: 

\begin{itemize}
  \item \textbf{Numerical Errors}: In real-world computations, numerical errors can accumulate during the iterative process, leading to deviations from the theoretical behavior. These errors can arise from finite precision arithmetic, round-off errors, and truncation errors, among others.
  
  \item \textbf{Ill-Conditioned Matrices}: Despite being designed for symmetric positive definite matrices, ill-conditioned matrices can pose challenges to convergence. In such cases, the number of iterations required for convergence may increase, and the convergence may not be as quick as for well-conditioned matrices.
  
  \item \textbf{Preconditioning Quality}: The effectiveness of preconditioning significantly affects the convergence behavior of iterative solvers like CG. If the preconditioner is not sufficiently effective in reducing the condition number of the matrix, the number of iterations needed for convergence may increase.
  
  \item \textbf{Convergence Criteria}: The choice of convergence criteria also impacts the number of iterations. If the convergence tolerance is set too tight, the iterative solver may require more iterations to meet the criteria, especially if the solution is close to the convergence threshold.
  
  \item \textbf{Matrix Structure and Eigenvalue Distribution}: The distribution of eigenvalues of the matrix $A$ can influence the convergence behavior of the CG method. Clustering of eigenvalues near zero can lead to slow convergence, while widely spread eigenvalues can accelerate convergence.
\end{itemize}


Preconditioners can be applied to accelerate the algorithm's convergence, transforming the original matrix into a new matrix with a lower condition number. The condition number of a matrix measures how sensitive its solution is to perturbations in its coefficients or right-hand side. A matrix with a high condition number is more sensitive to perturbations, and thus, it is more difficult to solve accurately.

The Conjugate Gradient method is efficient in solving large and sparse systems of equations as it only needs the matrix-vector product and can be computed in parallel using OpenCL. This makes the method suitable for scientific and engineering applications such as finite element analysis, image processing, and machine learning.

The \cg{} has been extensively studied and optimized over the years, and many variations and improvements have been proposed. These include preconditioned \cg{} methods, which use a preconditioner to improve the convergence rate of the algorithm, and restarted \cg{} methods, which periodically restart the algorithm with a new initial guess to enhance convergence of difficult problems. 

\subsubsection{\gmres{}}
One of the successful Krylov methods designed to solve large nonsymmetric sparse linear systems of equations of the form $Ax=b$ is \gmres{} (GMRES) method \cite{walker1994simpler}. GMRES iteratively constructs an orthogonal basis for the Krylov subspace generated by the matrix $A$ and the initial residual vector $r_0 = b - Ax_0$, where $x_0$ is the initial guess. GMRES iteratively improves the approximation to the solution by minimizing the residual error in this subspace. The number of iterations required for convergence depends on factors such as the condition number of the matrix $A$, the choice of preconditioner, and the convergence criterion.

Orthogonalization is a key aspect of GMRES. It ensures that search directions are linearly independent. GMRES makes new search directions orthogonal to old ones. It does this during the iterative process to keep orthogonality and numerical stability. GMRES commonly uses Gram-Schmidt orthogonalization. It is used to make the basis vectors of the Krylov subspace orthogonal. However, this process can lead to loss of orthogonality. This is especially true for ill-conditioned matrices. In Gram-Schmidt orthogonalization \cite{balabanov2022randomized}, the method stores the residual vector at each iteration. This allows for cheaper re-initialization after each cycle of the restarted method GMRES(m). This is especially helpful for the restarted method GMRES(m). In this method, the process iterates. It is periodically restarted with a new guess. We need to re-initialize because we lose orthogonality. This loss can happen over long iterations. To fix the loss of orthogonality and save memory, GMRES often uses restarting. The iterative process is periodically restarted with a new initial guess.

One of the disadvantages of the GMRES method is that it might not perform well when applied to linear systems of ill-posed problems. Ill-posed problems contain error-contaminated data in the right-hand side vector \cite{morikuni2014fgmres}. The problem is called an ill-posed problem if one of the following criteria is not met: the existence of the solution, the uniqueness of the solution, or the stability of the solution. Flexible GMRES performs better in such problems by allowing dynamic preconditioning. Another disadvantage of the GMRES method is memory usage. As all the previous residual vectors are stored to maintain the orthogonalization, it might be impractical for large systems. Despite some weaknesses in the GMRES method, it is still very well-suited and robust for non-symmetric matrices and for systems where the direct solvers are not effective.

\subsubsection{\fgmres{}}
Flexible GMRES (FGMRES), also introduced by Yousef Saad in 1993, is an extension of the GMRES (Generalized Minimal Residual) method.
It is designed to handle the cases where the preconditioner may vary during the iterations. This characteristic of the FGMRES method makes it suitable for problems where a fixed preconditioner cannot be easily maintained. The algorithm follows the same processes as GMRES. The difference is that, like GMRES, FGMRES builds a Krylov subspace incrementally but requires additional storage for the changing intermediate-step residuals. FGMRES stores each new search vector and the vector resulting from applying a preconditioner. This allows the solution to be restored to the original unpreconditioned system at the end of the algorithm. Because of high memory requirements, restarting the algorithm might be required for large-scale problems.

\subsection{Domain Decomposition Methods (DDM)}
Choosing between a direct solver or an iterative solver for large problems is not trivial. Iterative solvers are regarded as highly memory efficient, especially for large-scale problems; however, it is observed that the choice of an iterative solver is highly specific to the problem. Iterative solvers do not work very well for problems with highly ill-conditioned matrices. On the other hand, direct solvers are considered very robust and independent of the problem; however, they are seriously limited due to the large memory requirement \cite{raju2009domain}. In such cases, domain decomposition methods come into play to improve the following properties:
\begin{itemize}
  \item \textbf{Parallel Computation}: Domain decomposition methods divide a large computational domain into smaller subdomains. This allows different processors to handle separate subdomains concurrently, reducing computation time significantly.
  \item \textbf{Memory Management}:
  Large problems typically require substantial memory, which might exceed the capacity of a single processor or node. By dividing the domain, each processor handles a smaller part, making the memory requirements more manageable.
  \item \textbf{Scalability}: These methods enhance the scalability of algorithms to large systems, allowing them to efficiently utilize the increasing number of processors in modern computing environments.
  \item \textbf{Convergence Rate Improvement}: In iterative methods, domain decomposition can be used as a preconditioner to improve the convergence rates of iterative solvers.
\end{itemize}

Domain Decomposition Methods are a class of numerical techniques used to solve PDEs on large and complex computational domains. The idea behind these methods is to divide the computational domain into smaller regions called subdomains. Each subdomain is then solved separately, and the solutions from each subdomain are combined to obtain the solution for the entire domain. Doing this would enable the preconditioner to improve the properties of the matrix system globally, in particular, eigenvalue distribution, and hence contribute toward speeding up the convergence of the global iterative method. The most common types of domain decomposition methods are the Schwarz method, the substructuring method, and the mortar method.

In the context of Krylov subspace methods, domain decomposition methods can be used to improve the efficiency of the solver for large linear systems arising from PDEs. The GMRES, FGMRES, and CG mentioned above are iterative methods that require the solution of a linear system at each iteration. Domain decomposition methods can be used to parallelize the solution of the linear system, thereby reducing the time required to solve the problem.

There are two main types of Domain Decomposition methods: overlapping and non-overlapping \cite{toselli2004domain}. The non-overlapping domain decomposition sees a large computational domain divided into a number of subdomains. All the same, the subdivisions share the boundaries but don't overlap. Interactions between the subdomains are so important that the interfaces require special treatment to ensure continuity and transfer of information in the right manner. It gives us high efficiency in parallel computation as subdomains can be solved independently. Also reduces complexity per subdomain, making it manageable to use direct solvers for each. Non-overlapping domain decomposition methods demonstrate slower convergence than overlapping domain decomposition methods.

\begin{figure}[ht]
    \centering
    \includegraphics[scale=0.2]{images/Overlapping DD.png}
    \caption{Overlapping domains}
    \label{fig:overlapping-dd}
\end{figure}

On the other hand, overlapping methods, also known as Schwarz methods, involve subdomains overlapping by some margin\cite{cai2003overlapping}. The overlap region provides a better correction of errors and a smoother transfer of information for the boundaries of each subdomain. Every subdomain has part of its neighboring subdomains. In particular, the width of the overlap significantly impacts the convergence and stability of the method. In non-overlapping methods, the local problems are solved independently from other subdomains. However, the overlapping regions ensure that each subdomain solver would have more accurate information on the boundary, which would further improve the effectiveness of the solvers. We benefit from improved convergence due to a better exchange of information in the overlapping regions and better robustness to problem parameter variations or discretization errors. 

A simple graphical visualization of overlapping domain decomposition methods is shown in Figure \ref{fig:overlapping-dd}. It shows a domain formed by the union of a circle and a rectangle, which is donated as $\Omega = \Omega_1 \cup \Omega_2$
The \as{} (AS), %Generalized \as{} (GAS), 
Restricted \as{} (RAS), and 
\raas{} (RAAS) Schwarz methods are all examples of overlapping domain decomposition methods. The next section will describe RAAS in more detail because of its usage in our practical implementation.


\subsubsection{The \raas{} method}
The standard Additive Schwarz method divides the computational domain into overlapping subdomains with minimal overlap. The minimal overlap consists of a single layer of nodes in the discretization grid. Having a larger extended overlap between subdomains in domain decomposition methods can improve the convergence rate of iterative solvers by reducing the number of iterations needed to reach a solution. In particular, this is known to be true in the case of the Helmholtz problem (see Subsection \ref{helmholtz}) due to the spectral properties of its stiffness matrix obtained through the discretization. The overlap allows for a smoother transition of information between adjacent subdomains, which can help to reduce the influence of interface errors. The minimal overlap is, therefore, extended. Then, each subdomain solves the local problem independently on the extended subdomain. However, the Restricted Additive Schwarz methods modify the Additive Schwarz method by introducing a restriction operator \cite{efstathiou2003restricted}. The restriction operator limits the updates of each subdomain after the subdomain solution to its interior nodes. The Additive Average Schwar (AAS) method updates the values on the overlap with averages from each neighboring subdomain instead of adding up the separate contributions. With the \raas{} (RAAS) method, averaging is done only on the minimal overlap. Outside of the minimal overlap -- the extended overlap values -- are used only as an input for the extended subdomain problem solution. 

This adjustment leads to improved performance in terms of both convergence and stability, particularly when dealing with certain large-scale parallel computations, like the discretization from the Helmholtz equation \cite{graham2017domain}.

Figure \ref{fig:overlappingsubdomains} illustrates the partitioning of the unit square into four overlapping subdomains. The symbol $\Delta$ represents the extent of overlap between these subdomains, while red dashed lines demonstrate the minimal overlaps. Consequently, each extended subdomain will have a dimension of $0.5 + \Delta$. The solver will then be employed within these expanded regions. Upon completion of the subdomain solutions, the \raas{} method will compute an average from the solutions within the minimal overlap areas while discarding the solutions from nodes located in the extended overlaps.

\begin{figure}[!h]
    \centering
    \includegraphics[scale=0.6]{images/overlapping-subdomains.png}
    \caption{An example of partitioning the unit square into four overlapping subdomains with an extended overlap of width $\Delta$.}
    \label{fig:overlappingsubdomains}
\end{figure}