Parallel Algorithm for Solving the Inverse Two-Dimensional Fractional Diffusion Problem of Identifying the Source Term

Abstract

This paper is devoted to the development of a parallel algorithm for solving the inverse problem of identifying the space-dependent source term in the two-dimensional fractional diffusion equation. For solving the inverse problem, the regularized iterative conjugate gradient method is used. At each iteration of the method, we need to solve the auxilliary direct initial-boundary value problem. By using the finite difference scheme, this problem is reduced to solving a large system of a linear algebraic equation with a block-tridiagonal matrix at each time step. Solving the system takes almost the entire computation time. To solve this system, we construct and implement the direct parallel matrix sweep algorithm. We establish stability and correctness for this algorithm. The parallel implementations are developed for the multicore CPU using the OpenMP technology. The numerical experiments are performed to study the performance of parallel implementations.

1. Introduction

Numerous physical phenomena exhibit the characteristics of memory retention and nonlocal effects [1]. The formulation of mathematical frameworks to describe these phenomena often involves the utilization of fractional calculus [2]. Various fractional derivative operators possess distinct definitions and inherent properties. The fractional differential equation may be applied to a wide range of fields such as anomalous diffusion [3,4,5], viscoelasticity [6], ferroelectric media [7], fractional multi-pole and neuron modelling [8], and fractional Lèvy motion [9]. The review paper [10] presents a comprehensive survey on real-world applications of fractional calculus in various fields, namely, physics; control, signal and image processing; mechanics and dynamic systems; biology; environmental science; material studies; economics; and engineering.

The solution to direct and inverse problems for differential equations with fractional derivatives typically incurs substantial computational costs due to their nonlocal characteristics. Different numerical techniques are available for solving approximate initial-boundary problems for fractional differential equations [11,12,13,14], for example, the finite difference method. An established approach to enhance computational efficiency involves the utilization of parallel computing [15,16,17]. The problems for fractional differential equations may be solved using various original parallel algorithms [18,19].

Additionally, while forward problems are usually well developed, the inverse problems may present significant challenges to achieve the stable solutions [20]. Here, by direct problems, we mean the classical initial boundary problem of finding the unknown function from equation and additional boundary and initial conditions. Inverse problems include identifying unknown coefficients of an equation or unknown boundary or initial conditions [21]. In this case, a priori information is used to ensure the uniqueness of the solution. Such tasks are often incorrect. Regularization methods are used to achieve correctness. Hence, the development and implementation of parallel algorithms aimed at solving inverse problems emerging from fractional differential equations is of great importance.

In [22], an algorithm for solving the inverse problem of identifying the space-dependent source is constructed on the base of regularized Landweber iterative method. An iterative algorithm on the basis of the conjugate gradient method is constructed in [23]. For smoothing the values, the authors used the Savitzky–Golay filter.

In [24], existence, uniqueness, and stability estimates are established for the inverse source problem of the recovery of a space-dependent source term for a generalized subdiffusion equation. The Tikhonov regularization method was proposed for solving the problem of reconstruction of a time-dependent source term in a time-fractional diffusion-wave equation [25]. Convergence estimates and parameted choice rules for Tikhonov regularization for the inverse source problem for a time fractional diffusion equation were proposed in [26].

In our previous works [27,28] parallel algorithms for solving direct and inverse problems for one-dimensional time-fractional diffusion equation are constructed on the basis of the parallel sweep method for solving the systems of linear algebraic equations with tridiagonal matrices.

The goal of our work is to construct and implement an efficient algorithm for solving the inverse problem of identifying the stationary source term for the two-dimensional time-fractional diffusion equation. To solve the inverse problem, we apply the regularized iterative conjugate gradient method. It requires solving the auxilliary direct subproblem at each iteration. By using the finite difference scheme, we reduce the direct problem to solving a series of systems of linear algebraic equations (SLAE) with block tridiagonal matrices. Note that solving these SLAEs for the two-dimensional fractional diffusion problem takes most of the computation time (up to 99% of the total computational time). We construct and establish the stability and correctness of the direct parallel matrix sweep method for solving such systems. We have developed a parallel algorithm and a parallel code for solving the inverse problem. The code is intended for multicore processors and is implemented using C++17 and OpenMP 4.0 extension. In order to evaluate the validity of developed numerical methods and the performance of parallel algorithms, numerical experiments are conducted. In the future, we plan to use our algorithm and code for solving applied problems.

Below is a breakdown of the article’s structure. In Section 2, we formulate the direct and inverse problems and introduce the discretization and difference scheme for solving the direct initial-boundary value problem. We demonstrate the block-tridiagonal structure of the coefficient matrix of the SLAE. In Section 3, we describe the direct methods for solving SLAEs with the block-tridiagonal matrix. We establish the stability and correctness of the parallel matrix sweep algorithm. In Section 4, we describe the conjugate gradient algorithm for solving the inverse problem. Section 5 describes the development of the parallel code that implements the numerical algorithms described above. Section 6 presents the results of the performed numerical experiments. In Section 7, we discuss these results. Section 8 concludes our work.

2. Problem

2.1. Statement of the Problem

Consider the basis time-fractional elliptic partial differential equation in the following form:

$$\frac{{\partial }^{\alpha }U}{\partial t^{\alpha }}+\mathcal{L}U=\psi \left(\mathbf{x}\right)\eta (\mathbf{x},t),\mathbf{x}\in \Omega ,0<t\le T.$$

(1)

Here, $\mathbf{x}=(x_1,x_2,\dots ,x_d)\in \overline{\Omega }=\prod _{i=1}^d[0,{\gamma }_i],0<\alpha <1,q(x,t)\ne 0$ and $\mathcal{L}$ is an elliptic operator

$$\mathcal{L}U=-\sum _{i=1}^d\frac{\partial }{\partial x_i}\left(k_i(\mathbf{x},t)\frac{\partial U}{\partial x_i}\right),x_i\in (0,{\gamma }_i),0<t\le T.$$

The boundary condition is

$$U(\mathbf{x},t)=0,\mathbf{x}\in \delta \Omega .$$

(2)

The initial condition is

$$U(\mathbf{x},0)=g_0\left(\mathbf{x}\right),\mathbf{x}\in \Omega .$$

(3)

where $g_0\left(\mathbf{x}\right)$ is the given functions.

For simplicity, we have formulated our fractional diffusion equation with homogeneous Dirichlet boundary conditions. Note that the problem with inhomogeneous boundary conditions may be reduced to a problem with homogeneous conditions. To do this, we need to represent function U as the sum of a some function $U_0$ that satisfies the inhomogeneous boundary conditions (for example, found by solving the Dirichlet problem for the Laplace’s equation) plus a remainder function V that will satisfy the homogeneous conditions $U(\mathbf{x},t)=U_0(\mathbf{x},t)+V(\mathbf{x},t)$. The algorithms and approaches presented below may be utilized for more general formulations with the Neumann or mixed boundary conditions.

For this study, we consider the Caputo fractional derivative with order $\alpha$ in the form [29]

$$D^{\alpha }f\left(x\right)=\frac{1}{\Gamma (m-\alpha )}\underset{0}{\overset{x}{\int }}\frac{f^{\left(m\right)}\left(t\right)}{{(x-t)}^{\alpha -m+1}}dt,$$

with $\alpha \in (m-1,m)$, $m\in N,x>0$.

Assuming that the solution $U(\mathbf{x},t)$ exists and satisfies the Dirichlet boundary conditions, for the case of $0<\alpha <1$ ($m=1$) we can consider the following formula for the Caputo fractional partial derivative:

$$\frac{{\partial }^{\alpha }U(\mathbf{x},t)}{\partial t^{\alpha }}=\frac{1}{\Gamma (1-\alpha )}\underset{0}{\overset{t}{\int }}\frac{\partial U(\mathbf{x},s)}{\partial s}{(t-s)}^{-\alpha }ds.$$

(4)

The direct initial boundary problem consists in finding the unknown function $U(\mathbf{x},t)$ when all other components of Equation (1) are known.

In the present work, we study the inverse problem of restoring the space-dependent right-hand part $\psi \left(\mathbf{x}\right)$. Thus, the problem consists in finding the pair of unknown functions $\left[U\right(\mathbf{x},t),\psi (\mathbf{x}\left)\right]$. Additional information for the inverse problem is given in the form of final overdetermination

$$U(\mathbf{x},T)=\phi \left(\mathbf{x}\right),\mathbf{x}\in \Omega .$$

(5)

Conditions for the uniqueness of the solution of this inverse problem for a general multi-dimensional equation are formulated in [30].

2.2. Discretization of Equation and Difference Scheme

In this paper, for simplicity, we will consider the case of $d=2$ and $k_1=k_2=1$. For an arbitrary elliptic operator, the general structure of the coefficient matrix of the SLAE remains the same, and the methods and algorithms presented below will be applicable.

Equation (1) will take the following form:

