# Solving the Caputo Fractional Reaction-Diffusion  Equation on GPU

Authors: Jie Liu; Chunye Gong; Weimin Bao; Guojian Tang; Yuewen Jiang
Journal: Discrete Dynamics in Nature and Society
 (2014)

---

## Abstract

We present a parallel GPU solution of the Caputo fractional reaction-diffusion equation in one spatial dimension with explicit finite difference approximation. The parallel solution, which is implemented with CUDA programming model, consists of three procedures: preprocessing, parallel solver, and postprocessing. The parallel solver involves the parallel tridiagonal matrix vector multiplication, vector-vector addition, and constant vector multiplication. The most time consuming loop of vector-vector addition and constant vector multiplication is optimized and impressive performance improvement is got. The experimental results show that the GPU solution compares well with the exact solution. The optimized GPU solution on NVIDIA Quadro FX 5800 is 2.26
times faster than the optimized parallel CPU solution on multicore Intel Xeon E5540 CPU.


---

## Body

## 1. Introduction

The idea of fractional derivatives can be dated back to the 17th century. A fractional differential equation is a kind of equation which uses fractional derivatives. Fractional equations can be used to describe some physical phenomenons more accurately than the classical integer order differential equation [1]. The reaction-diffusion equations play an important role in dynamical systems of mathematics, physics, chemistry, bioinformatics, finance, and other research areas.Some analytical methods were proposed for fractional differential equations [2, 3]. The stability of Cauchy fractional differential equations was studied [4, 5] and more attention should be paid to the interesting Ulam’s type stability [6]. There have been a wide variety of numerical approximation methods proposed for fractional equations [7, 8], for example, finite difference method [9], finite element method, and spectral techniques [10]. Interest in fractional reaction-diffusion equations has increased [11]. In 2000, Henry and Wearne [12] derived a fractional reaction-diffusion equation from a continuous-time random walk model with temporal memory and sources. The fractional reaction-diffusion system with activator and inhibitor variables was studied by Gafiychuk et al. [13]. Haubold et al. [14] developed a solution in terms of the H-function for a unified reaction-diffusion equation. The generalized differential transform method [15] was presented for fractional reaction-diffusion equations. Saxena et al [16] gave investigation of a closed form solution of a generalized fractional reaction-diffusion equation.Parallel computing is used to solve computation intensive applications simultaneously [17–19]. In recent years, the computing accelerators such as graphics processing unit (GPU) provided a new parallel method of accelerating computation intensive simulations [20–22]. The use of general purpose GPU is possible by the advance of programming models and hardware resources. The GPU programming models such as NVIDIA’s compute unified device architecture (CUDA) [23] become more mature than before and simplify the development of nongraphic applications. GPU presents an energy efficient architecture for computation intensive domains like particle transport [24, 25] and molecular dynamics [26].It is time consuming to numerically solve fractional differential equations for high spatial dimension or big time integration. Short memory principle [27, 28] and parallel computing [29–32] can be used to overcome this difficulty. The parallel algorithms of one- and two- dimensional time fractional equations are studied and good parallel scalability is got [31, 32]. Optimization of the sum of constant vector multiplication is presented and 2-time speedup can be got [31]. The parallel implicit iterative algorithm was studied for two-dimensional time fractional problem at the first time [32].Gong et al. [29] presented a parallel algorithm for Riesz space fractional equations. The parallel efficiency of the presented parallel algorithm of 64 processes is up to 79.39% compared with 8 processes on a distributed memory cluster system. Diethelm [30] implemented the fractional version of the second-order Adams-Bashforth-Moulton method for fractional ordinary equations on a parallel computer. Domain decomposition method is regarded as the basic mathematical background for many parallel applications. A domain decomposition algorithm for time fractional reaction-diffusion equation with implicit finite difference method was proposed [33]. The domain decomposition algorithm keeps the same parallelism but needs much fewer iterations, compared with Jacobi iteration in each time step, until nothing has been recorded on accelerating the numerical solution of Caputo fractional reaction-diffusion equation on GPU.This paper focuses on the Caputo fractional reaction-diffusion equation:(1)
0
D
t
α
u
(
x
,
t
)
+
μ
u
(
x
,
t
)
=
∂
2
u
(
x
,
t
)
∂
x
2
+
K
f
(
x
,
t
)

(
0
<
α
<
1
)
u
(
x
,
0
)
=
ϕ
(
x
)
,
x
∈
[
0
,
x
R
]
u
(
0
,
t
)
=
u
(
x
R
,
t
)
=
0
,
x
∈
[
0
,
T
]
on a finite domain 0
≤
x
≤
x
R and 0
≤
t
≤
T. The μ and K are constants. If α equals 1, (1) is the classical reaction-diffusion equation. The fractional derivative is in the Caputo form.


## 2. Background

### 2.1. Numerical Solution

