MPI: Parallel Pectoral Inversion
Jumanazarov Mardonbek
Plan:
The importance of the Message Passing Interface (MPI) standard is that it allows a program to access additional computing nodes and, therefore, perform increasingly larger tasks by adding more nodes to the simulation. The term "message passing" refers to the ability to easily send messages from one process to another. MPI is widely used in high-performance computing. In many scientific fields, the use of supercomputers necessitates the implementation of MPI.
MPI was launched as an open standard in 1994 and within months became the dominant library-oriented parallel computing language. Since 1994, MPI has led to scientific breakthroughs from physics to machine learning and self-driving cars! Several MPI implementations are currently in widespread use. The two most common are MPICH from Argonne National Laboratories and OpenMPI. Hardware vendors often have customized versions of one of these two implementations for their platforms. The MPI standard, currently at version 3.1 in 2015, continues to evolve and change.
In this chapter, we'll show you how to implement MPI in your application. We'll start with a simple MPI program and then move on to a more complex example of how to link together separate computational grids on separate processes by passing boundary information. We'll touch on several advanced techniques that are essential for well-written MPI programs, such as creating application-specific MPI data types and using MPI's Cartesian topology features. Finally, we'll introduce a combination of MPI with OpenMP (MPI plus OpenMPI) and vectorization to achieve multiple levels of parallelism.
NOTE: We recommend referencing the examples in this chapter at https://github.com/EssentialsofParallelComputing/Chapter8.
1. MPI Program Basics
In this section, we'll cover the basics needed for a minimal MPI program. Some of these core requirements are defined by the MPI standard, while others are common across most MPI implementations. The basic structure and operation of MPI have remained remarkably consistent since the first standard.
To begin with, MPI is a fully library-oriented language. It requires no special compiler or operating system support. All MPI programs have a basic structure and flow, as shown in Figure 8.1. MPI always begins with a call to MPI_Init at the beginning of the program and MPI_Finalize at program exit. This contrasts with OpenMP, as discussed in Chapter 7, which doesn't require special startup and shutdown commands and simply places parallel directives around key loops.
After a parallel MPI program is written, it is compiled with an include file and library. It is then executed using a special launcher that creates parallel processes between and within nodes.
1. MPI Program Basics
Fig. 8.1 The MPI approach is library-based. Simply compile, linking to the MPI library, and run it using a dedicated parallel execution program.
1. MPI Program Basics
Basic MPI Function Calls for Every MPI Program
The basic MPI function calls include MPI_Init and MPI_Finalize. The MPI_Init call must occur immediately after program startup, and arguments from the main procedure must be passed to the initialization call. Typical calls look like this and may or may not have a return variable:
iret = MPI_Init(&argc, &argv);
iret = MPI_Finalize();
Most programs will require the number of processes and the process rank within a group capable of communicating, called a communicator. One of the main functions of the MPI interface is to start remote processes and link them so that messages can be passed between them. By default, the communicator is MPI_COMM_WORLD, which is configured at the beginning of each parallel job using the MPI_Init function. Let's pause and look at a few definitions:
The calls to obtain these important variables are:
iret = MPI_Comm_rank(MPI_COMM_WORLD, &rank);
iret = MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
1. MPI Program Basics
Compiler Wrappers for Simpler MPI Programs
Although MPI is a library, we can treat it as a compiler through MPI compiler wrappers. This simplifies building MPI applications because you don't need to know which libraries are required or where they are located. They are especially convenient for small MPI applications. There are compiler wrappers for each programming language:
Using these wrappers is optional. If you don't use compiler wrappers, they can still be useful for identifying the compiler flags required to build your application. The mpicc command has options that provide this information. The specified options for your MPI can be found using the man mpicc command. For the two most popular MPI implementations listed below, we list the console options for mpicc, mpicxx, and mpifort.
– --showme;
– --showme:compile;
– --showme:link.
– -show;
– -compile_info;
– -link_info.
1. MPI Program Basics
Using Parallel Startup Commands
Starting parallel processes for MPI is a complex operation handled by a dedicated command. Initially, this command was often mpirun. However, with the release of the MPI 2.0 standard in 1997, mpiexec was recommended as the startup command in an attempt to provide greater portability. However, this standardization attempt was not entirely successful, and today several names are used for the startup command:
Most MPI startup commands use the -n option, which specifies the number of processes, but others can take the value -np. Given the complexity of modern computer node architectures, startup commands have myriad options for affinity, placement, and environment (some of which we discuss in Chapter 14). These options vary across MPI implementations and even across releases of their MPI libraries. The simplicity of the options available in the initial startup commands has devolved into an endless morass of options that haven't yet fully stabilized. Fortunately, most of these options can be ignored for a novice MPI user, but they are essential for advanced use and fine-tuning.
1. MPI Program Basics
A Minimal Working Example of an MPI Program
Now that we've mastered all the basic components, we can combine them into a minimal working example, shown in Listing 8.1: we run a parallel job and print the rank and number of processes for each process. In the call to obtain the rank and size, we use the MPI_COMM_WORLD variable, which represents the group of all MPI processes and is predefined in the MPI header file. Note that the output can be printed in any order; the MPI program leaves it to the operating system to decide when and how it is printed.
1. MPI Program Basics
Listing 8.2 defines a simple makefile for building this example using MPI compiler wrappers. In this case, we use the mpicc wrapper to specify the location of the mpi.h include file and the MPI library.
For more complex builds on different systems, you may prefer CMake. The following listing shows the CMakeLists.txt file for this program.
1. MPI Program Basics
1. MPI Program Basics
Now, using the CMake build system, let's configure, build, and then run the test using the following commands:
cmake
make
make test
The printf command prints the output in any order. Finally, to clean up after the run, use the following commands:
make clean
make distclean
2. Send and receive commands for process-to-process data exchange
The core of the message-passing approach is sending messages from point to point, or perhaps more accurately, from process to process. The whole point of parallel processing is to coordinate work. To do this, you need to send messages for control or to distribute work. We'll show you how these messages are constructed and sent correctly. There are many variations of point-to-point procedures; we'll cover those recommended for most situations.
Figure 8.2 shows the components of a message. At either end of the system, there must be a mailbox. The size of the mailbox is important. The sending side knows the size of the message, but the receiving side doesn't. To ensure storage space for the message, it's usually best to post the received message first. This avoids message latency, which occurs when the receiving process needs to allocate temporary storage space until the received message is sent and can copy it to the proper location. By analogy, if a received message (mailbox) hasn't been sent (it's not there), the postman must linger, constantly checking until someone provides it. Sending the received message from the start avoids the possibility of running out of memory space on the receiving end for the temporary buffer allocated for storing the message.
2. Send and receive commands for process-to-process data exchange
Fig. 8.2 A message in MPI always consists of a memory pointer, a counter, and a type. The envelope contains an address consisting of a rank, tag, and communication group, as well as an internal MPI context.
2. Send and receive commands for process-to-process data exchange
The message itself always consists of a triplet at both ends: a pointer to a memory buffer, a counter, and a type. The send type and receive type can have different types and counters. The rationale for using types and counters is to allow type conversion between processes at the source and destination. This allows the message to be converted to a different form at the receiving end. In a heterogeneous environment, this may mean converting from a right-handed encoding to a left-handed encoding[1], which accounts for low-level differences in the byte order of data stored by different hardware vendors. Furthermore, the receive size can be larger than the send size. This allows the receiver to request the size of the sent data to properly manipulate the message. However, the received size cannot be smaller than the sent size, as this would result in writing past the end of the buffer.
The envelope also consists of a triplet. It identifies the message's sender, recipient, and message identifier, preventing the confusion of multiple messages. The triplet consists of a rank, a tag, and a communication group. A rank is assigned to the specified communication group. The tag helps the programmer and MPI identify the message itself and the reception to which it belongs. In MPI, the tag is a convenience. It can be set to MPI_ANY_TAG if an explicit tag number is not required. MPI uses a context created within the library to properly separate messages. A message is complete if both the communicator and the tag match.
2. Send and receive commands for process-to-process data exchange
NOTE: One of the strengths of the message-passing approach is the memory model. Each process has clear ownership of its data, as well as control and synchronization over changes to that data. You are guaranteed that another process cannot modify your memory while your back is turned.
Now let's try an MPI program with a simple send/receive. We must send data on one process and receive data on another. These cross-process calls can be made in different ways (Figure 8.3). Some combinations of basic blocking sends and receives are unsafe and can hang, such as the two combinations on the left in Figure 8.3. The third combination requires careful programming using conditional blocks. The method on the far right is one of several safe methods for scheduling data exchanges using non-blocking sends and receives. These are also called asynchronous or immediate calls, which explains the I (for "immediate") preceding the send and receive keywords (this case is shown on the far right of the figure).
[1] That is, from little-endian encoding, in which the least significant bit is on the right, to big-endian encoding, in which the least significant bit is on the left. – Translator's note.
2. Send and receive commands for process-to-process data exchange
Fig. 8.3 Blocking sending and receiving operations is difficult to get right. It is much safer and faster to use non-blocking or immediate forms of sending and receiving operations and then wait for completion.
2. Send and receive commands for process-to-process data exchange
The most basic MPI send and receive functions are MPI_Send and MPI_Recv. The basic send and receive functions have the following prototypes:
MPI_Send(void *data, int count, MPI_Datatype datatype, int dest, int tag, MPI_COMM comm)
MPI_Recv(void *data, int count, MPI_Datatype datatype, int source, int tag, MPI_COMM comm, MPI_Status *status)
Now let's walk through each of the four cases in Figure 8.3 to understand why some hang and others work fine. We'll start with MPI_Send and MPI_Receive, which were shown in the previous function prototypes and in the leftmost example in the figure. Both of these procedures are blocking. Blocking means they don't return until a specific condition is met. In the case of these two calls, the return condition is that the buffer is safe to use again. When sending, the buffer must be read and no longer needed. When receiving, the buffer must be full. If both processes in the message are blocking, a situation called a hang can occur. A hang occurs when one or more processes are waiting for an event that can never occur.
2. Send and receive commands for process-to-process data exchange
Example: A Hanging Blocking Send/Receive Program
This example highlights a common problem in concurrent programming. You must always be vigilant to avoid situations that could cause a deadlock. To avoid this, the following listing clearly illustrates how this can happen.
2. Send and receive commands for process-to-process data exchange
2. Send and receive commands for process-to-process data exchange
The tag and rank of a communication partner are calculated using integer and modular arithmetic, which pairs the tags for each send and receive and obtains the rank of the other member of the pair. Receives are then sent for each process with its partner. These receives are blocking and do not complete (return) until the buffer is full. Since send is called only after the receives are complete, the program hangs. Note that we wrote the send and receive calls without if statements (conditional blocks), based on rank. In parallel code, conditional blocks are the source of many errors, so they should generally be avoided.
2. Send and receive commands for process-to-process data exchange
Let's try reversing the order of sending and receiving. We'll show the modified lines of code from the original listing of the previous example in the listing below.
Will this example fail? It depends. The send call returns after the data buffer from sending is used. In most MPI implementations, data will be copied to pre-allocated buffers on the sender or receiver if the size is small enough. In this case, the send completes, and receive is called. If the message is large, the sender waits until receive allocates a buffer to hold the message before returning. But receive is never called, so the program hangs. We could alternate sending and receiving by rank to avoid hangs. For this option, we would use a conditional block, as shown in the listing below.
2. Send and receive commands for process-to-process data exchange
But in more complex data exchanges, this is difficult to do correctly and requires careful use of conditional blocks. A more appropriate implementation is to use the MPI_Sendrecv function call, as shown in the listing below. By using this call, you delegate responsibility for properly executing data exchange to the MPI library. For the programmer, this is a rather favorable trade-off.
2. Send and receive commands for process-to-process data exchange
The MPI_Sendrecv call is a good example of the benefits of using collective communication calls, which we'll present in Section 8.3. It's recommended to use collective communication calls whenever possible, as they delegate responsibility for preventing deadlocks and stalls, as well as ensuring good performance, to the MPI library.
As an alternative to the blocking communication calls from previous examples, we'll consider using MPI_Isend and MPI_Irecv in Listing 8.7. These are called immediate (I) versions because they return immediately. They're often referred to as asynchronous or non-blocking calls. Asynchrony means the call initiates an operation but doesn't wait for the work to complete.
2. Send and receive commands for process-to-process data exchange
Each process waits for the message to complete in MPI_Waitall on line 31 of the listing. You should also see a noticeable performance improvement by reducing the number of blocks for each send and receive call to a single MPI_Waitall. However, you should be careful not to modify the send buffer or access the receive buffer until the operation completes. There are other combinations that work. Let's look at the following listing, which uses one of these possibilities.
2. Send and receive commands for process-to-process data exchange
2. Send and receive commands for process-to-process data exchange
We begin the data exchange with an asynchronous send and then block with a blocking receive. Once the blocking receive completes, the process can continue even if the send is not complete. You must still release the request handle with MPI_Request_free or as a side effect of calling MPI_Wait or MPI_Test. This is done to prevent memory leaks. You can also call MPI_Request_free immediately after sending MPI_Isend.
Other send/receive modes are useful in special situations. The modes are designated by a one- or two-letter prefix similar to that found in the immediate mode, as follows:
2. Send and receive commands for process-to-process data exchange
The list of predefined MPI data types for C is extensive; the data types map to almost all C language types. MPI also has types corresponding to Fortran data types. We will list only the most common ones for C:
MPI_PACKED and MPI_BYTE are special types and map to any other type. MPI_BYTE expresses an untyped value, and count expresses a number of bytes. It bypasses any data conversion operations in heterogeneous data transfer systems. MPI_PACKED is used with the MPI_PACK procedure, as shown in the ghost data exchange example in Section 8.4.3. You can also define your own data type for use in these calls. This is also demonstrated in the ghost data exchange example. Additionally, there are many procedures for testing data exchange completion, including:
2. Send and receive commands for process-to-process data exchange
int MPI_Test(MPI_Request *request, int *flag, MPI_Status *status)
int MPI_Testany(int count, MPI_Request requests[], int *index, int *flag, MPI_Status *status)
int MPI_Testall(int count, MPI_Request requests[], int *flag, MPI_Status statuses[])
int MPI_Testsome(int incount, MPI_Request requests[], int *outcount, int indices[], MPI_Status statuses[])
int MPI_Wait(MPI_Request *request, MPI_Status *status)
int MPI_Waitany(int count, MPI_Request requests[], int *index, MPI_Status *status)
int MPI_Waitall(int count, MPI_Request requests[], MPI_Status statuses[]) int MPI_Waitsome(int incount, MPI_Request requests[], int *outcount,
int indices[], MPI_Status statuses[])
int MPI_Probe(int source, int tag, MPI_Comm comm, MPI_Status *status)
Additional variants of MPI_Probe are not listed here. The MPI_Waitall procedure is shown in several examples in this chapter. Other procedures are useful in more specialized situations. The names of the procedures give a good indication of the capabilities they provide.
3. Collective Data Exchange: A Powerful Component of MPI
In this section, we'll explore the rich set of collective communication calls in MPI. Collective communication calls operate on a group of processes contained within an MPI communicator. To operate on a partial set of processes, you can create your own MPI communicator for a subset of MPI_COMM_WORLD, like any other process. You can then use your communicator in collective communication calls instead of MPI_COMM_WORLD. Most collective communication routines operate on data. Figure 8.4 provides a visual representation of what each collective operation does.
We'll provide examples of how to use the most common collective operations, as they may be useful in your application. The first example (in Section 8.3.1) demonstrates the use of a barrier. This is the only collective procedure that doesn't operate on data. We'll then show some examples of broadcasting (Section 8.3.2), reducing (Section 8.3.3), and finally, scatter/gather operations (Sections 8.3.4 and 8.3.5). MPI also has many "all-to-all" routines. However, they are expensive and rarely used, so we won't describe them here. All of these collective operations operate on a group of processes, represented by a communication group. All members of the communication group must call the collective, or your program will hang.
3. Collective Data Exchange: A Powerful Component of MPI
Figure 8.4: Data movement within the most common MPI collective routines provides important functionality for parallel programs. The additional MPI_Scatterv, MPI_Gatherv, and MPI_Allgatherv options allow for sending or receiving variable amounts of data between processes. Not shown are some additional routines, such as MPI_Alltoall and similar functions.
3. Collective Data Exchange: A Powerful Component of MPI
Using a Barrier to Synchronize Timers
The simplest collective communication call is MPI_Barrier. It is used to synchronize all processes in the MPI communicator. This isn't necessary in most programs, but it is often used for debugging and synchronizing timers. Let's look at using MPI_Barrier to synchronize timers in the following listing. We also use the MPI_Wtime function to get the current time.
3. Collective Data Exchange: A Powerful Component of MPI
3. Collective Data Exchange: A Powerful Component of MPI
A barrier is inserted before the timer starts and then immediately before the timer stops. This forces the timers of all processes to start at approximately the same time. By inserting a barrier before the timer stops, we obtain the maximum time across all processes. Sometimes using a synchronized timer provides a less confusing measure of time, but in other cases, an unsynchronized timer is better.
NOTE: Synchronized timers and barriers should not be used in production runs; they can seriously slow down the application.
Using Broadcasting to Manipulate Small Input File Data
Broadcasting (or broadcasting) sends data from one processor to all other processors. This operation is shown in Figure 8.4 in the upper left corner. One use of broadcasting, MPI_Bcast, is to send values read from an input file to all other processes. If each process attempts to open the file when there are a large number of processes, it may take several minutes for the file to complete opening. This is because file systems are inherently sequential and are one of the slowest components of a computer system. For these reasons, for a small input file, it is recommended to open and read the file from only one process. The following listing shows how this is done.
3. Collective Data Exchange: A Powerful Component of MPI
3. Collective Data Exchange: A Powerful Component of MPI
3. Collective Data Exchange: A Powerful Component of MPI
It's better to transmit large chunks of data than to transmit many small individual values. Therefore, we broadcast the entire file. To do this, we first broadcast the size so that each process can allocate an input buffer, and then broadcast the data. File reading and broadcasting are performed from rank 0, usually referred to as the master process.
MPI_Bcast takes a pointer for the first argument, so when sending a scalar variable, we send a reference using the ampersand (&) operator to obtain the variable's address. Next comes the count and type to fully define the data being sent. The next argument specifies the originating process. In both of these calls, it is 0 because that is the rank where the data is stored. After this, all other processes in the MPI_COMM_WORLD exchange receive the data. This method is intended for small input files. For input or output of larger files, there are methods for performing parallel file operations. The complex world of parallel input and output is discussed in Chapter 16.
3. Collective Data Exchange: A Powerful Component of MPI
Using Reduce to Extract a Single Value from All Processes
The reduce pattern, discussed in Section 5.7, is one of the most important patterns in parallel computing. The reduce operation is shown in Figure 8.4, top middle. An example of reduce in Fortran syntax for arrays is xsum = sum(x(:)), where the Fortran intrinsic sum sums the array x and places it in the scalar variable xsum. MPI reduce calls take an array or multidimensional array and reduce the values to a scalar result. Many operations can be performed during reduce. The most common are:
The following listing shows the use of MPI_Reduce to obtain the minimum, maximum, and average value of a variable for each process.
3. Collective Data Exchange: A Powerful Component of MPI
3. Collective Data Exchange: A Powerful Component of MPI
3. Collective Data Exchange: A Powerful Component of MPI
The reduction result, in this case the maximum, is stored at rank 0 (argument 6 in the MPI_Reduce call), i.e., in this case, the main process. If we simply wanted to print it in the main process, that would be fine. But if we wanted all processes to be relevant, we would use the MPI_Allreduce procedure.
Alternatively, you can define your own operator. We'll use the example of the enhanced-precision Kahan summation, which we've already worked with and first introduced in Section 5.7. The task in a parallel distributed-memory environment is to perform the Kahan summation over all process ranks. We'll begin by accessing the main program in the following listing, and then move on to the other two parts of the program in Listings 8.13 and 8.14.
3. Collective Data Exchange: A Powerful Component of MPI
3. Collective Data Exchange: A Powerful Component of MPI
The main program (main) shows that the new MPI data type is created once at the beginning of the program and deallocated at the end, before MPI_Finalize. The call to perform a global Kahan summation is executed multiple times in a loop, doubling the data size. Now let's look at the following listing to see what needs to be done to initialize the new data type and operator.
3. Collective Data Exchange: A Powerful Component of MPI
3. Collective Data Exchange: A Powerful Component of MPI
First, we create a new EPSUM_TWO_DOUBLES data type by combining the two base MPI_DOUBLE data types on line 33. We must declare the type outside the procedure on line 19 so that it is available for use by the summation procedure. To create the new operator, we first write a function to use as the operator on lines 22–30. Then, we use esum_type to pass both double values back and forth. We also need to pass the length and data type it will operate on as the new EPSUM_TWO_DOUBLES type.
In creating the Kahan sum reduction operator, we showed you how to create a new MPI data type and a new MPI reduction operator. Now let's move on to actually calculating the global sum of an array over all MPI ranks, as shown in the following listing.
3. Collective Data Exchange: A Powerful Component of MPI
3. Collective Data Exchange: A Powerful Component of MPI
Computing the global Kahan sum is now relatively straightforward. We can perform the local Kahan sum, as shown in Section 5.7. However, we must add MPI_Allreduce on line 52 to obtain the global result. Here, we specify that the allreduce operation will produce a result on all processors, as shown in Figure 8.4 in the upper right corner.
Using the Collect Operation to Organize Debug Prints
The collect operation can be described as a sort operation, where the data from all processors is collected and stored in a single array, as shown in Figure 8.4 in the lower center. This collective communication call can be used to tidy up the console output from your program. By now, you should have noticed that the output printed from multiple ranks of an MPI program is randomly ordered, creating a jumbled, confusing mess. Let's consider a more efficient approach to resolving this situation, ensuring that only the main process's output is retrieved. If we print only the main process's output, the order will be correct. The following listing shows an example program that retrieves data from all processes and prints it in a convenient, organized format.
3. Collective Data Exchange: A Powerful Component of MPI
3. Collective Data Exchange: A Powerful Component of MPI
3. Collective Data Exchange: A Powerful Component of MPI
MPI_Gather accepts a standard triplet describing the data source. We need to use an ampersand to obtain the address of the scalar variable total_time. The destination is also a triplet with the target array times. The array is already an address, so the ampersand is not required. The collection is performed for process 0 of the global communication group MPI_COMM_WORLD. From there, a loop is required to print the time for each process. We prepend a number in the #: format to each string literal to make it clear which process the result applies to.
Using Scatter and Gather to Send Data to Processes for Work
The scatter operation, shown in the lower left corner of Figure 8.4, is the opposite of the gather operation. In this operation, data is transferred from one process to all others in the communication group. The most common use of the scatter operation is in a parallel strategy for distributing data arrays among other processes for work. This is provided by the routines MPI_Scatter and MPI_Scatterv. The following listing shows the implementation.
3. Collective Data Exchange: A Powerful Component of MPI
3. Collective Data Exchange: A Powerful Component of MPI
3. Collective Data Exchange: A Powerful Component of MPI
3. Collective Data Exchange: A Powerful Component of MPI
First, we need to calculate the data size for each process. The desired distribution should be as equal as possible. A simple way to calculate the size is shown in lines 13–15 using simple integer arithmetic. Now we need a global array, but we only need it for the master process. Therefore, we allocate and configure it on this process in lines 18–23. To distribute or collect data, the sizes and offsets for all processes must be known. We see a typical calculation of this data in lines 25–30. The actual scattering is performed using MPI_Scatterv in lines 32–34. The data source is described by the arguments buffer, counts, offsets, and data type. The destination is handled by a standard triple. Then, the rank of the source that will send the data is specified as rank 0. Finally, the last argument is comm, the communication group that will receive the data.
3. Collective Data Exchange: A Powerful Component of MPI
MPI_Gatherv performs the opposite operation, as shown in Figure 8.4. We only need a global array for the main process, so it's allocated only in lines 40–42. The arguments to MPI_Gatherv begin with the source description using a standard triplet. Then, the destination is described using the same four arguments used in scatter. The next argument is the target rank, followed by the communication group.
It should be noted that all sizes and offsets used in the MPI_Gatherv call are integer types. This limits the size of the data being manipulated. An attempt was made to change the data type to long to accommodate larger data in version 3 of the MPI standard. This was not adopted because it would have broken too many applications. Watch for new calls that will support the long integer type in a future MPI standard.
4. Data Parallelism Examples
The data parallelism strategy defined in Section 1.5 is the most common approach in parallel applications. In this section, we'll look at several examples of this approach. First, we'll consider a simple case of a stream triad, where no data sharing is required. Then, we'll look at more typical ghost cell data sharing techniques used to link subdivided domains allocated to each process.
Stream Triad for Measuring Throughput on a Node
STREAM Triad is the source code for the bandwidth benchmark presented in Section 3.2.4. This version uses MPI to allow more processes to run on a node, and potentially across multiple nodes. The goal of increasing the number of processes is to see the maximum throughput of a node while utilizing all processors. This yields a target throughput for more complex applications. As shown in Listing 8.17, the source code is simple because no data sharing is required between ranks. Timing is tracked only on the master process. You can run it first on a single processor, and then across all processors on your node. Do you get the full parallel acceleration you'd expect from increasing the number of processors? To what extent does system memory bandwidth limit your acceleration?
4. Data Parallelism Examples
4. Data Parallelism Examples
4. Data Parallelism Examples
Ghost Cell Communication in a 2D Computational Grid
Ghost cells are a mechanism we use to link computational grids on adjacent processors. They are used to cache values from adjacent processors, reducing the need for data exchanges. The ghost cell technique is the single most important method for achieving distributed memory parallelism in MPI.
Let's briefly discuss the terms "halo" and "ghost cell“. Even before the era of parallel processing, a region of cells around a grid was often used to implement boundary conditions. These boundary conditions could be reflective, incoming, outgoing, or periodic. To improve efficiency, programmers wanted to avoid if statements in the main computation loop. To achieve this, they added cells around the grid and set their appropriate values before the main computation loop. These cells resembled a halo, hence the name. Halo cells are any set of cells surrounding a computational grid, regardless of their purpose. Then the domain boundary halo is the halo cells used to impose a specific set of boundary conditions.
4. Data Parallelism Examples
After parallelizing applications, a similar outer region of cells was added to store values from adjacent cells. These cells are not real cells, but exist only to reduce data exchange costs. Because they are not real cells, they soon became known as "ghost cells." The real data in ghost cells resides on an adjacent processor, and the local copy is simply a ghost value. Ghost cells also look like halos and are also called halo cells. Updates or data exchanges in ghost cells are referred to as ghost cell updates and are only necessary for parallel multiprocess runs when you need updates to real values from adjacent processes.
Boundary conditions must be satisfied for both serial and parallel runs. Confusion arises because these operations are often referred to as halo updates, although it is unclear what exactly is meant. In our terminology, halo updates refer to both domain boundary updates and ghost cell updates. To optimize data exchange within MPI, we only need to look at updates or phantom cell exchanges and put boundary condition calculations aside for now.
4. Data Parallelism Examples
Let's now consider setting up ghost cells for the local grid boundaries in each process and executing communications between subdomains. Using ghost cells, the necessary messages are grouped into fewer messages than if one message were executed each time a cell value was needed from another process. This method is the most common and enables efficient use of data parallelism. In the implementations of ghost cell updates, we demonstrate the use of the MPI_Pack procedure and load the communication buffer using a simple cell-by-cell assignment in an array. In subsequent sections, we will also address executing the same communications using MPI data types, using MPI topology calls for setup and data exchange.
Implementing ghost cell updates in the data parallelism source code handles most of the necessary communications. This isolates the source code that provides parallelism into a small section of the application. This small section of code is important for optimizing for parallel efficiency. Let's look at some implementations of this functionality, starting with the setup in Listing 8.18 and the work performed by stencil loops in Listing 8.19. You may want to view the full source code for this chapter's example code in the GhostExchange/GhostExchange_Pack directory at https://github.com/EssentialsOfParallelComputing/Chapter8.
4. Data Parallelism Examples
4. Data Parallelism Examples
We allocate memory for the local size plus space for halos in each process. To make indexing a little easier, we shift the memory indexing so it starts at -nhalo and ends at isize+nhalo. This way, the actual cells will always be between 0 and isize-1, regardless of the halo width.
The following lines show a special malloc2D call with two additional arguments that shift the array addressing so that the actual portion of the array is between 0,0 and jsize,isize. This is done using some pointer arithmetic that moves the starting location of each pointer.
64 double** x = malloc2D(jsize+2*nhalo, isize+2*nhalo, nhalo, nhalo);
65 double** xnew = malloc2D(jsize+2*nhalo, isize+2*nhalo, nhalo, nhalo);
To provide the work, we use a simple stencil calculation from the blur operator, shown in Figure 1.10. Many applications have much more complex calculations that take much longer. The following listing shows the stencil calculation loops.
4. Data Parallelism Examples
4. Data Parallelism Examples
Now we can look at the critical source code for updating ghost cells. Figure 8.5 shows the required operation. The width of a ghost cell region can be one, two, or more cells deep. Some applications may also require corner cells. Each of the four processes (or ranks) needs data from that rank: left, right, top, and bottom. Each of these processes requires a separate data exchange and a separate data buffer. The width of the ghost region depends on the application and the need for corner cells.
Figure 8.5 shows an example of ghost cell exchange for a 4x4 grid on nine processes with a one-cell-wide halo and corners. The outer boundary halos are updated first, followed by horizontal data exchange, synchronization, and vertical data exchange. If corners are not needed, horizontal and vertical exchanges can be performed simultaneously. If corners are needed, synchronization between the horizontal and vertical exchanges is required.
A key observation when updating ghost cell data is that in C, row-wise data is continuous, while column-wise data is separated by increments equal to the row size. Sending individual values for columns is expensive, so we need to somehow group them together.
4. Data Parallelism Examples
There are several ways to update ghost cells using MPI. In this first version, Listing 8.20, we'll consider an implementation that uses the MPI_Pack call to pack column-wise data. Row-wise data is sent only using the standard MPI_Isend call. The width of the ghost cell region is specified by the nhalo variable, and angles can be requested as input, as appropriate.
4. Data Parallelism Examples
Figure 8.5. In the corner cell version of ghost cell updating, left and right data are exchanged first (in the top half of the figure), followed by top and bottom data (in the bottom half of the figure). With care, the left and right exchange can be reduced to just real cells plus cells along the outer boundary, although there's no harm in doing it across the entire vertical dimension of the grid. Updating the boundary cells around the grid is done separately.
4. Data Parallelism Examples
4. Data Parallelism Examples
4. Data Parallelism Examples
4. Data Parallelism Examples
Calling MPI_Pack is especially useful when multiple data types need to be transferred during ghost updates. Values are packed into a type-independent buffer and then unpacked from the other side. Exchanges with neighbors in the vertical direction are performed using continuous row data. When corners are taken into account, a single buffer works well. Without corners, individual halo rows are sent. Typically, there are only one or two halo cells, so this approach is reasonable.
Another way to load buffers for data exchange is with array assignments. This approach works well when there is a single simple data type, such as the double-precision floating-point type used in this example. The following listing shows the source code for replacing MPI_Pack loops with array assignments.
4. Data Parallelism Examples
4. Data Parallelism Examples
4. Data Parallelism Examples
The MPI_Irecv and MPI_Isend calls now use a counter (count) and the MPI_DOUBLE data type, rather than the generic byte type of the MPI_Pack function. We also need to know the data type for copying data to and from the communication buffer.
Communicating with Ghost Cells in 3D Stencil Calculation
Communicating with ghost cells is also possible for 3D stencil calculation. We do this in Listing 8.22. However, the organization of this process is slightly more complex. First, the process map is calculated as xcoord, ycoord, and zcoord values. Then, neighbors are determined and the data sizes on each processor are calculated.
4. Data Parallelism Examples
Updating ghost cells, including copying arrays to buffers, swapping, and copying, takes a couple hundred lines and is impossible to demonstrate here. For a detailed implementation, please refer to the source code examples (https://github.com/EssentialsofParallelComputing/Chapter8) included with this chapter. We will demonstrate a version of updating ghost cells with the MPI data type in Section 8.5.1.
5. Advanced MPI functionality to simplify source code and provide optimizations
MPI's superior design becomes evident when we see how MPI's core components can be combined to create higher-level functionality. We saw this in Section 8.3.3, when we created the new double-double type and the new reduce operator. This extensibility gives MPI important capabilities. We'll look at a couple of these advanced features, which are useful in common data parallelism applications. These include:
5. Advanced MPI functionality to simplify source code and provide optimizations
Using MPI Application-Specific Data Types to Improve Performance and Simplify Code
MPI has a rich set of functions for creating new MPI application-specific data types from MPI base types. This allows complex data to be encapsulated into a single application-specific data type that can be used in communication calls. As a result, a single communication call can send or receive many smaller chunks of data as a single unit. Here is a list of some of the MPI data type creation functions:
Below is a visual illustration to help you understand some of these data types. Figure 8.6 shows several simpler and more commonly used functions, including MPI_Type_contiguous, MPI_Type_vector, and MPI_Type_create_subarray.
5. Advanced MPI functionality to simplify source code and provide optimizations
Figure 8.6 Three application-specific MPI data types with illustrations of the arguments used to create them
5. Advanced MPI functionality to simplify source code and provide optimizations
Once a data type is declared and converted to a new data type, it must be initialized before it can be used. Several additional routines exist for committing and freeing types for this purpose. A type must be committed before it can be used, and it must be freed to avoid memory leaks. These routines are:
We can greatly simplify communication with ghost cells by defining an MPI application-specific data type, as shown in Figure 8.6, to represent a column of data and avoid calls to MPI_Pack. By defining an MPI data type, we can avoid the extra data copy. Data can be copied from its normal location directly into the MPI send buffers. Let's see how this is done in Listing 8.23. Listing 8.24 shows the second part of the program.
First, we set up application-specific data types. We use MPI_Type_vector for array access sets with index stride. In the case of continuous data for the vertical type, when we use angles, we use MPI_Type_contiguous, and in lines 139 and 140, we free the data type at the end before MPI_Finalize.
5. Advanced MPI functionality to simplify source code and provide optimizations
5. Advanced MPI functionality to simplify source code and provide optimizations
You can then write ghostcell_update more concisely and with better performance using MPI data types, as shown in the following listing. If we need to update the angles, synchronization between the two data exchange passes is necessary.
5. Advanced MPI functionality to simplify source code and provide optimizations
5. Advanced MPI functionality to simplify source code and provide optimizations
Higher performance is often cited as the reason for using MPI data types. This does allow MPI implementations to be implemented, in some cases avoiding an additional copy. However, in our view, the most important reason for using MPI data types is cleaner, simpler source code and fewer opportunities for errors.
The three-dimensional version using MPI data types is a bit more complex. In the following listing, we use MPI_Type_create_subarray to create three application-specific MPI data types that will be used for data exchange.
5. Advanced MPI functionality to simplify source code and provide optimizations
5. Advanced MPI functionality to simplify source code and provide optimizations
The following listing shows that the data exchange procedure using these three MPI data types is quite concise.
5. Advanced MPI functionality to simplify source code and provide optimizations
5. Advanced MPI functionality to simplify source code and provide optimizations
5. Advanced MPI functionality to simplify source code and provide optimizations
5. Advanced MPI functionality to simplify source code and provide optimizations
Cartesian Topology Support in MPI
In this section, we'll demonstrate how topology functions work in MPI. The operation is still a ghost cell exchange, as shown in Figure 8.5, but we can simplify the coding using Cartesian functions. General graph functions for unstructured applications are not covered. We'll start with setup routines and then move on to communication routines.
In the setup routines, you need to specify the values for assignments in the process grid, and then specify the neighbors, as was done in Listings 8.18 and 8.22. As shown in Listing 8.24 for two dimensions and Listing 8.25 for three dimensions, a process sets the dims array to the number of processes used in each dimension. If any of the values in the dims array is zero, the MPI_Dims_create function calculates some values that will be used. Note that the number of processes in each direction doesn't take into account the grid size and may not yield good results for long, narrow tasks. Consider the case of an 8x8x1000 grid and assign it 8 processors; the process grid will be 2x2x2, resulting in a 4x4x500 grid domain for each process.
5. Advanced MPI functionality to simplify source code and provide optimizations
The MPI_Cart_create procedure takes the resulting array, dims, and the input array, periodic, which declares whether the boundary will cycle to the opposite side and vice versa. The last argument is the reordering argument, which allows MPI to change the order of processes. In this example, it is zero (false). Now we have a new communicator containing topology information. Obtaining a process mesh diagram is a simple call to MPI_Cart_coords. Neighbors can be obtained by calling MPI_Cart_shift, where the second argument specifies the direction, and the third argument is the offset or number of processes in that direction. The output is the ranks of adjacent processors.
5. Advanced MPI functionality to simplify source code and provide optimizations
The 3D Cartesian topology setup is similar but has three dimensions as shown in the following listing.
5. Advanced MPI functionality to simplify source code and provide optimizations
5. Advanced MPI functionality to simplify source code and provide optimizations
Comparing this source code with the versions in Listings 8.19 and 8.23, we see that the topology functions in this relatively simple setup example don't save many lines of code or significantly reduce programming complexity. We can also effectively use the Cartesian communicator created in line 70 of Listing 8.28 for neighbor communication. This is where the greatest reduction in lines of code is observed. The MPI function has the following arguments:
int MPI_Neighbor_alltoallw(const void *sendbuf,
const int sendcounts[], const MPI_Aint sdispls[],
const MPI_Datatype sendtypes[],�
void *recvbuf,
const int recvcounts[], const MPI_Aint rdispls[],
const MPI_Datatype recvtypes[], MPI_Comm comm)
5. Advanced MPI functionality to simplify source code and provide optimizations
A neighbor call has many arguments, but once we've set it all up, the exchange is concise and performed in a single statement. We'll go over all the arguments in detail because they can be difficult to understand.
A neighbor communication call can either use a filled send and receive buffer or perform the operation in-place. We'll demonstrate the in-place method. The send and receive buffers are a two-dimensional array. We'll use the MPI data type to describe the data block, so the counters will be an array with a value of one for all four Cartesian sides in the case of two dimensions, or six sides in the case of three dimensions. The data exchange order for the sides is: bottom, top, left, right for two dimensions, and front, back, bottom, top, left, right for three dimensions, and is the same for both the send and receive types.
The data block is different for each direction: horizontal, vertical, and depth. We use the drawings in the generally accepted standard orientation, where x goes to the right, y goes up, and z (depth) goes back into the page. However, the data block is the same in each direction, but with different offsets to the beginning of the data block. The offsets are specified in bytes, so you'll see the offsets multiplied by 8, the size of a double-precision data type. Now let's look at how all this fits into the source code for setting up data exchange for the two-dimensional case in the following listing.
5. Advanced MPI functionality to simplify source code and provide optimizations
5. Advanced MPI functionality to simplify source code and provide optimizations
5. Advanced MPI functionality to simplify source code and provide optimizations
The configuration for three-dimensional Cartesian neighbor communication uses the MPI data types from Listing 8.25. The data types define the block of data to be transferred, but we need to define the byte offset to the starting location of the block of data for sending and receiving. We also need to define arrays for the sendtype and recvtype types in the proper order, as in the following listing.
5. Advanced MPI functionality to simplify source code and provide optimizations
5. Advanced MPI functionality to simplify source code and provide optimizations
The actual exchange is accomplished with a single call to MPI_Neighbor_alltoallw, as shown in Listing 8.31. There's also a second block of code for corner cases, which requires multiple calls with synchronization between them to ensure the corners are filled properly. The first call executes only the horizontal direction and then waits for completion before executing the vertical direction.
5. Advanced MPI functionality to simplify source code and provide optimizations
Three-dimensional Cartesian swapping with neighbors is similar, but with the addition of the z-coordinate (depth). Depth occupies the first position in the count and type arrays. In the phase-based swap for corners, depth is determined after horizontal and vertical swapping with ghost cells, as shown in the listing below.
5. Advanced MPI functionality to simplify source code and provide optimizations
5. Advanced MPI functionality to simplify source code and provide optimizations
5. Advanced MPI functionality to simplify source code and provide optimizations
Performance Tests for Ghost Cell Exchange Variations
Let's test these ghost cell exchange variants on a test system. We'll use two Broadwell nodes (Intel® Xeon® E5-2695 v4 processors running at 2.10 GHz) with 72 virtual cores each. We could run this test on more compute nodes with different MPI library implementations, halo sizes, cell sizes, and higher-performance interconnects to better understand how each ghost cell exchange variant performs. The source code is below:
mpirun -n 144 --bind-to hwthread ./GhostExchange -x 12 -y 12 -i 20000 \
-j 20000 -h 2 -t -c
mpirun -n 144 --bind-to hwthread ./GhostExchange -x 6 -y 4 -z 6 -i 700 \
-j 700 -k 700 -h 2 -t -c
5. Advanced MPI functionality to simplify source code and provide optimizations
The GhostExchange program options are:
-x (processes in the x-direction);
-y (processes in the y-direction);
-z (processes in the z-direction);
-i (grid size in the i- or x-direction);
-j (grid size in the j- or y-direction);
-k (grid size in the k- or z-direction);
-h (halo cell width, typically 1 or 2);
-c (include corner cells).
Exercise: Ghost Cell Tests
The included source code is configured to run the full range of ghost cell exchange cases. The batch.sh file can be used to change the halo size and whether or not to include corner cells. The file is configured to run all test cases 11 times on two Skylake Gold nodes with a total of 144 processes.
5. Advanced MPI functionality to simplify source code and provide optimizations
cd GhostExchange
./build.sh
./batch.sh |& tee results.txt
./get_stats.sh > stats.out
You can then generate graphs using the provided scripts. For Python plotting scripts, you need the matplotlib library. Results are shown for average execution times.
python plottimebytype.py
python plottimeby3Dtype.py
The following figure shows graphs for small test cases. The MPI data type versions are slightly faster even at this small scale, suggesting that data copying may be worth avoiding. For MPI Cartesian topology calls and MPI data types, the bigger benefit is that the ghost communication code is significantly simplified. Using these more advanced MPI calls requires additional setup, but this is only done once per run.
5. Advanced MPI functionality to simplify source code and provide optimizations
Relative performance of 2D and 3D ghost exchanges on two nodes with a total of 144 processes. MPI data types in MPI and CNeighbor are slightly faster than a buffer explicitly filled with array assignment loops. In 2D ghost exchanges, packing procedures are slower, although in this case, an explicitly filled buffer is faster than MPI data types.
6. Hybrid MPI plus OpenMP technique for maximum scalability
A combination of two or more parallelization techniques is called hybrid parallelization, in contrast to all MPI implementations, which are also called pure MPI or MPI-everywhere. In this section, we consider a hybrid MPI plus OpenMP technique, where MPI and OpenMP are used together in an application. This typically involves replacing some MPI ranks with OpenMP threads. For larger parallel applications spanning thousands of processes, replacing MPI ranks with OpenMP threads potentially reduces the overall size of the MPI domain and the memory required at extreme scale. However, the added performance of a thread-level parallelism layer may not always justify the added complexity and development time. For this reason, hybrid MPI plus OpenMP implementations are typically the domain of extreme applications in both size and performance.
Benefits of a Hybrid MPI and OpenMP Technique
When performance becomes critical enough to justify the added complexity of hybrid parallelism, there may be several advantages to adding an OpenMP parallel layer to MPI-based code. For example, these benefits may include:
6. Hybrid MPI plus OpenMP technique for maximum scalability
Spatially decomposed parallel applications that utilize subdomains with ghost cells (halos) will have fewer total ghost cells per node when adding thread-level parallelism. This reduces both memory requirements and communication costs, particularly on multi-core architectures such as Intel Knights Landing (KNL). Using shared-memory parallelism can also improve performance by reducing contention for the network interface card (NIC), avoiding unnecessary data copying used by MPI for per-node communications. Furthermore, many MPI algorithms are tree-based and scale as log2n. Reducing the execution time by 2n threads reduces the tree depth and gradually improves performance. Although the remaining work must still be performed by threads, this impacts performance by reducing synchronization costs and communication latency. Threads can also be used to improve load balancing within a NUMA region or compute node.
In some cases, a hybrid parallel approach is not only advantageous but also necessary to access the full performance potential of the hardware. For example, some hardware, and perhaps memory controller functionality, is only accessible to threads and not processes (MPI ranks). These issues were observed on the multi-core Intel Knights Corner and Knights Landing architectures. In MPI + X + Y, where X is threading and Y is GPU parlance, we often correlate ranks with the number of GPU processors. OpenMP allows an application to continue accessing other processors for CPU execution. Other solutions exist for this purpose, such as MPI_COMM groups and shared memory functionality for MPI, or simple GPU management from multiple MPI ranks.
6. Hybrid MPI plus OpenMP technique for maximum scalability
Thus, while executing source code using MPI-Everywhere technology is attractive on modern multi-core systems, there are concerns about scalability as the number of cores increases. If you are looking for extreme scalability, you will need an efficient OpenMP implementation in your application. We discussed our much more efficient high-level OpenMP design in the previous chapter, in sections 7.2.2 and 7.6.
An Example of MPI and OpenMP
The first steps to implementing a hybrid MPI and OpenMP technique are to provide information about what you will be doing. This is done in the call to MPI_Init at the very beginning of the program. You should replace the call to MPI_Init with a call to MPI_Init_thread as follows:
MPI_Init_thread(&argc, &argv, int thread_model required, int *thread_model_provided);
6. Hybrid MPI plus OpenMP technique for maximum scalability
The MPI standard defines four thread models. These models provide different levels of thread safety for MPI calls. In order of increasing thread safety:
Many applications perform communication at the main loop level, and OpenMP threads are used for key computational cycles. For this pattern, MPI_THREAD_FUNNELED works perfectly.
NOTE: It's best to use the lowest thread safety level you need. Each higher level imposes a performance penalty, since the MPI library must place mutexes or critical blocks around the send and receive queues and other core parts of MPI.
6. Hybrid MPI plus OpenMP technique for maximum scalability
Now let's look at the changes required in our stencil example to add OpenMP threading. For this exercise, we'll be modifying the CartExchange_Neighbor example. The following listing shows that the first change is a modification to the MPI initialization.
6. Hybrid MPI plus OpenMP technique for maximum scalability
A mandatory change is the use of MPI_Init_thread instead of MPI_Init on line 27. The additional code checks whether the requested thread safety level is available and terminates if it is not. We also print the number of threads in the main rank zero thread. Now let's move on to the changes to the computation loop, shown in the following listing.
6. Hybrid MPI plus OpenMP technique for maximum scalability
The changes required to add OpenMP threading consist of adding one pragma at line 157. As a bonus, we'll show how to add vectorization to the inner loop with another pragma inserted at line 159.
Now you can try running this MPI/OpenMP/Vectorization hybrid example on your system. However, to achieve good performance, you'll need to control the placement of MPI ranks and OpenMP threads. This is done by setting affinity—we'll cover this topic in more detail in Chapter 14.
DEFINITION: Affinity favors a particular hardware component when scheduling a process, rank, or thread. It's also known as pinning or affinity.
Tuning the affinity of your ranks and threads becomes more important as the node complexity increases and in hybrid parallel applications. In the previous examples, we used
--bind-to core and --bind-to hwthread to improve performance and reduce runtime performance variability caused by rank migration from one core to another. In OpenMP, we used environment variables to specify placements and affinities. An example is:
export OMP_PLACES=cores
export OMP_CPU_BIND=true
6. Hybrid MPI plus OpenMP technique for maximum scalability
For now, let's start by pinning MPI ranks in sockets so that threads can be distributed to other cores, as we demonstrated in our ghost cell benchmark example for the Skylake Gold process. Here's how it's done:
export OMP_NUM_THREADS=22
mpirun -n 4 --bind-to socket ./CartExchange -x 2 -y 2 -i 20000 -j 20000 \
-h 2 -t –c
We execute four MPI ranks, each spawning 22 threads, as specified by the OMP_NUM_THREADS environment variable, for a total of 88 processes. The --bind-to socket option for mpirun instructs it to bind processes to the socket they are hosted on.
7. Exercises
1. Why can't we simply block on receives, as we did for send/receive in the ghost exchange, using the packing buffer or array methods, respectively, in Listings 8.20 and 8.21?
2. Is it safe to block on receive, as shown in Listing 8.8 in the vector-type version of the ghost exchange? What are the advantages of only blocking on receive?
3. Modify the vector-type ghost cell exchange example in Listing 8.21 to use receive locking instead of wait. Will this be faster? Will it always work?
4. Try replacing explicit tags in one of the ghost exchange routines with MPI_ANY_TAG. Will this work? Will it be even slightly faster? What advantage do you see in using explicit tags?
5. Remove barriers for synchronized timers in one of the ghost exchange examples. Run the source code with the original synchronized timers and unsynchronized timers.
6. Add the timer statistics from Listing 8.11 to the thread triad throughput measurement source code in Listing 8.17.
7. Apply the steps for converting high-level OpenMP to the hybrid MPI plus OpenMP example in the source code accompanying this chapter (HybridMPIPlusOpenMP directory). Experiment with vectorization, number of threads, and MPI ranks on your platform.