Parallel Algorithm For Constructing a Cubic Spline on Multi-Core Processors in a Cluster

Parallel Algorithm For Constructing a Cubic Spline

on Multi-Core Processors in a Cluster

1

st

Oybek Mallaev

Computer engineering. TUIT

Tashkent, Uzbekistan

https://orcid.org/0000-0002-1918-0641

2

nd

Javohir Nurmurodov

Computer engineering. TUIT

Tashkent, Uzbekistan

[e-mail]

Abstract—The relevance of using parallel computing systems

is reflected, the main approaches to parallelizing processes and

data processing methods are reflected, the principles of parallel

programming technologies are described, the main parameters

of parallel algorithms are studied using cubic splines calculation

example, a comparative analysis of the effectiveness of using par-

allelism technologies is given. The parallel algorithm considered

for constructing the cubic spline of defect 1 as p to n leads to

the construction of a local cubic spline on each grid interval ω.

Index Terms—Parallel computing, UMA, NUMA, SMP, data

parallelization, task parallelization, data processing, MPI.

I. I

NTRODUCTION

An alternative to a common physical address space can

be separate physical address spaces with a common virtual

address space and assignment of data exchange work to the

operating system. This approach was tested, but we had a very

high level of costs in providing programmers with a workable

abstraction of shared memory [1], [2]. Several attempts have

been made to create high-performance computers based on

high-speed messaging networks, and they have shown better

absolute data exchange performance than clusters based on

local area networks [3]. The problem was that they were

much more expensive. Only a few applications could justify

a significantly higher cost [4]. Therefore, it is clusters that

have become the most common example of parallel computers

based on messaging today [6], [8]. Clusters in general are a

set of mass-produced computers combined by I/O systems and

standard network switches and cables. Each of the computers

is running a separate copy of the operating system [9], [11].

Almost every Internet service is based on clusters assembled

from mass servers and switches.

II. B

UILDING A

C

UBIC

S

PLINE

Let on the segment [a, b] of the real axis Ox on the grid

sentence, as in:

ω : a = x

0

< x

1

< ... < x

n

= b, h

i

= x

i

− x

i−1

, i = 1, n,

given the values of some grid function f

i

, i=0,1,,n. For

the interpolation cubic spline S(x) of defect 1, the following

relations hold:

Identify applicable funding agency here. If none, delete this.

S(x) =

3

X

k=0

a

i

k

(x − x

i

, x ∈ [x

i

, x

x+1

]), i = 1, n

(1)

S(x) = f

i

, i = 1, n.

(2)

S(x) ∈ C

2

[a, b],

(3)

where the superscript i denotes the spline unit number. It

follows from (1) that on ω : the spline S(x) depends on 4 n

coefficients

a

i

k

