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.

The purpose of the Distribute_Reshape directive is to achieve the desired data distribution without page-granularity limitations. You can accomplish the same thing through the use of dynamic memory allocation and pointers. The method follows this outline:

However, this approach can easily add great complexity to what ought to be a readable algorithm. As with most optimizations, it is usually better to leave the source code simple and let the compiler do the hard work behind the scenes.

When you use the Distribute_Reshape directive, there are some restrictions you should remember:

In addition, some care must be used in passing reshaped arrays to subroutines. As long as you pass the entire array, and the declared size of the array in the subroutine matches the size of array that was passed to it, problems are unlikely. Example 8-21 illustrates this.

      real a(n1,n2) 
!$sgi distribute_reshape a(block,block) 
      call sub(a,n1,n2) 
               . . .
      subroutine sub(a,n1,n2) 
      real a(n1,n2) 
                . . .

In Example 8-22, the formal parameter is declared differently from the actual reshaped array. Treating the two-dimensional array as a vector amounts to an assumption that consecutive elements are adjacent in memory, a false assumption. The compiler will catch this change of dimensionality and issue a warning message if you compile with the -mp option.

      real a(n1,n2) 
!$sgi distribute_reshape a(block,block) 
      call sub(a,n1*n2) 
                . . .
      subroutine sub(a,n) 
      real a(n) 
                . . .

Inside a subroutine, you don't need to declare how an array has been distributed. In fact, the subroutine is more general if you do not declare the distribution. The compiler will generate versions of the subroutine necessary to handle all distributions that are passed to it. Example 8-23 shows a subroutine receiving two arrays that are reshaped differently.

      real a(n1,n2) 
!$sgi distribute_reshape a(block,block) 
      real b(n1,n2) 
!$sgi distribute_reshape b(cyclic,cyclic) 
      call sub(a,n1,n2) 
      call sub(b,n1,n2) 
                . . .
      subroutine sub(a,n1,n2) 
      real a(n1,n2) 
                . . .

The compiler actually generates code to handle each distribution passed to the subroutine. As long as the data calculations are efficient for both distributions, this will achieve good performance. If a particular algorithm works only for a specific data distribution, you can declare the required distribution inside the subroutine by using a Distribute_Reshape directive there. Then, calls to the subroutine that pass a mismatched distribution cause compile- or link-time errors.

Most errors in accessing reshaped arrays are caught at compile time or link time. However, some errors, such as passing a portion of a reshaped array that spans a distributed dimension, can be caught only at run time. You can instruct the compiler to do runtime checks for these with the -MP:check_reshape=on flag. You should use this flag during the development and debugging of programs which use reshaped arrays. In addition, the Fortran manuals listed in “Compiler Manuals” have more information about reshaped arrays.

If you are using the data distribution directives, you can set the environment variable _DSM_VERBOSE. When it is set, the program will print out messages at run time that tell you whether the distribution directives are being used, which physical CPUs the threads are running on, and what the page sizes are. A report might look like Example 8-24.

% setenv _DSM_VERBOSE 
% a.out 
[ 0]: 16 CPUs, using 4 threads. 
        1 CPUs per memory 
        Migration: OFF 
        MACHINE: IP27 --- is NUMA --- 
        MACHINE: IP27 --- is NUMA --- 
Created 4 MLDs 
Created and placed MLD-set. Topology-free, Advisory. 
MLD attachments are: 
        MLD  0. Node /hw/module/3/slot/n4/node 
        MLD  1. Node /hw/module/3/slot/n3/node 
        MLD  2. Node /hw/module/3/slot/n1/node 
        MLD  3. Node /hw/module/3/slot/n2/node 
[ 0]: process_mldlink: thread  0 (pid  3832) to memory  0. -advisory-. 
Pagesize: stack 16384 data 16384 text 16384 (bytes) 
--- finished MLD initialization --- 
                                . . . 
