# Presentation: 183799

## CS 267: Applications of Parallel ComputersLecture 9:Sources of Parallelism and Locality in Simulation

- Horst D. Simon
**http://www.cs.berkeley.edu/~strive/cs267**

- CS267, Yelick

## Parallelism and Locality in Simulation

- Real world problems have parallelism and locality:
- Many objects operate independently of others.
- Objects often depend much more on nearby than distant objects.
- Dependence on distant objects can often be simplified.
- Scientific models may introduce more parallelism:
- When a continuous problem is discretized, temporal domain dependencies are generally limited to adjacent time steps.
- Far-field effects may be ignored or approximated in many cases.
- Many problems exhibit parallelism at multiple levels
- Example: circuits can be simulated at many levels, and within each there may be parallelism within and between subcircuits.

- CS267, Yelick

**Notes:**
Not all problems solved on parallel machines involve simulation, but many do

## Example: Circuit Simulation 

- Circuits are simulated at many different levels

| Level | Primitives | Examples |
| --- | --- | --- |
| Instruction level | Instructions | SimOS, SPIM |
| Cycle level | Functional units |           VIRAM-p |
| Register Transfer Level (RTL) | Register, counter, MUX | VHDL |
| Gate Level | Gate, flip-flop, memory cell |           Thor |
| Switch level | Ideal transistor | Cosmos |
| Circuit level | Resistors, capacitors, etc. | Spice |
| Device level | Electrons, silicon |  |

- CS267, Yelick

## Basic Kinds of Simulation

- Discrete event systems:
- Examples: “Game of Life,” logic level circuit simulation.
- Particle systems:
- Examples: billiard balls, semiconductor device simulation, galaxies.
- Lumped variables depending on continuous parameters:
- ODEs, e.g., circuit simulation (Spice), structural mechanics, chemical kinetics.
- Continuous variables depending on continuous parameters:
- PDEs, e.g., heat, elasticity, electrostatics.
- A given phenomenon can be modeled at multiple levels.
- Many simulations combine more than one of these techniques.

- CS267, Yelick

## Outline

- Discrete event systems
- Time and space are discrete
- Particle systems
- Important special case of lumped systems
- Ordinary Differential Equations (ODEs)
- Lumped systems
- Location/entities are discrete, time is continuous
- Partial Different Equations (PDEs)
- Time and space are continuous
- Next lecture

- discrete

- continuous

- CS267, Yelick

## Discrete Event Systems

- CS267, Yelick

## Discrete Event Systems

- Systems are represented as:
- finite set of variables.
- the set of all variable values at a given time is called the state.
- each variable is updated by computing a transition function depending on the other variables.
- System may be:
- synchronous: at each discrete timestep evaluate all transition functions; also called a finite state machine.
- asynchronous: transition functions are evaluated only if the inputs change, based on an “event” from another part of the system; also called event driven simulation.
- Example: The “game of life:”
- Also known as Sharks and Fish
- http://www.cs.berkeley.edu/~strive/cs267/hw2.html
- Space divided into cells, rules govern cell contents at each step

- CS267, Yelick

## A Model Problem: Sharks and Fish

- Illustration of parallel programming
- Original version (discrete event only) proposed by Geoffrey Fox
- Called WATOR
- Sharks and fish living in a 2D toroidal ocean
- We can imagine several variation to show different physical phenomenon
- Basic idea: sharks and fish living in an ocean
- rules for movement
- breeding, eating, and death
- forces in the ocean
- forces between sea creatures

- CS267, Yelick

## Sharks and Fish as Discrete Event System

- Ocean modeled as a 2D toroidal grid
- Each cell occupied by at most one sea creature

- CS267, Yelick

## Fish-only: the Game of Life

- An new fish is born if
- a cell is empty
- exactly 3 (of 8) neighbors contain fish
- A fish dies (of overcrowding) if
- cell contains a fish
- 4 or more neighboring cells are full
- A fish dies (of loneliness) if
- cell contains a fish
- less than 2 neighboring cells are full
- Other configurations are stable

- The original Wator problem adds fish-eating sharks

- CS267, Yelick

## Parallelism in Sharks and Fish

- The activities in this system are discrete events
- The simulation is synchronous
- use two copies of the grid (old and new)
- the value of each new grid cell in new depends only on the 9 cells (itself plus neighbors) in old grid
- Each grid cell update is independent: reordering or parallelism OK
- simulation proceeds in timesteps, where (logically) each cell is evaluated at every timestep

- old ocean

- new ocean

- CS267, Yelick

## Parallelism in Sharks and Fish

- Parallelism is straightforward
- ocean is regular data structure
- even decomposition across processors gives load balance
- Locality is achieved by using large patches of the ocean
- boundary values from neighboring patches are needed
- although, there isn’t much reuse...
- Advanced optimization: visit only occupied cells (and neighbors)  load balance is more difficult

