A Three-Level Parallelisation Scheme and Application to the Nelder-Mead   Algorithm

Rima Kriauzien\.e; Andrej Bugajev; and Raimondas \v{C}iegis

arXiv:1904.05208·cs.DC·September 24, 2019

A Three-Level Parallelisation Scheme and Application to the Nelder-Mead Algorithm

Rima Kriauzien\.e, Andrej Bugajev, and Raimondas \v{C}iegis

PDF

Open Access

TL;DR

This paper introduces a three-level parallelisation scheme to enhance the efficiency of parallel algorithms, exemplified by modified Nelder-Mead methods for solving PDEs and linear systems, with a focus on load balancing and scalability.

Contribution

It proposes a novel three-level parallelisation framework that improves scalability by incorporating less efficient algorithms at an additional level, addressing efficiency drops in traditional two-level schemes.

Findings

01

Enhanced parallel efficiency with three-level scheme.

02

Effective load balancing heuristic for large processor counts.

03

Successful application to PDEs and linear systems using modified Nelder-Mead.

Abstract

We consider a three-level parallelisation scheme. The second and third levels define a classical two-level parallelisation scheme and some load balancing algorithm is used to distribute tasks among processes. It is well-known that for many applications the efficiency of parallel algorithms of the second and third level starts to drop down after some critical parallelisation degree is reached. This weakness of the two-level template is addressed by introduction of one additional parallelisation level. As an alternative to the basic solver some new or modified algorithms are considered on this level. The idea of the proposed methodology is to increase the parallelisation degree by using less efficient algorithms in comparison with the basic solver. As an example we investigate two modified Nelder-Mead methods. For the selected application, a few partial differential equations are solved…

Tables7

Table 1. Table 1: Benchmarks with different sizes J m × N m subscript 𝐽 𝑚 subscript 𝑁 𝑚 J_{m}\times N_{m} of the discrete problem ( 9 )

Eq.	Sizes	Sizes	Sizes
	Benchmark 1	Benchmark 2	Benchmark 3
1	$8000 \times 40000$	$8000 \times 20000$	$8000 \times 10000$
2	$4000 \times 20000$	$4000 \times 20000$	$2000 \times 20000$
3	$2000 \times 20000$	$4000 \times 10000$	$2000 \times 10000$
4	$2000 \times 10000$	$2000 \times 10000$	$1000 \times 20000$

Table 2. Table 2: The results for Benchmark 1. T p subscript 𝑇 𝑝 T_{p} is the CPU time in seconds required to compute one useful point ( 11 )

p	16	32	64
	$k = 1$
Eq. 1	10	22	50
Eq. 2	3	5	8
Eq. 3	2	3	4
Eq. 4	1	2	2
Total number of $p$	16	32	64
Model $T_{M p}$	11.145	5.784	3.614
$T_{p}$	11.003	5.394	3.608
Speed-up	12.679	25.862	38.664

Table 3. Table 3: The results for Benchmark 2. T p subscript 𝑇 𝑝 T_{p} is the CPU time in seconds required to compute one useful point ( 11 ).

p	16	32	64	96	128	128
	$k = 1$			$k = 2$		$K = 1$
Eq. 1	9	18	37	26	37	56
Eq. 2	4	8	15	12	15	18
Eq. 3	2	4	8	7	8	8
Eq. 4	1	2	4	3	4	4
Total number of $p$	16	32	64	96	128	86
Model time	6.59	3.36	2.01	1.65	1.34	1.8
$T_{p}$	6.69	3.37	1.98	1.62	1.33	1.86
Speed-up	13.6	27.03	46.03	56.24	68.25	49.03

Table 4. Table 4: The results for Benchmark 3. T p subscript 𝑇 𝑝 T_{p} is the CPU time in seconds required to compute one useful point ( 11 ).

p	16	32	64	96	128	128
	$k = 1$			$k = 2$		$K = 1$
Eq. 1	8	16	32	24	32	56
Eq. 2	4	8	16	12	16	31
Eq. 3	2	4	8	6	8	8
Eq. 4	2	4	8	6	8	9
Total number of $p$	16	32	64	96	128	104
Model time	3.33	1.76	1.05	0.87	0.7	0.9
$T_{p}$	3.38	1.76	1.06	0.86	0.7	0.95
Speed-up	14.33	27.55	45.96	56.72	69.08	51.23

Table 5. Table 5: The results for Benchmark 1 with E m i n > 0 subscript 𝐸 𝑚 𝑖 𝑛 0 E_{min}>0 . T p subscript 𝑇 𝑝 T_{p} is the CPU time in seconds required to compute one useful point ( 11 ).

p		128
$E_{m i n}$	$0.75$	$0.8$	$0.6$
	$k = 3$	$k = 3$	$K = 1$
Eq. 1	26	19	42
Eq. 2	7	5	8
Eq. 3	4	3	4
Eq. 4	2	1	2
Total number of $p$	117	84	56
Model $T_{M p}$	2.45	3.17	3.84
$T_{p}$	2.49	3.08	3.76
Speed-up	56.03	45.37	37.11

Table 6. Table 6: The γ k subscript 𝛾 𝑘 \gamma_{k} values for direct Nelder-Mead parallelisation

k	$d = 3$	$d = 6$	$d = 7$
2	0.603	0.604	0.606
3	0.584	0.517	0.502

Table 7. Table 7: The γ k subscript 𝛾 𝑘 \gamma_{k} values for the parallel Nelder-Mead algorithm from (Lee and Wiswall, 2007)

k	$d = 3$	$d = 6$	$d = 7$
2	0.668	0.685	0.714
3	–	0.436	0.454
4	–	0.023	0.104
5	–	0.001	0.002
6	–	–	0.001

Equations55

A = {V_{1} ≺ V_{2} ≺ \dots ≺ V_{K}},

A = {V_{1} ≺ V_{2} ≺ \dots ≺ V_{K}},

V_{k} = {v_{1} (X_{k}), v_{2} (X_{k}), \dots, v_{M} (X_{k})}, k = 1, \dots, K,

V_{k} = {v_{1} (X_{k}), v_{2} (X_{k}), \dots, v_{M} (X_{k})}, k = 1, \dots, K,

t_{m} (p_{2}) < t_{m} (p_{1}), \mbox f or p_{1} < p_{2} \leq P_{m} .

t_{m} (p_{2}) < t_{m} (p_{1}), \mbox f or p_{1} < p_{2} \leq P_{m} .

t_{m} (p) \geq t_{m} (P_{m}), p > P_{m} .

t_{m} (p) \geq t_{m} (P_{m}), p > P_{m} .

ar g (p_{1}, \dots, p_{M}) \in Q min 1 \leq m \leq M max t_{m} (p_{m}),

ar g (p_{1}, \dots, p_{M}) \in Q min 1 \leq m \leq M max t_{m} (p_{m}),