The fractional derivative off
(
t
) in the Caputo sense is defined as [27]
(2)
0
C
D
t
α
f
(
t
)
=
1
Γ
(
1
-
α
)
∫
0
t
‍
f
′
(
ξ
)
(
t
-
ξ
)
α
d
ξ
(
0
<
α
<
1
)
.Iff
′
(
t
) is continuous bounded derivatives in [
0
,
T
] for every T
>
0, we can get
(3)
0
C
D
t
α
f
(
t
)
=
lim
⁡
ξ
→
0
,
n
ξ
=
t
ξ
α
∑
i
=
0
n
‍
(
-
1
)
i
(
α
i
)
=
f
(
0
)
t
-
α
Γ
(
1
-
α
)
+
1
Γ
(
1
-
α
)
∫
0
t
‍
f
′
(
ξ
)
(
t
-
ξ
)
α
d
ξ
.Defineτ
=
T
/
N, h
=
x
R
/
(
M
+
1
), t
n
=
n
τ, and x
i
=
0
+
i
h for 0
≤
n
≤
N
,
0
≤
i
≤
M
+
1. Define u
i
n, φ
i
n, and ϕ
i as the numerical approximation to u
(
x
i
,
t
n
), f
(
x
i
,
t
n
), and ϕ
(
x
i
). We can get [11]
(4)
0
C
D
t
α
u
(
x
,
t
)
|
x
i
t
n
=
1
τ
Γ
(
1
-
α
)
[
b
0
u
i
n
-
∑
k
=
1
n
-
1
‍
(
b
n
-
k
-
1
-
b
n
-
k
)
u
i
k
-
b
n
-
1
u
i
0
]
+
◯
(
τ
2
-
α
)
,
where 1
≤
i
≤
M, n
≥
1, and
(5)
b
l
=
τ
1
-
α
1
-
α
[
(
l
+
1
)
1
-
α
-
l
1
-
α
]
,
l
≥
0
.Using explicit center difference scheme for∂
2
u
(
x
,
t
)
/
∂
x
2, we can get
(6)
∂
2
u
(
x
,
t
)
∂
x
2
|
x
i
t
n
=
1
h
2
(
u
i
+
1
n
-
1
-
2
u
i
n
-
1
+
u
i
-
1
n
-
1
)
+
◯
(
h
2
)
.The explicit finite difference approximation for (1) is
(7)
1
τ
Γ
(
1
-
α
)
[
b
0
u
i
n
-
∑
k
=
1
n
-
1
‍
(
b
n
-
k
-
1
-
b
n
-
k
)
u
i
k
-
b
n
-
1
u
i
0
]
+
μ
u
i
n
=
u
i
+
1
n
-
1
-
2
u
i
n
-
1
+
u
i
-
1
n
-
1
h
2
+
K
φ
i
n
.Defines
=
b
0
+
μ
τ
Γ
(
1
-
α
), a
1
=
τ
Γ
(
1
-
α
)
/
(
s
h
2
), a
2
=
K
τ
Γ
(
1
-
α
)
/
s, U
n
=
(
u
1
n
,
u
2
n
,
…
,
u
M
n
)
T, F
n
=
(
f
1
n
,
f
2
n
,
…
,
f
M
n
)
T, f
i
n
=
τ
Γ
(
1
-
α
)
K
φ
i
n, and r
l as
(8)
r
l
=
b
l
-
b
l
+
1
s
.Equation (7) evolves as
(9)
U
n
=
∑
k
=
1
n
-
1
‍
r
n
-
1
-
k
U
k
+
b
n
-
1
U
0
+
A
U
n
-
1
+
a
2
F
n
,
where matrix A is a tridiagonal matrix, defined by
(10)
A
M
×
M
=
(
-
2
a
1
a
1
a
1
-
2
a
1
a
1
•
•
•
•
a
1
a
1
-
2
a
1
)
.


### 2.2. GPU Architecture and CUDA Programming Model

The architecture of GPU is optimized for rendering real-time graphics, a computation and memory access intensive problem domain with enormous inherent parallelism. Not like CPU, a much larger portion of GPUs resources is devoted to data processing rather than to caching or control flow. The NVIDIA Quadro FX 5800 GPU has 30 multiprocessor units which are called the streaming multiprocessors (SMs). Each SM contains 8 SIMD CUDA cores. Each CUDA core runs at 1.30 GHz. The multiprocessors create, manage, and execute concurrent threads in hardware with near zero scheduling overhead. The single instruction multiple thread (SIMT) unit, which is akin to SIMD vector organizations, creates, manages, schedules, and executes threads in groups of 32 parallel threads called warp [23].The programming model is a bridge between hardware and application. CUDA comes with a software environment that allows developers to use C as a high-level programming language. There are three key abstractions in CUDA: a hierarchy of thread execution model (grid or kernel, thread block, and thread), shared memory, and barrier synchronization. These abstractions provide fine-grained data level parallelism and thread parallelism, nested within coarse-grained data level parallelism and task level parallelism. Each CUDA kernel creates a single grid, which consists of many thread blocks. Each thread block schedules groups of threads that can share data through on-chip shared memory and synchronize with one another using barrier synchronization. Threads within a block are organized into warps which run in SIMT fashion. CUDA threads may access data from multiple memory spaces during their execution. The memory spaces include global, texture, and constant memory for grid, on-chip shared memory for thread block, and private register files for thread. The memory access time to different memory spaces varies from several to hundreds of cycles. These memory accesses must be coalesced to improve the performance of global memory access.


## 2.1. Numerical Solution