- CS267, Yelick

## Synchronous Circuit Simulation

- Circuit is a graph made up of subcircuits connected by wires
- Component simulations need to interact if they share a wire.
- Data structure is irregular (graph) of subcircuits.
- Parallel algorithm is timing-driven or synchronous:
- Evaluate all components at every timestep (determined by known circuit delay)
- Graph partitioning assigns subgraphs to processors (NP-complete)
- Determines parallelism and locality.
- Attempts to evenly distribute subgraphs to nodes  (load balance).
- Attempts to minimize edge crossing (minimize communication).

- CS267, Yelick

**Notes:**
Nodes and edges in graph may be weighted
Optimal graph partitioning is NP-complete (define for non-CS audience)
Good graph partitioning techniques exist
More on this in a later lecture

## Asynchronous Simulation

- Synchronous simulations may waste time:
- Simulate even when the inputs do not change,.
- Asynchronous simulations update only when an event arrives from another component:
- No global time steps, but individual events contain time stamp.
- Example: Game of life in loosely connected ponds (don’t simulate empty ponds).
- Example: Circuit simulation with delays (events are gates changing).
- Example: Traffic simulation (events are cars changing lanes, etc.).
- Asynchronous is more efficient, but harder to parallelize
- In MPI, events are naturally implemented as messages, but how do you know when to execute a “receive”?

- CS267, Yelick

## Scheduling Asynchronous Circuit Simulation

- Conservative:
- Only simulate up to (and including) the minimum time stamp of inputs.
- May need deadlock detection if there are cycles in graph, or else “null messages”.
- Example: Pthor circuit simulator in Splash1 from Stanford.
- Speculative (or Optimistic):
- Assume no new inputs will arrive and keep simulating.
- May need to backup if assumption wrong.
- Example: Timewarp [D. Jefferson], Parswec [Wen,Yelick].
- Optimizing load balance and locality is difficult:
- Locality means putting tightly coupled subcircuit on one processor.
- Since “active” part of circuit likely to be in a tightly coupled subcircuit, this may be bad for load balance.

- CS267, Yelick

## Summary of Discrete Even Simulations

- Model of the world is discrete
- Both time and space
- Approach
- Decompose domain, i.e., set of objects
- Run each component ahead using
- Synchronous: communicate at end of each timestep
- Asynchronous: communicate on-demand
- Conservative scheduling – wait for inputs
- Speculative scheduling – assume no inputs, roll back if necessary

- CS267, Yelick

## Particle Systems

- CS267, Yelick

## Particle Systems

- A particle system has
- a finite number of particles.
- moving in space according to Newton’s Laws (i.e. F = ma).
- time is continuous.
- Examples:
- stars in space with laws of gravity.
- electron beam and ion beam semiconductor manufacturing.
- atoms in a molecule with electrostatic forces.
- neutrons in a fission reactor.
- cars on a freeway with Newton’s laws plus model of driver and engine.
- Many simulations combine particle simulation techniques with some discrete event techniques (e.g., Sharks and Fish).

- CS267, Yelick

## Forces in Particle Systems

- Force on each particle decomposed into near and far:

- External force
- ocean current to sharks and fish world (S&F 1).
- externally imposed electric field in electron beam.
- Nearby force
- sharks attracted to eat nearby fish (S&F 5).
- balls on a billiard table bounce off of each other.
- Van der Wals forces in fluid (1/r6).
- Far-field force
- fish attract other fish by gravity-like (1/r2 ) force (S&F 2).
- gravity, electrostatics
- forces governed by elliptic PDE.

- force = external_force + nearby_force + far_field_force

- CS267, Yelick

## Parallelism in External Forces

- External forces are the simplest to implement.
- The force on each particle is independent of other particles.
- Called “embarrassingly parallel”.
- Evenly distribute particles on processors
- Any even distribution works.
- Locality is not an issue, no communication.
- For each particle on processor, apply the external force.

- CS267, Yelick

## Parallelism in Nearby Forces

- Nearby forces require interaction and therefore communication.
- Force may depend on other nearby particles:
- Example: collisions.
- simplest algorithm is O(n2): look at all pairs to see if they collide.
- Usual parallel model is decomposition* of physical domain:
-  O(n2/p) particles per processor if evenly distributed.

- Need to check for collisions between regions

- *often called “domain decomposition,” but the term also refers to a numerical technique.

- CS267, Yelick

## Parallelism in Nearby Forces

**Parallelism in Nearby Forces**

- Challenge 1: interactions of particles near processor boundary:
- need to communicate particles near boundary to neighboring processors.
- surface to volume effect means low communication.
- Which communicates less: squares (as below) or slabs?

- Communicate particles in boundary region to neighbors

- CS267, Yelick

## Parallelism in Nearby Forces

