Billiards and related systems on the bulk-synchronous parallel model.
With two examples we show the suitability of the bulk-synchronous parallel (BSP) model for discrete-event simulation of homogeneous large-scale systems. This model provides a unifying approach for general purpose parallel computing which in addition to efficient and scalable computation, ensures portability across different parallel architectures. A valuable feature of this approach is a simple cost model that enables precise performance prediction of BSP algorithms. We show both theoretically and empirically that systems with uniform event occurrence among their components, such as colliding hard-spheres and ising-spin models, can be efficiently simulated in practice on current parallel computers supporting the BSP model.
Introduction
Parallel discrete-event simulation of billiards and
related systems is considered a non-obvious algorithmic
problem, and has deserved attention in the literature
[1, 5, 7, 8, 9, 11, 13, 18, 24, 23, 25]. Currently an
important class of applications for these simulations
is in computational physics [6, 7, 10, 14, 15, 20, 21]
(e.g. hard-particle fluids, ising-spin models, disk-
packing problems). However this kind of systems
viewed through the more general setting of "many
moving objects" [3, 16], are present everywhere in
real life (e.g. big cities, transport problems, navigation
systems, computer games, and combat models!).
On the other hand, these systems have been considered
sufficiently general and computationally intensive
enough to be used as a sort of benchmark for Time
Warp simulation [5, 23, 25], whereas different simulation
techniques have been shown to be more efficient
when dealing with large systems [8, 9, 11, 13, 18].
Similar to most of the parallel software development
in the last few decades, the prevalent approach to the
simulation of these systems has followed a machine
dependent exploitation of the inherent parallelism associated
with the problem. Currently, however, one of
the greatest challenges in parallel computing is: "to
establish a solid foundation to guide the rapid process
of convergence observed in the field of parallel computer
systems and to enable architecture independent
software to be developed for the emerging range of
scalable parallel systems" [17]. The bulk synchronous
parallel (BSP) model has been proposed to provide
such a foundation [22] and, for a wide range of ap-
plications, this model has already been shown to be
successful in this bridging role (i.e. a bridge between
hardware and software in direct analogy with the role
played by the von Neumman model in sequential computing
over the last fifty years). At present, the BSP
model has been implemented in different parallel architectures
shared memory multi-processors, distributed
memory systems, and networks of workstations
enabling portable, efficient and scalable parallel
software to be developed for those machines [19, 4].
A first step in the BSP implementation of conservative
and optimistic parallel simulation algorithms
has so far been given in [12]. In this paper we follow
a different approach by using conservative algorithms
designed on purely BSP concepts, and evaluating their
performance under two examples: an ising-spin model
and a hard-particle fluid. Note that these potentially
large-scale systems have the property of a very random
and even distribution of events among their constituent
elements. We believe, however, that these two
examples exhibit sufficient generality and complexity
as to be representative of a wide range of other related
asynchronous systems (e.g. some instances of
the multiple-loop networks described in [8] and the
systems there mentioned). Note that because of the
synchronous nature of the BSP model, our algorithms
are reminiscent to those proposed in [8, 18].
2 The BSP model
For a detailed description of the BSP model the
reader is referred to [22, 17]. A bulk-synchronous
parallel (BSP) computer consists of: (i) a set of
processor-memory pairs, (ii) a communication net-work
that delivers messages in a point-to-point man-
ner, and (iii) a mechanism for the efficient barrier
synchronization of all, or a subset, of the processors.
There are no specialized broadcasting or combining
facilities.
If we define a time step to be the time required for
a single local operation, i.e. a basic operation such as
addition or multiplication on locally held data values,
then the performance of any BSP computer can be
characterized by the following four parameters: (i) p
the number of processors, (ii) s the processor speed,
i.e. number of time step per second, (iii) l the synchronization
periodicity, i.e. minimal number of time steps
Network l g
Ring O(p) O(p)
2D Array O( p p) O( p p)
Butterfly O(logp) O(logp)
Hypercube O(logp) O(1)

Table