[ 0]: process_mldlink: thread  1 (pid  3797) to memory  1. -advisory-. 
[ 0]: process_mldlink: thread  2 (pid  3828) to memory  2. -advisory-. 
[ 0]: process_mldlink: thread  3 (pid  3843) to memory  3. -advisory-. 
[ 0]: process_mldlink: thread  0 (pid  3832) to memory  0. -advisory-. 

These messages are useful in verifying that the requested data placements are actually being performed.

The utility dplace (see “Non-MP Library Programs and Dplace”) disables the data distribution directives. It should not normally be used to run MP library programs. If you do use it, however, it displays messages reminding you that “DSM” is disabled.

You can obtain information about data placement within the program with the dsm_home_threadnum() intrinsic. It takes an address as an argument, and returns the number of the CPU in whose local memory the page containing that address is stored. It is used in Fortran as follows:

integer dsm_home_threadnum 
numthread = dsm_home_threadnum(array(i)) 

Two CPUs are connected to each node and they share the same memory, so the function returns the lowest CPU number of the two running on the node with the data.

The numbering of CPUs is relative to the program. You can also determine which absolute physical node a page of memory is stored in. (Note, though, that this number will typically change from run to run.) The information is gotten via the syssgi() system call using the SGI_PHYSP command (see the syssgi(2) man page). The routine va2pa() translates a virtual address to a physical address using this system call (see Example C-9). You can translate a physical page address to a node number with the following macro:

#define ADDR2NODE(A) ((int) (va2pa(A) >> 32)) 

Note that the node number is generally half the lowest numbered CPU on the node, so if you would prefer to use CPU numbers, multiply the return value by two.

To explore this information, a test program was written that allocated two vectors, a and b, of size 256 KB, initialized them sequentially, and then printed out the locations of their pages. The default first-touch policy for a run using five CPUs resulted in the page placements shown in Example 8-25.

Distribution for array "a" 
address       byte index  thread  proc 
-------------------------------------- 
0xfffffb4760  0x       0       0    12 
0xfffffb8000  0x    38a0       0    12 
0xfffffbc000  0x    78a0       0    12 
0xfffffc0000  0x    b8a0       0    12 
0xfffffc4000  0x    f8a0       0    12 
0xfffffc8000  0x   138a0       0    12 
0xfffffcc000  0x   178a0       0    12 
0xfffffd0000  0x   1b8a0       0    12 
0xfffffd4000  0x   1f8a0       0    12 
0xfffffd8000  0x   238a0       0    12 
0xfffffdc000  0x   278a0       0    12 
0xfffffe0000  0x   2b8a0       0    12 
0xfffffe4000  0x   2f8a0       0    12 
0xfffffe8000  0x   338a0       0    12 
0xfffffec000  0x   378a0       0    12 
0xffffff0000  0x   3b8a0       0    12 
0xffffff4000  0x   3f8a0       0    12 
Distribution for array "b" 
address       byte index  thread  proc 
-------------------------------------- 
0xfffff6e0f8  0x       0       0    12 
0xfffff70000  0x    1f08       0    12 
0xfffff74000  0x    5f08       0    12 
0xfffff78000  0x    9f08       0    12 
0xfffff7c000  0x    df08       0    12 
0xfffff80000  0x   11f08       0    12 
0xfffff84000  0x   15f08       0    12 
0xfffff88000  0x   19f08       0    12 
0xfffff8c000  0x   1df08       0    12 
0xfffff90000  0x   21f08       0    12 
0xfffff94000  0x   25f08       0    12 
0xfffff98000  0x   29f08       0    12 
0xfffff9c000  0x   2df08       0    12 
0xfffffa0000  0x   31f08       0    12 
0xfffffa4000  0x   35f08       0    12 
0xfffffa8000  0x   39f08       0    12 
0xfffffac000  0x   3df08       0    12 

The sequential initialization combined with the first-touch policy caused all the pages to be placed in the memory of the first thread which, in this run, was CPU 12, or node 6.

