23. Welcome to the Machine¶
PDF pages 1109–1120
Chapter 23 Welcome to the Machine 23.1 Finite Representations of Real Numbers Almost all digital computers store numbers as binary strings of zeros and ones. Representing an integer in binary form is straightforward assuming you understand binary vs. decimal counting. The only subtle point is whether the integer is signed or unsigned. The only difference is in the interpretation of the sign bit. For example, a 32-bit unsigned integer can range from 0 to 232 −1 = 4 294 967 295. By contrast, in a 32-bit signed integer, which ranges from −231 = −2 147 483 648 to 231−1 = 2 147 483 647, the idea is the same, but the first bit is interpreted as the sign bit (with a one representing a negative integer). That is, any integer greater than 231 −1 is simply wrapped by subtracting 232. Indeed, integer arithmetic is always performed modulo the range of the integer type, so that in unsigned 32-bit arithmetic, 4 294 967 295+1 = 0. Note that modern Fortran has its own processor-independent integer model that is somewhat more restrictive than this (i.e., integers are always signed), but Fortran integers can effectively always be used in the same ways as C integers, with the right interpretation. Fixed-point arithmetic is the simplest way to represent real numbers. Basically, fixed-point numbers are integers with a decimal place stuck in some fixed position. Thus, they are in some sense a reinterpretation of integers. Addition is straightforward, whereas multiplication is a bit more complicated. Even though fixed- point arithmetic is typically very fast, the limited dynamic range of fixed-point numbers usually leads to serious problems with accuracy, and consequently they are only used in specialized applications. The floating-point number is the standard on modern computing machines. A floating point number is decomposed into two parts, a mantissa and an exponent. For example, the real number 328.925 is represented in floating-point arithmetic in usual scientific notation as 3.28925 \times 102, (23.1) or in computer-style notation as 3.28925E2 (in Fortran the E can be changed to a D to denote a double- precision number, which we will discuss below). The mantissa and exponent are then simply represented as integers, with a fixed decimal place associated with the mantissa; note that for a binary representation, the exponent is chosen such that the mantissa is in the range [1, 2). There are many ways to decide how to associate binary data for floating-point numbers, but by far the most common standard is dictated by the IEEE Standard for Binary Floating-Point Arithmetic (ANSI/IEEE Std 754-1985). IEEE single precision allocates 32 bits for a floating-point number. Without going into too much detail, one bit represents the sign, 8 bits represent the exponent, and the rest represent the mantissa. IEEE single precision is characterized by the following values: • ‘‘machine epsilon’’ (epsilon(1.0) in Fortran 90) of 1.1920929E-7 • smallest positive number (tiny(1.0) in Fortran 90) of 1.1754944E-38 • largest number (huge(1.0) in Fortran 90) of 3.4028235E+38
23.2.1 Memory Hierarchy¶
Chapter 23. Welcome to the Machine Here machine epsilon is a measure of the precision of the representation, and is defined as the smallest number that, when added to 1, returns a number different from 1. Thus, single precision gets about 7 significant digits, and ranges through roughly 10\pm38. (Actually, nonzero values closer to zero can be represented, but are not usually accessible in Fortran.) IEEE double precision extends this range considerably, and is characterized by the following values: • machine epsilon (epsilon(1.0D0) in Fortran 90) of 2.220446049250313E-16 • smallest positive number (tiny(1.0D0) in Fortran 90) of 2.2250738585072014E-308 • largest number (huge(1.0D0) in Fortran 90) of 1.7976931348623157E+308 Thus, we get about 16 significant digits and a range through about 10\pm308. There are two perfectly reasonable but opposite philosophies considering the choice of precision. The first goes: ‘‘modern computers are so fast that you may as well just use double precision just to be safe.’’ Actually, as 64-bit processors become more common, the speed penalty for using double precision instead of single precision is shrinking. However, the other philosophy is ‘‘if you need double precision, your algorithm is probably flawed; single precision is almost always sufficient.’’ You’ll just have to decide for yourself, but it is generally useful to write your code such that it is easy to switch precision. In Fortran 90, you do this by defining a parameter wp (for ‘‘working precision’’ as
to select double precision. Then, whenever you declare a variable, use this parameter, as in declaring
14 digits of precision’’ and ‘‘at least 4 digits of precision, please’’) as the argument to selected_real_kind above. 23.2 Machine Structure and Optimization Here we will discuss just a few of the basic features of modern processors that relate to high-performance numerical computing. This is an area that is complicated, varies widely among processor vendors, and evolves rapidly, and so a general and introductory discussion of this sort must necessarily sacrifice detail and completeness.1 However, knowledge of some basic concepts is invaluable in tuning your codes for high performance. 23.2.1 Memory Hierarchy The first thing to understand is how the processor accesses, stores, and manipulates information. There is a hierarchy of locations in which information can be stored. Here, we list these in decreasing order of access speed, and conversely, increasing order of size. • Registers. The registers are the memory locations for the data on which the processor is currently working. There are typically of the order of 10 registers on a processor, and obviously they must operate at the nominal execution speed of the processor. A typical machine instruction for a floating- point operation might involve, say, taking the contents of registers 0 and 1, multiplying them together, and storing the result in register 0. In the past, a special register called the accumulator was typically the register that was always involved in any particular operation and also the destination for the result of any operation. Now the term is less common, but one or two registers are still typically of more importance than the rest. For numerical calculations, it is useful to note that on some processors, the registers handle more data than regular memory addresses. For example, even when processing 64-bit data, the registers 1For more detailed discussions, see Kevin Dowd and Charles Severance, High-Performance Computing, 2nd ed. (O’Reilly, 1998). While this reference is somewhat dated, it is clear and readable, and still contains a great deal of relevant information.
23.2 Machine Structure and Optimization might be designed be 80-bit or 128-bit ‘‘extended precision’’ registers (as happens, for example, on modern Intel processors). Thus, if you compute the sum of many numbers, the processor can keep the intermediate sums with extra precision to reduce roundoff error, particularly in cases where the sum is sensitive to the order of addition due to roundoff errors. This also helps reduce unexpected effects due to unintended reordering of operations at the compiler or processor levels. • Cache. Next comes a small area of high-speed memory called cache memory. The idea is that the main memory, where large quantities of data can be stored is quite slow, but the processor needs to take in and push out data very quickly. Cache is the intermediate area that the processor can use on short time scales. Think of it as a small pile of papers on your desk that you want quick access to, compared to main memory, which is more like the filing cabinet. Because this memory is so fast, it is quite expensive, and thus its size is quite limited. Modern designs also use multiple levels of cache: L1 (‘‘level 1’’) cache is the fastest and smallest, typically being on the order of a few to a few hundred KB, and talks directly to the cpu; L2 cache is larger and slower, typically on the order of a few hundred KB or larger; and some designs even incorporate an L3 cache. Cache is a particularly important concept in modern computers, because the processors must not be kept waiting for the information they need. Of course, the cache can only keep a small subset of the total data in main memory on hand, but if it doesn’t have a particular piece of data when the processor requests it, the computation stalls while the data are fetched from a higher-level cache or even main memory. One fairly obvious strategy that helps here is to process data in fairly small-sized chunks such that ‘‘blocks’’ of the calculation can fit entirely in cache. Of course, most calculations where you care about speed will not fit in cache, and manually breaking up calculations into cache-sized chunks is difficult and guaranteed to render your code unreadable or at least ugly. The other strategy requires a bit more understanding about how cache works. Essentially, the various elements of cache are copies of various elements of main memory. However, if each location in cache were a copy of a completely independent location in main memory, we would need another bank of fast memory, the same size as the cache, just so we would know which main memory entry each cache entry referred to. To reduce this memory overhead, cache elements are grouped into cache lines, so that when one datum from memory is needed, a whole line’s worth of data are actually fetched from memory. Thus, as long as you have to fetch all the data together, you may as well make the best possible use of them. The basic strategy is this: stick to unit-stride access as much as possible. That is, if you are processing long arrays of data, try to access the elements only sequentially, if at possible. The canonical example here is in computing the sum of a matrix. In Fortran, a two-dimensional array A is stored in memory by column, or such that A(1,1) and A(2,1) are adjacent in memory, while A(1,1) and A(1,2) are separated in memory by at least the length of the first array dimension. Thus, for example, what might seem a reasonable method for computing the array sum, s = 0 do j = 1, m do k = 1, n
end do end do is actually a bad idea, because the access in memory in the inner loop has a stride of (at least) m. In the worst, case, an entire line of cache must be fetched for each addition operation, slowing things down considerably. Fortunately, this problem is easily fixed by switching the order of the loops: s = 0 do k = 1, n do j = 1, m
end do end do
Chapter 23. Welcome to the Machine The access here is now unit-stride, and this makes optimal use of the cache since (almost) all the fetched data are used in the sum. Most compilers will detect the problem in the former code sample and change it to the latter, depending on the level of optimization you request (after all, you might want the former calculation due to some issue with roundoff errors, since the results of the two codes are not guaranteed to be identical). In a more modern approach, as in Fortran 90, you can simply use an intrinsic such as sum, as in
or
to accomplish the same calculation, but explicitly giving the compiler freedom to choose the best way to perform the sum. If an intrinsic does not exist to do what you want, there are other constructs such as the forall loop in Fortran 95, which can be used to indicate that there is no dependency among the iterations of multiple loops, so that the compiler can perform the operations in any (presumably optimal) order. The problems with cache can be even a bit more insidious than what we have indicated. To see why, consider a slightly different code that performs a sum over the second dimension of an array: dimension A(2048, 128) do j = 1, 2048 do k = 1, 128
end do end do If A is stored contiguously in memory, then the elements A(j,k) for the same j but different k are separated by powers of two. But to perform the calculation, lines of cache are fetched, corresponding to A(j,k) for the same k but different j, which is not so useful. In a fully associative cache, a line in cache can be associated with any line in memory, and many lines of cache can be fetched on the first iteration of the j loop, so that they will be reused later on subsequent iterations. However, to reduce complexity and cost, most cache is set-associative, which means that a line in cache can only map to particular lines in memory, and conversely a line in memory can only map to a few (say two or four) cache lines. Typically, a cache line maps to locations in memory that are widely spaced by some number of lines given by a power of two. The problem in the above example is that the full array itself has as its first dimension a large power of two. All of the A(j,k) for fixed j are needed at the same time, but many will overlap to the same few lines of cache. A line of cache has to be fetched for each k, but they all can’t be stored in cache. So the same lines need to be refetched on subsequent iterations of the j loop, and in the worst case, a line of cache must be fetched for each addition operation. So even if the relevant data could have fit in cache, the power-of-two associativity caused fetched data to be flushed before it was needed, so that the cache needed to fetch it again. This behavior is called cache thrashing and can be highly detrimental to the performance of your code. In addition to unit stride access, cache thrashing can sometimes be avoided (as in this example) by padding the first array dimensions to some larger (non-power-of-two) value, to avoid problems with associativity. (An appropriate change in the example would be 2048 −\rightarrow 2048 + 128.) In general, you should avoid accessing data with large, power-of-two strides (some traditional FFT algorithms are notorious for doing just this). • Main Memory. Most of your data reside in main memory during a calculation, which can be large (in the range of tens of GB on the best single-box machines). However, it is slow, which is the point of having cache. The strategies we discussed for cache also apply here, since main memory is often also banked. That is, sequential chunks in memory can come from different banks of memory to reduce
23.2.2 Pipeline¶
23.2 Machine Structure and Optimization latency when fetching data quickly. The number of banks is typically a power of two, so power-of-two strides are bad here as well, and clearly unit stride is best: after accessing one bank for data, you would then access the next bank for the next piece of data, giving the first some time to recover and prepare for the next fetch before you bother it again. • Storage. Of course, memory for even more data in the longer term comes in the form of disks, tapes, and so on. These are really slow, and if you have to use slow storage on the fly because your computation is so large then you’re really in trouble. 23.2.2 Pipeline One crucially important concept in understanding modern, ultrahigh-speed processors is the processor pipeline. Let’s consider an analogy to motivate the pipeline. Suppose you have a ‘‘clothes-cleaning ma- chine.’’ It’s quite nice, you just pop in a load of clothes, and then it proceeds to wash, dry, iron, and fold them for you. Suppose each of these tasks takes the machine 15 minutes. Then the rate at which the machine washes clothes is 1 load/hour. Of course, that’s not very efficient, since while the clothes are being dried, the machinery related to washing the clothes sits idle. But if you have a lot of laundry to do, you’d want a more clever design for a machine would divide the machine into 4 units, one for each task. Each task still takes 15 minutes, and after you put in your first load of laundry, it’s an hour before you see a set of spanking-fresh linens. But after the first load finishes washing and goes to the dryer, you can start the second load of clothes right away. In all you can have a total of 4 loads of laundry ‘‘in flight,’’ and after a delay of one hour, you effectively are finishing one load of laundry every 15 minutes, or four times faster than the first machine. This is the essence of pipelining: dividing the work up into stages, so that several different instructions can be processed in different parts of the ‘‘assembly line’’ at once. It should be apparent that dividing up tasks into smaller pieces allow for longer pipelines, and thus for faster processors. Indeed, the recent offerings from Intel bear this out, with pipelines of 10 stages for the Pentium III, 20 stages for the Pentium 4, 31 stages for the Xeon Nocona (with processor speeds currently in the high 3 GHz range). By contrast, the older 8086 used no pipeline at all, and executed one instruction per cycle with no latency. The crucial point here is this: the pipeline is your friend only when full. With a full 31-stage pipeline, it appears (after a delay of 31 clock cycles) that one instruction is being executed per cycle. However, suppose that successive instructions depend on each other, so that the next instruction can’t be started until the
to implement a random-number generator). In this case, it takes 31 clock cycles to execute each operation, because no pipelining is possible, and you’ve effectively just cut the processor speed down by a factor of 31. That’s bad. In fact, you can see why cache thrashing is even worse on a fast processor: not only does the processor have to wait for the memory fetch, but the problem could be compounded if no other instructions are ready to go, since the pipeline will drain. Modern processors implement a bunch of tricks, many of which are aimed at keeping the pipeline packed with instructions, to keep the performance high. We will discuss them below. But it is important to note that with speed limitations on memory (especially main memory), and with constraints on instruction sets, and so on, it is typically very hard to keep a processor working at its ‘‘theoretical peak.’’ For example, Intel Xeon processors can theoretically execute 1 floating-point operation (flop) per processor cycle, with register, cache, and other overhead it is often difficult to execute flops on more than, say, 30% of the processor cycles. You should definitely keep this in mind when tuning your codes with hardware performance counters (see below), so you don’t have unrealistic goals for your flop counts. There is a particular class of processor, the vector processor (as opposed to the above cache-based, or scalar, processors), that is optimized for just one task: take contiguous chunks of data, and perform the same mathematical operation on all the data, and do this quickly. These are very good for scientific computation, again provided you use unit strides in your calculation. In these processors, it is much more realistic to achieve flop rates nearing the 100% ideal. However, these days, such processors are expensive and normally relegated to the best supercomputers (i.e., computers that aren’t sitting on your desktop). So it’s important to learn to deal with the constraints of the cheaper (and ubiquitous) scalar processor.
Chapter 23. Welcome to the Machine However, it is difficult to do much in a code to avoid ‘‘bubbles’’ in the pipeline, beyond what you would already do to make a ‘‘cache-friendly’’ code. This is especially true since pipelines vary greatly among processors, with most processors having multiple pipelines (say, to handle integer and floating-point instructions separately). You have to rely heavily on good compilers to provide a mix of instructions to the processor without a lot of dependencies to keep the pipelines working. Some rather simple tricks, like Intel’s ‘‘hyperthreading technology,’’ rely on executing multiple codes at once, so the different codes fill in each others’ pipeline bubbles, so at least the processor is staying more busy overall, even if each code is not executing faster. 23.2.2.1 Out-of-Order Execution One of the tricks implemented by all of the fastest modern processors is out-of-order execution. The concept is fairly simple: as the machine instructions are being sent to the pipelines for execution, they are first held in a buffer (of something like 100 instructions), and then analyzed for dependencies. If the processor detects an instruction that depends on the result of another instruction (either also in the buffer or already in the pipeline), it is free to dynamically reorder the instructions so that the dependent instruction is moved back in the queue so the pipeline isn’t idle while the instruction waits. This helps the execution speeds greatly, but makes it difficult to analyze how a set of machine instructions will actually be executed (especially if you’re trying to hand-code an optimized routine in machine language, or trying to disassemble some compiled code to analyze its performance). Again, there isn’t much for you to do here, you have to hope your compiler is good enough to provide a good mix of instructions for the out-of-order buffer to work with. 23.2.2.2 Loop Unrolling Very often, you need to perform a repeated computation on multiple elements on an array: do j = 1, n
end do The key idea here is that at each loop iteration, the processor must execute a branch (if/then) instruction to decide if the next iteration should be performed. This is particularly bad if the processor waits for each branch instruction to be carried out before starting the next loop multiplication, since effectively all benefits of having a pipeline are lost. (Branch prediction, described below, helps this somewhat.) Thus, it would help to rewrite the loop as do j = 1, n, 4 a(j)
- b(j)
end do This loop is said to have been ‘‘unrolled four times.’’ Now there is no branch separating the instructions in the loop, and four iterations of the original loop can be pipelined right away without any branch overhead. Of course, the code here assumes that n is a multiple of 4; the unrolled loop is more complicated for arbitrary n. The tradeoff here is that the resulting code is larger than the original. The benefits also decrease as the loop is unrolled more, so excessive loop unrolling is not useful. Generally, this is handled by the compiler optimizer, and not in your code, so your code stays readable. However it is useful to know that a loop can be unrolled only if its iterations are independent. In certain constructions, for example with pointers, the compiler may not be able to ‘‘prove’’ that the iterations are independent, and thus not optimize it. For example, consider this loop: real, dimension(n) :: a, b integer, dimension(n) :: c
23.2 Machine Structure and Optimization do j = 1, n
end do If you know that the array c is a permutation of the set {1, . . . , n}, then the loop can still be unrolled. However, a compiler would likely assume that values in the c array could have been repeated, and thus not optimized the loop. In this case, a compiler directive would help by telling the compiler that it can assume the loop iterations to be independent. In the HPF (High-Performance Fortran) language, this would look like this: !hpf$ independent do j = 1, n
end do This also occurs in the original loop, if it is in a subroutine (this is valid Fortran 90): subroutine foo(a, b, n) integer :: n, j real, dimension(n), intent(inout) :: a, b do j = 1, n
end do end subroutine foo The problem could arise in some languages if the two arrays overlap in memory, for example if the two arguments are overlapping parts of the same array, as in call foo(c(2:6), c(1:5), 5). In this case, the results depend on the order in which the loop iterations are executed, since the value of elements of b are changing. Actually, this problematic call is explicitly disallowed in Fortran 90: array arguments must not overlap if they are defined or modified by the subroutine. Fortran 90 tends to make choices to favor optimization over flexibility, and the loop could be unrolled by the compiler in this example. However, in most other languages (like C), the call would be acceptable and thus the compiler would not unroll the loop. 23.2.2.3 Branch Prediction As we mentioned above, branch instructions (i.e., instructions to jump to different parts of a program depending on some condition, as in an if or case statement) are problematic for pipelined machines, since in principle instructions after the branch can’t be executed until the result from the branch condition is known. Actually, you could start executing instructions before the branch result is known: just pick one outcome, and start executing the appropriate instructions, hoping to win if the processor guesses the right outcome in advance. This trick is called speculative execution. The problem is that canceling (or retiring) the finished and in-flight instructions in the case of a wrong guess involves a lot of overhead, and unless the processor has a good way to accurately guess the result of the branch condition, speculative execution could actually slow things down. Thus enters the art of branch prediction. Again, the prediction in modern processors must be good, since in some processors the cost of a mispredicted branch is stalling and flushing the entire pipeline to clear the false branch. Branch-prediction algorithms are numerous and can be quite complex2 As a simple example, we’ll consider branch prediction with a dynamically updated Moore machine.3 We’ll do this for a three-bit machine, so consider the following table of ‘‘addresses’’ three bits long: 2For nice coverage see the Wikipedia entry ‘‘Branch Predictor,’’ http://en.wikipedia.org/wiki/Branch_predictor. 3Edward F. Moore, ‘‘Gedanken-experiments on Sequential Machines,’’ Automata Studies (Annals of Mathematical Studies) 34, 129 (1956).
Chapter 23. Welcome to the Machine b1 b2 b3 output x x x x x x x x The outputs are initially undetermined (corresponding to some default value of, say, 0). The 0’s correspond to ‘‘branch true,’’ and the 1’s correspond to ‘‘branch false.’’ Then the output is the prediction given the last three branch results. That is, if the last three branches were true, false, and false, then our prediction for the next branch would be the ‘‘100’’ output. Correspondingly, after the real branch result is known, the result is recorded at the same place in the table. Clearly, if the branch result is always the same thing, this algorithm will predict it perfectly after four possible mispredictions (settling into either the 000 or 111 output). Suppose now that branches come in some more complicated pattern like 1000100010001000 . . ., and suppose we always default with 0. Shortly, the table will settle down to this: b1 b2 b3 output x x x x That is, after a transient of 4 branches, some of which are default mispredictions, the predictor predicts the pattern perfectly. Larger predictors will obviously predict longer-period patterns, but may take longer to ‘‘lock on.’’ Basically, regular branch patterns are easy to predict, whereas any change in pattern (e.g., at the end of a loop) or branching on a random number can seriously hurt your performance if done too often. 23.2.2.4 Addition, Multiplication, and Division One more thing we can discuss are the basic types of floating-point operations that you want to do. Additions are typically the easiest, with multiplications next, and divisions are by far the hardest. Most modern processors are set up to churn out (with pipelining) one floating-point operation per cycle, such as one addition. Many also can do one floating-point multiplication per cycle (per pipeline). However, there are many variations on this theme. Modern Intel Xeon processors can only pipeline one multiplication every other cycle; however, on the cycles between multiplications, you are allowed to pipeline an addition with no extra cost. That’s a sense in which additions are ‘‘cheaper’’ than multiplications (additions can be pipelined on every cycle). Other processors are set up to do one multiplication per cycle, but they can also do an addition at the same time as a multiplication in a single ‘‘multadd’’ operation. Thus, mixing floating-point multiplications and additions together can take advantage of hardware capabilities and result in very good efficiency. Divisions are extremely bad: for example, a double-precision division on an Intel Xeon can be finished every 38 processor cycles, compared to 2 for multiplication. Many optimizing compilers have options to substitute computing b−1 and then multiplying by a to compute a/b, although the result may be less accurate than the divide. (Similarly, it is useful to know that it is usually possible to compute inverse square roots very quickly, and some processors have special instructions to compute the sine and cosine of a number at the same time.) It is usually best to keep your code readable and let the compiler make the appropriate transformations: most compilers can easily change a statement like a/2 to 0.5*a. However, compilers occasionally miss operations, and these concepts can be useful in speeding things up.
23.2.3 Avoiding Overhead¶
23.2 Machine Structure and Optimization 23.2.3 Avoiding Overhead As we mentioned above in loop unrolling, we can greatly increase performance by decreasing overhead operations. Here we will briefly discuss a couple of situations where it is possible to profitably decrease overhead. 23.2.3.1 Procedure Inlining One place where much overhead can be eliminated is in the calling of a procedure (a function or a subroutine). Whenever a procedure is called, the code must jump to a new section, with some provisions for where to return at the completion of the procedure, as well as possible setup of temporary variables and arrays. If the program spends a lot of time in the procedure, then the overhead may not be a big deal. However, consider the loop in this code example: do j = 1, n
end do If foo and bar are relatively small functions, then the function-call overhead can be quite substantial. Further- more, the computations in the two procedures cannot be rearranged to improve performance (i.e., instructions from one procedure could fill a pipeline bubble in the other), if the procedures are compiled separately. An easy solution is to simply take the contents of the two procedures and paste them directly into the loop to eliminate the function calls. This trick is called procedure inlining. Obviously, manually inlining procedures will make your code a whole lot less readable, and inlining almost always increases the size of your code. Inlining is best done at the compiler level, and even then it should be limited to certain cases, such as procedure calls in the innermost loops of your code. For example, if inlining expands the code for a loop to the point where it no longer fits into instruction cache (if there are many copies of a procedure), then inlining may slow things down. It is also worth noting that compilers have trouble inlining under certain conditions. For example, it is difficult for a compiler to inline a procedure when it is defined in a different file from which it is called (e.g., when it is in a library). To inline such separated procedures, it must defer most of the optimizations until the linking phase, and a number of compilers now do this. 23.2.3.2 Compiler Issues Compiler vendors must work very hard to get their codes to implement the tricks we have discussed here (as well as many, many more). If your primary goal is speed at any cost you should note this: simple, old features of a language will be the best supported in a compiler in terms of optimization or just plain working, while using new, advance, fancy-schmancy features in your code will tend to inhibit optimization or even cause you to send a bug report to the vendor! For example, suppose you are integrating the Schrödinger equation, and in your subroutine to evaluate the time derivative of the wave-function array psi, you try to keep things organized by using a function call to compute the result of applying an operator A on \psi in one of the terms of the Schrödinger equation \partial t\psi = A\psi + \cdot \cdot \cdot :
In Fortran 90, this is an allowed construction: a function can return an array as a result, which here is then added to the psidot array. However, even otherwise excellent compilers can miss the fact that a temporary array to store the result of the function call to A_on_psi can in fact be eliminated. If this procedure is called often (as is likely in a code that evolves the Schrödinger equation), it is likely to speed things up by using a much less elegant subroutine call, which adds the term A\psi to psidot: call add_A_on_psi(psidot) The point is that the procedure should be inlined in either case, but when the new feature (array-valued functions), the optimizer support won’t be as good. The same goes for user-defined data types, pointers, and virtually anything object-oriented. Of course, there are good reasons for the existence of these features: they
23.2.5 Tuning Your Code¶
Chapter 23. Welcome to the Machine make your code more readable and elegant. There is usually a tradeoff between speed and general elegance when programming. 23.2.4 Parallel Programming Optimizing parallel codes is an advanced topic, and way beyond the scope of this discussion, even as it becomes more relevant with multicore computers and wider availability of computing clusters. We’ll just stick to crude, obvious observations here. Basically, parallel programming involves communication between processors. In the best case of symmetric multiprocessor (SMP) computing, the processors are all part of the same computer, or ‘‘box,’’ and share the same memory. Communication is simple, since the processors have access to the same pool of memory. In the worst case of parallel processing, as in many clusters, the computers must talk over standard networks. In any case, communication is much slower than the crunching that happens on a single computer. Your goal is to minimize the communication, or at least organize it to make it as efficient as possible (e.g., into occasional, large ‘‘bursts’’ to minimize effects of network latency). Parallel programming is complex and subtle. The quantum optician who doesn’t wish to spend all of his/her time writing code would do well to learn a data-parallel language such as High-Performance Fortran (HPF), where in many cases a well-written Fortran 90 code can be changed to a parallel code only by adding some compiler directives as comments in the Fortran 90 file. 23.2.5 Tuning Your Code So, how do you increase the speed of your code? We’ll only say a couple of superficial things here to get you started. The most important thing to do is to get a good compiler, and investigate all of its options. Good compilers support the above optimizations and many more, but usually only if you explicitly enable them. Good compilers will also let you write nice-looking, readable code without worrying about speed. For example, in this loop, do j = 1, n
end do it is possible to save on almost half of the multiplications by precomputing the constant product 2*n (the loop invariant):
do j = 1, n
end do However, this slightly obfuscates things, and anyway almost any worthwhile compiler will take care of this for you. The basic strategy for writing a good, fast code is this: • Write a readable, working code, not worrying about performance. • Use a profiler, a program that runs your code and tells you where your code is spending its execution time. Usually it will give you something like a percentage of time in each subroutine. A free and commonly available program is gprof. • If you identify one or a few routines in which the code is spending most of its time, then concentrate on optimizing those. If you can’t identify any such routine, it may not be worth the effort. (If you work very hard to double the speed of a routine that accounts for 1% of the total execution time, the net speedup will be truly unimpressive!) • If you have to resort to obfuscating a routine to speed it up, consider maintaining two versions: a readable one and a fast one. Or at least put the readable version in the comments.
23.2 Machine Structure and Optimization • Make sure to test each change you make, to make sure it is really helping. If possible, you should also use hardware performance counters to diagnose and characterize the performance of your code. Most processors have counters for performance-related events such as floating- point operations, cache misses, and so on. Unfortunately, tools for accessing these counters are often less convenient and less available than their profiling counterparts (in some operating systems, accessing the counters requires kernel-level access and thus recompiling the kernel to enable them).