Step 1: Identify Computational Hotspots

We begin any parallelization process by identifying the parts of the program that consume the most run time.

The goal is to know which code should be parallelized and which code should be recycled from the serial program.

The downside of parallel programming:

• Why not parallelize everything?

• Why do we care which part of our program should be parallelized and which should be recycled?

• The parallelization process:
  • Is time-consuming.
  • Increases the complexity of your program.
  • Reduces the portability of your program.
  • Introduces many opportunities for user and system bugs.
  • Makes your program harder to modify, or reuse, at a later time.

The reason is obvious:

• Greatly optimizing a procedure which only takes 10 percent of the total run time can at most improve your performance by about 10 percent.

• Greatly optimizing a procedure which takes 90 percent of the total run time can improve your performance by a factor of 10.

• Abstract analysis of the algorithm

• Visual inspection of the code

• Profiling tools. Many tools exist which can profile run time behavior and quickly identify hotspots. See the Performance Analysis Tools tutorial for a description of a number of such tools.

In addition, the complexity and feasibility of the parallelization process should be considered when evaluating parallelization targets.

Step 1: Hotspot Identification Checklist

• Are you satisfied you have identified a code fragment (your parallelization target) that consumes a large enough percentage of your programs run time?

• Are you satisfied that your parallelization target does not include extensive amounts of code that would not significantly improve your performance?

• Is your parallelization target parallelizable?

Step 1: Hotspot Identification Example
Generic File I/O

Many programs begin with several pages of code that read input data from a file, create complex data structures, write this data into those data structures, and perform a small amount of preprocessing on that data.

Similarly, many programs perform complimentary tasks at the end which postprocess the data and output results to files.

Observation: These processes are usually not good candidates for parallelization:

• While such input and output files can be quite large, their file I/O time typically consumes a small amount of run time relative to the rest of the code.

• Such files typically reside on a single disk; any attempts to parallelize the file I/O might ultimately be serialized at the disk.

Step 2: Partitioning

In this step we partition our problem into smaller tasks that can execute in parallel with each other.

The intent is to expose multiple opportunities for parallel execution rather then to suggest the actual partition of tasks that will be executed in parallel.

Later, considerations of communication, memory, and task requirements will lead us to agglomerate smaller tasks into larger tasks.

For now, the goal is to partition the tasks as finely as possible so as to provide greater flexibility to make choices during later steps.

There are two primary methods for partitioning a problem.

Data Decomposition: partition the data first. Then partition the computation based on the data partition. Some general heuristics:

• Begin by focusing on the largest data structures, or the ones that are accessed most frequently.

• Divide the data into small pieces.

• If possible, divide the data into same size units.

• Strive for more aggressive partitioning then your target computer will allow.

• Use the data partitioning as a guideline for partitioning the computation into separate tasks; associate some computation with each data element.

• Take communication requirements into some account when partitioning data.

• Partition the computation into many small tasks.

• Strive for that tasks with uniform computational requirements.

• Associate data with each task.

Alternatively, combine the two strategies at different points in your problem.

• Some problems focus on different data structures at different times. Consider different partitionings for these different phases.

• Problems where no single data or functional partitioning yields satisfactory results.

Step 2: Partitioning Checklist

There are exponentially many ways to partition a problem; some may be significantly better or worse than others. The following checkpoints should help you determine if your partitioning is appropriate:

1. Maximize flexibility; does your partition define at least an order of magnitude more tasks than there are processors?

2. Scaling to larger problem sizes: does your partition avoid redundant storage requirements?

3. Load balancing: Are tasks of comparable size?

4. Problem scalability: does the number of tasks scale with problem size?

5. Maximize flexibility; have you identified several alternative partitions?

Step 2: Partitioning Example
Generic 3-D Grid

Many programs focus their computation on a 3-dimensional grid structure. Depending on the nature of the program, computation could be focused on planes, columns, individual points, or blocks of data in any combination of the 3 dimensions. Some example data partitions would be:

• The 1-D decomposition might be appropriate when you know that computation is focused on planes of the grid with much interaction between elements in the same plane, but little interaction between elements in different planes.

• Similarly for the 2-D decomposition.

• A 3-D decomposition should be used if interaction between each element of the grid is not known to favor elements in any particular dimension.

Step 2: Partitioning Example
Matrix Multiply

In the common matrix multiply problem, two 2-d matrices, A and B, are multiplied to form a third 2-d matrix C. If all matrices are NxN, this problem requires about NxNxN floating-point operations. Each element in C is computed from the dot product of a row of A with a column of B.

   for(i=0;i < NRA;i++)
      for(j=0;j < NCB;j++)
         for(k=0;k < NCA;k++)
            c[i][j]+= a[i][k] * b[k][j];