The fractional derivative off
(
t
) in the Caputo sense is defined as [27]
(2)
0
C
D
t
α
f
(
t
)
=
1
Γ
(
1
-
α
)
∫
0
t
‍
f
′
(
ξ
)
(
t
-
ξ
)
α
d
ξ
(
0
<
α
<
1
)
.Iff
′
(
t
) is continuous bounded derivatives in [
0
,
T
] for every T
>
0, we can get
(3)
0
C
D
t
α
f
(
t
)
=
lim
⁡
ξ
→
0
,
n
ξ
=
t
ξ
α
∑
i
=
0
n
‍
(
-
1
)
i
(
α
i
)
=
f
(
0
)
t
-
α
Γ
(
1
-
α
)
+
1
Γ
(
1
-
α
)
∫
0
t
‍
f
′
(
ξ
)
(
t
-
ξ
)
α
d
ξ
.Defineτ
=
T
/
N, h
=
x
R
/
(
M
+
1
), t
n
=
n
τ, and x
i
=
0
+
i
h for 0
≤
n
≤
N
,
0
≤
i
≤
M
+
1. Define u
i
n, φ
i
n, and ϕ
i as the numerical approximation to u
(
x
i
,
t
n
), f
(
x
i
,
t
n
), and ϕ
(
x
i
). We can get [11]
(4)
0
C
D
t
α
u
(
x
,
t
)
|
x
i
t
n
=
1
τ
Γ
(
1
-
α
)
[
b
0
u
i
n
-
∑
k
=
1
n
-
1
‍
(
b
n
-
k
-
1
-
b
n
-
k
)
u
i
k
-
b
n
-
1
u
i
0
]
+
◯
(
τ
2
-
α
)
,
where 1
≤
i
≤
M, n
≥
1, and
(5)
b
l
=
τ
1
-
α
1
-
α
[
(
l
+
1
)
1
-
α
-
l
1
-
α
]
,
l
≥
0
.Using explicit center difference scheme for∂
2
u
(
x
,
t
)
/
∂
x
2, we can get
(6)
∂
2
u
(
x
,
t
)
∂
x
2
|
x
i
t
n
=
1
h
2
(
u
i
+
1
n
-
1
-
2
u
i
n
-
1
+
u
i
-
1
n
-
1
)
+
◯
(
h
2
)
.The explicit finite difference approximation for (1) is
(7)
1
τ
Γ
(
1
-
α
)
[
b
0
u
i
n
-
∑
k
=
1
n
-
1
‍
(
b
n
-
k
-
1
-
b
n
-
k
)
u
i
k
-
b
n
-
1
u
i
0
]
+
μ
u
i
n
=
u
i
+
1
n
-
1
-
2
u
i
n
-
1
+
u
i
-
1
n
-
1
h
2
+
K
φ
i
n
.Defines
=
b
0
+
μ
τ
Γ
(
1
-
α
), a
1
=
τ
Γ
(
1
-
α
)
/
(
s
h
2
), a
2
=
K
τ
Γ
(
1
-
α
)
/
s, U
n
=
(
u
1
n
,
u
2
n
,
…
,
u
M
n
)
T, F
n
=
(
f
1
n
,
f
2
n
,
…
,
f
M
n
)
T, f
i
n
=
τ
Γ
(
1
-
α
)
K
φ
i
n, and r
l as
(8)
r
l
=
b
l
-
b
l
+
1
s
.Equation (7) evolves as
(9)
U
n
=
∑
k
=
1
n
-
1
‍
r
n
-
1
-
k
U
k
+
b
n
-
1
U
0
+
A
U
n
-
1
+
a
2
F
n
,
where matrix A is a tridiagonal matrix, defined by
(10)
A
M
×
M
=
(
-
2
a
1
a
1
a
1
-
2
a
1
a
1
•
•
•
•
a
1
a
1
-
2
a
1
)
.


## 2.2. GPU Architecture and CUDA Programming Model

The architecture of GPU is optimized for rendering real-time graphics, a computation and memory access intensive problem domain with enormous inherent parallelism. Not like CPU, a much larger portion of GPUs resources is devoted to data processing rather than to caching or control flow. The NVIDIA Quadro FX 5800 GPU has 30 multiprocessor units which are called the streaming multiprocessors (SMs). Each SM contains 8 SIMD CUDA cores. Each CUDA core runs at 1.30 GHz. The multiprocessors create, manage, and execute concurrent threads in hardware with near zero scheduling overhead. The single instruction multiple thread (SIMT) unit, which is akin to SIMD vector organizations, creates, manages, schedules, and executes threads in groups of 32 parallel threads called warp [23].The programming model is a bridge between hardware and application. CUDA comes with a software environment that allows developers to use C as a high-level programming language. There are three key abstractions in CUDA: a hierarchy of thread execution model (grid or kernel, thread block, and thread), shared memory, and barrier synchronization. These abstractions provide fine-grained data level parallelism and thread parallelism, nested within coarse-grained data level parallelism and task level parallelism. Each CUDA kernel creates a single grid, which consists of many thread blocks. Each thread block schedules groups of threads that can share data through on-chip shared memory and synchronize with one another using barrier synchronization. Threads within a block are organized into warps which run in SIMT fashion. CUDA threads may access data from multiple memory spaces during their execution. The memory spaces include global, texture, and constant memory for grid, on-chip shared memory for thread block, and private register files for thread. The memory access time to different memory spaces varies from several to hundreds of cycles. These memory accesses must be coalesced to improve the performance of global memory access.


## 3. Details of GPU Solution

