An Algorithm for Dynamic Load Balancing

of Synchronous Monte Carlo Simulations on Multiprocessor Systems

Peter Altevogt^{∗}, Andreas Linke Institute for Supercomputing and Applied Mathematics (ISAM) Heidelberg Scientific Center IBM Deutschland Informationssysteme GmbH Vangerowstr. 18 69115 Heidelberg Tel.: 06221–59–4471 Germany

Abstract

We describe an algorithm for dynamic load balancing of geometrically parallelized synchronous Monte Carlo simulations of physical models. This algorithm is designed for a (heterogeneous) multiprocessor system of the MIMD type with distributed memory. The algorithm is based on a dynamic partitioning of the domain of the algorithm, taking into account the actual processor resources of the various processors of the multiprocessor system.

Keywords: Monte Carlo Simulation; Geometric Parallelization; Synchronous Algorithms; Dynamic Load Balancing; Dynamic Resizing

IBM preprint 75.93.08, hep-lat/9310021

∗altevogt@dhdibm1.bitnet

1 Introduction

During the last years, Monte Carlo simulations of scientific problems have turned out to be of outstanding importance [1, ?, ?, ?, ?]. It is now a general belief within the community of computational scientists that multiprocessor systems of the MIMD^{1} type with distributed memory are the most appropriate computer systems to provide the computational resources needed to solve the most demanding problems e.g. in High Energy Physics, Material Sciences or Meteorology.

Implementing a synchronous parallel algorithm on a heterogeneous MIMD system with distributed memory (e.g. on a cluster of different workstations), a load balancing between the processors of the system (taking into account the actual resources being available on each node) turns out to be crucial, because the processor with the least resources determines the speed of the complete algorithm.

The heterogenity of the MIMD system may not only result because of heterogeneous hardware resources, but also due to a heterogeneous use of homogeneous hardware resources (e.g. on a workstation cluster, there may exist several serial tasks running on some of the workstations of the cluster for some time in addition to the parallel application; this results in a temporary heterogenity of the cluster, even if the workstations of the cluster are identical). This kind of heterogenity can in general only be detected during the runtime of the parallel algorithm.

Therefore, the usual approach of geometric parallelization [6, ?] to parallelize algorithms by a static decomposition of a domain into subdomains and associating each subdomain with a processor of the multiprocessor system is not appropriate for a heterogeneous multiprocessor system. Instead, on a heterogeneous system this geometric parallelization should be done dynamic.

In the sequel, we consider geometrically parallelized Monte Carlo simulations consisting of update algorithms (e.g. Metropolis or heatbath algorithms) defined on e.g. spins at the sites of a lattice (e.g. in Ising models), matrices defined on the links of the lattice (e.g. in Lattice Gauge Theories), etc.. In general, a synchronization of the parallel processes takes place after each sweep through the lattice (each iteration).

For this class of simulations, we will introduce an algorithm for dynamic load balancing. Implementing and testing the algorithm for the two–dimensional Ising model, it will be shown, that this algorithm may drastically improve the performance of the Monte Carlo simulation on a heterogeneous multiprocessor system.

The paper consist of basically three parts. In the first part we will introduce a simple model for analyzing the performance of synchronous Monte Carlo simulations on multiprocessor systems with distributed memory, in the second part we will present our algorithm for the dynamic load balancing and finally we will present our numerical results.

1Multiple Instruction stream, Multiple Data stream. For an excellent overview on the various models of computation like SISD, SIMD, etc., see [8].

2 The Performance Model

We consider a heterogeneous^{2} multiprocessor system consisting of n processors. For a parallelized Monte Carlo simulation we measure^{3} at times t_{MC} the times Δt_{i}^{t_{MC}} the simulation has taken for a fixed number of sweeps on each of the processors^{4}. The parallelization is done by geometric parallelization, associating a sublattice i with a characteristic scale L_{i}^{t_{MC}} (e.g. a characteristic side of the sublattice, its volume, etc..) with each of the processors. These scales are choosen such that

L_{i}^{t_{MC}} = c_{i}^{t_{MC}}L  (1)

holds. Here L denotes the scale of the complete lattice and the c_{i}^{t_{MC}} are real numbers with 0 < c_{i}^{t_{MC}} < 1 and ∑_{i} c_{i}^{t_{MC}} = 1^{5}. Using these parameters, we can calculate quantities P_{i}^{t_{MC}}, characterizing the computing resources of processor i at time t_{MC}:

P_{i}^{t_{MC}} := \frac{c_{i}^{t_{MC}}}{Δt_{i}^{t_{MC}}}  (2)