• The A matrix is partitioned by rows; since rows of A are accessed together, this partition would minimize communication during dot product calculation.

• The B matrix is partitioned by columns; see above.

• Functional partitioning on each element of C.

Some potentially bad partitions:

• The A matrix is partitioned by columns; since rows of A are accessed together, all dot produces will involve communication with each other.

• The B matrix is partitioned by rows: see above.

Step 2: Partitioning Example
Mandelbrot Image Calculation

In a Mandelbrot image calculation, a 2-d image is generated in which each pixel's value is calculated solely as a function of its coordinates.

for (i=0;i < IMAGE_HEIGHT;i++)
{
   p_imag = ystart + ((i + (taskid * (IMAGE_HEIGHT/numtasks))) * yinc);
   for (j=0;j < IMAGE_WIDTH;j++)
   {
      p_real = xstart + (j * xinc); z_real = p_real; z_imag = p_imag;
      count = 0;
      t_real = z_real * z_real; t_imag = z_imag * z_imag; z_size = t_real + t_imag;
      while ((count < MAX_COUNT) && (z_size < MAX_SIZE))
      {
         z_imag = ((z_real + z_real) * z_imag) + p_imag;
         z_real = (t_real - t_imag) + p_real; t_real = z_real * z_real;
         t_imag = z_imag * z_imag; z_size = t_real + t_imag;
          count++;
      }
      output_image[i][j] = count * 2;
   }
}

Some example partitions would be:

• Partition by individual output pixels; most aggressive partitioning.

• Partition by rows.

• Partition by columns.

• Partition by 2-D blocks.

Step 3: Communication

Ideally, the tasks defined in the previous step could execute completely in parallel, but in most real world programs these tasks require data computed or used by other tasks; communication is necessary.

Data communication is one of the main obstacles to efficient parallel programming.

In the communication step, we identify the data communication requirements of your partition by drawing lines between tasks that either produce or consume data from other tasks.

The intermediate goal is to identify tasks which most frequently communicate with another. In the next step, those tasks will be agglomerated into larger tasks.

Another goal is to identify those data structures which are needed by many tasks and to make decisions on whether to replicate them.

The process: draw lines between those tasks that must communicate.

A small number of communication patterns typically arise:

• Local Communication: tasks communicate with a small number of local neighbor tasks.

• Global Communication: tasks communicate with potentially many, non-local tasks.
  • Broadcast: a data structure must be broadcast from one task to all other tasks. Consider computing the data structure locally at every task to avoid the broadcast: SP2 tip.

  • Reduction: an associative reduction operation must be performed on all data elements producing a single result: SP2 tip.

  • Unstructured: no identifiable pattern.

• Static Communication: communication links between tasks are fixed at compile- or load-time.

• Dynamic Communication: communication links between tasks are unknown at compile- or load-time.

Step 3: Communication Checklist

The communication step requires few decisions.

The purpose is to make you understand the relationships between different data structures and task in your program.

This will allow you to better map those data structures and tasks to processors.

Step 3: Communication Example
Jacobi Finite Difference Method

The communication patterns of the Jacobi finite difference method are a good example of an algorithm with local communication requirements. In this example, a two-dimensional grid is repeatedly updated by assigning a function of neighbor values to each element in a 2-d grid X:

The communication structure for a single element in such a grid would look like this:

Step 3: Communication Example
Global Reduction Operations

Many programs contain phases where a single value must be summed among each task. The communication structure of such an operation would look like this:

This type of communication operation is inefficient because all receive operations, and all additions, must be serialized.

A better communication strategy would be to use the divide-and-conquer strategy to perform the global sum with a binary tree:

The running time of this code fragment would be reduced from O(N) to O(LogN).

Step 3: Communication Example
Affine Image Warping

In Affine image warping with LaGrange interpolation, an input image is rotated, translated and scaled to produce an output image; each pixel in the output image is overwritten by some pixel in the input image. A small (4 by 64 floating-point numbers) interpolation table is used to interpolate input image values that fall in between pixels.

A code fragment is available here.

Regardless of our data/task partitioning scheme, each task will need to access unknown entries of the interpolation table.

If that table were stored at a single task, excessive communication would occur.

Observation: the interpolation table is computed only once.

Suggestion: replicating the table at each processor will reduce total communication costs and not severely consume task memory.

 

Step 4: Agglomeration

