You expect Speedup(n) to be less than n, reflecting the fact that not all parts of a program benefit from parallel execution. However, it is possible, in rare situations, for Speedup(n) to be larger than n. This is called a superlinear speedup—the program has been sped up by more than the increase of CPUs.

A superlinear speedup does not really result from parallel execution. It comes about because each CPU is now working on a smaller set of memory. The problem data handled by any one CPU fits better in cache, so each CPU executes faster than the single CPU could do. A superlinear speedup is welcome, but it indicates that the sequential program was being held back by cache effects.

If a properly parallelized program does not scale as well as expected, there are several potential causes. The first thing to check is whether some form of cache contention is slowing the program. You have presumably done this as part of single-CPU tuning (“Identifying Cache Problems with Perfex and SpeedShop” in Chapter 6), but new cache contention problems can appear when a program starts executing in parallel.

New problems can be an issue, however, only for data that are frequently updated or written. Data that are mostly read and rarely written do not cause cache coherency contention for parallel programs.

The mechanism used in SN0 to maintain cache coherence is described under “Understanding Directory-Based Coherency” in Chapter 1. When one CPU modifies a cache line, any other CPU that has a copy of the same cache line is notified, and discards its copy. If that CPU needs that cache line again, it fetches a new copy. This arrangement can cause performance problems in two cases:

In some cases, these problems arise because all parallel threads in the program are genuinely contending for use of the same key global variables. In that case, the only cure is an algorithmic change.

More often, cache contention arises because the CPUs are using and updating unrelated variables that happen to fall in the same cache line. This is false sharing. Fortunately, IRIX tools can help you identify these problems. Cache contention is revealed when perfex shows a high number of cache invalidation events: counter events 31, 12, 13, 28, and 29. A CPU that repeatedly updates the same cache line shows a high count of stores to shared cache lines: counter 31. (See “Cache Coherency Events” in Appendix B.)

False sharing can be demonstrated with code like that in Example 8-5.

      subroutine sum85 (a,s,m,n)
      integer m, n, i, j 
      real    a(m,n), s(m) 
!$omp parallel do private(i,j), shared(s,a)
      do i = 1, m 
         s(i) = 0.0 
         do j = 1, n 
            s(i) = s(i) + a(i,j) 
         enddo 
      enddo
      return
      end

This code calculates the sums of the rows of a matrix. For simplicity, assume m=4 and that the code is run on up to four CPUs. What you observe is that the time for the parallel runs is longer than when just one CPU is used. To understand what causes this, look at what happens when this loop is run in parallel.

The following is a time line of the operations that are carried out (more or less) simultaneously:

At each stage of the calculation, all four CPUs attempt to update one element of the sum array, s(i). For a CPU to update one element of s, it needs to gain exclusive access to the cache line holding that element, but the four words of s are probably contained in a single cache line, so only one CPU at a time can update an element of s. Instead of operating in parallel, the calculation is serialized.

Actually, it's a bit worse than merely serialized. For a CPU to gain exclusive access to a cache line, it first needs to invalidate cached copies that may reside in the other caches. Then it needs to read a fresh copy of the cache line from memory, because the invalidations will have caused data in some other CPU's cache to be written back to main memory. In a sequential version of the program, the element being updated can be kept in a register, but in the parallel version, false sharing forces the value to be continually reread from memory, in addition to serializing the updates.

This serialization is purely a result of the unfortunate accident that the different elements of s ended up in the same cache line. If each element were in a separate cache line, each CPU could keep a copy of the appropriate line in its cache, and the calculations could be done perfectly in parallel. A possible way to fix the problem is to spread the elements of s out so that each updated element resides in its own cache line, as shown in Example 8-6.

      subroutine sum86 (a,s,m,n)
      integer m, n, i, j 
      real    a(m,n), s(32,m)
!$omp parallel do private(i,j), shared(s,a)
      do i = 1, m 
         s(1,i) = 0.0 
         do j = 1, n 
            s(1,i) = s(1,i) + a(i,j) 
         enddo 
      enddo
      return
      end

The elements s(1,1), s(1,2), s(1,3) and s(1,4) are separated by at least 32 × 4 = 128 bytes, and so are guaranteed to fall in separate cache lines. Implemented this way, the code achieves perfect parallel speedup.