Q = {(p_{1}, \dots, p_{M}) : p_{1} + \dots + p_{M} \leq P} .

Q = {(p_{1}, \dots, p_{M}) : p_{1} + \dots + p_{M} \leq P} .

t_{m} (p_{m}) = i = 1 \sum n t_{l_{i}} (p_{l_{i}}), p_{l_{i}} = min (p_{m}, P_{l_{i}}) .

t_{m} (p_{m}) = i = 1 \sum n t_{l_{i}} (p_{l_{i}}), p_{l_{i}} = min (p_{m}, P_{l_{i}}) .

E_{p} (V_{k}) \geq E_{min}, \mbox f or p \leq P_{m},

E_{p} (V_{k}) \geq E_{min}, \mbox f or p \leq P_{m},

P_{m} = min (P_{m}, P_{m}) .

P_{m} = min (P_{m}, P_{m}) .

A_{j} = {V_{1}^{j} ≺ V_{2}^{j} ≺ \dots ≺ V_{K_{j}}^{j}}, j = 1, \dots, J .

A_{j} = {V_{1}^{j} ≺ V_{2}^{j} ≺ \dots ≺ V_{K_{j}}^{j}}, j = 1, \dots, J .

V_{k}^{j} = {v_{1}^{j} (X_{k}), v_{2}^{j} (X_{k}), \dots, v_{M_{j}}^{j} (X_{k})}, k = 1, \dots, K .

V_{k}^{j} = {v_{1}^{j} (X_{k}), v_{2}^{j} (X_{k}), \dots, v_{M_{j}}^{j} (X_{k})}, k = 1, \dots, K .

T_{P} (A_{j}) = T_{P} (V^{j}) K_{j}

T_{P} (A_{j}) = T_{P} (V^{j}) K_{j}

ar g 1 \leq j \leq J min T_{P} (A_{j}) .

ar g 1 \leq j \leq J min T_{P} (A_{j}) .

⎩ ⎨ ⎧ i \frac{\partial u}{\partial t} + \frac{\partial ^{2} u}{\partial x ^{2}} = 0, x \in (A, B), t \in (0, T], u (x, 0) = u_{0} (x), x \in [A, B], L_{l} u (A) = 0, L_{r} u (B) = 0, t \in (0, T],

⎩ ⎨ ⎧ i \frac{\partial u}{\partial t} + \frac{\partial ^{2} u}{\partial x ^{2}} = 0, x \in (A, B), t \in (0, T], u (x, 0) = u_{0} (x), x \in [A, B], L_{l} u (A) = 0, L_{r} u (B) = 0, t \in (0, T],

\displaystyle\omega_{h}=\big{\{}x_{j}:\;x_{0}=A,\,x_{J}=B,\,x_{k}=x_{k-1}+h,\,k=1,\ldots,J\big{\}},

\displaystyle\omega_{h}=\big{\{}x_{j}:\;x_{0}=A,\,x_{J}=B,\,x_{k}=x_{k-1}+h,\,k=1,\ldots,J\big{\}},

ω_{τ} = {t^{n} : t^{n} = n τ, n = 0, \dots, N, N τ = T} .

\partial_{x} U_{j}^{n} = (U_{j + 1}^{n} - U_{j}^{n}) / h, \partial_{\overset{x}{ˉ}} U_{j}^{n} = (U_{j}^{n} - U_{j - 1}^{n}) / h

\partial_{x} U_{j}^{n} = (U_{j + 1}^{n} - U_{j}^{n}) / h, \partial_{\overset{x}{ˉ}} U_{j}^{n} = (U_{j}^{n} - U_{j - 1}^{n}) / h

\partial_{\overset{ˉ}{t}} U_{j}^{n} = (U_{j}^{n} - U_{j}^{n - 1}) / τ, U_{j}^{n - 0.5} = 0.5 (U_{j}^{n} + U_{j}^{n - 1}) .

\partial_{\overset{ˉ}{t}} U_{j}^{n} = (U_{j}^{n} - U_{j}^{n - 1}) / τ, U_{j}^{n - 0.5} = 0.5 (U_{j}^{n} + U_{j}^{n - 1}) .

i \partial_{\overset{ˉ}{t}} U_{j}^{n} + \partial_{x} \partial_{\overset{x}{ˉ}} U_{j}^{n - 0.5} = 0, x_{j} \in ω_{h}, n > 0.

i \partial_{\overset{ˉ}{t}} U_{j}^{n} + \partial_{x} \partial_{\overset{x}{ˉ}} U_{j}^{n - 0.5} = 0, x_{j} \in ω_{h}, n > 0.

\partial_{n} u + e^{- i \frac{π}{4}} D_{t}^{1/2} u = 0

\partial_{n} u + e^{- i \frac{π}{4}} D_{t}^{1/2} u = 0

\displaystyle\partial_{n}u=-e^{-i\frac{\pi}{4}}\bigg{(}\bigg{(}\sum\limits_{k=0}^{l}a_{k}\bigg{)}u{-}\sum\limits_{k=1}^{l}a_{k}d_{k}\varphi_{k}\bigg{)},\;x=a,b,

\displaystyle\partial_{n}u=-e^{-i\frac{\pi}{4}}\bigg{(}\bigg{(}\sum\limits_{k=0}^{l}a_{k}\bigg{)}u{-}\sum\limits_{k=1}^{l}a_{k}d_{k}\varphi_{k}\bigg{)},\;x=a,b,

\frac{d φ _{k} ( x , t )}{d t} + d_{k} φ_{k} (x, t) = u (x, t), x = A, B, k = 1, \dots, l .

\frac{d φ _{k} ( x , t )}{d t} + d_{k} φ_{k} (x, t) = u (x, t), x = A, B, k = 1, \dots, l .

\begin{split}E_{\infty}^{c}=&\min\limits_{\{a_{k},d_{k}\}}V_{k},\quad V_{k}=\max\limits_{1\leq m\leq{\tilde{M}}}v_{m}(X_{k}),\\ v_{m}=&\max\limits_{j\in[0,J_{m}],\,n\in[0,N_{m}]}\big{|}u(x_{j,m},t^{n}_{m})-U_{j,m}^{n}\big{|},\end{split}

\begin{split}E_{\infty}^{c}=&\min\limits_{\{a_{k},d_{k}\}}V_{k},\quad V_{k}=\max\limits_{1\leq m\leq{\tilde{M}}}v_{m}(X_{k}),\\ v_{m}=&\max\limits_{j\in[0,J_{m}],\,n\in[0,N_{m}]}\big{|}u(x_{j,m},t^{n}_{m})-U_{j,m}^{n}\big{|},\end{split}

T_{W p} = 17 \frac{J}{p} + 8 p + T_{c 1} (p),

T_{W p} = 17 \frac{J}{p} + 8 p + T_{c 1} (p),

T_{O p} = c_{1} \frac{J}{p} + T_{c 2} (p) .