The parallel solution consists of three parts. The first part is preprocessing, which prepares the initial matrices, vectors, and so on. The second part is the parallel solver, which focuses on the iteration of time step with (9). The third part is postprocessing, which outputs the final results and so on.The preprocessing includes initialization of parallel environment, distribution of computing task, allocation of local memory space, and initialization of variables and arrays. MatricesA
M
×
M and F
M
×
N are prepared before the computation of (9). For example, matrix A can be got according to (10). The postprocessing is simple. The results of the exact solution are performed. The max absolute error of the exact and parallel solutions is computed and outputted. Both the results of the exact and parallel solution are saved in files which are necessary for plot. Other operations of postprocessing include free memory space and stop the parallel environment.In order to getU
n, the right-sided computation of (9) should be performed. There are mainly one tridiagonal matrix vector multiplication, many constant vector multiplications, and many vector-vector additions in the right-sided computation.(1)
The tridiagonal matrix vector multiplication isA
U
n.
(2)
The constant vector multiplications areV
k
=
r
n
-
1
-
k
U
k, a
2
F
n, and so on.
(3)
The vector-vector additions are∑
k
=
1
n
-
1
‍
V
k and so on.The parallel solution uses the data level parallelism of GPU architecture. The parallel solution with CUDA for (1) is described in Algorithm 1. The preprocessing involves lines 1 to 3. The parallel solver, which is the most time consuming procedure, involves lines 4 to 12. The postprocessing involves lines 13 to 14 and other additional operations are not shown in Algorithm 1. Because the time spent on the preprocessing and postprocessing is trivial when the number of time steps is big enough, the preprocessing and postprocessing time is omitted for the measured time. T
1 and T
2 are used to record the measured time of the parallel CPU and GPU solutions.Algorithm 1:Parallel solution for Caputo fractional reaction-diffusion equation with CUDA.
(1)    Init CUDA environment
(2)   Allocate GPU global memoryA
,
U
,
F
,
r
,
b
….
(3)   Init variables and arrays on GPU
(4)   record timeT
1
(5)   call kerneli
n
i
t
U
0
〈
〈
〈
M
/
BS
,
BS
〉
〉
〉
(
…
)
(6)for  n
=
1  to  N  by Step  1
do
(7)               call kernelt
r
i
m
v
m
C
U
〈
〈
〈
M
/
BS
,
BS
〉
〉
〉
(
…
)
(8)               call kernelc
v
m
_
v
v
a
C
U
〈
〈
〈
M
/
BS
,
BS
〉
〉
〉
(
…
)
(9)for  k
=
1  to  n  by Step  1
do
(10)                     call kernelc
v
m
a
C
U
〈
〈
〈
M
/
BS
,
BS
〉
〉
〉
(
…
)
(11)              call kernelc
v
m
C
U
〈
〈
〈
M
/
BS
,
BS
〉
〉
〉
(
…
)
(12) record timeT
2
(13) outputT
2
-
T
1 and U
N
…
(14) free GPU memory and stop CUDA environmentThe parallel solution uses the data level parallelism of GPU architecture. The parallel solution with CUDA for (1) is described in Algorithm 1. The preprocessing involves lines 1 to 3. The parallel solver, which is the most time consuming procedure, involves lines 4 to 12. The postprocessing involves lines 13 to 14 and other additional operations are not shown in Algorithm 1. BS stands for the CUDA thread block size and M/BS is the number of CUDA thread blocks. BS is the predefined constant like 16, 32, 64, and so forth. Because the time spent on the preprocessing and postprocessing is trivial when the number of time steps is big enough, the preprocessing and postprocessing time is omitted for the measured time. T
1 and T
2 are used to record the measured time of the serial and parallel solution.Except for the initialization of variables and arrays, there are four CUDA kernels. The first kernel isi
n
i
t
U
0, which computes the initial condition according to ϕ
(
x
) in (1). The second kernel is t
r
i
m
v
m
C
U, which performs the tridiagonal matrix vector multiplication. The third kernel is c
v
m
_
v
v
a
C
U, which stands for constant vector multiplication and vector-vector addition. The fourth kernel is c
v
m
a
C
U, which performs the constant vector multiplication and vector-vector addition. The last kernel is c
v
m
C
U, which stands for constant vector multiplication. The CUDA kernels i
n
i
t
U
1 and c
v
m
C
U are simple and will not be described in detail.


### 3.1. Implementation