1: BSP parameters for some parallel computers.
elapsed between two successive barrier synchronizations
of the processors, (iv) g the ratio total number of
local operations performed by all processors in one second
to total number of words delivered by the communication
network in one second, i.e. g is a normalized
measure of the time steps required to send/receive an
one-word message in a situation of continuous traffic
in the communication network. See Table 1 taken
from [17] which shows bounds for the values of g and
l for different communication networks.
A BSP computer operates in the following way. A
computation consists of a sequence of parallel super-
steps, where each superstep is a sequence of steps, followed
by a barrier synchronization of processors at
which point any remote memory accesses takes effect.
During a superstep each processor has to carry out a
set of programs or threads, and it can do the following:
(i) perform a number of computation steps, from its
set of threads, on values held locally at the start of the
send and receive a number of messages
corresponding to non-local read and write requests.
The complexity of a superstep S in a BSP algorithm
is determined as follows. Let the work w be the
maximumnumber of local computation steps executed
by any processor during S. Let h s be the maximum
number of messages sent by any processor during S,
and h r be the maximum number of messages received
by any processor during S. Then the cost of S is given
by time steps (or alternatively
g). The cost of a
BSP algorithm is simply the sum of the costs of its
supersteps.
The architecture independence is achieved in the
BSP model by designing algorithm which are parameterized
not only by n, the size of the problem, and p,
the number of processors, but also by l and g. The
resulting algorithms can then efficiently implemented
on a range of BSP architectures with widely differing l
and g values. For example, on a machine with large g
we must provide an algorithm with sufficient parallel
slackness (i.e. a v processor algorithm implemented on
a p processor machine with v ? p) to ensure that for
every non-local memory access at least g operations
on local data are performed.
3 Basic BSP simulation algorithms
The kind of systems relevant to this paper (i.e.
statistically homogeneous steady state systems with
event occurrences randomly and evenly distributed
among their constituent elements) can be simulated
on a BSP computer using two-phase conservative algorithms
as follows.
On a p-processor BSP computer the whole system
is divided into p equal-sized regions that are owned
by a unique processor. Events involving elements located
on the boundaries are called border zone events
(BZ events), and are used to synchronize the parallel
operation of the processors. The most conservative
(but less efficient) version of this algorithm works doing
iterations composed of two phases: (i) the parallel
phase where each processor is simultaneously allowed
to simulate sequentially and asynchronously its own
region, and (ii) the synchronization phase where the
occurrence of one border zone event is simulated by
only one processor while the other ones remain in an
idle state. [We further improve the efficiency of this
algorithm by exploiting opportunities to simulate at
most p border zone events in parallel during the synchronization
phase.]
The synchronization phase is used to cause the barrier
synchronization of the processors in the simulated
time, and to exchange state information among
neighboring regions. [In the system examples studied
below, this state information refer to the states of
particular atoms and particles located in neighboring
regions.] During the parallel phase every processor
simulates events whose times are less than the current
global next BZ event (i.e. the BZ event with the
least time among all of the local next BZ events held
in each region or processor). Thus, global processor
synchronization is issued periodically at variable time
intervals which are driven by the chronological occurrence
of the BZ events. See pseudo-code in Figure 1.
We assume there are n elements evenly distributed
throughout the whole system, with regions made up
of
n=p) \Theta (a =
n=p) elements. The goal
is to simulate the occurrence of M events which on
average are assumed to occur randomly and evenly
distributed among the elements (i.e. M=n per ele-
ment). This goal is achieved by the BSP algorithm in
I iterations, wherein each iteration simulates a total
of NPE events in the parallel phase plus one event in
the synchronization phase, namely
We define f I =M to be the fraction of BZ events
that occur during the simulation. As we show below,
for the kind of systems we are interested in we have
leading to
which shows that by choosing regions sufficiently large
it is always possible to achieve some degree of parallelism
with this strategy. However, the actual gain
in running time due to the parallel phase, where each
processor simulates about O(a=p) events sequentially,
crucially depends on the cost of the communication
and synchronization among the processors during the
synchronization phase (this cost depends on the particular
parallel computer).
The parallel prefix operation in Figure 1 is realized
as follows. A virtual t-ary tree is constructed among
k be the processor that owns the region where
the next border zone event (NBZE) is about to take
place. This event is scheduled to occur at time T bz .
The parallel prefix operation calculates the minimum
among a set of p local NBZEs distributed in the p
processors (the minimum is stored in each processor).]
Parallel Simulation [processor i]
Initialisation;
while( not end condition )
begin-superstep
Simulate events with time less than T bz ;
Processor k reads the state of neighboring
regions;
end-superstep
Processor k simulates the occurrence of
the NBZE;
endwhile