At this point we have defined a data partition and identified many tasks which can execute in parallel along with their communication requirements.

In the agglomeration step, we group many small tasks into larger tasks in such a way that we:

• reduce total communication (by grouping tasks that must communicate together).

• increase message-size granularity (by grouping pairs of tasks that frequently communicate together).

• This would increase total communication by failing to assign frequently communicating tasks to the same processor.

• This would sacrifice the improved communication performance inherent in large-grained messages.

• This would sacrifice communication locality.

How much should we agglomerate? Assume we start with T tasks that will run on P processors.

• P Agglomerated Tasks: If the task run times are knowable at compile time, wecan create fully-agglomerated tasks.

• ~10P Agglomerated Tasks: If tasks run times are unknowable, or their total number is unknown at load time, an alternative is to create semi-agglomerated tasks with an order of magnitude more tasks than processors and dynamically schedule them at run time.
  • Staticly combining the T smaller tasks into P larger tasks could result in gross load imbalances.

  • Dynamically created tasks must be scheduled at run time.

• T Non-Agglomerated Tasks: Generally this is a bad idea. Dynamic Load-balancing techniques impose a fixed overhead per task; less tasks reduce this overhead.

Step 4: Agglomeration:
How To Agglomerate

Whether we are creating fully- or semi-agglomerated tasks, the basic guidelines are the same (and an exercise in satisfying conflicting demands).

It's often useful to know the communication performance capabilities of your target machine while agglomerating tasks: SP2 Communication Performance Summary.

1. Tasks that can run concurrently should be placed in different groups.

2. Tasks that require the same input data should be grouped together.

3. Tasks that communicate with each other should be grouped together.

4. Pairs of tasks that frequently communicate should be grouped together (to increase message granularity).

5. Surface-to-Volume issues: for problems that only communicate with local neighbors, we can simultaneously increases the message granularity, decrease the total communication requirements, and increase the computation granularity by agglomerating neighbors.
   • Illustration
   • Mathematical explanation

Step 4: Agglomeration:
Special Case

Special Case: agglomerating T tasks with varying run times to P tasks. We must employ a static load balancing scheme to ensure processors have equal work load.

• Recursive Bisection

  • Used to partition unstructured data domains so that computational load per processor is even and communication between processors is minimized.

  • The divide-and-conquer approach is taken in agglomerating the domain; the entire domain is first divided in half, then the halves are divided in half, etc.

  • An illustration of a recursively bisected 2-D grid is available here.

• Local Algorithms
  • An initial agglomeration/mapping is specified.

  • Periodically, processors with neighbor data exchange their computational load and exchange data/tasks if an imbalance exists.

• Probabalistic Methods:
  • Tasks are mapped to processors at random.

  • If enough tasks exist, an even computational load per processor can be expected.

  • Cyclic mappings are a form of a probabalistic method.

Step 4: Agglomeration Checklist

1. Has agglomeration reduced total communication costs by grouping tightly coupled tasks together?

2. Has agglomeration yielded tasks with similar computation and communication costs? If we have created just one task per processor, then these tasks should have nearly identical costs.

3. Does the number of tasks still scale with problem size?

4. Has agglomeration preserved opportunities for concurrent execution?

5. If you expect to use a dynamic task scheduler, does your agglomeration have an order of magnitude more tasks than processors?

Step 4: Agglomeration Example
2-D Jacobi Finite Difference Method

We examine the 2-D Jacobi finite difference method again as an example for task agglomeration. Recall that each element in the 2-d grid is updated each iteration as a function of its 4 neighbors. By grouping together grid elements in 2-D blocks, we can reduce total communication, increase communication granularity and increase computational granularity.

Step 4: Agglomeration Example
Mandelbrot Image Calculation  

Recall that in Mandelbrot image calculation, a 2-d image is generated in which each pixel's value is calculated solely as a function of its coordinates. Example image

In Step 2 we suggested that partitioning by pixels would be most aggressive.

We now consider the agglomeration phase. Inspection of the code reveals that the computation time is not constant for each pixel. In fact, different pixel colors in the output image correspond to different amounts of computation.

Because computation is widely varying at each pixel, a potential load-balancing problem arises.

Agglomeration choices:

• Group pixels into continuous rows of the image; map continuous rows onto processors.
  • Advantage: this is the simplest choice both conceptually and implementation wise.

  • Disadvantage: processors which receive many "easy" rows will be left idle during much of the computation.

• Fully agglomerate pixels into 2-d blocks.
  • No difference from above; surface-to-volume benefits do not apply here.