Another possible way to fix the problem in Example 8-5 is to replace all of the s(i)s inside the loop on i with a new variable, stemp, that is declared private. Then add a final statement inside the loop on i to assign stemp to s(i).

Note that the problem of false sharing is not specific to the SN0 architecture. It occurs in any cached-coherent shared memory system.

To see how this works for a real code, see an example from a paper presented at the Supercomputing '96 conference. One example in this paper is a weather modeling program that shows poor parallel scaling (the red curve in Figure 8-2).

Figure 8-2. Performance of Weather Model Before and After Tuning

Running the program under perfex revealed that the number of secondary data cache misses (event 26) increased as the number of CPUs increased, as did the number of stores exclusive to a shared block (event 31). The large event 31 counts, increasing with the secondary cache misses, indicated a likely problem with false sharing. (The large number of secondary cache misses were a problem as well.)

The source of the problem was found using ssrun. There are several events that can be profiled to determine where false sharing occurs. The natural one to use is event 31, “store/prefetch exclusive to shared block in scache.” There is no explicitly named experiment type for this event, so profiling it requires setting the following environment variables:

% setenv _SPEEDSHOP_HWC_COUNTER_NUMBER 31 
% setenv _SPEEDSHOP_HWC_COUNTER_OVERFLOW 99 
% ssrun -exp prof_hwc a.out 
% prof a.out.prof_hwc.m* a.out.prof_hwc.p*

(See “Sampling Through Other Hardware Counters” in Chapter 4.) You could also profile external interventions (event 12) or external invalidations (event 13), but using event 31 shows where the source of the problem is in the program (some thread asking for exclusive access to a cache line) rather than the place a thread happens to be when an invalidation or intervention occurs. Another event that can be profiled is secondary cache misses (event 26). For this event, use the -dsc_hwc experiment type and don't set any environment variables.

For this program, profiling secondary cache misses was sufficient to locate the source of the false sharing. The profile showed that the majority of secondary cache misses occurred in an accumulation step, similar to the row summation in Example 8-5. Padding was used to move each CPU's accumulation variable into its own cache line. After doing this, the performance improved dramatically (the blue curve in Figure 8-2).  

You can often deal with cache contention by changing the layout of data in the source program, but sometimes you may have to make algorithmic changes as well. If profiling indicates that there is cache contention, examine the parallel regions identified; any assignment to memory in the parallel regions is a possible source of the contention. You need to determine if the assigned variable, or the data adjacent to it, is used by other CPUs at the same time. If so, the assignment forces the other CPUs to read a fresh copy of the cache line, and this is a source of contention.

To deal with cache contention, you have the following general strategies:

An update of one word (that is, a 4-byte quantity) invalidates all the 31 other words in the same L2 cache line. When those other words are not related to the new data, false sharing results. Use strategy 3 to eliminate the false sharing.

Be careful when your program defines a group of global status variables that is visible to all parallel threads. In the normal course of running the program, every CPU will cache a copy of most or all of this common area. Shared, read-only access does no harm. But if items in the block are volatile (frequently updated), those cache lines are invalidated often. For example, a global status area might contain the anchor for a LIFO queue. Each time a thread puts or takes an item from the queue, it updates the queue head, invalidating that cache line.

It is inevitable that a queue anchor field will be frequently invalidated. The time cost, however, can be isolated to the code that accesses the queue by applying strategy 2. Allocate the queue anchor in separate memory from the global status area. Put only a pointer to the queue anchor (a non-volatile item) in the global status block. Now the cost of fetching the queue anchor is born only by CPUs that access the queue. If there are other items that are updated along with the queue anchor—such as a lock that controls exclusive access to the queue—place them adjacent to the queue, aligned so that all are in the same cache line (strategy 4). However, if there are two queues that are updated at unrelated times, place the anchor of each in its own cache line (strategy 3).

Synchronization objects such as locks, semaphores, and message queues are global variables that must be updated by each CPU that uses them. You may as well assume that synchronization objects are always accessed at memory speeds, not cache speeds. You can do two things to reduce contention:

When you make a loop run in parallel, try to ensure that each CPU operates on its own distinct sections of the input and output arrays. Sometimes this falls out naturally, but there are also compiler directives for just this purpose. (These are described in “Using Data Distribution Directives”.)

