Physics Computer Science
ComputerScience,1970-2008-Computer Science-2005.jsonl.gz.gz
Computer Science
8f6a160e1361cea45f22d59d9c990a86c5b7170c
Parallel Programming with Matrix Distributed Processing

Matrix Distributed Processing (MDP) is a C++ library for fast development of efficient parallel algorithms. It constitues the core of FermiQCD. MDP enables programmers to focus on algorithms, while parallelization is dealt with automatically and transparently. Here we present a brief overview of MDP and examples of applications in Computer Science (Cellular Automata), Engineering (PDE Solver) and Physics (Ising Model).


Introduction
Matrix Distributed Processing (MDP) [1][2] is a collection of classes and functions written in C++ for fast development of parallel algorithms such as solvers for partial differential equations, mesh-like algorithms, and various types of graph-based problems. These algorithms find frequent application in many sectors of physics, engineering, electronics and computational finance.
MDP includes: • a natural syntax for the algorithms that is transparent to the underlying parallelization; • parallelization based optimization algorithms implemented in MPI; • functions for linear algebra computations with support of a Maple-like syntax; • statistical functions; • fitting functions; • a Parallel SIMulator (PSIM) which allows MDP algorithms to run on single processor machines without MPI (uses fork for creating the processes and socket pairs for communication).
MDP was originally developed for and still constitutes the core of FermiQCD [3][4], a library for Lattice QCD computations. Lattice QCD is a Monte Carlo based numerical approach to the physics of composite particles made of quarks (e.g., protons and neutrons) and it is considered one of the most computationally intensive projects of modern physics.
FermiQCD, was developed by the author at the University of Southampton (UK) and Fermilab (Department of Energy) and it is currently used by other physics departments around the world [5][6] [7][8] [9].
While Lattice QCD is currently one of the main applications of MDP, its range of applicability is not limited to it as we show in the following examples.