T_{O p} = c_{1} \frac{J}{p} + T_{c 2} (p) .

u (t, x) = \frac{exp ( - iπ /4 )}{4 t - i} exp (\frac{i x ^{2} - 6 x - 36 t}{4 t - i}),

u (t, x) = \frac{exp ( - iπ /4 )}{4 t - i} exp (\frac{i x ^{2} - 6 x - 36 t}{4 t - i}),

u (t, x) = \frac{1}{\root + \of 1 + i t / α} exp (ik (x - x^{(0)} - k t) - \frac{( x - x ^{(0)} - 2 k t ) ^{2}}{4 ( α + i t )}),

u (t, x) = \frac{1}{\root + \of 1 + i t / α} exp (ik (x - x^{(0)} - k t) - \frac{( x - x ^{(0)} - 2 k t ) ^{2}}{4 ( α + i t )}),

f (x_{1}, \dots, x_{d}) = i = 1 \sum d - 1 100 (x_{i + 1} - x_{i}^{2})^{2} + (1 - x_{i})^{2},

f (x_{1}, \dots, x_{d}) = i = 1 \sum d - 1 100 (x_{i + 1} - x_{i}^{2})^{2} + (1 - x_{i})^{2},

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Taxonomy

TopicsAdvanced Numerical Methods in Computational Mathematics · Electromagnetic Scattering and Analysis · Electromagnetic Simulation and Numerical Methods

Full text

ad2]Rima KRIAUZIENĖ ]Andrej BUGAJEV ]Raimondas ČIEGIS \correspondingauthorRima KRIAUZIENĖ

\RunauthorsR. Kriauzienė, A. Bugajev and Raimondas Čiegis

A Three-Level Parallelisation Scheme and

Application to the Nelder-Mead algorithm