Carefully review the design of any data collections that are used by parallel code. For example, the root and the first few branches of a binary tree are likely to be visited by every CPU that searches that tree, and they will be cached by every CPU. However, elements at higher levels in the tree may be visited by only a few CPUs. One option is to pre-build the top levels of the tree so that these levels never have to be updated once the program starts. Also, before you implement a balanced-tree algorithm, consider that tree-balancing can propagate modifications all the way to the root. It might be better to cut off balancing at a certain level and never disturb the lower levels of the tree. (Similar arguments apply to B-trees and other branching structures: the “thick” parts of the tree are widely cached and should be updated least often, while the twigs are less frequently used.)

Other classic data structures can cause memory contention, and algorithmic changes are needed to cure it:

Data placement policies used by the operating system are introduced under “Data Placement Policies” in Chapter 2. The default policy is called first-touch. Under this policy, the process that first touches (that is, writes to, or reads from) a page of memory causes that page to be allocated in the node on which the process is running. This policy works well for sequential programs and for many parallel programs as well. For example, this is just what you want for message-passing programs that run on SN0. In such programs each process has its own separate data space. Except for messages sent between the processes, all processes use memory that should be local. Each process initializes its own data, so memory is allocated from the node the process is running in, thus making the accesses local.

But for some parallel programs, the first-touch policy can have unintended side effects. As an example, consider a program parallelized using the MP library. In parallelizing such a program, the programmer worries only about parallelizing those parts of the program that take the most amount of time. Often, data initialization code takes little time, so they are not parallelized, and so are executed by the main thread of the program. Under the first-touch policy, all the program's memory ends up allocated in the node running the main thread. Having all the data concentrated in one node or within a small radius of it creates a bottleneck: all data accesses are satisfied by one hub, and this limits the memory bandwidth. If the program is run on only a few CPUs, you may not notice a bottleneck. But as you add more CPUs, the one hub that serves the memory becomes saturated, and the speed does not scale as predicted by Amdahl's law.

One easy way to test for this problem is to try other memory allocation policies. The first one you should try is round-robin allocation. Under this policy, data are allocated in a round-robin fashion from all the nodes the program runs on. Thus, even if the data are initialized sequentially, the memory holding them will not be allocated from a single node; it will be evenly spread out among all the nodes running the program. This may not place the data in their optimal locations—that is, the access times are not likely to be minimized—but there will be no bottlenecks, so scalability will not be limited.

You can enable round-robin data placement for an MP library program with the following environment variable setting (see pe_environ(5) man page) :

% setenv _DSM_ROUND_ROBIN 

The performance improvement this can make is demonstrated with an example. The plot in Figure 8-3 shows the performance achieved for three parallelized runs of the double-precision vector operation a(i) = b(i) + q*c(i) .

Figure 8-3. Calculated Bandwidth for Different Placement Policies

The bandwidth calculated assumes that 24 bytes are moved for each vector element. The red curve (“placement in one node”) shows the performance you attain using a first-touch placement policy and a sequential initialization so that all the data end up in the memory of one node. The memory bottleneck this creates severely limits the scaling. The green curve (“round-robin placement”) shows what happens when the data placement policy is changed to round-robin. The bottleneck disappears and the performance now scales with the number of CPUs. The performance isn't ideal—the blue line (“optimal placement”) shows what you measure if the data are placed optimally—but by spreading the data around, the round-robin policy makes the performance respectable. (Ideal data placement cannot be accomplished simply by initializing the data in parallel.)

If a round-robin data placement solves the scaling problem that led you to tune the data placement, then you're done! But if the performance is still not up to your expectations, the next experiment to perform is enabling dynamic page migration.

The IRIX automatic page migration facility is disabled by default because page migration is an expensive operation that impacts all CPUs, not just the ones used by the program whose data are being moved. You can enable dynamic page migration for a specific MP library program by setting the environment variable _DSM_MIGRATION. You can set it to either ON or ALL_ON:

% setenv _DSM_MIGRATION ON 
% setenv _DSM_MIGRATION ALL_ON 

When set to ALL_ON, all program data is subject to migration. When set only to ON, only those data structures that have not been explicitly placed via the compiler's data distribution directives (“Using Data Distribution Directives”) will be migrated.

Enabling migration is beneficial when a poor initial placement of the data is responsible for limited performance. However, it cannot be effective unless the program runs for at least half a minute, first because the operating system needs some time to move the data to the best layout, and second because the program needs to execute for some time with the pages in the best locations in order to recover the time cost of migrating.