Example 8-26 shows the results from the same code after setting the environment variable _DSM_ROUND_ROBIN to use a round-robin page placement.

Distribution for array "a" 
address       byte index  thread  proc 
-------------------------------------- 
0xfffffb4750  0x       0       2     2 
0xfffffb8000  0x    38b0       0     0 
0xfffffbc000  0x    78b0       4     4 
0xfffffc0000  0x    b8b0       2     2 
0xfffffc4000  0x    f8b0       0     0 
0xfffffc8000  0x   138b0       4     4 
0xfffffcc000  0x   178b0       2     2 
0xfffffd0000  0x   1b8b0       0     0 
0xfffffd4000  0x   1f8b0       4     4 
0xfffffd8000  0x   238b0       2     2 
0xfffffdc000  0x   278b0       0     0 
0xfffffe0000  0x   2b8b0       4     4 
0xfffffe4000  0x   2f8b0       2     2 
0xfffffe8000  0x   338b0       0     0 
0xfffffec000  0x   378b0       4     4 
0xffffff0000  0x   3b8b0       2     2 
0xffffff4000  0x   3f8b0       0     0 
Distribution for array "b" 
address       byte index  thread  proc 
-------------------------------------- 
0xfffff6e0e8  0x       0       0     0 
0xfffff70000  0x    1f18       4     4 
0xfffff74000  0x    5f18       2     2 
0xfffff78000  0x    9f18       0     0 
0xfffff7c000  0x    df18       4     4 
0xfffff80000  0x   11f18       2     2 
0xfffff84000  0x   15f18       0     0 
0xfffff88000  0x   19f18       4     4 
0xfffff8c000  0x   1df18       2     2 
0xfffff90000  0x   21f18       0     0 
0xfffff94000  0x   25f18       4     4 
0xfffff98000  0x   29f18       2     2 
0xfffff9c000  0x   2df18       0     0 
0xfffffa0000  0x   31f18       4     4 
0xfffffa4000  0x   35f18       2     2 
0xfffffa8000  0x   39f18       0     0 
0xfffffac000  0x   3df18       4     4 

Round-robin placement really did spread the data over all three nodes used by this 5-CPU run.

The placement file shown in Figure 8-14 has one limitation: it is not scalable. That is, if you want to change the number of threads, you need to change the placement file. A placement file can refer to environment variables, making it scalable. A file of this type is shown in Example 8-27.

# scalable placement_file 
memories $NP in topology none # set up memories which are close 
threads $NP                      # number of threads 
# run the last thread on the first memory etc. 
distribute threads $NP-1:0:-1 across memories 

Instead of writing the number of threads into the file, the environment variable $NP is used. (There is no significance to the name; any variable could be used.) When the placement file is read, the reference to the environment variable is replaced by its value taken from the shell. Thus, if NP is set to 2, the first two statements have the same effect as the version in Figure 8-14. To make the last two statements scalable, they are replaced with a distribute statement. This one statement specifies a mapping of all threads to all memories. Its use is covered under “Assigning Threads to Memories”.

You can use simple arithmetic on environment variables in placement files. A scalable placement file that uses a more conventional placement of two threads per memory is shown in Example 8-28. This example uses the MP_SET_NUMTHREADS environment variable used by the MP library (however, it is not recommended to apply dplace to MP library programs). A comparable variable is PT_ITC, the number of processes used by the POSIX Threads library to run a pthreaded program.

# Standard 2 thread per memory scalable placement_file 
memories ($MP_SET_NUMTHREADS + 1)/2 in topology none 
threads $MP_SET_NUMTHREADS 
distribute threads across memories 

The memories statement specifies the number of memories (SN0 nodes) to be allocated to the program and, optionally, their topology. The topologies supported are these:

The number of nodes in a hypercube is a power of two. As a result, cube topology can have strange side effects when used on machines that do not have a power of two number of nodes. In general, use the none (cluster) topology.