- Challenge 2: load imbalance, if particles cluster:
- galaxies, electrons hitting a device wall.
- To reduce load imbalance, divide space unevenly.
- Each region contains roughly equal number of particles.
- Quad-tree in 2D, oct-tree in 3D.

- Example: each square contains at most 3 particles

- See: _[http://njord.umiacs.umd.edu:1601/users/brabec/quadtree/points/prquad.html](http://njord.umiacs.umd.edu:1601/users/brabec/quadtree/points/prquad.html)_

- CS267, Yelick

## Parallelism in Far-Field Forces

- Far-field forces involve all-to-all interaction and therefore communication.
- Force depends on all other particles:
- Examples: gravity, protein folding
- Simplest algorithm is O(n2) as in S&F 2, 4, 5.
- Just decomposing space does not help since every particle needs to “visit” every other particle.
- Use more clever algorithms to beat O(n2).

- Implement by rotating particle sets.
- Keeps processors busy
- All processor eventually see all particles

- CS267, Yelick

## Far-field Forces: Particle-Mesh Methods

- Based on approximation:
- Superimpose a regular mesh.
- “Move” particles to nearest grid point.
- Exploit fact that the far-field force satisfies a PDE that is easy to solve on a regular mesh:
- FFT, multigrid (described in future lecture)
- Accuracy depends on the fineness of the grid is and the uniformity of the particle distribution.

- 1) Particles are moved to mesh (scatter)
- 2) Solve mesh problem
- 3) Forces are interpolated at particles (gather)

- CS267, Yelick

**Notes:**
Data structures: a 1D array of particles (each with a location)
                           a 2D (or 3D) array representing the mesh
Scatter: store particle value (say mass) at nearby mesh point
Gather: read forces from mesh, using particle location as index (as in scatter)

## Far-field forces: Tree Decomposition

- Based on approximation.
- Forces from  group of far-away particles “simplified” -- resembles a single large particle.
- Use tree; each node contains an approximation of descendants.
- O(n log n) or O(n) instead of O(n2).
- Several Algorithms
- Barnes-Hut.
- Fast multipole method (FMM)
-      of Greengard/Rohklin.
- Anderson’s method.
- Discussed in later lecture.

- CS267, Yelick

## Summary of Particle Methods

- Model contains discrete entities, namely, particles
- Time is continuous – is discretized to solve
- Simulation follows particles through timesteps
- All-pairs algorithm is simple, but inefficient, O(n2)
- Particle-mesh methods approximates by moving particles
- Tree-based algorithms approximate by treating set of particles as a group, when far away
- May think of this as a special case of a “lumped” system

- CS267, Yelick

## Lumped SystemsODEs

- CS267, Yelick

## System of Lumped Variables

- Many systems are approximated by
- System of “lumped” variables.
- Each depends on continuous parameter (usually time).
- Example -- circuit:
- approximate as graph.
- wires are edges.
- nodes are connections between 2 or more wires.
- each edge has resistor, capacitor, inductor or voltage source.
- system is “lumped” because we are not computing the voltage/current at every point in space along a wire, just endpoints.
- Variables related by Ohm’s Law, Kirchoff’s Laws, etc.
- Forms a system of ordinary differential equations (ODEs).
- Differentiated with respect to time

- CS267, Yelick

## Circuit Example

- State of the system is represented by
- v_n(t)  node voltages
- i_b(t) branch currents	all at time t
- v_b(t) branch voltages
- Equations include
- Kirchoff’s current
- Kirchoff’s voltage
- Ohm’s law
- Capacitance
- Inductance
- Write as single large system of ODEs (possibly with constraints).

- 0	A	0	v_n	0
- A’	0	-I     *	i_b    =	S
- 0	R	-I	v_b	0
- 0	-I       C*d/dt		0
- 0	L*d/dt	I		0

- CS267, Yelick

## Structural Analysis Example

- Another example is structural analysis in civil engineering:
- Variables are displacement of points in a building.
- Newton’s and Hook’s (spring) laws apply.
- Static modeling: exert force and determine displacement.
- Dynamic modeling: apply continuous force (earthquake).
- Eigenvalue problem: do the resonant modes of the building match an earthquake

- CS267, Yelick

## Solving ODEs

- In these examples, and most others, the matrices are sparse:
- i.e., most array elements are 0.
- neither store nor compute on these 0’s.
- Given a set of ODEs, two kinds of questions are:
- Compute the values of the variables at some time t
- Explicit methods
- Implicit methods
- Compute modes of vibration
- Eigenvalue problems

- CS267, Yelick

**Notes:**
Reason is a locality property from physics.

## Solving ODEs: Explicit Methods