Assuming, that the resources of the processors vary slowly compared with the time the simulation takes for one sweep

P_{i}^{t_{MC}+1} ∼ P_{i}^{t_{MC}},  (3)

we set

P_{i} = P_{i}^{t_{MC}} = const.  (4)

Now we reinterpret formula (2): For fixed P_{i} we want to calculate a set of {c_{i}^{t_{MC}+1}}, such that the time for the next sweep (excluding the time spent for communication)^{6}

Δt(\{c_{i}\}) := max_{i} Δt_{i}(c_{i}) with Δt_{i}(c_{i}) := \frac{c_{i}}{P_{i}}  (5)

for i = 1, ..., n has a minimal value

2In the sense of the introduction.

3Using e.g. a library routine provided by the operating system.

4Here the times Δt_{i} denote “wall clock times” measured in seconds and t_{MC} denotes the “internal”

time of the simulation, measured in numbers of Monte Carlo sweeps.

5The initial choice of the parameters L_{i}^{t_{MC}=0} respectively of the c_{i}^{t_{MC}=0} is quite arbitrary, one could choose e.g. c_{i}^{t_{MC}=0} = \frac{1}{n} for all i.

6From now on throughout this section we always consider the times Δt, Δt_{i}, etc. and the coefficients {c_{i}} to be taken at t_{MC} +1. For the sake of clarity we therefore drop this index throughout this section.

Δt_{min} := min_{{c_{i}}} Δt({c_{i}}).  (6)

A necessary condition for this solution is obviously, that all Δt_{i} must be equal^{7}. Remembering the normalization condition on the constants c_{i}, we arrive at

c_{i} = \frac{P_{i}}{∑_{i} P_{i}} ,  (7)

with i = 1, ..., n. Using (5) this results in

Δt_{min} = \frac{1}{∑_{i} P_{i}} .  (8)

For a homogenous system (all P_{i} equal) (7) would give

c_{i} = \frac{1}{n}  (9)

For a heterogenous system this choice of {c_{i}} results in (using (5))

Δt(\{c_{i} = \frac{1}{n}\}) = \frac{1}{n} \frac{1}{min_{i} P_{i}} .  (10)

As the homogenity H of the multiprocessor system we define therefore the ratio of Δt_{min} with (10):

H = n\frac{min_{i} P_{i}}{∑_{i} P_{i}}  (11)

with 0 ≤ H ≤ 1. The speedup S that can be obtained by the dynamic resizing of the sublattices is the inverse of H:

S = \frac{1}{H} .  (12)

Rewriting (8) in terms of H, we arrive at

Δt_{min} = \frac{1}{n} \frac{H}{min_{i} P_{i}} ,  (13)

For later comparison with our experimental data, let us look at the special case P := P_{1} = ... = P_{n−1} > P_{n} =: P_{min}. In this case we have for Δt_{min}:

7Let us assume that e.g. Δt_{1} > Δt_{2}. Then we could easily make Δt_{1} smaller by a redefinition of c_{1} and c_{2} with c_{1} + c_{2} = const.

Δt_{min} = Δt^{P} \frac{n − H}{n − 1} ,  (14)

with Δt^{P} = \frac{1}{nP} ^{8}. Without load balancing, the time for the total simulation will be determined by the time spent on processor n:

Δt_{max} = \frac{1}{nP_{n}}
= Δt^{P} \frac{n − H}{(n − 1)H}
= Δt_{min} \frac{1}{H} .  (15)

Figure 1 shows the “optimal” curve (using dynamic load balancing):

\frac{Δt^{P}}{Δt_{min}} = \frac{n − 1}{n − H}  (16)

and the one obtained without any load balancing:

\frac{Δt^{P}}{Δt_{max}} = \frac{n − 1}{n − H} H  (17)

for n = 4. These curves will be compared with our experimental results.

3 The Algorithm

In this section we describe our algorithm to perform the dynamic load balancing, based on the performance model described above.

3.1 The Input:

- A characteristic scale L of the lattice (e.g. a side length of the lattice).

- The number n of processors of the multiprocessor system.

- The total number of iterations n_{iter} to be done by the simulation.

- The number of iterations n_{resize} after which a resizing of the sublattices may take place.

- A control parameter ϵ with 0 < ϵ < 1 to determine if a resizing should be done.

8Using P_{min} = \frac{(n−1)H}{n−H} P.

3.2 The Output:

- A dynamic resizing of the domains associated with each processor of the multiprocessor system, taking into account the actual resources of the processors.

3.3 Formal Steps:

1. Read the input.