The CUDA kernelc
v
m
a
C
U for constant vector multiplication and vector-vector addition is shown in Algorithm 2. Algorithm 2 computes U
c
o
l
O
u
t
+
a
2
U
c
o
l
I
n and saves the final vector into U
c
o
l
O
u
t.Algorithm 2:CUDA kernel for constant vector multiplication and vector vector addition.
input: c
o
l
O
u
t
,
c
o
l
I
n
,
a
2
,
U
output: U
(1)for  all threads in each thread block
do in parallel
(2)         local variableg
i
d
,
ξ
(3)g
i
d
←
t
h
r
e
a
d
I
d
.
x
+
b
l
o
c
k
I
d
.
x
·
b
l
o
c
k
D
i
m
.
x
(4)ξ
←
a
2
×
U
[
M
·
c
o
l
I
n
+
g
i
d
]
(5)U
[
M
·
c
o
l
O
u
t
+
g
i
d
]
←
U
[
M
·
c
o
l
O
u
t
+
g
i
d
]
+
ξMost elements of tridiagonal matrixA are zero. The most common data structure used to store a sparse matrix for sparse matrix vector multiplication computations is compressed sparse row (CSR) format [34] shown in
(11)
A
3
×
M
=
(
0
a
1
⋯
a
1
a
1
-
2
a
1
-
2
a
1
⋯
-
2
a
1
-
2
a
1
a
1
a
1
⋯
a
1
0
)
.
So in the following parts of this paper, matrix A stands for the format of (11) not the format of (10). With the format of (11), the global memory is coalesced and can improve the performance of global memory access.The CUDA kernelt
r
i
m
v
m
C
U for tridiagonal matrix vector multiplication is shown in Algorithm 3. One thread block deals with the multiplication of one row of matrix A and one column of U. Algorithm 3 computes A
U
n
-
1 and saves the final vector into U
n. The shared memory is used to improve the memory access speed. The synchronization function _
_
s
y
n
c
t
h
r
e
a
d
s
(
) is used to ensure the correctness of the logic of the algorithm.Algorithm 3:CUDA kernel for tridiagonal matrix vector multiplication.
input: n
,
M
,
BS
,
A
,
U
output: U
(1)    shared memorys
m
[
BS
+
2
]
(2)for  all threads in every thread block
do in parallel
(3)                local variablest
i
d
,
g
i
d
,
BS
,
ξ
(4)t
i
d
←
t
h
r
e
a
d
I
d
.
x
(5)g
i
d
←
t
h
r
e
a
d
I
d
.
x
+
b
l
o
c
k
I
d
.
x
·
b
l
o
c
k
D
i
m
.
x
(6)s
m
[
t
i
d
+
1
]
←
U
[
M
·
(
n
-
1
)
+
g
i
d
]
(7)if  0
=
=
t
i
d  then
(8)if  0
=
=
g
i
d  then
(9)s
m
[
0
]
←
0
(10)else
(11)s
m
[
0
]
←
U
[
M
·
(
n
-
1
)
+
g
i
d
-
1
]
(12)if  BS
-
1
=
=
t
i
d  then
(13)if  M
=
=
g
i
d  then
(14)s
m
[
BS
+
1
]
←
0
(15)else
(16)s
m
[
BS
+
1
]
←
U
[
M
·
(
n
-
1
)
+
g
i
d
+
1
]
(17)_
_
s
y
n
c
t
h
r
e
a
d
s
(
)
(18)ξ
←
0.0
(19)ξ
←
ξ
+
s
m
[
t
i
d
+
0
]
·
A
[
0
*
M
+
g
i
d
]
(20)ξ
←
ξ
+
s
m
[
t
i
d
+
1
]
·
A
[
1
*
M
+
g
i
d
]
(21)ξ
←
ξ
+
s
m
[
t
i
d
+
2
]
·
A
[
2
*
M
+
g
i
d
]
(22)U
[
M
·
n
+
g
i
d
]
←
U
[
M
·
n
+
g
i
d
]
+
ξIn Algorithm3, each GPU thread deals with the multiplication of one row of tridiagonal matrix A and vector U
n
-
1. Each thread needs to use three elements of vector U
n. The nearest two threads use the same two elements of vector U
n
-
1. We can use the shared memory to improve the memory access performance. So the elements of vector U
n
-
1 which will be used by threads in a thread block are stored into shared memory, shown between lines 6 and 16 of Algorithm 3. The real computation of tridiagonal matrix multiplication is shown between lines 18 and 21. Finally, the results are stored into the global memory of U
n.


### 3.2. Optimization

In Algorithm1, the kernels c
v
m
_
v
v
a
C
U, c
v
m
C
U, and t
r
i
m
v
m
C
U are invoked N times. In each time step, kernel c
v
m
a
C
U is invoked n (1,2
,
…
,
N) times. Because
(12)
I
(
I
-
1
)
2
=
1
+
2
+
⋯
I
the total number of the invocations of kernel c
v
m
a
C
U is N
2
(
N
-
1
)
/
2. The most time consuming part of Algorithm 1 is the loop of line 9. Loop 9 can be combined into one CUDA kernel c
v
m
a
O
p
t
C
U as shown in Algorithm 4. The array g
p
u
R is the coefficient of (8) in global memory.Algorithm 4:Optimized CUDA kernel for constant vector multiplication and vector-vector addition.
input: n
,
g
p
u
R
,
U
output: U
(1)for  all threads in every thread block
do in parallel
(2)           local variablesg
i
d
,
ξ
,
k
(3)g
i
d
←
t
h
r
e
a
d
I
d
.
x
+
b
l
o
c
k
I
d
.
x
·
b
l
o
c
k
D
i
m
.
x
(4)ξ
←
0
(5)for  k
=
1  to  n
-
1  by Step  1
do
(6)ξ
←
ξ
+
g
p
u
R
[
k
]
·
U
[
M
·
(
n
-
1
-
k
)
+
g
i
d
]
(7)U
[
M
·
n
+
g
i
d
]
←
U
[
M
·
n
+
g
i
d
]
+
ξSo the optimized parallel solution for Caputo fractional reaction-diffusion equation is similar to Algorithm1 except that the loop (lines 9-10) in Algorithm 1 is replaced with the optimized CUDA kernel of Algorithm 4.


## 3.1. Implementation