(x − x

i

, i = 0, n − 1, k = 0, 3.

To determine them, there are 3 (n-1) continuity conditions (3)

at the internal nodes of the grid , as well as n+1 interpolation

condition (2). Thus, to determine the spline coefficients, two

more conditions must be specified. Four types of additional

conditions are usually considered:

S

0

(x) = f

0

(a), S

0

(b) = f

0

(b).

(4)

S

00

(x) = f

00

(a), S

00

(b) = f

00

(b).

(5)

S

0

(a) = S

0

(b), S

00

(a) = S

00

(b).

(6)

S

000

(x

1

− 0) = S

000

(x

1

+ 0),

S

000

(x

n−1

− 0) = S

000

(x

n−1

+ 0).

(7)

Additional conditions (6) are used for periodic splines with a

period (b - a). If for non-periodic spline there is no additional

information about the value of its first or second derivative

at x=a or x=b, then it is recommended to use additional

conditions (7). In this case

a

0

k

= a

1

k

, a

n−1

k

= a

n−2

k

, k = 0, 3.

We introduce the following notation:

M

i

= S

00

(x

i

), i = 0, n.

The quantities M

i

is called the moments of the spline S(x).

As follows from the definition of a cubic spline, its second

derivative is a linear function on each elementary interval of

the grid ω. For a link of such a broken line, the formula

S

00

(x) =

x

i

− x

h

i

M

i−1

+

x − x

i−1

h

i

M

i

, h

i

=

= x − x

i−1

, x ∈ [x

i−1

− x

i

].

(8)

Integrating (8) twice over x, we successively obtain

S

0

(x) =

(x

i

− x)

2

2h

i

M

i−1

+

(x − x

i−1

)

2

2h

i

M

i

+ C

i−1

1

, (9)

S(x) =

(x

i

− x)

3

6h

i

M

i−1

+

(x − x

i−1

)

3

6h

i

M

i

+ C

i−1

1

x + C

i−1

2

(10)

where C

i−1

1

, C

i−1

2

are the integration constants for the (i–

1) grid interval ω Substituting (10) and the interpolation

conditions (2), we obtain

S

0

(x) =

(x

i

− x)

2

2h

i

M

i−1

+

(x − x

i−1

)

2

2h

i

M

i

+

f − f

i−1

2h

i

−

h

i

6

(M

i

− M

i−1

),

(11)

S(x) =

(x

i

− x)

3

6h

i

M

i−1

+

(x − x

i−1

)

3

6h

i

M

i

+

(f

i−1

−

h

2

i

6

M

i−1

)

x

i

− x

h

i

+ (f

i

−

h

i

6

(M

i

)

x − x

i−1

h

i

,

∈ [x

i

, x].

(12)

Note that (11), (12) is a representation of the cubic spline

and its first derivative on the (i–1) interval of the grid ω. In

this interval, S

00

(x) is determined by function (8), and

s

00

(x) =

1

h

i

(M

i

− M

i−1

), x ∈ [x

i−1

− x

i

].

If in expressions (11), (12) we replace the index i by i+1, then

we obtain the representation S(x), S

0

(x) on the ith interval of

the grid ω:

S

0

(x) =

(x

i+1

− x)

2

2h

i+1

M

i

+

(x − x

i

)

2

2h

i+1

M

i+1

+

f

i+1

− f

i

h

i+1

−

h

i+1

6

(M

i+1

− M

i

),

(13)

S(x) =

(x

i+1

− x)

3

6h

i+1

M

i

+

(x − x

i

)

3

6h

i+1

M

i+1

+

(f

i

h

2

i+1

6

M

i

)

x

i+1

− x

h

i+1

+ (f

i+1

−

h

i+1

6

M

i+1

)

x − x

i

h

i+1

, x ∈ [x

i

, x

i+1

].

(14)

From (11) and (13) for x=x

i

we write down the one–sided

limits of the derivative:

S

0

(x

i

− 0) =

h

i

6

M

i−1

+

h

i

3

M

i

+

f

i

− f

i−1

h

i

, (15)

S

0

(x

i

− 0) =

h

i+1

6

M

i+1

+

h

i+1

3

M

i

+

f

i+1

− f

i

h

i+1

, (16)

From the continuity condition for the first derivative of the

spline at the internal nodes of the grid ω, we obtain a system

of equations for determining the moments of the cubic spline,

for example, in the case of additional conditions of the first

type

2M

0

+ M

1

=

6

h

i



f

0

(a) −

f

1

− f

0

h

1



,

a

i

M

i−1

+ aM

i

+ β

i

M

i+1

= 3f

¯

xx

, i = 1, n − 1,

2M

n

+ M

n−1

=

6

h

n



f

0

(b) −

f

n

− f

n−1

h

n



,

(17)

where

h

i

=

h

i

+ h

i+1

2

, α

i

=

h

i

2h

i

, β

i

=

h

i+1

2h

i

,

f

xx

=

1

h

i

 f

i+1

− f

i

h

i+1

−

f

i

− f

i−1

h

i



.

Having determined the array of moments by the sweep

method from this system and substituting , M

i

into expression

(14), we obtain the representation S(x) for X. Note that if,

when calculating the integral,

I =

Z

b

a

f (x)dx

replace the integrate f(x) with the cubic spline S(x), then we

obtain the formula for numerical integration

Z

b

a

S(x)dx =

1

2

n−1

X

i=0

h

i+1

(f

i

+ f

i+1

) −

1

24

n−1

X

i=0

h

3

i+1

(M

i

+ M

i+1

),

which has a higher accuracy than the trapezoid formula. We

now introduce the expanded grid

ω : x

−3

< x

−2

< ... < x

0

... < x

n

< x

n+1

< ... < x

n+3

.

and normalized basis functions B

i

(x), i = –1.0, ..., n, n+1,

where x

i

is the center of the carrier. Considering that these

functions form a basis in the space of cubic splines of defect

1, the spline S(x) from this set can be represented as

S(x) =

n+1

X

i=−1

b

i

B

i

(x),

(18)

where the b

i

coefficients are to be determined. Then, for ex-

ample, the system of equations for determining the coefficients

of the interpolation spline (18) for the additional conditions of

the first type (4) has the form







b

−1

B

0

−1

(x

0

) + b

0

B

0

0

(x

0

) + b

1

B

0

1

(x

0

) = f

0

0

,

b

i−1

B

i−1

(x

i

) + b

i

B

i

(x

i

) + b

i+1

B

i+1

(x

i

) = f

i

, i = 0, n,

b

n−1

B

0

n−1

(x

n

) + b

n

B

0

n

(x

n

) + b

n+1

B

0

n+1

(x

n

) = f

0

n

,

After eliminating b and b

n+i

, such a system can be resolved

by a sweep method if

max

| {z }

|i=j|−1

h

i

h

j

< (1 +

√

13)/2 ≈ 2.3,

which corresponds to the condition of diagonal dominance of

the matrix. When calculating B

i

(x) for x ∈ (x

i

, x

i+1

), one

can use the formulas:

B

i−1

(x) =

1

6

(1 − t)

3

,

B

i

(x) =

1

6

1 + 3(1 − t) + 3t(1 − t)

2

 ,

B

i+1

(x) =

1

6

1 + 3t + 3t

2

(1 − t)

 , B

i+2

(x) =

1

6

t

3

,

t =

x − x

i

h

, t ∈ [0, 1], h = (b − a)/n.

Note that representing a spline in terms of B–functions (18) is

convenient for large n, since its calculation requires storage of

two arrays of information (nodes x

i

and coefficients b

i

). On the

other hand, a change in one of the coefficients in (18) implies

a change in S(x) only on the support of the corresponding B–

function, which can be used to change the spline locally.

A significant number of applied problems are associated with

finding eigenvalues of large–dimensional matrices. Let A =

[a

ij

]

n

i

– is a real symmetric matrix. Then, using Gershgorin’s

theorem, we can determine the interval [a,b] containing all real

eigenvalues λ

i

of this matrix. If on [a,b] introduce a uniform

grid ω and associate a grid function with it

f

i

= det(A − ξ

i

E), ξ

i

= a + i · h, h = (b − a)/n, i = 0, n,

then the problem of determining the eigenvalues by the

spline method reduces to finding the roots of the cubic spline

(polynomial of the third degree) on each interval of the grid ω.

Using the decomposition of the grid domain by the number

of processor elements p, the analysis of the roots of cubic

equations can be performed in parallel mode after loading the

corresponding values of M

i

and f

i

into the processor memory.

To control the calculations, you can use the known relations

Sp(A) =

n

X

i=1

a

ij

=

n

X

i=1

λ

i

, det(A) = ±

n

Y

i=1

λ

i

.

Note that when determining the real roots of the equation

using a cubic spline, the grid step h must be selected so that

there are no more than three roots on each interval.

III. P

ARALLEL

A

LGORITHM

F

OR

C

ONSTRUCTING A

C

UBIC

S

PLINE

Suppose that values of some grid function f

0

, f

1

, ..., f

n

and

n >> p >> 1 are given on a uniform grid ω : a = x

0

< ... <

x

n

= b. Consider one of the options for a parallel algorithm

for constructing a cubic spline with an error of O(h

2

). Let n =

p x q, q > 3, p – the number of processors in whose memory

the data are respectively allocated:

M

0

, f

0

, f

1

, ..., f

q

, f

q+1

, f

q−1

, f

q

, f

q+1

, ..., f

2q

, f

2q+1

, ...,

f

n−q−1

, f

n−q

, f

n−q+1

, ..., f

n−1

, f

n

, M

n

This arrangement of the data corresponds to the decompo-

sition of the region with overlapping by two steps of the

grid, a fragment of which is shown in fig.1. The spline

building process consists of three stages: 1. The distribution

Fig. 1. Decomposition of area ω with overlapping.

of the source data in each processor in accordance with the

decomposition of the grid ω. 2. Calculation according to the

grid function with an error O(h

2

) of the moments of the spline

at the boundaries of the intervals of decomposition of the

region:

M

k

= (F

k−1

− 2F

k

+ F

k+1

)/F

2

,

M

k−q

= (F

k+q−1

− 2F

k+q

+ F

k+q+1

)/h

2

,

k = q, 2q, 3q, ..., (p − 1)q.

(19)

3. The solution by sweep method in each processor of the

three–diagonal system to determine the moments. In this case,

it will be necessary to perform the order of O(q) operations on

each processor. The degree of parallelism can be increased by

using one of the parallelization methods of the sweep method.

To calculate the spline value at some point x from [a,b], it is

necessary to determine the grid interval to which this point

belongs and use the formula (12) or (14). The determination

of the grid interval can be carried out in parallel, using the

decomposition of the grid region without overlapping.

If the spline is restored through basis functions, then to the

system (19) at the boundaries of the decomposition intervals,

it is necessary to add the equations































b

k−1

B

00

k−1

(x

k

) + b

k

B

00

k

(x

k

) + b

k+1

B

00

k+1

(x

k

) =

(f

0

k−1

− 2f

k

+ f

k+1

)/h

2

;

b

k+q−1

B

00

k+q−1

(x

k+q

)+

b

k+q

B

00

k+q

(x

k+q

) + b

k+q+1

B

00

k+q+1

(x

k+q

) =

= (f

k+q+1

− 2f

k+q

+ f

k+q+1

)/h

2

;

k = q, 2q, 3q, ..., (p − 1)q.

After bringing this system to a tridiagonal form, the solution is

found by the sweep method. Taking into account the properties

of the supports of B–splines from (18), it is convenient to use

the formula for calculating

S(z), z ∈ (x

i

, x

x+1

), i = 0, n − 1,

S(z) = b

i

B

i−1

(z) + b

i

B

i

(z) + b

i+1

B

i+1

(z) + b

i+2

B

i+2

(z).

The calculation of the spline using this formula can be carried

out on a p–processor system by decomposing the grid region

in accordance with fig.1. The parallel algorithm considered

above for constructing the cubic spline of defect 1 as p− > n

leads to the construction of a local cubic spline on each grid

interval ω. If conditions (20) are not used, then the solution

of spline systems can be found, for example, by the cyclic

reduction method [5], [7], [10], [12]. The problem statement

is such that it is required to calculate cubic spline of degree

m=3 using various parallel programming technologies and

analyze the results of calculations by the main efficiency

parameters by comparison. To implement the parallel algo-

TABLE I

P

ERFORMANCE

P

ARAMETERS OF

O

PEN

MP

AND

MPI

Num.

Number of OpenMP Threads

cal.

1

4

8

T1,mc

T4,mc

S

E

T8,mc

S

E

8X10

4

2,0

1

2

1.00

0.72

2.78

0.69

8X10

5

30,0

14,6

2,5

1.03

4.6

6.52

10.63

8X10

6

220,0

89,2

2,47

1.23

3,3

6.67

1.67

8X10

7

1710,0

526,8

3,25

1.62

269,1

7.15

1.79

8X10

8

16800,0

4845,5

3,47

1.73

2200,0

7.64

1.91

Num.

Number of MPI Processes

cal.

1

4

8

T1,mc

T4,mc

S

E

T8,mc

S

E

8X10

4

2,0

1.2

1,67

0.83

1

2.00

0.50

8X10

5

30,0

15,7

1,91

1.96

5.6

5.36

1.34

8X10

6

220,0

90,8

2,42

1.21

36,2

6.08

1.52

8X10

7

1710,0

546,6

3,13

1.56

279,1

6.12

1.53

8X10

8

16800,0

4995,5

3,36

1.68

2378,0

7.06

1.77

rithm for computing the cubic spline, we chose the integrated

development environment MS Visual Studio 2017 and the

C++ programming language, and the TUIT training cluster,

including 4 physical cores Intel (R) Core (TM) i7, was used

as the technological platform for conducting computational

experiments –7700 CPU 3.60 GHz (8 threads). In the process

of computing, parallelism based on data was applied to the

parallel algorithm, that is, each process calculated a certain

part of the spline function.

C

ONCLUSION

So, the use of parallel computing systems today is relevant.

The development of various fields of science and technology

expands the range of applied tasks more laborious than the

previous ones, for the solution of which it is necessary to use

powerful computers with the ability to perform parallel com-

puting. Determining the main parameters of parallel computing

systems allows you to determine the degree of its effectiveness.

The calculation of the cubic spline using various parallel

programming technologies has shown that MPI technology is

the most effective. It is the most relevant in modern times, but

it requires significant costs.

R

EFERENCES

[1] A.R. Aydinya, “Apparatnyye sredstva vychislitelnoy tekhniki,” Direct-

media, 2016, pp. 125–145.

[2] A.C. Antonov, “ Parallelnoye programmirovaniye s ispolzovaniyem

tekhnologii MPI,” textbook, MGU, 2004. p. 71.

[3] V.A. Billing, “Parallelnyye vychisleniya i mnogopotochnoye program-

mirovaniye,” textbook, NOU TUIT, volume 2, 2016, p. 311.

[4] I.G. Burova, and Yu.K. Demyanovich, “Minimalnyye splayny i ikh

prilozheniya,” textbook, SPbGU, 2010, p. 364.

[5] V.P. Gergel, “ Vysokoproizvoditelnyye vychisleniya dlya mnogoy-

adernykh mnogoprotsessornykh sistem,” textbook, NNGU, 2010, p. 364.

[6] O.U. Mallaev, “Algorithm and software for speeding up computer

memory using OpenMP technology,” Descendants of Mohammed al-

Khwarizmi TUIT, pp. 71–75.

[7] H. Zayniddinov, S. Madhusudan, and S. Dhananjay, “Polynomial Splines

for Digital Signal and Systems,” Monograph, LAMBERT Academic

Press, Germany, 2016, p. 208.

[8] H.N. Zaynidinov, G.O. Tojiboev, and O.U. Mallaev, “Parallel algorithms

for processing seismic signals on multicore processors,” Automation and

software engineering, 2018, volume 1, pp. 89–95.

[9] H.N. Zaynidinov, and O.U. Mallayev, “Paralleling of calculations and

vectorization of processes in digital treatment of seismic signals by cubic

spline,” IOP Conf. Series: Materials Science and Engineering 537, 2019,

doi:10.1088/1757-899X/537/3/032002, pp. 68–69.

[10] S. Dhananjay, S. Madhusudan, and H. Zaynidinov, “Parabolic Splines

based One-Dimensional Polynomial,” Applied Sciences and Technology,

2019, volume 1, pp. 1-10.

[11] H.N. Zaynidinov, and O.U. Mallayev, “Vectorization of parallel process-

ing of seismic signals using technology of Parallel Studio,” Descendants

of Mohammed al-Khwarizmi TUIT, 2018, volume 2, pp. 14-25.

[12] H.N. Zaynidinov, and O.U. Mallayev, “Definition of synchronization

processes during parallel signal processing in multicore processors,”

2019 International Conference on Information Science and Communi-

cations Technologies (ICISCT), 2019, 4–6 Nov, Publisher: IEEE, DOI:

10.1109/ICISCT47635.2019.9012006.