An afternoon discussion at work among some of our interns revolved around a question that one (A) posed to another (B): "How many molecules of Abraham Lincoln are in a glass of water?" After sufficient simplifying assumptions, B (a mathematician) was able to solve the problem with a simple computation. The process resulted in a debate between A and B over the virtues of being able to derive formulae and constants, rather than looking them up with a search engine. I sided with B, of course! as having and exercising the ability to carry out that process is much more useful than demanding that a search engine do all the work. Simple units and order-of-magnitude calculations are things my PhD advisor does very rapidly, and when I was his student, I always found myself lagging behind and feeling not so intelligent. Exercising that skill at least saved me a little embarrassment!
Later that evening, A was explaining to me how one writes applications for a cell phone. I found myself feeling perfectly happy that I didn't have to do all that work. Cell phones are things I want to use to make phone calls, and maybe look up directions in an emergency; I struggled to consider why I would want to write my own app. However, then I realized I was doing the same thing that A had done above: not valuing the process over the answer. Learning something about writing a cell phone application might serve me well at some future time, even though I can't imagine the value of this time investment now. I suppose I was being a mathematical chauvinist!
21 July 2011
11 October 2010
Git: recovering from an incorrect e-mail address
Our big software project uses Git (the version control system) in a CVS mode, with a single central repository to which we push commits (after they pass the check-in test suite). This can be good, because it forces us to fix small but annoying things -- like setting the correct e-mail address! I was testing on a different machine than my usual development workstation; I made some commits there, and then pushed to my workstation's local repository (Git shines for stuff like this). Then, I activated the check-in tests, but the push failed: my e-mail address wasn't set! (The tests check for that -- handy!)
While I had set my e-mail address for Git correctly on my workstation, I hadn't set it on the other machine. There were seven commits sitting on my workstation (where I run the check-in test suite) with the wrong e-mail address. Here's how I started the fix. First, I ran
$ git rebase -i HEAD~7
which let me fix up the last 7 local commits. "-i" means "interactive," so Git fired up a text editor and let me decide which commits to change (and how). I then repeated the following three tasks seven times:
$ git commit --amend --reset-author
(Edit the commit message -- no author e-mail address here, but it gets fixed due to "--reset-author")
$ git rebase --continue
That fixed it! There's probably an easier way, but I was pleased.
While I had set my e-mail address for Git correctly on my workstation, I hadn't set it on the other machine. There were seven commits sitting on my workstation (where I run the check-in test suite) with the wrong e-mail address. Here's how I started the fix. First, I ran
$ git rebase -i HEAD~7
which let me fix up the last 7 local commits. "-i" means "interactive," so Git fired up a text editor and let me decide which commits to change (and how). I then repeated the following three tasks seven times:
$ git commit --amend --reset-author
(Edit the commit message -- no author e-mail address here, but it gets fixed due to "--reset-author")
$ git rebase --continue
That fixed it! There's probably an easier way, but I was pleased.
04 October 2010
Larval sparse matrices?
John Rose's post on "larval objects in the JVM" suggests a type distinction between "larval" and "adult" objects. Larval objects are mutable, suitable for in-place updates and for use as scratch computation, but are not safe for sharing by multiple threads in the same shared-memory space. Adult objects are immutable and safe to share in parallel, but "changing" them requires copying, which can be expensive. His examples mainly involve small objects, like complex numbers or dates. However, he does suggest a "large-scale" example of a database, where the default state is immutable ("adult"), but a mutable API can be made available by request.
This pattern of a "larval" initialization stage -- too complex and incremental to handle in a constructor -- followed by a "publish" command (we call it "fillComplete") that transforms the data structure into its "adult" state -- is more or less how sparse matrices are used in practice. The data for the sparse matrix is processed in "element" form (more or less equation by equation), and then "assembled" into the optimized data structure. Changing the optimized final data structure can be expensive and may require re-optimization. Thus, we like to think of it as immutable when optimized.
The problem is, the sparse matrix API contains operations for the "larval" stage, which are either invalid or expensive when the matrix is in its "adult" stage. Rather than discouraging users from doing the wrong thing, it would be nice for the API to forbid invalid operations syntactically (not by throwing an exception and frustrating new users with a mysterious error message), and also make clear which operations are allowed but undesirable due to their cost. Making it syntactically harder to do suboptimal things is a good way to intimidate nonexpert programmers from doing them!
John mentions the library solution to this problem: "builder" objects or functions, like Java's StringBuilder. They take as input some number of updates, applied incrementally. When the user signals that all updates are done, they produce a "final" immutable object. In the case of sparse matrices, that would mean users can only call mutating methods on the "SparseMatrixBuilder" object. To mutate a "finished" sparse matrix, they would have to request a "builder."
The problem with this approach is the ability to request a builder, after the original completion of the immutable object. It means that immutable objects can change state. It also means that the immutable object (which is now only immutable in terms of its API) and the mutable "view" coinhabit the same scope. Now the state of the "immutable" object depends on the semantics of view updates. Do they happen as soon as the update methods are called? (Simultaneity doesn't mean much in the context of multiple threads.) Do they get accumulated and applied in batches? Do they involve an asynchronous call to some remote compute device, like a GPU? Programmers shouldn't be depending on these sorts of semantics. It's better to avoid the mess and forbid changing immutable objects back into mutable ones, or return to the original solution and provide only a single object type with immutable and mutable states.
What's missing here is a "rebinding syntax" that prevents programmers from accessing the supposedly immutable object while it's in a mutable state. This is a good use case for the "WITH-" Lisp idiom, or for Python's context guards: the mutable view created by the "with" isn't allowed to escape its scope.
This pattern of a "larval" initialization stage -- too complex and incremental to handle in a constructor -- followed by a "publish" command (we call it "fillComplete") that transforms the data structure into its "adult" state -- is more or less how sparse matrices are used in practice. The data for the sparse matrix is processed in "element" form (more or less equation by equation), and then "assembled" into the optimized data structure. Changing the optimized final data structure can be expensive and may require re-optimization. Thus, we like to think of it as immutable when optimized.
The problem is, the sparse matrix API contains operations for the "larval" stage, which are either invalid or expensive when the matrix is in its "adult" stage. Rather than discouraging users from doing the wrong thing, it would be nice for the API to forbid invalid operations syntactically (not by throwing an exception and frustrating new users with a mysterious error message), and also make clear which operations are allowed but undesirable due to their cost. Making it syntactically harder to do suboptimal things is a good way to intimidate nonexpert programmers from doing them!
John mentions the library solution to this problem: "builder" objects or functions, like Java's StringBuilder. They take as input some number of updates, applied incrementally. When the user signals that all updates are done, they produce a "final" immutable object. In the case of sparse matrices, that would mean users can only call mutating methods on the "SparseMatrixBuilder" object. To mutate a "finished" sparse matrix, they would have to request a "builder."
The problem with this approach is the ability to request a builder, after the original completion of the immutable object. It means that immutable objects can change state. It also means that the immutable object (which is now only immutable in terms of its API) and the mutable "view" coinhabit the same scope. Now the state of the "immutable" object depends on the semantics of view updates. Do they happen as soon as the update methods are called? (Simultaneity doesn't mean much in the context of multiple threads.) Do they get accumulated and applied in batches? Do they involve an asynchronous call to some remote compute device, like a GPU? Programmers shouldn't be depending on these sorts of semantics. It's better to avoid the mess and forbid changing immutable objects back into mutable ones, or return to the original solution and provide only a single object type with immutable and mutable states.
What's missing here is a "rebinding syntax" that prevents programmers from accessing the supposedly immutable object while it's in a mutable state. This is a good use case for the "WITH-" Lisp idiom, or for Python's context guards: the mutable view created by the "with" isn't allowed to escape its scope.
24 September 2010
Using ScaLAPACK: not so painful!
I finally got an excuse this week to try out ScaLAPACK, a library for distributed-memory parallel dense linear algebra computations. I've mainly been working on iterative methods for sparse problems, so it's not common for me to want to factor some giant dense matrix. However, part of those iterative methods involves an orthogonalization method for small groups of dense vectors. Comparing against ScaLAPACK's QR factorization PDGEQRF should make a good sanity check for accuracy and performance. (PDGEQRF should be slower, since it uses an algorithm that requires more communication, and also since it likely doesn't use multicore parallelism as effectively for the panel factorization.)
It took me two and a half days to get a prototype accuracy test and benchmark up and running, which is a lot less than I had thought. The Intel MKL Link Line Advisor was a big help in linking the right libraries. I've tried building ScaLAPACK from source before -- while that was some four years ago, it was a real pain and I was happy not to have to try that. The last few hours of those two-and-a-half days were spent trying to figure out why the matrix distribution wasn't working. The routine in question, PDGEQRF, either refused to accept the input data (it helpfully identified the first argument of the routine which was off in some way, but didn't say why it was wrong), or only factored part of the matrix. Staring at the ScaLAPACK Users' Guide -- in particular, at the diagrams illustrating global matrix layout -- helped me decode the mysterious parameters describing how the matrix is distributed among processors.
Using ScaLAPACK turned out not to be so difficult as I had thought. The Users' Guide is well-written, with many helpful diagrams. The routines have great documentation (and helpfully repeat documentation for certain library-standard interfaces, like the matrix descriptor array). Having a pre-built library with the linker flags ready for copy-and-paste into my build configuration system was a huge help. I imagine much of the pain people have reported came from trying to build everything themselves. Note also that Matlab now has support for distributed-memory parallel computation, and distributed-memory data layouts. It's based on ScaLAPACK, which means the MathWorks' developers have gone through a lot of this pain for you ;-)
A final interesting point: the ScaLAPACK Users' Guide was written when High Performance Fortran (HPF) was under serious consideration as the way to write distributed-memory parallel code. ScaLAPACK doesn't use HPF, but it does use HPF as an alternate notation to describe different ways to distribute matrices over processors. HPF proved unpopular because it was too hard to build an efficient compiler for it, but parts of HPF made their way into modern versions of the Fortran standard.
It took me two and a half days to get a prototype accuracy test and benchmark up and running, which is a lot less than I had thought. The Intel MKL Link Line Advisor was a big help in linking the right libraries. I've tried building ScaLAPACK from source before -- while that was some four years ago, it was a real pain and I was happy not to have to try that. The last few hours of those two-and-a-half days were spent trying to figure out why the matrix distribution wasn't working. The routine in question, PDGEQRF, either refused to accept the input data (it helpfully identified the first argument of the routine which was off in some way, but didn't say why it was wrong), or only factored part of the matrix. Staring at the ScaLAPACK Users' Guide -- in particular, at the diagrams illustrating global matrix layout -- helped me decode the mysterious parameters describing how the matrix is distributed among processors.
Using ScaLAPACK turned out not to be so difficult as I had thought. The Users' Guide is well-written, with many helpful diagrams. The routines have great documentation (and helpfully repeat documentation for certain library-standard interfaces, like the matrix descriptor array). Having a pre-built library with the linker flags ready for copy-and-paste into my build configuration system was a huge help. I imagine much of the pain people have reported came from trying to build everything themselves. Note also that Matlab now has support for distributed-memory parallel computation, and distributed-memory data layouts. It's based on ScaLAPACK, which means the MathWorks' developers have gone through a lot of this pain for you ;-)
A final interesting point: the ScaLAPACK Users' Guide was written when High Performance Fortran (HPF) was under serious consideration as the way to write distributed-memory parallel code. ScaLAPACK doesn't use HPF, but it does use HPF as an alternate notation to describe different ways to distribute matrices over processors. HPF proved unpopular because it was too hard to build an efficient compiler for it, but parts of HPF made their way into modern versions of the Fortran standard.
15 September 2010
The user experience of lazy automatic tuning
I've been doing some performance experiments with an MPI-parallel linear algebra kernel. I've implemented the kernel in two different ways: a "reduce" (not quite MPI_Reduce(), though it has a similar communication pattern) and "broadcast," and a "butterfly." The latter is more general (it can be used to compute things other than what I'm benchmarking). However, it might be slow in certain cases (e.g., if the network can't handle too many messages at once, or if there are many MPI ranks per node and the network card serializes on processing all those messages going in and out of the node). Each run involves testing two different implementation alternatives, and four different scalar data types (float, double, complex float, and complex double).
When I ran the benchmarks, I noticed that the very first run (butterfly implementation strategy, "float" data type) of each invocation of the benchmark is 10 - 100x slower than for other data types with the same implementation strategy. When I rebuild the executable using a different MPI library, this problem went away. What was going on? I noticed that for the "slow" MPI library, the minimum benchmark run time was reasonable, but the maximum run time was a lot more. After reading through some messages in the e-mail archives of that MPI library, I found out that the library delays some setup costs until execution time. That means the first few calls of MPI collectives might be slow, but overall they will be fast. The other MPI library doesn't seem to be doing this.
I was reminded of some troubles that the MathWorks had with integrating FFTW into Matlab, a few years ago. FFTW does automatic performance tuning as a function of problem size. The tuning phase takes time, but once it's finished, running any problem of that size will likely be much faster. However, users didn't realize this. They ran an FFT once, noticed that it was much slower than before (because FFTW was busy doing its thing), and then complained. One can imagine all sorts of fixes for this, but the point is that something happened to performance which was not adequately communicated to users.
Now, MPI users are likely more sophisticated programmers than most Matlab users. They likely have heard of automatic run-time performance tuning, and maybe even use it themselves. I'm fairly familiar with such things and have been programming in MPI for a while. Yet, this phenomenon still puzzled me until I poked around through the e-mail list archives. What I was missing was an obvious pointer to documentation describing performance phenomena such as this one. The performance of this MPI library is not transparent. This is not a bad thing in this case, because the opaqueness is hiding a performance tuning process that will make my code faster overall. However, I want to know this right away, so I don't have to think about what I could be doing wrong.
See, I'm not just calling MPI directly. I've wrapped up some other wrapper over MPI. I didn't know if my wrapper was broken, the wrapper underneath was broken, the MPI library was broken, the cluster job execution software was broken, or if nothing was wrong and this was just expected behavior. That's why I want to know right away whether I should expect significantly nonuniform performance behavior from some library.
When I ran the benchmarks, I noticed that the very first run (butterfly implementation strategy, "float" data type) of each invocation of the benchmark is 10 - 100x slower than for other data types with the same implementation strategy. When I rebuild the executable using a different MPI library, this problem went away. What was going on? I noticed that for the "slow" MPI library, the minimum benchmark run time was reasonable, but the maximum run time was a lot more. After reading through some messages in the e-mail archives of that MPI library, I found out that the library delays some setup costs until execution time. That means the first few calls of MPI collectives might be slow, but overall they will be fast. The other MPI library doesn't seem to be doing this.
I was reminded of some troubles that the MathWorks had with integrating FFTW into Matlab, a few years ago. FFTW does automatic performance tuning as a function of problem size. The tuning phase takes time, but once it's finished, running any problem of that size will likely be much faster. However, users didn't realize this. They ran an FFT once, noticed that it was much slower than before (because FFTW was busy doing its thing), and then complained. One can imagine all sorts of fixes for this, but the point is that something happened to performance which was not adequately communicated to users.
Now, MPI users are likely more sophisticated programmers than most Matlab users. They likely have heard of automatic run-time performance tuning, and maybe even use it themselves. I'm fairly familiar with such things and have been programming in MPI for a while. Yet, this phenomenon still puzzled me until I poked around through the e-mail list archives. What I was missing was an obvious pointer to documentation describing performance phenomena such as this one. The performance of this MPI library is not transparent. This is not a bad thing in this case, because the opaqueness is hiding a performance tuning process that will make my code faster overall. However, I want to know this right away, so I don't have to think about what I could be doing wrong.
See, I'm not just calling MPI directly. I've wrapped up some other wrapper over MPI. I didn't know if my wrapper was broken, the wrapper underneath was broken, the MPI library was broken, the cluster job execution software was broken, or if nothing was wrong and this was just expected behavior. That's why I want to know right away whether I should expect significantly nonuniform performance behavior from some library.
08 September 2010
CPUs and GPUs, part 3: memory and programming models
My last post got a good comment which pointed out that I hadn't been covering another important class of programming models: annotations. C and C++ programmers would call these "pragmas," since their compilers understand lines starting with "#pragma" as directives to the compiler that affect one or more lines of code that follow. They would likely be familiar with OpenMP, which uses annotations to tell the compiler to parallelize certain loops or fragments of code. Other programming languages, like Common Lisp, use annotations as hints for optimization -- in particular, as type hints.
The commenter pointed out that the Hybrid Multicore Parallel Programming (HMPP) system of annotations, developed by CAPS and supported by the PathScale ENZO compiler suit, offers another programming model for heterogeneous node computing. You write a kernel in a subset of C(++) or Fortran, annotated with operations that the compute device (e.g., GPU) understands (like when to transfer data to / from the host CPU, or when to synchronize device threads), and also label the kernel with an annotation to identify it for deployment to the device. I was just introduced to HMPP last night, so please take a look at the comment to read the manual for yourself :-) However, a brief read of the manual reminds me of Copperhead, a Python annotation system (based on decorators) that accomplishes similar goals.
Annotations require the programmer to (1) know the annotation language as well as the main programming language being used, (2) understand what subset of the programming language can be sent to the device, and (3) understand when and how the annotation changes program semantics. Compilers should help with #2, but may not; if given invalid code for the device, they may generate incorrect device code, which may crash or silently give the wrong answer. This is true for other annotation systems as well. OpenMP, for example, doesn't help you avoid race conditions if multiple threads modify shared state. You have to identify and annotate critical sections and / or atomic updates yourself. #3 can also be a challenge, especially if annotations are sprinkled throughout the code.
Before talking about HMPP, I'll first use OpenMP as an example of the difficulties of popular annotation-based approaches. Many OpenMP performance problems and code-writing challenges relate to the memory model. You've taken a sequential code, where memory is laid out for sequential operation, and introduced parallelism to the computation but not the memory layout. On nonuniform memory access (NUMA) hardware, performance tanks, because the data lives in the wrong place for most of the threads. At least modern OpenMP implementations try not to let threads wander around to different NUMA regions from where they started, which means that the operating system could migrate pages closer to the cores that are working on them. That suggests a "preconditioning" approach, where you first use OpenMP parallelism to read through the whole data set (thereby letting the operating system move data to where it needs to be), and then use OpenMP to do the computations you wanted to do in the first place. Operating systems in common use don't do that yet! Another approach is just syntactic sugar for explicitly threaded programming: each thread allocates its own data, and then you program in something like a distributed-memory model. (It would be easier and less buggy just to run MPI with one process per core.)
Both of these OpenMP approaches are awkward. The problem is that OpenMP separates computation and memory layout. The annotations only govern computation; they don't say where the data should live in memory. That often means throwing away useful information, which many OpenMP programmers know, about which threads will operate on which data. If you know how to allocate an array, you likely know how you would divide up simple parallel for loops or reductions on the array among the CPU cores. If you know that, you almost certainly know how to lay out the array in the various NUMA regions.
Heterogeneous programming has its own memory layout problem: data either lives on the host CPU or the device. Those are two separate memory spaces. Even if they weren't separate, it still might be expensive to compute on remote data, or even to copy data from one space to another. This is exactly the same problem as on multicore CPUs with NUMA. When I look at HMPP, I see annotations for when to perform data transfers between the host and the device, and even ways to overlap data transfers with useful computation. However, I don't see annotations for where to allocate the data. The HMPP model assumes (at least from what I can tell from a brief look -- please correct me if I'm wrong) that data always lives on the host, is always transferred to the device for computations, and is always transferred back to the host afterwards. This is like allocating all your data in one NUMA region and operating on it using CPU cores in all the NUMA regions.
One could augment the annotation approach to include a memory layout model, but it seems like this is just making the language more baroque. Plenty of programming languages reify data layout without the complication of additional annotations. Co-Array Fortran is a good example. This need not be done in the language; libraries can also hide the details of layout and allocation. Kokkos in Trilinos does this for GPUs, for example. CUDA and OpenCL are languages, but their special routines for allocating memory on the device are just plain C.
I like the Array Building Blocks (and Kokkos) approach, because it also gives a good model for interactions between legacy code (with legacy memory layout) and new parallel code (that expects or should have data in a particular nonlegacy layout). Kokkos does this with smart pointers and copy routines: you aren't allowed to call parallel kernels on legacy memory. Array Building Blocks does this with "bind", which links up legacy memory and compute buffers. The general idea is that "you don't get the pointer": memory on which heterogeneous or node-level parallel computation is done, is encapsulated and organized optimally by the language or library system. This is also a great way to unify different kinds of nodes with different memory space(s), which is why Kokkos does it. Pretty much any kind of node has parallel computing resources, associated to one or more memory spaces, for which data transfer costs are high and data layout matters for performance.
The commenter pointed out that the Hybrid Multicore Parallel Programming (HMPP) system of annotations, developed by CAPS and supported by the PathScale ENZO compiler suit, offers another programming model for heterogeneous node computing. You write a kernel in a subset of C(++) or Fortran, annotated with operations that the compute device (e.g., GPU) understands (like when to transfer data to / from the host CPU, or when to synchronize device threads), and also label the kernel with an annotation to identify it for deployment to the device. I was just introduced to HMPP last night, so please take a look at the comment to read the manual for yourself :-) However, a brief read of the manual reminds me of Copperhead, a Python annotation system (based on decorators) that accomplishes similar goals.
Annotations require the programmer to (1) know the annotation language as well as the main programming language being used, (2) understand what subset of the programming language can be sent to the device, and (3) understand when and how the annotation changes program semantics. Compilers should help with #2, but may not; if given invalid code for the device, they may generate incorrect device code, which may crash or silently give the wrong answer. This is true for other annotation systems as well. OpenMP, for example, doesn't help you avoid race conditions if multiple threads modify shared state. You have to identify and annotate critical sections and / or atomic updates yourself. #3 can also be a challenge, especially if annotations are sprinkled throughout the code.
Before talking about HMPP, I'll first use OpenMP as an example of the difficulties of popular annotation-based approaches. Many OpenMP performance problems and code-writing challenges relate to the memory model. You've taken a sequential code, where memory is laid out for sequential operation, and introduced parallelism to the computation but not the memory layout. On nonuniform memory access (NUMA) hardware, performance tanks, because the data lives in the wrong place for most of the threads. At least modern OpenMP implementations try not to let threads wander around to different NUMA regions from where they started, which means that the operating system could migrate pages closer to the cores that are working on them. That suggests a "preconditioning" approach, where you first use OpenMP parallelism to read through the whole data set (thereby letting the operating system move data to where it needs to be), and then use OpenMP to do the computations you wanted to do in the first place. Operating systems in common use don't do that yet! Another approach is just syntactic sugar for explicitly threaded programming: each thread allocates its own data, and then you program in something like a distributed-memory model. (It would be easier and less buggy just to run MPI with one process per core.)
Both of these OpenMP approaches are awkward. The problem is that OpenMP separates computation and memory layout. The annotations only govern computation; they don't say where the data should live in memory. That often means throwing away useful information, which many OpenMP programmers know, about which threads will operate on which data. If you know how to allocate an array, you likely know how you would divide up simple parallel for loops or reductions on the array among the CPU cores. If you know that, you almost certainly know how to lay out the array in the various NUMA regions.
Heterogeneous programming has its own memory layout problem: data either lives on the host CPU or the device. Those are two separate memory spaces. Even if they weren't separate, it still might be expensive to compute on remote data, or even to copy data from one space to another. This is exactly the same problem as on multicore CPUs with NUMA. When I look at HMPP, I see annotations for when to perform data transfers between the host and the device, and even ways to overlap data transfers with useful computation. However, I don't see annotations for where to allocate the data. The HMPP model assumes (at least from what I can tell from a brief look -- please correct me if I'm wrong) that data always lives on the host, is always transferred to the device for computations, and is always transferred back to the host afterwards. This is like allocating all your data in one NUMA region and operating on it using CPU cores in all the NUMA regions.
One could augment the annotation approach to include a memory layout model, but it seems like this is just making the language more baroque. Plenty of programming languages reify data layout without the complication of additional annotations. Co-Array Fortran is a good example. This need not be done in the language; libraries can also hide the details of layout and allocation. Kokkos in Trilinos does this for GPUs, for example. CUDA and OpenCL are languages, but their special routines for allocating memory on the device are just plain C.
I like the Array Building Blocks (and Kokkos) approach, because it also gives a good model for interactions between legacy code (with legacy memory layout) and new parallel code (that expects or should have data in a particular nonlegacy layout). Kokkos does this with smart pointers and copy routines: you aren't allowed to call parallel kernels on legacy memory. Array Building Blocks does this with "bind", which links up legacy memory and compute buffers. The general idea is that "you don't get the pointer": memory on which heterogeneous or node-level parallel computation is done, is encapsulated and organized optimally by the language or library system. This is also a great way to unify different kinds of nodes with different memory space(s), which is why Kokkos does it. Pretty much any kind of node has parallel computing resources, associated to one or more memory spaces, for which data transfer costs are high and data layout matters for performance.
05 September 2010
GPUs and CPUs, part 2: programming models
After finishing last night's post, I was thinking about CPU and GPU programming models -- in particular, about short vector instructions. On CPUs these are called SIMD (single instruction, multiple data) instructions. NVIDIA markets them instead as "SIMT" (single instruction, multiple thread), not because of differences in the hardware implementation, but because the programming model is different. The "CUDA C Programming Guide" version 3.1 says: "... SIMD vector organizations expose the SIMD width to the software, whereas SIMT instructions specify the execution and branching behavior of a single thread." Here, "thread" corresponds to an element in an array, operated on by a vector instruction.
CPU vendors have historically done a poor job presenting programmers with a usable interface to SIMD instructions. Compilers are supposed to vectorize code automatically, so that you write ordinary scalar loops and vectorized code comes out. However, my experience a few years ago writing SIMD code is that compilers have a hard time doing this for anything but the simplest code. Often, they will use SIMD instructions, but only operate on one element in the short vector, just because the vector pipes are faster than the scalar ones (even if the extra element(s) in the vector are wasted). The reason they struggle at vectorization -- more so than compilers for traditional vector machines like the Crays -- is because the SIMD instructions are more restrictive than traditional vector instructions. They demand contiguous memory, aligned on wide boundaries (like 128 bits). If the compiler can't prove that your memory accesses have this nice friendly pattern, it won't vectorize your code. SIMD branching is expensive (it means you lose the vector parallelism), so compilers might refuse to vectorize branchy code.
Note that GPU hardware has the same properties! They also like contiguous memory accesses and nonbranchy code. The difference is the programming model: "thread" means "array element" and consecutive threads are contiguous, whether you like it or not. The model encourages coding in a SIMD-friendly way. In contrast, CPU SIMD instructions didn't historically have much of a programming model at all, other than trusting in the compiler. The fact that compilers make the instructions directly available via "intrinsics" (only one step removed from inline assembler code) already indicates that coders had no other satisfactory interface to these instructions, yet didn't trust the compiler to vectorize. Experts like Sam Williams, a master performance tuner of computational kernels like sparse matrix-vector multiply and stencils, used the intrinsics to vectorize. This made their codes dependent on the particular instruction set, as well as on details of the instruction set implementation. (For example, older AMD CPUs used a "half-piped" implementation of SIMD instructions on 64-bit floating-point words. This meant that the implementation wasn't parallel, even though the instructions were. Using the x87 scalar instructions instead offered comparable performance, was easier to program, and even offered better accuracy for certain computations, since temporaries are stored with extra precision.) Using SIMD instructions complicated their code in other ways as well. For example, allocated memory has to be aligned to a particular size, such as 128 bits. That bit value depends on the SIMD vector width, which again decreases portability. SIMD instruction widths increase over time, and continue to do so. Furthermore, these codes are brittle, because feeding nonaligned memory to SIMD instructions often results in errors that can crash one's program.
CPU vendors are finally starting to think about programming models that make it easier to exploit SIMD instructions. OpenCL, while as low-level a model as CUDA, also lets programmers reason about vector instructions in a less hardware-dependent way. One of the most promising programming models is Intel's Array Building Blocks, featured in a recent Dr. Dobb's Journal article. I'm excited about Array Building Blocks for the following reasons:
CPU vendors have historically done a poor job presenting programmers with a usable interface to SIMD instructions. Compilers are supposed to vectorize code automatically, so that you write ordinary scalar loops and vectorized code comes out. However, my experience a few years ago writing SIMD code is that compilers have a hard time doing this for anything but the simplest code. Often, they will use SIMD instructions, but only operate on one element in the short vector, just because the vector pipes are faster than the scalar ones (even if the extra element(s) in the vector are wasted). The reason they struggle at vectorization -- more so than compilers for traditional vector machines like the Crays -- is because the SIMD instructions are more restrictive than traditional vector instructions. They demand contiguous memory, aligned on wide boundaries (like 128 bits). If the compiler can't prove that your memory accesses have this nice friendly pattern, it won't vectorize your code. SIMD branching is expensive (it means you lose the vector parallelism), so compilers might refuse to vectorize branchy code.
Note that GPU hardware has the same properties! They also like contiguous memory accesses and nonbranchy code. The difference is the programming model: "thread" means "array element" and consecutive threads are contiguous, whether you like it or not. The model encourages coding in a SIMD-friendly way. In contrast, CPU SIMD instructions didn't historically have much of a programming model at all, other than trusting in the compiler. The fact that compilers make the instructions directly available via "intrinsics" (only one step removed from inline assembler code) already indicates that coders had no other satisfactory interface to these instructions, yet didn't trust the compiler to vectorize. Experts like Sam Williams, a master performance tuner of computational kernels like sparse matrix-vector multiply and stencils, used the intrinsics to vectorize. This made their codes dependent on the particular instruction set, as well as on details of the instruction set implementation. (For example, older AMD CPUs used a "half-piped" implementation of SIMD instructions on 64-bit floating-point words. This meant that the implementation wasn't parallel, even though the instructions were. Using the x87 scalar instructions instead offered comparable performance, was easier to program, and even offered better accuracy for certain computations, since temporaries are stored with extra precision.) Using SIMD instructions complicated their code in other ways as well. For example, allocated memory has to be aligned to a particular size, such as 128 bits. That bit value depends on the SIMD vector width, which again decreases portability. SIMD instruction widths increase over time, and continue to do so. Furthermore, these codes are brittle, because feeding nonaligned memory to SIMD instructions often results in errors that can crash one's program.
CPU vendors are finally starting to think about programming models that make it easier to exploit SIMD instructions. OpenCL, while as low-level a model as CUDA, also lets programmers reason about vector instructions in a less hardware-dependent way. One of the most promising programming models is Intel's Array Building Blocks, featured in a recent Dr. Dobb's Journal article. I'm excited about Array Building Blocks for the following reasons:
- It includes a memory model with a segregated memory space. This can cover all kinds of complicated hardware details (like alignment and NUMA affinity). It's also validation for the work being done in libraries like Trilinos' Kokkos to hide array memory management from the programmer ("you don't get a pointer"), thus freeing the library to place and manage memory as it sees fit. All of this will help future-proof code from the details of allocation, and also make code safer (you don't get the pointer, so you can't give nonaligned memory to instructions that want aligned memory, or give GPU device memory to a CPU host routine).
- It's a programming language (embedded in C++, with a run time compilation model) -- which will expose mainstream programmers to the ideas of embedded special-purpose programming languages, and run time code generation and compilation. Scientific computing folks tend to be conservative about generating and compiling code at run time, in part because the systems on which they have to run often don't support it or only support it weakly. If Array Building Blocks gives the promised performance benefits and helps "future-proof" code, HPC system vendors will be driven to support these features.
- Its parallel operators are deterministic; programmers won't have to reason about shared-memory nightmares like race conditions or memory models. (Even writing a shared-memory barrier is a challenging task, and the best-performing barrier implementation depends on the application as well as the definition of "performance" (speed or energy?).)
Subscribe to:
Posts (Atom)