2. Introduce characteristic scales L_{i}^{t_{MC}} of the sublattices with i = 1, .., n and t_{MC} = 1, .., n_{iter}/n_{resize}, where i denotes the processors and t_{MC} counts the number of resizings that have been done.

3. Calculate the initial characteristic sizes of the sublattices L_{i}^{t_{MC}=0} for all processors according to L_{i}^{t_{MC}=0} = \frac{L}{n}.

4. Associate each of the sublattices with one of the processors.

5. Do on each processor i = 1, .., n (in parallel)

{

- Set t_{MC} = 0.

- For m = 1, .., n_{iter}:

{

- Perform iteration of the Monte Carlo update algorithm on the sublattice^{9}.

- If (m mod n_{resize}) = 0 then {

- Measure the wall–clock time Δt_{i}^{t_{MC}} spent on processor i for doing the calculations excluding the time spent for communications.

- Calculate

P_{i}^{t_{MC}} := \frac{L_{i}^{t_{MC}}}{Δt_{i}^{t_{MC}}}  (18)

to measure the actual resources of each node of the multiprocessor system.

- Communicate the results to all processors of the multiprocessor system.

- Calculate the new characteristic sizes L_{i}^{t_{MC}+1} of the sublattices with

L_{i}^{t_{MC}+1} = \frac{P_{i}^{t_{MC}}}{∑_{n=1}^{n} P_{n}^{t_{MC}}} L (i = 1, .., n − 1)  (19)

and

L_{n}^{t_{MC}+1} = L − ∑_{n=1}^{n−1} L_{n}^{t_{MC}+1}.  (20)

9Possibly including communication with other processors.

- Resize the sublattices if

|L_{i}^{t_{MC}+1} − L_{i}^{t_{MC}}| > ϵL.  (21)

This step may include the communication of parts of the sublattices between the processors and is certainly the critical part of the algorithm. We will introduce an algorithm for this resizing for a special case below.

- Set t_{MC} = t_{MC} + 1.

}

}
}

In the sequel we present an algorithm for resizing the sublattices for the special case that the splitting of the sublattices takes place only in one dimension. We use the host–node (respectively client–server) parallel programming paradigm (see [7]), associating each sublattice with a server process and leaving the more administration oriented tasks (like reading the global parameters of the simulation, starting the server processes, etc.) to the host process^{10}. Let us assume the size of the lattice in the direction of the splitting to be L and that the host process holds arrays a[i] with (i = 1, …, n + 1) for the “Monte Carlo times” t_{MC} and t_{MC} − 1 containing the first coordinate in that direction of the “slice” of the lattice associated with each processor:

a[1] = 1 ≤ a[2] ≤ ... ≤ a[n + 1] = L + 1.  (22)

(In terms of the constants \{c_{i}\} this would mean c_{i} = \frac{a[i+1]−a[i]}{L}.) Now the host process sends messages containing instructions to the node processes in two passes:

1. For i = 2, …, n

{

if ((d := a_{t_{MC}}[i] − a_{t_{MC}−1}[i]) > 0)

{ send message to server i telling it to send its “first” d slices to processor i − 1

}

}

2. For i = n − 1, n − 2, …, 1

{

if (d := a_{t_{MC}}[i + 1] − a_{t_{MC}−1}[i + 1]) < 0

{

10Of course these tasks could also be done by the server processes, resulting in the hostless paradigm of parallel programming. Therefore, our algorithm could also be implemented in a hostless model and our limitation to the host–node model is not a loss of generality.

send message to server i telling it to send its “last” d slices to processor i + 1

}

}

The node processes wait for messages from either the host process or from neighbouring node processes. If there are not enough slices available on a node process to be sent, the node process waits for a message from a neighbour node process to receive additional slices. The two pass algorithm prevents deadlocks.

If the resources of the processors of the multiprocessor system change very rapidly, a multiple communication of data may be necessary and will drastically reduce the efficiency of this algorithm. But this is consistent with the fact, that our complete approach to dynamic load balancing is anyhow only valid for systems with moderatly varying resources, as was already pointed out at the beginning of section 2, see (3).

4 Results for the Two–Dimensional Ising Model

The above described algorithm has been implemented for the parallelized simulation of the two–dimensional Ising model on a cluster of four IBM RISC System/6000 – 550 workstations [7], using the PVM programming environment [9, ?]. Here we have a two–dimensional lattice which is divided into stripes. The objects defined on the lattice sites are spins (i.e. binary variables) and an iteration defined on these objects consists e.g. of a Metropolis algorithm to generate a new spin configuration on the lattice. Each stripe is associated with one workstation. The characteristic scales of the stripes are their widths and the characteristic scale of the lattice is the sum of all widths.