The CUDA kernelc
v
m
a
C
U for constant vector multiplication and vector-vector addition is shown in Algorithm 2. Algorithm 2 computes U
c
o
l
O
u
t
+
a
2
U
c
o
l
I
n and saves the final vector into U
c
o
l
O
u
t.Algorithm 2:CUDA kernel for constant vector multiplication and vector vector addition.
input: c
o
l
O
u
t
,
c
o
l
I
n
,
a
2
,
U
output: U
(1)for  all threads in each thread block
do in parallel
(2)         local variableg
i
d
,
ξ
(3)g
i
d
←
t
h
r
e
a
d
I
d
.
x
+
b
l
o
c
k
I
d
.
x
·
b
l
o
c
k
D
i
m
.
x
(4)ξ
←
a
2
×
U
[
M
·
c
o
l
I
n
+
g
i
d
]
(5)U
[
M
·
c
o
l
O
u
t
+
g
i
d
]
←
U
[
M
·
c
o
l
O
u
t
+
g
i
d
]
+
ξMost elements of tridiagonal matrixA are zero. The most common data structure used to store a sparse matrix for sparse matrix vector multiplication computations is compressed sparse row (CSR) format [34] shown in
(11)
A
3
×
M
=
(
0
a
1
⋯
a
1
a
1
-
2
a
1
-
2
a
1
⋯
-
2
a
1
-
2
a
1
a
1
a
1
⋯
a
1
0
)
.
So in the following parts of this paper, matrix A stands for the format of (11) not the format of (10). With the format of (11), the global memory is coalesced and can improve the performance of global memory access.The CUDA kernelt
r
i
m
v
m
C
U for tridiagonal matrix vector multiplication is shown in Algorithm 3. One thread block deals with the multiplication of one row of matrix A and one column of U. Algorithm 3 computes A
U
n
-
1 and saves the final vector into U
n. The shared memory is used to improve the memory access speed. The synchronization function _
_
s
y
n
c
t
h
r
e
a
d
s
(
) is used to ensure the correctness of the logic of the algorithm.Algorithm 3:CUDA kernel for tridiagonal matrix vector multiplication.
input: n
,
M
,
BS
,
A
,
U
output: U
(1)    shared memorys
m
[
BS
+
2
]
(2)for  all threads in every thread block
do in parallel
(3)                local variablest
i
d
,
g
i
d
,
BS
,
ξ
(4)t
i
d
←
t
h
r
e
a
d
I
d
.
x
(5)g
i
d
←
t
h
r
e
a
d
I
d
.
x
+
b
l
o
c
k
I
d
.
x
·
b
l
o
c
k
D
i
m
.
x
(6)s
m
[
t
i
d
+
1
]
←
U
[
M
·
(
n
-
1
)
+
g
i
d
]
(7)if  0
=
=
t
i
d  then
(8)if  0
=
=
g
i
d  then
(9)s
m
[
0
]
←
0
(10)else
(11)s
m
[
0
]
←
U
[
M
·
(
n
-
1
)
+
g
i
d
-
1
]
(12)if  BS
-
1
=
=
t
i
d  then
(13)if  M
=
=
g
i
d  then
(14)s
m
[
BS
+
1
]
←
0
(15)else
(16)s
m
[
BS
+
1
]
←
U
[
M
·
(
n
-
1
)
+
g
i
d
+
1
]
(17)_
_
s
y
n
c
t
h
r
e
a
d
s
(
)
(18)ξ
←
0.0
(19)ξ
←
ξ
+
s
m
[
t
i
d
+
0
]
·
A
[
0
*
M
+
g
i
d
]
(20)ξ
←
ξ
+
s
m
[
t
i
d
+
1
]
·
A
[
1
*
M
+
g
i
d
]
(21)ξ
←
ξ
+
s
m
[
t
i
d
+
2
]
·
A
[
2
*
M
+
g
i
d
]
(22)U
[
M
·
n
+
g
i
d
]
←
U
[
M
·
n
+
g
i
d
]
+
ξIn Algorithm3, each GPU thread deals with the multiplication of one row of tridiagonal matrix A and vector U
n
-
1. Each thread needs to use three elements of vector U
n. The nearest two threads use the same two elements of vector U
n
-
1. We can use the shared memory to improve the memory access performance. So the elements of vector U
n
-
1 which will be used by threads in a thread block are stored into shared memory, shown between lines 6 and 16 of Algorithm 3. The real computation of tridiagonal matrix multiplication is shown between lines 18 and 21. Finally, the results are stored into the global memory of U
n.


## 3.2. Optimization

In Algorithm1, the kernels c
v
m
_
v
v
a
C
U, c
v
m
C
U, and t
r
i
m
v
m
C
U are invoked N times. In each time step, kernel c
v
m
a
C
U is invoked n (1,2
,
…
,
N) times. Because
(12)
I
(
I
-
1
)
2
=
1
+
2
+
⋯
I
the total number of the invocations of kernel c
v
m
a
C
U is N
2
(
N
-
1
)
/
2. The most time consuming part of Algorithm 1 is the loop of line 9. Loop 9 can be combined into one CUDA kernel c
v
m
a
O
p
t
C
U as shown in Algorithm 4. The array g
p
u
R is the coefficient of (8) in global memory.Algorithm 4:Optimized CUDA kernel for constant vector multiplication and vector-vector addition.
input: n
,
g
p
u
R
,
U
output: U
(1)for  all threads in every thread block
do in parallel
(2)           local variablesg
i
d
,
ξ
,
k
(3)g
i
d
←
t
h
r
e
a
d
I
d
.
x
+
b
l
o
c
k
I
d
.
x
·
b
l
o
c
k
D
i
m
.
x
(4)ξ
←
0
(5)for  k
=
1  to  n
-
1  by Step  1
do
(6)ξ
←
ξ
+
g
p
u
R
[
k
]
·
U
[
M
·
(
n
-
1
-
k
)
+
g
i
d
]
(7)U
[
M
·
n
+
g
i
d
]
←
U
[
M
·
n
+
g
i
d
]
+
ξSo the optimized parallel solution for Caputo fractional reaction-diffusion equation is similar to Algorithm1 except that the loop (lines 9-10) in Algorithm 1 is replaced with the optimized CUDA kernel of Algorithm 4.


## 4. Experimental Results

### 4.1. Experiment Platforms