- Assume ODE is x’(t) = f(x) = A*x, where A is a sparse matrix
- Compute x(i*dt) = x[i]
-          at i=0,1,2,...
- Approximate x’(i*dt)
-       x[i+1]=x[i] + dt*slope
- Explicit methods, e.g., (Forward) Euler’s method.
- Approximate x’(t)=A*x by (x[i+1] - x[i] )/dt = A*x[i].
- x[i+1] = x[i]+dt*A*x[i],  i.e. sparse matrix-vector multiplication.
- Tradeoffs:
- Simple algorithm: sparse matrix vector multiply.
- Stability problems: May need to take very small time steps, especially if system is stiff (i.e. can change rapidly).

- Use slope at x[i]

- CS267, Yelick

## Solving ODEs: Implicit Methods

- Assume ODE is x’(t) = f(x) = A*x, where A is a sparse matrix
- Compute x(i*dt) = x[i]
-          at i=0,1,2,...
- Approximate x’(i*dt)
-        x[i+1]=x[i] + dt*slope
- Implicit method, e.g., Backward Euler solve:
- Approximate x’(t)=A*x by (x[i+1] - x[i] )/dt = A*x[i+1].
- (I - dt*A)*x[i+1] = x[i],  i.e. we need to solve a sparse linear system of equations.
- Trade-offs:
- Larger timestep possible: especially for stiff problems
- More difficult algorithm: need to do a sparse solve at each step

- Use slope at x[i+1]

- CS267, Yelick

## Solving ODEs: Eigensolvers

- Computing modes of vibration: finding eigenvalues and eigenvectors.
- Seek solution of x’’(t) = A*x of form x(t) = sin(f*t)*x0, where x0 is a constant vector.
- Plug in to get  -f2 *x0 = A*x0, so that –f2  is an eigenvalue and x0 is an eigenvector of A.
- Solution schemes reduce either to sparse-matrix multiplication, or solving sparse linear systems.

- CS267, Yelick

## ODEs and Sparse Matrices

- All these reduce to sparse matrix problems
- Explicit: sparse matrix-vector multiplication.
- Implicit: solve a sparse linear system
- direct solvers (Gaussian elimination).
- iterative solvers (use sparse matrix-vector multiplication).
- Eigenvalue/vector algorithms may also be explicit or implicit.

- CS267, Yelick

## Parallel Sparse Matrix-vector multiplication

- y = A*x, where A is a sparse  n x n matrix
- Questions
- which processors store
- y[i], x[i], and A[i,j]
- which processors compute
- y[i] = sum (from 1 to n) A[i,j] * x[j]
-             = (row i of A) * x          ... a sparse dot product
- Partitioning
- Partition index set {1,...,n} = N1 + N2 + ... + Np.
- For all i in Nk, Processor k stores y[i], x[i], and row i of A
- For all i in Nk, Processor k computes y[i] = (row i of A) * x
- “owner computes” rule: Processor k compute the y[i]s it owns.

- i: [j1,v1], [j2,v2],...

- x

- y

- Po
- P1
- P2
- P3

- Most problematic

- CS267, Yelick

## Matrix Reordering via Graph Partitioning

- “Ideal” matrix structure for parallelism: block diagonal
- p (number of processors) blocks, can all be computed locally.
- few non-zeros outside these blocks, which require communication.
- Can we reorder the rows/columns to achieve this?

- P0
- P1
- P2
- P3
- P4

**=**

*****

- P0    P1   P2   P3  P4

- CS267, Yelick

## Goals of Reordering

- Performance goals
- balance load (how is load measured?).
- balance storage (how much does each processor store?).
- minimize communication (how much is communicated?).
- Some algorithms reorder for other reasons
- Reduce # nonzeros in answer (fill)
- Improve numerical properties

- CS267, Yelick

## Graph Partitioning and Sparse Matrices 

- 1    1     1                      1
- 2    1     1             1               1
- 3                   1     1               1
- 4           1      1     1       1
- 5    1                    1       1      1
- 6            1     1              1      1

-   1     2      3      4      5      6

- Relationship between matrix and graph

- A “good” partition of the graph has
- equal (weighted) number of nodes in each part (load and storage balance).
- minimum number of edges crossing between (minimize communication).
- Reorder the rows/columns by putting all nodes in one partition together.

- CS267, Yelick

## Implicit Methods and Eigenproblems

- Direct methods (Gaussian elimination)
- Called LU Decomposition, because we factor A = L*U.
- Future lectures will consider both dense and sparse cases.
- More complicated than sparse-matrix vector multiplication.
- Iterative solvers
- Will discuss several of these in future.
- Jacobi, Successive over-relaxation (SOR) , Conjugate Gradient (CG), Multigrid,...
- Most have sparse-matrix-vector multiplication in kernel.
- Eigenproblems
- Future lectures will discuss dense and sparse cases.
- Also depend on sparse-matrix-vector multiplication,  direct methods.

- CS267, Yelick