Example: Parallel Game of Life
As a first example we show here how to write, in a few lines, a parallel program to find stable configurations for the Game Of Life. The game works as follows: • It is defined on a board (N × N cells) with periodic boundary conditions.
• The game consists of a series of iterations and, at each iteration, if a cell is dead and the number of beighbor cells who are alive is greater or equal 3, the cell dies; if a cell is alive and the number of beighbor cells who are alive is 2 or 3, the cell remains alive; otherwise the cell will be dead.
• Given a random starting configuration for the board, we want to determine its evolution and whether it reaches a stable configuration.
• Line 14 opens the parallel communication channels, line 37 closes them.
• Line 15 declares a 2D array containing the board size, L.
• Line 16 declares a lattice object (board) that represents the set of board sites and its topology. By default a mesh topology with toroid boundary conditions. Two arguments are passed to the constructor: the dimension, and size of the board. For each board site, it is possible to specify on which parallel process it is going to be allocated as well as alternative topologies.
• Line 17 declares a field of cells (C) on the board (the cells that live on it).
• Line 18 declares an auxiliary field that will be necessary for the computation.
• Line 19 declares a x a variable that represents a generic site of the board and will be used for looping over the sites and therefore the cells that live on them.
• Line 21 sets the initial configuration of the board by looping over all sites x and setting the corresponding cell C(x) to a random alive or dead. Each process loops only over the sites that are allocated locally by the process. Notice how MDP provides a local random number generator for each site of the board, board.random(x), which is vital in order to be able to reproduce results of stochastic algorithms.
• Lines 23-35 loop over the iterations of the game.
• Line 24 is performs a critical operation; it informs MDP that the value of the field C has changed and parallel communication is needed to synchronize those buffers that contain copies of non-local sites (synchronization).
• Lines 25-28 apply the rules of the game at each site. The new values of the states for the cells are stored in the auxiliary field newC. Each process loops over the local sites only. For each local site (x), x+0 represents the neighbor site when coordinate 0 is incremented by one, x-0 represents the neighbor site when coordinate 0 is decremented by one, x+1 represents the neighbor site when coordinate 1 is incremented by one, etc. 1 .
• Lines 29-33 perform two operations in a single parallel loop: the new states for the cells (newC) are copied back into C; the number of cells that have changed state are counted and the number is stored in diff. 1 This notation may appear bizzarre but it is nothing more than a sum ( x is a vector that represents a site on the board and the integer − → i represents a unit verson in direction i.
• Line 34 terminates the iteration if no cell has changed its state.
• Line 36 saves the board configuration in a file cells.dat.
Most of the parallel work is done by the constructor (mdp lattice) of the board which, from the board topology, determines how to partition it and determines the optimal communication patterns; by the constructor (mdp field) of the field of cells (C) which allocates the local cells and the buffers to store copies of non-local cells; and by the method update which performs communication to copy remote cell values into the local buffers.
We'll show in the next few examples that is equally easy to implement fields of any class of objects. In the next example we'll consider a field of matrices.

Example: Parallel PDE Solver
Consider here, as a different example that presents many similarities with the one above, the following Laplace equation: The initial conditions are ϕ initial (x) = 0. We will also assume that x i + L i = x i (torus topology). Solution: In order to solve eq. (1) we first discretize the Laplacian (∇ 2 = ∂ 2 0 + ∂ 2 1 + ∂ 2 2 ) and rewrite it as where µ is a unit vector in the discretized space in direction µ. Hence we solve for ϕ(x) and obtain the following recurrence relation The following is a typical MDP program that solves eq. (1) by recursively iterating eq. (4). Notice how the program is parallel but there are no explicit calls to communication functions: phi.save("field_phi.mdp"); // save field 24 mdp.close_wormholes(); // close communications 25 }

Notes:
• Line 04 declares the size of the box used to approximate the space L = {L 0 , L 1 , L 2 }.
• Line 05 declares a 3-dimensional lattice, called space, on the box L. MDP supports up to 10dimensional lattices. By default, a lattice object is a mesh with torus topology.
• Line 06 declares a site variable site x that will be used to loop over the lattice.
• Line 07 declares a field of 2 × 2 matrices, called phi, over the lattice space.
• Lines 11-12 initialize the field phi. Notice that phi is distributed over the parallel processes and forallsites is a parallel loop.
• Line 16 loops over all sites in parallel.
• Line 23 saves the field. Notice that all fields, including user-defined ones, inherit save and load methods from the basic mdp field class.

Example: Parallel Ising Model
As one more example of usage of MDP we report here a simple program for the Ising model. • Lines 3-4 declare a 1D lattice of 100 points (line).
• Line 5 declares a field of integers (spin) on the lattice.
• Line 7 sets the total magnetization M for this initial spin configuration.
• Line 9 sets the initial configuration: all field variables equal to 1.
• Line 14 computes the energy variation (dE) of each Markov Chain Monte Carlo (MCMC) step.
• Lines 15-16 perform the Monte Carlo accept-reject. If a change is accepted the spin at site x is flipped and the total magnetization M changes (line 16).
Note how at each MCMC step, first the code tries to flip the spins at even locations then, after it updates the lattice sites, it tries to flip the spins at odd locations (line 13). This guarantees computation results are independent on parallelization of the lattice line.
Since this even-odd distinction is common in many lattice algorithms, MDP stores all even lattice sites and all odd lattice sites close together. This speeds up loops over one of the two subsets and also speeds up communication. In fact, in this example, we are able to limit the synchronization (update) to the site of a given parity (line 18).

Linear Algebra and Other Tools
MDP includes a Linear Algebra package and other tools. Some of the most important classes are: • class mdp real, that should be used in place of float or double; • class mdp complex, for complex numbers; • class mdp array, for vectors and/or multidimensional tensors; • class mdp matrix, for any kind of complex rectangular matrix; • class mdp measure, for error propagation; • class mdp jackboot, a container for jackknife and bootstrap algorithms.
The most notable difference between our linear algebra package and other existing packages is its natural syntax.
For example: A=Random.SU (7); B=exp(A)+inv(A)*hermitian(A)+5; reads like A and B are matrices A is a random SU (7) matrix B = e A + A −1 A H + 5 · 1 Note that each matrix can be resized at will and is resized automatically when a value is assigned. MDP includes functions for fitting such as the Levenberger-Marquardt algorithm.

Lattice, Site, and Field
An mdp lattice is the class that describes the space on which fields are defined; it stores the topology of the space (by default that of a torus in d dimensions) and information about partioning of the space over the parallel processes. In MDP, a lattice is a graph, defined as a collection of points (lattice sites) connected by links (they specify the topology). Each site is uniquely mapped to one of the parallel processes.
The only restriction is that the graph must have a degree less than 20. From now on we will assume the default topology of a torus; therefore the lattice should be thought of as a mesh in d ≤ 10 dimensions.
The constructor class mdp lattice determines on which process to store each site, determines the neighbors of each site, and determines the sizes of the buffers where each process keeps copies of those sites that are non-local but are neighbors of the local sites.
The constructor also allocates a parallel random number generator so that each site of the lattice has its own independent random number generator. This is important for parallel Monte Carlo applications of MDP and ensures reproducibility of computations on different architectures. On each lattice it is possible to allocate fields. Some fields are built-in, for example This code declares a 10 × 10 × 10 lattice (cube) and a field (psi), that lives on the cube. The site variables of psi, psi(x) belong to class S (assuming x is an mdp site on the cube).
User-defined fields can be saved, loaded, and synchronized psi.save("filename"); psi.load("filename"); psi.update(); Synchronization means that all processes will perform MPI communications to make sure all buffers that contain copies of non-local site variables are updated with the proper values. The method update should be called immediately after the local site variables of a field have been changed.
Once a field object is declared in the field constructor, each process dynamically allocates memory for the buffers that store the copies of those sites that are non-local, but are neighbors of the local sites. These buffers are created in such as way to ensure optimal communication patterns.
Every time a field changes, for example in a parallel loop such as the program notifies the field that its values have been changed by calling phi.update(); The method update performs all required communication to copy site variables that need to be synchronized between each couple of overlapping processes.
Lattice sites are represented by objects of class mdp site. Site objects can be looped over in parallel loops (such as forallsites) but it also possible to explicitly address a specific site by specifying the site coordinates. Obviously, only the process that stores a site locally should address a specific site. Class mdp site has methods to check if a site is local, if it is non-local but a local copy is present, and which process stores the site locally.

Optimal Communication Patterns
In MDP, the lattice objects, according to the lattice topology and the parallel partitioning, determine the optimal way to store site variables in memory and performing parallel communication. This information is then used by the field method update that performs the synchronization of the field variables.
Note that MDP does not attempt to overlap computation and communication. By "optimal communication" pattern we mean that, under the assumptions below, the method update minimizes network traffic and data copies.
Current optimizations are based on the assumption that each processing node has one and only one network card and that the network is isotropic (latency and bandwidth for each couple of nodes is the same). This assumption is generally true for Ethernet and Myrinet clusters.
Communications are optimal in the sense that: • Each process retrieves all non-local site variables in a single send/recv for each process that contains sites which are neighbors of local sites.
• Two processes that do not store neighbor sites do not communicate with each other.
• No process is involved in more than a single send and a single receive at one time.
• Each process stores close in memory those copies of non-local site variables which are local to the same process. In this way synchronization does not require the use of additional buffers receiving buffers.
This technque is particularly efficient for algorithms that only require next-neighbor communication and run on a all-to-all network topology such as Ethernet or Myrinet.
It is possible, in principle, to change the above communication patterns to optimize communication for other network topologies.
Although communications are currently based on MPI, they do not make use of communication tags. It is therefore possible, in principle, to speed-up communication by using a faster tagless and bufferless protocol such as Myricom GM.
Our communication patterns have the effect of making communication almost insensitive to network latency, and communication speed is dominated by network bandwidth. Benchmarks are very much application dependent since parallel efficiency is greatly affected by the lattice size, by the amount of computation performed per site, processor speed and type of interconnection. In many typical applications, like the one described in the preceding example, the drop in efficiency is less than 10% up to 8 nodes (processes) and less than 20% up to 32 (our tests are usually performed on a cluster of Pentium 4 PCs (2.2GHz) running Linux and connected by Myrinet).

MDP and PSIM
For portability reasons, MDP is based on MPI. Nevertheless it is desirable to be able to run, test and debug MDP programs on a single node (with single or multi processor architecture) without having to install MPI. The latest version of MDP includes a Parallel SIMulator (PSIM). Despite the name, this is not quite a simulator but an emulator, i.e., a message passing library that uses local (unix/posix) socket pairs. PSIM is Objected Oriented and is not based on MPI.
When compiling with PSIM, the parallel processes are created at start-up by forking. The number of parallel processes is specified at runtime by passing the following command line argument to any MDP executable program, -PSIM_NPROCS=4 (this makes 4 parallel processes). PSIM also creates a communication log that can be used for debugging. MDP with PSIM has been tested on Linux, Mac and Windows (with cygwin).
For single processor node, using PSIM does not introduce any speed-up, but, for a small number of processes (2-16) it does not slow down the code either. For multi-process shared memory architecture, parallelization of PSIM should produce a speed-up comparable with MPI. We have not been yet performed such tests.
Moreover, PSIM should perform well on openMosix clusters as soon as openMosix starts supporting migratable sockets since all communication between the parallel processes will be done by the operating system. Unfortunately openMosix does not support migratable sockets yet.

Conclusions
MDP is an easy, powerful, and reliable tool for developing efficient parallel numerical applications. We have shown here examples in Computer Science (Cellular Automata), Mathematics (PDE solver) and Physics (Ising model).
MDP enables the programmer to focus on algorithm design while parallelization issues are dealt with automatically in a way transparent to the programmer.
The underlying communication functions are written in MPI but it is possible also to compile it without MPI using the built-in PSIM emulator which enables one to run parallel processes on a single node and/or single processor architectures, such as a PC. This is useful for testing and debugging purposes.
All of the features here described are fully functional and have been tested in real-life applications such as FermiQCD, developed by the University of Southampton (UK) and Fermilab (Department of Energy).