The experiment platforms consist of one GPU and one CPU listed in Table1. For the purpose of fair comparisons, we measure the performance provided by GPU compared to the MPI code running on multicore CPU [31]. Both codes run on double precision floating point operations. The CPU code is compiled by the mpif90 compiler with level three optimization. The GPU code is compiled by the NVCC compiler provided by CUDA version 3.0 with level three optimization too.Table 1
Technical specifications of experiment platforms.
CPU
Intel Xeon E5540,4 cores, 2.53 GHz
GPU
NVIDIA Quadro FX 5800,240 SPs, 1.30 GHz
Operating system
Kylin server version 3.1
CPU compilerGPU compiler
mpif90, Intel Fortran version 11.1NVCC, CDUA version 3.0
Communication
MPICH2 version 1.3rc2


### 4.2. Numerical Example

The following Caputo fractional reaction-diffusion equation [11] was considered:
(13)
0
D
t
α
u
(
x
,
t
)
+
μ
u
(
x
,
t
)
=
∂
2
u
(
x
,
t
)
∂
x
2
+
K
f
(
x
,
t
)

(
0
<
α
<
1
)
u
(
x
,
0
)
=
0
,
x
∈
(
0,2
)
u
(
0
,
t
)
=
u
(
2
,
t
)
=
0
with μ
=
1
,
K
=
1, and
(14)
f
(
x
,
t
)
=
2
Γ
(
2.3
)
x
(
2
-
x
)
t
1.3
+
x
(
2
-
x
)
t
2
+
2
t
2
.
The exact solution of (13) is
(15)
u
(
x
,
t
)
=
x
(
2
-
x
)
t
2
.


### 4.3. Accuracy of the GPU Implementation

The GPU solution compares well with the exact solution to the time fractional diffusion equation in this test case of (13), shown in Figure 1. The τ and h for the GPU solution are T
/
2048 and 2.0
/
16. The maximum absolute errors for T
=
0.3, 0.5, and 0.7 are 3.29
×
1
0
-
5, 1.07
×
1
0
-
4, and 2.30
×
1
0
-
4. In fact, the accuracy and convergence of the GPU solution are the same as the serial and parallel MPI solution on CPU [31].Figure 1
Comparison of exact solution to the parallel GPU solution atT
=
0.3
,
0.5
,
0.7.


### 4.4. Total Performance Improvement

In this section, the performance of the optimized GPU solution presented in this paper is compared with the performance of the parallel CPU solution [31]. The parallel CPU solution makes full use of four cores of E5540. The optimized GPU solution is presented in Section 3.2.For fixedN
=
128, the performance comparison between GPU and multicore CPU is shown in Table 2. The thread block size is 32.Table 2
Performance comparison between optimized GPU solution on Quadro FX 5800 and parallel CPU solution on E5540 with fixedN
=
128.
M
CPU
GPU
Speedup
245760
0.58
0.46
1.26
491520
1.23
0.90
1.36
737280
1.95
1.27
1.54
983040
3.27
1.89
1.73
1228800
4.66
2.06
2.26For fixedM
=
122880, the performance comparison between optimized GPU solution and parallel CPU solution is shown in Table 3. The thread block size is 32.Table 3
Performance comparison between optimized GPU solution on Quadro FX 5800 and parallel CPU solution on E5540 with fixedM
=
122880.
N
CPU
GPU
Speedup
128
0.29
0.25
1.14
256
1.06
0.84
1.26
512
3.96
2.82
1.40
1024
15.20
9.33
1.63
2048
60.16
33.59
1.79


### 4.5. Performance Issues of GPU Solution

WithM
=
491520, N
=
128, and thread block size 64, the runtime of Algorithm 1 on Quadro FX 5800 is 1.228 seconds. Without the loop of line 9 in Algorithm 1, the runtime is only 0.032 seconds. That is to say, about 97.39% of runtime is spent on the loop of line 9. This is the reason why we develop the optimized GPU solution with an optimized CUDA kernel of Algorithm 4.The impact of the optimized CUDA kernel on constant vector multiplication and vector-vector addition with fixedN
=
128 is shown in Table 4. The performance improvement with fixed M
=
491520 is shown in Table 5. All thread block sizes are 64. The basic GPU solution is Algorithm 1 and the optimized GPU solution uses the optimized CUDA kernel of Algorithm 4.Table 4
Performance improvement for fixedN with the optimization of constant vector multiplication and vector-vector addition.
M
Original
Optimization
Speedup
245760
0.64
0.53
1.21
491520
1.26
1.05
1.20
737280
1.85
1.56
1.19
983040
2.53
2.10
1.20
1228800
3.06
2.59
1.18Table 5
Performance improvement for fixedM with the optimization of constant vector multiplication and vector-vector addition.
N
Original
Optimization
Speedup
256
0.72
0.53
1.36
512
2.86
1.75
1.63
1024
11.36
5.60
2.03
2048
45.32
18.65
2.43
4096
181.01
66.71
2.71The thread block size (BS) is a key parameter for parallel GPU algorithms. The impact of BS is shown in Table6. From Table 6, we can see that thread block size 32 is the best choice.Table 6
Impact of thread block size (BS).
BS
Runtime
BS
Runtime
4
13.5
64
3.35
8
6.84
128
3.61
16
3.59
256
3.55
32
2.82
512
3.47


## 4.1. Experiment Platforms

The experiment platforms consist of one GPU and one CPU listed in Table1. For the purpose of fair comparisons, we measure the performance provided by GPU compared to the MPI code running on multicore CPU [31]. Both codes run on double precision floating point operations. The CPU code is compiled by the mpif90 compiler with level three optimization. The GPU code is compiled by the NVCC compiler provided by CUDA version 3.0 with level three optimization too.Table 1
Technical specifications of experiment platforms.
CPU
Intel Xeon E5540,4 cores, 2.53 GHz
GPU
NVIDIA Quadro FX 5800,240 SPs, 1.30 GHz
Operating system
Kylin server version 3.1
CPU compilerGPU compiler
mpif90, Intel Fortran version 11.1NVCC, CDUA version 3.0
Communication
MPICH2 version 1.3rc2