The cluster being completely homogeneous, the heterogeneous situation has been simulated by starting independent processes on one or several nodes of the cluster. This allows the heterogenity of the multiprocessor system to be introduced in a controlled manner^{11}, i.e. to vary the homogenity H and measure (16) resp. (17) as functions of H. Our results are presented in figure 2 for a 1000 × 1000 and a 2000 × 2000 lattice. One clearly sees the qualitative agreement with the prediction of our performance model, see figure 1^{12}. A different point of view consists of looking at the (mega) updates done by the Metropolis algorithm on each spin per second (“MUps”)^{13}. These are presented for a 1000×1000 and a 2000 × 2000 lattice as a function of H in figures 3 and 4 with the dynamic load balancing being done after a certain number of sweeps. It turns out, that the optimal number of sweeps between the load balancing to be performed depends on the size of the problem.

11During the measurements cited below, the cluster has been dedicated to our application.

12Considering the fact, that we have not included the time spent for communication in our model, a quantitative agreement between the theoretical and measured performance cannot be expected. An inclusion of the communication in our model would be very difficult and highly system dependend, e.g. because system parameters like latency and bandwidth may be be complicated functions of the homogenity H.

13These “MUps” constitute a benchmark for spin models.

Figure 1: Performance predicted by our model for 4 processors with and without load balancing.

Figure 2: Performance measured for 4 processors with and without load balancing. 11 Figure 3: Performance measured in MUps for 4 processors for a 1000 × 1000 lattice with and without load balancing being done after a certain number of sweeps.

Figure 4: Performance measured in MUps for 4 processors for a 2000 × 2000 lattice with and without load balancing being done after a certain number of sweeps.

5 Summary

We have introduced an algorithm for dynamic load balancing for synchronous Monte Carlo simulations on a heterogenous multiprocessor system with distributed memory. Implementing this algorithm for the two–dimensional Ising model, we have shown, that it may result in a speedup of a factor 5 - 6 for the above described class of geometrically parallelized algorithms. In many cases, the implementation of the algorithm is straight forward with only little overhead in calculation and communication. For homogeneous systems, almost no performance is lost because the algorithm detects that no resizing is necessary by applying (21). For systems with slowly changing heterogenity^{14}, the algorithm converges very fast and the requirements of the algorithm concerning the computing environment are minimal: the system only has to provide a routine to measure the wall–clock time; such a routine should be available on all operating systems.

Considering the generality of the algorithm introduced above, it may also be useful applied to problems other than Monte Carlo simulations, e.g. in parallel iterative methods for solving linear or nonlinear equations appearing in engineering problems^{15}.

14compared to the time needed for one iteration (sweep)

15Here the domain consists of a lattice, with a matrix element being associated with each of the nodes of the lattice.

References

[1] M.Creutz (ed.), Quantum Fields on the Computer (World Scientific Publishing Co. Pte. Ltd., Singapore 1992).

[2] K. Binder (ed.), The Monte Carlo Method in Condensed Matter Physics, Topics in Applied Physics, Vol. 71, (Springer–Verlag, Berlin, Heidelberg 1992).

[3] K. Binder (ed.), Monte Carlo Methods in Statistical Physics, Topics in Current Physics, Vol. 7, 2nd edition (Springer–Verlag, Berlin, Heidelberg 1986).

[4] K. Binder, D.W.Heermann, Monte Carlo Simulation in Statistical Physics: An Introduction, Springer Series in Solid–State Sciences, Vol. 80 (Springer–Verlag, Berlin, Heidelberg 1988).

[5] J.P.Mesirov (ed.), Very Large Scale Computation in the 21st Century

[6] D.W. Heermann and A.N. Burkitt, Parallel Algorithms in Computational Science, Springer Series in Information Sciences (Springer–Verlag, Berlin, Heidelberg 1991).

[7] P. Altevogt, A. Linke, Parallelization of the Two–Dimensional Ising Model on a Cluster of IBM RISC System/6000 Workstations, Parallel Comp. 19, Vol.9 (1993).

[8] S.G. Akl, The Design and Analysis of Parallel Algorithms, Prentice–Hall International Editions (Prentice–Hall, Inc. 1989).

[9] A. Beguelin, J.J. Dongarra, G.A. Geist, R.Manchek and V.S. Sunderam, A user’s guide to PVM parallel virtual machine. Technical Report ORNL/TM-118 26, Oak Ridge National Laboratory, July 1991.

[10] V.S. Sunderam, PVM: A Framework for Parallel Distributed Computing, Concurrency: Practice&Experience Vol.2 No.4, Dec. 1990.