In addition to the MP library environment variables, migration can also be enabled as follows:

You can try to reduce the cost of migrating from a poor initial location to the optimal one by combining a round-robin initial placement with migration. You won't know whether round-robin, migration, or both together will produce the best results unless you try the different combinations, so tuning the data placement requires experimentation. Fortunately, the environment variables make these experiments easy to perform.

Figure 8-4 shows the results of several experiments on the vector operation

a(i) = b(i) + q*c(i) 

In addition to combinations of the round-robin policy and migration, the effect of serial and parallel initializations are shown. For the results presented in the diagram, the vector operation was iterated 100 times and the time per iteration was plotted to see the performance improvement as the data moved toward optimal layout.

Figure 8-4. Calculated Iteration Times for Different Placement Policies

The curve at the top (“fixed placement in one node”) shows the performance for the poor initial placement caused by a sequential data initialization with a first-touch policy. No migration is used, so this performance is constant over all iterations. The line near the bottom (“fixed round-robin placement”) shows that a round-robin initial placement fixes most of the problems with the serial initialization: the time per iteration is five times faster than with the first-touch policy. Ultimately, though, the performance is still a factor of two slower than if an optimal data placement had been used. Once again, migration is not enabled, so the performance is constant over all iterations.

The flat curve just below the curve shows the performance achieved when a parallel initialization is used in conjunction with a round-robin policy. Its per-iteration time is constant and nearly the same as when round-robin is used with a sequential initialization. This is to be expected: the round-robin policy spreads the data evenly among the CPUs, so it matters little how the data are initialized.