• Semi-agglomerate into full rows. Dynamically map individual rows to processors at run time.
  • Advantage: processors with easy rows will receive more useful work once they finish their initial computation.

  • Advantage: task scheduling time is small; the master merely has to send a worker the row number.

  • Advantage: output rows are relatively easy to assemble into final image.

  • Disadvantage: requires significantly more programming work.

  • Disadvantage: will not scale well if the number of rows is not at least an order of magnitude greater than the number of processors.

• Do not agglomerate at all; the master sends pixel-sized tasks to processors at run time.
  • Advantage: processor load balance will almost surely be even.

  • Disadvantage: the overhead of creating a separate task for each pixel, and assembling the output image from byte-size messages will almost surely be greater then the pixel computation time itself.

Step 5: Mapping tasks to Processors

At this point we should have our program partitioned into reasonable sized tasks, but not yet mapped to processors.

In the case of static problems where we fully agglomerated T smaller tasks into P larger tasks, the mapping task is straightforward; one task per processor; SP2 tip.

In the tougher problems, the number of tasks is unknown; some form of dynamic scheduling and dynamic load balancing must be used.

Step 5: Mapping: Task-Scheduling Algorithms

Task-scheduling algorithms, where tasks are dynamically mapped to processors as computation proceeds, work best when each task has small communication requirements:

• Master/Worker: The simplest task scheduling algorithm.
  • A master process assigns tasks to workers; when a worker finishes one task, the master assigns it a new one.

  • Works best if there are an order of magnitude more tasks than processors.

  • Works best if the cost of assigning tasks to workers is small relative to their computation time.

  • Performance can be improved by prefetching a new task when the current one is about to complete.

• Hierarchical Master/Worker:
  • Similar to generic master/worker; a hierarchy of masters and workers are designed.

  • Can improve performance in problems where a single master would be overloaded with scheduling tasks to processors.

• Decentralized Task Pool:
  • Each processor maintains its own task pool.

  • If a processor exhausts its task pool, it requests tasks from neighbor processors.

  • Requires no central task manager.

Step 5: Mapping Checklist

1. If using a centralized load-balancing scheme, have you verified that the manager will not become a bottleneck?

2. If using probabilistic or cyclic methods, do you have a large enough number of tasks to ensure reasonable load balance?

3. If using a dynamic load-balancing scheme, have you evaluated the relative costs of different strategies (ie. probabilistic mappings)?

Step 5: Mapping Examples
Recursive Bisection Example

Step 5: Mapping Examples
Mandelbrot Image Calculation

Discussion

 

Case Study: 2-D FFT
Overview

The problem is described in a serial context, then the above steps are applied to yield a parallel program.

In this problem a series of one-dimensional FFT's are applied to a two-dimensional image.

The input and output of the algorithm is a 2-D image. A brief pseudo-code description of the algorithm is:

   do i = 1, IMAGE_HEIGHT
      row(tmp_image, i) = 1-d_fft(row(input_image, i))
   end do

   do i = 1, IMAGE_WIDTH
      column(output_image, i) = 1-d_fft(column(tmp_image, i))
   end do

2-D FFT Parallelization

Identify the Computation Hotspots: intuitive and visual inspection of the code reveals that the bulk of the run time is spent performing the N row FFT's and the N column FFT's. Running the program with prof confirms this.

Partitioning: for 2-D FFT, partitioning is the most important step and there are several choices:

• Partition all images by pixels: although this is the most aggressive partitioning possible, later steps would show it is unsatisfactory.
  • Visual examination of the code reveals that every output pixel in a 1-D FFT is a function of every input pixel to that 1-D FFT.

  • Communication and agglomeration steps would show that it's best to perform each 1-D FFT entirely in the same task.

• Partition all images by rows. This works perfectly for the row FFT's but is poor for the column FFT's.

• Partition all images by columns. See above.

• Duplicate the input image across all processors. No benefit arises because row-partitioning already provides communication-free row FFT's and the input image is not used for the column FFT's.

• Duplicate the tmp image on all processors; this would require either replicating the 1-D row FFT computations in all tasks (effectively eliminating parallelism for this phase) or gather/broadcasting the tmp image to all tasks (which could take longer than the row FFT computation itself): SP2 Communication Performance Summary

• Functional decomposition: define a task for each row FFT and each column FFT.
  • Use row partitioning for the row FFT's.

  • Use column partitioning for the column FFT's.

  • Repartition the tmp image in between 1-D FFT phases.

  • Advantage: both the row and column FFT tasks can proceed without communication.

  • Disadvantage: the repartitioning of the tmp image will introduce an AlltoAll communication overhead:SP2 Communication Performance Summary