$$\frac{{\partial }^{\alpha }U(x_1,x_2,t)}{\partial t^{\alpha }}=\frac{{\partial }^{2}U(x_1,x_2,t)}{\partial x_1^2}+\frac{{\partial }^{2}U(x_1,x_2,t)}{\partial x_2^2}+\psi (x_1,x_2)\eta (x_1,x_2,t),$$

(6)

where $U(x_1,x_2,t)$ and $\psi (x_1,x_2)$ are the sought functions; $a(x_1,x_2),b(x_1,x_2),c(x_1,x_2)$, $\eta (x_1,x_2,t)$ are the given functions; and $0<\alpha <1$ is the the order of the fractional derivative.

The problem is considered for the area $\Omega :0\le x_1\le {\gamma }_1,0\le x_2\le {\gamma }_2$ and time interval $0\le t\le T$.

To discretize Equation (6), we construct the partitioning on the solution domain. On intervals $[0,{\gamma }_1],[0,{\gamma }_2]$, we introduce the grid with n and N points. The steps are $h_1=\Delta x_1={\gamma }_1/n$ and $h_2=\Delta x_2={\gamma }_2/N$. For the time interval $[0,T]$, we use the uniform grid of $\tilde{N}$ points. The step is $\tau =\Delta t=T/\tilde{N}$. We denote the nodes as $x_{1,i_1}=i_1{h}_1,x_{2,i_2}=i_2{h}_2,i_1\in \left\{0,1,\dots ,n\right\},i_2\in \left\{0,1,\dots ,N\right\}$, and $t_j=j\tau ,j\in \left\{0,1,\dots ,\tilde{N}\right\}$. The values of discretized functions U and others are denoted as $U_{i_1,i_2,j}=U(x_{1,i_1},x_{2,i_2},t_j)$.

In this work, the Grunwald–Letnikov formula [2] is used for approximating the Caputo fractional derivative in Equation (6)

$$\begin{array}{cc}\hfill & D_t^{\alpha }{U}_{i_1,i_2,j}\cong {\sigma }_{\alpha ,\tau }\sum _{\mathsf{\ell }=1}^j{w}_{\mathsf{\ell }}^{\left(\alpha \right)}(U_{i_1,i_2,j-\mathsf{\ell }+1}-U_{i_1,i_2,j-\mathsf{\ell }}),\hfill \\ \hfill & {\sigma }_{\alpha ,\tau }=\frac{1}{\Gamma (1-\alpha )(1-\alpha ){\tau }^{\alpha }},w_{\mathsf{\ell }}^{\left(\alpha \right)}={\mathsf{\ell }}^{1-\alpha }-{(\mathsf{\ell }-1)}^{1-\alpha }.\hfill \end{array}$$

(7)

The implicit two-step finite difference scheme is used for approximating Equation (6). For the grid point $(x_{i_1},x_{i_2})$ at time layer $t_j$, the difference equation has the form

$$\begin{array}{cc}\hfill & {\sigma }_{\alpha ,\tau }\sum _{\mathsf{\ell }=1}^j{w}_{\mathsf{\ell }}^{\left(\alpha \right)}(U_{i_1,i_2,j-\mathsf{\ell }+1}-U_{i_1,i_2,j-\mathsf{\ell }})=\hfill \\ \hfill & =\frac{U_{i_1-1,i_2,j}-2U_{i_1,i_2,j}+U_{i_1+1,i_2,j}}{h_1^2}+\frac{U_{i_1,i_2-1,j}-2U_{i_1,i_2,j}+U_{i_1,i_2+1,j}}{h_2^2}+{\psi }_{i_1,i_2}{\eta }_{i_1,i_2,j}.\hfill \end{array}$$

(8)

2.3. Constructing the SLAE

Let us apply the following transformations:

$$\begin{array}{cc}\hfill & {\sigma }_{\alpha ,\tau }(U_{i_1,i_2,j}-U_{i_1,i_2,j-1})+{\sigma }_{\alpha ,\tau }\sum _{\mathsf{\ell }=2}^j{w}_{\mathsf{\ell }}^{\left(\alpha \right)}(U_{i_1,i_2,j-\mathsf{\ell }+1}-U_{i_1,i_2,j-\mathsf{\ell }})=\hfill \\ \hfill & =r_{i_1,i_2}{U}_{i_1,i_2-1,j}+p_{i_1,i_2}{U}_{i_1-1,i_2,j}+(-2p_{i_1,i_2}-2q_{i_1,i_2})U_{i_1,i_2,j}+\hfill \\ \hfill & +p_{i_1,i_2}{U}_{i_1+1,i_2,j}+r_{i_1,i_2}{U}_{i_1,i_2+1,j}+{\psi }_{i_1,i_2}{\eta }_{i_1,i_2,j};\hfill \end{array}$$

$$\begin{array}{cc}\hfill & -r_{i_1,i_2}{U}_{i_1,i_2-1,j}-p_{i_1,i_2}{U}_{i_1-1,i_2,j}+q_{i_1,i_2}{U}_{i_1,i_2,j}-p_{i_1,i_2}{U}_{i_1+1,i_2,j}-r_{i_1,i_2}{U}_{i_1,i_2+1,j}=\hfill \\ \hfill & ={\sigma }_{\alpha ,\tau }\left(U_{i_1,i_2,j-\mathsf{\ell }}-\sum _{\mathsf{\ell }=2}^j{w}_{\mathsf{\ell }}^{\left(\alpha \right)}(U_{i_1,i_2,j-\mathsf{\ell }+1}-U_{i_1,i_2,j-\mathsf{\ell }})\right)+{\psi }_{i_1,i_2}{\eta }_{i_1,i_2,j},\hfill \end{array}$$

where

$$p_{i_1,i_2}=\frac{1}{h_1^2},r_{i_1,i_2}=\frac{1}{h_2^2},q_{i_1,i_2}={\sigma }_{\alpha ,\tau }+2p_{i_1,i_2}+2q_{i_1,i_2}.$$

Let us denote

$$\begin{array}{cc}\hfill & f_{i_1,i_2,j}={\sigma }_{\alpha ,\tau }\left(U_{i_1,i_2,j-\mathsf{\ell }}-\sum _{\mathsf{\ell }=2}^j{w}_{\mathsf{\ell }}^{\left(\alpha \right)}(U_{i_1,i_2,j-\mathsf{\ell }+1}-U_{i_1,i_2,j-\mathsf{\ell }})\right)+{\psi }_{i_1,i_2}{\eta }_{i_1,i_2,j},j>1,\hfill \\ \hfill & f_{i_1,i_2,1}={\sigma }_{\alpha ,\tau }{U}_{i_1,i_2,0}+{\psi }_{i_1,i_2}{\eta }_{i_1,i_2,0}.\hfill \end{array}$$

(9)

The difference equation will take the form

$$-r_{i_1,i_2}{U}_{i_1,i_2-1,j}-p_{i_1,i_2}{U}_{i_1-1,i_2,j}+q_{i_1,i_2}{U}_{i_1,i_2,j}-p_{i_1,i_2}{U}_{i_1+1,i_2,j}-r_{i_1,i_2}{U}_{i_1,i_2+1,j}=f_{i_1,i_2,j}.$$

(10)

Note that the values $U_{0,i_2,j},U_{n,i_2,j},i_2\in \left\{0,1,\dots ,N\right\}$ and $U_{i_1,0,j},U_{i_1,N,j}$, $i_1\in \left\{0,1,\dots ,n\right\}$ at the boundaries are given. Then, for all inner points $(x_{i_1},x_{i_2})$, $i_1\in \left\{1,2,\dots ,n-1\right\}$,$i_2\in \left\{1,2,\dots ,N-1\right\}$, difference Equation (10) constitutes an SLAE

$$AY=F,$$

(11)

where A is the block-tridiagonal matrix of dimension $(n-2)(N-2)\times (n-2)(N-2)$

$$\begin{array}{cc}\hfill & A=\left[\begin{array}{ccccccc}{C}_1& -B_1& & & & & \\ -A_2& C_2& -B_2& & \\ & \ddots & \ddots & \ddots & & \\ & & -A_{i_2}& C_{i_2}& -B_{i_2}\\ & & & \ddots & \ddots & \ddots & & \\ & & & & -A_{N-2}& C_{N-2}& -B_{N-2}\\ & & & & & -A_{N-1}& C_{N-1}\end{array}\right],\hfill \\ \hfill & i_2\in \left\{1,2,\dots ,N-1\right\}.\hfill \end{array}$$

Blocks $A_{i_2},B_{i_2},C_{i_2}$ of dimension $(n-2)\times (n-2)$ are defined as

$$\begin{array}{cc}\hfill & A_{i_2}=B_{i_2}=\left[\begin{array}{ccccc}{r}_{1,i_2}& & & & \\ & r_{2,i_2}& & & \\ & & \ddots & & & \\ & & & r_{i_1,i_2}& \\ & & & & \ddots & & & \\ & & & & & r_{n-2,i_2}& \\ & & & & & & r_{n-1,i_2}\end{array}\right],\hfill \\ \hfill & i_1\in \left\{1,2,\dots ,n-1\right\},\hfill \end{array}$$