The remaining five curves all asymptote to the optimal time per iteration. Four of these eventually achieve the optimal time due to the use of migration. The teal curve (at the bottom of the figure) uses a parallel initialization with a first-touch policy to achieve the optimal data placement from the outset. The orange curve just above it starts out the same, but for it, migration is enabled. The result of using migration on perfectly placed data is only to scramble the pages around before finally settling down on the ideal time. Above this curve are a magenta and a black curve (round-robin initial placement). They show the effect of combining migration and a round-robin policy. The only difference is that the magenta curve used a serial initialization while the black curve used a parallel initialization. For these runs, the serial initialization took longer to reach the steady-state performance (the black curve doesn't show up well in the figure), but you should not conclude that this is a general property; when a round-robin policy is used, it doesn't matter how the data are initialized.

Finally, the blue curve (“initial placement in one node”) shows how migration improves a poor initial placement in which all the data began in one node. It takes longer to migrate to the optimal placement than the other cases, indicating that by combining a round-robin policy with migration, one can do better than by using just one of the remedies alone. Dramatic evidence of this is seen when the cumulative program run time is plotted, as shown in Figure 8-5.   

Figure 8-5. Cumulative Run Time for Different Placement Policies

These results show clearly that, for this program, migration eventually gets the data to their optimal location and that migrating from an initial round-robin placement is faster than migrating from an initial placement in which all the data are on one node. Note, though, that migration does take time: if only a small number of iterations are going to be performed, it is better to use a round-robin policy without migration.

The results for data migration were generated using this environment variable setting:

% setenv _DSM_MIGRATION ON 

(This program has no data distribution directives in it, so in this case, there is no difference between a setting of ON and one of ALL_ON.) But you can also adjust the migration level, that is, control how aggressively the operating system migrates data from one memory to another. The migration level is controlled with the following environment variable setting (see pe_environ(5) man page):

% setenv _DSM_MIGRATION_LEVEL level 

A high level (100 maximum) means aggressive migration, and a low level means nonaggressive migration. A level of 0 means no migration. The migration level can be used to experiment with how aggressively migration is performed. The plot in Figure 8-6 shows the improvement in time per iteration for various migration levels as memory is migrated from an initial placement in one node. As you would expect, with migration enabled, the data will eventually end up in an optimal location, while the more aggressive migration is, the faster the pages get to where they should be.

Figure 8-6. Effect of Migration Level on Iteration Time

The best level for a particular program will depend on how much time is available to move the data to their optimal location. A low migration level is fine for a long-running program with a poor initial data placement. But if a program can only afford a small amount of time for the data to redistribute, a more aggressive migration level is needed. Keep in mind that migration has an impact on system-wide performance, and the more aggressive the migration, the greater the impact. For most programs the default migration level setting is 100 (agressive).

In general, you should try the following experiments when you need to tune data placement:

Note that these data placement tuning experiments apply only to MP library programs:

If there is a single placement of data that is optimal for your program, you can use the default first-touch policy to cause each CPU's share of the data to be allocated from memory local to its node. As a simple example, consider parallelizing the vector operation in Example 8-7.

integer i, j, n, niters 
      parameter (n = 8*1024*1024, ndim = n+35, niters = 100) 
      real a(ndim), b(ndim), q 
c initialization 
      do i = 1, n 
         a(i) = 1.0 - 0.5*i 
         b(i) = -10.0 + 0.01*(i*i) 
      enddo 
c real work 
      do it = 1, niters 
         q = 0.01*it 
         do i = 1, n 
            a(i) = a(i) + q*b(i) 
         enddo 
      enddo
      print *, a(1), a(n), q
      end

This vector operation is easy to parallelize; the work can be divided among the CPUs of a shared memory computer any way you would like. For example, if p is the number of CPUs, the first CPU can carry out the first n/p iterations, the second CPU the next n/p iterations, and so on. (This is called SCHEDULE (STATIC).) Alternatively, each thread can perform one of the first p iterations, then one of the next p iterations and so on. (This is called SCHEDULE (STATIC,1).) If Example 8-7 is compiled with option -O2 and Example 8-8 is compiled with options -mp -O2, the speedup is approximately 3.2 when four processors are specified for the code using the OpenMP directives.

In a cache-based machine, not all divisions of work produce the same performance. The reason is that if a CPU accesses the element a(i), the entire cache line containing a(i) is moved into its cache. If the same CPU works on the following elements, they will be in cache. But if different CPUs work on the following elements, the cache line will have to be loaded into each one's cache. Even worse, false sharing is likely to occur. Thus performance is best for work allocations in which each CPU is responsible for blocks of consecutive elements.

Example 8-8 shows the above vector operation parallelized for SN0.

      integer i, j, n, niters 
      parameter (n = 8*1024*1024, niters = 1000) 
      real a(n), b(n), q 
c initialization 
c$doacross local(i), shared(a,b) 
      do i = 1, n 
         a(i) = 1.0 - 0.5*i 
         b(i) = -10.0 + 0.01*(i*i) 
      enddo 
c real work 
      do it = 1, niters 
         q = 0.01*it 
c$doacross local(i), shared(a,b,q) 
         do i = 1, n 
            a(i) = a(i) + q*b(i) 
         enddo 
      enddo 

Because the schedule type is not specified, it defaults to simple: that is., process 0 performs iterations 1 to n/p, process 1 performs iterations 1 + (n/p) to 2 × n/p, and so on. Each process accesses blocks of memory with stride 1.

Because the initialization takes a small amount of time compared with the “real work,” parallelizing it doesn't reduce the sequential time of this code by much. Some programmers wouldn't bother to parallelize the first loop for a traditional shared memory computer. However, if you are relying on the first-touch policy to ensure a good data placement, it is critical to parallelize the initialization code in exactly the same way as the processing loop.

Due to the correspondence of iteration number with data element, the parallelization of the “real work” loop means that elements 1 to n/p of the vectors a and b are accessed by process 0. To minimize memory access times, you would like these data elements to be allocated from the memory of the node running process 0. Similarly, elements 1 + (n/p) to 2 × n/p are accessed by process 1, and you would like them allocated from the memory of the node running it, and so on. This is accomplished automatically by the first-touch policy. The CPU that first touches a data element causes the page holding that data element to be allocated from its local memory. Thus, if the data are to be allocated so that each CPU can make local accesses during the “real work” section of the code, each CPU must be the one to initialize its share of the data. This means that the initialization loop must be parallelized the same way as the “real work” loop.

Now consider why the simple schedule type is important to Example 8-8. Data are placed in units of one page, so they will end up in their optimal location only if the same CPU processes all the data elements in a page. The default page size is 16 KB, or 4096 REAL*4 data elements; this is a fairly large number. Since the static schedule blocks together as many elements as possible for a single CPU to work on (n/p), it will create more optimally-allocated pages than any other work distribution.

Figure 8-7. Effect of Page Granularity in First-Touch Allocation

For Example 8-8, n = 8 × 1024 × 1024 = 8388608. If the program is run on 128 CPUs, n/p = 65536, which means that each CPU's share of each array fills sixteen 16 KB pages ((65536 elements × 4) ÷ 16 KB). However, it is unlikely that any array begins exactly on a page boundary, so you would expect 15 of a CPU's 16 pages to contain only elements it processes, and two additional pages, each of which will contain some elements of the appropriate CPU along with some elements belonging to a neighboring CPU. This effect is diagrammed in Figure 8-7.

In a perfect data layout, no pages would share elements from multiple CPUs, but this small imperfection has a negligible effect on performance.

On the other hand, if you used a SCHEDULE(STATIC,1) clause, all 128 CPUs repeatedly attempt to concurrently touch 128 adjacent data elements. Because 128 adjacent data elements are almost always on the same page, the resulting data placement could be anything from a random distribution of the pages, to one in which all pages end up in the memory of a single CPU. This initial data placement will certainly affect the performance.

In general, you should try to arrange it so that each CPU's share of a data structure exceeds the size of a page. For finer granularity, you need the directives discussed under “Using Reshaped Distribution Directives”.

Programming for the default first-touch policy is effective if the application requires only one distribution of data. When different data placements are required at different points in the program, more sophisticated techniques are required. The data directives (see “Using Data Distribution Directives”) allow data structures to be redistributed at execution time, so they provide the most general approach to handling multiple data placements.

In many cases, you don't need to go to this extra effort. For example, the initial data placement can be optimized with first-touch, and migration can redistribute data during the run of the program. In other cases, copying can be used to handle multiple data placements. As an example, consider the FT kernel of the NAS FT benchmark (see the link under “Third-Party Resources”). This is a three-dimensional Fast Fourier Transform (FFT) carried out multiple times. It is implemented by performing one-dimensional FFTs in each of the three dimensions. As long as you don't need more than 128-way parallelism, the x- and y-calculations can be parallelized by partitioning data along the z-axis, as sketched in Figure 8-8.

Figure 8-8. Data Partition for NAS FT Kernel

For SN0, first-touch is used to place nz/p, xy-planes in the local memory of each CPU. The CPUs are then responsible for performing the FFTs in their shares of planes.

Once the FFTs are complete, z-transforms need to be done. On a conventional distributed memory computer, this presents a problem. The xy-planes have been distributed, no CPU has access to all the data along the lines in the z-dimension, so the z-FFTs cannot be performed. The conventional solution is to redistribute the data by transposing the array, so that each CPU holds a set of zy-planes. One-dimensional FFTs can then be carried out along the first dimension in each plane. A transpose is convenient since it is a well-defined operation that can be optimized and packaged into a library routine. In addition, it moves a lot of data together, so it is more efficient than moving individual elements to perform z-FFTs one at a time. (This technique was discussed earlier at “Understanding Transpositions” in Chapter 6.)

SN0, however, is a shared memory computer, so explicit data redistribution is not needed. Instead, split up the z-FFTs among the CPUs by partitioning work along the y-axis, as sketched in Figure 8-9.

Figure 8-9. NAS FT Kernel Data Redistributed

Each CPU then copies several z-lines at a time into a scratch array, performs the z-FFTs on the copied data, and completes the operation by copying back the results.

Copying has several advantages over explicitly redistributing the data:

Combining first-touch placement with copying solves a problem that might otherwise require two data distributions.    

In Example 8-9, Distribute specifies a BLOCK mapping for both arrays. BLOCK means that, when running with p CPUs, each array is to be divided into p blocks of size ceiling(n/p), with the first block assigned to the first CPU, the second block assigned to the second CPU, and so on. The intent of assigning blocks to CPUs is to allow each block to be stored in one CPUs' local memory.

Only whole pages are distributed, so when a page straddles blocks belonging to different CPUs, it is stored entirely in the memory of a single node. As a result, some CPUs will use nonlocal accesses to some of the data on one or two pages per array. (The situation is diagrammed in Figure 8-7.) An imperfect data distribution has a negligible effect on performance because, as long as a “block” comprises at least a few pages, the ratio of nonlocal to local accesses is small. When the block size is less than a page, you must live with a larger fraction of nonlocal accesses, or use the reshaped directives. (Data placement is rarely important for arrays smaller than a page, because if they are used heavily, they fit entirely in cache.)

Now look at the Parallel Do directive in Example 8-9. It instructs the compiler to run this loop in parallel. But instead of distributing iterations to CPUs by specifying a schedule type such as SCHEDULE(STATIC) or SCHEDULE(STATIC,1), the AFFINITY clause is used. This new clause tells the compiler to execute iteration i on the CPU that is responsible for data element a(i) under data distribution. Thus work is divided among the CPUs according to their access to local data. This is the situation achieved using first-touch allocation and default scheduling, but now it is specified explicitly. The AFFINITY clause is placed on a separate line beginning with !$SGI+ so that the clause winon-SGI ll be ignored if the program is run on a platform other than SGI.

Page granularity is not considered in assigning work; the first CPU carries out iterations 1 to ceiling(n/p), the second CPU performs iterations ceiling(n/p) + 1 to 2 × ceiling(n/p), and so on. This ensures a proper balance of work among CPUs, at the possible expense of a few nonlocal memory accesses.

For the default BLOCK distribution used in this example, an  AFFINITY clause assigns work to CPUs identically to the default SCHEDULE(STATIC). The  AFFINITY clause, however, has advantages. First, for some complicated data distributions, there are no defined schedule types that can achieve optimal work distribution. Second, the  AFFINITY clause makes the code easier to maintain, because you can change the data distribution without having to modify the Distribute directive to realign the work with the data.

The Distribute directive allows you specify a different mapping for each dimension of each distributed array. The possible mappings are these:

For example, suppose a Fortran program has a two-dimensional array A and uses p CPUs. Then !$SGI DISTRIBUTE A(*,CYCLIC) assigns columns 1, p+1, 2p+1, ... of A to the first CPU; columns 2, p+2, 2p+2, ... to the second CPU; and so on. CYCLIC(x) specifies a block-cyclic mapping, in which blocks, rather than individual elements, are cyclically assigned to CPUs, with the block size given by the compile-time expression x. For example, !$SGI DISTRIBUTE A(*,CYCLIC(2)) assigns columns 1, 2, 2p+1, 2p+2, ... to the first CPU; 3, 4, 2p+3, 2p+4, ... to the second CPU; and so on.

Combinations of the mappings on different dimensions can produce a variety of distributions, as shown in Figure 8-10, in which each color represents the intended assignment of data elements to one CPU out of a total of four CPUs.

Figure 8-10. Some Possible Regular Distributions for Four Processors

The distributions illustrated in Figure 8-10 are ideals, and do not take into consideration the restriction to whole pages. The intended distribution of data is achieved only when at least a page of data is assigned to a CPU's local memory. In the figure, mappings (a) and (d) produce the desired data distributions for arrays of moderate size, while mappings (e) and (f) require large arrays for the data placements to be near-optimal. For the cyclic mappings (b, c, and g), reshaped data distribution directives should be used to achieve the intended results (“Using Reshaped Distribution Directives”).

When an array is distributed in multiple dimensions, data is apportioned to CPUs as equally as possible across each dimension. For example, if an array has two distributed dimensions, execution on six CPUs assigns the first dimension to two CPUs and the second to three (2x 3= 6). The optional ONTO clause allows you to override this default and explicitly control the number of CPUs in each dimension. The clause ONTO(3,2) assigns three CPUs to the first dimension and two to the second. Some possible arrangements are sketched in Figure 8-11.

Figure 8-11. Possible Outcomes of Distribute ONTO Clause

The arguments of ONTO specify the aspect ratio desired. If the available CPUs cannot exactly match the specified aspect ratio, the best approximation is used. In a six-CPU execution, ONTO(1,2) assigns two CPUs to the first dimension and three to the second. An argument of asterisk means that all remaining CPUs are to fill in that dimension: ONTO(2,*) assigns two CPUs to the first dimension and p/2 to the second.

The  AFFINITY clause extends the Parallel Do directive to data distribution. Affinity has two forms. The data form, in which iterations are assigned to CPUs to match a data distribution, is used in Example 8-9. In that example, the correspondence between iterations and data distribution is quite simple. The clause, however, allows more complicated associations. One is shown in Example 8-10.

parameter (n=800)
      real a(2*n+3,0:n)
!$sgi distribute a(block,cyclic(1))
      do j=0,n
        do i=1,2*n+3
          a(i,j)=1000*i+j
        enddo
      enddo
      print *,a(5,0),a(5,1)
!$omp parallel do private(i,j), shared(a)
!$sgi+     affinity(i) = data(a(2*i+3,j))
      do i=1,n
        do j=1,n
          a(2*i+3,j)=a(2*i+3,j-1)
        enddo
      enddo
      print *,a(5,0),a(5,1)
      end

Here the compiler and runtime try to execute loop iterations on i in those CPUs for which element A(2*i+3,j) is in local memory. The loop-index variable (i in the example) cannot appear in more than one dimension of the clause, and the expressions involving it are limited to the form m x i + k, where m and k are integer constants with m greater than zero.

A different kind of affinity relates iterations to threads of execution, without regard to data location. This form of the clause is shown in Example 8-11.

integer n, p, i 
      parameter (n = 8*1024*1024) 
      real a(n) 
      p = 1 
!$    p = mp_numthreads() 
!$omp parallel do private(i), shared (a,p) 
!$sgi+     affinity(i) = thread(i/((n+p-1)/p)) 
      do i = 1, n 
         a(i) = 0.0 
      enddo
      print *, n, a(n)
      end

This form of affinity has no direct relationship to data distribution. Rather, it replaces the schedule(type) clause (that is, #pragma omp for schedule(...), or !$OMP DO SCHEDULE) with a schedule based on an expression in an index variable. The code in Example 8-11 executes iteration i on the thread given by an expression, modulo the number of threads that exist. The expression may need to be evaluated in each iteration of the loop, so variables (other than the loop index) must be declared shared and not changed during the execution of the loop.

Multi-dimensional mappings (for example, mappings f and g in Figure 8-10) often benefit from parallelization over more than just the one loop to which a standard directive applies. The NEST clause permits such parallelizations. It specifies that the full set of iterations in the loop nest may be executed concurrently. Typical syntax is shown in Example 8-12.

parameter (n=1000)
      real a(n,n) 
!$omp parallel do private(i,j), shared(a) 
!$sgi+ nest(i,j)
      do j = 1, n 
         do i = 1, n 
            a(i,j) = 0.0 
         enddo 
      enddo 
      print *, n, a(n,n)
      end

The normal and default parallelization of the loop in Example 8-12 would schedule the iterations of each column of the array—iterations in i for each value of j—on parallel CPUs. The work of the outer loop, which is only index incrementing, is done serially. In the example, the NEST clause specifies that all n x n iterations can be executed in parallel. To use this directive, the loops to be parallelized must be perfectly nested, that is, no code is allowed except within the innermost loop. Using NEST, when there are enough CPUs, every combination of i and j executes in parallel. In any event, when there are more than n CPUs, all will have work to do.

parameter (n=1000)
      real a(n,n) 
!$sgi distribute a(block,block)
!$omp parallel do private(i,j), shared(a) 
!$sgi+ nest(i,j), affinity(i,j) = data(a(i,j))
      do j = 1, n 
         do i = 1, n 
            a(i,j) = 0.0 
         enddo 
      enddo 
      print *, n, a(n,n)
      end

Example 8-13 shows the combination of the NEST and AFFINITY (to data) clauses. All iterations of the loop are treated independently, but each CPU operates on its share of distributed data.

parameter (m=2000,n=1000)
      real a(m,n) 
!$sgi distribute a(block,block) onto (2,1)
!$omp parallel do private(i,j), shared(a) 
!$sgi+ nest(i,j), affinity(i,j) = data(a(i,j))
      do j = 1, n 
         do i = 1, m 
            a(i,j) = 0.0 
         enddo 
      enddo 
      print *, m, n, a(m,n)
      end

Example 8-14 shows the combination of all three clauses: NEST to parallelize all iterations; AFFINITY to causes CPUs to work on local data; and an ONTO clause to specify the aspect ratio of the parallelization. In this example, twice as many CPUs are used to parallelize the i-dimension as the j-dimension.

The Distribute directive is commonly used to specify a single arrangement of data that is constant throughout the program. In some programs you need to distribute an array differently during different phases of the program. Two directives allow you to accomplish this easily:

Redistribute allows you to change the distribution defined by a Distribute directive to another mapping. However, redistribution should not be seen as a general performance tool. For example, the FFT algorithm is better with data copying, as discussed under “First-Touch Placement with Multiple Data Distributions”.

Despite its name, the Distribute_Reshape directive does not perform redistribution. It tells the compiler that it can reshape the distributed data into nonstandard layouts.