Agglomeration/Mapping: This step, and the mapping step, are simple for this example.

• The computational run time is the same for each row or column FFT.

• No communication is required during the row and column FFT's.

• We fully agglomerate N/P row and column FFT's per task.

Source code for serial and parallel versions in C:

• mpi_2dfft.c
• mpi_2dfft.h
• ser_2dfft.c

Source code for serial and parallel versions in Fortran:

• mpi_2dfft.f
• ser_fft.f
• ser_2dfft.f
• timing_fgettod.c

Case Study: 2-D CFD
Overview

The problem is described in a serial context, then the above steps are applied to yield a parallel program.

This case study is a computational fluid dynamics example in two dimensions. In this example, we'll be simulating airflow over a wing. More specifically, this code calculates the solution of a compressible, inviscid, non-conducting flow over a body in two dimensions.

The code solves the 2D, compressible Euler equations on unstructured triangular grids:

• U_t + AU_x + BU_y = 0

The Euler equations, where U is the vector of conserved variables. A and B are the Jacobian matrices defined in the program.

The algorithm approximately solves the Euler equations over the domain by solving the discretized equations on the unstructured grid, or mesh. The values for density, x-momentum, y-momentum, and total energy are stored at the vertices of the mesh.

Basically, the field values at a given timestep are checked to see if they satisfy the discretized Euler equations. If not, then they are adjusted and checked again. This iterative process continues until the field values satisfy the discretized Euler equations to a specified tolerance.

2-D CFD: Computational Elements

The wing is divided into a triangular mesh. Each triangle is called a cell. Each cell has 3 vertices or nodes. Click here for an illustration of the 2-D mesh.

The "residual" is calculated at each cell using field values held at the vertices of that cell.

This fluctuation is then split-up and distributed to the vertices of the cell.

The sum of the fluctuation contributions from all the cells that surround a vertex make up the nodal fluctuation for that vertex.

Once the code has looped through all of the cells, the nodal residuals will have been calculated for each vertex, and the vertex values can be updated for the next timestep.

2-D CFD: Useful Information

• The grid is represented by a vector of N elements; each element contains the X and Y coordinates of the vertices, along with other data. The actual vertices appear in no particular order in this vector.

• CFD grids commonly have from thousands to millions of cells.

• Convergence is typically reached after 1,000 to 100,000 time iterations.

• Each cell needs approximately 400 floating-point operations per time iteration.

CFD Parallelization

Identify the Computation Hotspots: visual inspection of the code shows that the main loop in which the values (density, x-dir momentum, y-dir momentum and total energy) are calculated for each cell are the computationally expensive parts of the algorithm. File I/O, pre-, and post-processing steps will remain serial.

Partitioning: The cell is a natural element for data-partitioning:

• Each cell is essentially the same as every other cell.

• The computation performed at each cell is essentially the same as at every other cell.

• No obvious larger, or smaller, data partition elements appear to exist.

Communication: each cell only needs data from its vertices. Each vertex only needs data from vertices that it is directly connected to. Therefore, the 2-D grid itself already identifies communication links.

Agglomeration: at this point we have a separate task for each cell and far more cells than processors. We must agglomerate to larger tasks.

• We attempt to fully agglomerate cells/vertices to P separate tasks.
  • The total number of tasks is fixed.

  • The computation requirements of each task is the same.

  • The communication requirements of each task is the same.

• Since cells/vertices only need to communicate with cells/vertices they are directly adjacent to, we group physically close cells/vertices into the same agglomeration (simply putting the first N/P cells/vertices in task 0, etc. would not necessarily agglomerate physically close cells/vertices into the same task).

• We expect to reap surface-to-volume benefits due to the local nature of communication requirements.

• Dividing the physical 2-D space into equal sized blocks would generate an uneven load balance; the amount of cells/vertices per unit area of physical space is widely varying.

• We need to agglomerate cells/vertices such that the amount of cells per task is the same, and communication requirements between tasks is minimized.

• We choose a commercial package, METIS, which uses a recursive bisection technique to divide the grid into halves, halves into quarters, etc.

• The agglomerated tasks are illustrated here.

Source code (compressed tar format) for serial and parallel versions in C:

• Serial C files
• Parallel C files
• Grid input file
• Viewing utility

References and More Information

• Written September 1996 by George Gusciora

• This tutorial draws heavily from the book Designing and Building Parallel Programs by Ian Foster; the 5-step parallelization method, and many of the figures were borrowed from this book. Additional information can be obtained here.