$$\begin{array}{cc}\hfill & C_{i_2}=\left[\begin{array}{ccccc}{q}_{1,i_2}& -p_{1,i_2}& & & \\ -p_{2,i_2}& q_{2,i_2}& -p_{2,i_2}& & \\ & \ddots & \ddots & \ddots & & \\ & & -p_{i_1,i_2}& q_{i_1,i_2}& -p_{i_1,i_2}\\ & & & \ddots & \ddots & \ddots & & \\ & & & & -p_{n-2,i_2}& q_{n-2,i_2}& -p_{n-2,i_2}\\ & & & & & -p_{n-1,i_2}& q_{n-1,i_2}\end{array}\right],\hfill \\ \hfill & i_1\in \left\{1,2,\dots ,n-1\right\}.\hfill \end{array}$$

The sought vector consists of values $U_{i_1,i_2,j}$ collapsed in a row-major order

$$\begin{array}{cc}\hfill & Y={\left[Y_1,Y_2,...,Y_{N-1}\right]}^{⊺},\hfill \\ \hfill & Y_1=\left[U_{1,1,j},U_{2,1,j},\dots ,U_{n-2,1,j},U_{n-1,1,j}\right];\hfill \\ \hfill & Y_2=\left[U_{1,2,j},U_{2,2,j},\dots ,U_{n-2,2,j},U_{n-1,2,j}\right];\hfill \\ \hfill & \dots \hfill \\ \hfill & Y_{N-1}=\left[U_{1,N-1,j},U_{2,N-1,j},\dots ,U_{n-1,N-1,j}\right].\hfill \end{array}$$

Right-hand vector is constructed similarly, taking into account the boundary points

$$\begin{array}{cc}\hfill F=[& F_1,F_2,\dots ,F_{N-1}{]}^{⊺};\hfill \\ \hfill F_1=[& F_{1,1,j}+r_{1,1}{U}_{1,0,j}+p_{1,1}{U}_{0,1,j},F_{2,1,j}+r_{2,1}{U}_{2,0,j},\dots ,\hfill \\ \hfill & F_{n-2,1,j}+r_{n-2,1}{U}_{n-2,0,j},F_{n-1,1,j}+r_{n-1,0}{U}_{n-2,0,j}+p_{n-1,1}{U}_{n,1,j}];\hfill \\ \hfill F_2=[& F_{1,2,j}+p_{1,2}{U}_{0,2,j},F_{2,2,j},\dots ,F_{n-1,2,j},F_{n-1,2,j}+p_{n-1,2}{U}_{n,2,j}];\hfill \\ \hfill \dots & \\ \hfill F_{N-1}=[& F_{1,N-1,j}+r_{1,N-1}{U}_{1,N,j}+p_{1,N-1}{U}_{0,N-1,j},F_{2,N-1,j}+r_{2,N-1}{U}_{2,N,j},\dots ,\hfill \\ \hfill & F_{n-1,N-1,j}+r_{n-1,N-1}{U}_{n-1,N,j}+p_{n-1,N-1}{U}_{n,N-1,j}].\hfill \end{array}$$

To numerically solve the initial boundary problem, we need to solve system (11) at each subsequent time level $j=1,...,\tilde{N}$.

3. Numerical Methods for Solving the SLAE

To solve the SLAEs with various matrix structures, different numerical methods may be used. In this work, for solving the block-tridiagonal SLAE (11) we will use the block-elimination method [31] and parallel matrix sweep method [32].

3.1. Block-Elimination Method

Let us consider an SLAE with a block tridiagonal matrix in the following form:

$$\left\{\begin{array}{cccc}\hfill C_0{\overline{Y}}_0-B_0{\overline{Y}}_1& =& {\overline{F}}_0,& i=0,\hfill \\ \hfill -A_i{\overline{Y}}_{i-1}+C_i{\overline{Y}}_i-B_i{\overline{Y}}_{i+1}& =& {\overline{F}}_i,& i=1,2,\cdots ,N-1,\hfill \\ \hfill -A_N{\overline{Y}}_{N-1}+C_N{\overline{Y}}_N& =& {\overline{F}}_N,& i=N,\hfill \end{array}\right.$$

(12)

where ${\overline{Y}}_i$ are the sought vectors of the n dimension; ${\overline{F}}_i$ are the given right-hand vectors of the n dimension; and $A_i,B_i,C_i$ are the square matrices of the $n\times n$ dimension.

The direct block elimination (or matrix sweep) method is intended for solving the SLAE with block tridiagonal matrix. The auxilliary coefficients ${\alpha }_i$ (matrices of $n\times n$ dimension) and ${\beta }_i$ (vectors of n dimension) are found by the reccurent formulae (the forward elimination phase)

$$\begin{array}{cc}\hfill & {\alpha }_0=C_0^{-1}{B}_0,{\alpha }_i={(C_i-A_i{\alpha }_i)}^{-1}{B}_i,i=1,2,\dots ,N-1,\hfill \\ \hfill & {\beta }_0=C_0^{-1}{F}_0,{\beta }_i={(C_i-A_i{\alpha }_i)}^{-1}(\overline{F_i}+A_i{\beta }_{i-1}),i=1,2,\dots ,N.\hfill \end{array}$$

(13)

Solution vectors $\overline{Y_i}$ are found by formulae (the backward substitution phase)

$$\begin{array}{c}\hfill \overline{Y_N}={\beta }_N,\overline{Y_i}={\alpha }_i\overline{Y_{i+1}}+{\beta }_i,i=N-1,N-2,\dots ,1,0.\end{array}$$

(14)

Remark 1.

The algorithms (13) and (14) are correct if matrices $C_0$ and $(C_i-A_i{\alpha }_i)$ are nonsingular for $i=1,2,...,N$. The algorithm is stable if $∥{\alpha }_i∥\le 1$ for $i=1,2,...,N$.

The stability condition for this algorithm is the following.

Lemma 1.

If matrices $C_i,i=0,1,...,N$ are nonsingular, matrices $A_i,B_i,i=1,...,N-1$ are non-null, and conditions

$$\parallel C_0^{-1}{B}_0\parallel \le 1,\parallel C_N^{-1}{A}_N\parallel \le 1,\parallel C_i^{-1}{A}_i\parallel +\parallel C_i^{-1}{B}_i\parallel \le 1,i=1,\cdots ,N-1,$$

(15)

are satisfied where at least one of the inequalities is strict; then, the algorithms (13) and (14) are stable and correct.

Remark 2.

The coefficients ${\alpha }_i,{\beta }_i$ are dependent on the previous ones ${\alpha }_{i-1},{\beta }_{i-1}$. Thus, we cannot distribute these calculations to independent workers. Thus, the parallelization is limited to operations of matrix inversion and multiplication.

3.2. Parallel Matrix Sweep Method

To construct a direct parallel algorithm, let us split the interval $i=0,1,...,N$ into L subintervals of length M such as $N=L\times M$. Consider the unknown values ${\overline{Y}}_K,K=0,M,...,N$ as parameters.

Now, let us construct the reduced SLAE for ${\overline{Y}}_K$. To do this, consider the following problems for the interval $(K,K+M)$

$$\left\{\begin{array}{ccc}-A_i{\overline{U}}_{i-1}^1+C_i{\overline{U}}_i^1-B_i{\overline{U}}_{i+1}^1=0,& {\overline{U}}_K^1=(10...0),& {\overline{U}}_{K+M}^1=(00...0),\\ \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots & \dots \dots \dots \dots & \dots \dots \dots \dots .\\ \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots & \dots \dots \dots \dots & \dots \dots \dots \dots .\\ -A_i{\overline{U}}_{i-1}^n+C_i{\overline{U}}_i^n-B_i{\overline{U}}_{i+1}^n=0,& {\overline{U}}_K^n=(00...1),& {\overline{U}}_{K+M}^n=(00...0),\end{array}\right.$$

(16)

$$\left\{\begin{array}{ccc}-A_i{\overline{V}}_{i-1}^1+C_i{\overline{V}}_i^1-B_i{\overline{V}}_{i+1}^1=0,& {\overline{V}}_K^1=(00...0),& {\overline{V}}_{K+M}^1=(10...0),\\ \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots & \dots \dots \dots \dots & \dots \dots \dots \dots .\\ \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots & \dots \dots \dots \dots & \dots \dots \dots \dots .\\ -A_i{\overline{V}}_{i-1}^n+C_i{\overline{V}}_i^n-B_i{\overline{V}}_{i+1}^n=0,& {\overline{V}}_K^n=(00...0),& {\overline{V}}_{K+M}^n=(00...1),\end{array}\right.$$

(17)

$$\left\{\begin{array}{ccc}-A_i{\overline{W}}_{i-1}+C_i{\overline{W}}_i-B_i{\overline{W}}_{i+1}={\overline{F}}_i,& {\overline{W}}_K=(00...0),& {\overline{W}}_{K+M}=(00...0),\end{array}\right.$$

(18)

where $i=K+1,\cdots ,K+M-1$.

Theorem 1.

If ${\overline{U}}_i^1,\cdots ,{\overline{U}}_i^n$ are solutions of auxilliary problem (16); ${\overline{V}}_i^1,\cdots ,{\overline{V}}_i^n$ are solutions of problem (17); ${\overline{W}}_i$ are solutions of problem (18); and ${\overline{Y}}_i$ are solutions of the basic problem (12) for interval $(K,K+M)$, then

$${\overline{Y}}_i=({\overline{U}}_i^1{\overline{U}}_i^2\cdots {\overline{U}}_i^n){\overline{Y}}_K+({\overline{V}}_i^1{\overline{V}}_i^2\cdots {\overline{V}}_i^n){\overline{Y}}_{K+M}+{\overline{W}}_i.$$

(19)

Proof. 

Consider system (12) for the inner subinterval $(K,K+M)$:

$$-A_i{\overline{Y}}_{i-1}+C_i{\overline{Y}}_i-B_i{\overline{Y}}_{i+1}={\overline{F}}_i,i=K+1,\cdots ,K+M-1.$$

(20)

This system contains the parameters ${\overline{Y}}_K$ and ${\overline{Y}}_{K+M}$.

Let us rewrite (20) in the form

$$\left\{\begin{array}{c}{C}_{K+1}{\overline{Y}}_{K+1}-B_{K+1}{\overline{Y}}_{K+2}=A_{K+1}{\overline{Y}}_K+{\overline{F}}_{K+1},\hfill \\ -A_i{\overline{Y}}_{i-1}+C_i{\overline{Y}}_i-B_i{\overline{Y}}_{i+1}={\overline{F}}_i,i=K+2,\cdots ,K+M-2,\hfill \\ -A_{K+M-1}{\overline{Y}}_{K+M-2}+C_{K+M-1}{\overline{Y}}_{K+M-1}=B_{K+M-1}{\overline{Y}}_{K+M}+{\overline{F}}_{K+M-1}.\hfill \end{array}\right.$$

(21)

In this system (21), $A_{K+1}{\overline{Y}}_K$ and $B_{K+M-1}{\overline{Y}}_{K+M}$ have the following form:

$$A_{K+1}{\overline{Y}}_K=\left(\begin{array}{ccc}{A}_{K+1}^{11}& ...& A_{K+1}^{1N}\\ ...& ...& ...\\ A_{K+1}^{N1}& ...& A_{K+1}^{NN}\end{array}\right)\left(\begin{array}{c}{Y}_K^1\\ ...\\ Y_K^N\end{array}\right)=$$

$$=\left(\begin{array}{c}{A}_{K+1}^{11}\\ ...\\ A_{K+1}^{N1}\end{array}\right)Y_K^1+\cdots +\left(\begin{array}{c}{A}_{K+1}^{1N}\\ ...\\ A_{K+1}^{NN}\end{array}\right)Y_K^N=$$

$$={\overline{A}}_{K+1}^1{Y}_K^1+{\overline{A}}_{K+1}^2{Y}_K^2+\cdots +{\overline{A}}_{K+1}^N{Y}_K^N;$$

(22)

$$B_{K+M-1}{\overline{Y}}_{K+M}=\left(\begin{array}{ccc}{B}_{K+M-1}^{11}& ...& B_{K+M-1}^{1N}\\ ...& ...& ...\\ B_{K+M-1}^{N1}& ...& B_{K+M-1}^{NN}\end{array}\right)\left(\begin{array}{c}{Y}_{K+M}^1\\ ...\\ Y_{K+M}^N\end{array}\right)=$$

$$=\left(\begin{array}{c}{B}_{K+M-1}^{11}\\ ...\\ B_{K+M-1}^{N1}\end{array}\right)Y_{K+M}^1+\cdots +\left(\begin{array}{c}{B}_{K+M-1}^{1N}\\ ...\\ B_{K+M-1}^{NN}\end{array}\right)Y_{K+M}^N=$$

$$={\overline{B}}_{K+M-1}^1{Y}_{K+M}^1+{\overline{B}}_{K+M-1}^2{Y}_{K+M}^2+\cdots +{\overline{B}}_{K+M-1}^N{Y}_{K+M}^N.$$

Here, the line over the symbol denotes the vector column of corresponding matrix.

Taking into account the formulae (22), SLAEs (21) have the following form:

$$\left(\begin{array}{ccccc}{C}_{K+1}& -B_{K+1}& & & \\ -A_{K+2}& C_{K+2}& -B_{K+2}\\ & \ddots & \ddots & \ddots & \\ & & -A_{K+M-1}& C_{K+M-1}\end{array}\right)\left(\begin{array}{c}{\overline{Y}}_{K+1}\\ {\overline{Y}}_{K+2}\\ \dots \\ {\overline{Y}}_{K+M-1}\end{array}\right)=$$

$$=\left(\begin{array}{c}{\overline{F}}_{K+1}\\ {\overline{F}}_{K+2}\\ ...\\ {\overline{F}}_{K+M-1}\end{array}\right)+\left(\begin{array}{c}{\overline{A}}_{K+1}^1\\ 0\\ ...\\ 0\end{array}\right)Y_K^1+\cdots +\left(\begin{array}{c}{\overline{A}}_{K+1}^N\\ 0\\ ...\\ 0\end{array}\right)Y_K^N+$$

(23)

$$+\left(\begin{array}{c}0\\ ...\\ 0\\ {\overline{B}}_{K+M-1}^1\end{array}\right)Y_{K+M}^1+\cdots +\left(\begin{array}{c}0\\ ...\\ 0\\ {\overline{B}}_{K+M-1}^N\end{array}\right)Y_{K+M}^N.$$

System (23) is equivalent to the following one:

$$\Lambda \overline{Y}=\overline{F}+{\overline{A}}^1{Y}_K^1+\cdots +{\overline{A}}^N{Y}_K^N+{\overline{B}}^1{Y}_{K+M}^1+\cdots +{\overline{B}}^N{Y}_{K+M}^N,$$

(24)

where $\Lambda$ is the submatrix of the basic block tridiagonal matrix of system (12) that corresponds to the interval $(K,K+M)$.

If matrix $\Lambda$ is inversible, then

$$\overline{Y}={\Lambda }^{-1}\overline{F}+{\Lambda }^{-1}{\overline{A}}^1{Y}_K^1+\cdots +{\Lambda }^{-1}{\overline{A}}^N{Y}_K^N+{\Lambda }^{-1}{\overline{B}}^1{Y}_{K+M}^1+\cdots +{\Lambda }^{-1}{\overline{B}}^N{Y}_{K+M}^N.$$

From this, it follows that

$$\overline{Y}={\overline{U}}^1{Y}_K^1+\cdots +{\overline{U}}^N{Y}_K^N+{\overline{V}}^1{Y}_{K+M}^1+\cdots +{\overline{V}}^N{Y}_{K+M}^N+\overline{W},$$

(25)

where

$${\overline{U}}^1\mathrm{is}a\mathrm{solution}\mathrm{of}\Lambda {\overline{U}}^1={\overline{A}}^1,\cdots ,$$

$${\overline{U}}^N\mathrm{is}a\mathrm{solution}\mathrm{of}\Lambda {\overline{U}}^N={\overline{A}}^N,$$

$${\overline{V}}^1\mathrm{is}a\mathrm{solution}\mathrm{of}\Lambda {\overline{V}}^1={\overline{B}}^1,\cdots ,$$

(26)

$${\overline{V}}^N\mathrm{is}a\mathrm{solution}\mathrm{of}\Lambda {\overline{V}}^N={\overline{B}}^N,$$

$$\overline{W}\mathrm{is}a\mathrm{solution}\mathrm{of}\Lambda \overline{W}=\overline{F}.$$

Let us rewrite ${\overline{U}}^1,\cdots ,{\overline{U}}^N,{\overline{V}}^1,\cdots ,{\overline{V}}^N,\overline{W}$ in a more detailed manner.

$${\overline{U}}^1=\left(\begin{array}{c}{\overline{U}}_{K+1}^1\\ ...\\ {\overline{U}}_{K+M-1}^1\end{array}\right),\cdots ,{\overline{U}}^N=\left(\begin{array}{c}{\overline{U}}_{K+1}^N\\ ...\\ {\overline{U}}_{K+M-1}^N\end{array}\right),{\overline{V}}^1=\left(\begin{array}{c}{\overline{V}}_{K+1}^1\\ ...\\ {\overline{V}}_{K+M-1}^1\end{array}\right),\cdots ,$$

$${\overline{V}}^N=\left(\begin{array}{c}{\overline{V}}_{K+1}^N\\ ...\\ {\overline{V}}_{K+M-1}^N\end{array}\right),\overline{W}=\left(\begin{array}{c}{\overline{W}}_{K+1}\\ ...\\ {\overline{W}}_{K+M-1}\end{array}\right).$$

(27)

Taking into account (23) and (27), problems (26) are equivalent to following:

$$\left\{\begin{array}{c}{C}_{K+1}{\overline{U}}_{K+1}^1-B_{K+1}{\overline{U}}_{K+2}^1={\overline{A}}_{K+1}^1,\hfill \\ ...\hfill \\ -A_{K+M-1}{\overline{U}}_{K+M-2}^1+C_{K+M-1}{\overline{U}}_{K+M-1}^1=0.\hfill \end{array}\right.$$

$$\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots .$$

$$\left\{\begin{array}{c}{C}_{K+1}{\overline{U}}_{K+1}^N-B_{K+1}{\overline{U}}_{K+2}^N={\overline{A}}_{K+1}^N;\hfill \\ ...\hfill \\ -A_{K+M-1}{\overline{U}}_{K+M-2}^N+C_{K+M-1}{\overline{U}}_{K+M-1}^N=0;\hfill \end{array}\right.$$

$$\left\{\begin{array}{c}{C}_{K+1}{\overline{V}}_{K+1}^1-B_{K+1}{\overline{V}}_{K+2}^1=0;\hfill \\ ...\hfill \\ -A_{K+M-1}{\overline{V}}_{K+M-2}^1+C_{K+M-1}{\overline{V}}_{K+M-1}^1={\overline{B}}_{K+M-1}^1;\hfill \end{array}\right.$$

(28)

$$\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots .$$

$$\left\{\begin{array}{c}{C}_{K+1}{\overline{V}}_{K+1}^N-B_{K+1}{\overline{V}}_{K+2}^N=0;\hfill \\ ...\hfill \\ -A_{K+M-1}{\overline{V}}_{K+M-2}^N+C_{K+M-1}{\overline{V}}_{K+M-1}^N={\overline{B}}_{K+M-1}^1;\hfill \end{array}\right.$$

$$\left\{\begin{array}{c}{C}_{K+1}{\overline{W}}_{K+1}-B_{K+1}{\overline{W}}_{K+2}={\overline{F}}_{K+1};\hfill \\ ...\hfill \\ -A_{K+M-1}{\overline{W}}_{K+M-2}+C_{K+M-1}{\overline{W}}_{K+M-1}={\overline{F}}_{K+M-1}.\hfill \end{array}\right.$$

Systems (28) for intervals $(K,K+M),K=0,M,\cdots ,N$ are equivalent to problems (16)–(18).

Thus, the original solution (25) on interval $(K,K+M)$ has the form

$${\overline{Y}}_i=({\overline{U}}_i^1{\overline{U}}_i^2\cdots {\overline{U}}_i^n){\overline{Y}}_K+({\overline{V}}_i^1{\overline{V}}_i^2\cdots {\overline{V}}_i^n){\overline{Y}}_{K+M}+{\overline{W}}_i,i=K+1,\cdots ,K+M-1.$$

   □

Remark 3.

Matrices $U_i,V_i$ and vector ${\overline{W}}_i$ can be obtained by the block elimination methods (13) and (14). For arbitrary L block elimination formulae for intervals $(K,K+M)$, $K=0,M,\cdots ,N$ have the following form.

The forward phase for Equations (16)–(18) is as follows:

$$\begin{array}{cc}\hfill & {\alpha }_{K+1}=C_{K+1}^{-1}{B}_{K+1},{\beta }_{K+1}=C_{K+1}^{-1}{A}_{K+1},{\overline{\gamma }}_{K+1}=C_{K+1}^{-1}{\overline{F}}_{K+1},\hfill \\ \hfill & {\alpha }_i={[C_i-A_i{\alpha }_{i-1}]}^{-1}{B}_i,{\beta }_i={[C_i-A_i{\alpha }_{i-1}]}^{-1}{A}_i{\beta }_{i-1},\hfill \\ \hfill & {\overline{\gamma }}_i={[C_i-A_i{\alpha }_{i-1}]}^{-1}({\overline{F}}_i+A_i{\overline{\gamma }}_{i-1}),i=K+2,\cdots ,K+M-1.\hfill \end{array}$$

(29)

The backward phase for Equations (16)–(18) is as follows:

$$\begin{array}{cc}\hfill & U_{K+M-1}={\beta }_{K+M-1},V_{K+M-1}={\alpha }_{K+M-1},{\overline{W}}_{K+M-1}={\overline{\gamma }}_{K+M-1},\hfill \\ \hfill & U_i={\alpha }_i{U}_{i+1}+{\beta }_i,V_i={\alpha }_i{V}_{i+1},{\overline{W}}_i={\alpha }_i{\overline{W}}_{i+1}+{\overline{\gamma }}_i,i=K+M-2,\cdots ,K+1.\hfill \end{array}$$

(30)

When we substitute expression (19) for indices $K=0,M,\cdots ,N$ into the basic system (12), we will obtain a reduced SLAE for the parametric unknown values (vectors ${\overline{Y}}_K$). This system has a similar structure to (12) but has a smaller size.

$$\left\{\begin{array}{c}[C_0-B_0{U}_1]{\overline{Y}}_0-\left[B_0{V}_1\right]{\overline{Y}}_M={\overline{F}}_0+B_0{\overline{W}}_1,K=0,\hfill \\ -\left[A_K{U}_{K-1}\right]{\overline{Y}}_{K-M}+[C_K-A_K{V}_{K-1}-B_K{U}_{K+1}]{\overline{Y}}_K-\left[B_K{V}_{K+1}\right]{\overline{Y}}_{K+M}=\hfill \\ ={\overline{F}}_K+A_K{\overline{W}}_{K-1}+B_K{\overline{W}}_{K+1},K=M,\cdots ,N-M,\hfill \\ -\left[A_N{U}_{N-1}\right]{\overline{Y}}_{N-M}+[C_N-A_N{V}_{N-1}]{\overline{Y}}_N={\overline{F}}_N+A_N{\overline{W}}_{N-1},K=N,\hfill \end{array}\right.$$

(31)

where $U_K$ and $V_K$ are matrices of dimension $n\times n$.

Problem (31) is solved by the block-elimination method (13) and (14) in the single-threaded mode or on a single node of a cluster. Auxilliary problems (16)–(18) are solved independently for each of the L intervals. Thus, this workload can be distributed between L threads or processes. After obtaining ${\overline{Y}}_K$, the rest of the unknown values are found by formula (19). This work also can be performed independently for each of the L intervals.

The parallel matrix sweep algorithm for solving system (12) is presented in Listing 1.

Listing 1. Parallel matrix sweep algorithm for solving SLAE.

1. Find values $U_i,V_i,{\overline{W}}_i$ from problems (16)–(18) for inner points $i\in (K,K+M)$ of each subinterval $K=0,M,2M,...,N$. Methods (29) and (30) are used to solve the subproblems independently on each of the subintervals. 2. Calculate the coefficients for reduced system (31). The coefficients may be calculated independently for each K, but to solve the resulting reduced system, we need to transfer these coefficients to a single process. This requires synchronization or gather-type communication. 3. Find values $\overline{Y_K}$ from the reduced system (31). Compared to the basic system (12), its dimension is much smaller. It is solved by the block elimination algorithms (13) and (14) in serial mode. The computed values ${\overline{Y}}_K$ must be transmitted to processors. This step requires synchronization or communication. 4. Use formula (19) to calculate the sought values ${\overline{Y}}_i,i\in (K,K+M)$. These computations may be performed independently for each subinterval K.

Thus, steps 1, 2, and 4 of this algorithm may be parallelized, while step 3 must be performed in serial mode.

To establish the stability (see Remark 1) of the parallel matrix sweep algorithm, let us prove the following theorems.

Theorem 2.

If original system (12) satisfies the condition

$$\parallel C_i\parallel \ge \parallel A_i\parallel +\parallel B_i\parallel +\delta ,\delta >0,$$

then reduced system (31) also satisfies this condition in the form

$$\parallel C_K-A_K{V}_{K-1}-B_K{U}_{K+1}\parallel \ge \parallel A_K{U}_{K-1}\parallel +\parallel B_K{V}_{K+1}\parallel +\delta .$$

Proof. 

$$\begin{array}{cc}\hfill & \parallel C_K-A_K{V}_{K-1}-B_K{U}_{K+1}\parallel \ge \parallel C_K\parallel -\parallel A_K{V}_{K-1}\parallel -\parallel B_K{U}_{K+1}\parallel \ge \hfill \\ \hfill \ge & \parallel C_K\parallel -\parallel A_K(I-U_{K-1})\parallel -\parallel B_K(I-V_{K+1})\parallel \ge \hfill \\ \hfill \ge & \parallel C_K\parallel -\parallel A_K\parallel +\parallel A_K{U}_{K-1}\parallel -\parallel B_K\parallel +\parallel B_K{V}_{K+1}\parallel \ge \hfill \\ \hfill \ge & \parallel A_K{U}_{K-1}\parallel +\parallel B_K{V}_{K+1}\parallel +\parallel C_K\parallel -\parallel A_K\parallel -\parallel B_K\parallel \ge \hfill \\ \hfill \ge & \parallel A_K{U}_{K-1}\parallel +\parallel B_K{V}_{K+1}\parallel +\delta ,\hfill \end{array}$$

since $\parallel U_K\parallel +\parallel V_K\parallel \le 1$.    □

Theorem 3.

If basic system (12) satisfies the stability conditions of the matrix sweep method (Lemma 1), then these conditions are sufficient for the stability of the matrix sweep method for reduced system (31) for ${\overline{Y}}_K$.

Proof. 

Let us construct a proof by the mathematical induction method.

We will utilize the following statement [31]. If square matrix S satisfies $\parallel S\parallel \le q\le 1$, then matrices ${(E-S)}^{-1}$ and ${\parallel (E-S)}^{-1}\parallel \le 1/(1-q)$ must exist.

Let ${\tilde{\alpha }}_1,\cdots ,{\tilde{\alpha }}_K,\cdots ,{\tilde{\alpha }}_N$ be the elimination coefficients for the matrix sweep method for system (31).

1. Let us demonstrate that $\parallel {\tilde{\alpha }}_1\parallel \le 1.$

$$\parallel C_0^{-1}{B}_0{U}_1\parallel \le \parallel C_0^{-1}{B}_0\parallel ·\parallel U_1\parallel <1,$$

therefore, there are ${(E-C_0^{-1}{B}_0{U}_1)}^{-1}$ and

$$\parallel {\tilde{\alpha }}_1\parallel =\parallel {(C_0-B_0{U}_1)}^{-1}{B}_0{V}_1\parallel \le \parallel C_0^{-1}{(E-C_0^{-1}{B}_0{U}_1)}^{-1}{B}_0{V}_1\parallel \le$$

$$\le \parallel {(E-C_0^{-1}{B}_0{U}_1)}^{-1}\parallel ·\parallel C_0^{-1}{B}_0\parallel ·\parallel V_1\parallel \le \frac{1}{1-\parallel C_0^{-1}{B}_0\parallel ·\parallel U_1\parallel }·\parallel C_0^{-1}{B}_0\parallel ·\parallel V_1\parallel \le$$

$$\le \frac{\parallel V_1\parallel }{1-\parallel U_1\parallel }\le \frac{\parallel V_1\parallel }{\parallel V_1\parallel }=1,\mathrm{since}\parallel U_1\parallel +\parallel V_1\parallel \le 1.$$

2. Assume $\parallel {\tilde{\alpha }}_K\parallel \le 1.$ Let us demonstrate that $\parallel {\tilde{\alpha }}_{K+1}\parallel \le 1.$

$$\parallel {\tilde{\alpha }}_{K+1}\parallel =\parallel {(C_K-A_K{V}_{K-1}-B_K{U}_{K+1}-A_K{U}_{K-1}{\tilde{\alpha }}_K)}^{-1}·B_K{V}_{K+1}\parallel .$$

Consider,

$$\parallel {(C_K^{-1}{A}_K{V}_{K-1}+C_K^{-1}{B}_K{U}_{K+1}+C_K^{-1}{A}_K{U}_{K-1}{\tilde{\alpha }}_K)}^{-1}\parallel \le$$

$$\le \parallel C_K^{-1}{A}_K\parallel ·\parallel V_{K-1}\parallel +\parallel C_K^{-1}{B}_K\parallel ·\parallel U_{K+1}\parallel +\parallel C_K^{-1}{A}_K\parallel ·\parallel U_{K-1}\parallel \le$$

$$\le \parallel C_K^{-1}{A}_K\parallel +\parallel C_K^{-1}{B}_K\parallel ·\parallel U_{K+1}\parallel \le 1-\parallel C_K^{-1}{B}_K\parallel +\parallel C_K^{-1}{B}_K\parallel ·\parallel U_{K+1}\parallel \le$$

$$\le 1-\parallel C_K^{-1}{B}_K\parallel ·(1-\parallel U_{K+1}\parallel )\le 1-\parallel C_K^{-1}{B}_K\parallel ·\parallel V_{K+1}\parallel <1,$$

since

$$\parallel U_{K+1}\parallel +\parallel V_{K+1}\parallel \le 1.$$

Therefore, there are

$${(E-C_K^{-1}{A}_K{V}_{K-1}+C_K^{-1}{B}_K{U}_{K+1}+C_K^{-1}{A}_K{U}_{K-1}{\tilde{\alpha }}_K)}^{-1}\mathrm{and}$$

$$\parallel {\tilde{\alpha }}_{K+1}\parallel \le \parallel {(E-C_K^{-1}{A}_K{V}_{K-1}+C_K^{-1}{B}_K{U}_{K+1}+C_K^{-1}{A}_K{U}_{K-1}{\tilde{\alpha }}_K)}^{-1}\parallel \times$$

$$\times \parallel C_K^{-1}{B}_K{V}_{K+1}\parallel \le \frac{1}{\parallel C_K^{-1}{B}_K\parallel ·\parallel V_{K+1}\parallel }·\parallel C_K^{-1}{B}_K\parallel ·\parallel V_{K+1}\parallel =1.$$

  3. Let us demonstrate that $\parallel {\tilde{\alpha }}_N\parallel \le 1.$

$$\parallel C_N^{-1}{A}_N{V}_{N-1}\parallel \le \parallel C_N^{-1}{A}_N\parallel ·\parallel V_{N-1}\parallel <1,$$

therefore, there are

$${(E-C_N^{-1}{A}_N{V}_{N-1})}^{-1}.$$

Consider

$$\parallel {\tilde{\alpha }}_N\parallel =\parallel {(C-A_N{V}_{N-1})}^{-1}{A}_N{U}_{N-1}\parallel \le$$

$$\le \parallel {(E-C_N^{-1}{A}_N{V}_{N-1})}^{-1}\parallel ·\parallel C_N^{-1}{A}_N\parallel ·\parallel U_{N-1}\parallel \le$$

$$\le \frac{1}{\parallel C_N^{-1}{A}_N\parallel ·\parallel V_{N-1}\parallel }·\parallel C_N^{-1}{A}_N\parallel ·\parallel U_{N-1}\parallel \le \frac{\parallel U_{N-1}\parallel }{1-\parallel V_{N-1}\parallel }\le 1,$$

since $\parallel U_{N-1}\parallel +\parallel V_{N-1}\parallel \le 1.$    □

4. Numerical Method for Solving the Inverse Problems

To solve the inverse problems (1)–(3) and (5), we use the iterative conjugate gradient method [23,33]. Consider that additional information $\phi$ may contain a random perturbation ${\phi }^{\delta }=\phi ·(1+rand(-\delta ,\delta ))$. To overcome this, we will regularize the inverse problems using the Lavrentyev scheme [34].

The resulting algorithm for solving the inverse problems is presented in Listing 2, where $\epsilon >0$ is the regularization parameter.

Listing 2. Regularized conjugate gradient algorithm for solving the inverse problems.

Initialization: 1. Set $s=0$ as the iterative step. 2. Set the initial approximation ${\psi }_0\left(\mathbf{x}\right)$, for example, ${\psi }_0\left(\mathbf{x}\right)=0$. 3. Solve the initial boundary problems (1)–(3) by substituting the right-hand part $\psi$ with ${\psi }_0$; obtain $U_T^0=U(\mathbf{x},T){|}_{{\psi }_0}$. 4. Calculate the initial residual $r_0\left(\mathbf{x}\right)=\phi \left(\mathbf{x}\right)-(U_T^0+\epsilon {\psi }_0)$ and initial estimation $p_0\left(\mathbf{x}\right)=r_0\left(\mathbf{x}\right)$, where $\epsilon >0$ is the regularization parameter. Iterations: 5. Set $s=s+1$. 6. Solve the initial boundary problems by substituting the right-hand part with $p_s\left(\mathbf{x}\right)$; obtain $U_T^s=U(\mathbf{x},T){|}_{p_s}$. 7. Calculate the coefficient ${\alpha }_s=(r_s,r_s)/(p_s,(U_T^s+\epsilon p_s))$. 8. Calculate the estimation and residual for next step ${\psi }_{s+1}={\psi }_s+{\alpha }_s{p}_s,r_{s+1}=r_s-{\alpha }_s(U_T^s+\epsilon p_s)$. 9. Calculate the coefficient ${\beta }_s=(r_{s+1},r_{s+1})/(r_s,r_s)$. 10. Calculate $p_{s+1}=r_{s+1}+{\beta }_s{p}_s.$ 11. Check the stopping rule $∥r_s∥/∥\phi ∥<\mu$, $0<\mu <1$. If not met, go to step 5.

5. Parallel Implementation of Algorithms for Solving the Inverse Problem

The numerical solution to the problems related to the fractional differential equation is an expensive task that requires a lot of computing time.

The most time-consuming subroutine of the regularized conjugate gradient algorithm (see Listing 2) is solving the auxiliary initial boundary problem at each iteration. In turn, this procedure consists in forming and solving SLAEs (11) at each subsequent time step.

5.1. Efficient Computation of the Right-Hand Parts

Forming SLAE requires the calculation of the right-hand parts using formula (9). In our earlier work [27], when solving one-dimensional problems, the fraction of time spent to calculate the right-hand part was up to 70% of the total time. To optimize this procedure, we implemented the logarithmic memory approach. It consists of using the non-uniform time grid when computing the approximation of the fractional derivative. The fine time step is used for the latest history part. For the more distant history, successively larger time steps are used. This approach allowed us to reduce the computing time for the one-dimensional case by up to 1.5 times.

This approach may also be utilized for the ctwo-dimensional problem. For solving system (11), a modified variant of formula (9) takes the form

$$\begin{array}{cc}\hfill & f_{i_1,i_2,1}={\sigma }_{\alpha ,\tau }{U}_{i_1,i_2,0}+{\psi }_{i_1,i_2}{\eta }_{i_1,i_2,0},\hfill \\ \hfill & f_{i_1,i_2,j}={\sigma }_{\alpha ,\tau }{U}_{i_1,i_2,j-1}-\sum _{(\mathsf{\ell },k)}{\sigma }_{\alpha ,\tau }^{(\mathsf{\ell },k)}{w}_{\mathsf{\ell },k}^{\left(\alpha \right)}(U_{i_1,i_2,j-\mathsf{\ell }+1}-U_{i_1,i_2,j-k})+{\psi }_{i_1,i_2}{\eta }_{i_1,i_2,j},j>1,\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill & (\mathsf{\ell },k)\in \left\{(2,2+{\theta }^0),(2+{\theta }^0,2+{\theta }^1),(2+{\theta }^1,2+{\theta }^2),\dots ,(2+{\theta }^{\lfloor {log}_{\theta }j\rfloor },j)\right\},\hfill \\ \hfill & {\sigma }_{\alpha ,\tau }^{(\mathsf{\ell },k)}=\frac{1}{\Gamma (1-\alpha )(1-\alpha ){\theta }^{k-\mathsf{\ell }}{\tau }^{\alpha }},w_{\mathsf{\ell },k}^{\left(\alpha \right)}={\left(k\right)}^{1-\alpha }-{(\mathsf{\ell }-1)}^{1-\alpha },0<\alpha <1,\hfill \end{array}$$

where $\theta \in \mathbb{N}$ is the stretching coefficient and $\lfloor {log}_{\theta }n\rfloor$ is the floor function (integer part). Note that approach has complexity $O(\tilde{N}·log\tilde{N})$ in contrast of $O\left({\tilde{N}}^2\right)$ of the uniform time grid.

In the next section, we will explore the usefulness of this approach for the case of a two-dimensional problem.

5.2. Parallel Implementation of the SLAE Solver

To speed up SLAE solving, we implement the parallel matrix sweep Algorithm (see Listing 1). For comparison, we also implemented the serial block-elimination methods (13) and (14).

The parallel algorithm for solving the inverse problem of finding the source term of the time-fractional diffusion equation was implemented for the multicore processor using OpenMP technology [35] and the Intel MKL library [36]. The parallelization is performed as follows.

• The workload of calculating the right-hand parts (9) and (32) is distributed to OpenMP threads utilizing to the same subinterval decomposition that is used for parallel matrix sweep algorithm.
• The matrix operations and the inversion of the matrix blocks are performed with MKL routines (gemv, gemm, getrf, and getri).
• In the parallel matrix sweep algorithm (Listing 1), steps 1, 2, and 4 are performed by the individual OpenMP threads on their corresponding subintervals. To perform step 3, thread synchronization is required. This is performed by the ‘#pragma omp barrier’ directive. Note that this synchronization requires additional time.

6. Numerical Experiments

In this section, we present the numerical experiments of solving the direct and inverse problems for the two-dimensional fractional diffusion equation. The experiments were performed using the developed code on the Intel i9-12900k CPU, which has 8 P-cores. The goal is to study the validity of the proposed numerical methods, as well as the efficiency of the parallel code.

6.1. Problem 1

Consider the two-dimensional equation

$$\frac{{\partial }^{\alpha }U(x_1,x_2,t)}{\partial t^{\alpha }}=\frac{{\partial }^{2}U(x_1,x_2,t)}{\partial x_1^2}+\frac{{\partial }^{2}U(x_1,x_2,t)}{\partial x_2^2}+\left(\frac{2}{\Gamma (3-\alpha )}{t}^{2-\alpha }+2t^2\right)sin\left(x_1\right)sin\left(x_2\right)$$

(33)

with initial and boundary conditions

$$U(x_1,x_2,0)=0,$$

$$U(x_1,0,0)=0,U(x_1,\pi ,0)=0,$$

$$U(0,x_2,0)=0,U(\pi ,x_2,0)=0,$$

and area

$$0\le x_1,x_2\le {\gamma }_1={\gamma }_2=\pi ,0\le t\le T=1,$$

for order

$$0<\alpha <1.$$

Paper [37] presents the exact solution for this equation

$$U(x_1,x_2,t)=t^{2}sin\left(x_1\right)sin\left(x_2\right).$$

6.1.1. Experiment 1

Experiment 1 consists of solving the forward (initial boundary) problem for Equation (33) on the various grids $n=N=\{128;256;512\}$, $\tilde{N}=\{64;128;256\}$ and various parameters $\alpha =\{0.5;0.8;0.95\}$. It was solved using the difference scheme described in Section 2.2. For solving the SLAE, two methods were applied, namely, the classical serial block-elimination method (13) and (14) and parallel matrix sweep method (see Listing 1).

Figure 1 shows the exact solution $U(x_1,x_2,1)=t^{2}sin\left(x_1\right)sin\left(x_2\right)$ and approximate solution $\tilde{U}(x_1,x_2,1)$ for Problem 1 obtained by the parallel matrix sweep algorithm for grid size $n=N=512$, $\tilde{N}=256$, and the order of fractional derivative $\alpha =0.5$. The approximate solutions obtained by the matrix sweep and parallel matrix sweep methods coincide with each other up to the machine precision (${10}^{-15}$, as we used the double precision format).

Table 1. Results of Experiment 1 for Problem 1: relative error of solving the direct problem.

	Grid Size	$\tilde{\mathit{N}}=64$	$\tilde{\mathit{N}}=128$	$\tilde{\mathit{N}}=256$
$\alpha =0.5$	$n=N=128$	$3.64\times {10}^{-4}$	$1.48\times {10}^{-4}$	$7.07\times {10}^{-5}$
	$n=N=256$	$3.45\times {10}^{-4}$	$1.26\times {10}^{-4}$	$4.93\times {10}^{-5}$
	$n=N=512$	$3.40\times {10}^{-4}$	$1.21\times {10}^{-4}$	$1.98\times {10}^{-5}$
$\alpha =0.8$	$n=N=128$	$2.25\times {10}^{-3}$	$9.44\times {10}^{-4}$	$4.23\times {10}^{-4}$
	$n=N=256$	$2.13\times {10}^{-3}$	$9.25\times {10}^{-4}$	$4.04\times {10}^{-4}$
	$n=N=512$	$2.12\times {10}^{-3}$	$9.20\times {10}^{-4}$	$4.02\times {10}^{-4}$
$\alpha =0.95$	$n=N=128$	$5.14\times {10}^{-3}$	$2.48\times {10}^{-3}$	$1.20\times {10}^{-3}$
	$n=N=256$	$5.12\times {10}^{-3}$	$2.46\times {10}^{-3}$	$1.18\times {10}^{-3}$
	$n=N=512$	$5.11\times {10}^{-3}$	$2.45\times {10}^{-3}$	$1.18\times {10}^{-3}$

contains the relative error of the solutions $∥U-\tilde{U}∥/∥U∥$ for various grid sizes and parameters $\alpha$. The experiments show that taking a finer grid either for space or time reduces the relative error of the resulting solution. For the time grid, the rate of convergence is close to linear (increasing the number of grid points twofold reduces the error approximately by two times). As parameter $\alpha$ approaches 1, the error of the solution increases. To achieve higher accuracy, we need to use a finer grid.

Table 2. Results of Experiment 1 for Problem 1: computing time of solving the direct problem.

Method	Number L of OpenMP Threads	${\mathit{T}}_{\mathit{L}}$ (Minutes)	${\mathit{T}}_{\mathbf{SLAE}}$	${\mathit{T}}_{\mathbf{Right}}$
Matrix Sweep (13) and (14)	Serial	28.6	28.5	0.13
Parallel Matrix Sweep (Listing 1)	2	30.7	30.63	0.06
Parallel Matrix Sweep	4	16.3	16.26	0.03
Parallel Matrix Sweep	8	10.1	10.08	0.016
Parallel Matrix Sweep	16	10.5	10.48	0.017

presents the total computing time $T_L$ of solving the direct problem. For the parallel matrix sweep method, the computing time for various numbers L of OpenMP threads are presented. Total time $T_L$ consists of time $T_{SLAE}$ spent on solving the SLAEs plus time $T_{Right}$ spent on computing the right-hand part for these SLAEs using formula (9). The table shows that for the case of two-dimensional problem solving SLAE with a more complex matrix structure (block-tridiagonal), the time spent on solving the SLAE is up to 600 times larger than the time spent on computing the right-hand part. Thus, utilizing the optimized approach for computing the fractional derivative is less relevant as it would bring a miniscule speedup.

6.1.2. Experiment 2

Experiment 2 consists of solving the inverse problem for Equation (33). We assume that $\eta \left(t\right)=\left(\frac{2}{\Gamma (3-\alpha )}{t}^{2-\alpha }+2t^2\right)$ and $\phi (x_1,x_2)=sin\left(x_1\right)sin\left(x_2\right)$. Thus, we need to solve the inverse problem of finding unknown $\left[U(x_1,x_2,T),\psi (x_1,x_2)\right]$. For this experiment, we introduce the varying level $\delta$ of random perturbation to the a priori data ${\phi }^{\delta }=\phi ·(1+rand(-\delta ,\delta ))$.

Remark 4.

Note that this level corresponds to an error in the infinity norm, i.e., $\delta \approx {∥\phi -{\phi }^{\delta }∥}_{\infty }/{∥\phi ∥}_{\infty }$. In the rest of the paper, the norm is implied to be $L_2$-norm. In tables, we provide the corresponding ${\delta }_2=∥\phi -{\phi }^{\delta }∥/∥\phi ∥$.

The inverse problem was solved by the regularized conjugate gradient method Algorithm (see Listing 2). For solving the SLAEs, the parallel matrix sweep method was used. The grid size was $n=N=256$, $\tilde{N}=256$.

Figure 2 shows the exact solution $\psi (x_1,x_2)=sin\left(x_1\right)sin\left(x_2\right)$ and the approximate solution $\tilde{\psi }(x_1,x_2)$ for Problem 1 for the noise level $\delta =0.02$.

Table 3. Results of Experiment 2 for Problem 1: solving the inverse problem.

Noise Level $\mathit{\delta }$	Noise Level ${\mathit{\delta }}_2$	Regularization Parameter $\mathit{\epsilon }$	Stopping Rule $\mathsf{\mu }$	Number of Iterations S	Error of Solution $∥\mathit{\psi }-\tilde{\mathit{\psi }}∥/∥\mathit{\psi }∥$
0	0	0	0.001	3	5$\times {10}^{-5}$
0.01	0.006	0.1	0.006	2	0.1
0.02	0.01	0.1	0.01	2	0.14
0.05	0.03	0.2	0.03	1	0.16

presents the results of Experiment 2 for varying levels of noise $\delta ,{\delta }_2$. It contains the values of regularization parameter $\epsilon$, the threshold $\mu$ for the stopping criterion (we used $\mu ={\delta }_2$ for experiments), number of iterations S, and the relative error of the resulting solution.

6.2. Problem 2

Consider the two-dimensional equation [22]

$$\frac{{\partial }^{\alpha }U(x_1,x_2,t)}{\partial t^{\alpha }}=\frac{{\partial }^{2}U(x_1,x_2,t)}{\partial x_1^2}+\frac{{\partial }^{2}U(x_1,x_2,t)}{\partial x_2^2}+sin\left(x_1\right)sin\left(x_2\right)+sin\left(2x_1\right)sin\left(2x_2\right)$$

(34)

with initial and boundary conditions

$$U(x_1,x_2,0)=0,$$

$$U(x_1,0,0)=0,U(x_1,\pi ,0)=0,$$

$$U(0,x_2,0)=0,U(\pi ,x_2,0)=0,$$

and area

$$0\le x_1,x_2\le \pi ,0\le t\le 1,$$

for order

$$0<\alpha <1.$$

The numerical experiment consists of solving the inverse problem for Equation (34). We assume that $\eta (x_1,x_2,t)=1$. Additional data $\phi (x_1,x_2)=U(x_1,x_2,1)$ are obtained by solving the direct problem substituting exact $\psi (x_1,x_2)=sin\left(x_1\right)sin\left(x_2\right)+sin\left(2x_1\right)sin\left(2x_2\right)$. They are shown in Figure 3.

Table 4. Results of experiments for Problem 2: solving the inverse problem.

Noise Level $\mathit{\delta }$	Noise Level ${\mathit{\delta }}_2$	Regularization Parameter $\mathit{\epsilon }$	Stopping Rule $\mathsf{\mu }$	Number of Iterations S	Error of Solution $∥\mathit{\psi }-\tilde{\mathit{\psi }}∥/∥\mathit{\psi }∥$
0	0	0	0.001	3	3.2$\times {10}^{-4}$
0.01	0.006	0.02	0.006	3	0.13
0.02	0.01	0.02	0.01	3	0.19
0.05	0.03	0.05	0.03	1	0.23

presents the results of experiments for Problem 2 for varying levels of noise $\delta ,{\delta }_2$. It contains the values of the regularization parameter $\epsilon$, the threshold $\mu$ for the stopping criterion (we used $\mu ={\delta }_2$ for experiments), the number of iterations S, and the relative error of the resulting solution. The grid size was $n=N=256$, $\tilde{N}=256$, and order $\alpha =0.5$.

Figure 4 shows the approximate solutions $\tilde{\psi }(x_1,x_2)$ for Problem 2 for various noise levels.

7. Discussion

According to the experiments, the relative error of the direct problem solution decreases with finer grid size. This indicates the experimental confirmation of convergence of the finite difference scheme.

In the case of two-dimensional problem solving, SLAE takes a significantly larger time (up to 600) than computing the right-hands part for SLAE. This makes the parallel implementation of the SLAE solver more important than the optimization of the procedures for computing the fractional derivative.

Experiments show that the parallel matrix sweep method for solving the SLAE has good parallel efficiency. The minimal computing time is achieved by using eight OpenMP threads on an eight-core processor.

The parallel code performance is mainly limited by memory bandwidth. Adding more than eight threads does not reduce the computing time. The largest speed up is only three-fold for $512\times 512$ spatial grid. Figure 5 presents the roofline analysis performed by the Intel Advisor tool. Most subroutines of the parallel code lie primarily below and near the slanted line that represents DRAM bandwidth. This indicates that the code is memory-bound.

Several approaches can be used to overcome this limitation. Computing hardware with large memory bandwidth, such as graphics processors (GPU), can be used. Central processors with DDR4 or DDR5 RAM can achieve up to 100 GB/s, while modern GPUs have a bandwidth of 1000 GB/s or higher. Another option is to use massive distributed memory systems. Since each node works independently, the memory speed of individual node is effectively summed. Moreover, it enables larger problems to be solved with data that cannot be accommodated in the memory of a single computing node.

The experiments in Table 3 and Table 4 show that for the model problems, the regularized conjugate gradient method allows us to solve the inverse problem even with noised data. The results are comparable with other works in terms of accuracy (for example, see [23], Table 3).

We also note that while this work is devoted to the case of the two-dimensional equation, the results may be extended to three-dimensional elliptic equations. The structure of matrix A in Equation (11) will remain block tridiagonal, but the inner structure of the blocks will be more complex.

8. Conclusions

In this work, we construct the parallel algorithm for solving the inverse problem of finding the space-dependent component of a source term in a two-dimensional fractional diffusion equation. The considered inverse problem is solved by the iterative conjugate gradient method. At each iteration, it is necessary to solve an auxiliary direct initial-boundary value problem. Applying the finite difference scheme, we reduce the initial-boundary value problem to solving an SLAE with block tridiagonal matrices at each subsequent time level. For the efficient solution of such SLAEs, we construct and implement the direct parallel matrix sweep method. Stability and correctness for the parallel matrix sweep method are established. In the two-dimensional case, computing the fractional derivative (the right-hand part of the SLAE) takes little time in comparison with solving the SLAE.

The algorithm is implemented for the multicore processors using the OpenMP technology. In the numerical experiments, we investigated the validity of numerical methods and the efficiency and speedup of the parallel algorithm. The utilization of the parallel sweep algorithm reduces the computing time by up to three times on a eight-core processor. Using Lavrentyev regularization method allows us to solve the inverse problem with a disturbed data.

In future, the authors plan to implement a similar approach to solving the retrospective inverse problem (identifying the initial value) for a fractional differential equation. The developed algorithms may be utilized for real applications. The parallel algorithms will be implemented on graphics processors.