The memories are identified by number, starting from 0, in the order that they are allocated. Using cube or cluster topology, there is no way to relate a memory number to a specific SN0 node; memory 0 could be identified with any node in the system on any run of the program.

The threads statement specifies the number of IRIX processes the program will use. The MP library uses IRIX processes as does MPI, PVM, and standard UNIX software. POSIX threads programs are implemented using some number of processes; the number is available within the program using the pthread_get_concurrency() function.

The threads statement does not create the processes; creation happens as a result of the program's own actions. The threads statement tells dplace in advance how many threads the program will create. Threads are identified by number, starting from 0, in the sequence in which the program creates them. Thus the first process, the parent process of the program, is thread 0 to dplace.

There are two statements used to assign threads to memories, run and distribute. The run statement tells dplace to assign a specific thread, by number, to a specific memory by number. Optionally, it can assign the thread to a specific CPU in that node, although there is no particular advantage in doing so. The run statement deals with only one thread, so a placement file that uses it cannot be scalable.

The distribute statement specifies a mapping of all threads to all memories. In Example 8-27, a distribute statement is used to map threads in reverse order to memories, as follows:

distribute threads $NP-1:0:-1 across memories 

The three-part expression, reminiscent of a Fortran DO loop, specifies the threads in the program as first:last:stride. In this example, the last thread ($NP-1) is mapped to memory 0, the next-to-last thread is mapped to memory 1, and so on.

This example maps threads to memories in reverse order. If the normal order were desired, the statement could be written in any of the ways shown in Example 8-29.

# Explicitly indicate order of threads 
distribute threads 0:$NP-1:1 across memories 
# No first:last:stride means 0:last:1 
distribute threads across memories 
# The stride can be omitted when it is 1 
distribute threads 0:$NP-1 across memories 
# first:last:stride notation applies to memories too
distribute threads across memories 0:$NP-1:1 
distribute threads $NP-1:0:-1 across memories $NP-1:0:-1 

The distribute statement also supports a clause telling how threads are assigned to memories, when there are more threads than memories. The clause can be either block or cyclic. In a block mapping with a block size of two, the first two threads are assigned to the first memory, the second two to the second, and so on. In a cyclic mapping, threads are dealt to memories as you deal a deck of cards. The default is to use a block mapping, and the default block size is ceiling (number of threads ÷ number of memories); that is, 2 if there are twice as many threads as memories, and 1 if the number of threads and memories are equal.

The MPI library (version 3.1 or higher) is aware of, and compatible with, the MP library. As a result, you rarely need to use dplace with MPI programs. However, in the rare cases when you do want this combination, the key to success is to apply the mpirun command to dplace, not the reverse. In addition, you must set an environment variable, MPI_DSM_OFF, as in the following two-line example:

% setenv MPI_DSM_OFF 
% mpirun -np 4 dplace -place pfile a.out

This example starts copies of dplace on four CPUs. Each of them, in turn, invokes the program a.out, applying the placement file pfile.

When using a tool such as dplace that produces output on the standard output or standard error streams, it can be difficult to capture the tool's output. The correct approach, as shown in Example 8-31, is to write a small script to do the redirection, and to use mpirun to launch the script.

> cat myscript 
#!/bin/sh 
setenv MPI_DSM_OFF 
dplace -verbose a.out 2> outfile 
> mpirun -np 4 myscript 
hello world from process 0 
hello world from process 1 
hello world from process 2 
hello world from process 3 
> cat outfile 
there are now 1 threads 
Setting up policies and initial thread. 
Migration is off. 
Data placement policy is PlacementDefault. 
Creating data PM. 
Data pagesize is 16k. 
Setting data PM. 
Creating stack PM. 
Stack pagesize is 16k. 
Stack placement policy is PlacementDefault. 
Setting stack PM. 
there are now 2 threads 
there are now 3 threads 
there are now 4 threads 
there are now 5 threads