Figure

1: Hyper-conservative simulation algorithm.
the p processors: from the leaves to the root the partial
are calculated, and then the absolute minimum
is distributed among the processors going from
the root to the leaves. The cost of this operation is
where the value of t depends on the parameters g and
l (e.g. for a small number of processors p it could be
more convenient to set
The efficiency of the algorithm in Figure 1 is improved
by attempting to simulate in parallel at most
border zone events per iteration. We explain this
procedure with an example. Let us assume a situation
with next BZ events e a
regions R a and R b respectively. In our
bg is the identifier of the element in
region R i which has scheduled the next BZ event e i
to occur at time t i . In addition, we define t
to be
the time at which an element i 0 (i
(R i
has scheduled a BZ event. We assume that
the elements are related due to the topology
of the system being simulated (e.g. neighboring atoms
in the ising-spin model described below). Note that t i 0
is not necessarily the time of the next BZ event in region
However, the simulation of both e i and e i 0
is restricted by the order relation between their respective
scheduled times t i and t i 0 . If t
must simulate e i before e i 0 , otherwise we first simulate
. Thus we simulate in parallel the two next
BZ events e a and e b only if t
Otherwise, we must process sequentially
more BZ events in the region with lesser t i until the
above condition is reached. For each new BZ event
processed in a region R i the non-BZ events in the time
interval between two consecutive BZ events have to
Parallel Simulation [region R a
Initialisation;
while( not end condition )
begin-superstep
end-superstep
begin-superstep
Simulate events e k with time t k
so that t k  t a and t k  t a 0
if (t a  t a 0 ) then
Simulate next BZ event e a ;
endif
end-superstep
endwhile

Figure

2: Conservative BSP simulation algorithm.
be simulated as well. This is described in the pseudo-code
for region R a shown in Figure 2. The operation
reads the value t a 0 of the element
a 0 stored in region R b .
4 Ising-spin systems
The ising-spin system is modeled as a
n \Theta
toroidal network. Every node i of the network is an
atom with magnetic spin value \Gamma1 or +1. Each atom
i attempts to change its spin value at discrete times
given by t is the time at
which the atom i has been currently simulated, and X
is a random variable with negative exponential distri-
bution. The new spin value of i is decided considering
the current spin values of its four neighbors. The goal
of the simulation is to process the occurrence of M
spin changes (events).
The sequential simulation of this system is trivial
since it is only necessary to deal with one type of
event and to use an efficient event-list to administer
the times t i(k+1) . Then the cost C 1 of processing each
event that takes place in the sequential algorithm is
O(logn) or even O(1) if a calendar queue were used [in
[2] it has been conjectured that the calendar queue has
O(1) cost under a work-load very similar to the one
produced by the ising-spin system]. The cost of the
sequential simulation of the whole system of n atoms
is then
In the case of the parallel simulation, the
toroidal network is divided into p p \Theta p p regions with
n=p) \Theta (a =
atoms each. For each
region there are a total of 4 (a \Gamma 1) atoms in the border
zone, i.e., f In each
region the same sequential event-list algorithm is applied
during the parallel phase, although it is executed
on a smaller number of atoms (n=p). The cost C p
of processing every event during the parallel phase is
event-list is used.
For each iteration, the cost of the parallel phase is
determined by the maximum number of events simulated
in any processor during that period. This number
is hard to determine. We optimistically assume
that on average a very similar number of events are
simulated by each processor. We are going to assume
that from the total of M simulated
in all of the parallel phases executed during the simu-
lation, a total of M are simulated by each
processor (this introduces a constant error since the
average maximum per iteration should be considered
here). Also we assume in our analysis that the M f bz
BZ events that take place are simulated sequentially
(hyper-conservative algorithm of Figure 1). Thus the
cost TP of the parallel simulation is given by
where TCS (p; g; l) is the cost in communication (g) and
synchronization (l) among the (p) processors generated
in each iteration.
To predict the performance of a BSP algorithm we
need to compare it with the fastest sequential algorithm
for the same problem. With this aim define the
speed-up S to be
shows that there exists a BSP algorithm
with total cost smaller than the
cost TS of the sequential alternative (i.e. in a BSP
algorithm we not only consider its computation cost
but also its cost in communication and synchronization
among processors).
For the case of the ising-spin model we have
Since C p  C 1 , we can replace C p by C 1 to obtain
an upper bound for TCS required to achieve S ? 1 ,
which expressed as
shows that the effect of the cost TCS is essentially absorbed
with f bz and C 1 . That is, given a particular
machine (characterized by its parameters g and l) we
can always achieve a speedup S ? 1 for a sufficiently
large problem (characterized by its parameters f bz and
For example, in an extreme situation of a system
with very low C 1 to be simulated on an inefficient ma-
chine, say very high g and l, the only way to achieve
increasing the parallel slackness (by
increasing a and/or reducing p) enough to reduce the
effect of TCS (p;
For the hyper-conservative algorithm given in Figure
1 the cost TCS is dominated by the parallel prefix
operation, i.e., p). For
Ising-spin system
f exp

Table

2: Results on an 8-processor IBM/SP2.
the more efficient algorithm shown in Figure 2 this
cost depends on the number q of different processors
to which every processor has to communicate with in
order to decide whether or not to simulate its next BZ
event, namely (for the
case of the 2D ising-spin model we have q  2).
For the ising-spin model we can estimate bounds for
TCS which ensure S ? 1. To this end we substitute in
to obtain
O( a log p )
for the hyper-conservative algorithm. On the other
hand, if for the less conservative algorithm we assume
that p BZ events are processed in each iteration,
namely
and
then we obtain the better bound
Note that the more restricted case C leads to
the bounds O(a) and O(a p) respectively. Given the
bounds for g and l shown in Table 1 we can see that
the restriction (upper bound) for TCS is possible to
satisfy in practice. For example, running the hyper-
conservative algorithm on a 2D array computer would
require one to adjust a so that a = O( p p).
In

Table

2 we show empirical results for S. We
obtained S using the running time of the O(log n) sequential
algorithm and the less conservative parallel
algorithm in Figure 2. In the column 4 (t) we show the
fraction of BZ events per iteration, where a value 1.0
means that one BZ event is simulated in each processor
(the best case). This column (t) shows the average
among all the iterations.
5 Hard-disk fluids
Our second example is more complex, it consists of
a two-dimensional box of size L \Theta L which contains
hard-disks evenly distributed in it. After a random
assignment of velocities, the (non-obvious) problem
consists of simulating a total of NDDC elastic disk-disk
collisions (DDC events) in a running time as small as
possible. In this section we show that similar bounds
for TCS (although with much higher constant factors)
are required to simulate these systems on a BSP computer
efficiently.
To achieve an efficient sequential running time, the
whole box is divided into p n c \Theta p n c cells of size oe \Theta oe
with
d, where d is the diameter of each
disk. The box is periodical in the sense that every
time a disk runs out of the box through a boundary
wall, it re-enters the box at the opposite point. The
neighborhood of a disk i whose center is located in the
cell c, is composed of the cell c itself and the eight cells
immediately (periodical) adjacent to c. We define m to
be the average number of disks per cell. Since oe  d,
a disk i can only collide with the 9 located
in the neighborhood of i. This reduces from O(n)
to O(logn) the cost associated with the simulation of
every DDC event that takes place since
regarded as a constant and we use a O(log n) event-list
to administer the pending events (collisions).
As the disks move freely between DDC events they
will eventually cross into neighboring cells. We regard
the instant when a disk i crosses from a cell c to a
neighboring cell c   as a virtual wall collision (VWC)
event. Then each time a VWC event takes place it is
necessary to consider the possible collisions among i
and the disks located in the new cells that become part
of the neighborhood of i, i.e., the cells immediately
adjacent to c   which are not adjacent to c (3 m disks
should be considered here). To consider the effect of
these events, we define  to be the average number of
VWC events that take place between two consecutive
DDC events. So the goal of simulating NDDC DDC
events actually involves processing the occurrence of
events. Note that 1
represents
the probability that the next event to take place, in
a given instant of the simulated time, is a DDC event
whereas
is the probability that the next event is
a VWC event.
To perform the simulation it is necessary to maintain
for each disk i, updated information of the time t
of all the possible collisions between i and the disks j
located in the neighborhood of i. It is also necessary
to periodically update the time when i crosses to a
neighboring cell through a virtual wall w. These computations
are done in a pair-wise manner by considering
only the positions and velocities of the two objects
involved in the event being calculated. The outcome
is a dynamic set of event-tuples
represents a disk j or a virtual wall w and e indicates
a DDC or VWC event. At initialization, the
first future events (event-tuples) are predicted for the
n disks of the system and then new future events are
successively calculated as the simulation advances to
the end, namely every time a disk suffer a DDC or
event.
Notice that only a subset of all the events calculated
for each disk i are the ones that really occur during the
simulation, and it is not obvious how to identify these
events in advance. Different methods to cope with this
problem have been proposed in the literature [6, 10,
14, 20]. However the common principle is to use an
efficient data structure to maintain an event-list where
the future events are stored until they are removed to
take place or they are invalidated by earlier events;
a DDC event E(t 1 ; stored in the event-list
becomes invalidated if another DDC event E(t
place during the simulation.
We assume here that only one event E is actually
maintained for each disk i (the one with minimal time
E(t)) in the event-list, and if this event E becomes invalidated
a new event E 0 for i is calculated considering
the complete neighborhood of i. In other words, after
every DDC and VWC, and when an invalid event is
retrieved from the event-list, new collisions are calculated
considering the 9 located in the neighborhood
(this implies a fairly slower sequential simulation
but also simplifies its implementation and analysis).
Note that the fraction of invalid events which are retrieved
as the "next event" is less than 15% [14], so
we neglect this effect in our analysis as well.
After initialization, the simulation enters a basic
cycle essentially composed of the following operations:
(i) picking the chronological next event from the event-
list, (ii) updating the state(s) of the disk(s) involved
in the current event, (iii) calculating new events for
this (these) disk(s) (one VWC event and several DDC
events) and (iv) inserting one event in the event-list
per disk involved in the current event. These operations
are cyclically performed until some end condition
is reached (i.e., the occurrence of a border zone event
in the case of the parallel simulation).
The running time of the sequential algorithm can be
estimated as follows. Constant factors are neglected
by considering that each of the operations of updating
the position or velocity of a disk i, and calculating
one DDC event for i, are all single operations with
cost O(1). Also, for every disk i it is necessary to consider
the disks located in the neighborhood of
while calculating new DDC events for i. Calculating
a VWC takes time O(1) as well. The cost associated
with the event-list is log n per event insertion whereas
retrieving the next event is negligible. Selecting the
event with minimal time for a disk i is also negligible
since this can be done as the new events are calcu-
lated. This gives the costs 2 (3 log n) and
for the simulation of each DDC and
VWC event that takes place respectively. Then the
overall cost of the simulation of each DDC that takes
place is
Notice that C 1 includes the cost of the  VWC events
that take place between two consecutive DDC events.
The total running T S of the sequential algorithm is
then
Using theory of hard-disk fluids, we have derivated
the following expression for  [13],
is the disk-area density of the
system. Calculating d TDDC
we obtained the ex-
pressionm
which even for extreme conditions (e.g., ae = 0:01 and
On the other
hand, the restriction oe  d imposes a lower bound for
m opt . Replacing
where in practice choosing
to an efficient simulation in terms of the total running
time and space used by the cells.
During the parallel phase every processor simulates
the evolution of the disks located in its own region. If
there are a total of p processors and an average of n=p
disks in every region, then by logarithmic property we
have
which is the average running time spent by each processor
computing the occurrence of two consecutive
DDC events in its region. We emphasize here that
hard-disk systems are by far more difficult to simulate
in parallel than ising-spin models. In particular,
it is necessary to cope with the problem that an event
scheduled for a particular disk may not occur at the
predicted time; this disk can be hit by a neighboring
disk in an earlier simulated time. This necessarily
leads one to deal with the possibility of "rollbacks"
where the whole simulation is re-started at some check
point passed without error.
For example, in the algorithm of Figure 2 after that
processor simulating region R a fetches time t a 0
from
region R b in superstep s, it might occur that processor
simulating region R b changes the value t a 0
in the next
invalidating in this way all the work
made by processor simulating R a during its parallel
phase (i.e., superstep s 1). To cope with this prob-
lem, we maintain an additional copy of the whole state
of the simulation. This state is an array with one entry
per disk. Each entry keeps disk's information such
as position, velocity and local simulation time. We
also use a single linked list to register each position of
the main state array that is modified during a complete
iteration (i.e., parallel phase plus synchronization
phase). If the above described problem occurs,
then we use the linked list to make the two state arrays
identical and repeat the iteration (s; s
a better estimation of t a 0 . See [13] for specific details
on the BSP simulation of hard-disk fluids (e.g.,
since a BSP computer is a distributed memory system
we maintain in each region a copy of the disks
located in the border zone of its neighboring regions,
thus the synchronization phase also involves the transference
of information among neighboring regions to
properly updated the states of these disk copies).
The regions to be simulated by each processor are
made up of a \Theta a cells with
Also we
define the BZ cells to be the 4 (a \Gamma 1) cells located
in the boundaries of every region. By studying the
probabilities of all the cases when a BZ event takes
place we can calculate f bz , which is a function of a
and other parameters of the hard-disk system. The
general expression for f bz is given by
where PDDCBZ and PVWCBZ are the probabilities of a
DDC and VWC event taking place in border zone res-
pectively. These probabilities can be calculated considering
that a given disk has the same probability of
being in an arbitrary cell and that its direction also
has the same probability. These calculations are a bit
involved because of the many cases to be considered.
Briefly, the expressions given below were obtained
by studying the two types of BZ events (DDC and
VWC) and the positions of the disk E(i) involved
in them. With probability
disk is located in a BZ cell whereas with probability
the disk is located in a cell
neighboring to a BZ cell. For a VWC event the disk
crosses to any cell with probability 1/4 whereas for
a DDC event the disk collides with a disk located in
any neighboring cell with probability 1/8. If the disk
is located in a cell neighboring to a BZ cell then with
probability the disk is located in a corner of
the box. In this case the probability of a BZ event is
5/8 for a DDC and 2/4 for a VWC. When the disk is
located out of the corner these last probabilities are
3/8 and 1/4 respectively. Similar considerations are
used when the disk is located inside a BZ cell for a
VWC. For DDC the probability of a BZ event is just
bz . By doing the weighted sum of all the cases we
have obtained
a 2
and
a
a 2
where Pm represents the probability of a DDC between
two disks located in the same cell. This probability
depends on the size of the cells oe, but for the purpose
of our analysis it is enough to say
f bz !a
O
'a
The running time TP of the parallel algorithm is
given by
where TPP is the time spent simulating NPE events
(DDC and VWC) during the parallel phase, and T SP
is the time spent in the synchronization phase simulating
one event plus the cost TCS associated with the
communication and synchronization among the pro-
cessors. Note that from the NPE events a total of
(1=(1+)) NPE are DDC events, and these events are
evenly distributed among the p processors, namely p
processors simultaneously simulate (1=(1 +)) NPE =p
DDC events. Then TPP is given by
''f bz
and TSP is given by
namely
and therefore
This last expression can be made more exigent for TCS
assuming to obtain
log n TCS
which shows that for a practical simulation with
1=a, the bound for TCS is similar
to those of the ising-spin system. On the other hand, if
we assume that p border zone events are simulated in
each iteration of the less conservative algorithm, then
we have
f bz
and
which leads to a bound similar to the one for the ising-
spin model as well. It is important to note that in the
calculations involved in the derivation of the speedup
S, we have been very conservative in the sense that
we are mixing BSP cost units with the ones defined
Hard-disk fluid
f exp

Table

3: Empirical results on the IBM/SP2.
by ourselves. Our basic unit of cost (updating disk
state or calculating a new event for a disk) is much
higher than the cost of each time step assumed for g
and l in the BSP model.
In

Table

3 we show empirical results for the hard-disk
fluid simulated with the less conservative algorithm
running on a IBM/SP2 parallel computer.
6 Final comments
In this paper we have derived upper bounds for
the cost of communication and synchronization among
processors in order to perform the efficient conservative
simulation of two system examples. We conclude
that it is possible to satisfy such bounds on current
parallel computers. Our empirical results confirm this
conclusion. We believe that the examples analyzed in
this paper exhibit sufficient generality and complexity
to be considered as representatives of a wide class of
systems where the events takes place randomly and
evenly distributed among their constituent elements.
The first example is a very simple system where
the links among neighboring regions (processors) are
maintained fixed during the whole simulation. How-
ever, the cost of each event processed in this system is
extremely low. This imposes harder requirements on
the cost of communication and synchronization (upper
bounds with much lower constant factors). The second
example is noticeably more complex because of the
dynamic nature of the system. Here the links among
regions change randomly during the simulation. Thus,
even for the hyper-conservative algorithm of Figure 1,
it is necessary to cope with the problem of roll-backs
since the scheduled events associated with each disk do
not necessarily occur at the predicted time. However
the simulated time progresses statistically at the same
rate in each region and the upper bounds for communication
and synchronization are similar to those of the
first simpler example (notice that the constant factors
are much higher in this second system example which
relaxes the requirements of these bounds).
All of our empirical results were obtained with the
less conservative algorithm shown in Figure 2 running
on an IBM/SP2. We could not obtain speedup S ? 1
with the hyper-conservative algorithm of Figure 1 running
under similar conditions. We emphasize, how-
ever, that these results were obtained for a particular
machine with fairly high g and l values. Only by increasing
the slackness to O(n 5 ) disks per processor
we obtained with the algorithm of

Figure

1 under the experiments described in Table 3.

Acknowledgements

The author has been supported by University of
Magallanes (Chile) and a Chilean scholarship.



--R

"Distributed simulation and time warp Part 1: Design of Colliding Pucks"
"Calendar queues: A fast O(1) priority queue implementation for the simulation event set problem"
"Discrete event simulation of object movement and interactions"
"The green BSP library"
"Per- formance of the colliding pucks simulation on the time warp operating system part 2: Asynchronous behavior & sectoring"
"An efficient algorithm for the hard-sphere problem"
"Efficient parallel simulations of dynamic ising spin systems"
"Efficient distributed event-driven simulations of multiple-loop networks"
"Simulating colliding rigid disks in parallel using bounded lag without time warp"
"How to simulate billiards and similars systems"
"Simulating billiards: Serially and in parallel"
"Direct BSP algorithms for parallel discrete-event simulation"
"Event-driven hard-particle molecular dynamics using bulk synchronous parallelism"
"Ef- ficient algorithms for many-body hard particle molecular dynamics"
"An empirical assessment of priority queues in event-driven molecular dynamics simulation"
"An object oriented C++ approach for discrete event simulation of complex and large systems of many moving ob- jects"
"General purpose parallel comput- ing"
"Parallel simulation of billiard balls using shared variables"

"The event scheduling problem in molecular dynamics simulation"
"Reduction of the event-list for molecular dynamic simulation"
"A bridging model for parallel com- putation"
"Distributed combat simulation and time warp: The model and its performance"
"Im- plementing a distributed combat simulation on the time warp operating system"
"Case studies in serial and parallel simulation"
--TR
Efficient parallel simulations of dynamic Ising spin systems
Calendar queues: a fast 0(1) priority queue implementation for the simulation event set problem
Implementing a distributed combat simulation on the Time Warp operating system
Efficient distributed event-driven simulations of multiple-loop networks
A bridging model for parallel computation
How to simulate billiards and similar systems
Discrete event simulation of object movement and interactions
General purpose parallel computing
Efficient algorithms for many-body hard particle molecular dynamics
An efficient algorithm for the hard-sphere problem
Parallel simulation of billiard balls using shared variables
An object oriented C++ approach for discrete event simulation of complex and large systems of many moving objects

--CTR
Wentong Cai , Emmanuelle Letertre , Stephen J. Turner, Dag consistent parallel simulation: a predictable and robust conservative algorithm, ACM SIGSIM Simulation Digest, v.27 n.1, p.178-181, July 1997
billiards;general purpose parallel computing;discrete-event simulation;homogeneous large-scale systems;discrete event simulation;ising-spin models;scalable computation;colliding hard-spheres;bulk-synchronous parallel model