[

Department of Mathematical modelling

Vilnius Gediminas Technical University, Lithuania, Sauletekio ave. 11, LT-10223 Vilnius

e-mails: [email protected], [email protected], [email protected]

Institute of Data Science and Digital Technologies

Vilnius university, Akademijos str. 4, LT-08663 Vilnius.

Abstract

We consider a three-level parallelisation scheme. The second and third levels define a classical two-level parallelisation scheme and some load balancing algorithm is used to distribute tasks among processes. It is well-known that for many applications the efficiency of parallel algorithms of the second and third level starts to drop down after some critical parallelisation degree is reached. This weakness of the two-level template is addressed by introduction of one additional parallelisation level. As an alternative to the basic solver some new or modified algorithms are considered on this level. The idea of the proposed methodology is to increase the parallelisation degree by using less efficient algorithms in comparison with the basic solver. As an example we investigate two modified Nelder-Mead methods. For the selected application, a few partial differential equations are solved numerically on the second level, and on the third level the parallel Wang’s algorithm is used to solve systems of linear equations with tridiagonal matrices. A greedy workload balancing heuristic is proposed, which is oriented to the case of a large number of available processors. The complexity estimates of the computational tasks are model-based, i.e. they use empirical computational data.

keywords:

multi-level parallelisation, load balancing and task assignment, parallel optimisation, Nelder-Mead algorithm, Wang’s algorithm, model-based parallelisation, finite difference methods.

1 Introduction

Current trends in supercomputing show that in order to accumulate high computing power, computers with more, but not faster, processors are used. This trend induces changes in the development of parallel algorithms. The important challenge is to develop parallelization techniques which enable exploitation of substantially more computational resources than the standard existing methods.

This paper deals with problems that can be split into a collection of independent subproblems and this splitting step is repeated iteratively. The solutions of subproblems define the solution of an initial problem. Thus, an additional splitting step increases the potential parallelisation degree of a parallel algorithm.

Any multi-level parallelisation can be considered as a way to generate a pool of tasks. After the pool of tasks is obtained, it is not important how many parallelisation levels were used. However, often such final simplification of the template leads to a loss of an important information and as a consequence to degraded efficiency of the parallel algorithm. Especially this is true if different levels of the scheme are characterised by different properties of an algorithm that should be properly addressed.

In this paper, we consider a special case of a three level parallelisation. The template of this approach is given in Fig. 1:

•

At the first level of parallelisation we assume that there exist a few parallel alternatives $A_{j}$ (see Figure 1) to the original modelling algorithm. The first level of parallelisation becomes a part of a new parallel algorithm and the degree of the first level parallelism can be selected dynamically during the computations – a selection of the best algorithm is performed. In this paper as an example we consider two new parallel modifications of the Nelder-Mead method [nelder1965simplex].

•

On the second level, a set of computational tasks $V^{j}=\{v_{1}^{j},v_{2}^{j},\ldots,v_{M_{j}}^{j}\}$ (see Figure 1) with different computational complexities is defined. These tasks are solved in parallel. As an example we investigate the case when computation of one value of the objective function requires to solve numerically $M$ partial differential equations. The computational complexities of tasks are non-equal because different discretisation steps must be used for different equations in order to achieve the same accuracy for each equation.

•

The third level defines parallel algorithms to solve tasks from the second level. As an example we use Wang’s algorithm to parallelise the solution of systems of linear equations with tridiagonal matrices [wang1981parallel].

The second and the third levels define a well-investigated two-level parallelisation template. We note that load balancing techniques for two-level parallelisation are widely used in applications, see, e.g., [Baravykaite], [Huismann].

The scheduling problem can be formulated representing a parallel algorithm by a directed acyclic graph (DAG). The vertices define computational tasks, the edges define connections/order among tasks. Then a set of partially ordered computational tasks is scheduled onto a multiprocessors system to minimise the computational time (or to optimise some other performance criteria). It is well known that the scheduling problem is NP complete. Many interesting heuristics are proposed to solve it, we mention greedy algorithms [silko_LNCS2002], genetic algorithms [sharma2015], [Singh], simulated annealing and tabu search algorithms [Kirkpatrick671], [tabu1], [tabu2]. Such algorithms include a possibility of dynamic scheduling and allow for tasks to arrive continuously and they can consider variable in time computational resources.

A scheduling task can be very challenging due to specificity of a given application problem and the necessity to parallelise it on modern parallel architectures. As an example we mention the particle simulation which is solved by appropriate domain decomposition techniques [Furuichi]. Another example is the dynamic load balancing on heterogeneous clusters for parallel ant colony optimisation [Llanes]. In the recent work [datta2019exploiting] it is concentrated on the problem of high-dimensionality of the data while solving subspace clustering problem.

In this article we focus on the scheduling problem, when all tasks in the set are independent and can be solved in parallel. It is well known that the given optimisation problem can be redefined as a problem to equalise the computational times of all processes. The simplest load balancing algorithm is based on the assumption that the computation time is proportional to sizes of sub-tasks. Then the domain decomposition algorithm is applied to guarantee that the sizes of subtasks scheduled for each group of processors are equal [Baravykaite].

The quasi-optimal distributions of tasks can be obtained using the greedy strategy to distribute the work on demand, i.e. to apply dynamic load balancing techniques such as work-stealing [Imam], self-organising process rescheduling [Righi].

However, the efficiency of two-level approach is limited due to a typical saturation of the speed up of parallel algorithms for increased numbers of processors and fixed sizes of tasks. Exactly this situation has motivated us to introduce an additional level of parallelisation template. In most cases the usage of it leads to a less efficient algorithms than the initial state-of-the-art algorithm. But the additional degree of parallelism on the second level gives a large overall speed-up, if the number of available resources is large.

Recent developments of new architectures of parallel processors make even more challenging the task to build accurate theoretical performance models. The empirical data shows that for some advanced algorithms the efficiency of parallel computations can depend non-monotonically on the size of a task. Thus the model-based load balancing method starts to become the main tool in developing efficient and accurate task scheduling algorithms. In our work we build the model for prediction of computation time empirically by solving the specialised benchmarks for a wide range of problem sizes and numbers of processors. In fact this analysis resemblance the classical experimental strong scalability analysis of a given parallel algorithm. We note, that these measurements are always done for all processes working simultaneously in order to reflect their actual performance during the execution of real applications (see, also [Lastovetsky2, Lastovetsky1]).

Here we mention two interesting papers, where the model-based task scheduling algorithms are considered. In [Lastovetsky1], it is concentrated on multicore co-processors Xeon Phi, where the empirical computation time curves are used to find optimal parameters for a workload distribution. The obtained model predicts non-monotonic dependence of computation speed on the sizes of problems. The authors call their approach ”load imbalancing”, however, it can be considered as an advanced balancing which adapts the scheduling algorithm to the specificity of Xeon Phi processors. Obviously in this case the assumption that computation time is proportional to the task size is not valid. In a similar research [Lastovetsky2], computations were performed on non-uniform memory access (NUMA) parallel platform with various shared on-chip resources such as Last Level Cache. Again the model-based approach enables to take into account the specific properties of the algorithm and processors. The matrix multiplication and Fast Fourier Transform are used as benchmark problems. It is interesting to note that, according to the presented results, the globally optimal solutions may not load-balance the sizes of sub-tasks. The authors pay a special attention to the energy efficiency of calculations. We note, that there are some papers that are specifically dedicated to load balancing of energy efficiency [Perez]. In our work we formulate some restrictions that are connected to energy efficiency as well – we do not use additional available computational resources if the parallelisation efficiency drops below some specific level. The other work [Lastovetsky3] is dedicated to model-based optimisation on hybrid heterogeneous systems composed of CPUs and accelerators. In that research authors investigate the problem of communications costs due to uneven workload distribution between accelerators and CPUs. They propose to generalise the $\tau$ -Lop [tauLop] model for heterogeneous computations.

In this paper we propose a general methodology for parallelisation of algorithms. As an example we use it to solve some applied optimisations problems. is shown The superiority of the three level parallelisation scheme is shown, comparing it with two level paralleisation scheme. On the second level a set of different-size tasks is defined, which is a typical situation for computation of one value of a black box objective function. In most cases these tasks (or groups of tasks) are independent but computationally costly. Thus each task also should be solved in parallel. This fact leads to a necessity of the third level. The second and third levels of the template define a set of tasks solved in parallel and some load balancing algorithm should be used to take into account the different sizes of subtasks. The necessity of the additional first level comes from the assumption of having more computational resources than can be utilised by the two-level parallelisation approach. It is a consequence of the efficiency saturation for parallel algorithms when the size of the problem is fixed and the number of processes is increased. We select a different optimization method (or a modification of the basic solver) which gives additional degrees of parallelisation thus enabling the possibility to use more processors. At the first level of the template the optimal algorithm is selected. This part requires to find a compromise between the increased parallelisation degree and the decreased convergence rate of the modified parallel optimization algorithm.

In this work we are also interested to address some green computing (GC) challenges. In a broader sense GC is the practices and procedures of designing, manufacturing, using of computing resources in an environment friendly way while maintaining overall computing performance and finally disposing in a way that reduces their environmental impact [Saha2018]. The research in green computing is done in many areas [nemalikanti2011achieving]: Energy Consumption; E-Waste Recycling; Data Center Consolidation and Optimization; Virtualization; I. T Products and Eco-labeling. One of approaches for optimisation of energy consumption on the software level is the autotuning software, which is able to optimise its own execution parameters with respect to a specific objective function (usually, it is execution time) [carretero2015energy]. Well known examples of autotuning software are: FFTW [frigo2005design] (fast Fourier transformations); ATLAS [whaley1998automatically], PHiPAC [bilmes1997optimizing] (dense matrix computations); OSKI [vuduc2005oski], SPARSITY [im2001optimizing] (sparse matrix computations).

Usually, the goal for any autotuning software is to achieve the same result with the same resources, however, reducing the computation time – in terms of parallelisation it means to increase the parallelisation efficiency. Another way to decrease the power consumption is to increase the efficiency by avoiding inefficient calculations; this may slightly increase the execution time, however will give a reasonable increase of parallel efficiency, which leads to the energy savings. We propose to control the efficiency of the parallel algorithm on the load balancing stage of the parallelization template. In many cases this strategy reduces the amount of computational resources used in computations. This analysis is done a priori, meaning that the user knows how many cores should be used for solving a specific parallel task even before starting real computations.

This paper makes the following contributions:

We propose to extend the typical two level parallelisation, which is usually accompanied by some load balancing technique, by adding one additional level. Also, we investigate the possibility to limit the number of used processors to sustain the parallelisation efficiency at the selected level. This approach let us to avoid the inefficient calculations, supporting the green computing technology.

As an example two different families of parallel Nelder-Mead methods were investigated: the family of the generalised parallel Nelder-Mead method [lee2007parallel] and the parallel versions of the classical Nelder-Mead method. In order to perform the load balancing on the second and third levels of the proposed template, we use the complexity model of tasks which is based on the computational data (also known as model-based), as it is done in recent state-of-the-art works [Lastovetsky1],[Lastovetsky2],[Lastovetsky3]. We demonstrate a big potential of this new technique. 2. 2.

A parallel version of the Nelder-Mead method is proposed, which does not change the convergence properties of the sequential optimisation algorithm. We note, that there were some attempts to parallelise this optimisation method before [lee2007parallel], [Klein2014]. However, in these papers the convergence properties are changed and these changes are not studied comprehensively enough. Moreover, it is questionable whether these parallel algorithms are applicable in the case of small-dimension problems.

Our parallel algorithm leads to an increasement of the parallelisation degree up to factor three. However, the introduced changes do not affect the convergence of the sequential optimisation algorithm. The experimental comparison of this new parallelization algorithm with the state-of-the-art technique [lee2007parallel] is provided. The obtained experimental results show that in the case of the Rosenbrock function the convergence properties of the parallel algorithm [lee2007parallel] are much worse than of the classical sequential Nelder-Mead algorithm.

The rest of this paper is organised as follows. In Section 2 the workload balancing problem is formulated, the selection of the optimal algorithm is provided and a general strategy for workload distribution is presented along with the efficient workload distribution algorithm. In Section 3 the detailed description of three parallelisation levels are given for the studied case. We consider the approximation of boundary conditions of Schrödinger equation. The modified Nelder-Mead method is used to solve local optimisation problems on the first level, on the second level a set of partial differential equations are solved numerically, and on the third level Wang’s algorithm is used to solve systems of linear equations in parallel. In Section 4 the results of computational experiments are provided and the efficiency of the proposed three-level parallelisation template is analysed. In Section 5 the comparison of different Nelder-Mead parallelisation methods is presented. The final conclusions are done in Section 6.

2 Workload balancing problem

In this section we formulate the workload balancing problem for the two level parallelisation. Also we present a greedy scheduling algorithm to distribute the processes among tasks. Next, we introduce the additional level – the first and second levels of the two level parallelisation technique become the second and the third levels, accordingly and the first level is a new parallelisation level. On the first level the selection of the optimal algorithm is performed.

First, we will present two level parallelisation template. Assume that we solve a given problem by using the basic method $A$ . The solution process consists of $K$ blocks of tasks (a simple DAG)

[TABLE]

and all blocks must be solved sequentially one after another. Each block consists of $M$ tasks

[TABLE]

where $X_{k}$ defines a set of parameters for the $V_{k}$ block. $V_{k}$ defines the first level of two level parallelisation scheme. Each task $v_{m}$ can be solved by parallel algorithm – this is the second level of the scheme.

The complexities of tasks $v_{m}$ are different, however, they are known in advance and do not depend on $k$ . For each task $v_{m}$ the prediction of computation time $t_{m}(p)$ , $p\leq P$ , $m=1,\ldots,M$ is given – it is based on the modelling results, $P$ is the number of processors in a parallel system. We assume that up to $P_{m}$ processes the computation time monotonically decreases:

[TABLE]

For $P_{m}$ the predicted computation time function $t_{m}(p)$ reaches the minimum value:

[TABLE]

Such a model of computation time $t_{m}(p)$ is important for algorithms with limited scalability such as Wang’s algorithm. In Fig. 2 we present speed-ups of this algorithm for different sizes of linear systems. It is important to mention that the provided results include some additional costs for computation of the objective function along with Wang’s algorithm computational costs. These additional calculations slightly increase the overall parallelisation scalability, thus the provided figure represents the optimistic scenario for general Wang’s algorithm and the realistic scenario for actual computations, that were done in this paper.

In our specific case this data was derived from a simple benchmark implementing Wang’s algorithm. This benchmark performs computations using different numbers of processes and different problem complexity parameters $J$ . It is important to note, that nodes were artificially loaded with calculations to imitate the real situation. For example, with the number of processes $p=4$ there were 32 tasks that were solved by 128 processes at the same time. Thus this benchmark must be run once, using all processes available.

From Figure 2 it follows that the computation time monotonically decrease till some critical number of processes and therefore the efficient usage of processes is limited to this number of processes. Even for large size systems, when the number of equations is $J=16000$ , the maximum number of processes $P_{m}$ does not exceed 80. This analysis justifies our motivation to use the multi-level approach in order to solve the given applied problem.

In the two-level parallelisation scheme for each block of tasks $V_{k}$ we select the number of processes such that the overall solution time is minimised:

[TABLE]

where a set of feasible processors distributions $Q$ is defined as

[TABLE]

Remark 1

In the case when we solve only few large size tasks and the remaining tasks are much smaller and the number of processes $P$ is not very big, the optimal scheduling is obtained when a few smaller tasks are combined into one group $\widetilde{v}_{m}$ . Then sub-task $\widetilde{v}_{m}$ consists of tasks $v_{l_{1}},\ldots,v_{l_{n}}$ . The computation time for this combined task is predicted by the model:

[TABLE]

In this work we are interested to solve the scheduling problem, when the number of processes is large, so the aggregation step is not used.

Next, we propose a simple greedy partitioning algorithm, which is described in Algorithm 1. It aims to find a near-optimal distribution of $M$ tasks of different sizes between homogeneous $P$ processes by using the model-based complexity model $t_{m}(p)$ (similar ideas are also used in [Lastovetsky2]). We assume that $P\geq M$ . The interesting feature of the presented algorithm is that for a given number of processes $P$ the number of active processes can be taken less than $P$ to minimise the overall execution time of the parallel algorithm.

The algorithm starts from the initial distribution when one process is assigned for each task and the predictions of parallel execution times are calculated using the selected performance model. Then, the greedy iterative procedure is applied to distribute the remaining processes. At each iteration, one additional process is assigned to the task which has the largest predicted computation time. Then its parallel execution time is updated. Iterations are repeated until all processes are distributed or the number of processes for some task reaches the limit $P_{m}$ .

Note, that before $t_{m}(p)$ has reached the minimum, value starts to decrease slowly, thus the parallelisation efficiency drops. Therefore, it may be wise to restrict the number of processes by taking into account the efficieny value.

We define the maximum number of processes $\widetilde{P}_{k}$ for which the efficiency condition is still satisfied

[TABLE]

where $E_{p}(V_{k})=S_{p}(V_{k})/p$ ir $S_{p}(V_{k})=t_{k}(1)/t_{k}(p)$ , $E_{min}\in[0,1]$ is a given efficiency lower bound. Estimate (4) is used to modify the limit of the maximum number of processes (3) that can be used to solve the $j$ -th task

[TABLE]

Therefore, in the presented technique $P_{m}$ includes two restrictions:

•

The number of processes cannot exceed the number after which the speed-up drops down (see Fig. 2).

•

The number of processes is limited by efficiency requirement (4), which states: the number of processes per block of tasks $V_{k}$ is not allowed to be increased if the efficiency of the parallel algorithm on the third level reaches the critical value $E_{min}$ .

In fact the second level of the two-level scheme can be used alone, however, it is limited due to Amdahl’s law [amdahl], i.e. the efficiency begin to drop as the number of processes increases for a fixed size of problem. Two-level approach let us to solve this issue up to some point.

Exactly this situation has motivated us to introduce an additional level of parallelisation template.

In the new three-level parallelisation scheme, the second and third levels represent the two-level scheme part described before. Additionally, we add new first level of the template. We assume that there exist parallel alternative algorithms $A_{j}$ :

[TABLE]

Each block $V_{k}^{j}$ consists of $M_{j}$ independent tasks

[TABLE]

The numbers of blocks of tasks $K_{j}$ , the numbers of tasks per block $M_{j}$ , the sizes of tasks $|v_{m}^{j}|$ may be different for different $j$ .

Next, we select the optimal algorithm according to the number of resources available. We denote

[TABLE]

the total solution time for algorithm $A_{j}$ . The block of tasks $V^{j}$ is solved by using the heuristic proposed above. Then the optimal algorithm is defined as

[TABLE]

The usage of $j>1$ may lead to a less efficient algorithm than the initial basic algorithm. But the additional degree of parallelism gives a large overall speed-up.

3 Application of the three-level parallelisation scheme

First, we briefly present the problem which is used to test our methodology. We solve an initial-boundary value Schrödinger problem formulated in a finite space domain [ShrodABC]:

[TABLE]

where operators $L_{l},L_{r}$ define the nonlocal/transparent boundary conditions.

Let $\omega_{h}$ and $\omega_{\tau}$ be discrete uniform grids with space and time steps $h$ , $\tau$ :

[TABLE]

Let $U^{n}_{j}$ be a numerical approximation of the exact solution $u_{j}^{n}=u(x_{j},t^{n})$ at the grid points $(x_{j},t^{n})\in\omega_{h}\times\omega_{\tau}$ . For functions defined on the grid we introduce the forward and backward difference quotients with respect to $x$

[TABLE]

and similarly the backward difference quotient and the averaging operator with respect to $t$

[TABLE]

We approximate the differential equation (6) by the Crank-Nicolson finite difference scheme [radziunas]

[TABLE]

A very interesting approach to construct the approximate local artificial boundary conditions is based on approximation of the transparent boundary condition

[TABLE]

by rational functions. The discrete boundary conditions can be written as:

[TABLE]

where $\partial_{n}u$ is the normal derivative, $\varphi_{k}$ are solutions of the initial value problem for ODEs [ShrodABC]:

[TABLE]

Our aim is to find optimal values of parameters $\{a_{0},a_{1},\dots a_{l},d_{1},a_{2},\dots d_{l}\}$ , when the following minimisation problem is solved

[TABLE]

and $\widetilde{M}$ specially selected benchmark PDEs are solved.

In all examples we use $l=3$ , i.e., the dimensionality of the optimization problem (11) is equal to 7. Here discrete approximations of PDEs represent the tasks $v_{m}$ in (1). To solve $v_{m}$ we must find solutions of $N$ systems of linear equation with tridiagonal matrix [ShrodABC]. According to our three-level parallelisation scheme, the calculations of a single point in minimisation problem (11) define the block of tasks $V_{k}$ .

The systems of linear equations with tridiagonal matrices are solved using Wang’s algorithm. It is well known that if the size of a system is $J$ and $p$ processes are used then the computation time can be estimated as

[TABLE]

where $T_{c1}(p)$ defines communication costs. The time to compute a value of the objective function $f$ for the specified equation can be estimated as

[TABLE]

In this work instead of theoretical complexity models (12) and (13) we use $t_{m}(p)$ , $m=1,\ldots,M$ , based on empirical computations for a selected set of benchmark problems. Such an approach takes into account all specific details of the parallel algorithm and the computer system.

It is interesting to note that the complexity of computational task $v_{m}$ depends on both parameters: the number of linear equations $J_{m}$ of the system and the number of integration in time steps $N_{m}$ . The computation time $T_{mp}$ is equal to $N_{m}t_{m}(p)$ , but the scalability of the parallel algorithm depends on $J_{m}$ only, since the integration in time is done sequentially step by step.

Next, we present an example with $M=4$ , where four different benchmark PDE problems (6) with explicit solutions [szeftel, Zlotnik] are defined as:

[TABLE]

$x\in[-5,5]$ , $t\in[0,0.8]$ . The problem is approximated on the uniform grid $J\times N=8000\times 4000$ . 2. 2.

[TABLE]

where $k=100,\alpha=1/120,x^{(0)}=0.8$ . $x\in[0,1.5]$ , $t\in[0,0.04]$ . We use the uniform discretisation grid $J\times N=12000\times 4000$ . 3. 3.

The solution is defined by (14), $x\in[-10,10]$ , $t\in[0,2]$ . We use the uniform discretisation grid $J\times N=16000\times 10000$ . 4. 4.

The solution if defined by (15), where $k=100,\alpha=1/120,x^{(0)}=0.8$ . $x\in[0,2]$ , $t\in[0,0.08]$ . We use the uniform discretisation grid $J\times N=16000\times 8000$ .

Next, we consider the problem (11) as a local optimisation problem, which can be solved using an iterative algorithm with a given initial starting point. As a local optimiser Nelder-Mead algorithm is used [nelder1965simplex].

We propose a family of modifications of the original Nelder-Mead algorithm in order to increase the parallelisation degree of it.

At each iteration the following four different scenarios can be obtained:

•

Reflection – compute the value $f_{R}$ of the objective function at the point $X_{R}$ . Depending on the value $f_{R}$ this can be the end of the iteration.

•

Expansion – depending on the $f_{R}$ , an additional computation of the objective function at the point $X_{E}$ is done, meaning the total computation of two objective function values: $f_{R},f_{E}$ .

•

Contraction – depending on the $f_{R}$ , an additional computation of the objective function at point $X_{C}$ is done, meaning the total computation of two objective function values: $f_{R},f_{C}$ .

•

Compression – compute $m$ objective function values, as well as $f_{R}$ and $f_{C}$ . Here $m$ is the number of simplex dimensions.

The first three scenarios require to compute one or two values of the objective function from the set: $f_{R}$ , $f_{E}$ , $f_{C}$ . We can neglect the last scenario, because it occurs very rarely. For the first three scenarios we propose to compute two or three points simultaneously. Algorithmically this means that we change the order of computations, which let us to parallelise the Nelder-Mead method. In most cases only two of three points will be used. Therefore, some redundant calculations will be performed, however, this modification gives an additional parallelisation of computations.

Thus, two modifications of the sequential ( $A_{1}$ ) Nelder-Mead method are defined. For $A_{2}$ we compute in parallel values $f_{R},f_{E}$ and for $A_{3}=3$ we compute in parallel all three values $f_{R},f_{E},f_{C}$ . As a test case we assume that the first scenario is relatively rare, the extension step is done with probability 2/3 and contraction steps occurs with probability 1/3. Then we get that the algorithmic efficiency of the proposed parallel modifications are equal to $\gamma_{2}=0.75$ and $\gamma_{3}=2/3$ , respectively. We note, that these values can be estimated more precisely for specific applications, and one example is given for the computational experiments with the Rosenbrock objective function in Section 5.

On the first level different parallel algorithms can be used, however, the proposed approach is oriented to the cases when the increased degree of parallelisation gives the speed-up at the cost of efficiency which is a typical situation in parallel algorithms theory (Amdahl’s law). As one more example we mention new algorithms developed to solve the global optimisation problems. The modification of the well-known DIRECT method [finkel2003direct] was presented in [Stripinis2018], it is called DIRECT-GL. The new modification is based on the idea at each iteration to analyse more potential optimal rectangles. This approach increases the global sensitivity of the method but in many cases this property is achieved at the cost of additional computations. The potential parallelisation degree of DIRECT-GL algorithm can increase up to 2-3 times. But the results of computational experiments in [Stripinis2018] show that for many benchmark problems (in [Stripinis2018] these cases are numbered 1,2,5,6,20,21,22,24,35,37,38,47,48,49,52) the DIRECT-GL algorithm increases the computational costs to achieve the same accuracy of approximations as DIRECT algorithm. Thus, the classical DIRECT algorithm and its modification DIRECT-GL fit well into the proposed three-level parallelisation template. Then the degree of parallelisation should be increased only if this increasement compensates the reduced efficiency of the modified algorithm. Thus we state, that in order to apply the proposed three level parallelisation scheme, first the computations of one point should be parallelised by a two-level parallelisation approach. Then alternative cases of parallel algorithms with additional degrees of parallelisation should be identified and the optimal algorithm should be selected.

4 Experimental results

In this section we present results of the parallel scalability tests. All parallel numerical tests in this work were performed on the computer cluster “HPC Sauletekis” at the High Performance Computing Center of Vilnius University, Faculty of Physics. We have used up to 8 nodes with Intel® Xeon® processors E5-2670 with 16 cores (2.60 GHz) and 128 GB of RAM per node. Computational nodes are interconnected via the InfiniBand network.

Our main goal is to investigate the efficiency of the proposed three level template of workload distribution between processes. First, we have selected three specific benchmarks with different discretizations (7), when $M=4$ discrete approximations of PDEs (9) are solved numerically to compute one value of the objective function. The sizes $(J_{m}\times N_{m})$ , $m=1,\ldots,4$ of discrete problems are given in Table 1.

In the first benchmark the size of one task $v_{1}$ is much bigger than sizes of the remaining three tasks. In the second benchmark two changes are done. They make this set of tasks more suited for parallelisation on large number of processes: the size of task $v_{1}$ is reduced twice by taking a smaller number of time steps $N_{1}$ ; the size of task $v_{3}$ remains the same, but the number of points $J_{3}$ is increased twice, therefore the scalability of Wang’s algorithm is improved for this task. In the third benchmark the relative sizes of tasks $v_{m}$ are more homogeneous than in the first benchmark, but this result is achieved by reducing the number of space grid points $J_{2},J_{4}$ , therefore the scalability of Wang’s algorithm is decreased for these two tasks, especially for $v_{4}$ .

First, we exclude the efficiency condition from the load balancing algorithm by taking $E_{min}=0$ in (4). The distribution of processors between tasks are presented in Tables 2–4. We also provide the actual computation time $T_{p}$ along with $T_{Mp}$ that were predicted by the theoretical complexity model. As we can see from Table 2 the model and experimental times are close to each other. The experimental time is smaller in cases when there is no interpolation error. Also it is smaller than the model time – it is expected result, model times (see Figure 2) are based on benchmark, that imitate pessimistic scenario – as it was mentioned before, all nodes were artificially loaded at the same time. The prediction accuracy depends on many parameters such as cluster architecture, network loads during computations.

For comparison purposes we provide the results obtained by using the two-level parallelisation template. $K=1$ , then the first level of the three-level template is not used.

It is important to note, that in Tables 2-5 we present the CPU time needed to compute one useful point (11), i.e., the actual time is divided by $\gamma_{k}\,k$ , which represents the usefulness of computations. Optimal algorithm $A_{k}$ is selected automatically using the approach that was described above.

As it follows from Table 2, the usage of the first level with $k=3$ and $P=128$ processes increases the potential speed-up from 38.75 to 60.44. If $P=128$ and $k=1$ then only 70 processes are used. However the result is very similar to the case when $P=64$ processes are used, which means that these additional resources are used very inefficiently.

In the Fig. 3 the Gantt charts show theoretical model time $t_{m}(p)$ , that is needed to obtain the solutions of different equations. The workload distribution becomes closer to uniform as the number of processes is increased.

4.1 The control of efficiency

The reduction of the energy consumption is an important goal, especially when increasment of computation speed-up are small for additional processes. The presented results indicate that in some cases there is a highly inefficient usage of computational resources.

For the purposes of controlling the efficiency of calculations the condition (4) was introduced in Algorithm 1. This condition guarantees that the efficiency of the numerical solution of each block of tasks will be at least $E_{min}$ . It is important to note, that we are not attempting to generate optimal mappings of processors – we have developed an heuristic that provides the quality of distribution of tasks, that is sufficient for the most practical purposes. The quality of the algorithm is improved when more processors are available.

Next, a more detailed analysis of the Benchmark 1 is provided. In Table 5 the results for $E_{min}>0$ are presented. Comparing the results in Table 5 with the results in Table 2 we see that for $K=1$ and $E_{min}=0.6$ the number of processes for the first equation is decreased by $14$ , however, the computation times are almost the same as it was in the case of $E_{min}=0$ . Also, for $K=3$ the efficiency requirement begins to limit the number of processes for $E_{min}=0.75$ and it decreases further with $E_{min}=0.8$ .

However, even then a three level approach with $K=3$ is superior to the standard two-level approach in terms of the final speed-up. The results in Table 5 indicate that even for the efficiency limitation $E_{min}=0.75$ the proposed three-level approach lets to maintain a big number of parallel processes active, this number is equal to $(26+7+4+2)\times 3=117$ . The speed-up is 56 and the efficiency of the parallel algorithm is $56/117\approx 0.48$ . The last column in Table 5 with $K=1$ presents the results for the two-level approach (without the first level). A straightforward two-level parallelisation approach would have the limited parallelisation possibility especially for problems of the size $J=2000$ . For such small subproblems it would be possible to utilise only up to 32 processes (Fig. 2), the speed-up would be quite limited as well.

Note, that all previous results represents the analysis based on a single Nelder-Mead iteration. Next, we solve the actual real-world optimisation problem (11). The maximum number of processes $P=128$ the load balancing algorithm has selected $k=1$ . The number of Nelder-Mead method iterations was fixed to 1000. The parallel and sequentional versions gave the same results the minimum value of the error $E_{\infty}^{C}=0.0806$ . The sequentional version of computations took $180286$ seconds, the parallel version computations took $2232$ seconds. Thus, a speed-up factor of $81.8$ was achieved. The selection of $k=1$ indicates that the number of processes can be greatly increased – the algorithm has selected $k=1$ automatically for a given number of processes.

5 The comparison of different Nelder-Mead parallelisation methods

Here we present the analysis of the convergence properties of different modifications of the Nelder-Mead method. As it was mentioned before, the convergence rate of the selected algorithm directly affects the parallelisation efficiency, which is represented by $\gamma_{k}$ , where $k$ is the parallelisation degree. In this section we measure $\gamma_{k}$ by measuring the experimental parallel efficiency of algorithms.

The detailed analysis of convergence behaviour for different objective functions is out of the scope of this research. However, the objective function from the previous sections is suitable for a narrow class of applications. Thus, to perform a comparison of different parallel versions of Nelder-Mead method we minimise the Rosenbrock objective function that is widely used by researchers in the field of optimisation theory [fajfar2018nelder], [Stripinis2018].

We show that in the case of the Rosenbrock function the real experimental $\gamma_{k}$ values are different than were assumed to be in the experiments of the previous sections. The reason is that the significant number of iterations require to compute only one point $F_{R}$ .

We compare the results of our parallel modification of the Nelder-Mead method with the state-of-art technique proposed in [lee2007parallel]. As a benchmark we use the Rosenbrock function

[TABLE]

which makes the optimisation problem challenging. It should be noted that the parallel algorithm [lee2007parallel] can achieve the parallelisation degree $K$ that is equal to the optimisation problem dimension $d$ . Thus potentially this algorithm is well suited for parallel computers with a big number of processes.

In the Table 6 we compare three cases $d=3,6,7$ : $d=3$ – the minimum, that is needed for parallelilsation with both methods, $d=7$ – the case that was studied in previous section, $d=6$ – to show the tendency for smaller $d$ . We provide results obtained when the Rosenbrock function of different dimensions $d=3,6,7$ was minimized by using our parallel modification of the Nelder-Mead method. The values of the efficiency coefficients $\gamma_{k}$ are presented. They show that this parallel algorithm is quite stable and it is well-suited to be used in the three-level template solver for small dimension objective functions.

Table 7 presents results obtained by using the state-of-the-art parallel Nelder-Mead algorithm from [lee2007parallel]. It follows, that in all investigated cases the parallelisation degree is very limited, since the convergence drops significantly when the parallelization degree is increased. This method is mainly targeted to solve problems when the dimension of the objective function is big (e.g. problems in financial mathematics, when $d\approx 100$ ).

6 Conclusions

In this paper we introduced a three-level parallelisation template which utilises a new model-based load balancing which is based on experimental data. The technique was tested for three benchmarks. The experimental results confirmed the good accuracy of the new time prediction model.

Comparing the three-level template to the classical two-level scheme, the proposed scheme looks more promising for development of efficient parallel algorithms in the case when a big number is computational resources is available.

The possibilities of the three-level parallelization template are demonstrated for solving local optimization problems. On the first level a well-known Nelder-Mead algorithm was used. We proposed a family of parallel versions of this method, which increases the parallelisation degree up to the factor three. The proposed load balancing algorithm chooses the optimal version of the parallel Nelder-Mead algorithm. It dynamically increases the parallelisation degree on the first level when the speed-up of the second and third levels begins to saturate.

For the considered test problem on the second level $M$ PDEs were solved numerically and on the third parallelisation level Wang’s algorithm was used to solve systems of linear equations. It was shown that there exists a limit for the speed-up that can be achieved due to limitations of Wang’s algorithm. The proposed approach extends the parallelisation degree allowing to achieve an additional speed-up.

The proposed load balancing algorithm limits the size of computational resources to preserve the efficiency requirement which can be controlled by selecting the parameter $E_{min}$ .

{acknowledgment}

Computations were performed using resources at the High Performance Computing Centre ”HPC Sauletekis” in Vilnius University Faculty of Physics.

\harvardpreambledefs\harvardpreambletext

Bibliography39

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1[1] \harvarditem Amdahl 1967 amdahl Amdahl, G. \harvardyearleft 1967 \harvardyearright . Validity of the single processor approach to achieving large-scale computing capabilities, AFIPS Conference Proceedings 30 : 483–485.
2[2] \harvarditem [Bilmes et al .]Bilmes, Asanovic, Chin \harvardand Demmel 1997 bilmes 1997 optimizing Bilmes, J., Asanovic, K., Chin, C.-W. \harvardand Demmel, J. \harvardyearleft 1997 \harvardyearright . Optimizing matrix multiply using phipac: a portable, high-performance, ansi c coding methodology, Proceedings of the 11th international conference on Supercomputing , ACM, pp. 340–347.
3[3] \harvarditem [Bugajev et al .]Bugajev, Ciegis, Kriauziene, Leonaviciene \harvardand Zilinskas 2017 Shrod ABC Bugajev, A., Ciegis, R., Kriauziene, R., Leonaviciene, T. \harvardand Zilinskas, J. \harvardyearleft 2017 \harvardyearright . On the accuracy of some absorbing boundary conditions for the schrodinger equation, Mathematical Modelling and Analysis 22 (3): 408–423.
4[4] \harvarditem [Carretero et al .]Carretero, Distefano, Petcu, Pop, Rauber, Rünger \harvardand Singh 2015 carretero 2015 energy Carretero, J., Distefano, S., Petcu, D., Pop, D., Rauber, T., Rünger, G. \harvardand Singh, D. E. \harvardyearleft 2015 \harvardyearright . Energy-efficient algorithms for ultrascale systems, Supercomputing frontiers and innovations 2 (2): 77–104.
5[5] \harvarditem Ciegis \harvardand Baravykaite 2007 Baravykaite Ciegis, R. \harvardand Baravykaite, M. \harvardyearleft 2007 \harvardyearright . Implementation of a black-box global optimization algorithm with a parallel branch and bound template, Applied Parallel Computing: State of the Art in Scientific Computing 4699 : 1115–1125.
6[6] \harvarditem Čiegis \harvardand Šilko 2002 silko_LNCS 2002 Čiegis, R. \harvardand Šilko, G. \harvardyearleft 2002 \harvardyearright . A scheme for partitioning regular graphs, in R. Wyrzykowski, E. Deelman, J. Dongarra, K. Karczewski, J. Kitowski \harvardand K. Wiatr (Eds), Proc. 4th International Conference on Parallel Processing and Applied Mathematics (PPAM 2001, Naleczsow, Poland, September 9-12, 2001) , Vol. 2328 of Lecture Notes in Computer Science , Springer, Berlin, Germa
7[7] \harvarditem [Datta et al .]Datta, Kaur, Lauer \harvardand Chabbouh 2019 datta 2019 exploiting Datta, A., Kaur, A., Lauer, T. \harvardand Chabbouh, S. \harvardyearleft 2019 \harvardyearright . Exploiting multi–core and many–core parallelism for subspace clustering, International Journal of Applied Mathematics and Computer Science 29 (1): 81–91.
8[8] \harvarditem [Fajfar et al .]Fajfar, Bűrmen \harvardand Puhan 2018 fajfar 2018 nelder Fajfar, I., Bűrmen, Á. \harvardand Puhan, J. \harvardyearleft 2018 \harvardyearright . The nelder–mead simplex algorithm with perturbed centroid for high-dimensional function optimization, Optimization Letters pp. 1–15.