## 4.2. Numerical Example

The following Caputo fractional reaction-diffusion equation [11] was considered:
(13)
0
D
t
α
u
(
x
,
t
)
+
μ
u
(
x
,
t
)
=
∂
2
u
(
x
,
t
)
∂
x
2
+
K
f
(
x
,
t
)

(
0
<
α
<
1
)
u
(
x
,
0
)
=
0
,
x
∈
(
0,2
)
u
(
0
,
t
)
=
u
(
2
,
t
)
=
0
with μ
=
1
,
K
=
1, and
(14)
f
(
x
,
t
)
=
2
Γ
(
2.3
)
x
(
2
-
x
)
t
1.3
+
x
(
2
-
x
)
t
2
+
2
t
2
.
The exact solution of (13) is
(15)
u
(
x
,
t
)
=
x
(
2
-
x
)
t
2
.


## 4.3. Accuracy of the GPU Implementation

The GPU solution compares well with the exact solution to the time fractional diffusion equation in this test case of (13), shown in Figure 1. The τ and h for the GPU solution are T
/
2048 and 2.0
/
16. The maximum absolute errors for T
=
0.3, 0.5, and 0.7 are 3.29
×
1
0
-
5, 1.07
×
1
0
-
4, and 2.30
×
1
0
-
4. In fact, the accuracy and convergence of the GPU solution are the same as the serial and parallel MPI solution on CPU [31].Figure 1
Comparison of exact solution to the parallel GPU solution atT
=
0.3
,
0.5
,
0.7.


## 4.4. Total Performance Improvement

In this section, the performance of the optimized GPU solution presented in this paper is compared with the performance of the parallel CPU solution [31]. The parallel CPU solution makes full use of four cores of E5540. The optimized GPU solution is presented in Section 3.2.For fixedN
=
128, the performance comparison between GPU and multicore CPU is shown in Table 2. The thread block size is 32.Table 2
Performance comparison between optimized GPU solution on Quadro FX 5800 and parallel CPU solution on E5540 with fixedN
=
128.
M
CPU
GPU
Speedup
245760
0.58
0.46
1.26
491520
1.23
0.90
1.36
737280
1.95
1.27
1.54
983040
3.27
1.89
1.73
1228800
4.66
2.06
2.26For fixedM
=
122880, the performance comparison between optimized GPU solution and parallel CPU solution is shown in Table 3. The thread block size is 32.Table 3
Performance comparison between optimized GPU solution on Quadro FX 5800 and parallel CPU solution on E5540 with fixedM
=
122880.
N
CPU
GPU
Speedup
128
0.29
0.25
1.14
256
1.06
0.84
1.26
512
3.96
2.82
1.40
1024
15.20
9.33
1.63
2048
60.16
33.59
1.79


## 4.5. Performance Issues of GPU Solution

WithM
=
491520, N
=
128, and thread block size 64, the runtime of Algorithm 1 on Quadro FX 5800 is 1.228 seconds. Without the loop of line 9 in Algorithm 1, the runtime is only 0.032 seconds. That is to say, about 97.39% of runtime is spent on the loop of line 9. This is the reason why we develop the optimized GPU solution with an optimized CUDA kernel of Algorithm 4.The impact of the optimized CUDA kernel on constant vector multiplication and vector-vector addition with fixedN
=
128 is shown in Table 4. The performance improvement with fixed M
=
491520 is shown in Table 5. All thread block sizes are 64. The basic GPU solution is Algorithm 1 and the optimized GPU solution uses the optimized CUDA kernel of Algorithm 4.Table 4
Performance improvement for fixedN with the optimization of constant vector multiplication and vector-vector addition.
M
Original
Optimization
Speedup
245760
0.64
0.53
1.21
491520
1.26
1.05
1.20
737280
1.85
1.56
1.19
983040
2.53
2.10
1.20
1228800
3.06
2.59
1.18Table 5
Performance improvement for fixedM with the optimization of constant vector multiplication and vector-vector addition.
N
Original
Optimization
Speedup
256
0.72
0.53
1.36
512
2.86
1.75
1.63
1024
11.36
5.60
2.03
2048
45.32
18.65
2.43
4096
181.01
66.71
2.71The thread block size (BS) is a key parameter for parallel GPU algorithms. The impact of BS is shown in Table6. From Table 6, we can see that thread block size 32 is the best choice.Table 6
Impact of thread block size (BS).
BS
Runtime
BS
Runtime
4
13.5
64
3.35
8
6.84
128
3.61
16
3.59
256
3.55
32
2.82
512
3.47


## 5. Conclusions and Future Work

In this paper, the numerical solution of Caputo fractional reaction-diffusion equation with explicit finite difference method is accelerated on GPU. The iteration of time step consists of tridiagonal matrix vector multiplication, constant vector multiplication, and vector-vector addition. The details of the GPU solution and some basic CUDA kernels are presented. The most time consuming loop (constant vector multiplication and vector-vector addition) is optimized. The experimental results show the GPU solution compares well with the exact analytic solution and is much faster than parallel CPU solution. So the power of parallel computing on GPU for solving fractional applications should be recognized.As a part of the future work, first, the stability, like Ulam’s type, of different fractional equations should be paid attention to [35–37]. Second, parallelizing the implicit numerical method of fractional differential equations is challenging.



---

*Source: 820162-2014-06-17.xml*