High Performance Computing 5.2
High Performance Computing 5.2
By:
Charles Severance
Kevin Dowd
High Performance Computing
By:
Charles Severance
Kevin Dowd
Online:
< http://cnx.org/content/col11136/1.5/ >
CONNEXIONS
1
2
and the ability to print locally can make this book available in any country and any school in the world.
Like Wikipedia, those of us who use the book can become the volunteers who will help improve the book
and become co-authors of the book.
I need to thank Kevin Dowd who wrote the rst edition and graciously let me alter it from cover to cover
in the second edition. Mike Loukides of O'Reilly was the editor of both the rst and second editions and we
talk from time to time about a possible future edition of the book. Mike was also instrumental in helping
to release the book from O'Reilly under Creative Commons Attribution. The team at Connexions has been
wonderful to work with. We share a passion for High Performance Computing and new forms of publishing
so that the knowledge reaches as many people as possible. I want to thank Jan Odegard and Kathi Fletcher
for encouraging, supporting and helping me through the re-publishing process. Daniel Williamson did an
amazing job of converting the materials from the O'Reilly formats to the Connexions formats.
I truly look forward to seeing how far this book will go now that we can have an unlimited number of
co-authors to invest and then use the book. I look forward to work with you all.
Charles Severance - November 12, 2009
3
4
High performance RISC processors are designed to be easily inserted into a multiple-processor system
with 2 to 64 CPUs accessing a single memory using symmetric multi processing (SMP). Programming
multiple processors to solve a single problem adds its own set of additional challenges for the programmer.
The programmer must be aware of how multiple processors operate together, and how work can be eciently
divided among those processors.
Even though each processor is very powerful, and small numbers of processors can be put into a single
enclosure, often there will be applications that are so large they need to span multiple enclosures. In order to
cooperate to solve the larger application, these enclosures are linked with a high-speed network to function
as a network of workstations (NOW). A NOW can be used individually through a batch queuing system or
can be used as a large multicomputer using a message passing tool such as parallel virtual machine
(PVM)
or message-passing interface (MPI).
For the largest problems with more data interactions and those users with compute budgets in the millions
of dollars, there is still the top end of the high performance computing spectrum, the scalable parallel
processing systems with hundreds to thousands of processors. These systems come in two avors. One type
is programmed using message passing. Instead of using a standard local area network, these systems are
connected using a proprietary, scalable, high-bandwidth, low-latency interconnect (how is that for marketing
speak?). Because of the high performance interconnect, these systems can scale to the thousands of processors
while keeping the time spent (wasted) performing overhead communications to a minimum.
The second type of large parallel processing system is the scalable non-uniform memory access (NUMA)
systems. These systems also use a high performance inter-connect to connect the processors, but instead of
exchanging messages, these systems use the interconnect to implement a distributed shared memory that can
be accessed from any processor using a load/store paradigm. This is similar to programming SMP systems
except that some areas of memory have slower access than others.
Measuring Performance
When a computer is being purchased for computationally intensive applications, it is important to determine
how well the system will actually perform this function. One way to choose among a set of competing systems
is to have each vendor loan you a system for a period of time to test your applications. At the end of the
evaluation period, you could send back the systems that did not make the grade and pay for your favorite
system. Unfortunately, most vendors won't lend you a system for such an extended period of time unless
there is some assurance you will eventually purchase the system.
More often we evaluate the system's potential performance using benchmarks
. There are industry bench-
marks and your own locally developed benchmarks. Both types of benchmarks require some careful thought
and planning for them to be an eective tool in determining the best system for your application.
• A basic understanding of modern computer architecture. You don't need an advanced degree in
computer engineering, but you do need to understand the basic terminology.
• A basic understanding of benchmarking, or performance measurement, so you can quantify your own
successes and failures and use that information to improve the performance of your application.
This book is intended to be an easily understood introduction and overview of high performance computing.
It is an interesting eld, and one that will become more important as we make even greater demands on
our most common personal computers. In the high performance computer eld, there is always a tradeo
between the single CPU performance and the performance of a multiple processor system. Multiple processor
systems are generally more expensive and dicult to program (unless you have this book).
Some people claim we eventually will have single CPUs so fast we won't need to understand any type of
advanced architectures that require some skill to program.
So far in this eld of computing, even as performance of a single inexpensive microprocessor has increased
over a thousandfold, there seems to be no less interest in lashing a thousand of these processors together to
get a millionfold increase in power. The cheaper the building blocks of high performance computing become,
the greater the benet for using many processors. If at some point in the future, we have a single processor
that is faster than any of the 512-processor scalable systems of today, think how much we could do when we
connect 512 of those new processors together in a single system.
That's what this book is all about. If you're interested, read on.
1.1 Memory
1.1.1 Introduction1
1.1.1.1 Memory
Let's say that you are fast asleep some night and begin dreaming. In your dream, you have a time machine
and a few 500-MHz four-way superscalar processors. You turn the time machine back to 1981. Once you
arrive back in time, you go out and purchase an IBM PC with an Intel 8088 microprocessor running at 4.77
MHz. For much of the rest of the night, you toss and turn as you try to adapt the 500-MHz processor to the
Intel 8088 socket using a soldering iron and Swiss Army knife. Just before you wake up, the new computer
nally works, and you turn it on to run the Linpack2 benchmark and issue a press release. Would you expect
this to turn out to be a dream or a nightmare? Chances are good that it would turn out to be a nightmare,
just like the previous night where you went back to the Middle Ages and put a jet engine on a horse. (You
have got to stop eating double pepperoni pizzas so late at night.)
Even if you can speed up the computational aspects of a processor innitely fast, you still must load and
store the data and instructions to and from a memory. Today's processors continue to creep ever closer to
innitely fast processing. Memory performance is increasing at a much slower rate (it will take longer for
memory to become innitely fast). Many of the interesting problems in high performance computing use a
large amount of memory. As computers are getting faster, the size of problems they tend to operate on also
goes up. The trouble is that when you want to solve these problems at high speeds, you need a memory
system that is large, yet at the same time fasta big challenge. Possible approaches include the following:
• Every memory system component can be made individually fast enough to respond to every memory
access request.
• Slow memory can be accessed in a round-robin fashion (hopefully) to give the eect of a faster memory
system.
• The memory system design can be made wide so that each transfer contains many bytes of informa-
tion.
• The system can be divided into faster and slower portions and arranged so that the fast portion is used
more often than the slow one.
Again, economics are the dominant force in the computer business. A cheap, statistically optimized memory
system will be a better seller than a prohibitively expensive, blazingly fast one, so the rst choice is not much
of a choice at all. But these choices, used in combination, can attain a good fraction of the performance
1 This content is available online at <http://cnx.org/content/m32733/1.3/>.
2 See Chapter 15, Using Published Benchmarks, for details on the Linpack benchmark.
7
8 CHAPTER 1. MODERN COMPUTER ARCHITECTURES
you would get if every component were fast. Chances are very good that your high performance workstation
incorporates several or all of them.
Once the memory system has been decided upon, there are things we can do in software to see that it
is used eciently. A compiler that has some knowledge of the way memory is arranged and the details of
the caches can optimize their use to some extent. The other place for optimizations is in user applications,
as we'll see later in the book. A good pattern of memory access will work with, rather than against, the
components of the system.
In this chapter we discuss how the pieces of a memory system work. We look at how patterns of data
and instruction access factor into your overall runtime, especially as CPU speeds increase. We also talk a
bit about the performance implications of running in a virtual memory environment.
AT models were introduced in the mid-1980s with CPUs that clocked more quickly than the access times
of available commodity memory. Faster memory was available for a price, but vendors punted by selling
computers with wait statesadded to the memory access cycle. Wait states are articial delays that slow
down references so that memory appears to match the speed of a faster CPU at a penalty. However, the
technique of adding wait states begins to signicantly impact performance around 25?33MHz. Today, CPU
speeds are even farther ahead of DRAM speeds.
The clock time for commodity home computers has gone from 210 ns for the XT to around 3 ns for a
300-MHz Pentium-II, but the access time for commodity DRAM has decreased disproportionately less
from 200 ns to around 50 ns. Processor performance doubles every 18 months, while memory performance
doubles roughly every seven years.
The CPU/memory speed gap is even larger in workstations. Some models clock at intervals as short as
1.6 ns. How do vendors make up the dierence between CPU speeds and memory speeds? The memory in
the Cray-1 supercomputer used SRAM that was capable of keeping up with the 12.5-ns clock cycle. Using
SRAM for its main memory system was one of the reasons that most Cray systems needed liquid cooling.
Unfortunately, it's not practical for a moderately priced system to rely exclusively on SRAM for storage.
It's also not practical to manufacture inexpensive systems with enough storage using exclusively SRAM.
The solution is a hierarchy of memories using processor registers, one to three levels of SRAM cache,
DRAM main memory, and virtual memory stored on media such as disk. At each point in the memory
hierarchy, tricks are employed to make the best use of the available technology. For the remainder of this
chapter, we will examine the memory hierarchy and its impact on performance.
In a sense, with today's high performance microprocessor performing computations so quickly, the task
of the high performance programmer becomes the careful management of the memory hierarchy. In some
sense it's a useful intellectual exercise to view the simple computations such as addition and multiplication
as innitely fast in order to get the programmer to focus on the impact of memory operations on the overall
performance of the program.
1.1.3 Registers5
At least the top layer of the memory hierarchy, the CPU registers, operate as fast as the rest of the processor.
The goal is to keep operands in the registers as much as possible. This is especially important for intermediate
values used in a long computation such as:
X = G * 2.41 + A / W - W * M
While computing the value of A divided by W, we must store the result of multiplying G by 2.41. It would
be a shame to have to store this intermediate result in memory and then reload it a few instructions later.
On any modern processor with moderate optimization, the intermediate result is stored in a register. Also,
the value W is used in two computations, and so it can be loaded once and used twice to eliminate a wasted
load.
Compilers have been very good at detecting these types of optimizations and eciently making use of
the available registers since the 1970s. Adding more registers to the processor has some performance benet.
It's not practical to add enough registers to the processor to store the entire problem data. So we must still
use the slower memory technology.
1.1.4 Caches6
Once we go beyond the registers in the memory hierarchy, we encounter caches. Caches are small amounts
of SRAM that store a subset of the contents of the memory. The hope is that the cache will have the right
subset of main memory at the right time.
5 This content is available online at <http://cnx.org/content/m32681/1.3/>.
6 This content is available online at <http://cnx.org/content/m32725/1.3/>.
Registers 2 ns
L1 On-Chip 4 ns
L2 On-Chip 5 ns
L3 O-Chip 30 ns
Memory 220 ns
When every reference can be found in a cache, you say that you have a 100% hit rate. Generally, a hit
rate of 90% or better is considered good for a level-one (L1) cache. In level-two (L2) cache, a hit rate of
above 50% is considered acceptable. Below that, application performance can drop o steeply.
One can characterize the average read performance of the memory hierarchy by examining the probability
that a particular load will be satised at a particular level of the hierarchy. For example, assume a memory
architecture with an L1 cache speed of 10 ns, L2 speed of 30 ns, and memory speed of 300 ns. If a memory
reference were satised from L1 cache 75% of the time, L2 cache 20% of the time, and main memory 5% of
the time, the average memory performance would be:
You can easily see why it's important to have an L1 cache hit rate of 90% or higher.
Given that a cache holds only a subset of the main memory at any time, it's important to keep an index of
which areas of the main memory are currently stored in the cache. To reduce the amount of space that must
be dedicated to tracking which memory areas are in cache, the cache is divided into a number of equal sized
slots known as lines . Each line contains some number of sequential main memory locations, generally four
to sixteen integers or real numbers. Whereas the data within a line comes from the same part of memory,
other lines can contain data that is far separated within your program, or perhaps data from somebody
else's program, as in Figure 1.1 (Cache lines can come from dierent parts of memory). When you ask for
something from memory, the computer checks to see if the data is available within one of these cache lines.
If it is, the data is returned with a minimal delay. If it's not, your program may be delayed while a new line
is fetched from main memory. Of course, if a new line is brought in, another has to be thrown out. If you're
lucky, it won't be the one containing the data you are just about to need.
Figure 1.1
On multiprocessors (computers with several CPUs), written data must be returned to main memory so
the rest of the processors can see it, or all other processors must be made aware of local cache activity.
Perhaps they need to be told to invalidate old lines containing the previous value of the written variable so
that they don't accidentally use stale data. This is known as maintaining coherency
between the dierent
caches. The problem can become very complex in a multiprocessor system.7
Caches are eective because programs often exhibit characteristics that help kep the hit rate high. These
characteristics are called spatial
and temporal locality of reference
; programs often make use of instructions
and data that are near to other instructions and data, both in space and time. When a cache line is
retrieved from main memory, it contains not only the information that caused the cache miss, but also some
neighboring information. Chances are good that the next time your program needs data, it will be in the
cache line just fetched or another one recently fetched.
Caches work best when a program is reading sequentially through the memory. Assume a program is
reading 32-bit integers with a cache line size of 256 bits. When the program references the rst word in
the cache line, it waits while the cache line is loaded from main memory. Then the next seven references to
memory are satised quickly from the cache. This is called unit stride
because the address of each successive
data element is incremented by one and all the data retrieved into the cache is used. The following loop is
a unit-stride loop:
DO I=1,1000000
SUM = SUM + A(I)
END DO
When a program accesses a large data structure using non-unit stride, performance suers because data is
loaded into cache that is not used. For example:
7 Section 3.2.1 describes cache coherency in more detail.
DO I=1,1000000, 8
SUM = SUM + A(I)
END DO
This code would experience the same number of cache misses as the previous loop, and the same amount of
data would be loaded into the cache. However, the program needs only one of the eight 32-bit words loaded
into cache. Even though this program performs one-eighth the additions of the previous loop, its elapsed
time is roughly the same as the previous loop because the memory operations dominate performance.
While this example may seem a bit contrived, there are several situations in which non-unit strides occur
quite often. First, when a FORTRAN two-dimensional array is stored in memory, successive elements in the
rst column are stored sequentially followed by the elements of the second column. If the array is processed
with the row iteration as the inner loop, it produces a unit-stride reference pattern as follows:
REAL*4 A(200,200)
DO J = 1,200
DO I = 1,200
SUM = SUM + A(I,J)
END DO
END DO
Interestingly, a FORTRAN programmer would most likely write the loop (in alphabetical order) as follows,
producing a non-unit stride of 800 bytes between successive load operations:
REAL*4 A(200,200)
DO I = 1,200
DO J = 1,200
SUM = SUM + A(I,J)
END DO
END DO
Because of this, some compilers can detect this suboptimal loop order and reverse the order of the loops to
make best use of the memory system. As we will see in Section 1.2.1, however, this code transformation may
produce dierent results, and so you may have to give the compiler permission to interchange these loops
in this particular example (or, after reading this book, you could just code it properly in the rst place).
while ( ptr != NULL ) ptr = ptr->next;
The next element that is retrieved is based on the contents of the current element. This type of loop bounces
all around memory in no particular pattern. This is called pointer chasing
and there are no good ways to
improve the performance of this code.
A third pattern often found in certain types of codes is called gather (or scatter) and occurs in loops
such as:
Figure 1.2
The arrays A and B both take up exactly 4 KB of storage, and their inclusion together in COMMON assures
that the arrays start exactly 4 KB apart in memory. In a 4-KB direct mapped cache, the same line that is
used for A(1) is used for B(1), and likewise for A(2) and B(2), etc., so alternating references cause repeated
cache misses. To x it, you could either adjust the size of the array A, or put some other variables into
COMMON, between them. For this reason one should generally avoid array dimensions that are close to powers
of two.
When the processor goes looking for a piece of data, the cache lines are asked all at once whether any of
them has it. The cache line containing the data holds up its hand and says I have it; if none of them do,
there is a cache miss. It then becomes a question of which cache line will be replaced with the new data.
Rather than map memory locations to cache lines via an algorithm, like a direct- mapped cache, the memory
system can ask the fully associative cache lines to choose among themselves which memory locations they
will represent. Usually the least recently used line is the one that gets overwritten with new data. The
assumption is that if the data hasn't been used in quite a while, it is least likely to be used in the future.
Fully associative caches have superior utilization when compared to direct mapped caches. It's dicult
to nd real-world examples of programs that will cause thrashing in a fully associative cache. The expense
of fully associative caches is very high, in terms of size, price, and speed. The associative caches that do
exist tend to be small.
Like the previous cache thrasher program, this forces repeated accesses to the same cache lines, except that
now there are three variables contending for the choose set same mapping instead of two. Again, the way
to x it would be to change the size of the arrays or insert something in between them, in COMMON. By the
way, if you accidentally arranged a program to thrash like this, it would be hard for you to detect it aside
from a feeling that the program runs a little slow. Few vendors provide tools for measuring cache misses.
Figure 1.3
four-way set-associative L1 caches for instruction and data and a combined L2 cache.
Figure 1.4
The operating system stores the page-table addresses virtually, so it's going to take a virtual-to-physical
translation to locate the table in memory. One more virtual-to- physical translation, and we nally have
the true address of location 1000. The memory reference can complete, and the processor can return to
executing your program.
easiest case to construct is one where every memory reference your program makes causes a TLB miss:
REAL X(10000000)
COMMON X
DO I=0,9999
DO J=1,10000000,10000
SUM = SUM + X(J+I)
END DO
END DO
Assume that the TLB page size for your computer is less than 40 KB. Every time through the inner loop
in the above example code, the program asks for data that is 4 bytes*10,000 = 40,000 bytes away from the
last reference. That is, each reference falls on a dierent memory page. This causes 1000 TLB misses in the
inner loop, taken 1001 times, for a total of at least one million TLB misses. To add insult to injury, each
reference is guaranteed to cause a data cache miss as well. Admittedly, no one would start with a loop like
the one above. But presuming that the loop was any good to you at all, the restructured version in the code
below would cruise through memory like a warm knife through butter:
REAL X(10000000)
COMMON X
DO I=1,10000000
SUM = SUM + X(I)
END DO
The revised loop has unit stride, and TLB misses occur only every so often. Usually it is not necessary to
explicitly tune programs to make good use of the TLB. Once a program is tuned to be cache-friendly, it
nearly always is tuned to be TLB friendly.
Because there is a performance benet to keeping the TLB very small, the TLB entry often contains a
length eld. A single TLB entry can be over a megabyte in length and can be used to translate addresses
stored in multiple virtual memory pages.
on the memory access patterns. Interestingly, an increase in cache size on the part of vendors can render a
benchmark obsolete.
Simple Memory System
Figure 1.5
Up to 1992, the Linpack 100×100 benchmark was probably the single most- respected benchmark to
determine the average performance across a wide range of applications. In 1992, IBM introduced the IBM
RS-6000 which had a cache large enough to contain the entire 100×100 matrix for the duration of the
benchmark. For the rst time, a workstation had performance on this benchmark on the same order of
supercomputers. In a sense, with the entire data structure in a SRAM cache, the RS-6000 was operating like
a Cray vector supercomputer. The problem was that the Cray could maintain and improve the performance
for a 120×120 matrix, whereas the RS-6000 suered a signicant performance loss at this increased matrix
size. Soon, all the other workstation vendors introduced similarly large caches, and the 100×100 Linpack
benchmark ceased to be useful as an indicator of average application performance.
Figure 1.6
One way to make the cache-line ll operation faster is to widen the memory system as shown in Figure 1.7
(Wide memory system). Instead of having two rows of DRAMs, we create multiple rows of DRAMs. Now
on every 100-ns cycle, we get 32 contiguous bits, and our cache-line lls are four times faster.
Figure 1.7
We can improve the performance of a memory system by increasing the width of the memory system up
to the length of the cache line, at which time we can ll the entire line in a single memory cycle. On the
SGI Power Challenge series of systems, the memory width is 256 bits. The downside of a wider memory
system is that DRAMs must be added in multiples. In many modern workstations and personal computers,
memory is expanded in the form of single inline memory modules (SIMMs). SIMMs currently are either
30-, 72-, or 168-pin modules, each of which is made up of several DRAM chips ready to be installed into a
memory sub-system.
Figure 1.8
Figure 1.9
Dierent access patterns are subject to bank stalls of varying severity. For instance, accesses to every
fourth word in an eight-bank memory system would also be subject to bank stalls, though the recovery would
occur sooner. References to every second word might not experience bank stalls at all; each bank may have
recovered by the time its next reference comes around; it depends on the relative speeds of the processor
and memory system. Irregular access patterns are sure to encounter some bank stalls.
In addition to the bank stall hazard, single-word references made directly to a multibanked memory
system carry a greater latency than those of (successfully) cached memory accesses. This is because references
are going out to memory that is slower than cache, and there may be additional address translation steps
as well. However, banked memory references are pipelined. As long as references are started well enough in
advance, several pipelined, multibanked references can be in ight at one time, giving you good throughput.
The CDC-205 system performed vector operations in a memory-to-memory fashion using a set of explicit
memory pipelines. This system had superior performance for very long unit-stride vector computations. A
single instruction could perform 65,000 computations using three memory pipes.
This is not the actual FORTRAN. Prefetching is usually done in the assembly code generated by the compiler
when it detects that you are stepping through the array using a xed stride. The compiler typically estimate
how far ahead you should be prefetching. In the above example, if the cache-lls were particularly slow, the
value 8 in I+8 could be changed to 16 or 32 while the other values changed accordingly.
In a processor that could only issue one instruction per cycle, there might be no payback to a prefetch
instruction; it would take up valuable time in the instruction stream in exchange for an uncertain benet.
On a superscalar processor, however, a cache hint could be mixed in with the rest of the instruction stream
and issued alongside other, real instructions. If it saved your program from suering extra cache misses, it
would be worth having.
Unlike a vector processor or a prefetch instruction, the post-RISC processor does not need to anticipate
the precise pattern of memory references so it can carefully control the memory subsystem. As a result, the
post-RISC processor can achieve peak performance in a far-wider range of code sequences than either vector
processors or in-order RISC processors with prefetch capability.
This implicit tolerance to memory latency makes the post-RISC processors ideal for use in the scalable
shared-memory processors of the future, where the memory hierarchy will become even more complex than
current processors with three levels of cache and a main memory.
Unfortunately, the one code segment that doesn't benet signicantly from the post-RISC architecture is
the linked-list traversal. This is because the next address is never known until the previous load is completed
so all loads are fundamentally serialized.
Fast page mode DRAM saves time by allowing a mode in which the entire address doesn't have to be re-
clocked into the chip for each memory operation. Instead, there is an assumption that the memory will be
accessed sequentially (as in a cache-line ll), and only the low-order bits of the address are clocked in for
successive reads or writes.
EDO RAM is a modication to output buering on page mode RAM that allows it to operate roughly
twice as quickly for operations other than refresh.
Synchronous DRAM is synchronized using an external clock that allows the cache and the DRAM to
coordinate their operations. Also, SDRAM can pipeline the retrieval of multiple memory bits to improve
overall throughput.
RAMBUS is a proprietary technology capable of 500 MB/sec data transfer. RAMBUS uses signicant
logic within the chip and operates at higher power levels than typical DRAM.
Cached DRAM combines a SRAM cache on the same chip as the DRAM. This tightly couples the
SRAM and DRAM and provides performance similar to SRAM devices with all the limitations of any cache
architecture. One advantage of the CDRAM approach is that the amount of cache is increased as the amount
of DRAM is increased. Also when dealing with memory systems with a large number of interleaves, each
interleave has its own SRAM to reduce latency, assuming the data requested was in the SRAM.
An even more advanced approach is to integrate the processor, SRAM, and DRAM onto a single chip
clocked at say 5 GHz, containing 128 MB of data. Understandably, there is a wide range of technical
problems to solve before this type of component is widely available for $200 but it's not out of the
question. The manufacturing processes for DRAM and processors are already beginning to converge in some
ways (RAMBUS). The biggest performance problem when we have this type of system will be, What to do
if you need 160 MB?
1.1.9 Exercises16
Exercise 1.1.9.1
The following code segment traverses a pointer chain:
while ((p = (char *) *p) != NULL);
How will such a code interact with the cache if all the references fall within a small portion of mem-
ory? How will the code interact with the cache if references are stretched across many megabytes?
Exercise 1.1.9.2
How would the code in Exercise 1.1.9.1 behave on a multibanked memory system that has no
cache?
Exercise 1.1.9.3
A long time ago, people regularly wrote self-modifying code programs that wrote into instruction
memory and changed their own behavior. What would be the implications of self-modifying code
on a machine with a Harvard memory architecture?
Exercise 1.1.9.4
Assume a memory architecture with an L1 cache speed of 10 ns, L2 speed of 30 ns, and memory
speed of 200 ns. Compare the average memory system performance with (1) L1 80%, L2 10%, and
memory 10%; and (2) L1 85% and memory 15%.
Exercise 1.1.9.5
On a computer system, run loops that process arrays of varying length from 16 to 16 million:
ARRAY(I) = ARRAY(I) + 3
How does the number of additions per second change as the array length changes? Experiment
with REAL*4, REAL*8, INTEGER*4, and INTEGER*8.
Which has more signicant impact on performance: larger array elements or integer versus
oating-point? Try this on a range of dierent computers.
Exercise 1.1.9.6
Create a two-dimensional array of 1024×1024. Loop through the array with rows as the inner loop
and then again with columns as the inner loop. Perform a simple operation on each element. Do
the loops perform dierently? Why? Experiment with dierent dimensions for the array and see
the performance impact.
Exercise 1.1.9.7
Write a program that repeatedly executes timed loops of dierent sizes to determine the cache size
for your system.
1.2.2 Reality18
The real world is full of real numbers. Quantities such as distances, velocities, masses, angles, and other
quantities are all real numbers.19 A wonderful property of real numbers is that they have unlimited accuracy.
For example, when considering the ratio of the circumference of a circle to its diameter, we arrive at a value
of 3.141592.... The decimal value for pi
does not terminate. Because real numbers have unlimited accuracy,
pi
even though we can't write it down, is still a real number. Some real numbers are rational numbers because
they can be represented as the ratio of two integers, such as 1/3. Not all real numbers are rational numbers.
Not surprisingly, those real numbers that aren't rational numbers are called irrational. You probably would
not want to start an argument with an irrational number unless you have a lot of free time on your hands.
Unfortunately, on a piece of paper, or in a computer, we don't have enough space to keep writing the
pi
digits of . So what do we do? We decide that we only need so much accuracy and round real numbers to
a certain number of digits. For example, if we decide on four digits of accuracy, our approximation of is pi
3.142. Some state legislature attempted to pass a law that pi
was to be three. While this is often cited as
evidence for the IQ of governmental entities, perhaps the legislature was just suggesting that we only need
pi
one digit of accuracy for . Perhaps they foresaw the need to save precious memory space on computers
when representing real numbers.
1.2.3 Representation20
Given that we cannot perfectly represent real numbers on digital computers, we must come up with a
compromise that allows us to approximate real numbers.21 There are a number of dierent ways that have
been used to represent real numbers. The challenge in selecting a representation is the trade-o between
space and accuracy and the tradeo between speed and accuracy. In the eld of high performance computing
we generally expect our processors to produce a oating- point result every 600-MHz clock cycle. It is pretty
clear that in most applications we aren't willing to drop this by a factor of 100 just for a little more accuracy.
Before we discuss the format used by most high performance computers, we discuss some alternative (albeit
slower) techniques for representing real numbers.
17 This content is available online at <http://cnx.org/content/m32739/1.3/>.
18 This content is available online at <http://cnx.org/content/m32741/1.3/>.
19 In high performance computing we often simulate the real world, so it is somewhat ironic that we use simulated real numbers
(oating-point) in those simulations of the real world.
20 This content is available online at <http://cnx.org/content/m32772/1.3/>.
21 Interestingly, analog computers have an easier time representing real numbers. Imagine a water- adding analog computer
which consists of two glasses of water and an empty glass. The amount of water in the two glasses are perfectly represented
real numbers. By pouring the two glasses into a third, we are adding the two real numbers perfectly (unless we spill some),
and we wind up with a real number amount of water in the third glass. The problem with analog computers is knowing just
how much water is in the glasses when we are all done. It is also problematic to perform 600 million additions per second using
this technique without getting pretty wet. Try to resist the temptation to start an argument over whether quantum mechanics
would cause the real numbers to be rational numbers. And don't point out the fact that even digital computers are really
analog computers at their core. I am trying to keep the focus on oating-point values, and you keep drifting away!
123.45
0001 0010 0011 0100 0101
This format allows the programmer to choose the precision required for each variable. Unfortunately, it is
dicult to build extremely high-speed hardware to perform arithmetic operations on these numbers. Because
each number may be far longer than 32 or 64 bits, they did not t nicely in a register. Much of the oating-
point operations for BCD were done using loops in microcode. Even with the exibility of accuracy on BCD
representation, there was still a need to round real numbers to t into a limited amount of space.
Another limitation of the BCD approach is that we store a value from 09 in a four-bit eld. This eld
is capable of storing values from 015 so some of the space is wasted.
Figure 1.10
The limitation that occurs when using rational numbers to represent real numbers is that the size of the
numerators and denominators tends to grow. For each addition, a common denominator must be found. To
keep the numbers from becoming extremely large, during each operation, it is important to nd the greatest
common divisor (GCD) to reduce fractions to their most compact representation. When the values grow
and there are no common divisors, either the large integer values must be stored using dynamic memory or
some form of approximation must be used, thus losing the primary advantage of rational numbers.
For mathematical packages such as Maple or Mathematica that need to produce exact results on smaller
data sets, the use of rational numbers to represent real numbers is at times a useful technique. The perfor-
mance and storage cost is less signicant than the need to produce exact results in some instances.
1.2.3.4 Mantissa/Exponent
The oating-point format that is most prevalent in high performance computing is a variation on scientic
notation. In scientic notation the real number is represented using a mantissa, base, and exponent: 6.02 ×
1023 .
The mantissa typically has some xed number of places of accuracy. The mantissa can be represented in
base 2, base 16, or BCD. There is generally a limited range of exponents, and the exponent can be expressed
as a power of 2, 10, or 16.
The primary advantage of this representation is that it provides a wide overall range of values while using
a xed-length storage representation. The primary limitation of this format is that the dierence between
two successive values is not uniform. For example, assume that you can represent three base-10 digits,
and your exponent can range from 10 to 10. For numbers close to zero, the distance between successive
numbers is very small. For the number 1.72 × 10−10 , the next larger number is 1.73 × 10−10 . The distance
between these two close small numbers is 0.000000000001. For the number 6.33 × 1010 , the next larger
number is 6.34 × 1010 . The distance between these close large numbers is 100 million.
In Figure 1.11 (Distance between successive oating-point numbers), we use two base-2 digits with an
exponent ranging from 1 to 1.
22 Perhaps banks round this instead of truncating, knowing that they will always make it up in teller machine fees.
Figure 1.11
There are multiple equivalent representations of a number when using scientic notation:
6.00 × 105
0.60 × 106
0.06 × 107
By convention, we shift the mantissa (adjust the exponent) until there is exactly one nonzero digit to
the left of the decimal point. When a number is expressed this way, it is said to be normalized. In the
above list, only 6.00 × 105 is normalized. Figure 1.12 (Normalized oating-point numbers) shows how some
of the oating-point numbers from Figure 1.11 (Distance between successive oating-point numbers) are not
normalized.
While the mantissa/exponent has been the dominant oating-point approach for high performance com-
puting, there were a wide variety of specic formats in use by computer vendors. Historically, each computer
vendor had their own particular format for oating-point numbers. Because of this, a program executed on
several dierent brands of computer would generally produce dierent answers. This invariably led to heated
discussions about which system provided the right answer and which system(s) were generating meaningless
results.23
Figure 1.12
23 Interestingly, there was an easy answer to the question for many programmers. Generally they trusted the results from the
computer they used to debug the code and dismissed the results from other computers as garbage.
When storing oating-point numbers in digital computers, typically the mantissa is normalized, and then
the mantissa and exponent are converted to base-2 and packed into a 32- or 64-bit word. If more bits were
allocated to the exponent, the overall range of the format would be increased, and the number of digits of
accuracy would be decreased. Also the base of the exponent could be base-2 or base-16. Using 16 as the
base for the exponent increases the overall range of exponents, but because normalization must occur on
four-bit boundaries, the available digits of accuracy are reduced on the average. Later we will see how the
IEEE 754 standard for oating-point format represents numbers.
REAL*4 X,Y
X = 0.1
Y = 0
DO I=1,10
Y = Y + X
ENDDO
IF ( Y .EQ. 1.0 ) THEN
PRINT *,'Algebra is truth'
ELSE
PRINT *,'Not here'
ENDIF
PRINT *,1.0-Y
END
At rst glance, this appears simple enough. Mathematics tells us ten times 0.1 should be one. Unfortunately,
because 0.1 cannot be represented exactly as a base-2 decimal, it must be rounded. It ends up being rounded
down to the last bit. When ten of these slightly smaller numbers are added together, it does not quite add
up to 1.0. When X and Y are REAL*4, the dierence is about 10-7 , and when they are REAL*8, the dierence
is about 10-16 .
One possible method for comparing computed values to constants is to subtract the values and test to
see how close the two values become. For example, one can rewrite the test in the above code to be:
X = 1.25E8
Y = X + 7.5E-3
IF ( X.EQ.Y ) THEN
PRINT *,'Am I nuts or what?'
ENDIF
While both of these numbers are precisely representable in oating-point, adding them is problematic. Prior
to adding these numbers together, their decimal points must be aligned as in Figure 1.13 (Figure 4-4: Loss
of accuracy while aligning decimal points).
Figure 1.13
Unfortunately, while we have computed the exact result, it cannot t back into a REAL*4 variable (7
digits of accuracy) without truncating the 0.0075. So after the addition, the value in Y is exactly 1.25E8.
Even sadder, the addition could be performed millions
of times, and the value for Y would still be 1.25E8.
Because of the limitation on precision, not all algebraic laws apply all the time. For instance, the answer
you obtain from X+Y will be the same as Y+X, as per the commutative law for addition. Whichever operand
you pick rst, the operation yields the same result; they are mathematically equivalent. It also means that
you can choose either of the following two forms and get the same answer:
(X + Y) + Z
(Y + X) + Z
(Y + Z) + X
The third version isn't equivalent to the rst two because the order of the calculations has changed. Again, the
rearrangement is equivalent algebraically, but not computationally. By changing the order of the calculations,
we have taken advantage of the associativity of the operations; we have made an associative transformation
of the original code.
To understand why the order of the calculations matters, imagine that your computer can perform
arithmetic signicant to only ve decimal places.
Also assume that the values of X, Y, and Z are .00005, .00005, and 1.0000, respectively. This means that:
but:
The two versions give slightly dierent answers. When adding Y+Z+X, the sum of the smaller numbers was
insignicant when added to the larger number. But when computing X+Y+Z, we add the two small numbers
rst, and their combined sum is large enough to inuence the nal answer. For this reason, compilers
that rearrange operations for the sake of performance generally only do so after the user has requested
optimizations beyond the defaults.
For these reasons, the FORTRAN language is very strict about the exact order of evaluation of ex-
pressions. To be compliant, the compiler must ensure that the operations occur exactly as you express
them.26
26 Often even if you didn't mean it.
a = x + (y + z);
However, the C compiler is free to ignore you, and combine X, Y, and Z in any order it pleases.
Now armed with this knowledge, view the following harmless-looking code segment:
REAL*4 SUM,A(1000000)
SUM = 0.0
DO I=1,1000000
SUM = SUM + A(I)
ENDDO
Begins to look like a nightmare waiting to happen. The accuracy of this sum depends of the relative
magnitudes and order of the values in the array A. If we sort the array from smallest to largest and then
perform the additions, we have a more accurate value. There are other algorithms for computing the sum
of an array that reduce the error without requiring a full sort of the data. Consult a good textbook on
numerical analysis for the details on these algorithms.
If the range of magnitudes of the values in the array is relatively small, the straight- forward computation
of the sum is probably sucient.
Figure 1.14
To perform this computation and round it correctly, we do not need to increase the number of signicant
digits for stored
values. We do, however, need additional digits of precision while performing the computation.
The solution is to add extra guard digits
which are maintained during the interim steps of the compu-
tation. In our case, if we maintained six digits of accuracy while aligning operands, and rounded before
normalizing and assigning the nal value, we would get the proper result. The guard digits only need to be
present as part of the oating-point execution unit in the CPU. It is not necessary to add guard digits to
the registers or to the values stored in memory.
It is not necessary to have an extremely large number of guard digits. At some point, the dierence in
the magnitude between the operands becomes so great that lost digits do not aect the addition or rounding
results.
• Storage formats
• Precise specications of the results of operations
• Special values
• Specied runtime behavior on illegal operations
Specifying the oating-point format to this level of detail insures that when a computer system is compliant
with the standard, users can expect repeatable execution from one hardware platform to another when
operations are executed in the same order.
Table 1.2
In FORTRAN, the 32-bit format is usually called REAL, and the 64-bit format is usually called DOUBLE.
However, some FORTRAN compilers double the sizes for these data types. For that reason, it is safest to
declare your FORTRAN variables as REAL*4 or REAL*8. The double-extended format is not as well supported
in compilers and hardware as the single- and double-precision formats. The bit arrangement for the single
and double formats are shown in Figure 1.15 (IEEE754 oating-point formats).
Based on the storage layouts in Table 1.2: Parameters of IEEE 32- and 64-Bit Formats, we can derive
the ranges and accuracy of these formats, as shown in Table 1.3.
Figure 1.15
Table 1.3 : Range and Accuracy of IEEE 32- and 64-Bit Formats
Figure 1.16
The 64-bit format is similar, except the exponent is 11 bits long, biased by adding 1023 to the exponent,
and the signicand is 54 bits long.
• Addition
• Subtraction
• Multiplication
• Division
• Square root
• Remainder (modulo)
• Conversion to/from integer
• Conversion to/from printed base-10
These operations are specied in a machine-independent manner, giving exibility to the CPU designers to
implement the operations as eciently as possible while maintaining compliance with the standard. During
operations, the IEEE standard requires the maintenance of two guard digits and a sticky bit for intermediate
values. The guard digits above and the sticky bit are used to indicate if any of the bits beyond the second
guard digit is nonzero.
29 This content is available online at <http://cnx.org/content/m32756/1.3/>.
Figure 1.17
In Figure 1.17 (Computation using guard and sticky bits), we have ve bits of normal precision, two
guard digits, and a sticky bit. Guard bits simply operate as normal bits as if the signicand were 25
bits. Guard bits participate in rounding as the extended operands are added. The sticky bit is set to 1 if
any of the bits beyond the guard bits is nonzero in either operand.30 Once the extended sum is computed,
it is rounded so that the value stored in memory is the closest possible value to the extended sum including
the guard digits. Table 1.4 shows all eight possible values of the two guard digits and the sticky bit and the
resulting stored value with an explanation as to why.
30 If you are somewhat hardware-inclined and you think about it for a moment, you will soon come up with a way to properly
maintain the sticky bit without ever computing the full innite precision sum. You just have to keep track as things get
shifted around.
The rst priority is to check the guard digits. Never forget that the sticky bit is just a hint, not a real
digit. So if we can make a decision without looking at the sticky bit, that is good. The only decision we
are making is to round the last storable bit up or down. When that stored value is retrieved for the next
computation, its guard digits are set to zeros. It is sometimes helpful to think of the stored value as having
the guard digits, but set to zero.
Two guard digits and the sticky bit in the IEEE format insures that operations yield the same rounding
as if the intermediate result were computed using unlimited precision and then rounded to t within the
limits of precision of the nal computed value.
At this point, you might be asking, Why do I care about this minutiae? At some level, unless you are a
hardware designer, you don't care. But when you examine details like this, you can be assured of one thing:
when they developed the IEEE oating-point standard, they looked at the details very
carefully. The goal
was to produce the most accurate possible oating-point standard within the constraints of a xed-length
32- or 64-bit format. Because they did such a good job, it's one less thing you have to worry about. Besides,
this stu makes great exam questions.
Table 1.5
The value of the exponent and signicand determines which type of special value this particular oating-
point number represents. Zero is designed such that integer zero and oating-point zero are the same bit
pattern.
Denormalized numbers can occur at some point as a number continues to get smaller, and the exponent
has reached the minimum value. We could declare that minimum to be the smallest representable value.
However, with denormalized values, we can continue by setting the exponent bits to zero and shifting the
signicand bits to the right, rst adding the leading 1 that was dropped, then continuing to add leading
zeros to indicate even smaller values. At some point the last nonzero digit is shifted o to the right, and the
value becomes zero. This approach is called gradual underow
where the value keeps approaching zero and
then eventually becomes zero. Not all implementations support denormalized numbers in hardware; they
might trap to a software routine to handle these numbers at a signicant performance cost.
At the top end of the biased exponent value, an exponent of all 1s can represent the Not a Number
(NaN) value or innity. Innity occurs in computations roughly according to the principles of mathematics.
If you continue to increase the magnitude of a number beyond the range of the oating-point format, once
the range has been exceeded, the value becomes innity. Once a value is innity, further additions won't
increase it, and subtractions won't decrease it. You can also produce the value innity by dividing a nonzero
value by zero. If you divide a nonzero value by innity, you get zero as a result.
The NaN value indicates a number that is not mathematically dened. You can generate a NaN by
dividing zero by zero, dividing innity by innity, or taking the square root of -1. The dierence between
innity and NaN is that the NaN value has a nonzero signicand. The NaN value is very sticky. Any
operation that has a NaN as one of its inputs always produces a NaN result.
• Overow to innity
• Underow to zero
• Division by zero
• Invalid operation
• Inexact operation
According to the standard, these traps are under the control of the user. In most cases, the compiler runtime
library manages these traps under the direction from the user through compiler ags or runtime library calls.
Traps generally have signicant overhead compared to a single oating-point instruction, and if a program
is continually executing trap code, it can signicantly impact performance.
In some cases it's appropriate to ignore traps on certain operations. A commonly ignored trap is the
underow trap. In many iterative programs, it's quite natural for a value to keep reducing to the point where
it disappears. Depending on the application, this may or may not be an error situation so this exception
can be safely ignored.
If you run a program and then it terminates, you see a message such as:
Overflow handler called 10,000,000 times
It probably means that you need to gure out why your code is exceeding the range of the oating-point
format. It probably also means that your code is executing more slowly because it is spending too much
time in its error handlers.
32 This content is available online at <http://cnx.org/content/m32760/1.3/>.
• The compiler is too conservative in trying to generate IEEE-compliant code and produces code that
doesn't operate at the peak speed of the processor. On some processors, to fully support gradual under-
ow, extra instructions must be generated for certain instructions. If your code will never underow,
these instructions are unnecessary overhead.
• The optimizer takes liberties rewriting your code to improve its performance, eliminating some neces-
sary steps. For example, if you have the following code:
Z = X + 500
Y = Z - 200
The optimizer may replace it with Y = X + 300. However, in the case of a value for X that is close to
overow, the two sequences may not produce the same result.
Sometimes a user prefers fast code that loosely conforms to the IEEE standard, and at other times the
user will be writing a numerical library routine and need total control over each oating-point operation.
Compilers have a challenge supporting the needs of both of these types of users. Because of the nature of
the high performance computing market and benchmarks, often the fast and loose approach prevails in
many compilers.
• Look for compiler options that relax or enforce strict IEEE compliance and choose the appropriate
option for your program. You may even want to change these options for dierent portions of your
program.
• Use REAL*8 for computations unless you are sure REAL*4 has sucient precision. Given that REAL*4
has roughly 7 digits of precision, if the bottom digits become meaningless due to rounding and compu-
tations, you are in some danger of seeing the eect of the errors in your results. REAL*8 with 13 digits
makes this much less likely to happen.
• Be aware of the relative magnitude of numbers when you are performing additions.
• When summing up numbers, if there is a wide range, sum from smallest to largest.
• Perform multiplications before divisions whenever possible.
• When performing a comparison with a computed value, check to see if the values are close rather
than identical.
• Make sure that you are not performing any unnecessary type conversions during the critical portions
of your code.
33 This content is available online at <http://cnx.org/content/m32762/1.3/>.
34 This content is available online at <http://cnx.org/content/m32768/1.3/>.
An excellent reference on oating-point issues and the IEEE format is What Every Computer Scientist
Should Know About Floating-Point Arithmetic, written by David Goldberg, in ACM Computing Surveys
magazine (March 1991). This article gives examples of the most common problems with oating-point and
outlines the solutions. It also covers the IEEE oating-point format very thoroughly. I also recommend
you consult Dr. William Kahan's home page (http://www.cs.berkeley.edu/∼wkahan/35 ) for some excellent
materials on the IEEE format and challenges using oating-point arithmetic. Dr. Kahan was one of the
original designers of the Intel i8087 and the IEEE 754 oating-point format.
1.2.13 Exercises36
Exercise 1.2.13.1
Run the following code to count the number of inverses that are not perfectly accurate:
REAL*4 X,Y,Z
INTEGER I
I = 0
DO X=1.0,1000.0,1.0
Y = 1.0 / X
Z = Y * X
IF ( Z .NE. 1.0 ) THEN
I = I + 1
ENDIF
ENDDO
PRINT *,'Found ',I
END
Exercise 1.2.13.2
Change the type of the variables to REAL*8 and repeat. Make sure to keep the optimization at a
suciently low level (-00) to keep the compiler from eliminating the computations.
Exercise 1.2.13.3
Write a program to determine the number of digits of precision for REAL*4 and REAL*8.
Exercise 1.2.13.4
Write a program to demonstrate how summing an array forward to backward and backward to
forward can yield a dierent result.
Exercise 1.2.13.5
Assuming your compiler supports varying levels of IEEE compliance, take a signicant computa-
tional code and test its overall performance under the various IEEE compliance options. Do the
results of the program change?
35 http://www.cs.berkeley.edu/∼wkahan/
36 This content is available online at <http://cnx.org/content/m32765/1.3/>.
47
48 CHAPTER 2. PROGRAMMING AND TUNING SOFTWARE
2.1.2 History of Compilers3
If you have been in high performance computing since its beginning in the 1950s, you have programmed
in several languages during that time. During the 1950s and early 1960s, you programmed in assembly
language. The constraint on memory and slow clock rates made every instruction precious. With small
memories, overall program size was typically small, so assembly language was sucient. Toward the end
of the 1960s, programmers began writing more of their code in a high-level language such as FORTRAN.
Writing in a high-level language made your work much more portable, reliable, and maintainable. Given
the increasing speed and capacity of computers, the cost of using a high-level language was something most
programmers were willing to accept. In the 1970s if a program spent a particularly large amount of time
in a particular routine, or the routine was part of the operating system or it was a commonly used library,
most likely it was written in assembly language.
During the late 1970s and early 1980s, optimizing compilers
continued to improve to the point that all
but the most critical portions of general-purpose programs were written in high-level languages. On the
average, the compilers generate better code than most assembly language programmers. This was often
because a compiler could make better use of hardware resources such as registers. In a processor with 16
registers, a programmer might adopt a convention regarding the use of registers to help keep track of what
value is in what register. A compiler can use each register as much as it likes because it can precisely track
when a register is available for another use.
However, during that time, high performance computer architecture was also evolving. Cray Research
was developing vector processors at the very top end of the computing spectrum. Compilers were not quite
ready to determine when these new vector instructions could be used. Programmers were forced to write
assembly language or create highly hand-tuned FORTRAN that called the appropriate vector routines in
their code. In a sense, vector processors turned back the clock when it came to trusting the compiler for a
while. Programmers never lapsed completely into assembly language, but some of their FORTRAN started
looking rather un-FORTRAN like. As the vector computers matured, their compilers became increasingly
able to detect when vectorization could be performed. At some point, the compilers again became better
than programmers on these architectures. These new compilers reduced the need for extensive directives or
language extensions.4
The RISC revolution led to an increasing dependence on the compiler. Programming early RISC proces-
sors such as the Intel i860 was painful compared to CISC processors. Subtle dierences in the way a program
was coded in machine language could have a signicant impact on the overall performance of the program.
For example, a programmer might have to count the instruction cycles between a load instruction and the
use of the results of the load in a computational instruction. As superscalar processors were developed,
certain pairs of instructions could be issued simultaneously, and others had to be issued serially. Because
there were a large number of dierent RISC processors produced, programmers did not have time to learn
the nuances of wringing the last bit of performance out of each processor. It was much easier to lock the
processor designer and the compiler writer together (hopefully they work for the same company) and have
them hash out the best way to generate the machine code. Then everyone would use the compiler and get
code that made reasonably good use of the hardware.
The compiler became an important tool in the processor design cycle. Processor designers had much
greater exibility in the types of changes they could make. For example, it would be a good design in the
next revision of a processor to execute existing codes 10% slower than a new revision, but by recompiling
the code, it would perform 65% faster. Of course it was important to actually provide that compiler when
the new processor was shipped and have the compiler give that level of performance across a wide range of
codes rather than just one particular benchmark suite.
3 This content is available online at <http://cnx.org/content/m33686/1.3/>.
4 The Livermore Loops was a benchmark that specically tested the capability of a compiler to eectively optimize a set of
loops. In addition to being a performance benchmark, it was also a compiler benchmark.
Figure 2.1
The compilation process is typically broken down into a number of identiable steps, as shown in Fig-
ure 2.1 (Basic compiler processes). While not all compilers are implemented in exactly this way, it helps to
understand the dierent functions a compiler must perform:
1. A precompiler or preprocessor phase is where some simple textual manipulation of the source code
is performed. The preprocessing step can be processing of include les and making simple string
substitutions throughout the code.
2. The lexical analysis phase is where the incoming source statements are decomposed into tokens such
as variables, constants, comments, or language elements.
3. The parsing phase is where the input is checked for syntax, and the compiler translates the incoming
program into an intermediate language that is ready for optimization.
7 This content is available online at <http://cnx.org/content/m33694/1.3/>.
As compilers become more and more sophisticated in order to wring the last bit of performance from the
processor, some of these steps (especially the optimization and code-generation steps) become more and
more blurred. In this chapter, we focus on the traditional optimizing compiler, and in later chapters we will
look more closely at how modern compilers do more sophisticated optimizations.
A = -B + C * D / E
Taken all at once, this statement has four operators and four operands: /, *, +, and - (negate), and B, C, D,
and E. This is clearly too much to t into one quadruple. We need a form with exactly one operator and, at
most, two operands per statement. The recast version that follows manages to do this, employing temporary
variables to hold the intermediate results:
T1 = D / E
T2 = C * T1
T3 = -B
A = T3 + T2
A workable intermediate language would, of course, need some other features, like pointers. We're going to
suggest that we create our own intermediate language to investigate how optimizations work. To begin, we
need to establish a few rules:
• Instructions consist of one opcode, two operands, and a result. Depending on the instruction, the
operands may be empty.
• Assignments are of the form X := Y op Z, meaning X gets the result of op applied to Y and Z.
• All memory references are explicit load from or store to temporaries t . n
• Logical values used in branches are calculated separately from the actual branch.
8 By denitions, we mean the assignment of values: not declarations.
9 More generally, code can be cast as n-tuples. It depends on the level of the intermediate language.
If we were building a compiler, we'd need to be a little more specic. For our purposes, this will do. Consider
the following bit of C code:
while (j < n) {
k = k + j * 2;
m = j * 2;
j++;
}
This loop translates into the intermediate language representation shown here:
A:: t1 := j
t2 := n
t3 := t1 < t2
jmp (B) t3
jmp (C) TRUE
B:: t4 := k
t5 := j
t6 := t5 * 2
t7 := t4 + t6
k := t7
t8 := j
t9 := t8 * 2
m := t9
t10 := j
t11 := t10 + 1
j := t11
jmp (A) TRUE
C::
Each C source line is represented by several IL statements. On many RISC processors, our IL code is so
close to machine language that we could turn it directly into object code.10 Often the lowest optimization
level does a literal translation from the intermediate language to machine code. When this is done, the code
generally is very large and performs very poorly. Looking at it, you can see places to save a few instructions.
For instance, j gets loaded into temporaries in four places; surely we can reduce that. We have to do some
analysis and make some optimizations.
way, each basic block has one entrance (at the top) and one exit (at the bottom). Figure 2.2 (Intermediate
language divided into basic blocks) represents our IL code as a group of three basic blocks. Basic blocks
make code easier to analyze. By restricting ow of control within a basic block from top to bottom and
eliminating all the branches, we can be sure that if the rst statement gets executed, the second one does
too, and so on. Of course, the branches haven't disappeared, but we have forced them outside the blocks in
the form of the connecting arrows the ow graph
.
Figure 2.2
We are now free to extract information from the blocks themselves. For instance, we can say with
certainty which variables a given block uses and which variables it denes (sets the value of ). We might
not be able to do that if the block contained a branch. We can also gather the same kind of information
about the calculations it performs. After we have analyzed the blocks so that we know what goes in and
what comes out, we can modify them to improve performance and just worry about the interaction between
blocks.
• No optimization : Generates machine code directly from the intermediate language, which can be very
large and slow code. The primary uses of no optimization are for debuggers and establishing the correct
program output. Because every operation is done precisely as the user specied, it must be right.
• Basic optimizations :Similar to those described in this chapter. They generally work to minimize the
intermediate language and generate fast compact code.
• Interprocedural analysis : Looks beyond the boundaries of a single routine for optimization opportu-
nities. This optimization level might include extending a basic optimization such as copy propagation
across multiple routines. Another result of this technique is procedure inlining where it will improve
performance.
• Runtime prole analysis : It is possible to use runtime proles to help the compiler generate improved
code based on its knowledge of the patterns of runtime execution gathered from prole information.
• Floating-point optimizations : The IEEE oating-point standard (IEEE 754) species precisely how
oating- point operations are performed and the precise side eects of these operations. The compiler
may identify certain algebraic transformations that increase the speed of the program (such as replac-
ing a division with a reciprocal and a multiplication) but might change the output results from the
unoptimized code.
• Data ow analysis : Identies potential parallelism between instructions, blocks, or even successive loop
iterations.
• Advanced optimization : May include automatic vectorization, parallelization, or data decomposition
on advanced architecture computers.
These optimizations might be controlled by several dierent compiler options. It often takes some time
to gure out the best combination of compiler ags for a particular code or set of codes. In some cases,
programmers compile dierent routines using dierent optimization settings for best overall performance.
X = Y
12 This content is available online at <http://cnx.org/content/m33696/1.3/>.
Z = 1.0 + X
As written, the second statement requires the results of the rst before it can proceed you need X
to calculate Z. Unnecessary dependencies could translate into a delay at runtime.13 With a little bit of
rearrangement we can make the second statement independent of the rst, by propagating
a copy of Y. The
new calculation for Z uses the value of Y directly:
X = Y
Z = 1.0 + Y
Notice that we left the rst statement, X=Y, intact. You may ask, Why keep it? The problem is that we
can't tell whether the value of X is needed elsewhere. That is something for another analysis to decide. If it
turns out that no other statement needs the new value of X, the assignment is eliminated later by dead code
removal.
PROGRAM MAIN
INTEGER I,K
PARAMETER (I = 100)
K = 200
J = I + K
END
Because I and K are constant individually, the combination I+K is constant, which means that J is a constant
too. The compiler reduces constant expressions like I+K into constants with a technique called constant
folding.
How does constant folding work? You can see that it is possible to examine every path along which a
given variable could be dened en route to a particular basic block. If you discover that all paths lead back
to the same value, that is a constant; you can replace all references to that variable with that constant. This
replacement has a ripple-through eect. If the compiler nds itself looking at an expression that is made
up solely of constants, it can evaluate the expression at compile time and replace it with a constant. After
several iterations, the compiler will have located most of the expressions that are candidates for constant
folding.
A programmer can sometimes improve performance by making the compiler aware of the constant values
in your application. For example, in the following code segment:
X = X * Y
13 This code is an example of a ow dependence. I describe dependencies in detail in Section 3.1.1.
DO I = 1,10000
DO J=1,IDIM
.....
ENDDO
ENDDO
After looking at the code, it's clear that IDIM was either 1, 2, or 3, depending on the data set in use. Clearly
if the compiler knew that IDIM was 1, it could generate much simpler and faster code.
You can easily write some unreachable code into a program by directing the ow of control around it
permanently. If the compiler can tell it's unreachable, it will eliminate it. For example, it's impossible to
reach the statement I = 4 in this program:
PROGRAM MAIN
I = 2
WRITE (*,*) I
STOP
I = 4
WRITE (*,*) I
END
The compiler throws out everything after the STOP statement and probably gives you a warning. Unreachable
code produced by the compiler during optimization will be quietly whisked away.
Computations with local variables can produce results that are never used. By analyzing a variable's
denitions and uses, the compiler can see whether any other part of the routine references it. Of course the
compiler can't tell the ultimate fate of variables that are passed between routines, external or common, so
those computations are always kept (as long as they are reachable).14 In the following program, computations
involving k contribute nothing to the nal answer and are good candidates for dead code elimination:
14 If a compiler does sucient interprocedural analysis, it can even optimize variables across routine boundaries. Interproce-
dural analysis can be the bane of benchmark codes trying to time a computation without using the results of the computation.
main ()
{
int i,k;
i = k = 1;
i += 1;
k += 2;
printf ("%d\n",i);
}
Dead code elimination has often produced some amazing benchmark results from poorly written benchmarks.
See for an example of this type of code.
REAL X,Y
Y = X**2
J = K*2
For the exponentiation operation on the rst line, the compiler generally makes an embedded mathematical
subroutine library call. In the library routine, X is converted to a logarithm, multiplied, then converted back.
Overall, raising X to a power is expensive taking perhaps hundreds of machine cycles. The key is to notice
that X is being raised to a small integer power. A much cheaper alternative would be to express it as X*X,
and pay only the cost of multiplication. The second statement shows integer multiplication of a variable K
by 2. Adding K+K yields the same answer, but takes less time.
There are many opportunities for compiler-generated strength reductions; these are just a couple of them.
We will see an important special case when we look at induction variable simplication. Another example
of a strength reduction is replacing multiplications by integer powers of two by logical shifts.
x = y * z;
q = r + x + x;
When the compiler recognizes that a variable is being recycled, and that its current and former uses are
independent, it can substitute a new variable to keep the calculations separate:
x0 = y * z;
q = r + x0 + x0;
x = a + b;
Variable renaming is an important technique because it claries that calculations are independent of each
other, which increases the number of things that can be done in parallel.
D = C * (A + B)
E = (A + B)/2.
Rather than calculate A + B twice, the compiler can generate a temporary variable and use it wherever A +
B is required:
temp = A + B
D = C * temp
E = temp/2.
Dierent compilers go to dierent lengths to nd common subexpressions. Most pairs, such as A+B, are
recognized. Some can recognize reuse of intrinsics, such as SIN(X). Don't expect the compiler to go too far
though. Subexpressions like A+B+C are not computationally equivalent to reassociated forms like B+C+A, even
though they are algebraically the same. In order to provide predictable results on computations, FORTRAN
must either perform operations in the order specied by the user or reorder them in a way to guarantee
exactly the same result. Sometimes the user doesn't care which way A+B+C associates, but the compiler
cannot assume the user does not care.
Address calculations provide a particularly rich opportunity for common subexpression elimination. You
don't see the calculations in the source code; they're generated by the compiler. For instance, a reference to
an array element A(I,J) may translate into an intermediate language expression such as:
address(A) + (I-1)*sizeof_datatype(A)
+ (J-1)*sizeof_datatype(A) * column_dimension(A)
If A(I,J) is used more than once, we have multiple copies of the same address computation. Common
subexpression elimination will (hopefully) discover the redundant computations and group them together.
DO I=1,N
A(I) = B(I) + C * D
E = G(K)
ENDDO
Below, we have modied the expressions to show how they can be moved to the outside:
temp = C * D
DO I=1,N
A(I) = B(I) + temp
ENDDO
E = G(K)
It is possible to move code before or after the loop body. As with common subexpression elimination, address
arithmetic is a particularly important target for loop- invariant code motion. Slowly changing portions of
index calculations can be pushed into the suburbs, to be executed only when needed.
DO I=1,N
K = I*4 + M
...
ENDDO
Induction variable simplication replaces calculations for variables like K with simpler ones. Given a starting
point and the expression's rst derivative, you can arrive at K's value for the nth iteration by stepping
through the n-1 intervening iterations:
K = M
DO I=1,N
K = K + 4
...
ENDDO
The two forms of the loop aren't equivalent; the second won't give you the value of K given any value of I.
n
Because you can't jump into the middle of the loop on the th iteration, K always takes on the same values
it would have if we had kept the original expression.
Induction variable simplication probably wouldn't be a very important optimization, except that array
address calculations look very much like the calculation for K in the example above. For instance, the address
calculation for A(I) within a loop iterating on the variable I looks like this:
Performing all that math is unnecessary. The compiler can create a new induction variable for references to
A and simplify the address calculations:
Induction variable simplication is especially useful on processors that can automatically increment a register
each time it is used as a pointer for a memory reference. While stepping through a loop, the memory reference
and the address arithmetic can both be squeezed into a single instructiona great savings.
Some instructions in the repertoire also save your compiler from having to issue others. Examples are
auto-increment for registers being used as array indices or conditional assignments in lieu of branches. These
both save the processor from extra calculations and make the instruction stream more compact.
Lastly, there are opportunities for increased parallelism. Programmers generally think serially, specifying
steps in logical succession. Unfortunately, serial source code makes serial object code. A compiler that hopes
to eciently use the parallelism of the processor will have to be able to move instructions around and nd
operations that can be issued side by side. This is one of the biggest challenges for compiler writers today. As
superscalar and very long instruction word(VLIW) designs become capable of executing more instructions
per clock cycle, the compiler will have to dig deeper for operations that can execute at the same time.
2.1.8 Exercises16
Exercise 2.1.8.1
Does your compiler recognize dead code in the program below? How can you be sure? Does the
compiler give you a warning?
main()
{
int k=1;
if (k == 0)
printf ("This statement is never executed.\n");
}
Exercise 2.1.8.2
Compile the following code and execute it under various optimization levels.
Try to guess the dierent types of optimizations that are being performed to improve the
performance as the optimization is increased.
REAL*8 A(1000000)
DO I=1,1000000
A(I) = 3.1415927
15 This content is available online at <http://cnx.org/content/m33699/1.3/>.
16 This content is available online at <http://cnx.org/content/m33700/1.3/>.
Exercise 2.1.8.3
Take the following code segment and compile it at various optimization levels. Look at the gener-
ated assembly language code (S option on some compilers) and nd the eects of each optimization
level on the machine language. Time the program to see the performance at the dierent optimiza-
tion levels. If you have access to multiple architectures, look at the code generated using the same
optimization levels on dierent architectures.
REAL*8 A(1000000)
COMMON/BLK/A
.... Call Time
DO I=1,1000000
A(I) = A(I) + 1.234
ENDDO
.... Call Time
END
It's clear why you might care about the performance of your program if the workload increases. Trying to
cram 25 hours of computing time into a 24-hour day is an administrative nightmare. But why should people
who are considering a new machine care about the runtime? After all, the new machine is presumably faster
than the old one, so everything should take less time. The reason is that when people are evaluating new
machines, they need a basis of comparisona benchmark. People often use familiar programs as benchmarks.
It makes sense: you want a benchmark to be representative of the kind of work you do, and nothing is more
representative of the work you do than the work you do!
Benchmarking sounds easy enough, provided you have timing tools. And you already know the meaning
of time.18 You just want to be sure that what those tools are reporting is the same as what you think you're
17 This content is available online at <http://cnx.org/content/m33704/1.3/>.
18 Time is money.
getting; especially if you have never used the tools before. To illustrate, imagine if someone took your watch
and replaced it with another that expressed time in some funny units or three overlapping sets of hands. It
would be very confusing; you might have a problem reading it at all. You would also be justiably nervous
about conducting your aairs by a watch you don't understand.
UNIX timing tools are like the six-handed watch, reporting three dierent kinds of time measurements.
They aren't giving conicting information they just present more information than you can jam into a
single number. Again, the trick is learning to read the watch. That's what the rst part of this chapter is
about. We'll investigate the dierent types of measurements that determine how a program is doing.
If you plan to tune a program, you need more than timing information. Where is time being spent
in a single loop, subroutine call overhead, or with memory problems? For tuners, the latter sections of this
chapter discuss how to prole code at the procedural and statement levels. We also discuss what proles
mean and how they predict the approach you have to take when, and if, you decide to tweak the code for
performance, and what your chances for success will be.
2.2.2 Timing19
We assume that your program runs correctly. It would be rather ridiculous to time a program that's not
running right, though this doesn't mean it doesn't happen. Depending on what you are doing, you may
be interested in knowing how much time is spent overall, or you may be looking at just a portion of the
program. We show you how to time the whole program rst, and then talk about timing individual loops or
subroutines.
• User time
• System time
• Elapsed time
These timing gures are easier to understand with a little background. As your program runs, it switches
back and forth between two fundamentally dierent modes: user mode
and kernel mode
. The normal
operating state is user mode. It is in user mode that the instructions the compiler generated on your behalf
get executed, in addition to any subroutine library calls linked with your program.20 It might be enough to
run in user mode forever, except that programs generally need other services, such as I/O, and these require
the intervention of the operating system the kernel. A kernel service request made by your program, or
perhaps an event from outside your program, causes a switch from user mode into kernel mode.
Time spent executing in the two modes is accounted for separately. The user time
gure describes time
spent in user mode. Similarly, system timeis a measure of the time spent in kernel mode. As far as user time
goes, each program on the machine is accounted for separately. That is, you won't be charged for activity in
somebody else's application. System time accounting works the same way, for the most part; however, you
can, in some instances, be charged for some system services performed on other people's behalf, in addition
to your own. Incorrect charging occurs because your program may be executing at the moment some outside
activity causes an interrupt. This seems unfair, but take consolation in the fact that it works both ways:
other users may be charged for your system activity too, for the same reason.
19 This content is available online at <http://cnx.org/content/m33706/1.3/>.
20 Cache miss time is buried in here too.
People often record the CPU time and use it as an estimate for elapsed time. Using CPU time is okay on
a single CPU machine, provided you have seen the program run when the machine was quiet and noticed
the two numbers were very close together. But for multiprocessors, the total CPU time can be far dierent
from the elapsed time. Whenever there is a doubt, wait until you have the machine to your- self and time
your program then, using elapsed time. It is very important to produce timing results that can be veried
using another run when the results are being used to make important purchasing decisions.
If you are running on a Berkeley UNIX derivative, the C shell's built-in time command can report a
number of other useful statistics. The default form of the output is shown in Figure 2.3 (The built-in csh
time function). Check with your csh
manual page for more possibilities.
In addition to gures for CPU and elapsed time, csh
time command produces information about CPU
utilization, page faults, swaps, blocked I/O operations (usually disk activity), and some measures of how
much physical memory our pro- gram occupied when it ran. We describe each of them in turn.
Figure 2.3
The second average memory utilization measurement, unshared-memory space, describes the averagereal
storage dedicated to your program's data structures as it ran. This storage includes saved local variables
and COMMON for FORTRAN, and static and external variables for C. We stress the word real here and above
because these numbers talk about physical memory usage, taken over time. It may be that you have allocated
arrays with 1 trillion elements (virtual space), but if your program only crawls into a corner of that space,
your runtime memory requirements will be pretty low.
What the unshared-memory space measurement doesn't tell you, unfortunately, is your program's demand
for memory at its greediest. An application that requires 100 MB 1/10th of the time and 1 KB the rest
of the time appears to need only 10 MB on average not a revealing picture of the program's memory
requirements.
If, for instance, X's primary job is to calculate particle positions, divide by the total time to obtain a number
for particle positions/second. You have to be careful though; too many calls to the timing routines, and
the observer becomes part of the experiment. The timing routines take time too, and their very presence
can increase instruction cache miss or paging. Furthermore, you want X to take a signicant amount of
time so that the measurements are meaningful. Paying attention to the time between timer calls is really
important because the clock used by the timing functions has a limited resolution. An event that occurs
within a fraction of a second is hard to measure with any accuracy.
start = etime(tarray)
finish = etime(tarray)
Not every vendor supplies an etime function; in fact, one doesn't provide a timing routine for FORTRAN at
all. Try it rst. If it shows up as an undened symbol when the program is linked, you can use the following
C routine. It provides the same functionality as etime
:
#include <sys/times.h>
#define TICKS 100.
{
struct tms local;
times (&local);
parts->user= (float) local.tms_utime/TICKS;
parts->system = (float) local.tms_stime/TICKS;
return (parts->user + parts->system);
}
There are a couple of things you might have to tweak to make it work. First of all, linking C routines
with FORTRAN routines on your computer may require you to add an underscore (_) after the function
name. This changes the entry to float etime_ (parts). Furthermore, you might have to adjust the TICKS
parameter. We assumed that the system clock had a resolution of 1/100 of a second (true for the Hewlett-
Packard machines that this version of etime was written for). 1/60 is very common. On an RS-6000 the
number would be 1000. You may nd the value in a le named /usr/include/sys/param.h on your machine,
or you can determine it empirically.
A C routine for retrieving the wall time using calling gettimeofday
is shown below. It is suitable for use
with either C or FORTRAN programs as it uses call-by-value parameter passing:
#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
gettimeofday(&tp, &tzp);
Given that you will often need both CPU and wall time, and you will be continu- ally computing the
dierence between successive calls to these routines, you may want to write a routine to return the elapsed
wall and CPU time upon each call as follows:
SUBROUTINE HPCTIM(WTIME,CTIME)
Figure 2.4
A sharp prole
says that most of the time is spent in one or two procedures, and if you want to improve
the program's performance you should focus your eorts on tuning those procedures. A minor optimization
in a heavily executed line of code can sometimes have a great eect on the overall runtime, given the right
opportunity. A at prole
,23 on the other hand, tells you that the runtime is spread across many routines,
and eort spent optimizing any one or two will have little benet in speeding up the program. Of course,
there are also programs whose execution prole falls somewhere in the middle.
23 The term at prole is a little overloaded. We are using it to describe a prole that shows an even distribution of time
throughout the program. You will also see the label at prole used to draw distinction from a call graph prole, as described
below.
Figure 2.5
We cannot predict with absolute certainty what you are likely to nd when you prole your programs, but
there are some general trends. For instance, engineering and scientic codes built around matrix solutions
often exhibit very sharp proles. The runtime is dominated by the work performed in a handful of routines.
To tune the code, you need to focus your eorts on those routines to make them more ecient. It may
involve restructuring loops to expose parallelism, providing hints to the compiler, or rearranging memory
references. In any case, the challenge is tangible; you can see the problems you have to x.
There are limits to how much tuning one or two routines will improve your runtime, of course. An often
quoted rule of thumb is Amdahl's Law , derived from remarks made in 1967 by one of the designers of the
IBM 360 series, and founder of Amdahl Computer, Gene Amdahl. Strictly speaking, his remarks were about
the performance potential of parallel computers, but people have adapted Amdahl's Law to describe other
things too. For our purposes, it goes like this: Say you have a program with two parts, one that can be
optimized so that it goes innitely fast and another that can't be optimized at all. Even if the optimizable
portion makes up 50% of the initial runtime, at best you will be able to cut the total runtime in half. That
is, your runtime will eventually be dominated by the portion that can't be optimized. This puts an upper
limit on your expectations when tuning.
Even given the nite return on eort suggested by Amdahl's Law, tuning a program with a sharp prole
can be rewarding. Programs with at proles are much more dicult to tune. These are often system codes,
nonnumeric applications, and varieties of numerical codes without matrix solutions. It takes a global tuning
approach to reduce, to any justiable degree, the runtime of a program with a at prole. For instance,
you can sometimes optimize instruction cache usage, which is complicated because of the program's equal
distribution of activity among a large number of routines. It can also help to reduce subroutine call overhead
by folding callees into callers. Occasionally, you can nd a memory reference problem that is endemic to the
whole program and one that can be xed all at once.
When you look at a prole, you might nd an unusually large percentage of time spent in the library
routines such as log, exp, or sin. Often these functions are done in software routines rather than inline.
You may be able to rewrite your code to eliminate some of these operations. Another important pattern to
look for is when a routine takes far longer than you expect. Unexpected execution time may indicate you
are accessing memory in a pattern that is bad for performance or that some aspect of the code cannot be
optimized properly.
In any case, to get a prole, you need a proler. One or two subroutine prolers
come standard with
the software development environments on all UNIX machines. We discuss two of them: andprof . In gprof
addition, we mention a few line-by-line prolers. Subroutine prolers can give you a general overall view of
where time is being spent. You probably should start with prof
, if you have it (most machines do). Otherwise,
use gprof
. After that, you can move to a line-by- line proler if you need to know which statements take the
most time.
2.2.3.1 prof
prof is the most common of the UNIX proling tools. In a sense, it is an extension of the compiler, linker,
and object libraries, plus a few extra utilities, so it is hard to look at any one thing and say this proles
your code. prof works by periodically sampling the program counter as your application runs. To enable
proling, you must recompile and relink using the p ag. For example, if your program has two modules,
stu.c and junk.c, you need to compile and link according to the following code:
% cc stuff.c -p -O -c
% cc junk.c -p -O -c
% cc stuff.o junk.o -p -o stuff
This creates a stu binary that is ready for proling. You don't need to do anything special to run it. Just
treat it normally by entering stuff. Because runtime statistics are being gathered, it takes a little longer
than usual to execute.24 At completion, there is a new le called mon.out
in the directory where you ran it.
This le contains the history of stu
in binary form, so you can't look at it directly. Use the prof
utility to
read mon.outand create a prole of stu
. By default, the information is written to your screen on standard
output, though you can easily redirect it to a le:
main () {
int l;
24 Remember: code with proling enabled takes longer to run. You should recompile and relink the whole thing without the
p ag when you have nished proling.
Again, you need to compile and link loops with the p ag, run the program, and then run the prof utility
to extract a prole, as follows:
% cc loops.c -p -o loops
% ./loops
% prof loops > loops.prof
The following example shows what a loops.prof should look like. There are six columns.
• msec/call Seconds divided by number of calls giving the average length of time taken by each invo-
cation of the routine
• Name The name of this routine
2.2.3.2 gprof
Just as it's important to know how time is distributed when your program runs, it's also valuable to be
able to tell who called who in the list of routines. Imagine, for instance, if something labeled _exp showed
up high in the list in the prof output. You might say: Hmmm, I don't remember calling anything named
exp(). I wonder where that came from. A call tree helps you nd it.
Subroutines and functions can be thought of as members of a family tree. The top of the tree, or root, is
actually a routine that precedes the main routine you coded for the application. It calls your main routine,
which in turn calls others, and so on, all the way down to the leaf nodes of the tree. This tree is properly
known as a call graph
.25 The relationship between routines and nodes in the graph is one of parents and
children. Nodes separated by more than one hop are referred to as ancestors and descendants.
Figure 6-4 graphically depicts the kind of call graph you might see in a small application. main is the
parent or ancestor of most of the rest of the routines. G has two parents, E and C. Another routine, A,
doesn't appear to have any ancestors or descendants at all. This problem can happen when routines are not
compiled with proling enabled, or when they aren't invoked with a subroutine call such as would be the
case if A were an exception handler.
The UNIX proler that can extract this kind of information is called gprof
. It replicates the abilities of
prof , plus it gives a call graph prole so you can see who calls whom, and how often. The call graph prole
is handy if you are trying to gure out how a piece of code works or where an unknown routine came from,
or if you are looking for candidates for subroutine inlining.
To use call graph proling you need go through the same steps as with prof
, except that a pg ag is
substituted for the p ag.26 Additionally, when it comes time to produce the actual prole, you use the
gprof utility instead of prof
. One other dierence is that the name of the statistics le is gmon.out
instead
ofmon.out :
% cc -pg stuff.c -c
% cc stuff.o -pg -o stuff
% stuff
% gprof stuff > stuff.gprof
25 It doesn't have to be a tree. Any subroutine can have more than one parent. Furthermore, recursive subroutine calls
introduce cycles into the graph, in which a child calls one of its parents.
26 On HP machines, the ag is G.
Figure 2.6
The rst section textually maps out the call graph. The second section lists routines, the percentage of
time devoted to each, the number of calls, etc. (similar to prof
). The third section is a cross reference
so that you can locate routines by number, rather than by name. This section is especially useful for large
applications because routines are sorted based on the amount of time they use, and it can be dicult to
locate a particular routine by scanning for its name. Let's invent another trivial application to illustrate
how gprof
works. Figure 2.7 (FORTRAN example) shows a short piece of FORTRAN code, along with a
diagram of how the routines are connected together. Subroutines A and B are both called by MAIN, and, in
turn, each calls C. The following example shows a section of the output from gprof
's call graph prole:27
27 In the interest of conserving space, we clipped out the section most relevant to our discussion and included it in this
example. There was a lot more to it, including calls of setup and system routines, the likes of which you will see when you run
gprof.
FORTRAN example
Figure 2.7
called/total parents
index %time self descendants called+self name index
called/total children
-----------------------------------------------
-----------------------------------------------
Sandwiched between each set of dashed lines is information describing a given routine and its relationship to
parents and children. It is easy to tell which routine the block represents because the name is shifted farther
to the left than the others. Parents are listed above, children below. As with prof, underscores are tacked
onto the labels.28 A description of each of the columns follows:
• index You will notice that each routine name is associated with a number in brackets ([n]). This is
a cross-reference for locating the routine elsewhere in the prole. If, for example, you were looking at
the block describing _MAIN_ and wanted to know more about one of its children, say _a_, you could
nd it by scanning down the left side of the page for its index, [5].
• %time The meaning of the %time eld is a little dierent than it was for prof
. In this case it describes
the percentage of time spent in this routine plus
the time spent in all of its children. It gives you a
quick way to determine where the busiest sections of the call graph can be found.
• self Listed in seconds, the self column has dierent meanings for parents, the routine in question,
and its children. Starting with the middle entry the routine itself the self gure shows how
much overall time was dedicated to the routine. In the case _b_, for instance, this amounts to 3.23
seconds.
Each self column entry shows the amount of time that can be attributed to calls from the parents. If
you look at routine _c_, for example, you will see that it consumed a total time of 3.23 seconds. But
note that it had two parents: 1.62 seconds of the time was attributable to calls from _a_, and 1.62
seconds to _b_.
For the children, the self gure shows how much time was spent executing each child due to calls
from this routine. The children may have consumed more time overall, but the only time accounted
for is time-attributable to calls from this routine. For example, _c_ accumulated 3.23 seconds overall,
but if you look at the block describing _b_, you see _c_ listed as a child with only 1.62 seconds. That's
the total time spent executing _c_ on behalf of _b_.
• descendants As with the self column, gures in the descendants column have dierent meanings for
the routine, its parents, and children. For the routine itself, it shows the number of seconds spent in
all of its descendants.
For the routine's parents, the descendants gure describes how much time spent in the routine can be
traced back to calls by each parent. Looking at routine _c_ again, you can see that of its total time,
3.23 seconds, 1.62 seconds were attributable to each of its two parents, _a_ and _b_.
For the children, the descendants column shows how much of the child's time can be attributed to calls
from this routine. The child may have accumulated more time overall, but the only time displayed is
time associated with calls from this routine.
• calls The calls column shows the number of times each routine was invoked, as well as the distri-
bution of those calls associated with both parents and children. Starting with the routine itself, the
28 You may have noticed that there are two main routines: _MAIN_ and _main. In a FORTRAN program, _MAIN_ is the actual
FORTRAN main routine. It's called as a subroutine by _main, provided from a system library at link time. When you're
proling C code, you won't see _MAIN_.
gure in the calls column shows the total number of entries into the routine. In situations where the
routine called itself, you will also see a +n n
immediately appended, showing that additional calls were
made recursively.
Parent and child gures are expressed as ratios. For the parents, the ratio m/n n
says of the times
the routine was called, m n
of those calls came from this parent. For the child, it says of the times
m
this child was called, of those calls came from this routine.
• %time Again, we see a eld that describes the runtime for each routine as a percent- age of the overall
time taken by the program. As you might expect, all the entries in this column should total 100%
(nearly).
• cumulative seconds For any given routine, the column called cumulative seconds tallies a running
sum of the time taken by all the preceding routines plus its own time. As you scan towards the bottom,
the numbers asymptotically approach the total runtime for the program.
• self seconds Each routine's individual contribution to the runtime.
• calls The number of times this particular routine was called.
• self ms/call Seconds spent inside the routine, divided by the number of calls. This gives the average
length of time taken by each invocation of the routine. The gure is presented in milliseconds.
• total ms/call Seconds spent inside the routine plus its descendants, divided by the number of calls.
• name The name of the routine. Notice that the cross-reference number appears here too.
In the example prole, each run along the way creates a new gmon.outle that we renamed to make room
for the next one. At the end, gprof combines the infor- mation from each of the data les to produce a
summary prole of bar in the legprof.summary.out . Additionally (you don't see it here), gprof
creates a
le named gmon.sum that contains the merged data from the original three data les. gmon.sum
has the
same format as gmon.out , so you can use it as input for other merged proles down the road.
In form, the output from a merged prole looks exactly the same as for an individual run. There are a
couple of interesting things you will note, however. For one thing, the main routine appears to have been
invoked more than once one time for each run, in fact. Furthermore, depending on the application,
multiple runs tend to either smooth the contour of the prole or exaggerate its features. You can imagine
how this might happen. If a single routine is consistently called while others come and go as the input data
changes, it takes on increasing importance in your tuning eorts.
Figure 2.8
BAR and FOO take turns running. In fact, BAR takes more time than FOO. But because the sampling interval
closely matches the frequency at which the two subroutines alternate, we get a quantizing error: most of the
samples happen to be taken while FOO is running. Therefore, the prole tells us that FOO took more CPU
time than BAR.
We have described the tried and true UNIX subroutine prolers that have been available for years. In
many cases, vendors have much better tools available for the asking or for a fee. If you are doing some
serious tuning, ask your vendor representative to look into other tools for you.
2.2.4.1 tcov
tcov, available on Sun workstations and other SPARC machines that run SunOS, gives execution statistics
that describe the number of times each source statement was executed. It is very easy to use. Assume for
illustration that we have a source program called foo.c. The following steps create a basic block prole:
% cc -a foo.c -o foo
% foo
% tcov foo.c
The -a option tells the compiler to include the necessary support for tcov
.31 Several les are created in the
process. One called foo.d
accumulates a history of the exe- cution frequencies within the program . That foo
is, old data is updated with new data each time foo
is run, so you can get an overall picture of what happens
inside foo
, given a variety of data sets. Just remember to clean out the old data if you want to start over.
The prole itself goes into a le called .foo.tcov
Let's look at an illustration. Below is a short C program that performs a bubble sort of 10 integers:
29 This content is available online at <http://cnx.org/content/m33710/1.3/>.
30 A basic block is a section of code with only one entrance and one exit. If you know how many times the block was entered,
you know how many times each of the statements in the block was executed, which gives you a line-by-line prole. The concept
of a basic block is explained in detail in Section 2.1.1
31 On Sun Solaris systems, the xa option is used.
tcov produces a basic block prole that contains execution counts for each source line, plus some summary
statistics (not shown):
The numbers to the left tell you the number of times each block was entered. For instance, you can see that
the routine was entered just once, and that the highest count occurs at the test n[j] < n[j+1]. tcov
shows
more than one count on a line in places where the compiler has created more than one block.
2.2.4.2 pixie
pixie is a little dierent from tcov
. Rather than reporting the number of times each source line was executed,
pixie reports the number of machine clock cycles devoted to executing each line. In theory, you could use
this to calculate the amount of time spent per statement, although anomalies like cache misses are not
represented.
pixieworks by pixifying an executable le that has been compiled and linked in the normal way. Below
we run pixie foo
on to create a new executable called foo.pixie
:
% cc foo.c -o foo
% pixie foo
% foo.pixie
% prof -pixie foo
Below, we have listed the output of the third section for the bubble sort:
Here you can see three entries for the main routine from foo.c
, plus a number of system library routines.
The entries show the associated line number and the number of machine cycles dedicated to executing that
line as the program ran. For instance, line 7 of foo.c
took 605 cycles (12% of the runtime).
2.2.5.1 Gauging the Size of Your Program and the Machine's Memory
How can you tell if you are running out-of-core? There are ways to check for pag- ing on the machine,
but perhaps the most straightforward check is to compare the size of your program against the amount of
available memory. You do this with the size
command:
% size myprogram
The rst three elds describe the amount of memory required for three dierent portions of your program.
The rst, text, accounts for the machine instructions that make up your program. The second, data, includes
initialized values in your pro- gram such as the contents of data statements, common blocks, externals,
character strings, etc. The third component, bss, (block started by symbol), is usually the largest. It
describes an uninitialized data area in your program. This area would be made of common blocks that are
not set by a block data. The last eld is a total for all three sections added together, in bytes.34
Next, you need to know how much memory you have in your system. Unfortunately, there isn't a standard
UNIX command for this. On the RS/6000, /etc/lscfg
tells you. On an SGI machine, /etc/hinv
does it. Many
System V UNIX implementations have an /etc/memsize command. On any Berkeley derivative, you can
type:
% ps aux
This command gives you a listing of all the processes running on the machine. Find the process with the
largest value in the %MEM. Divide the value in the RSS eld by the percentage of memory used to get a rough
gure for how much memory your machine has:
memory = RSS/(%MEM/100)
For instance, if the largest process shows 5% memory usage and a resident set size (RSS) of 840 KB, your
machine has 840000/(5/100) = 16 MB of memory.35 If the answer from the size command shows a total
that is anywhere near the amount of memory you have, you stand a good chance of paging when you run
especially if you are doing other things on the machine at the same time.
% vmstat 5
This command produces output every ve seconds.
% sar -r 5 5
This command shows you the amount of free memory and swap space presently available. If the free memory
gure is low, you can assume that your program is paging:
35 You could also reboot the machine! It will tell you how much memory is available when it comes up.
As we mentioned earlier, if you must run a job larger than the size of the memory on your machine, the
same sort of advice that applied to conserving cache activity applies to paging activity.36 Try to minimize
the stride in your code, and where you can't, blocking memory references helps a whole lot.
A note on memory performance monitoring tools: you should check with your workstation vendor to see
what they have available beyond vmstat sar
or . There may be much more sophisticated (and often graphical)
tools that can help you understand how your program is using memory.
2.2.7 Exercises38
Exercise 2.2.7.1
Prole the following program using gprof. Is there any way to tell how much of the time spent in
routine c was due to recursive calls?
36 By the way, are you getting the message Out of memory? If you are running csh, try typing unlimit to see if the message
goes away. Otherwise, it may mean that you don't have enough swap space available to run the job.
37 This content is available online at <http://cnx.org/content/m33714/1.3/>.
38 This content is available online at <http://cnx.org/content/m33718/1.3/>.
main()
{
int i, n=10;
for (i=0; i<1000; i++) {
c(n);
a(n);
}
}
c(n)
int n;
{
if (n > 0) {
a(n-1);
c(n-1);
}
}
a(n)
int n;
{
c(n);
}
Exercise 2.2.7.2
Prole an engineering code (oating-point intensive) with full optimization on and o. How does
the prole change? Can you explain the change?
Exercise 2.2.7.3
Write a program to determine the overhead of the getrusage and the etime calls. Other than
consuming processor time, how can making a system call to check the time too often alter the
application performance?
DO I=1,N
A(I) = A(I) + B(I) * C
ENDDO
The code below performs the same calculations, but look at what we have done:
DO I=1,N
CALL MADD (A(I), B(I), C)
ENDDO
SUBROUTINE MADD (A,B,C)
A = A + B * C
RETURN
END
Each iteration calls a subroutine to do a small amount of work that was formerly within the loop. This is a
particularly painful example because it involves oating- point calculations. The resulting loss of parallelism,
coupled with the procedure call overhead, might produce code that runs 100 times slower. Remember, these
operations are pipelined, and it takes a certain amount of wind-up time before the throughput reaches one
operation per clock cycle. If there are few oating-point operations to perform between subroutine calls, the
time spent winding up and winding down pipelines gures prominently.
Subroutine and function calls complicate the compiler's ability to eciently man- age COMMON and
external variables, delaying until the last possible moment actually storing them in memory. The compiler
uses registers to hold the live values of many variables. When you make a call, the compiler cannot tell
whether the subroutine will be changing variables that are declared as external or COMMON. Therefore, it's
forced to store any modied external or COMMON variables back into memory so that the callee can nd
them. Likewise, after the call has returned, the same variables have to be reloaded into registers because the
compiler can no longer trust the old, register-resident copies. The penalty for saving and restoring variables
can be substantial, especially if you are using lots of them. It can also be unwarranted if variables that ought
to be local are specied as external or COMMON, as in the following code:
COMMON /USELESS/ K
DO K=1,1000
IF (K .EQ. 1) CALL AUX
ENDDO
In this example, K has been declared as a COMMON variable. It is used only as a do-loop counter, so there
really is no reason for it to be anything but local. However, because it is in a COMMON block, the call to
AUX forces the compiler to store and reload K each iteration. This is because the side eects of the call are
unknown.
So far, it looks as if we are preparing a case for huge main programs without any subroutines or functions!
Not at all. Modularity is important for keeping source code compact and understandable. And frankly, the
need for maintainability and modularity is always more important than the need for small
performance
improvements. However, there are a few approaches for streamlining subroutine calls that don't require you
to scrap modular coding techniques: macros and procedure inlining.
Remember, if the function or subroutine does a reasonable amount of work, procedure call overhead isn't
going to matter very much. However, if one small routine appears as a leaf node in one of the busiest sections
of the call graph, you might want to think about inserting it in appropriate places in the program.
2.3.2.1 Macros
Macros are little procedures that are substituted inline at compile time. Unlike subroutines or functions,
which are included once during the link, macros are replicated every place they are used. When the compiler
makes its rst pass through your program, it looks for patterns that match previous macro denitions and
expands them inline. In fact, in later stages, the compiler sees an expanded macro as source code.
The rst compilation step for a C program is a pass through the C preprocessor, cpp
. This happens
automatically when you invoke the compiler. cpp
expands #define statements inline, replacing the pattern
matched by the macro denition. In the program above, the statement:
a = average(p,q);
a = ((p+q)/2);
You have to be careful how you dene the macro because it literally replaces the pattern located by cpp.
For instance, if the macro denition said:
c = multiply(x+t,y+v);
the resulting expansion would be x+t*y+v probably not what you intended.
If you are a C programmer you may be using macros without being conscious of it. Many C header les
.h
( ) contain macro denitions. In fact, some standard C library functions are really dened as macros in
the header les. For instance, the function getchar
can be linked in when you build your program. If you
have a statement:
#include <stdio.h>
Without a little preparation, the #define statement is rejected by the FORTRAN compiler. The program
rst has to be preprocessed through cpp
to replace the use of AVERAG with its macro denition. It makes
compilation a two-step procedure, but that shouldn't be too much of a burden, especially if you are building
your programs under the control of the make
utility. We would also suggest you store FORTRAN programs
containing cppdirectives under lename.F
to distinguish them from unadorned FORTRAN. Just be sure
you make your changes only to the .F
les and not to the output from cpp
. This is how you would preprocess
FORTRAN .F
les by hand:
The FORTRAN compiler never sees the original code. Instead, the macro denition is substituted inline as
if you had typed it yourself:
C
PROGRAM MAIN
REAL A,P,Q
42 Some programmers use the standard UNIX m4 preprocessor for FORTRAN
The directives and compile line options are not standard, so you have to check your compiler documentation.
Unfortunately, you may learn that there is no such feature (yet, always yet), or that it's an expensive extra.
The third form of inlining in the list, automatic, is available from just a few vendors. Automatic inlining
depends on a sophisticated compiler that can view the denitions of several modules at once.
There are some words of caution with regard to procedure inlining. You can easily do too much of it. If
everything and anything is ingested into the body of its parents, the resulting executable may be so large that
it repeatedly spills out of the instruction cache and becomes a net performance loss. Our advice is that you
use the caller/callee information prolers give you and make some intelligent decisions about inlining, rather
than trying to inline every subroutine available. Again, small routines that are called often are generally the
best candidates for inlining.
2.3.3 Branches43
People sometimes take a week to make a decision, so we can't fault a computer if it takes a few tens of
nanoseconds. However, if an if-statement appears in some heavily traveled section of the code, you might
get tired of the delay. There are two basic approaches to reducing the impact of branches:
• Streamline them.
• Move them out to the computational suburbs. Particularly, get them out of loops.
In Section 2.3.4 we show you some easy ways to reorganize conditionals so they execute more quickly.
43 This content is available online at <http://cnx.org/content/m33722/1.3/>.
The idea was that if the multiplier, A(I), were reasonably small, there would be no reason to perform the
math in the center of the loop. Because oating-point operations weren't pipelined on many machines, a
comparison and a branch was cheaper; the test would save time. On an older CISC or early RISC processor,
a comparison and branch is probably still a savings. But on other architectures, it costs a lot less to just
perform the math and skip the test. Eliminating the branch eliminates a control dependency and allows
the compiler to pipeline more arithmetic operations. Of course, the answer could change slightly if the test
is eliminated. It then becomes a question of whether the dierence is signicant. Here's another example
where a branch isn't necessary. The loop nds the absolute value of each element in an array:
DO I=1,N
IF (A(I) .LT. 0.) A(I) = -A(I)
ENDDO
But why perform the test at all? On most machines, it's quicker to perform the abs() operation on every
element of the array.
We do have to give you a warning, though: if you are coding in C, the absolute value, fabs(), is a
subroutine call. In this particular case, you are better o leaving the conditional in the loop.45
When you can't always throw out the conditional, there are things you can do to minimize negative
performance. First, we have to learn to recognize which conditionals within loops can be restructured and
which cannot. Conditionals in loops fall into several categories:
DO I=1,K
IF (N .EQ. 0) THEN
A(I) = A(I) + B(I) * C
ELSE
A(I) = 0.
ENDIF
ENDDO
Invariant means that the outcome is always the same. Regardless of what happens to the variables A, B,
C, and I, the value of N won't change, so neither will the outcome of the test.
You can recast the loop by making the test outside and replicating the loop body twice once for when
the test is true, and once for when it is false, as in the following example:
IF (N .EQ. 0) THEN
DO I=1,K
A(I) = A(I) + B(I) * C
ENDDO
ELSE
DO I=1,K
A(I) = 0
ENDDO
ENDIF
The eect on the runtime is dramatic. Not only have we eliminated K-1 copies of the test, we have also
assured that the computations in the middle of the loop are not control-dependent on the if-statement, and
are therefore much easier for the compiler to pipeline.
We remember helping someone optimize a program with loops containing similar conditionals. They were
checking to see whether debug output should be printed each iteration inside an otherwise highly optimizable
loop. We can't fault the person for not realizing how much this slowed the program down. Performance
wasn't important at the time. The programmer was just trying to get the code to produce good answers.
But later on, when performance mattered, by cleaning up invariant conditionals, we were able to speed up
the program by a factor of 100.
DO I=1,N
DO J=1,N
IF (J .LT. I)
A(J,I) = A(J,I) + B(J,I) * C
ELSE
A(J,I) = 0.0
ENDIF
ENDDO
ENDDO
Notice how the if-statement partitions the iterations into distinct sets: those for which it is true and those
for which it is false. You can take advantage of the predictability of the test to restructure the loop into
several loops each custom-made for a dierent partition:
DO I=1,N
DO J=1,I-1
A(J,I) = A(J,I) + B(J,I) * C
ENDDO
DO J=I,N
A(J,I) = 0.0
ENDDO
ENDDO
The new version will almost always be faster. A possible exception is when N is a small value, like 3, in
which case we have created more clutter. But then, the loop probably has such a small impact on the total
runtime that it won't matter which way it's coded.
DO I=1,N
DO J=1,N
IF (B(J,I) .GT. 1.0) A(J,I) = A(J,I) + B(J,I) * C
ENDDO
ENDDO
There is not much you can do about this type of conditional. But because every iteration is independent,
the loop can be unrolled or can be performed in parallel.
DO I=1,N
IF (X .LT. A(I)) X = X + B(I)*2.
ENDDO
You can't know which way the branch will go for the next iteration until you are done with the current
iteration. To recognize the dependency, try to unroll the loop slightly by hand. If you can't start the second
test until the rst has nished, you have a dependent loop conditional. You may want to look at these types
of loops to see if you can eliminate the iteration-to-iteration value.
2.3.4.5 Reductions
Keep an eye out for loops in which the if-statement is performing a max or min function on a array. This is
areduction , so called because it reduces a array to a scalar result (the previous example was a reduction too,
by the way). Again, we are getting a little bit ahead of ourselves, but since we are talking about if-statements
in loops, I want to introduce a trick for restructuring reductions max and min to expose more parallelism.
The following loop searches for the maximum value, z, in the array a by going through the elements one at
a time:
As written, it's recursive like the loop from the previous section. You need the result of a given iteration
before you can proceed to the next. However, since we are looking for the greatest element in the whole
array, and since that will be the same element (essentially) no matter how we go about looking for it, we
can restructure the loop to check several elements at a time (we assume n is evenly divisible by 2 and do not
include the preconditioning loop):
z0 = 0.;
z1 = 0.;
for (i=0; i< n-1; i+=2) {
z0 = z0 < a[i] ? a[i] : z0;
z1 = z1 < a[i+1] ? a[i+1] : z1;
}
z = z0 < z1 ? z1 : z0;
Do you see how the new loop calculates two new maximum values each iteration? These maximums are
then compared with one another, and the winner becomes the new ocial max
. It's analogous to a play-o
arrangement in a Ping-Pong tournament. Whereas the old loop was more like two players competing at a
time while the rest sat around, the new loop runs several matches side by side. In general this particular
optimization is not a good one to code by hand. On parallel processors, the compiler performs the reduction
in its own way. If you hand-code similar to this example, you may inadvertently limit the compiler's exibility
on a parallel system.
DO I=1,N
DO J=1,N
IF (B(J,I) .EQ. 0 ) THEN
PRINT *,I,J
STOP
ENDIF
A(J,I) = A(J,I) / B(J,I)
ENDDO
ENDDO
Avoiding these tests is one of the reasons that the designers of the IEEE oating- point standard added the
trap feature for operations such as dividing by zero. These traps allow the programmer in a performance-
critical section of the code to achieve maximum performance yet still detect when an error occurs.
INTEGER NUMEL, I
PARAMETER (NUMEL = 1000)
REAL*8 A(NUMEL)
REAL*4 B(NUMEL)
DO I=1,NUMEL
A(I) = A(I) + B(I)
ENDDO
In each iteration, B(I) has to be promoted to double precision before the addition can occur. You don't see
the promotion in the source code, but it's there, and it takes time.
C programmers beware: in Kernighan and Ritchie (K&R) C, all oating-point calculations in C programs
take place in double precision even if all the variables involved are declared as oat
. It is possible for you
to write a whole K+R application in one precision, yet suer the penalty of many type conversions.
Another data typerelated mistake is to use character operations in IF tests. On many systems, character
operations have poorer performance than integer operations since they may be done via procedure calls. Also,
the optimizers may not look at code using character variables as a good candidate for optimization. For
example, the following code:
DO I=1,10000
IF ( CHVAR(I) .EQ. 'Y' ) THEN
A(I) = A(I) + B(I)*C
ENDIF
ENDDO
might be better written using an integer variable to indicate whether or not a computation should be
performed:
DO I=1,10000
IF ( IFLAG(I) .EQ. 1 ) THEN
A(I) = A(I) + B(I)*C
ENDIF
ENDDO
Another way to write the code, assuming the IFLAG variable was 0 or 1, would be as follows:
DO I=1,10000
The last approach might actually perform slower on some computer systems than the approach using the IF
and the integer variable.
c = a + b + d
e = q + a + b
becomes:
temp = a + b
c = temp + d
e = q + temp
Substituting for a+b eliminates some of the arithmetic. If the expression is reused many times, the savings can
be signicant. However, a compiler's ability to recognize common subexpressions is limited, especially when
there are multiple components, or their order is permuted. A compiler might not recognize that a+b+c and
c+b+a are equivalent.48 For important parts of the program, you might consider doing common subexpression
elimination of complicated expressions by hand. This guarantees that it gets done. It compromises beauty
somewhat, but there are some situations where it is worth it.
Here's another example in which the function sin is called twice with the same argument:
x = r*sin(a)*cos(b);
y = r*sin(a)*sin(b);
z = r*cos(a);
becomes:
temp = r*sin(a);
48 And because of overow and round-o errors in oating-point, in some situations they might not be equivalent.
We have replaced one of the calls with a temporary variable. We agree, the savings for eliminating one
transcendental function call out of ve won't win you a Nobel prize, but it does call attention to an important
point: compilers typically do not perform common subexpression elimination over subroutine or function
calls. The compiler can't be sure that the subroutine call doesn't change the state of the argument or some
other variables that it can't see.
The only time a compiler might eliminate common subexpressions containing function calls is when they
are intrinsics, as in FORTRAN. This can be done because the compiler can assume some things about their
side eects. You, on the other hand, can see into subroutines, which means you are better qualied than the
compiler to group together common subexpressions involving subroutines or functions.
DO I=1,N
A(I) = A(I) / SQRT(X*X + Y*Y)
ENDDO
becomes:
We hoisted an expensive, invariant operation out of the loop and assigned the result to a temporary variable.
Notice, too, that we made an algebraic simplication when we exchanged a division for multiplication by an
inverse. The multiplication will execute much more quickly. Your compiler might be smart enough to make
these transformations itself, assuming you have instructed the compiler that these are legal transformations;
but without crawling through the assembly language, you can't be positive. Of course, if you rearrange
code by hand and the runtime for the loop suddenly goes down, you will know that the compiler has been
sandbagging all along.
Sometimes you want to sink an operation below the loop. Usually, it's some calculation performed each
iteration but whose result is only needed for the last. To illustrate, here's a sort of loop that is dierent from
the ones we have been looking at. It searches for the nal character in a character string:
becomes:
The new version of the loop moves the assignment of c beyond the last iteration. Admittedly, this transfor-
mation would be a reach for a compiler and the savings wouldn't even be that great. But it illustrates the
notion of sinking an operation very well.
Again, hoisting or sinking instructions to get them out of loops is something your compiler should be
capable of doing. But often you can slightly restructure the calculations yourself when you move them to
get an even greater benet.
DO I=1,N
XOLD(I) = X(I)
X(I)= X(I) + XINC(I)
ENDDO
In reality, the steps that go into retrieving X(I) are just additional common subex- pressions: an address
calculation (possibly) and a memory load operation. You can see that the operation is repeated by rewriting
the loop slightly:
DO I=1,N
TEMP= X(I)
XOLD(I) = TEMP
X(I)= TEMP + XINC(I)
ENDDO
FORTRAN compilers should recognize that the same X(I) is being used twice and that it only needs to
be loaded once, but compilers aren't always so smart. You sometimes have to create a temporary scalar
Unless the compiler can see the denitions of x, xinc, and xold, it has to assume that they are pointers
leading back to the same storage, and repeat the loads and stores. In this case, introducing temporary
variables to hold the values x, xinc, and xold is an optimization the compiler wasn't free to make.
Interestingly, while putting scalar temporaries in the loop is useful for RISC and superscalar machines,
it doesn't help code that runs on parallel hardware. A parallel compiler looks for opportunities to eliminate
the scalars or, at the very least, to replace them with temporary vectors. If you run your code on a parallel
machine from time to time, you might want to be careful about introducing scalar temporary variables into
a loop. A dubious performance gain in one instance could be a real performance loss in another.
Sometimes this means we make changes that are not beautiful. However, they are often quick.
2.3.7 Exercises50
Exercise 2.3.7.1
How would you simplify the following loop conditional?
DO I=1,N
A(I) = A(I) * B
IF (I .EQ. N/2) A(I) = 0.
ENDDO
Exercise 2.3.7.2
Time this loop on your computer, both with and without the test. Run it with three sets of data:
one with all A(I)s less than SMALL, one with all A(I)s greater than SMALL, and one with an even
split. When is it better to leave the test in the loop, if ever?
Exercise 2.3.7.3
Write a simple program that calls a simple subroutine in its inner loop. Time the program
execution. Then tell the compiler to inline the routine and test the performance again. Finally,
modify the code to perform the operations in the body of the loop and time the code. Which option
ran faster? You may have to look at the generated machine code to gure out why.
• Loop unrolling
• Nested loop optimization
• Loop interchange
• Memory reference optimization
• Blocking
• Out-of-core solutions
Someday, it may be possible for a compiler to perform all these loop optimizations automatically. Typically
loop unrolling is performed as part of the normal compiler optimizations. Other optimizations may have to
51 This content is available online at <http://cnx.org/content/m33728/1.3/>.
DO I=1,N
A(I,J,K) = A(I,J,K) + B(J,I,K)
ENDDO
This loop contains one oating-point addition and three memory references (two loads and a store). There
are some complicated array index expressions, but these will probably be simplied by the compiler and
executed in the same cycle as the memory and oating-point operations. For each iteration of the loop, we
must increment the index variable and test to determine if the loop has completed.
A 3:1 ratio of memory references to oating-point operations suggests that we can hope for no more than
1/3 peak oating-point performance from the loop unless we have more than one path to memory. That's
bad news, but good information. The ratio tells us that we ought to consider memory reference optimizations
rst.
The loop below contains one oating-point addition and two memory operations a load and a store.
Operand B(J) is loop-invariant, so its value only needs to be loaded once, upon entry to the loop:
DO I=1,N
A(I) = A(I) + B(J)
ENDDO
Again, our oating-point throughput is limited, though not as severely as in the previous loop. The ratio of
memory references to oating-point operations is 2:1.
52 This content is available online at <http://cnx.org/content/m33729/1.3/>.
53 Take a look at the assembly language output to be sure, which may be going a bit overboard. To get an assembly language
listing on most machines, compile with the S ag. On an RS/6000, use the qlist ag.
54 The compiler reduces the complexity of loop index expressions with a technique called induction variable simplication.
See Section 2.1.1.
The next example shows a loop with better prospects. It performs element-wise multiplication of two
vectors of complex numbers and assigns the results back to the rst. There are six memory operations (four
loads and two stores) and six oating-point operations (two additions and four multiplications):
It appears that this loop is roughly balanced for a processor that can perform the same number of memory
operations and oating-point operations per cycle. However, it might not be. Many processors perform a
oating-point multiply and add in a single instruction. If the compiler is good enough to recognize that the
multiply-add is appropriate, this loop may also be limited by memory references; each iteration would be
compiled into two multiplications and two multiply-adds.
Again, operation counting is a simple way to estimate how well the requirements of a loop will map onto
the capabilities of the machine. For many loops, you often nd the performance of the loops dominated by
memory references, as we have seen in the last three examples. This suggests that memory reference tuning
is very important.
DO I=1,N
A(I) = A(I) + B(I) * C
ENDDO
You can unroll the loop, as we have below, giving you the same operations in fewer iterations with less loop
overhead. You can imagine how this would help on any computer. Because the computations in one iteration
do not depend on the computations in other iterations, calculations from dierent iterations can be executed
together. On a superscalar processor, portions of these four statements may actually execute in parallel:
55 This content is available online at <http://cnx.org/content/m33732/1.3/>.
DO I=1,N,4
A(I) = A(I) + B(I) * C
A(I+1) = A(I+1) + B(I+1) * C
A(I+2) = A(I+2) + B(I+2) * C
A(I+3) = A(I+3) + B(I+3) * C
ENDDO
II = IMOD (N,4)
DO I=1,II
A(I) = A(I) + B(I) * C
ENDDO
DO I=1+II,N,4
A(I) = A(I) + B(I) * C
A(I+1) = A(I+1) + B(I+1) * C
A(I+2) = A(I+2) + B(I+2) * C
A(I+3) = A(I+3) + B(I+3) * C
ENDDO
The number of iterations needed in the preconditioning loop is the total iteration count modulo for this
unrolling amount. If, at runtime, N turns out to be divisible by 4, there are no spare iterations, and the
preconditioning loop isn't executed.
Speculative execution in the post-RISC architecture can reduce or eliminate the need for unrolling a loop
that will operate on values that must be retrieved from main memory. Because the load operations take
such a long time relative to the computations, the loop is naturally unrolled. While the processor is waiting
for the rst load to nish, it may speculatively execute three to four iterations of the loop ahead of the rst
load, eectively unrolling the loop in the Instruction Reorder Buer.
PARAMETER (NITER = 3)
DO I=1,NITER
A(I) = B(I) * C
ENDDO
Because NITER is hardwired to 3, you can safely unroll to a depth of 3 without worrying about a precon-
ditioning loop. In fact, you can throw out the loop structure altogether and leave just the unrolled loop
innards:
PARAMETER (NITER = 3)
A(1) = B(1) * C
A(2) = B(2) * C
A(3) = A(3) * C
Of course, if a loop's trip count is low, it probably won't contribute signicantly to the overall runtime,
unless you nd such a loop at the center of a larger loop. Then you either want to unroll it completely or
leave it alone.
DO I=1,N
DO J=1,N
IF (B(J,I) .GT. 1.0) A(J,I) = A(J,I) + B(J,I) * C
ENDDO
ENDDO
Each iteration is independent of every other, so unrolling it won't be a problem. We'll just leave the outer
loop undisturbed:
II = IMOD (N,4)
DO I=1,N
DO J=1,II
IF (B(J,I) .GT. 1.0)
+ A(J,I) = A(J,I) + B(J,I) * C
ENDDO
DO J=II+1,N,4
IF (B(J,I) .GT. 1.0)
+ A(J,I) = A(J,I) + B(J,I) * C
IF (B(J+1,I) .GT. 1.0)
+ A(J+1,I) = A(J+1,I) + B(J+1,I) * C
IF (B(J+2,I) .GT. 1.0)
+ A(J+2,I) = A(J+2,I) + B(J+2,I) * C
IF (B(J+3,I) .GT. 1.0)
+ A(J+3,I) = A(J+3,I) + B(J+3,I) * C
ENDDO
ENDDO
This approach works particularly well if the processor you are using supports conditional execution. As
described earlier, conditional execution can replace a branch and an operation with a single conditionally
executed assignment. On a superscalar processor with conditional execution, this unrolled loop executes
quite nicely.
To unroll an outer loop, you pick one of the outer loop index variables and replicate the innermost loop
body so that several iterations are performed at the same time, just like we saw in the Section 2.4.4. The
dierence is in the index variable for which you unroll. In the code below, we have unrolled the middle (j)
loop twice:
We left the k loop untouched; however, we could unroll that one, too. That would give us outer and inner
loop unrolling at the same time:
57 This content is available online at <http://cnx.org/content/m33734/1.3/>.
We could even unroll the i loop too, leaving eight copies of the loop innards. (Notice that we completely
ignored preconditioning; in a real application, of course, we couldn't.)
DO I=1,N
DO J=1,M
A(J,I) = B(J,I) + C(J,I) * D
ENDDO
ENDDO
Unrolling the I loop gives you lots of oating-point operations that can be overlapped:
II = IMOD (N,4)
DO I=1,II
DO J=1,M
A(J,I) = B(J,I) + C(J,I) * D
ENDDO
ENDDO
DO I=II,N,4
DO J=1,M
A(J,I) = B(J,I) + C(J,I) * D
A(J,I+1) = B(J,I+1) + C(J,I+1) * D
A(J,I+2) = B(J,I+2) + C(J,I+2) * D
A(J,I+3) = B(J,I+3) + C(J,I+3) * D
ENDDO
ENDDO
In this particular case, there is bad news to go with the good news: unrolling the outer loop causes strided
memory references on A, B, and C. However, it probably won't be too much of a problem because the inner
loop trip count is small, so it naturally groups references to conserve cache entries.
Outer loop unrolling can also be helpful when you have a nest with recursion in the inner loop, but not
in the outer loops. In this next example, there is a rst- order linear recursion in the inner loop:
DO J=1,M
DO I=2,N
A(I,J) = A(I,J) + A(I-1,J) * B
ENDDO
ENDDO
Because of the recursion, we can't unroll the inner loop, but we can work on several copies of the outer loop
at the same time. When unrolled, it looks like this:
JJ = IMOD (M,4)
DO J=1,JJ
DO I=2,N
A(I,J) = A(I,J) + A(I-1,J) * B
ENDDO
ENDDO
DO J=1+JJ,M,4
DO I=2,N
A(I,J) = A(I,J) + A(I-1,J) * B
A(I,J+1) = A(I,J+1) + A(I-1,J+1) * B
A(I,J+2) = A(I,J+2) + A(I-1,J+2) * B
A(I,J+3) = A(I,J+3) + A(I-1,J+3) * B
ENDDO
ENDDO
You can see the recursion still exists in the I loop, but we have succeeded in nding lots of work to do
anyway.
Sometimes the reason for unrolling the outer loop is to get a hold of much larger chunks of things that
can be done in parallel. If the outer loop iterations are independent, and the inner loop trip count is high,
then each outer loop iteration represents a signicant, parallel chunk of work. On a single CPU that doesn't
matter much, but on a tightly coupled multiprocessor, it can translate into a tremendous increase in speeds.
In practice, KDIM is probably equal to 2 or 3, where J or I, representing the number of points, may be in the
thousands. The way it is written, the inner loop has a very low trip count, making it a poor candidate for
unrolling.
By interchanging the loops, you update one quantity at a time, across all of the points. For tuning
purposes, this moves larger trip counts into the inner loop and allows you to do some strategic unrolling:
DO K=1,KDIM
DO J=1,JDIM
DO I=1,IDIM
D(K,J,I) = D(K,J,I) + V(K,J,I) * DT
ENDDO
ENDDO
ENDDO
This example is straightforward; it's easy to see that there are no inter-iteration dependencies. But how
can you tell, in general, when two loops can be inter- changed? Interchanging loops might violate some
dependency, or worse, only violate it occasionally, meaning you might not catch it when optimizing. Can we
interchange the loops below?
DO I=1,N-1
DO J=2,N
A(I,J) = A(I+1,J-1) * B(I,J)
C(I,J) = B(J,I)
ENDDO
ENDDO
While it is possible to examine the loops by hand and determine the dependencies, it is much better if
the compiler can make the determination. Very few single-processor compilers automatically perform loop
interchange. However, the compilers for high-end vector and parallel computers generally interchange loops
if there is some benet and if interchanging the loops won't alter the program results.60
DO J=1,N
DO I=1,N
A(I,J) = B(I,J) + C(I,J) * D
ENDDO
ENDDO
In contrast, the next loop is slower because its stride is N (which, we assume, is greater than 1). As N
increases from one to the length of the cache line (adjusting for the length of each element), the performance
worsens. Once N is longer than the length of the cache line (again adjusted for element size), the performance
won't decrease:
DO J=1,N
DO I=1,N
A(J,I) = B(J,I) + C(J,I) * D
ENDDO
ENDDO
60 When the compiler performs automatic parallel optimization, it prefers to run the outermost loop in parallel to minimize
overhead and unroll the innermost loop to make best use of a superscalar or vector processor. For this reason, the compiler
needs to have some exibility in ordering the loops in a loop nest.
61 This content is available online at <http://cnx.org/content/m33738/1.3/>.
DO J=1,N
DO I=1,N
A(J,I) = B(J,I) + C(J,I) * D
ENDDO
ENDDO
After interchange, A, B, and C are referenced with the leftmost subscript varying most quickly. This modi-
cation can make an important dierence in performance. We traded three N-strided memory references for
unit strides:
DO I=1,N
DO J=1,N
A(J,I) = B(J,I) + C(J,I) * D
ENDDO
ENDDO
DO I=1,N
DO J=1,N
SUM = 0
DO K=1,N
SUM = SUM + A(I,K) * B(K,J)
ENDDO
C(I,J) = SUM
ENDDO
ENDDO
The problem with this loop is that the A(I,K) will be non-unit stride. Each iteration in the inner loop
consists of two loads (one non-unit stride), a multiplication, and an addition.
Given the nature of the matrix multiplication, it might appear that you can't eliminate the non-unit
stride. However, with a simple rewrite of the loops all the memory accesses can be made unit stride:
DO J=1,N
DO I=1,N
C(I,J) = 0.0
ENDDO
ENDDO
DO K=1,N
DO J=1,N
SCALE = B(K,J)
DO I=1,N
C(I,J) = C(I,J) + A(I,K) * SCALE
ENDDO
ENDDO
ENDDO
Now, the inner loop accesses memory using unit stride. Each iteration performs two loads, one store, a
multiplication, and an addition. When comparing this to the previous loop, the non-unit stride loads have
been eliminated, but there is an additional store operation. Assuming that we are operating on a cache-based
system, and the matrix is larger than the cache, this extra store won't add much to the execution time. The
store is to the location in C(I,J) that was used in the load. In most cases, the store is to a line that is
already in the in the cache. The B(K,J) becomes a constant scaling factor within the inner loop.
DO I=1,N DO 20 J=1,M
DO J=1,M DO 10 I=1,N
A(J,I) = B(I,J) A(J,I) = B(I,J)
ENDDO ENDDO
ENDDO ENDDO
Whichever way you interchange them, you will break the memory access pattern for either A or B. Even
more interesting, you have to make a choice between strided loads vs. strided stores: which will it be?65
We really need a general method for improving the memory access patterns for both
A and B, not one or the
other. We'll show you such a method in Section 2.4.9.
DO I=1,N
DO J=1,N
A(J,I) = A(J,I) + B(I,J)
ENDDO
ENDDO
This loop involves two vectors. One is referenced with unit stride, the other with a stride of N. We can
interchange the loops, but one way or another we still have N-strided array references on either A or B, either
of which is undesirable. The trick is to block
references so that you grab a few elements of A, and then a
few of B, and then a few of A, and so on in neighborhoods. We make this happen by combining inner and
outer loop unrolling:
DO I=1,N,2
DO J=1,N,2
A(J,I) = A(J,I) + B(I,J)
A(J+1,I) = A(J+1,I) + B(I,J+1)
A(J,I+1) = A(J,I+1) + B(I+1,J)
A(J+1,I+1) = A(J+1,I+1) + B(I+1,J+1)
ENDDO
ENDDO
Use your imagination so we can show why this helps. Usually, when we think of a two-dimensional array, we
think of a rectangle or a square (see Figure 2.9 (Arrays A and B)). Remember, to make programming easier,
the compiler provides the illusion that two-dimensional arrays A and B are rectangular plots of memory as
in Figure 2.9 (Arrays A and B). Actually, memory is sequential storage. In FORTRAN, a two-dimensional
array is constructed in memory by logically lining memory strips up against each other, like the pickets
of a cedar fence. (It's the other way around in C: rows are stacked on top of one another.) Array storage
starts at the upper left, proceeds down to the bottom, and then starts over at the top of the next column.
Stepping through the array with unit stride traces out the shape of a backwards N, repeated over and over,
moving to the right.
Arrays A and B
Figure 2.9
Figure 2.10
Imagine that the thin horizontal lines of Figure 2.10 (How array elements are stored) cut memory storage
into pieces the size of individual cache entries. Picture how the loop will traverse them. Because of their
index expressions, references to A go from top to bottom (in the backwards N shape), consuming every bit
of each cache line, but references to B dash o to the right, using one piece of each cache entry and discarding
the rest (see Figure 2.11 (2×2 squares), top). This low usage of cache entries will result in a high number of
cache misses.
If we could somehow rearrange the loop so that it consumed the arrays in small rectangles, rather than
strips, we could conserve some of the cache entries that are being discarded. This is exactly what we
accomplished by unrolling both the inner and outer loops, as in the following example. Array A is referenced
in several strips side by side, from top to bottom, while B is referenced in several strips side by side, from
left to right (see Figure 2.11 (2×2 squares), bottom). This improves cache performance and lowers runtime.
For really big problems, more than cache entries are at stake. On virtual memory machines, memory
references have to be translated through a TLB. If you are dealing with large arrays, TLB misses, in addition
to cache misses, are going to add to your runtime.
2×2 squares
Figure 2.11
Here's something that may surprise you. In the code below, we rewrite this loop yet again, this time
blocking references at two dierent levels: in 2×2 squares to save cache entries, and by cutting the original
loop in two parts to save TLB entries:
DO I=1,N,2
DO J=1,N/2,2
A(J,I) = A(J,I) + B(I,J)
A(J+1,I) = A(J+1,I) + B(I+1,J)
A(J,I+1) = A(J,I+1) + B(I+1,J)
A(J+1,I+1) = A(J+1,I+1) + B(I+1,J+1)
You might guess that adding more loops would be the wrong thing to do. But if you work with a reasonably
large value of N, say 512, you will see a signicant increase in performance. This is because the two arrays
A and B are each 256 KB × 8 bytes = 2 MB when N is equal to 512 larger than can be handled by the
TLBs and caches of most processors.
The two boxes in Figure 2.12 (Picture of unblocked versus blocked references) illustrate how the rst few
references to A and B look superimposed upon one another in the blocked and unblocked cases. Unblocked
references to B zing o through memory, eating through cache and TLB entries. Blocked references are more
sparing with the memory system.
Figure 2.12
You can take blocking even further for larger problems. This code shows another method that limits the
size of the inner loop and visits it repeatedly:
II = MOD (N,16)
JJ = MOD (N,4)
DO I=1,N
DO J=1,JJ
A(J,I) = A(J,I) + B(J,I)
ENDDO
ENDDO
DO I=1,II
DO J=JJ+1,N
A(J,I) = A(J,I) + B(J,I)
A(J,I) = A(J,I) + 1.0D0
ENDDO
ENDDO
DO I=II+1,N,16
DO J=JJ+1,N,4
DO K=I,I+15
A(J,K) = A(J,K) + B(K,J)
A(J+1,K) = A(J+1,K) + B(K,J+1)
A(J+2,K) = A(J+2,K) + B(K,J+2)
A(J+3,K) = A(J+3,K) + B(K,J+3)
ENDDO
ENDDO
ENDDO
Where the inner I loop used to execute N iterations at a time, the new K loop executes only 16 iterations.
This divides and conquers a large memory address space by cutting it into little pieces.
While these blocking techniques begin to have diminishing returns on single-processor systems, on large
multiprocessor systems with nonuniform memory access (NUMA), there can be signicant benet in carefully
arranging memory accesses to maximize reuse of both cache lines and main memory pages.
Again, the combined unrolling and blocking techniques we just showed you are for loops with mixed
stride expressions. They work very well for loop nests like the one we have been looking at. However, if all
array references are strided the same way, you will want to try loop unrolling or loop interchange rst.
With a software-managed approach, the programmer has recognized that the problem is too big and has
modied the source code to move sections of the data out to disk for retrieval at a later time. The other
method depends on the computer's memory system handling the secondary storage requirements on its own,
some- times at a great cost in runtime.
67 This content is available online at <http://cnx.org/content/m33770/1.3/>.
2.4.12 Exercises69
Exercise 2.4.12.1
Why is an unrolling amount of three or four iterations generally sucient for simple vector loops
on a RISC processor? What relationship does the unrolling amount have to oating-point pipeline
depths?
Exercise 2.4.12.2
On a processor that can execute one oating-point multiply, one oating-point addi-
tion/subtraction, and one memory reference per cycle, what's the best performance you could
expect from the following loop?
DO I = 1,10000
A(I) = B(I) * C(I) - D(I) * E(I)
ENDDO
Exercise 2.4.12.3
Try unrolling, interchanging, or blocking the loop in subroutine BAZFAZ to increase the performance.
What method or combination of methods works best? Look at the assembly language created by
the compiler to see what its approach is at the highest level of optimization.
note: Compile the main routine and BAZFAZ separately; adjust NTIMES so that the untuned run
takes about one minute; and use the compiler's default optimization level.
PROGRAM MAIN
IMPLICIT NONE
INTEGER M,N,I,J
PARAMETER (N = 512, M = 640, NTIMES = 500)
DOUBLE PRECISION Q(N,M), R(M,N)
C
DO I=1,M
DO J=1,N
Q(J,I) = 1.0D0
R(I,J) = 1.0D0
ENDDO
ENDDO
C
DO I=1,NTIMES
CALL BAZFAZ (Q,R,N,M)
ENDDO
END
Exercise 2.4.12.4
Code the matrix multiplication algorithm in the straightforward manner and compile it with
various optimization levels. See if the compiler performs any type of loop interchange.
Try the same experiment with the following code:
DO I=1,N
Do you see a dierence in the compiler's ability to optimize these two loops? If you see a dierence,
explain it.
Exercise 2.4.12.5
Code the matrix multiplication algorithm both the ways shown in this chapter. Execute the
program for a range of values for N. Graph the execution time divided by N3 for values of N
ranging from 50×50 to 500×500. Explain the performance you see.
• When performing a 32-bit integer addition, using a carry lookahead adder, you can partially add bits
0 and 1 at the same time as bits 2 and 3.
• On a pipelined processor, while decoding one instruction, you can fetch the next instruction.
• On a two-way superscalar processor, you can execute any combination of an integer and a oating-point
instruction in a single cycle.
• On a multiprocessor, you can divide the iterations of a loop among the four processors of the system.
• You can split a large array across four workstations attached to a network. Each workstation can
operate on its local information and then exchange boundary values at the end of each time step.
DO I=1,16000
A(I) = B(I) * 3.14159
ENDDO
123
124 CHAPTER 3. SHARED-MEMORY PARALLEL PROCESSORS
We have expressed the loop in a way that would imply that A(1) must be computed rst, followed by
A(2), and so on. However, once the loop was completed, it would not have mattered if A(16000), were
computed rst followed by A(15999), and so on. The loop could have computed the even values of I and
then computed the odd values of I. It would not even make a dierence if all 16,000 of the iterations were
computed simultaneously using a 16,000-way superscalar processor.2 If the compiler has exibility in the
order in which it can execute the instructions that make up your program, it can execute those instructions
simultaneously when parallel hardware is available.
One technique that computer scientists use to formally analyze the potential parallelism in an algorithm
is to characterize how quickly it would execute with an innite-way superscalar processor.
Not all loops contain as much parallelism as this simple loop. We need to identify the things that limit
the parallelism in our codes and remove them whenever possible. In previous chapters we have already
looked at removing clutter and rewriting loops to simplify the body of the loop.
This chapter also supplements Section 2.1.1, in many ways. We looked at the mechanics of compiling
code, all of which apply here, but we didn't answer all of the whys. Basic block analysis techniques form
the basis for the work the compiler does when looking for more parallelism. Looking at two pieces of data,
instructions, or data and instructions, a modern compiler asks the question, Do these things depend on
each other? The three possible answers are yes, no, and we don't know. The third answer is eectively the
same as a yes, because a compiler has to be conservative whenever it is unsure whether it is safe to tweak
the ordering of instructions.
Helping the compiler recognize parallelism is one of the basic approaches specialists take in tuning code.
A slight rewording of a loop or some supplementary information supplied to the compiler can change a we
don't know answer into an opportunity for parallelism. To be certain, there are other facets to tuning
as well, such as optimizing memory access patterns so that they best suit the hardware, or recasting an
algorithm. And there is no single best approach to every problem; any tuning eort has to be a combination
of techniques.
3.1.2 Dependencies3
Imagine a symphony orchestra where each musician plays without regard to the conductor or the other
musicians. At the rst tap of the conductor's baton, each musician goes through all of his or her sheet
music. Some nish far ahead of others, leave the stage, and go home. The cacophony wouldn't resemble
music (come to think of it, it would resemble experimental jazz) because it would be totally uncoordinated.
Of course this isn't how music is played. A computer program, like a musical piece, is woven on a fabric
that unfolds in time (though perhaps woven more loosely). Certain things must happen before or along with
others, and there is a rate to the whole process.
With computer programs, whenever event A must occur before event B can, we say that B is dependent
on A. We call the relationship between them a dependency. Sometimes dependencies exist because of
calculations or memory operations; we call these data dependencies
. Other times, we are waiting for a
branch or do-loop exit to take place; this is called a control dependency
. Each is present in every program
to varying degrees. The goal is to eliminate as many dependencies as possible. Rearranging a program so
that two chunks of the computation are less dependent exposes parallelism
, or opportunities to do several
things at once.
When calculations occur as a consequence of the ow of control, we say there is a control dependency
,
as in the code below and shown graphically in Figure 3.1 (Control dependency). The assignment located
inside the block-if may or may not be executed, depending on the outcome of the test X .NE. 0. In other
words, the value of Y depends on the ow of control in the code around it. Again, this may sound to you
like a concern for compiler designers, not programmers, and that's mostly true. But there are times when
you might want to move control-dependent instructions around to get expensive calculations out of the way
(provided your compiler isn't smart enough to do it for you). For example, say that Figure 3.2 (A little
section of your program) represents a little section of your program. Flow of control enters at the top and
goes through two branch decisions. Furthermore, say that there is a square root operation at the entry point,
and that the ow of control almost always goes from the top, down to the leg containing the statement A=0.0.
This means that the results of the calculation A=SQRT(B) are almost always discarded because A gets a new
value of 0.0 each time through. A square root operation is always expensive because it takes a lot of time
to execute. The trouble is that you can't just get rid of it; occasionally it's needed. However, you could
move it out of the way and continue to observe the control dependencies by making two copies of the square
root operation along the less traveled branches, as shown in Figure 3.3 (Expensive operation moved so that
it's rarely executed). This way the SQRT would execute only along those paths where it was actually needed.
Control dependency
Figure 3.1
Figure 3.2
This kind of instruction scheduling will be appearing in compilers (and even hardware) more and more
as time goes on. A variation on this technique is to calculate results that might be needed at times when
there is a gap in the instruction stream (because of dependencies), thus using some spare cycles that might
otherwise be wasted.
Figure 3.3
A = X + Y + COS(Z)
B = A * C
This dependency is easy to recognize, but others are not so simple. At other times, you must be careful
not to rewrite a variable with a new value before every other computation has nished using the old value.
We can group all data dependencies into three categories: (1) ow dependencies, (2) antidependencies,
and (3) output dependencies. Figure 3.4 (Types of data dependencies) contains some simple examples to
demonstrate each type of dependency. In each example, we use an arrow that starts at the source of the
dependency and ends at the statement that must be delayed by the dependency. The key problem in each
of these dependencies is that the second statement can't execute until the rst has completed. Obviously in
the particular output dependency example, the rst computation is dead code and can be eliminated unless
there is some intervening code that needs the values. There are other techniques to eliminate either output
or antidependencies. The following example contains a ow dependency followed by an output dependency.
Figure 3.4
X = A / B
Y = X + 2.0
X = D - E
While we can't eliminate the ow dependency, the output dependency can be eliminated by using a scratch
variable:
Xtemp = A/B
Y = Xtemp + 2.0
X = D - E
As the number of statements and the interactions between those statements increase, we need a better way
to identify and process these dependencies. Figure 3.5 (Multiple dependencies) shows four statements with
four dependencies.
Multiple dependencies
Figure 3.5
None of the second through fourth instructions can be started before the rst instruction completes.
Figure 3.6
For a basic block of code, we build our DAG in the order of the instructions. The DAG for the previous
four instructions is shown in Figure 3.7 (A more complex data ow graph). This particular example has
many dependencies, so there is not much opportunity for parallelism. Figure 3.8 (Extracting parallelism from
a DAG) shows a more straightforward example shows how constructing a DAG can identify parallelism.
From this DAG, we can determine that instructions 1 and 2 can be executed in parallel. Because we
see the computations that operate on the values A and B while processing instruction 4, we can eliminate a
common subexpression during the construction of the DAG. If we can determine that Z is the only variable
that is used outside this small block of code, we can assume the Y computation is dead code.
Figure 3.7
By constructing the DAG, we take a sequence of instructions and determine which must be executed in a
particular order and which can be executed in parallel. This type of data ow analysis is very important in
the codegeneration phase on super-scalar processors. We have introduced the concept of dependencies and
how to use data ow to nd opportunities for parallelism in code sequences within a basic block. We can
also use data ow analysis to identify dependencies, opportunities for parallelism, and dead code between
basic blocks.
Figure 3.8
To illustrate, suppose that we have the ow graph in Figure 3.9 (Flow graph for data ow analysis).
Beside each basic block we've listed the variables it uses and the variables it denes. What can data ow
analysis tell us?
Notice that a value for A is dened in block X but only used in block Y. That means that A is dead upon
exit from block Y or immediately upon taking the right-hand branch leaving X; none of the other basic blocks
uses the value of A. That tells us that any associated resources, such as a register, can be freed for other
uses.
Looking at Figure 3.9 (Flow graph for data ow analysis) we can see that D is dened in basic block X,
but never used. This means that the calculations dening D can be discarded.
Something interesting is happening with the variable G. Blocks X and W both use it, but if you look
closely you'll see that the two uses are distinct from one another, meaning that they can be treated as two
independent variables.
A compiler featuring advanced instruction scheduling techniques might notice that W is the only block
that uses the value for E, and so move the calculations dening E out of block Y and into W, where they are
needed.
Figure 3.9
In addition to gathering data about variables, the compiler can also keep information about subexpres-
sions. Examining both together, it can recognize cases where redundant calculations are being made (across
basic blocks), and substitute previously computed values in their place. If, for instance, the expression H*I
appears in blocks X, Y, and W, it could be calculated just once in block X and propagated to the others that
use it.
3.1.3 Loops5
Loops are the center of activity for many applications, so there is often a high payback for simplifying or
moving calculations outside, into the computational suburbs. Early compilers for parallel architectures used
pattern matching to identify the bounds of their loops. This limitation meant that a hand-constructed
loop using if-statements and goto-statements would not be correctly identied as a loop. Because modern
compilers use data ow graphs, it's practical to identify loops as a particular subset of nodes in the ow graph.
To a data ow graph, a hand constructed loop looks the same as a compiler-generated loop. Optimizations
can therefore be applied to either type of loop.
Once we have identied the loops, we can apply the same kinds of data-ow analysis we applied above.
Among the things we are looking for are calculations that are unchanging within the loop and variables that
change in a predictable (linear) fashion from iteration to iteration.
5 This content is available online at <http://cnx.org/content/m32784/1.3/>.
• A given node has to dominate all other nodes within the suspected loop. This means that all paths to
any node in the loop have to pass through one particular node, the dominator. The dominator node
forms the header at the top of the loop.
• There has to be a cycle in the graph. Given a dominator, if we can nd a path back to it from one of
the nodes it dominates, we have a loop. This path back is known as the back edge
of the loop.
The ow graph in Figure 3.10 (Flowgraph with a loop in it) contains one loop and one red herring. You can
see that node B dominates every node below it in the subset of the ow graph. That satises Condition 1
(list, p. 134) and makes it a candidate for a loop header. There is a path from E to B, and B dominates E,
so that makes it a back edge, satisfying Condition 2 (list, p. 134). Therefore, the nodes B, C, D, and E form
a loop. The loop goes through an array of linked list start pointers and traverses the lists to determine the
total number of nodes in all lists. Letters to the extreme right correspond to the basic block numbers in the
ow graph.
Figure 3.10
At rst glance, it appears that the nodes C and D form a loop too. The problem is that C doesn't dominate
D (and vice versa), because entry to either can be made from B, so condition 1 (list, p. 134) isn't satised.
Generally, the ow graphs that come from code segments written with even the weakest appreciation for a
structured design oer better loop candidates.
After identifying a loop, the compiler can concentrate on that portion of the ow graph, looking for
instructions to remove or push to the outside. Certain types of subexpressions, such as those found in array
index expressions, can be simplied if they change in a predictable fashion from one iteration to the next.
In the continuing quest for parallelism, loops are generally our best sources for large amounts of paral-
lelism. However, loops also provide new opportunities for those parallelism-killing dependencies.
DO I=1,N
A(I) = A(I) + B(I)
ENDDO
For any two values of I and K, can we calculate the value of A(I) and A(K) at the same time? Below, we
have manually unrolled several iterations of the previous loop, so they can be executed together:
You can see that none of the results are used as an operand for another calculation. For instance, the
calculation for A(I+1) can occur at the same time as the calculation for A(I) because the calculations are
independent; you don't need the results of the rst to determine the second. In fact, mixing up the order of
the calculations won't change the results in the least. Relaxing the serial order imposed on these calculations
makes it possible to execute this loop very quickly on parallel hardware.
DO I=2,N
A(I) = A(I-1) + B(I)
ENDDO
In this case, there is a dependency problem. The value of A(I+1) depends on the value of A(I), the value of
A(I+2) depends on A(I+1), and so on; every iteration depends on the result of a previous one. Dependencies
that extend back to a previous calculation and perhaps a previous iteration (like this one), are loop carried
ow dependencies backward dependencies
or . You often see such dependencies in applications that perform
Gaussian elimination on certain types of matrices, or numerical solutions to systems of dierential equations.
However, it is impossible to run such a loop in parallel (as written); the processor must wait for intermediate
results before it can proceed.
In some cases, ow dependencies are impossible to x; calculations are so dependent upon one another
that we have no choice but to wait for previous ones to complete. Other times, dependencies are a function
of the way the calculations are expressed. For instance, the loop above can be changed to reduce the
dependency. By replicating some of the arithmetic, we can make it so that the second and third iterations
depend on the rst, but not on one another. The operation count goes up we have an extra addition that
we didn't have before but we have reduced the dependency between iterations:
DO I=2,N,2
A(I) = A(I-1) + B(I)
A(I+1) = A(I-1) + B(I) + B(I+1)
ENDDO
The speed increase on a workstation won't be great (most machines run the recast loop more slowly).
However, some parallel computers can trade o additional calculations for reduced dependency and chalk
up a net win.
3.1.4.2 Antidependencies
It's a dierent story when there is a loop-carried antidependency, as in the code below:
DO I=1,N
A(I) = B(I) * E
B(I) = A(I+2) * C
ENDDO
In this loop, there is an antidependency between the variable A(I) and the variable A(I+2). That is, you
must be sure that the instruction that uses A(I+2) does so before the previous one redenes it. Clearly, this
is not a problem if the loop is executed serially, but remember, we are looking for opportunities to overlap
instructions. Again, it helps to pull the loop apart and look at several iterations together. We have recast
the loop by making many copies of the rst statement, followed by copies of the second:
A(I) = B(I) * E
A(I+1) = B(I+1) * E
A(I+2) = B(I+2) * E
...
B(I) = A(I+2) * C ← assignment makes use of the new
B(I+1) = A(I+3) * C value of A(I+2) incorrect.
B(I+2) = A(I+4) * C
The reference to A(I+2) needs to access an old value, rather than one of the new ones being calculated.
If you perform all of the rst statement followed by all of the second statement, the answers will be wrong.
If you perform all of the second statement followed by all of the rst statement, the answers will also be
wrong. In a sense, to run the iterations in parallel, you must either save the A values to use for the second
statement or store all of the B value in a temporary area until the loop completes.
We can also directly unroll the loop and nd some
parallelism:
1 A(I) = B(I) * E
2 B(I) = A(I+2) * C →
3 A(I+1) = B(I+1) * E | Output dependency
4 B(I+1) = A(I+3) * C |
5 A(I+2) = B(I+2) * E ←
6 B(I+2) = A(I+4) * C
Statements 14 could all be executed simultaneously. Once those statements completed execution, statements
58 could execute in parallel. Using this approach, there are sucient intervening statements between the
dependent statements that we can see some parallel performance improvement from a superscalar RISC
processor.
DO I=1,N
A(I) = C(I) * 2.
A(I+2) = D(I) + E
ENDDO
A(I) = C(I) * 2.
A(I+1) = C(I+1) * 2.
A(I+2) = C(I+2) * 2.
A(I+2) = D(I) + E ← Output dependency violated
A(I+3) = D(I+1) + E
A(I+4) = D(I+2) + E
Whether or not you have to worry about output dependencies depends on whether you are actually paral-
lelizing the code. Your compiler will be conscious of the danger, and will be able to generate legal code
and possibly even fast code, if it's clever enough. But output dependencies occasionally become a problem
for programmers.
DO I = 1,N
D = B(I) * 17
A(I) = D + 14
ENDDO
When we look at the loop, the variable D has a ow dependency. The second statement cannot start until
the rst statement has completed. At rst glance this might appear to limit parallelism signicantly. When
we look closer and manually unroll several iterations of the loop, the situation gets worse:
D = B(I) * 17
A(I) = D + 14
D = B(I+1) * 17
A(I+1) = D + 14
D = B(I+2) * 17
A(I+2) = D + 14
Now, the variable D has ow, output, and antidependencies. It looks like this loop has no hope of running in
parallel. However, there is a simple solution to this problem at the cost of some extra memory space, using
a technique called promoting a scalar to a vector . We dene D as an array withN elements and rewrite the
code as follows:
DO I = 1,N
D(I) = B(I) * 17
A(I) = D(I) + 14
ENDDO
Now the iterations are all independent and can be run in parallel. Within each iteration, the rst statement
must run before the second statement.
3.1.4.5 Reductions
The sum of an array of numbers is one example of a reduction
so called because it reduces a vector to a
scalar. The following loop to determine the total of the values in an array certainly looks as though it might
be able to be run in parallel:
SUM = 0.0
DO I=1,N
SUM = SUM + A(I)
ENDDO
This loop also has all three types of dependencies and looks impossible to parallelize. If we are willing to
accept the potential eect of rounding, we can add some parallelism to this loop as follows (again we did not
add the preconditioning loop):
SUM0 = 0.0
SUM1 = 0.0
SUM2 = 0.0
SUM3 = 0.0
DO I=1,N,4
SUM0 = SUM0 + A(I)
SUM1 = SUM1 + A(I+1)
SUM2 = SUM2 + A(I+2)
SUM3 = SUM3 + A(I+3)
ENDDO
Again, this is not precisely the same computation, but all four partial sums can be computed independently.
The partial sums are combined at the end of the loop.
Loops that look for the maximum or minimum elements in an array, or multiply all the elements of an
array, are also reductions. Likewise, some of these can be reorganized into partial results, as with the sum, to
expose more computations. Note that the maximum and minimum are associative operators, so the results
of the reorganized loop are identical to the sequential loop.
DO I=1,N
A(I) = B(I) * E
B(I) = A(I+2) * C
ENDDO
Because each variable reference is solely a function of the index, I, it's clear what kind of dependency we
are dealing with. Furthermore, we can describe how far apart (in iterations) a variable reference is from
its denition. This is called the dependency distance
. A negative value represents a ow dependency; a
positive value means there is an antidependency. A value of zero says that no dependency exists between
the reference and the denition. In this loop, the dependency distance for A is +2 iterations.
However, array subscripts may be functions of other variables besides the loop index. It may be dicult
to tell the distance between the use and denition of a particular element. It may even be impossible to tell
whether the dependency is a ow dependency or an antidependency, or whether a dependency exists at all.
Consequently, it may be impossible to determine if it's safe to overlap execution of dierent statements, as
in the following loop:
DO I=1,N
A(I) = B(I) * E
B(I) = A(I+K) * C ← K unknown
ENDDO
If the loop made use of A(I+K), where the value of K was unknown, we wouldn't be able to tell (at least
by looking at the code) anything about the kind of dependency we might be facing. If K is zero, we have a
dependency within the iteration and no loop-carried dependencies. If K is positive we have an antidependency
with distance K. Depending on the value for K, we might have enough parallelism for a superscalar processor.
If K is negative, we have a loop-carried ow dependency, and we may have to execute the loop serially.
Ambiguous references, like A(I+K) above, have an eect on the parallelism we can detect in a loop. From
the compiler perspective, it may be that this loop does contain two independent calculations that the author
7 This content is available online at <http://cnx.org/content/m32788/1.3/>.
whimsically decided to throw into a single loop. But when they appear together, the compiler has to treat
them conservatively, as if they were interrelated. This has a big eect on performance. If the compiler has
to assume that consecutive memory references may ultimately access the same location, the instructions
involved cannot be overlapped. One other option is for the compiler to generate two versions of the loop
and check the value for K at runtime to determine which version of the loop to execute.
A similar situation occurs when we use integer index arrays in a loop. The loop below contains only a
single statement, but you can't be sure that any iteration is independent without knowing the contents of
the K and J arrays:
DO I=1,N
A(K(I)) = A(K(I)) + B(J(I)) * C
ENDDO
For instance, what if all of the values for K(I) were the same? This causes the same element of the array A
to be rereferenced with each iteration! That may seem ridiculous to you, but the compiler can't tell.
With code like this, it's common for every value of K(I) to be unique. This is called a permutation. If
you can tell a compiler that it is dealing with a permutation, the penalty is lessened in some cases. Even so,
there is insult being added to injury. Indirect references require more memory activity than direct references,
and this slows you down.
C compilers don't enjoy the same restrictions on aliasing. In fact, there are cases where aliasing could be
desirable. Additionally, C is blessed with pointer types, increasing the opportunities for aliasing to occur.
This means that a C compiler has to approach operations through pointers more conservatively than a
FORTRAN compiler would. Let's look at some examples to see why.
The following loop nest looks like a FORTRAN loop cast in C. The arrays are declared or allocated all at
once at the top of the routine, and the starting address and leading dimensions are visible to the compiler.
This is important because it means that the storage relationship between the array elements is well known.
Hence, you could expect good performance:
#define N ...
double *a[N][N], c[N][N], d;
for (i=0; i<N; i++)
Now imagine what happens if you allocate the rows dynamically. This makes the address calculations more
complicated. The loop nest hasn't changed; however, there is no guaranteed stride that can get you from
one row to the next. This is because the storage relationship between the rows is unknown:
#define N ...
double *a[N], *c[N], d;
for (i=0; i<N; i++) {
a[i] = (double *) malloc (N*sizeof(double));
c[i] = (double *) malloc (N*sizeof(double));
}
for (i=0; i<N; i++)
for (j=0; j<N; j++)
a[i][j] = a[i][j] + c[j][i] * d;
In fact, your compiler knows even less than you might expect about the storage relationship. For instance,
how can it be sure that references to a and c aren't aliases? It may be obvious to you that they're not. You
might point out that malloc never overlaps storage. But the compiler isn't free to assume that. Who knows?
You may be substituting your own version of malloc
!
Let's look at a dierent example, where storage is allocated all at once, though the declarations are not
visible to all routines that are using it. The following subroutine bob performs the same computation as
our previous example. However, because the compiler can't see the declarations for a and c (they're in the
main routine), it doesn't have enough information to be able to overlap memory references from successive
iterations; the references could be aliases:
#define N...
main()
{
double a[N][N], c[N][N], d;
...
bob (a,c,d,N);
}
bob (double *a,double *c,double d,int n)
{
int i,j;
double *ap, *cp;
for (i=0;i<n;i++) {
ap = a + (i*n);
cp = c + i;
for (j=0; j<n; j++)
*(ap+j) = *(ap+j) + *(cp+(j*n)) * d;
}
}
To get the best performance, make available to the compiler as many details about the size and shape of your
data structures as possible. Pointers, whether in the form of formal arguments to a subroutine or explicitly
declared, can hide important facts about how you are using memory. The more information the compiler
has, the more it can overlap memory references. This information can come from compiler directives or from
making declarations visible in the routines where performance is most critical.
In the coming chapters, we will begin to learn more about executing our programs on parallel multiprocessors.
At some point we will escape the bonds of compiler automatic optimization and begin to explicitly code the
parallel portions of our code.
To learn more about compilers and dataow, read The Art of Compiler Design: Theory and Practice by
Thomas Pittman and James Peters (Prentice-Hall).
3.1.7 Exercises9
Exercise 3.1.7.1
Identify the dependencies (if there are any) in the following loops. Can you think of ways to
organize each loop for more parallelism?
a.
DO I=1,N-2
A(I+2) = A(I) + 1.
ENDDO
b.
DO I=1,N-1,2
A(I+1) = A(I) + 1.
ENDDO
c.
8 This content is available online at <http://cnx.org/content/m32789/1.3/>.
9 This content is available online at <http://cnx.org/content/m32792/1.3/>.
DO I=2,N
A(I) = A(I-1) * 2.
B = A(I-1)
ENDDO
d.
DO I=1,N
IF(N .GT. M)
A(I) = 1.
ENDDO
e.
DO I=1,N
A(I,J) = A(I,K) + B
ENDDO
f.
DO I=1,N-1
A(I+1,J) = A(I,K) + B
ENDDO
g.
Exercise 3.1.7.2
Imagine that you are a parallelizing compiler, trying to generate code for the loop below. Why are
references to A a challenge? Why would it help to know that K is equal to zero? Explain how you
could partially vectorize the statements involving A if you knew that K had an absolute value of
at least 8.
DO I=1,N
Exercise 3.1.7.3
The following three statements contain a ow dependency, an antidependency and an output
dependency. Can you identify each? Given that you are allowed to reorder the statements, can you
nd a permutation that produces the same values for the variables C and B? Show how you can
reduce the dependencies by combining or rearranging calculations and using temporary variables.
B = A + C
B = C + D
C = B + D
A shared-memory multiprocessor
Figure 3.11
Figure 3.12
A crossbar is a hardware approach to eliminate the bottleneck caused by a single bus. A crossbar is
like several buses running side by side with attachments to each of the modules on the machine CPU,
memory, and peripherals. Any module can get to any other by a path through the crossbar, and multiple
paths may be active simultaneously. In the 4×5 crossbar of Figure 3.13 (A crossbar), for instance, there can
be four active data transfers in progress at one time. In the diagram it looks like a patchwork of wires, but
there is actually quite a bit of hardware that goes into constructing a crossbar. Not only does the crossbar
connect parties that wish to communicate, but it must also actively arbitrate between two or more CPUs
that want access to the same memory or peripheral. In the event that one module is too popular, it's the
crossbar that decides who gets access and who doesn't. Crossbars have the best performance because there
is no single shared bus. However, they are more expensive to build, and their cost increases as the number
of ports is increased. Because of their cost, crossbars typically are only found at the high end of the price
and performance spectrum.
Whether the system uses a bus or crossbar, there is only so much memory bandwidth to go around; four
or eight processors drawing from one memory system can quickly saturate all available bandwidth. All of
the techniques that improve memory performance (as described in ) also apply here in the design of the
memory subsystems attached to these buses or crossbars.
Figure 3.13
Figure 3.14
In actuality, on some of the fastest bus-based systems, the memory bus is suciently fast that up to
20 processors can access memory using unit stride with very little conict. If the processors are accessing
memory using non-unit stride, bus and memory bank conict becomes apparent, with fewer processors.
This bus architecture combined with local caches is very popular for general-purpose multiprocessing
loads. The memory reference patterns for database or Internet servers generally consist of a combination of
time periods with a small working set, and time periods that access large data structures using unit stride.
Scientic codes tend to perform more non-unit-stride access than general-purpose codes. For this reason,
the most expensive parallel-processing systems targeted at scientic codes tend to use crossbars connected
to multibanked memory systems.
The main memory system is better shielded when a larger cache is used. For this reason, multiprocessors
sometimes incorporate a two-tier cache system, where each processor uses its own small on-chip local cache,
backed up by a larger second board-level cache with as much as 4 MB of memory. Only when neither can
satisfy a memory request, or when data has to be written back to main memory, does a request go out over
the bus or crossbar.
3.2.2.2 Coherency
Now, what happens when one CPU of a multiprocessor running a single program in parallel changes the
value of a variable, and another CPU tries to read it? Where does the value come from? These questions
are interesting because there can be multiple copies of each variable, and some of them can hold old or stale
values.
For illustration, say that you are running a program with a shared variable A. Processor 1 changes the
value of A and Processor 2 goes to read it.
Figure 3.15
In Figure 3.15 (Multiple copies of variable A), if Processor 1 is keeping A as a register-resident variable,
then Processor 2 doesn't stand a chance of getting the correct value when it goes to look for it. There is no
way that 2 can know the contents of 1's registers; so assume, at the very least, that Processor 1 writes the
new value back out. Now the question is, where does the new value get stored? Does it remain in Processor
1's cache? Is it written to main memory? Does it get updated in Processor 2's cache?
Really, we are asking what kind of cache coherency protocol the vendor uses to assure that all processors
see a uniform view of the values in memory. It generally isn't something that the programmer has to
worry about, except that in some cases, it can aect performance. The approaches used in these systems
are similar to those used in single-processor systems with some extensions. The most straight-forward cache
coherency approach is called a write-through policy : variables written into cache are simultaneously written
into main memory. As the update takes place, other caches in the system see the main memory reference
being performed. This can be done because all of the caches continuously monitor (also known as snooping
) the trac on the bus, checking to see if each address is in their cache. If a cache notices that it contains
a copy of the data from the locations being written, it may either invalidate
its copy of the variable or
obtain new values (depending on the policy). One thing to note is that a write-through cache demands a
fair amount of main memory bandwidth since each write goes out over the main memory bus. Furthermore,
successive writes to the same location or bank are subject to the main memory cycle time and can slow the
machine down.
A more sophisticated cache coherency protocol is called copyback writeback
or . The idea is that you
write values back out to main memory only when the cache housing them needs the space for something else.
Updates of cached data are coordinated between the caches, by the caches, without help from the processor.
Copyback caching also uses hardware that can monitor (snoop) and respond to the memory transactions of
the other caches in the system. The benet of this method over the write-through method is that memory
trac is reduced considerably. Let's walk through it to see how it works.
Exclusive: There are no other caches that have this cache line.
Shared: There are read-only copies of this line in two or more caches.
Empty/Invalid: This cache line doesn't contain any useful data.
This particular coherency protocol is often called MESI
. Other cache coherency protocols are more compli-
cated, but these states give you an idea how multiprocessor writeback cache coherency works.
We start where a particular cache line is in memory and in none of the writeback caches on the systems.
The rst cache to ask for data from a particular part of memory completes a normal memory access; the
main memory system returns data from the requested location in response to a cache miss. The associated
cache line is marked exclusive
, meaning that this is the only cache in the system containing a copy of the
data; it is the owner of the data. If another cache goes to main memory looking for the same thing, the
request is intercepted by the rst cache, and the data is returned from the rst cache not main memory.
Once an interception has occurred and the data is returned, the data is marked shared
in both of the caches.
When a particular line is marked shared, the caches have to treat it dierently than they would if they
were the exclusive owners of the data especially if any of them wants to modify it. In particular, a write to
a shared cache entry is preceded by a broadcast message to all the other caches in the system. It tells them
to invalidate their copies of the data. The one remaining cache line gets marked as modied
to signal that
it has been changed, and that it must be returned to main memory when the space is needed for something
else. By these mechanisms, you can maintain cache coherence across the multiprocessor without adding
tremendously to the memory trac.
By the way, even if a variable is not shared, it's possible for copies of it to show up in several caches.
On a symmetric multiprocessor, your program can bounce around from CPU to CPU. If you run for a little
while on this CPU, and then a little while on that, your program will have operated out of separate caches.
That means that there can be several copies of seemingly unshared variables scattered around the machine.
Operating systems often try to minimize how often a process is moved between physical CPUs during context
switches. This is one reason not to overload the available processors in a system.
% ps -a
PID TTY TIME CMD
28410 pts/34 0:00 tcsh
28213 pts/38 0:00 xterm
10488 pts/51 0:01 telnet
28411 pts/34 0:00 xbiff
11123 pts/25 0:00 pine
3805 pts/21 0:00 elm
6773 pts/44 5:48 ansys
...
% ps --a | grep ansys
6773 pts/44 6:00 ansys
For each process we see the process identier (PID), the terminal that is executing the command, the amount
of CPU time the command has used, and the name of the command. The PID is unique across the entire
system. Most UNIX commands are executed in a separate process. In the above example, most of the
processes are waiting for some type of event, so they are taking very few resources except for memory.
Process 677313 seems to be executing and using resources. Running ps
again conrms that the CPU time is
increasing for the ansysprocess:
% vmstat 5
procs memory page disk faults cpu
r b w swap free re mf pi po fr de sr f0 s0 -- -- in sy cs us sy id
3 0 0 353624 45432 0 0 1 0 0 0 0 0 0 0 0 461 5626 354 91 9 0
3 0 0 353248 43960 0 22 0 0 0 0 0 0 14 0 0 518 6227 385 89 11 0
Running the vmstat 5command tells us many things about the activity on the system. First, there are
three runnable processes. If we had one CPU, only one would actually be running at a given instant. To
allow all three jobs to progress, the operating system time-shares between the processes. Assuming equal
priority, each process executes about 1/3 of the time. However, this system is a two-processor system, so
each process executes about 2/3 of the time. Looking across the vmstat
output, we can see paging activity
pi po cs us sy
( , ), context switches ( ), overall user time ( ), system time ( ), and idle time ( ). id
Each process can execute a completely dierent program. While most processes are completely inde-
pendent, they can cooperate and share information using interprocess communication (pipes, sockets) or
various operating system-supported shared-memory areas. We generally don't use multiprocessing on these
shared-memory systems as a technique to increase single-application performance.
13 ANSYS is a commonly used structural-analysis package.
main () {
int pid,status,retval;
int stackvar; /* A stack variable */
globvar = 1;
stackvar = 1;
printf("Main - calling fork globvar=%d stackvar=%d\n",globvar,stackvar);
pid = fork();
printf("Main - fork returned pid=%d\n",pid);
if ( pid == 0 ) {
printf("Child - globvar=%d stackvar=%d\n",globvar,stackvar);
sleep(1);
printf("Child - woke up globvar=%d stackvar=%d\n",globvar,stackvar);
globvar = 100;
stackvar = 100;
printf("Child - modified globvar=%d stackvar=%d\n",globvar,stackvar);
retval = execl("/bin/date", (char *) 0 );
printf("Child - WHY ARE WE HERE retval=%d\n",retval);
} else {
printf("Parent - globvar=%d stackvar=%d\n",globvar,stackvar);
globvar = 5;
stackvar = 5;
printf("Parent - sleeping globvar=%d stackvar=%d\n",globvar,stackvar);
sleep(2);
printf("Parent - woke up globvar=%d stackvar=%d\n",globvar,stackvar);
printf("Parent - waiting for pid=%d\n",pid);
retval = wait(&status);
status = status 8; /* Return code in bits 15-8 */
printf("Parent - status=%d retval=%d\n",status,retval);
}
}
The key to understanding this code is to understand how the fork( ) function operates. The simple
summary is that the fork( ) function is called once in a process and returns twice, once in the original
process and once in a newly created process. The newly created process is an identical copy of the original
process. All the variables (local and global) have been duplicated. Both processes have access to all of the
open les of the original process. Figure 3.16 (How a fork operates) shows how the fork operation creates a
new process.
14 These examples are written in C using the POSIX 1003.1 application programming interface. This example runs on most
UNIX systems and on other POSIX-compliant systems including OpenNT, Open- VMS, and many others.
Tracing this through, rst the program sets the global and stack variable to one and then calls fork( ).
During the fork( ) call, the operating system suspends the process, makes an exact duplicate of the process,
and then restarts both processes. You can see two messages from the statement immediately after the fork.
The rst line is coming from the original process, and the second line is coming from the new process. If you
were to execute a ps command at this moment in time, you would see two processes running called fork.
One would have a process identier of 19336.
Figure 3.16
As both processes start, they execute an IF-THEN-ELSE and begin to perform dierent actions in the
parent and child. Notice that globvar
and stackvar
are set to 5 in the parent, and then the parent sleeps for
two seconds. At this point, the child begins executing. The values for andglobvar stackvar
are unchanged
in the child process. This is because these two processes are operating in completely independent memory
spaces. The child process sleeps for one second and sets its copies of the variables to 100. Next, the child
process calls the execl( ) function to overwrite its memory space with the UNIX date program. Note that
the execl( ) never returns; the date program takes over all of the resources of the child process. If you
were to do a ps
at this moment in time, you still see two processes on the system but process 19336 would
be called date. The date command executes, and you can see its output.15
The parent wakes up after a brief two-second sleep and notices that its copies of global and local variables
have not been changed by the action of the child process. The parent then calls the wait( ) function to
15 It's not uncommon for a human parent process to fork and create a human child process that initially seems to have the
same identity as the parent. It's also not uncommon for the child process to change its overall identity to be something very
dierent from the parent at some later point. Usually human children wait 13 years or so before this change occurs, but in
UNIX, this happens in a few microseconds. So, in some ways, in UNIX, there are many parent processes that are disappointed
because their children did not turn out like them!
#define THREAD_COUNT 3
void *TestFunc(void *);
int globvar; /* A global variable */
int index[THREAD_COUNT] /* Local zero-based thread index */
pthread_t thread_id[THREAD_COUNT]; /* POSIX Thread IDs */
main() {
int i,retval;
pthread_t tid;
globvar = 0;
printf("Main - globvar=%d\n",globvar);
for(i=0;i<THREAD_COUNT;i++) {
index[i] = i;
retval = pthread_create(&tid,NULL,TestFunc,(void *) index[i]);
printf("Main - creating i=%d tid=%d retval=%d\n",i,tid,retval);
thread_id[i] = tid;
}
16 This example uses the IEEE POSIX standard interface for a thread library. If your system supports POSIX threads, this
example should work. If not, there should be similar routines on your system for each of the thread functions.
Figure 3.17
The global shared areas in this case are those variables declared in the static area outside the main( ) code.
The local variables are any variables declared within a routine. When threads are added, each thread gets
its own function call stack. In C, the automatic variables that are declared at the beginning of each routine
are allocated on the stack. As each thread enters a function, these variables are separately allocated on that
particular thread's stack. So these are the thread-local variables.
Unlike the fork( ) function, the pthread_create( ) function creates a new thread, and then control is
returned to the calling thread. One of the parameters of the pthread_create( ) is the name of a function.
New threads begin execution in the function TestFunc( ) and the thread nishes when it returns from
this function. When this program is executed, it produces the following output:
Main - globvar=0
Main - creating i=0 tid=4 retval=0
Main - creating i=1 tid=5 retval=0
Main - creating i=2 tid=6 retval=0
Main thread - threads started globvar=0
Main - waiting for join 4
TestFunc me=0 - self=4 globvar=0
TestFunc me=0 - sleeping globvar=15
TestFunc me=1 - self=5 globvar=15
TestFunc me=1 - sleeping globvar=16
TestFunc me=2 - self=6 globvar=16
TestFunc me=2 - sleeping globvar=17
TestFunc me=2 - done param=6 globvar=17
TestFunc me=1 - done param=5 globvar=17
TestFunc me=0 - done param=4 globvar=17
Main - back from join 0 retval=0
Main - waiting for join 5
Main - back from join 1 retval=0
Main - waiting for join 6
Main - back from join 2 retval=0
Main thread -- threads completed globvar=17
recs %
You can see the threads getting created in the loop. The master thread completes the pthread_create( )
loop, executes the second loop, and calls the pthread_join( ) function. This function suspends the master
thread until the specied thread completes. The master thread is waiting for Thread 4 to complete. Once
the master thread suspends, one of the new threads is started. Thread 4 starts executing. Initially the
,
variable globvar is set to 0 from the main program. The self, me and param variables are thread-local
variables, so each thread has its own copy. Thread 4 sets globvar to 15 and goes to sleep. Then Thread 5
begins to execute and sees globvar set to 15 from Thread 4; Thread 5 sets globvar to 16, and goes to sleep.
This activates Thread 6, which sees the current value for globvar and sets it to 17. Then Threads 6, 5, and
4 wake up from their sleep, all notice the latest value of 17 in globvar, and return from the TestFunc( )
routine, ending the threads.
All this time, the master thread is in the middle of a pthread_join( ) waiting for Thread 4 to complete.
As Thread 4 completes, the pthread_join( ) returns. The master thread then calls pthread_join( )
repeatedly to ensure that all three threads have been completed. Finally, the master thread prints out the
value for globvar that contains the latest value of 17.
To summarize, when an application is executing with more than one thread, there are shared global areas
and thread private areas. Dierent threads execute at dierent times, and they can easily work together in
shared areas.
threads.
We can explore this eect by substituting the following SpinFunc( ) function, replacing TestFunc( )
function in the pthread_create( ) call in the previous example:
recs % ps
PID TTY TIME CMD
23921 pts/35 0:09 create2
recs % ps
PID TTY TIME CMD
17 The pthreads library supports both user-space threads and operating-system threads, as we shall soon see. Another popular
early threads package was called cthreads.
We run the program in the background18 and everything seems to run ne. All the threads go to sleep for
1, 2, and 3 seconds. The rst thread wakes up and starts the loop waiting for globvar to be incremented by
the other threads. Unfortunately, with user space threads, there is no automatic time sharing. Because we
are in a CPU loop that never makes a system call, the second and third threads never get scheduled so they
can complete their sleep( ) call. To x this problem, we need to make the following change to the code:
With this sleep19 call, Threads 2 and 3 get a chance to be scheduled. They then nish their sleep calls,
increment the globvar variable, and the program terminates properly.
You might ask the question, Then what is the point of user space threads? Well, when there is a high
performance database server or Internet server, the multiple logical threads can overlap network I/O with
database I/O and other background computations. This technique is not so useful when the threads all
want to perform simultaneous CPU-intensive computations. To do this, you need threads that are created,
managed, and scheduled by the operating system rather than a user library.
#define THREAD_COUNT 2
void *SpinFunc(void *);
int globvar; /* A global variable */
int index[THREAD_COUNT]; /* Local zero-based thread index */
pthread_t thread_id[THREAD_COUNT]; /* POSIX Thread IDs */
pthread_attr_t attr; /* Thread attributes NULL=use default */
18 Because we know it will hang and ignore interrupts.
19 Some thread libraries support a call to a routine sched_yield( ) that checks for runnable threads. If it nds a runnable
thread, it runs the thread. If no thread is runnable, it returns immediately to the calling thread. This routine allows a thread
that has the CPU to ensure that other threads make progress during CPU-intensive periods of its code.
globvar = 0;
pthread_attr_init(&attr); /* Initialize attr with defaults */
pthread_attr_setscope(&attr, PTHREAD_SCOPE_SYSTEM);
printf("Main - globvar=%d\n",globvar);
for(i=0;i<THREAD_COUNT;i++) {
index[i] = i;
retval = pthread_create(&tid,&attr,SpinFunc,(void *) index[i]);
printf("Main - creating i=%d tid=%d retval=%d\n",i,tid,retval);
thread_id[i] = tid;
}
printf("Main thread - threads started globvar=%d\n",globvar);
for(i=0;i<THREAD_COUNT;i++) {
printf("Main - waiting for join %d\n",thread_id[i]);
retval = pthread_join( thread_id[i], NULL ) ;
printf("Main - back from join %d retval=%d\n",i,retval);
}
printf("Main thread - threads completed globvar=%d\n",globvar);
}
The code executed by the master thread is modied slightly. We create an attribute data structure and
set the PTHREAD_SCOPE_SYSTEM attribute to indicate that we would like our new threads to be created and
scheduled by the operating system. We use the attribute information on the call to pthread_create( ).
None of the other code has been changed. The following is the execution output of this new program:
recs % create3
Main - globvar=0
Main - creating i=0 tid=4 retval=0
SpinFunc me=0 - sleeping 1 seconds ...
Main - creating i=1 tid=5 retval=0
Main thread - threads started globvar=0
Main - waiting for join 4
SpinFunc me=1 - sleeping 2 seconds ...
SpinFunc me=0 - wake globvar=0...
SpinFunc me=0 - spinning globvar=1...
SpinFunc me=1 - wake globvar=1...
SpinFunc me=1 - spinning globvar=2...
SpinFunc me=1 - done globvar=2...
SpinFunc me=0 - done globvar=2...
Main - back from join 0 retval=0
Main - waiting for join 5
Main - back from join 1 retval=0
Main thread - threads completed globvar=2
recs %
Now the program executes properly. When the rst thread starts spinning, the operating system is context
switching between all three threads. As the threads come out of their sleep( ), they increment their shared
variable, and when the nal thread increments the shared variable, the other two threads instantly notice
the new value (because of the cache coherency protocol) and nish the loop. If there are fewer than three
CPUs, a thread may have to wait for a time-sharing context switch to occur before it notices the updated
global variable.
With operating-system threads and multiple processors, a program can realistically break up a large
computation between several independent threads and compute the solution more quickly. Of course this
presupposes that the computation could be done in parallel in the rst place.
globvar++;
LOAD R1,globvar
ADD R1,1
STORE R1,globvar
What if globvar contained 0, Thread 1 was running, and, at the precise moment it completed the LOAD into
Register R1 and before it had completed the ADD or STORE instructions, the operating system interrupted the
thread and switched to Thread 2? Thread 2 catches up and executes all three instructions using its registers:
loading 0, adding 1 and storing the 1 back into globvar. Now Thread 2 goes to sleep and Thread 1 is
restarted at the ADD instruction. Register R1 for Thread 1 contains the previously loaded value of 0; Thread
1 adds 1 and then stores 1 into globvar. What is wrong with this picture? We meant to use this code to
count the number of threads that have passed this point. Two threads passed the point, but because of a
bad case of bad timing, our variable indicates only that one thread passed. This is because the increment of
a variable in memory is not atomic. That is, halfway through the increment, something else can happen.
Another way we can have a problem is on a multiprocessor when two processors execute these instructions
simultaneously. They both do the LOAD, getting 0. Then they both add 1 and store 1 back to memory.21
Which processor actually got the honor of storing their
1 back to memory is simply a race.
We must have some way of guaranteeing that only one thread can be in these three instructions at the
same time. If one thread has started these instructions, all other threads must wait to enter until the rst
thread has exited. These areas are called critical sections
. On single-CPU systems, there was a simple
solution to critical sections: you could turn o interrupts for a few instructions and then turn them back
on. This way you could guarantee that you would get all the way through before a timer or other interrupt
occurred:
However, this technique does not work for longer critical sections or when there is more than one CPU. In
these cases, you need a lock, a semaphore, or a mutex. Most thread libraries provide this type of routine.
To use a mutex, we have to make some modications to our example code:
21 Boy, this is getting pretty picky. How often will either of these events really happen? Well, if it crashes your airline
reservation system every 100,000 transactions or so, that would be way too often.
...
pthread_mutex_t my_mutex; /* MUTEX data structure */
...
main() {
...
pthread_attr_init(&attr); /* Initialize attr with defaults */
pthread_mutex_init (&my_mutex, NULL);
.... pthread_create( ... )
...
}
void *SpinFunc(void *parm) {
...
pthread_mutex_lock (&my_mutex);
globvar ++;
pthread_mutex_unlock (&my_mutex);
while(globvar < THREAD_COUNT ) ;
printf("SpinFunc me=%d -- done globvar=%d...\n", me, globvar);
...
}
The mutex data structure must be declared in the shared area of the program. Before the threads are
created, pthread_mutex_init must be called to initialize the mutex. Before globvar is incremented, we
must lock the mutex and after we nish updating globvar (three instructions later), we unlock the mutex.
With the code as shown above, there will never be more than one processor executing the globvar++ line
of code, and the code will never hang because an increment was missed. Semaphores and locks are used in
a similar way.
Interestingly, when using user space threads, an attempt to lock an already locked mutex, semaphore, or
lock can cause a thread context switch. This allows the thread that owns the lock a better chance to make
progress toward the point where they will unlock the critical section. Also, the act of unlocking a mutex can
cause the thread waiting for the mutex to be dispatched by the thread library.
3.2.4.3 Barriers
Barriers are dierent than critical sections. Sometimes in a multithreaded application, you need to have all
threads arrive at a point before allowing any threads to execute beyond that point. An example of this is
a time-based simulation . Each task processes its portion of the simulation but must wait until all of the
threads have completed the current time step before any thread can begin the next time step. Typically
threads are created, and then each thread executes a loop with one or more barriers in the loop. The rough
pseudocode for this type of approach is as follows:
main() {
for (ith=0;ith<NUM_THREADS;ith++) pthread_create(..,work_routine,..)
for (ith=0;ith<NUM_THREADS;ith++) pthread_join(...) /* Wait a long time */
exit()
}
work_routine() {
#define MAX_THREAD 4
void *SumFunc(void *);
int ThreadCount; /* Threads on this try */
double GlobSum; /* A global variable */
int index[MAX_THREAD]; /* Local zero-based thread index */
pthread_t thread_id[MAX_THREAD]; /* POSIX Thread IDs */
pthread_attr_t attr; /* Thread attributes NULL=use default */
pthread_mutex_t my_mutex; /* MUTEX data structure */
main() {
int i,retval;
pthread_t tid;
double single,multi,begtime,endtime;
/* Initialize things */
for (i=0; i<MAX_SIZE; i++) array[i] = drand48();
pthread_attr_init(&attr); /* Initialize attr with defaults */
22 This content is available online at <http://cnx.org/content/m32804/1.3/>.
recs % addup
Single sum=7999998000000.000000 time=0.256624
Threads=2
SumFunc me=0 start=0 end=2000000
SumFunc me=1 start=2000000 end=4000000
Sum=7999998000000.000000 time=0.133530
Efficiency = 0.960923
Threads=3
SumFunc me=0 start=0 end=1333333
SumFunc me=1 start=1333333 end=2666666
SumFunc me=2 start=2666666 end=4000000
Sum=7999998000000.000000 time=0.091018
Efficiency = 0.939829
Threads=4
SumFunc me=0 start=0 end=1000000
SumFunc me=1 start=1000000 end=2000000
SumFunc me=2 start=2000000 end=3000000
SumFunc me=3 start=3000000 end=4000000
Sum=7999998000000.000000 time=0.107473
Efficiency = 0.596950
recs %
There are some interesting patterns. Before you interpret the patterns, you must know that this system is
a three-processor Sun Enterprise 3000. Note that as we go from one to two threads, the time is reduced to
one-half. That is a good result given how much it costs for that extra CPU. We characterize how well the
additional resources have been used by computing an eciency factor that should be 1.0. This is computed
by multiplying the wall time by the number of threads. Then the time it takes on a single processor is divided
by this number. If you are using the extra processors well, this evaluates to 1.0. If the extra processors are
used pretty well, this would be about 0.9. If you had two threads, and the computation did not speed up at
all, you would get 0.5.
At two and three threads, wall time is dropping, and the eciency is well over 0.9. However, at four
threads, the wall time increases, and our eciency drops very dramatically. This is because we now have
more threads than processors. Even though we have four threads that could execute, they must be time-
sliced between three processors.23 This is even worse that it might seem. As threads are switched, they move
from processor to processor and their caches must also move from processor to processor, further slowing
performance. This cache-thrashing eect is not too apparent in this example because the data structure is
so large, most memory references are not to values previously in cache.
It's important to note that because of the nature of oating-point (see Section 1.2.1), the parallel sum
may not be the same as the serial sum. To perform a summation in parallel, you must be willing to tolerate
these slight variations in your results.
23 It is important to match the number of runnable threads to the available resources. In compute code, when there are
more threads than available processors, the threads compete among themselves, causing unnecessary overhead and reducing
the eciency of your computation.
3.2.7 Exercises25
Exercise 3.2.7.1
Experiment with the fork code in this chapter. Run the program multiple times and see how the
order of the messages changes. Explain the results.
Exercise 3.2.7.2
Experiment with the create1 and create3 codes in this chapter. Remove all of the sleep( )
calls. Execute the programs several times on single and multiprocessor systems. Can you explain
why the output changes from run to run in some situations and doesn't change in others?
Exercise 3.2.7.3
Experiment with the parallel sum code in this chapter. In the SumFunc( ) routine, change the
for-loop to:
Remove the three lines at the end that get the mutex and update the GlobSum. Execute the
code. Explain the dierence in values that you see for GlobSum. Are the patterns dierent on a
single processor and a multiprocessor? Explain the performance impact on a single processor and
a multiprocessor.
Exercise 3.2.7.4
Explain how the following code segment could cause deadlock two or more processes waiting
for a resource that can't be relinquished:
...
call lock (lword1)
call lock (lword2)
...
call unlock (lword1)
call unlock (lword2)
.
24 This content is available online at <http://cnx.org/content/m32807/1.3/>.
25 This content is available online at <http://cnx.org/content/m32810/1.3/>.
Exercise 3.2.7.5
If you were to code the functionality of a spin-lock in C, it might look like this:
while (!lockword);
lockword = !lockword;
As you know from the rst sections of the book, the same statements would be compiled into explicit
loads and stores, a comparison, and a branch. There's a danger that two processes could each load
lockword, nd it unset, and continue on as if they owned the lock (we have a race condition). This
suggests that spin-locks are implemented dierently that they're not merely the two lines of C
above. How do you suppose they are implemented?
have to do is turn on a compiler ag and buy a good parallel processor. For example, look at the following
code:
PARAMETER(NITER=300,N=1000000)
REAL*8 A(N),X(N),B(N),C
DO ITIME=1,NITER
DO I=1,N
A(I) = X(I) + B(I) * C
ENDDO
CALL WHATEVER(A,X,B,C)
ENDDO
Here we have an iterative code that satises all the criteria for a good parallel loop. On a good parallel
processor with a modern compiler, you are two ags away from executing in parallel. On Sun Solaris
systems, the autopar ag turns on the automatic parallelization, and the loopinfo ag causes the compiler
to describe the particular optimization performed for each loop. To compile this code under Solaris, you
simply add these ags to your f77 call:
real 30.9
user 30.7
sys 0.1
E6000:
If you simply run the code, it's executed using one thread. However, the code is enabled for parallel processing
for those loops that can be executed in parallel. To execute the code in parallel, you need to set the UNIX
environment to the number of parallel threads you wish to use to execute the code. On Solaris, this is done
using the PARALLEL variable:
real 30.9
user 30.7
sys 0.1
E6000: setenv PARALLEL 2
E6000: /bin/time daxpy
real 8.2
user 32.0
sys 0.5
E6000: setenv PARALLEL 8
E6000: /bin/time daxpy
real 4.3
user 33.0
sys 0.8
Speedup is the term used to capture how much faster the job runs using N processors compared to the
performance on one processor. It is computed by dividing the single processor time by the multiprocessor
time for each number of processors. Figure 3.18 (Improving performance by adding processors) shows the
wall time and speedup for this application.
Figure 3.18
Figure 3.19 (Ideal and actual performance improvement) shows this information graphically, plotting
speedup versus the number of processors.
Figure 3.19
Note that for a while we get nearly perfect speedup, but we begin to see a measurable drop in speedup at
four and eight processors. There are several causes for this. In all parallel applications, there is some portion
of the code that can't run in parallel. During those nonparallel times, the other processors are waiting for
work and aren't contributing to eciency. This nonparallel code begins to aect the overall performance as
more processors are added to the application.
So you say, this is more like it! and immediately try to run with 12 and 16 threads. Now, we see
the graph in Figure 3.21 (Diminishing returns) and the data from Figure 3.20 (Increasing the number of
threads).
Figure 3.20
Diminishing returns
Figure 3.21
What has happened here? Things were going so well, and then they slowed down. We are running this
program on a 16-processor system, and there are eight other active threads, as indicated below:
E6000:uptime
4:00pm up 19 day(s), 37 min(s), 5 users, load average: 8.00, 8.05, 8.14
E6000:
Once we pass eight threads, there are no available processors for our threads. So the threads must be time-
shared between the processors, signicantly slowing the overall operation. By the end, we are executing 16
threads on eight processors, and our performance is slower than with one thread. So it is important that
you don't create too many threads in these types of applications.
• Which loops can execute in parallel, producing the exact same results as the sequential executions of
the loops? This is done by checking for dependencies that span iterations. A loop with no interiteration
dependencies is called a DOALL loop.
• Which loops are worth executing in parallel? Generally very short loops gain no benet and may
execute more slowly when executing in parallel. As with loop unrolling, parallelism always has a cost.
It is best used when the benet far outweighs the cost.
• In a loop nest, which loop is the best candidate to be parallelized? Generally the best performance
occurs when we parallelize the outermost loop of a loop nest. This way the overhead associated with
beginning a parallel loop is amortized over a longer parallel loop duration.
• Can and should the loop nest be interchanged? The compiler may detect that the loops in a nest
can be done in any order. One order may work very well for parallel code while giving poor memory
performance. Another order may give unit stride but perform poorly with multiple threads. The
compiler must analyze the cost/benet of each approach and make the best choice.
• How do we break up the iterations among the threads executing a parallel loop? Are the iterations
short with uniform duration, or long with wide variation of execution time? We will see that there
are a number of dierent ways to accomplish this. When the programmer has given no guidance, the
compiler must make an educated guess.
Even though it seems complicated, the compiler can do a surprisingly good job on a wide variety of codes.
It is not magic, however. For example, in the following code we have a loop-carried ow dependency:
PROGRAM DEP
PARAMETER(NITER=300,N=1000000)
REAL*4 A(N)
DO ITIME=1,NITER
CALL WHATEVER(A)
DO I=2,N
When we compile the code, the compiler gives us the following message:
The compiler throws its hands up in despair, and lets you know that the loop at Line 8 had an unsafe
dependence, and so it won't automatically parallelize the loop. When the code is executed below, adding a
thread does not aect the execution performance:
E6000:setenv PARALLEL 1
E6000:/bin/time dep
real 18.1
user 18.1
sys 0.0
E6000:setenv PARALLEL 2
E6000:/bin/time dep
real 18.3
user 18.2
sys 0.0
E6000:
A typical application has many loops. Not all the loops are executed in parallel. It's a good idea to run a
prole of your application, and in the routines that use most of the CPU time, check to nd out which loops
are not being parallelized. Within a loop nest, the compiler generally chooses only one loop to execute in
parallel.
• You may have a compiler ag to enable the automatic parallelization of reduction operations. Because
the order of additions can aect the nal value when computing a sum of oating-point numbers, the
compiler needs permission to parallelize summation loops.
• Flags that relax the compliance with IEEE oating-point rules may also give the compiler more ex-
ibility when trying to parallelize a loop. However, you must be sure that it's not causing accuracy
problems in other areas of your code.
• Often a compiler has a ag called unsafe optimization or assume no dependencies. While this ag
may indeed enhance the performance of an application with loops that have dependencies, it almost
certainly produces incorrect results.
There is some value in experimenting with a compiler to see the particular combination that will yield good
performance across a variety of applications. Then that set of compiler options can be used as a starting
point when you encounter a new application.
3.3.3.1 Assertions
In a previous example, we compiled a program and received the following output:
An uneducated programmer who has not read this book (or has not looked at the code) might exclaim, What
unsafe dependence, I never put one of those in my code! and quickly add a no dependencies
assertion. This
is the essence of an assertion. Instead of telling the compiler to simply parallelize the loop, the programmer
is telling the compiler that its conclusion that there is a dependence is incorrect. Usually the net result is
that the compiler does indeed parallelize the loop.
We will briey review the types of assertions that are typically supported by these compilers. An assertion
is generally added to the code using a stylized comment.
29 This content is available online at <http://cnx.org/content/m32814/1.3/>.
DO I=1,N
A(I) = A(I+K) * B(I)
ENDDO
If you know that k is greater than -1 or less than -n, you can get the compiler to parallelize the loop:
C$ASSERT NO_DEPENDENCIES
DO I=1,N
A(I) = A(I+K) * B(I)
ENDDO
Of course, blindly telling the compiler that there are no dependencies is a prescription for disaster. If k
equals -1, the example above becomes a recursive loop.
3.3.3.1.2 Relations
You will often see loops that contain some potential dependencies, making them bad candidates for a no
dependencies directive. However, you may be able to supply some local facts about certain variables. This
allows partial parallelization without compromising the results. In the code below, there are two potential
dependencies because of subscripts involving k and j:
Perhaps we know that there are no conicts with references to a[i] and a[i+k]. But maybe we aren't so
sure about c[i] and c[i+j]. Therefore, we can't say in general that there are no dependencies. However, we
may be able to say something explicit about k (like k is always greater than -1), leaving j out of it. This
information about the relationship of one expression to another is called a relation assertion. Applying a
relation assertion allows the compiler to apply its optimization to the rst statement in the loop, giving us
partial parallelization.30
Again, if you supply inaccurate testimony that leads the compiler to make unsafe optimizations, your
answer may be wrong.
30 Notice that, if you were tuning by hand, you could split this loop into two: one parallelizable and one not.
3.3.3.1.3 Permutations
As we have seen elsewhere, when elements of an array are indirectly addressed, you have to worry about
whether or not some of the subscripts may be repeated. In the code below, are the values of K(I) all unique?
Or are there duplicates?
DO I=1,N
A(K(I)) = A(K(I)) + B(I) * C
END DO
If you know there are no duplicates in K (i.e., that A(K(I)) is a permutation), you can inform the compiler
so that iterations can execute in parallel. You supply the information using a permutation assertion
.
3.3.3.1.4 No equivalences
Equivalenced arrays in FORTRAN programs provide another challenge for the compiler. If any elements of
two equivalenced arrays appear in the same loop, most compilers assume that references could point to the
same memory storage location and optimize very conservatively. This may be true even if it is abundantly
apparent to you that there is no overlap whatsoever.
You inform the compiler that references to equivalenced arrays are safe with a no equivalences
assertion.
Of course, if you don't use equivalences, this assertion has no eect.
C$ASSERT TRIPCOUNT>100
DO I=L,N
A(I) = B(I) + C(I)
END DO
Your compiler is going to look at every loop as a candidate for unrolling or parallelization. It's working in
the dark, however, because it can't tell which loops are important and tries to optimize them all. This can
lead to the surprising experience of seeing your runtime go up after optimization!
A trip count assertion provides a clue to the compiler that helps it decide how much to unroll a loop
or when to parallelize a loop.31 Loops that aren't important can be identied with low or zero trip counts.
Important loops have high trip counts.
C$ASSERT NO_SIDE_EFFECTS
DO I=1,N
CALL BIGSTUFF (A,B,C,I,J,K)
END DO
Even if the compiler has all the source code, use of common variables or equivalences may mask call inde-
pendence.
PROGRAM ONE
EXTERNAL OMP_GET_THREAD_NUM, OMP_GET_MAX_THREADS
INTEGER OMP_GET_THREAD_NUM, OMP_GET_MAX_THREADS
IGLOB = OMP_GET_MAX_THREADS()
PRINT *,'Hello There'
C$OMP PARALLEL PRIVATE(IAM), SHARED(IGLOB)
IAM = OMP_GET_THREAD_NUM()
PRINT *, 'I am ', IAM, ' of ', IGLOB
C$OMP END PARALLEL
PRINT *,'All Done'
END
The C$OMP is the sentinel that indicates that this is a directive and not just another comment. The output
of the program when run looks as follows:
% setenv OMP_NUM_THREADS 4
% a.out
Hello There
I am 0 of 4
I am 3 of 4
I am 1 of 4
I am 2 of 4
All Done
%
Execution begins with a single thread. As the program encounters the PARALLEL directive, the other threads
are activated to join the computation. So in a sense, as execution passes the rst directive, one thread
becomes four. Four threads execute the two statements between the directives. As the threads are executing
independently, the order in which the print statements are displayed is somewhat random. The threads wait
at the END PARALLEL directive until all threads have arrived. Once all threads have completed the parallel
region, a single thread continues executing the remainder of the program.
In Figure 3.22 (Data interactions during a parallel region), the PRIVATE(IAM) indicates that the IAM
variable is not shared across all the threads but instead, each thread has its own private version of the
variable. The IGLOB variable is shared across all the threads. Any modication of IGLOB appears in all the
other threads instantly, within the limitations of the cache coherency.
Figure 3.22
During the parallel region, the programmer typically divides the work among the threads. This pattern of
going from single-threaded to multithreaded execution may be repeated many times throughout the execution
of an application.
Because input and output are generally not thread-safe, to be completely correct, we should indicate that
the print statement in the parallel section is only to be executed on one processor at any one time. We use a
directive to indicate that this section of code is a critical section. A lock or other synchronization mechanism
ensures that no more than one processor is executing the statements in the critical section at any one time:
C$OMP CRITICAL
PRINT *, 'I am ', IAM, ' of ', IGLOB
C$OMP END CRITICAL
DO I=1,1000000
TMP1 = ( A(I) ** 2 ) + ( B(I) ** 2 )
TMP2 = SQRT(TMP1)
B(I) = TMP2
ENDDO
To manually parallelize this loop, we insert a directive at the beginning of the loop:
C$OMP PARALLEL DO
DO I=1,1000000
TMP1 = ( A(I) ** 2 ) + ( B(I) ** 2 )
TMP2 = SQRT(TMP1)
B(I) = TMP2
ENDDO
C$OMP END PARALLEL DO
When this statement is encountered at runtime, the single thread again summons the other threads to join
the computation. However, before the threads can start working on the loop, there are a few details that
must be handled. The PARALLEL DO directive accepts the data classication and scoping clauses as in the
parallel section directive earlier. We must indicate which variables are shared across all threads and which
variables have a separate copy in each thread. It would be a disaster to have TMP1 and TMP2 shared across
threads. As one thread takes the square root of TMP1, another thread would be resetting the contents of
TMP1. A(I) and B(I) come from outside the loop, so they must be shared. We need to augment the directive
as follows:
The iteration variable I also must be a thread-private variable. As the dierent threads increment their way
through their particular subset of the arrays, they don't want to be modifying a global value for I.
There are a number of other options as to how data will be operated on across the threads. This
summarizes some of the other data semantics available:
Firstprivate: These are thread-private variables that take an initial value from the global variable of the
same name immediately before the loop begins executing.
Lastprivate: These are thread-private variables except that the thread that executes the last iteration of
the loop copies its value back into the global variable of the same name.
Each vendor may have dierent terms to indicate these data semantics, but most support all of these common
semantics. Figure 3.23 (Variables during a parallel region) shows how the dierent types of data semantics
operate.
Now that we have the data environment set up for the loop, the only remaining problem that must be
solved is which threads will perform which iterations. It turns out that this is not a trivial task, and a wrong
choice can have a signicant negative impact on our overall performance.
C VECTOR ADD
DO IPROB=1,10000
A(IPROB) = B(IPROB) + C(IPROB)
ENDDO
C PARTICLE TRACKING
DO IPROB=1,10000
RANVAL = RAND(IPROB)
CALL ITERATE_ENERGY(RANVAL) ENDDO
ENDDO
Figure 3.23
In both loops, all the computations are independent, so if there were 10,000 processors, each processor
could execute a single iteration. In the vector-add example, each iteration would be relatively short, and the
execution time would be relatively constant from iteration to iteration. In the particle tracking example, each
iteration chooses a random number for an initial particle position and iterates to nd the minimum energy.
Each iteration takes a relatively long time to complete, and there will be a wide variation of completion
times from iteration to iteration.
These two examples are eectively the ends of a continuous spectrum of the iteration scheduling challenges
facing the FORTRAN parallel runtime environment:
Static
At the beginning of a parallel loop, each thread takes a xed continuous portion of iterations of the loop
based on the number of threads executing the loop.
Dynamic
With dynamic scheduling, each thread processes a chunk of data and when it has completed processing, a
new chunk is processed. The chunk size can be varied by the programmer, but is xed for the duration of
the loop.
These two example loops can show how these iteration scheduling approaches might operate when ex-
Figure 3.24
Since the loop body (a single statement) is short with a consistent execution time, static scheduling
should result in roughly the same amount of overall work (and time if you assume a dedicated CPU for each
thread) assigned to each thread per loop execution.
An advantage of static scheduling may occur if the entire loop is executed repeatedly. If the same iterations
are assigned to the same threads that happen to be running on the same processors, the cache might actually
contain the values for A, B, and C from the previous loop execution.33 The runtime pseudo-code for static
scheduling in the rst loop might look as follows:
It's not always a good strategy to use the static approach of giving a xed number of iterations to each
thread. If this is used in the second loop example, long and varying iteration times would result in poor load
33 The operating system and runtime library actually go to some lengths to try to make this happen. This is another reason
not to have more threads than available processors, which causes unnecessary context switching.
balancing. A better approach is to have each processor simply get the next value for IPROB each time at the
top of the loop.
That approach is called dynamic scheduling , and it can adapt to widely varying iteration times. In
Figure 3.25 (Iteration assignment in dynamic scheduling), the mapping of iterations to processors using
dynamic scheduling is shown. As soon as a processor nishes one iteration, it processes the next available
iteration in order.
Iteration assignment in dynamic scheduling
Figure 3.25
If a loop is executed repeatedly, the assignment of iterations to threads may vary due to subtle timing
issues that aect threads. The pseudo-code for the dynamic scheduled loop at runtime is as follows:
ILOCAL is used so that each thread knows which iteration is currently processing. The IPROB value is altered
by the next thread executing the critical section.
While the dynamic iteration scheduling approach works well for this particular loop, there is a signicant
negative performance impact if the programmer were to use the wrong approach for a loop. For example, if
the dynamic approach were used for the vector-add loop, the time to process the critical section to determine
which iteration to process may be larger than the time to actually process the iteration. Furthermore, any
IPROB = 1
CHUNKSIZE = 100
WHILE (IPROB <= 10000 )
BEGIN_CRITICAL_SECTION
ISTART = IPROB
IPROB = IPROB + CHUNKSIZE
END_CRITICAL_SECTION
DO ILOCAL = ISTART,ISTART+CHUNKSIZE-1
RANVAL = RAND(ILOCAL)
CALL ITERATE_ENERGY(RANVAL)
ENDDO
ENDWHILE
The choice of chunk size is a compromise between overhead and termination imbalance. Typically the
programmer must get involved through directives in order to control chunk size.
Part of the challenge of iteration distribution is to balance the cost (or existence) of the critical section
against the amount of work done per invocation of the critical section. In the ideal world, the critical section
would be free, and all scheduling would be done dynamically. Parallel/vector supercomputers with hardware
assistance for load balancing can nearly achieve the ideal using dynamic approaches with relatively small
chunk size.
Because the choice of loop iteration approach is so important, the compiler relies on directives from the
programmer to specify which approach to use. The following example shows how we can request the proper
iteration scheduling for our loops:
C VECTOR ADD
C$OMP PARALLEL DO PRIVATE(IPROB) SHARED(A,B,C) SCHEDULE(STATIC)
DO IPROB=1,10000
A(IPROB) = B(IPROB) + C(IPROB)
ENDDO
C$OMP END PARALLEL DO
C PARTICLE TRACKING
C$OMP PARALLEL DO PRIVATE(IPROB,RANVAL) SCHEDULE(DYNAMIC)
DO IPROB=1,10000
RANVAL = RAND(IPROB)
CALL ITERATE_ENERGY(RANVAL)
ENDDO
C$OMP END PARALLEL DO
3.3.5 Exercises37
Exercise 3.3.5.1
Take a static, highly parallel program with a relative large inner loop. Compile the application for
parallel execution. Execute the application increasing the threads. Examine the behavior when the
number of threads exceed the available processors. See if dierent iteration scheduling approaches
make a dierence.
Exercise 3.3.5.2
Take the following loop and execute with several dierent iteration scheduling choices. For chunk-
based scheduling, use a large chunk size, perhaps 100,000. See if any approach performs better than
static scheduling:
DO I=1,4000000
A(I) = B(I) * 2.34
ENDDO
Exercise 3.3.5.3
Execute the following loop for a range of values for N from 1 to 16 million:
DO I=1,N
34 This content is available online at <http://cnx.org/content/m32820/1.3/>.
35 On the other hand, if the person is a computer scientist, improving the performance might result in anything from a poster
session at a conference to a journal article! This makes for lots of intra-departmental masters degree projects.
36 http://cnx.org/content/m32820/latest/www.openmp.org
37 This content is available online at <http://cnx.org/content/m32819/1.3/>.
Run the loop in a single processor. Then force the loop to run in parallel. At what point do you
get better performance on multiple processors? Do the number of threads aect your observations?
Exercise 3.3.5.4
Use an explicit parallelization directive to execute the following loop in parallel with a chunk size
of 1:
J = 0
C$OMP PARALLEL DO PRIVATE(I) SHARED(J) SCHEDULE(DYNAMIC)
DO I=1,1000000
J = J + 1
ENDDO
PRINT *, J
C$OMP END PARALLEL DO
Execute the loop with a varying number of threads, including one. Also compile and execute the
code in serial. Compare the output and execution times. What do the results tell you about cache
coherency? About the cost of moving data from one cache to another, and about critical section
costs?
• FORTRAN 90
• HPF: High Performance FORTRAN
These languages are designed for use on high-end computing systems. We will follow a simple program
through each of these languages, using a simple nite-dierence computation that roughly models heat ow.
It's a classic problem that contains a great deal of parallelism and is easily solved on a wide variety of parallel
architectures.
We introduce and discuss the concept of single program multiple data (SPMD) in that we treat MIMD
computers as SIMD computers. We write our applications as if a large SIMD system were going to solve the
problem. Instead of actually using a SIMD system, the resulting application is compiled for a MIMD system.
The implicit synchronization of the SIMD systems is replaced by explicit synchronization at runtime on the
MIMD systems.
191
192 CHAPTER 4. SCALABLE PARALLEL PROCESSING
To do this we break the rod into 10 segments and track the temperature over time for each segment.
Intuitively, within a time step, the next temperature of a portion of the plate is an average of the surrounding
temperatures. Given xed temperatures at some points in the rod, the temperatures eventually converge to
a steady state after sucient time steps. Figure 4.1 (Heat ow in a rod) shows the setup at the beginning
of the simulation.
Figure 4.1
PROGRAM HEATROD
PARAMETER(MAXTIME=200)
INTEGER TICKS,I,MAXTIME
REAL*4 ROD(10)
ROD(1) = 100.0
DO I=2,9
ROD(I) = 0.0
ENDDO
ROD(10) = 0.0
DO TICKS=1,MAXTIME
IF ( MOD(TICKS,20) .EQ. 1 ) PRINT 100,TICKS,(ROD(I),I=1,10)
DO I=2,9
ROD(I) = (ROD(I-1) + ROD(I+1) ) / 2
ENDDO
ENDDO
100 FORMAT(I4,10F7.2)
END
% f77 heatrod.f
heatrod.f:
MAIN heatrod:
% a.out
1 100.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
21 100.00 87.04 74.52 62.54 51.15 40.30 29.91 19.83 9.92 0.00
41 100.00 88.74 77.51 66.32 55.19 44.10 33.05 22.02 11.01 0.00
61 100.00 88.88 77.76 66.64 55.53 44.42 33.31 22.21 11.10 0.00
81 100.00 88.89 77.78 66.66 55.55 44.44 33.33 22.22 11.11 0.00
101 100.00 88.89 77.78 66.67 55.56 44.44 33.33 22.22 11.11 0.00
121 100.00 88.89 77.78 66.67 55.56 44.44 33.33 22.22 11.11 0.00
141 100.00 88.89 77.78 66.67 55.56 44.44 33.33 22.22 11.11 0.00
161 100.00 88.89 77.78 66.67 55.56 44.44 33.33 22.22 11.11 0.00
181 100.00 88.89 77.78 66.67 55.56 44.44 33.33 22.22 11.11 0.00
%
Clearly, by Time step 101, the simulation has converged to two decimal places of accuracy as the numbers
have stopped changing. This should be the steady-state approximation of the temperature at the center of
each segment of the bar.
Now, at this point, astute readers are saying to themselves, "Um, don't look now, but that loop has a
ow dependency." You would also claim that this won't even parallelize a little bit. It is so bad you can't
even unroll the loop for a little instruction-level parallelism!
A person familiar with the theory of heat ow will also point out that the above loop doesn't exactly
implement the heat ow model. The problem is that the values on the right side of the assignment in the
ROD loop are supposed to be from the previous time step, and that the value on the left side is the next
time step. Because of the way the loop is written, the ROD(I-1) value is from the next time step, as shown
in p. 193.
This can be solved using a technique called red-black , where we alternate between two arrays. Figure 4.3
(Using two arrays to eliminate a dependency) shows how the red-black version of the computation operates.
This kills two birds with one stone! Now the mathematics is precisely correct, and
there is no recurrence.
Sounds like a real win-win situation.
Computing the new value for a cell
Figure 4.2
Figure 4.3
The only downside to this approach is that it takes twice the memory storage and twice the memory
bandwidth.3 The modied code is as follows:
PROGRAM HEATRED
PARAMETER(MAXTIME=200)
INTEGER TICKS,I,MAXTIME
REAL*4 RED(10),BLACK(10)
RED(1) = 100.0
BLACK(1) = 100.0
DO I=2,9
RED(I) = 0.0
ENDDO
RED(10) = 0.0
BLACK(10) = 0.0
DO TICKS=1,MAXTIME,2
IF ( MOD(TICKS,20) .EQ. 1 ) PRINT 100,TICKS,(RED(I),I=1,10)
3 There is another red-black approach that computes rst the even elements and then the odd elements of the rod in two
passes. This approach has no data dependencies within each pass. The ROD array never has all the values from the same
time step. Either the odd or even values are one time step ahead of the other. It ends up with a stride of two and doubles the
bandwidth but does not double the memory storage required to solve the problem.
DO I=2,9
BLACK(I) = (RED(I-1) + RED(I+1) ) / 2
ENDDO
DO I=2,9
RED(I) = (BLACK(I-1) + BLACK(I+1) ) / 2
ENDDO
ENDDO
100 FORMAT(I4,10F7.2)
END
% f77 heatred.f
heatred.f:
MAIN heatred:
% a.out
1 100.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
21 100.00 82.38 66.34 50.30 38.18 26.06 18.20 10.35 5.18 0.00
41 100.00 87.04 74.52 61.99 50.56 39.13 28.94 18.75 9.38 0.00
61 100.00 88.36 76.84 65.32 54.12 42.91 32.07 21.22 10.61 0.00
81 100.00 88.74 77.51 66.28 55.14 44.00 32.97 21.93 10.97 0.00
101 100.00 88.84 77.70 66.55 55.44 44.32 33.23 22.14 11.07 0.00
121 100.00 88.88 77.76 66.63 55.52 44.41 33.30 22.20 11.10 0.00
141 100.00 88.89 77.77 66.66 55.55 44.43 33.32 22.22 11.11 0.00
161 100.00 88.89 77.78 66.66 55.55 44.44 33.33 22.22 11.11 0.00
181 100.00 88.89 77.78 66.67 55.55 44.44 33.33 22.22 11.11 0.00
%
Interestingly, the modied program takes longer to converge than the rst version. It converges at Time step
181 rather than 101. If you look at the rst version, because of the recurrence, the heat ended up owing up
faster from left to right because the left element of each average was the next-time-step value. It may seem
nifty, but it's wrong.4 Generally, in this problem, either approach converges to the same eventual values
within the limits of oating-point representation.
This heat ow problem is extremely simple, and in its red-black form, it's inherently very parallel with
very simple data interactions. It's a good model for a wide range of problems where we are discretizing
two-dimensional or three-dimensional space and performing some simple simulations in that space.
This problem can usually be scaled up by making a ner grid. Often, the benet of scalable processors is
to allow a ner grid rather than a faster time to solution. For example, you might be able to to a worldwide
weather simulation using a 200-mile grid in four hours on one processor. Using 100 processors, you may be
able to do the simulation using a 20-mile grid in four hours with much more accurate results. Or, using 400
processors, you can do the ner grid simulation in one hour.
4 There are other algorithmic approaches to solving partial dierential equations, such as the "fast multipole method" that
accelerates convergence "legally." Don't assume that the brute force approach used here is the only method to solve this
particular problem. Programmers should always look for the best available algorithm (parallel or not) before trying to scale up
the "wrong" algorithm. For folks other than computer scientists, time to solution is more important than linear speed-up.
DO I=1,N
C(I) = A(I) + B(I)
END DO
This seems reasonable, but look what happened. We imposed an order on the calculations! Wouldn't it be
enough to say "C gets A plus B"? That would free the compiler to add the vectors using any hardware
at its disposal, using any method it likes. This is what parallel languages are about. They seek to supply
primitives suitable for expressing parallel computations.
New parallel languages aren't being proposed as rapidly as they were in the mid-1980s. Developers have
realized that you can come up with a wonderful scheme, but if it isn't compatible with FORTRAN or C, few
people will care about it. The reason is simple: there are billions of lines of C and FORTRAN code, but only
a few lines of Fizgibbet
, or whatever it is you call your new parallel language. Because of the predominance
of C and FORTRAN, the most signicant parallel language activities today seek to extend those languages,
thus protecting the 20 or 30 years of investment in programs already written.6 It is too tempting for the
developers of a new language to test their language on the eight-queens problem and the game of life, get
good results, then declare it ready for prime time and begin waiting for the hordes of programmers converting
to their particular language.
was being developed, the dominant high performance computer architectures were scalable SIMD systems
such as the Connection Machine and shared-memory vector-parallel processor systems from companies like
Cray Research.
FORTRAN 90 does a surprisingly good job of meeting the needs of these very dierent architectures.
Its features also map reasonably well onto the new shared uniform memory multiprocessors. However, as
we will see later, FORTRAN 90 alone is not yet sucient to meet the needs of the scalable distributed and
nonuniform access memory systems that are becoming dominant at the high end of computing.
The FORTRAN 90 extensions to FORTRAN 77 include:
• Array constructs
• Dynamic memory allocation and automatic variables
• Pointers
• New data types, structures
• New intrinsic functions, including many that operate on vectors or matrices
• New control structures, such as a WHERE statement
• Enhanced procedure interfaces
A = A + B
DO I=1,N
A(I) = A(I) + B(I)
ENDDO
The code generated by the compiler on your workstation may not look any dierent, but for some of the
parallel machines available now and workstations just around the corner, the dierence are signicant. The
FORTRAN 90 version states explicitly that the computations can be performed in any order, including all
in parallel at the same time.
One important eect of this is that if the FORTRAN 90 version experienced a oating-point fault adding
element 17, and you were to look at the memory in a debugger, it would be perfectly legal for element 27 to
be already computed.
You are not limited to one-dimensional arrays. For instance, the element-wise addition of two two-
dimensional arrays could be stated like this:8
8 Just in case you are wondering, A*B gives you an element-wise multiplication of array members not matrix multiplication.
That is covered by a FORTRAN 90 intrinsic function.
A = A + B
in lieu of:
DO J=1,M
DO I=1,N
A(I,J) = A(I,J) + B(I,J)
END DO
END DO
Naturally, when you want to combine two arrays in an operation, their shapes have to be compatible. Adding
a seven-element vector to an eight-element vector doesn't make sense. Neither would multiplying a 2×4 array
by a 3×4 array. When the two arrays have compatible shapes, relative to the operation being performed
upon them, we say they are in shape conformance , as in the following code:
Scalars are always considered to be in shape conformance with arrays (and other scalars). In a binary
operation with an array, a scalar is treated as an array of the same size with a single element duplicated
throughout.
Still, we are limited. When you reference a particular array, A, for example, you reference the whole
thing, from the rst element to the last. You can imagine cases where you might be interested in specifying
a subset of an array. This could be either a group of consecutive elements or something like "every eighth
element" (i.e., a non-unit stride through the array). Parts of arrays, possibly noncontiguous, are calledarray
sections .
FORTRAN 90 array sections can be specied by replacing traditional subscripts with triplets of the form
a:b:c, meaning "elements a through b, taken with an increment of c." You can omit parts of the triplet,
provided the meaning remains clear. For example, a:b means "elements a through b;" a: means "elements
from a to the upper bound with an increment of 1." Remember that a triplet replaces a single subscript, so
n n
an -dimension array can have triplets.
You can use triplets in expressions, again making sure that the parts of the expression are in conformance.
Consider these statements:
The rst statement above assigns the last 10 elements of Y to the 10th row of X. The second statement
expresses the same thing slightly dierently. The lone " : " tells the compiler that the whole range (1
through 10) is implied.
REAL A(100,10,2)
...
A = SIN(A)
Each element of array A is replaced with its sine. FORTRAN 90 intrinsics work with array sections too, as
long as the variable receiving the result is in shape conformance with the one passed:
REAL A(100,10,2)
REAL B(10,10,100)
...
B(:,:,1) = COS(A(1:100:10,:,1))
Other intrinsics, such as SQRT, LOG, etc., have been extended as well. Among the new intrinsics are:
Reductions: FORTRAN 90 has vector reductions such as MAXVAL, MINVAL, and SUM. For higher-order arrays
(anything more than a vector) these functions can perform a reduction along a particular dimension.
Additionally, there is a DOT_PRODUCT function for the vectors.
Matrix manipulation: Intrinsics MATMUL and TRANSPOSE can manipulate whole matrices.
Constructing or reshaping arrays: RESHAPE allows you to create a new array from elements of an old
one with a dierent shape. SPREAD replicates an array along a new dimension. MERGE copies portions
of one array into another under control of a mask. CSHIFT allows an array to be shifted in one or more
dimensions.
Inquiry functions: SHAPE, SIZE, LBOUND, and UBOUND let you ask questions about how an array is con-
structed.
Parallel tests: Two other new reduction intrinsics, ANY and ALL, are for testing many array elements in
parallel.
In places where the logical expression is TRUE, A gets 1.0 and C gets B+1.0. In the ELSEWHERE clause, A gets
-1.0. The result of the operation above would be arrays A and C with the elements:
Again, no order is implied in these conditional assignments, meaning they can be done in parallel. This lack
of implied order is critical to allowing SIMD computer systems and SPMD environments to have exibility
in performing these computations.
SUBROUTINE RELAX(N,A)
INTEGER N
REAL, DIMENSION (N) :: A, B
Two arrays are declared: A, the dummy argument, and B, an automatic, explicit shape array. When the
subroutine returns, B ceases to exist. Notice that the size of B is taken from one of the arguments, N.
Allocatable arrays give you the ability to choose the size of an array after examining other variables in
the program. For example, you might want to determine the amount of input data before allocating the
arrays. This little program asks the user for the matrix's size before allocating storage:
INTEGER M,N
REAL, ALLOCATABLE, DIMENSION (:,:) :: X
...
WRITE (*,*) 'ENTER THE DIMENSIONS OF X'
READ (*,*) M,N
ALLOCATE (X(M,N))
...
do something with X
...
DEALLOCATE (X)
...
The ALLOCATE statement creates an M × N array that is later freed by the DEALLOCATE statement. As with C
programs, it's important to give back allocated memory when you are done with it; otherwise, your program
might consume all the virtual storage available.
PROGRAM HEATROD
PARAMETER(MAXTIME=200)
INTEGER TICKS,I,MAXTIME
REAL*4 ROD(10)
ROD(1) = 100.0
DO I=2,9
ROD(I) = 0.0
ENDDO
ROD(10) = 0.0
DO TICKS=1,MAXTIME
IF ( MOD(TICKS,20) .EQ. 1 ) PRINT 100,TICKS,(ROD(I),I=1,10)
ROD(2:9) = (ROD(1:8) + ROD(3:10) ) / 2
ENDDO
100 FORMAT(I4,10F7.2)
END
The program is identical, except the inner loop is now replaced by a single statement that computes the
"new" section by averaging a strip of the "left" elements and a strip of the "right" elements.
The output of this program is as follows:
If you look closely, this output is the same as the red-black implementation. That is because in FORTRAN
90:
is a single
assignment statement. As shown in Figure 4.4 (Data alignment and computations), the right side
is completely evaluated before the resulting array section is assigned into ROD(2:9). For a moment, that
might seem unnatural, but consider the following statement:
I = I + 1
We know that if I starts with 5, it's incremented up to six by this statement. That happens because the
right side (5+1) is evaluated before the assignment of 6 into I is performed. In FORTRAN 90, a variable
is
can be an entire array. So, this a red-black operation. There is an "old" ROD on the right side and a "new"
ROD on the left side!
To really "think" FORTRAN 90, it's good to pretend you are on an SIMD system with millions of little
CPUs. First we carefully align the data, sliding it around, and then wham in a single instruction, we
add all the aligned values in an instant. Figure 4.4 (Data alignment and computations) shows graphically
this act of "aligning" the values and then adding them. The data ow graph is extremely simple. The top
two rows are read-only, and the data ows from top to bottom. Using the temporary space eliminates the
seeming dependency. This approach of "thinking SIMD" is one of the ways to force ourselves to focus our
thoughts on the data rather than the control. SIMD may not be a good architecture for your problem but
if you can express it so that SIMD could work, a good SPMD environment can take advantage of the data
parallelism that you have identied.
The above example actually highlights one of the challenges in producing an ecient implementation of
FORTRAN 90. If these arrays contained 10 million elements, and the compiler used a simple approach, it
would need 30 million elements for the old "left" values, the old "right" values, and for the new values. Data
ow optimization is needed to determine just how much extra data must be maintained to give the proper
results. If the compiler is clever, the extra memory can be quite small:
Figure 4.4
SAVE1 = ROD(1)
DO I=2,9
SAVE2 = ROD(I)
ROD(I) = (SAVE1 + ROD(I+1) ) / 2
SAVE1 = SAVE2
ENDDO
This does not have the parallelism that the full red-black implementation has, but it does produce the correct
results with only two extra data elements. The trick is to save the old "left" value just before you wipe it
out. A good FORTRAN 90 compiler uses data ow analysis, looking at a template of how the computation
moves across the data to see if it can save a few elements for a short period of time to alleviate the need for
a complete extra copy of the data.
The advantage of the FORTRAN 90 language is that it's up to the compiler whether it uses a complete
copy of the array or a few data elements to insure that the program executes properly. Most importantly, it
can change its approach as you move from one architecture to another.
• There is a concern that the use of pointers and dynamic data structures would ruin performance and
lose the optimization advantages of FORTRAN over C. Some people would say that FORTRAN 90
is trying to be a better C than C. Others would say, "who wants to become more like the slower
language!" Whatever the reason, there was some controversy when FORTRAN 90 was implemented,
leading to some reluctance in adoption by programmers. Some vendors said, "You can use FORTRAN
90, but FORTRAN 77 will always be faster."
• Because vendors often implemented dierent subsets of FORTRAN 90, it was not as portable as
FORTRAN 77. Because of this, users who needed maximum portability stuck with FORTRAN 77.
• Sometimes vendors purchased their fully compliant FORTRAN 90 compilers from a third party who
demanded high license fees. So, you could get the free (and faster according to the vendor) FORTRAN
77 or pay for the slower (wink wink) FORTRAN 90 compiler.
• Because of these factors, the number of serious applications developed in FORTRAN 90 was small.
So the benchmarks used to purchase new systems were almost exclusively FORTRAN 77. This fur-
ther motivated the vendors to improve their FORTRAN 77 compilers instead of their FORTRAN 90
compilers.
• As the FORTRAN 77 compilers became more sophisticated using data ow analysis, it became rela-
tively easy to write portable "parallel" code in FORTRAN 77, using the techniques we have discussed
in this book.
• One of the greatest potential benets to FORTRAN 90 was portability between SIMD and the par-
allel/vector supercomputers. As both of these architectures were replaced with the shared uniform
memory multiprocessors, FORTRAN 77 became the language that aorded the maximum portability
across the computers typically used by high performance computing programmers.
• The FORTRAN 77 compilers supported directives that allowed programmers to ne-tune the perfor-
mance of their applications by taking full control of the parallelism. Certain dialects of FORTRAN
77 essentially became parallel programming "assembly language." Even highly tuned versions of these
codes were relatively portable across the dierent vendor shared uniform memory multiprocessors.
So, events conspired against FORTRAN 90 in the short run. However, FORTRAN 77 is not well suited for
the distributed memory systems because it does not lend itself well to data layout directives. As we need to
partition and distribute the data carefully on these new systems, we must give the compiler lots
of exibility.
FORTRAN 90 is the language best suited to this purpose.
Decomposing computations: We have already discussed this technique. When the decomposition is done
based on computations, we come up with some mechanism to divide the computations (such as the
iterations of a loop) evenly among our processors. The location of the data is generally ignored, and
the primary issues are iteration duration and uniformity. This is the preferred technique for the shared
uniform memory systems because the data can be equally accessed by any processor.
Decomposing data: When memory access is nonuniform, the tendency is to focus on the distribution of
the data rather than computations. The assumption is that retrieving "remote" data is costly and
should be minimized. The data is distributed among the memories. The processor that contains the
data performs the computations on that data after retrieving any other data necessary to perform the
computation.
Decomposing tasks: - When the operations that must be performed are very independent, and take some
time, a task decomposition can be performed. In this approach a master process/thread maintains a
queue of work units. When a processor has available resources, it retrieves the next "task" from the
queue and begins processing. This is a very attractive approach for embarrassingly parallel computa-
tions.10
In some sense, the rest of this chapter is primarily about data decomposition. In a distributed memory system,
the communication costs usually are the dominant performance factor. If your problem is so embarrassingly
parallel that it can be distributed as tasks, then nearly any technique will work. Data-parallel problems
occur in many disciplines. They vary from those that are extremely parallel to those that are just sort of
parallel. For example, fractal calculations are extremely parallel; each point is derived independently of the
rest. It's simple to divide fractal calculations among processors. Because the calculations are independent,
the processors don't have to coordinate or share data.
Our heat ow problem when expressed in its red-black (or FORTRAN 90) form is extremely parallel but
requires some sharing of data. A gravitational model of a galaxy is another kind of parallel program. Each
point exerts an inuence on every other. Therefore, unlike the fractal calculations, the processors do have
to share data.
In either case, you want to arrange calculations so that processors can say to one another, "you go over
there and work on that, and I'll work on this, and we'll get together when we are nished."
Problems that oer less independence between regions are still very good candidates for domain decompo-
sition. Finite dierence problems, short-range particle interaction simulations, and columns of matrices can
be treated similarly. If you can divide the domain evenly between the processors, they each do approximately
the same amount of work on their way to a solution.
Other physical systems are not so regular or involve long-range interactions. The nodes of an unstructured
grid may not be allocated in direct correspondence to their physical locations, for instance. Or perhaps the
model involves long-range forces, such as particle attractions. These problems, though more dicult, can
be structured for parallel machines as well. Sometimes various simplications, or "lumping" of intermediate
eects, are needed. For instance, the inuence of a group of distant particles upon another may be treated
as if there were one composite particle acting at a distance. This is done to spare the communications that
would be required if every processor had to talk to every other regarding each detail. In other cases, the
parallel architecture oers opportunities to express a physical system in dierent and clever ways that make
sense in the context of the machine. For instance, each particle could be assigned to its own processor, and
these could slide past one another, summing interactions and updating a time step.
Depending on the architecture of the parallel computer and problem, a choice for either dividing or
replicating (portions of ) the domain may add unacceptable overhead or cost to the whole project.
9 This content is available online at <http://cnx.org/content/m33762/1.3/>.
10 The distributed RC5 key-cracking eort was coordinated in this fashion. Each processor would check out a block of keys
and begin testing those keys. At some point, if the processor was not fast enough or had crashed, the central system would
reissue the block to another processor. This allowed the system to recover from problems on individual computers.
• Identify scalars and arrays that will be distributed across a parallel machine.
• Say how they will be distributed. Will they be strips, blocks, or something else?
• Specify how these variables will be aligned with respect to one another.
• Redistribute and realign data structures at runtime.
• Add a FORALL control construct for parallel assignments that are dicult or impossible to construct
using FORTRAN 90's array syntax.
• Make improvements to the FORTRAN 90 WHERE control construct.
• Add intrinsic functions for common parallel operations.
There were several sources of inspiration for the HPF eort. Layout directives were already part of the
FORTRAN 90 programming environment for some SIMD computers (i.e., the CM-2). Also, PVM, the rst
portable message-passing environment, had been released a year earlier, and users had a year of experi-
ence trying to decompose by hand programs. They had developed some basic usable techniques for data
decomposition that worked very well but required far too much bookkeeping.12
The HPF eort brought together a diverse set of interests from all the major high performance computing
vendors. Vendors representing all the major architectures were represented. As a result HPF was designed
to be implemented on nearly all types of architectures.
There is an eort underway to produce the next FORTRAN standard: FORTRAN 95. FORTRAN 95 is
expected to adopt some but not all of the HPF modications.
As the user adds directives to the program, the semantics of the program are not changed. If the
user completely misunderstands the application and inserts extremely ill-conceived directives, the program
produces correct results very slowly. An HPF compiler doesn't try to "improve on" the user's directives. It
assumes the programmer is omniscient.13
Once the user has determined how the data will be distributed across the processors, the HPF compiler
attempts to use the minimum communication necessary and overlaps communication with computation
whenever possible. HPF generally uses an "owner computes" rule for the placement of the computations. A
particular element in an array is computed on the processor that stores that array element.
All the necessary data to perform the computation is gathered from remote processors, if necessary, to
perform the computation. If the programmer is clever in decomposition and alignment, much of the data
needed will be from the local memory rather then a remote memory. The HPF compiler is also responsible
for allocating any temporary data structures needed to support communications at runtime.
In general, the HPF compiler is not magic - it simply does a very good job with the communication
details when the programmer can design a good data decomposition. At the same time, it retains portability
with the single CPU and shared uniform memory systems using FORTRAN 90.
REAL*4 ROD(10)
!HPF$ DISTRIBUTE ROD(BLOCK)
The !HPF$ prex would be a comment to a non-HPF compiler and can safely be ignored by a straight
FORTRAN 90 compiler. The DISTRIBUTE directive indicates that the ROD array is to be distributed across
multiple processors. If this directive is not used, the ROD array is allocated on one processor and com-
municated to the other processors as necessary. There are several distributions that can be done in each
dimension:
REAL*4 BOB(100,100,100),RICH(100,100,100)
!HPF$ DISTRIBUTE BOB(BLOCK,CYCLIC,*)
!HPF$ DISTRIBUTE RICH(CYCLIC(10))
BLOCK The array is distributed across the processors using contiguous blocks of the index value. The blocks
are made as large as possible.
CYCLIC The array is distributed across the processors, mapping each successive element to the "next"
processor, and when the last processor is reached, allocation starts again on the rst processor.
CYCLIC(n) The array is distributed the same as CYCLIC except that n successive elements are placed on
each processor before moving on to the next processor.
Figure 4.5
Figure 4.5 (Distributing array elements to processors) shows how the elements of a simple array would
be mapped onto three processors with dierent directives.
It must allocate four elements to Processors 1 and 2 because there is no Processor 4 available for the
leftover element if it allocated three elements to Processors 1 and 2. In Figure 4.5 (Distributing array
elements to processors), the elements are allocated on successive processors, wrapping around to Processor
1 after the last processor. In Figure 4.5 (Distributing array elements to processors), using a chunk size with
CYCLIC is a compromise between pure BLOCK and pure CYCLIC.
To explore the use of the *, we can look at a simple two-dimensional array mapped onto four processors. In
Figure 4.6 (Two-dimensional distributions), we show the array layout and each cell indicates which processor
will hold the data for that cell in the two-dimensional array. In Figure 4.6 (Two-dimensional distributions),
the directive decomposes in both dimensions simultaneously. This approach results in roughly square patches
in the array. However, this may not be the best approach. In the following example, we use the * to indicate
that we want all the elements of a particular column to be allocated on the same processor. So, the column
values equally distribute the columns across the processors. Then, all the rows in each column follow where
the column has been placed. This allows unit stride for the on-processor portions of the computation and is
benecial in some applications. The * syntax is also called on-processor distribution.
Two-dimensional distributions
Figure 4.6
When dealing with more than one data structure to perform a computation, you can either separately
distribute them or use the ALIGN directive to ensure that corresponding elements of the two data structures
are to be allocated together. In the following example, we have a plate array and a scaling factor that must
be applied to each column of the plate during the computation:
DIMENSION PLATE(200,200),SCALE(200)
!HPF$ DISTRIBUTE PLATE(*,BLOCK)
!HPF$ ALIGN SCALE(I) WITH PLATE(J,I)
Or:
DIMENSION PLATE(200,200),SCALE(200)
!HPF$ DISTRIBUTE PLATE(*,BLOCK)
!HPF$ ALIGN SCALE(:) WITH PLATE(*,:)
In both examples, the PLATE and the SCALE variables are allocated to the same processors as the corresponding
columns of PLATE. The * and : syntax communicate the same information. When * is used, that dimension
is collapsed, and it doesn't participate in the distribution. When the : is used, it means that dimension
follows the corresponding dimension in the variable that has already been distributed.
You could also specify the layout of the SCALE variable and have the PLATE variable "follow" the layout
of the SCALE variable:
DIMENSION PLATE(200,200),SCALE(200)
!HPF$ DISTRIBUTE SCALE(BLOCK)
!HPF$ ALIGN PLATE(J,I) WITH SCALE(I)
You can put simple arithmetic expressions into the ALIGN directive subject to some limitations. Other
directives include:
PROCESSORS Allows you to create a shape of the processor conguration that can be used to align other
data structures.
REDISTRIBUTE and REALIGN Allow you to dynamically reshape data structures at runtime as the commu-
nication patterns change during the course of the run.
TEMPLATE Allows you to create an array that uses no space. Instead of distributing one data structure and
aligning all the other data structures, some users will create and distribute a template and then align
all of the real data structures to that template.
The use of directives can range from very simple to very complex. In some situations, you distribute the one
large shared structure, align a few related structures and you are done. In other situations, programmers
attempt to optimize communications based on the topology of the interconnection network (hypercube, multi-
stage interconnection network, mesh, or toroid) using very detailed directives. They also might carefully
redistribute the data at the various phases of the computation.
Hopefully your application will yield good performance without too much eort.
This can be expressed in native FORTRAN 90 but it is rather ugly, counterintuitive, and prone to error.
Another control structure is the ability to declare a function as "PURE." A PURE function has no
side eects other than through its parameters. The programmer is guaranteeing that a PURE function
can execute simultaneously on many processors with no ill eects. This allows HPF to assume that it will
only operate on local data and does not need any data communication during the duration of the function
execution. The programmer can also declare which parameters of the function are input parameters, output
parameters, and input-output parameters.
development, a number of these operations had been identied and implemented. HPF took the opportunity
to dene standardized syntax for these operations.
A sample of these operations includes:
While there are a large number of these intrinsic functions, most applications use only a few of the operations.
INTEGER PLATESIZ,MAXTIME
PARAMETER(PLATESIZ=2000,MAXTIME=200)
!HPF$ DISTRIBUTE PLATE(*,BLOCK)
REAL*4 PLATE(PLATESIZ,PLATESIZ)
INTEGER TICK
PLATE = 0.0
* Add Boundaries
PLATE(1,:) = 100.0
PLATE(PLATESIZ,:) = -40.0
PLATE(:,PLATESIZ) = 35.23
PLATE(:,1) = 4.5
DO TICK = 1,MAXTIME
PLATE(2:PLATESIZ-1,2:PLATESIZ-1) = (
+ PLATE(1:PLATESIZ-2,2:PLATESIZ-1) +
+ PLATE(3:PLATESIZ-0,2:PLATESIZ-1) +
+ PLATE(2:PLATESIZ-1,1:PLATESIZ-2) +
+ PLATE(2:PLATESIZ-1,3:PLATESIZ-0) ) / 4.0
PRINT 1000,TICK, PLATE(2,2)
1000 FORMAT('TICK = ',I5, F13.8)
ENDDO
*
END
% cat hostfile
frodo.egr.msu.edu
gollum.egr.msu.edu
mordor.egr.msu.edu
%
After some nontrivial machinations with paths and environment variables, you can start the PVM console:
% pvm hostfile
pvmd already running.
15 This content is available online at <http://cnx.org/content/m33781/1.3/>.
16 Notice I said not that much more eort.
17 This content is available online at <http://cnx.org/content/m33779/1.3/>.
Many dierent users can be running virtual machines using the same pool of resources. Each user has their
own view of an empty machine. The only way you might detect other virtual machines using your resources
is in the percentage of the time your applications get the CPU.
ps
There is a wide range of commands you can issue at the PVM console. The command shows the running
processes in your virtual machine. It's quite possible to have more processes than computer systems. Each
process is time-shared on a system along with all the other load on the system. The reset
command performs
a soft reboot on your virtual machine. You are the virtual system administrator of the virtual machine you
have assembled.
To execute programs on your virtual computer, you must compile and link your programs with the PVM
library routines:18
When the rst PVM call is encountered, the application contacts your virtual machine and enrolls itself in
the virtual machine. At that point it should show up in the output of the ps
command issued at the PVM
console.
From that point on, your application issues PVM calls to create more processes and interact with those
processes. PVM takes the responsibility for distributing the processes on the dierent systems in the virtual
machine, based on the load and your assessment of each system's relative performance. Messages are moved
across the network using user datagram protocol
(UDP) and delivered to the appropriate process.
Typically, the PVM application starts up some additional PVM processes. These can be additional copies
of the same program or each PVM process can run a dierent PVM application. Then the work is distributed
among the processes, and results are gathered as necessary.
18 Note: the exact compilation may be dierent on your system.
There are several basic models of computing that are typically used when working with PVM:
Master/Slave: : When operating in this mode, one process (usually the initial process) is designated as
the master that spawns some number of worker processes. Work units are sent to each worker process,
and the results are returned to the master. Often the master maintains a queue of work to be done and
as a slave nishes, the master delivers a new work item to the slave. This approach works well when
there is little data interaction and each work unit is independent. This approach has the advantage
that the overall problem is naturally load-balanced even when there is some variation in the execution
time of individual processes.
Broadcast/Gather: : This type of application is typically characterized by the fact that the shared data
structure is relatively small and can be easily copied into every processor's node. At the beginning
of the time step, all the global data structures are broadcast from the master process to all of the
processes. Each process then operates on their portion of the data. Each process produces a partial
result that is sent back and gathered by the master process. This pattern is repeated for each time
step.
SPMD/Data decomposition: : When the overall data structure is too large to have a copy stored in
every process, it must be decomposed across multiple processes. Generally, at the beginning of a time
step, all processes must exchange some data with each of their neighboring processes. Then with their
local data augmented by the necessary subset of the remote data, they perform their computations.
At the end of the time step, necessary data is again exchanged between neighboring processes, and the
process is restarted.
The most complicated applications have nonuniform data ows and data that migrates around the system
as the application changes and the load changes on the system.
In this section, we have two example programs: one is a master-slave operation, and the other is a data
decomposition-style solution to the heat ow problem.
% cat mast.c
#include <stdio.h>
#include "pvm3.h"
#define MAXPROC 5
#define JOBS 20
main()
{
int mytid,info;
int tids[MAXPROC];
int tid,input,output,answers,work;
mytid = pvm_mytid();
info=pvm_spawn("slav", (char**)0, 0, "", MAXPROC, tids);
pvm_exit();
}
%
One of the interesting aspects of the PVM interface is the separation of calls to prepare a new message, pack
data into the message, and send the message. This is done for several reasons. PVM has the capability to
convert between dierent oating-point formats, byte orderings, and character formats. This also allows a
single message to have multiple data items with dierent types.
The purpose of the message type in each PVM send or receive is to allow the sender to wait for a particular
type of message. In this example, we use two message types. Type one is a message from the master to the
slave, and type two is the response.
When performing a receive, a process can either wait for a message from a specic process or a message
from any process.
In the second phase of the computation, the master waits for a response from any slave, prints the
response, and then doles out another work unit to the slave or tells the slave to terminate by sending a
message with a value of -1.
The slave code is quite simple it waits for a message, unpacks it, checks to see if it is a termination
message, returns a response, and repeats:
% cat slav.c
#include <stdio.h>
#include "pvm3.h"
while(1) {
pvm_recv( -1, 1 ); /* -1 = any task 1=msgtype */
pvm_upkint(&input, 1, 1);
if ( input == -1 ) break; /* All done */
output = input * 2;
pvm_initsend( PvmDataDefault );
pvm_pkint( &mytid, 1, 1 );
pvm_pkint( &input, 1, 1 );
pvm_pkint( &output, 1, 1 );
pvm_send( pvm_parent(), 2 );
}
pvm_exit();
}
%
% pheat
Thanks to 262204 2*0=0
Thanks to 262205 2*1=2
Thanks to 262206 2*2=4
Thanks to 262207 2*3=6
Thanks to 262204 2*5=10
Thanks to 262205 2*6=12
Thanks to 262206 2*7=14
Thanks to 262207 2*8=16
Thanks to 262204 2*9=18
Thanks to 262205 2*10=20
Thanks to 262206 2*11=22
Thanks to 262207 2*12=24
Thanks to 262205 2*14=28
Thanks to 262207 2*16=32
Thanks to 262205 2*17=34
Thanks to 262207 2*18=36
Thanks to 262204 2*13=26
Thanks to 262205 2*19=38
Thanks to 262206 2*15=30
Thanks to 262208 2*4=8
%
Figure 4.7
The data will be spread across all of the processes using a (*, BLOCK) distribution. Columns are distributed
to processes in contiguous blocks, and all the row elements in a column are stored on the same process. As
with HPF, the process that owns a data cell performs the computations for that cell after retrieving any
data necessary to perform the computation.
We use a red-black approach but for simplicity, we copy the data back at the end of each iteration. For
a true red-black, you would perform the computation in the opposite direction every other time step.
Note that instead of spawning slave process, the parent process spawns additional copies of itself. This
is typical of SPMD-style programs. Once the additional processes have been spawned, all the processes wait
at a barrier before they look for the process numbers of the members of the group. Once the processes have
arrived at the barrier, they all retrieve a list of the dierent process numbers:
% cat pheat.f
PROGRAM PHEAT
INCLUDE '../include/fpvm3.h'
INTEGER NPROC,ROWS,COLS,TOTCOLS,OFFSET
PARAMETER(NPROC=4,MAXTIME=200)
PARAMETER(ROWS=200,TOTCOLS=200)
PARAMETER(COLS=(TOTCOLS/NPROC)+3)
REAL*8 RED(0:ROWS+1,0:COLS+1), BLACK(0:ROWS+1,0:COLS+1)
LOGICAL IAMFIRST,IAMLAST
INTEGER INUM,INFO,TIDS(0:NPROC-1),IERR
INTEGER I,R,C
INTEGER TICK,MAXTIME
CHARACTER*30 FNAME
* Find my pals and get their TIDs - TIDS are necessary for sending
DO I=0,NPROC-1
CALL PVMFGETTID('pheat', I, TIDS(I))
ENDDO
At this point in the code, we have NPROC processes executing in an SPMD mode. The next step is to determine
which subset of the array each process will compute. This is driven by the INUM variable, which ranges from
0 to 3 and uniquely identies these processes.
We decompose the data and store only one quarter of the data on each process. Using the INUM variable,
we choose our continuous set of columns to store and compute. The OFFSET variable maps between a global
column in the entire array and a local column in our local subset of the array. Figure 4.8 (Assigning grid
elements to processors) shows a map that indicates which processors store which data elements. The values
marked with a B are boundary values and won't change during the simulation. They are all set to 0. This
code is often rather tricky to gure out. Performing a (BLOCK, BLOCK) distribution requires a two-dimensional
decomposition and exchanging data with the neighbors above and below, in addition to the neighbors to the
left and right:
Figure 4.8
* Start Cold
DO C=0,COLS+1
DO R=0,ROWS+1
BLACK(R,C) = 0.0
ENDDO
ENDDO
Now we run the time steps. The rst act in each time step is to reset the heat sources. In this simulation,
we have four heat sources placed near the middle of the plate. We must restore all the values each time
through the simulation as they are modied in the main loop:
Now we perform the exchange of the ghost values with our neighboring processes. For example, Process 0
contains the elements for global column 50. To compute the next time step values for column 50, we need
column 51, which is stored in Process 1. Similarly, before Process 1 can compute the new values for column
51, it needs Process 0's values for column 50.
Figure 4.9 (Pattern of communication for ghost values) shows how the data is transferred between pro-
cessors. Each process sends its leftmost column to the left and its rightmost column to the right. Because
the rst and last processes border unchanging boundary values on the left and right respectively, this is not
necessary for columns one and 200. If all is done properly, each process can receive its ghost values from
their left and right neighbors.
Figure 4.9
The net result of all of the transfers is that for each space that must be computed, it's surrounded by
one layer of either boundary values or ghost values from the right or left neighbors:
ENDIF
This next segment is the easy part. All the appropriate ghost values are in place, so we must simply perform
the computation in our subspace. At the end, we copy back from the RED to the BLACK array; in a real
simulation, we would perform two time steps, one from BLACK to RED and the other from RED to BLACK, to
save this extra copy:
* Copy back - Normally we would do a red and black version of the loop
DO C=1,MYLEN
DO R=1,ROWS
BLACK(R,C) = RED(R,C)
ENDDO
ENDDO
ENDDO
Now we nd the center cell and send to the master process (if necessary) so it can be printed out. We also
dump out the data into les for debugging or later visualization of the results. Each le is made unique by
appending the instance number to the lename. Then the program terminates:
CALL SENDCELL(RED,ROWS,COLS,OFFSET,MYLEN,INUM,TIDS(0),
+ ROWS/2,TOTCOLS/2)
The SENDCELL routine nds a particular cell and prints it out on the master process. This routine is called
in an SPMD style: all the processes enter this routine although all not at precisely the same time. Depending
on the INUM and the cell that we are looking for, each process may do something dierent.
If the cell in question is in the master process, and we are the master process, print it out. All other
processes do nothing. If the cell in question is stored in another process, the process with the cell sends it
to the master processes. The master process receives the value and prints it out. All the other processes do
nothing.
This is a simple example of the typical style of SPMD code. All the processes execute the code at roughly
the same time, but, based on information local to each process, the actions performed by dierent processes
may be quite dierent:
SUBROUTINE SENDCELL(RED,ROWS,COLS,OFFSET,MYLEN,INUM,PTID,R,C)
INCLUDE '../include/fpvm3.h'
INTEGER ROWS,COLS,OFFSET,MYLEN,INUM,PTID,R,C
REAL*8 RED(0:ROWS+1,0:COLS+1)
REAL*8 CENTER
Like the previous routine, the STORE routine is executed on all processes. The idea is to store a value into
a global row and column position. First, we must determine if the cell is even in our process. If the cell is
in our process, we must compute the local column (I) in our subset of the overall matrix and then store the
value:
SUBROUTINE STORE(RED,ROWS,COLS,OFFSET,MYLEN,R,C,VALUE,INUM)
REAL*8 RED(0:ROWS+1,0:COLS+1)
REAL VALUE
INTEGER ROWS,COLS,OFFSET,MYLEN,R,C,I,INUM
I = C - OFFSET
IF ( I .LT. 1 .OR. I .GT. MYLEN ) RETURN
RED(R,I) = VALUE
RETURN
END
% pheat
INUM: 0 Local 1 50 Global 1 50
Master Received 100 100 3.4722390023541D-07
%
We see two lines of print. The rst line indicates the values that Process 0 used in its geometry computation.
The second line is the output from the master process of the temperature at cell (100,100) after 200 time
steps.
One interesting technique that is useful for debugging this type of program is to change the number of
processes that are created. If the program is not quite moving its data properly, you usually get dierent
results when dierent numbers of processes are used. If you look closely, the above code performs correctly
with one process or 30 processes.
Notice that there is no barrier operation at the end of each time step. This is in contrast to the way
parallel loops operate on shared uniform memory multiprocessors that force a barrier at the end of each loop.
Because we have used an owner computes rule, and nothing is computed until all the required ghost data
is received, there is no need for a barrier. The receipt of the messages with the proper ghost values allows a
process to begin computing immediately without regard to what the other processes are currently doing.
This example can be used either as a framework for developing other grid-based computations, or as a
good excuse to use HPF and appreciate the hard work that the HPF compiler developers have done. A well-
done HPF implementation of this simulation should outperform the PVM implementation because HPF can
make tighter optimizations. Unlike us, the HPF compiler doesn't have to keep its generated code readable.
Communicators: : A communicator is a subset of the active processes that can be treated as a group
for collective operations such as broadcast, reduction, barriers, sending, or receiving. Within each
communicator, a process has a rank
that ranges from zero to the size of the group. A process may be a
member of more than one communicator and have a dierent rank within each communicator. There
is a default communicator that refers to all the MPI processes that is called MPI_COMM_WORLD.
Topologies: : A communicator can have a topology associated with it. This arranges the processes that
belong to a communicator into some layout. The most common layout is a Cartesian decomposition.
For example, 12 processes may be arranged into a 3×4 grid.22 Once these topologies are dened, they
can be queried to nd the neighboring processes in the topology. In addition to the Cartesian (grid)
topology, MPI also supports a graph-based topology.
Communication modes: : MPI supports multiple styles of communication, including blocking and non-
blocking. Users can also choose to use explicit buers for sending or allow MPI to manage the buers.
The nonblocking capabilities allow the overlap of communication and computation. MPI can support
a model in which there is no available memory space for buers and the data must be copied directly
from the address space of the sending process to the memory space of the receiving process. MPI
also supports a single call to perform a send and receive that is quite useful when processes need to
exchange data.
Single-call collective operations: : Some of the calls in MPI automate collective operations in a single
call. For example, the broadcast operation sends values from the master to the slaves and receives
the values on the slaves in the same operation. The net result is that the values are updated on all
20 This content is available online at <http://cnx.org/content/m33783/1.3/>.
21 One should not diminish the positive contributions of PVM, however. PVM was the rst widely avail- able portable message-
passing environment. PVM pioneered the idea of heterogeneous distributed computing with built-in format conversion.
22 Sounds a little like HPF, no?
processes. Similarly, there is a single call to sum a value across all of the processes to a single value.
By bundling all this functionality into a single call, systems that have support for collective operations
in hardware can make best use of this hardware. Also, when MPI is operating on a shared-memory
environment, the broadcast can be simplied as all the slaves simply make a local copy of a shared
variable.
Clearly, the developers of the MPI specication had signicant experience with developing message-passing
applications and added many widely used features to the message-passing library. Without these features,
each programmer needed to use more primitive operations to construct their own versions of the higher-level
operations.
PROGRAM MHEATC
INCLUDE 'mpif.h'
INCLUDE 'mpef.h'
INTEGER ROWS,COLS,TOTCOLS
PARAMETER(MAXTIME=200)
* This simulation can be run on MINPROC or greater processes.
* It is OK to set MINPROC to 1 for testing purposes
* For a large number of rows and columns, it is best to set MINPROC
* to the actual number of runtime processes
PARAMETER(MINPROC=2) PARAMETER(ROWS=200,TOTCOLS=200,COLS=TOTCOLS/MINPROC)
DOUBLE PRECISION RED(0:ROWS+1,0:COLS+1),BLACK(0:ROWS+1,0:COLS+1)
INTEGER S,E,MYLEN,R,C
INTEGER TICK,MAXTIME
CHARACTER*30 FNAME
The basic data structures are much the same as in the PVM example. We allocate a subset of the heat arrays
in each process. In this example, the amount of space allocated in each process is set by the compile-time
variable MINPROC. The simulation can execute on more than MINPROC processes (wasting some space in each
process), but it can't execute on less than MINPROC processes, or there won't be sucient total space across
all of the processes to hold the array:
INTEGER COMM1D,INUM,NPROC,IERR
INTEGER DIMS(1),COORDS(1)
LOGICAL PERIODS(1)
LOGICAL REORDER
INTEGER NDIM
INTEGER STATUS(MPI_STATUS_SIZE)
INTEGER RIGHTPROC, LEFTPROC
The call to MPI_INIT creates the appropriate number of processes. Note that in the output, the PRINT
statement before the call only appears once, but the second PRINT appears once for each process. We call
MPI_COMM_SIZE to determine the size of the global communicator MPI_COMM_WORLD. We use this value to set
up our Cartesian topology:
DIMS(1) = NPROC
PERIODS(1) = .FALSE.
REORDER = .TRUE.
NDIM = 1
Now we create a one-dimensional (NDIM=1) arrangement of all of our processes (MPI_COMM_WORLD). All of the
parameters on this call are input values except for COMM1D and IERR. COMM1D is an integer communicator
handle. If you print it out, it will be a value such as 134. It is not actually data, it is merely a handle that is
used in other calls. It is quite similar to a le descriptor or unit number used when performing input-output
to and from les.
The topology we use is a one-dimensional decomposition that isn't periodic. If we specied that we wanted
a periodic decomposition, the far-left and far-right processes would be neighbors in a wrapped-around fashion
making a ring. Given that it isn't periodic, the far-left and far-right processes have no neighbors.
In our PVM example above, we declared that Process 0 was the far-right process, Process NPROC-1 was
the far-left process, and the other processes were arranged linearly between those two. If we set REORDER
to .FALSE., MPI also chooses this arrangement. However, if we set REORDER to .TRUE., MPI may choose to
arrange the processes in some other fashion to achieve better performance, assuming that you are commu-
nicating with close neighbors.
Once the communicator is set up, we use it in all of our communication operations:
Within each communicator, each process has a rank from zero to the size of the communicator minus 1. The
MPI_COMM_RANK tells each process its rank within the communicator. A process may have a dierent rank in
the COMM1D communicator than in the MPI_COMM_WORLD communicator because of some reordering.
Given a Cartesian topology communicator,23 we can extract information from the communicator using
the MPI_CART_GET routine:
In this call, all of the parameters are output values rather than input values as in the MPI_CART_CREATE
call. The COORDS variable tells us our coordinates within the communicator. This is not so useful in our
one-dimensional example, but in a two-dimensional process decomposition, it would tell our current position
in that two-dimensional grid:
* Returns the left and right neighbors 1 unit away in the zeroth dimension
* of our Cartesian map - since we are not periodic, our neighbors may
* not always exist - MPI_CART_SHIFT handles this for us
We can use MPI_CART_SHIFT to determine the rank number of our left and right neighbors, so we can exchange
our common points with these neighbors. This is necessary because we can't simply send to INUM-1 and
INUM+1 if MPI has chosen to reorder our Cartesian decomposition. If we are the far-left or far-right process,
the neighbor that doesn't exist is set to MPI_PROC_NULL, which indicates that we have no neighbor. Later
when we are performing message sending, it checks this value and sends messages only to real processes. By
not sending the message to the null process, MPI has saved us an IF test.
To determine which strip of the global array we store and compute in this process, we call a utility
routine called MPE_DECOMP1D that simply does several calculations to evenly split our 200 columns among
our processes in contiguous strips. In the PVM version, we need to perform this computation by hand.
The MPE_DECOMP1D routine is an example of an extended MPI library call (hence the MPE prex). These
extensions include graphics support and logging tools in addition to some general utilities. The MPE library
23 Remember, each communicator may have a topology associated with it. A topology can be grid, graph, or none. Interest-
ingly, the MPI_COMM_WORLD communicator has no topology associated with it.
* Start Cold
DO C=0,COLS+1
DO R=0,ROWS+1
BLACK(R,C) = 0.0
ENDDO
ENDDO
As in the PVM example, we set the plate (including boundary values) to zero.
All processes begin the time step loop. Interestingly, like in PVM, there is no need for any synchronization.
The messages implicitly synchronize our loops.
The rst step is to store the permanent heat sources. We need to use a routine because we must make
the store operations relative to our strip of the global array:
All of the processes set these values independently depending on which process has which strip of the overall
array.
Now we exchange the data with our neighbors as determined by the Cartesian communicator. Note that
we don't need an IF test to determine if we are the far-left or far-right process. If we are at the edge, our
neighbor setting is MPI_PROC_NULL and the MPI_SEND and MPI_RECV calls do nothing when given this as a
source or destination value, thus saving us an IF test.
Note that we specify the communicator COMM1D because the rank values we are using in these calls are
relative to that communicator:
Just to show o, we use both the separate send and receive, and the combined send and receive. When
given a choice, it's probably a good idea to use the combined operations to give the runtime environment
more exibility in terms of buering. One downside to this that occurs on a network of workstations (or
any other high-latency interconnect) is that you can't do both send operations rst and then do both receive
operations to overlap some of the communication delay.
Once we have all of our ghost points from our neighbors, we can perform the algorithm on our subset of
the space:
* Copy back - Normally we would do a red and black version of the loop
DO C=1,MYLEN
DO R=1,ROWS
BLACK(R,C) = RED(R,C)
ENDDO
ENDDO
ENDDO
Again, for simplicity, we don't do the complete red-black computation.24 We have no synchronization at the
bottom of the loop because the messages implicitly synchronize the processes at the top of the next loop.
Again, we dump out the data for verication. As in the PVM example, one good test of basic correctness
is to make sure you get exactly the same results for varying numbers of processes:
As in the PVM example, we need a routine to store a value into the proper strip of the global array. This
routine simply checks to see if a particular global element is in this process and if so, computes the proper
location within its strip for the value. If the global element is not in this process, this routine simply returns
doing nothing:
SUBROUTINE STORE(RED,ROWS,COLS,S,E,R,C,VALUE,INUM)
REAL*8 RED(0:ROWS+1,0:COLS+1)
REAL VALUE
INTEGER ROWS,COLS,S,E,R,C,I,INUM
IF ( C .LT. S .OR. C .GT. E ) RETURN
I = ( C - S ) + 1
* PRINT *,'STORE, INUM,R,C,S,E,R,I',INUM,R,C,S,E,R,I,VALUE RED(R,I) = VALUE
RETURN
END
As you can see, we call MPI_INIT to activate the four processes. The PRINT statement immediately after
the MPI_INIT call appears four times, once for each of the activated processes. Then each process prints
out the strip of the array it will process. We can also see the neighbors of each process including -1 when a
process has no neighbor to the left or right. Notice that Process 0 has no left neighbor, and Process 3 has
no right neighbor. MPI has provided us the utilities to simplify message-passing code that we need to add
to implement this type of grid- based application.
When you compare this example with a PVM implementation of the same problem, you can see some of
the contrasts between the two approaches. Programmers who wrote the same six lines of code over and over
in PVM combined them into a single call in MPI. In MPI, you can think data parallel and express your
program in a more data-parallel fashion.
In some ways, MPI feels less like assembly language than PVM. However, MPI does take a little getting
used to when compared to PVM. The concept of a Cartesian communicator may seem foreign at rst, but
with understanding, it becomes a exible and powerful tool.
PROGRAM MHEAT
INCLUDE 'mpif.h'
INCLUDE 'mpef.h'
INTEGER ROWS,COLS
PARAMETER(MAXTIME=200)
PARAMETER(ROWS=200,COLS=200)
DOUBLE PRECISION RED(0:ROWS+1,0:COLS+1),BLACK(0:ROWS+1,0:COLS+1)
We need fewer variables for the MPI calls because we aren't creating a communicator. We simply use the
default communicator MPI_COMM_WORLD. We start up our processes, and nd the size and rank of our process
group:
INTEGER INUM,NPROC,IERR,SRC,DEST,TAG
Since we are broadcasting initial values to all of the processes, we only have to set things up on the master
process:
* Start Cold
IF ( INUM.EQ.0 ) THEN
DO C=0,COLS+1
DO R=0,ROWS+1
BLACK(R,C) = 0.0
ENDDO
ENDDO
ENDIF
As we run the time steps (again with no synchronization), we set the persistent heat sources directly. Since
the shape of the data structure is the same in the master and all other processes, we can use the real array
coordinates rather than mapping them as with the previous examples. We could skip the persistent settings
on the nonmaster processes, but it doesn't hurt to do it on all processes:
Now we broadcast the entire array from process rank zero to all of the other processes in the MPI_COMM_WORLD
communicator. Note that this call does the sending on rank zero process and receiving on the other processes.
The net result of this call is that all the processes have the values formerly in the master process in a single
call:
Now we perform the subset computation on each process. Note that we are using global coordinates because
the array has the same shape on each of the processes. All we need to do is make sure we set up our particular
strip of columns according to S and E:
DO C=S,E
DO R=1,ROWS
RED(R,C) = ( BLACK(R,C) +
+ BLACK(R,C-1) + BLACK(R-1,C) +
+ BLACK(R+1,C) + BLACK(R,C+1) ) / 5.0
ENDDO
ENDDO
Now we need to gather the appropriate strips from the processes into the appropriate strip in the master
array for rebroadcast in the next time step. We could change the loop in the master to receive the messages
in any order and check the STATUS variable to see which strip it received:
We use MPE_DECOMP1D to determine which strip we're receiving from each process.
In some applications, the value that must be gathered is a sum or another single value. To accomplish
this, you can use one of the MPI reduction routines that coalesce a set of distributed values into a single
value using a single call.
Again at the end, we dump out the data for testing. However, since it has all been gathered back onto
the master process, we only need to dump it on one process:
CALL MPI_FINALIZE(IERR)
END
When this program executes with four processes, it produces the following output:
% mpif77 -c mheat.f
mheat.f:
MAIN mheat:
% mpif77 -o mheat mheat.o -lmpe
% mheat -np 4
Calling MPI_INIT
My Share 1 4 51 100
My Share 0 4 1 50
My Share 3 4 151 200
My Share 2 4 101 150
%
The ranks of the processes and the subsets of the computations for each process are shown in the output.
So that is a somewhat contrived example of the broadcast/gather approach to parallelizing an applica-
tion. If the data structures are the right size and the amount of computation relative to communication is
appropriate, this can be a very eective approach that may require the smallest number of code modications
compared to a single-processor version of the code.
25 http://www.netlib.org/mpi/
26 This content is available online at <http://cnx.org/content/m33784/1.3/>.
239
240 APPENDIX
you have a problem, you can debug it at the source level. However, 30 years ago, respectable programmers
understood the machine's instruction set. High-level language compilers were commonly available, but
they didn't generate the fastest code, and they weren't terribly thrifty with memory. When programming,
you needed to save both space and time, which meant you knew how to program in assembly language.
Accordingly, you could develop an opinion about the machine's instruction set. A good instruction set
was both easy to use and powerful. In many ways these qualities were the same: powerful instructions
accomplished a lot, and saved the programmer from specifying many little steps which, in turn, made
them easy to use. But they had other, less apparent (though perhaps more important) features as well:
powerful instructions saved memory and time.
Back then, computers had very little storage by today's standards. An instruction that could roll all the
steps of a complex operation, such as a do-loop, into single opcode4 was a plus, because memory was precious.
To put some stakes in the ground, consider the last vacuum-tube computer that IBM built, the model 704
(1956). It had hardware oating-point, including a division operation, index registers, and instructions that
could operate directly on memory locations. For instance, you could add two numbers together and store the
result back into memory with a single command. The Philco 2000, an early transistorized machine (1959),
had an operation that could repeat a sequence of instructions until the contents of a counter was decremented
to zero very much like a do-loop. These were complex operations, even by today's standards. However,
both machines had a limited amount of memory 32-K words. The less memory your program took up,
the more you had available for data, and the less likely that you would have to resort to overlaying portions
of the program on top of one another.
Complex instructions saved time, too. Almost every large computer following the IBM 704 had a memory
system that was slower than its central processing unit (CPU). When a single instruction can perform several
operations, the overall number of instructions retrieved from memory can be reduced. Minimizing the number
of instructions was particularly important because, with few exceptions, the machines of the late 1950s were
very sequential; not until the current instruction was completed did the computer initiate the process of
going out to memory to get the next instruction.5 By contrast, modern machines form something of a
bucket brigade passing instructions in from memory and guring out what they do on the way so there
are fewer gaps in processing.
If the designers of early machines had had very fast and abundant instruction memory, sophisticated
compilers, and the wherewithal to build the instruction bucket brigade cheaply they might have
chosen to create machines with simple instruction sets. At the time, however, technological choices indicated
that instructions should be powerful and thrifty with memory.
As it turned out, assembly language programmers used the complicated machine instructions, but com-
pilers generally did not. It was dicult enough to get a compiler to recognize when a complicated instruction
could be used, but the real problem was one of optimizations: verbatim translation of source constructs isn't
very ecient. An optimizing compiler works by simplifying and eliminating redundant computations. After
a pass through an optimizing compiler, opportunities to use the complicated instructions tend to disappear.
• The number of transistors that could t on a single chip was increasing. It was clear that one would
eventually be able to t all the components from a processor board onto a single chip.
• Techniques such as pipelining were being explored to improve performance. Variable-length instructions
and variable-length instruction execution times (due to varying numbers of microcode steps) made
implementing pipelines more dicult.
• As compilers improved, they found that well-optimized sequences of stream- lined instructions often
outperformed the equivalent complicated multi-cycle instructions. (See Appendix A, Processor Archi-
tectures, and Appendix B, Looking at Assembly Language.)
The RISC designers sought to create a high performance single-chip processor with a fast clock rate. When
a CPU can t on a single chip, its cost is decreased, its reliability is increased, and its clock speed can be
increased. While not all RISC processors are single-chip implementation, most use a single chip.
To accomplish this task, it was necessary to discard the existing CISC instruction sets and develop a new
minimal instruction set that could t on a single chip. Hence the term reduced instruction set computer
. In
a sense reducing the instruction set was not an end but a means to an end.
For the rst generation of RISC chips, the restrictions on the number of components that could be
manufactured on a single chip were severe, forcing the designers to leave out hardware support for some
instructions. The earliest RISC processors had no oating-point support in hardware, and some did not
even support integer multiply in hardware. However, these instructions could be implemented using software
routines that combined other instructions (a microcode of sorts).
These earliest RISC processors (most severely reduced) were not overwhelming successes for four reasons:
• It took time for compilers, operating systems, and user software to be retuned to take advantage of
the new processors.
• If an application depended on the performance of one of the software-implemented instructions, its
performance suered dramatically.
• Because RISC instructions were simpler, more instructions were needed to accomplish the task.
• Because all the RISC instructions were 32 bits long, and commonly used CISC instructions were as
short as 8 bits, RISC program executables were often larger.
As a result of these last two issues, a RISC program may have to fetch more memory for its instructions
than a CISC program. This increased appetite for instructions actually clogged the memory bottleneck until
sucient caches were added to the RISC processors. In some sense, you could view the caches on RISC
processors as the microcode store in a CISC processor. Both reduced the overall appetite for instructions
that were loaded from memory.
6 This content is available online at <http://cnx.org/content/m33673/1.3/>.
• Instruction pipelining
• Pipelining oating-point execution
• Uniform instruction length
• Delayed branching
• Load/store architecture
• Simple addressing modes
This list highlights the dierences between RISC and CISC processors. Naturally, the two types of
instruction-set architectures have much in common; each uses registers, memory, etc. And many of these
techniques are used in CISC machines too, such as caches and instruction pipelines. It is the fundamental
dierences that give RISC its speed advantage: focusing on a smaller set of less powerful instructions makes
it possible to build a faster computer.
However, the notion that RISC machines are generally simpler than CISC machines isn't correct. Other
features, such as functional pipelines, sophisticated memory systems, and the ability to issue two or more
instructions per clock make the latest RISC processors the most complicated ever built. Furthermore, much
of the complexity that has been lifted from the instruction set has been driven into the compilers, making a
good optimizing compiler a prerequisite for machine performance.
Let's put ourselves in the role of computer architect again and look at each item in the list above to
understand why it's important.
7 And they did it without ever taking out a single instruction!
8 The typical CISC microprocessor in the 1980s supported oating-point operations in a separate coprocessor.
5.1.3.2 Pipelines
Everything within a digital computer (RISC or CISC) happens in step with a clock
: a signal that paces the
computer's circuitry. The rate of the clock, or clock speed
, determines the overall speed of the processor.
There is an upper limit to how fast you can clock a given computer.
A number of parameters place an upper limit on the clock speed, including the semiconductor technology,
packaging, the length of wires tying the pieces together, and the longest path in the processor. Although
it may be possible to reach blazing speed by optimizing all of the parameters, the cost can be prohibitive.
Furthermore, exotic computers don't make good oce mates; they can require too much power, produce too
much noise and heat, or be too large. There is incentive for manufacturers to stick with manufacturable and
marketable technologies.
Reducing the number of clock ticks it takes to execute an individual instruction is a good idea, though cost
and practicality become issues beyond a certain point. A greater benet comes from partially overlapping
instructions so that more than one can be in progress simultaneously. For instance, if you have two additions
to perform, it would be nice to execute them both at the same time. How do you do that? The rst,
and perhaps most obvious, approach, would be to start them simultaneously. Two additions would execute
together and complete together in the amount of time it takes to perform one. As a result, the throughput
would be eectively doubled. The downside is that you would need hardware for two adders in a situation
where space is usually at a premium (especially for the early RISC processors).
Other approaches for overlapping execution are more cost-eective than side-by-side execution. Imagine
what it would be like if, a moment after launching one operation, you could launch another without waiting
for the rst to complete. Perhaps you could start another of the same type right behind the rst one like
the two additions. This would give you nearly the performance of side-by-side execution without duplicated
hardware. Such a mechanism does exist to varying degrees in all computers CISC and RISC. It's called
a pipeline. A pipeline takes advantage of the fact that many operations are divided into identiable steps,
each of which uses dierent resources on the processor.9
A Pipeline
Figure 5.1
Figure 5.1 (A Pipeline) shows a conceptual diagram of a pipeline. An operation entering at the left
proceeds on its own for ve clock ticks before emerging at the right. Given that the pipeline stages are
independent of one another, up to ve operations can be in ight at a time as long as each instruction is
9 Here is a simple analogy: imagine a line at a fast-food drive up window. If there is only one window, one customer orders
and pays, and the food is bagged and delivered to the customer before the second customer orders. For busier restaurants,
there are three windows. First you order, then move ahead. Then at a second window, you pay and move ahead. At the third
window you pull up, grab the food and roar o into the distance. While your wait at the three-window (pipelined) drive-up
may have been slightly longer than your wait at the one-window (non-pipelined) restaurant, the pipeline solution is signicantly
better because multiple customers are being processed simultaneously.
Ideally, instruction
Our pipeline is ve stages deep, so it should be possible to get ve instructions in ight all at once. If
we could keep it up, we would see one instruction complete per clock cycle.
Simple as this illustration seems, instruction pipelining is complicated in real life. Each step must be able
to occur on dierent instructions simultaneously, and delays in any stage have to be coordinated with all
those that follow. In Figure 5.2 (Three instructions in ight through one pipeline) we see three instructions
being executed simultaneously by the processor, with each instruction in a dierent stage of execution.
Figure 5.2
For instance, if a complicated memory access occurs in stage three, the instruction needs to be delayed
before going on to stage four because it takes some time to calculate the operand's address and retrieve it
from memory. All the while, the rest of the pipeline is stalled. A simpler instruction, sitting in one of the
earlier stages, can't continue until the trac ahead clears up.
Now imagine how a jump to a new program address, perhaps caused by an if statement, could disrupt
the pipeline ow. The processor doesn't know an instruction is a branch until the decode stage. It usually
doesn't know whether a branch will be taken or not until the execute stage. As shown in Figure 5.3
(Detecting a branch), during the four cycles after the branch instruction was fetched, the processor blindly
fetches instructions sequentially and starts these instructions through the pipeline.
Figure 5.3
If the branch falls through, then everything is in great shape; the pipeline simply executes the next
instruction. It's as if the branch were a no-op instruction. However, if the branch jumps away, those three
partially processed instructions never get executed. The rst order of business is to discard these in-ight
instructions from the pipeline. It turns out that because none of these instructions was actually going to
do anything until its execute stage, we can throw them away without hurting anything (other than our
eciency). Somehow the processor has to be able to clear out the pipeline and restart the pipeline at the
branch destination.
Unfortunately, branch instructions occur every ve to ten instructions in many programs. If we executed
a branch every fth instruction and only half our branches fell through, the lost eciency due to restarting
the pipeline after the branches would be 20 percent.
You need optimal conditions to keep the pipeline moving. Even in less-than-optimal conditions, instruc-
tion pipelining is a big win especially for RISC processors. Interestingly, the idea dates back to the
late 1950s and early 1960s with the UNI- VAC LARC and the IBM Stretch. Instruction pipelining became
mainstreamed in 1964, when the CDC 6600 and the IBM S/360 families were introduced with pipelined
instruction units on machines that represented RISC-ish and CISC designs, respectively. To this day,
ever more sophisticated techniques are being applied to instruction pipelining, as machines that can overlap
instruction execution become commonplace.
several stages without delaying the rest of the processor. The result appears in a register at some point in
the future.
Some processors are limited in the amount of overlap their oating-point pipelines can support. Internal
components of the pipelines may be shared (for adding, multiplying, normalizing, and rounding intermediate
results), forcing restrictions on when and how often you can begin new operations. In other cases, oating-
point operations can be started every cycle regardless of the previous oating- point operations. We say that
such operations are fully pipelined .
The number of stages in oating-point pipelines for aordable computers has decreased over the last
10 years. More transistors and newer algorithms make it possible to perform a oating-point addition or
multiplication in just one to three cycles. Generally the most dicult instruction to perform in a single cycle
is the oating-point multiply. However, if you dedicate enough hardware to it, there are designs that can
operate in a single cycle at a moderate clock rate.
Figure 5.4
The processor has no way of knowing how long an instruction will be until it reaches the decode stage
and determines what it is. If it turns out to be a long instruction, the processor may have to go back to
memory and get the portion left behind; this stalls the pipeline. We could eliminate the problem by requiring
that all instructions be the same length, and that there be a limited number of instruction formats as shown
in Figure 5.5 (Variable-length CISC versus xed-length RISC instructions). This way, every instruction
entering the pipeline is known a priori to be complete not needing another memory access. It would
also be easier for the processor to locate the instruction elds that specify registers or constants. Altogether
because RISC can assume a xed instruction length, the pipeline ows much more smoothly.
Figure 5.5
branch and consult a table that kept the recent behavior of the branch; it would then make a guess. Based
on the guess, the CPU would immediately begin fetching at the predicted location. As long as the guesses
were correct, branches cost exactly the same as any other instruction.
If the prediction was wrong, the instructions that were in process had to be can- celled, resulting in
wasted time and eort. A simple branch prediction scheme is typically correct well over 90% of the time,
signicantly reducing the overall negative performance impact of pipeline stalls due to branches. All recent
RISC designs incorporate some type of branch prediction, making branch delay slots eectively unnecessary.
Another mechanism for reducing branch penalties is conditional execution. These are instructions that
look like branches in source code, but turn out to be a special type of instruction in the object code. They are
very useful because they replace test and branch sequences altogether. The following lines of code capture
the sense of a conditional branch:
IF ( B < C ) THEN
A = D
ELSE
A = E
ENDIF
Using branches, this would require at least two branches to ensure that the proper value ended up in A.
Using conditional execution, one might generate code that looks as follows:
COMPARE B < C
IF TRUE A = D conditional instruction
IF FALSE A = E conditional instruction
Explicit load and store instructions can kick o memory references in the pipeline's execute stage, to be
completed at a later time (they might complete immediately; it depends on the processor and the cache).
An operation downstream may require the result of the reference, but that's all right, as long as it is far
enough downstream that the reference has had time to complete.
• Improve the manufacturing processes to simply make the clock rate faster. Take a simple design; make
it smaller and faster. This approach was taken by the Alpha processors from DEC. Alpha processors
typically have had clock rates double those of the closest competitor.
• Add duplicate compute elements on the space available as we can manufacture chips with more transis-
tors. This could allow two instructions to be executed per cycle and could double performance without
increasing clock rate. This technique is called superscalar.
• Increase the number of stages in the pipeline above ve. If the instructions can truly be decom-
posed evenly into, say, ten stages, the clock rate could theoretically be doubled without requiring new
manufacturing processes. This technique was called superpipelining
. The MIPS processors used this
technique with some success.
or vice versa). You could also execute multiple xed-point instructions compares, integer additions, etc.
at the same time, provided that they, too, are independent. Another term used to describe superscalar
processors ismultiple instruction issueprocessors.
Figure 5.6
The number and variety of operations that can be run in parallel depends on both the program and
the processor. The program has to have enough usable parallelism so that there are multiple things to do,
and the processor has to have an appropriate assortment of functional units and the ability to keep them
busy. The idea is conceptually simple, but it can be a challenge for both hardware designers and compiler
writers. Every opportunity to do several things in parallel exposes the danger of violating some precedence
(i.e., performing computations in the wrong order).
Figure 5.7
Theoretically, if the reduced complexity allows the processor to clock faster, you can achieve nearly
the same performance as superscalar processors, yet without instruction mix preferences. For illustration,
picture a superscalar processor with two units xed- and oating-point executing a program that is
composed solely of xed-point calculations; the oating-point unit goes unused. This reduces the superscalar
performance by one half compared to its theoretical maximum. A superpipelined processor, on the other
hand, will be perfectly happy to handle an unbalanced instruction mix at full speed.
Superpipelines are not new; deep pipelines have been employed in the past, notably on the CDC 6600.
The label is a marketing creation to draw contrast to superscalar processing, and other forms of ecient,
high-speed computing.
Superpipelining can be combined with other approaches. You could have a superscalar machine with
deep pipelines (DEC AXP and MIPS R-8000 are examples). In fact, you should probably expect that
faster pipelines with more stages will become so commonplace that nobody will remember to call them
superpipelines after a while.
Interestingly, the reason that the rst two are feasible is that adder units take up so little space, it is possible
to put one adder into the decode unit and another into the load/store unit. Most visualization instruction
sets take up very little chip area. They often provide ganged 8-bit computations to allow a 64-bit register
to be used to perform eight 8-bit operations in a single instruction.
Assume that (1) we are executing the load instruction, (2) R5 and R6 are already loaded from earlier
instructions, (3) it takes 30 cycles to do a oating-point divide, and (4) there are no instructions that need
the divide unit between the LD and the FDIV. Why not start the divide unit computing
the FDIV right
now, storing the result in some temporary scratch area? It has nothing better to do. When or if we arrive
at the FDIV, we will know the result of the calculation, copy the scratch area into R4, and the FDIV will
appear to executein one cycle. Sound farfetched? Not for a post-RISC processor.
The post-RISC processor must be able to speculatively compute results before the processor knows
whether or not an instruction will actually execute. It accomplishes this by allowing instructions to start
that will never nish and allowing later instructions to start before earlier instructions nish.
To store these instructions that are in limbo between started and nished, the post-RISC processor needs
some space on the processor. This space for instructions is called the instruction reorder buer
(IRB).
Figure 5.8
The IRB holds up to 60 or so instructions that are waiting to execute for one reason or another. In a sense,
the fetch and decode/predict phases operate until the buer lls up. Each time the decode unit predicts a
branch, the following instructions are marked with a dierent indicator so they can be found easily if the
prediction turns out to be wrong. Within the buer, instructions are allowed to go to the computational
units when the instruction has all of its operand values. Because the instructions are computing results
without being executed, any instruction that has its input values and an available computation unit can be
computed. The results of these computations are stored in extra registers not visible to the programmer
calledrename registers . The processor allocates rename registers, as they are needed for instructions being
computed.
The execution units may have one or more pipeline stages, depending on the type of the instruction. This
part looks very much like traditional superscalar RISC processors. Typically up to four instructions can begin
computation from the IRB in any cycle, provided four instructions are available with input operands and
there are sucient computational units for those instructions.
Once the results for the instruction have been computed and stored in a rename register, the instruction
must wait until the preceding instructions nish so we know that the instruction actually executes. In
addition to the computed results, each instruction has ags associated with it, such as exceptions. For
example, you would not be happy if your program crashed with the following message: Error, divide by
zero. I was precomputing a divide in case you got to the instruction to save some time, but the branch
was mispredicted and it turned out that you were never going to execute that divide anyway. I still had to
blow you up though. No hard feelings? Signed, the post-RISC CPU. So when a speculatively computed
instruction divides by zero, the CPU must simply store that fact until it knows the instruction will execute
and at that moment, the program can be legitimately crashed.
If a branch does get mispredicted, a lot of bookkeeping must occur very quickly. A message is sent to all
the units to discard instructions that are part of all control ow paths beyond the incorrect branch.
Instead of calling the last phase of the pipeline writeback, it's called retire. The retire phase is what
executes the instructions that have already been computed. The retire phase keeps track of the instruction
execution order and retires the instructions in program order, posting results from the rename registers to
the actual registers and raising exceptions as necessary. Typically up to four instructions can be retired per
cycle.
So the post-RISC pipeline is actually three pipelines connected by two buers that allow instructions to
be processed out of order. However, even with all of this speculative computation going on, the retire unit
forces the processor to appear as a simple RISC processor with predictable execution and interrupts.
SUBROUTINE ADDEM(A,B,C,N)
REAL A(10000),B(10000),C(10000)
INTEGER N,I
DO 10 I=1,N
A(I) = B(I) + C(I)
ENDDO
END
We have gathered these examples over the years from a number of dierent compilers, and the results are not
particularly scientic. This is not intended to review a particular architecture or compiler version, but rather
just to show an example of the kinds of things you can learn from looking at the output of the compiler.
16 This content is available online at <http://cnx.org/content/m33682/1.3/>.
17 This content is available online at <http://cnx.org/content/m33787/1.3/>.
Because there are so few registers, the variable I is kept in memory and loaded several times throughout the
loop. The inc instruction at the end of the loop actually updates the value in memory. Interestingly, at the
top of the loop, the value is then reloaded from memory.
In this type of architecture, the available registers put such a strain on the exibility of the compiler,
there is often not much optimization that is practical.
The value for I is stored in the d0 register. Each time through the loop, it's incremented by 1. At the
same time, register d1 is initialized to the value for N and decremented each time through the loop. Each
time through the loop, I is stored into memory, so the proper value for I ends up in memory when the
loop terminates. Registers a1, a2, and a3 are preloaded to be the rst address of the arrays B, A, and C
respectively. However, since FORTRAN arrays begin at 1, we must subtract 4 from each of these addresses
before we can use I as the oset. The lea instructions are eectively subtracting 4 from one address register
and storing it in another.
The following instruction performs an address computation that is almost a one-to- one translation of an
array reference:
This instruction retrieves a oating-point value from the memory. The address is computed by rst multi-
plying d0 by 4 (because these are 32-bit oating-point numbers) and adding that value to a0. As a matter
of fact, the lea and fmoves instructions could have been combined as follows:
To compute its memory address, this instruction multiplies d0 by 4, adds the contents of a1, and then
subtracts 4. The resulting address is used to load 4 bytes into oating-point register fp0. This is almost a
literal translation of fetching B(I). You can see how the assembly is set up to track high-level constructs.
It is almost as if the compiler were trying to show o and make use of the nifty assembly language
instructions.
Like the Intel, this is not a load-store architecture. The fadds instruction adds a value from memory to
a value in a register (fp0) and leaves the result of the addition in the register. Unlike the Intel 8088, we
have enough registers to store quite a few of the values used throughout the loop (I, N, the address of A, B,
and C) in registers to save memory operations.
! d3 = I
! d1 = Address of A
! d2 = Address of B
! d0 = Address of C
! a6@(20) = N
moveq #0,d3 ! Initialize I
bras L5 ! Jump to End of the loop
L1: movl d3,a1 ! Make copy of I
movl a1,d4 ! Again
asll #2,d4 ! Multiply by 4 (word size)
movl d4,a1 ! Put back in an address register
fmoves a1@(0,d2:l),fp0 ! Load B(I)
movl a6@(16),d0 ! Get address of C
fadds a1@(0,d0:l),fp0 ! Add C(I)
fmoves fp0,a1@(0,d1:l) ! Store into A(I)
addql #1,d3 ! Increment I
L5:
cmpl a6@(20),d3
bits L1
We rst see the value of I being copied into several registers and multiplied by 4 (using a left shift of 2,
strength reduction). Interestingly, the value in register a1 is I multiplied by 4. Registers d0, d1, and d2
are the addresses of C, B, and A respectively. In the load, add, and store, a1 is the base of the address
computation and d0, d1, and d2 are added as an oset to a1 to compute each address.
This is a simplistic optimization that is primarily trying to maximize the values that are kept in registers
during loop execution. Overall, it's a relatively literal translation of the C language semantics from C to
assembly. In many ways, C was designed to generate relatively ecient code without requiring a highly
sophisticated optimizer.
! a0 = Address of C(I)
! a1 = Address of B(I)
! a2 = Address of A(I)
First o, the compiler is smart enough to do all of its address adjustment outside the loop and store the
adjusted addresses of A, B, and C in registers. We do the load, add, and store in quick succession. Then we
advance the array addresses by 4 and perform the subtraction to determine when the loop is complete.
This is very tight code and bears little resemblance to the original FORTRAN code.
This is some pretty poor code. We don't need to go through it line by line, but there are a few quick
observations we can make. The value for I is loaded from memory ve times in the loop. The address of I
is computed six times throughout the loop (each time takes two instructions). There are no tricky memory
addressing modes, so multiplying I by 4 to get a byte oset is done explicitly three times (at least they use
a shift). To add insult to injury, they even put a NO-OP in the branch delay slot.
One might ask, Why do they ever generate code this bad? Well, it's not because the compiler isn't
capable of generating ecient code, as we shall see below. One explanation is that in this optimization level,
it simply does a one-to-one translation of the tuples (intermediate code) into machine language. You can
almost draw lines in the above example and precisely identify which instructions came from which tuples.
One reason to generate the code using this simplistic approach is to guarantee that the program will
produce the correct results. Looking at the above code, it's pretty easy to argue that it indeed does exactly
what the FORTRAN code does. You can track every single assembly statement directly back to part of a
FORTRAN statement.
It's pretty clear that you don't want to execute this code in a high performance production environment
without some more optimization.
This is a signicant improvement from the previous example. Some loop constant computations (subtracting
4) were hoisted out of the loop. We only loaded I 4 times during a loop iteration. Strangely, the compiler
didn't choose to store the addresses of A(0), B(0), and C(0) in registers at all even though there were plenty
of registers. Even more perplexing is the fact that it loaded a value from memory immediately after it had
stored it from the exact same register!
But one bright spot is the branch delay slot. For the rst iteration, the load was done before the loop
started. For the successive iterations, the rst load was done in the branch delay slot at the bottom of the
loop.
Comparing this code to the moderate optimization code on the MC68020, you can begin to get a sense
of why RISC was not an overnight sensation. It turned out that an unsophisticated compiler could generate
much tighter code for a CISC processor than a RISC processor. RISC processors are always executing extra
instructions here and there to compensate for the lack of slick features in their instruction set. If a processor
has a faster clock rate but has to execute more instructions, it does not always have better performance than
a slower, more ecient processor.
But as we shall soon see, this CISC advantage is about to evaporate in this particular example.
addem_:
ld [%o3],%g2 ! Load N
cmp %g2,1 ! Check to see if it is <1
bl .L77000006 ! Check for zero trip loop
or %g0,1,%g1 ! Delay slot - Set I to 1
.L77000003:
ld [%o1],%f0 ! Load B(I) First time Only
.L900000109:
ld [%o2],%f1 ! Load C(I)
fadds %f0,%f1,%f0 ! Add
add %g1,1,%g1 ! Increment I
add %o1,4,%o1 ! Increment Address of B
add %o2,4,%o2 ! Increment Address of C
cmp %g1,%g2 ! Check Loop Termination
st %f0,[%o0] ! Store A(I)
add %o0,4,%o0 ! Increment Address of A
ble,a .L900000109 ! Branch w/ annul
ld [%o1],%f0 ! Load the B(I)
.L77000006:
retl ! Leaf Return (No window)
nop ! Branch Delay Slot
This is tight code. The registers o0, o1, and o2 contain the addresses of the rst elements of A, B, and C
respectively. They already point to the right value for the rst iteration of the loop. The value for I is never
stored in memory; it is kept in global register g1. Instead of multiplying I by 4, we simply advance the
three addresses by 4 bytes each iteration.
The branch delay slots are utilized for both branches. The branch at the bottom of the loop uses the
annul feature to cancel the following load if the branch falls through.
The most interesting observation regarding this code is the striking similarity to the code and the code
generated for the MC68020 at its top optimization level:
L3:
fmoves a1@,fp0 ! Load B(I)
fadds a0@,fp0 ! Add C(I)
fmoves fp0,a2@ ! Store A(I)
addql #4,a0 ! Advance by 4
addql #4,a1 ! Advance by 4
addql #4,a2 ! Advance by 4
subql #1,d0 ! Decrement I
tstl d0
bnes L3
The two code sequences are nearly identical! For the SPARC, it does an extra load because of its load-store
architecture. On the SPARC, I is incremented and compared to N, while on the MC68020, I is decremented
and compared to zero.
Initially, the vector length register is set to N. We assume that for the rst iteration, N is greater than 128.
The next instruction is a vector load instruction into register v0. This loads 128 32-bit elements into this
register. The next instruction also loads 128 elements, and the following instruction adds those two registers
and places the results into a third vector register. Then the 128 elements in Register v2 are stored back
into memory. After those elements have been processed, N is decremented by 128 (after all, we did process
128 elements). Then we add 512 to each of the addresses (4 bytes per element) and loop back up. At some
point, during the last iteration, if N is not an exact multiple of 128, the vector length register is less than
128, and the vector instructions only process those remaining elements up to N.
One of the challenges of vector processors is to allow an instruction to begin executing before the previous
instruction has completed. For example, once the load into Register v1 has partially completed, the processor
could actually begin adding the rst few elements of v0 and v1 while waiting for the rest of the elements of
v1 to arrive. This approach of starting the next vector instruction before the previous vector instruction has
completed is called chaining. Chaining is an important feature to get maximum performance from vector
processors.
The RS-6000 also supports a memory addressing mode that can add a value to its address register before using
the address register. Interestingly, these two features (branch on count and pre-increment load) eliminate
several instructions when compared to the more pure SPARC processor. The SPARC processor has 10
instructions in the body of its loop, while the RS-6000 has 6 instructions.
The advantage of the RS-6000 in this particular loop may be less signicant if both processors were
two-way superscalar. The instructions were eliminated on the RS-6000 were integer instructions. On a
two-way superscalar processor, those integer instructions may simply execute on the integer units while the
oating-point units are busy performing the oating-point computations.
5.2.1.6 Conclusion
In this section, we have attempted to give you some understanding of the variety of assembly language that
is produced by compilers at dierent optimization levels and on dierent computer architectures. At some
point during the tuning of your code, it can be quite instructive to take a look at the generated assembly
language to be sure that the compiler is not doing something really stupid that is slowing you down.
Please don't be tempted to rewrite portions in assembly language. Usually any problems can be solved
by cleaning up and streamlining your high-level source code and setting the proper compiler ags.
It is interesting that very few people actually learn assembly language any more. Most folks nd that
the compiler is the best teacher of assembly language. By adding the appropriate option (often -S), the
compiler starts giving you lessons. I suggest that you don't print out all of the code. There are many pages
of useless variable declarations, etc. For these examples, I cut out all of that useless information. It is best
T tag, 1.1.5(13)
W wait states, 1.1.2(8)
wordy tests, 2.3.1(85)
techniques, 3.2.4(163)
workload, 2.2.1(62)
temporal locality of reference, 1.1.4(9)
write-through policy, 3.2.2(146)
thrashing, 1.1.5(13)
writeback, 3.2.2(146)
thread, 3.2.3(151)
Attributions
Collection:High Performance Computing
Edited by: Charles Severance, Kevin Dowd
URL: http://cnx.org/content/col11136/1.5/
License: http://creativecommons.org/licenses/by/3.0/
Module: "1.0 Introduction to the Connexions Edition"
Used here as: "Introduction to the Connexions Edition"
By: Charles Severance
URL: http://cnx.org/content/m32709/1.1/
Pages: 1-2
Copyright: Charles Severance
License: http://creativecommons.org/licenses/by/3.0/
Module: "Introduction to High Performance Computing"
By: Charles Severance, Kevin Dowd
URL: http://cnx.org/content/m32676/1.3/
Pages: 3-5
Copyright: Charles Severance, Kevin Dowd
License: http://creativecommons.org/licenses/by/3.0/
Module: "Memory - Introduction"
Used here as: "Introduction"
By: Charles Severance, Kevin Dowd
URL: http://cnx.org/content/m32733/1.3/
Pages: 7-8
Copyright: Charles Severance, Kevin Dowd
License: http://creativecommons.org/licenses/by/3.0/
Module: "Memory - Memory Technology"
Used here as: "Memory Technology"
By: Charles Severance, Kevin Dowd
URL: http://cnx.org/content/m32716/1.3/
Pages: 8-9
Copyright: Charles Severance, Kevin Dowd
License: http://creativecommons.org/licenses/by/3.0/
Module: "Memory - Registers"
Used here as: "Registers"
By: Charles Severance, Kevin Dowd
URL: http://cnx.org/content/m32681/1.3/
Page: 9
Copyright: Charles Severance, Kevin Dowd
License: http://creativecommons.org/licenses/by/3.0/
Module: "Memory - Caches"
Used here as: "Caches"
By: Charles Severance, Kevin Dowd
URL: http://cnx.org/content/m32725/1.3/
Pages: 9-13
Copyright: Charles Severance, Kevin Dowd
License: http://creativecommons.org/licenses/by/3.0/
About Connexions
Since 1999, Connexions has been pioneering a global system where anyone can create course materials and
make them fully accessible and easily reusable free of charge. We are a Web-based authoring, teaching and
learning environment open to anyone interested in education, including students, teachers, professors and
lifelong learners. We connect ideas and facilitate educational communities.
Connexions's modular, interactive courses are in use worldwide by universities, community colleges, K-12
schools, distance learners, and lifelong learners. Connexions materials are in many languages, including
English, Spanish, Chinese, Japanese, Italian, Vietnamese, French, Portuguese, and Thai. Connexions is part
of an exciting new information distribution system that allows for Print on Demand Books. Connexions
has partnered with innovative on-demand publisher QOOP to accelerate the delivery of printed course
materials and textbooks into classrooms worldwide at lower prices than traditional academic publishers.