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.
Showing posts with label programming. Show all posts
Showing posts with label programming. Show all posts
04 October 2010
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.
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?).)
25 August 2010
"Archetype by prototype": a suggestion for making compile-time polymorphism more usable
A recent post on interface matching in ParaSail (a new parallel programming language under development) inspired me to think about compile-time polymorphism in languages like C++. The author of the article observed that "C++ templates define essentially no 'contract' so there is no easy way to find out what sort of type will be acceptable without trying the instantiation."
This stirred me to think about what C++ programmers can (and do) do to get around this. A standard technique is to use "concept checks." These test the required syntax by, well, trying the instantiation ;-) in a way that makes it easy to catch syntax errors. They can be tedious, because doing them right requires two things:
1. Concept checking code, that verifies the interface statically and as simply as possible (so that compiler error messages are easier to read)
2. An "archetype class": like a Java interface, except with (skeletal) implementations of the required interface
In the above, i'm using the terminology of the Boost Concept Check Library. #1 alone adds to the code maintenance burden, since as the "contract" evolves, you have to synchronize the checks with how you are actually using the object. #2 alone doesn't ensure synchronization of your idea of the contract, with the classes you think are implementing it. Nevertheless, there are still four separate bodies of code to synchronize: the archetype class, the concept checking code, the application class(es), and the application code that invokes these class(es). The most straightforward (but not easiest) way to avoid this burden would be to
1. Write an archetype class, and make the compiler check the concept, using application code
I call this approach "declared archetypes." It calls for a syntax for declaring that a particular concrete data type implements an archetype's interface. This approach has implications for avoiding error-prone tedium. A common case for templates is generic code for different kinds of numbers. I've written a lot of code templated on "Scalar", which could be anything representing a real or complex number, and Ordinal, which represents an array index type. It would be tedious to declare an archetype supporting all the different kinds of things one might like to do with a Scalar, for example. It would at least have to include all of arithmetic, and all kinds of transcendental functions (such as trigonometric functions and logarithms). That would be an awful lot of effort, especially since most people are only going to instantiate Scalar with "float" or "double". (A few unlucky folks will make Scalar a complex number and discover all the places where the supposedly generic code was using less-than on Scalar!)
To avoid this tedium, the "concepts" proposal considered for addition to the C++ standard (but rejected) offered a shorthand for "aggregate archetypes." You could say that Scalar is "FloatingPointLike" or that Ordinal is "SignedIntegralLike", and that would hopefully bring in the syntax that you want. You could also "inherit" from archetypes to create your own.
The trouble with this approach is that somebody has to write a huge library of arbitrarily named predefined archetypes. Huge, because pretty much anything you might want to do with a "plain old data" type has to have its own concept. Arbitrarily named, because somebody has to decide what "ArithmeticLike" means. (Does it mean integers or real numbers?) This path seems to call for abstract algebra, where you have semigroups and rings and fields of a certain characteristic, and all kinds of things that nonmathematicians don't want to understand. (I'm convinced the main reason why people don't use Haskell more is because it's so hard to explain what a monad is, and why you want to know.)
This is overkill because in most cases, programmers can accomplish their work without so much formality. The ParaSail blog post alludes to the reason why: "... when ad hoc matching is used, there is a presumption that the module formal (or target type) is very stable and new operations are not expected to be added to it." The typical use case of C++ templates (for example) is for things that "look like double" or "look like int." That suggests a second approach:
2. "Archetype by prototype" or "looks like type T"
If T is something simple like "double" or "int", then you save a lot of syntax and / or library writing. If T is complicated, this forces you to write an archetype class, which is probably a good thing if T is complicated! For the numerical linear algebra algorithms i write and maintain, this would help remind me whether an algorithm (that claims to be generic on a Scalar type) can handle complex numbers. (Complex arithmetic changes linear algebra in subtle ways that don't only have to do with syntax.)
This approach does not require special syntax for defining archetypes. However, it might still be nice to have such a syntax, and also to have a small library of archetype classes. I could see this being useful for iterators, for example.
"Archetype by prototype" would impose requirements "lazily," much like users of C++ template metaprogramming expect. The "actual archetype" would consist only of those operations with the concrete datatype that are actually written in code. For example, if you never ask for the cosine of a Scalar, you don't need to have an implementation of cosine that takes a Scalar argument. It would be nice to have some development environment support for "extracting" the "actual archetype." C++ compilers do some of this already in their error messages, if you have the patience to read them (after about a page of compiler backtrace, you get something like "no implementation of cos(Scalar) exists for Scalar = YourType").
In summary, i propose "archetype by prototype" as an alternative to "concepts" for making compile-time polymorphism more usable.
This stirred me to think about what C++ programmers can (and do) do to get around this. A standard technique is to use "concept checks." These test the required syntax by, well, trying the instantiation ;-) in a way that makes it easy to catch syntax errors. They can be tedious, because doing them right requires two things:
1. Concept checking code, that verifies the interface statically and as simply as possible (so that compiler error messages are easier to read)
2. An "archetype class": like a Java interface, except with (skeletal) implementations of the required interface
In the above, i'm using the terminology of the Boost Concept Check Library. #1 alone adds to the code maintenance burden, since as the "contract" evolves, you have to synchronize the checks with how you are actually using the object. #2 alone doesn't ensure synchronization of your idea of the contract, with the classes you think are implementing it. Nevertheless, there are still four separate bodies of code to synchronize: the archetype class, the concept checking code, the application class(es), and the application code that invokes these class(es). The most straightforward (but not easiest) way to avoid this burden would be to
1. Write an archetype class, and make the compiler check the concept, using application code
I call this approach "declared archetypes." It calls for a syntax for declaring that a particular concrete data type implements an archetype's interface. This approach has implications for avoiding error-prone tedium. A common case for templates is generic code for different kinds of numbers. I've written a lot of code templated on "Scalar", which could be anything representing a real or complex number, and Ordinal, which represents an array index type. It would be tedious to declare an archetype supporting all the different kinds of things one might like to do with a Scalar, for example. It would at least have to include all of arithmetic, and all kinds of transcendental functions (such as trigonometric functions and logarithms). That would be an awful lot of effort, especially since most people are only going to instantiate Scalar with "float" or "double". (A few unlucky folks will make Scalar a complex number and discover all the places where the supposedly generic code was using less-than on Scalar!)
To avoid this tedium, the "concepts" proposal considered for addition to the C++ standard (but rejected) offered a shorthand for "aggregate archetypes." You could say that Scalar is "FloatingPointLike" or that Ordinal is "SignedIntegralLike", and that would hopefully bring in the syntax that you want. You could also "inherit" from archetypes to create your own.
The trouble with this approach is that somebody has to write a huge library of arbitrarily named predefined archetypes. Huge, because pretty much anything you might want to do with a "plain old data" type has to have its own concept. Arbitrarily named, because somebody has to decide what "ArithmeticLike" means. (Does it mean integers or real numbers?) This path seems to call for abstract algebra, where you have semigroups and rings and fields of a certain characteristic, and all kinds of things that nonmathematicians don't want to understand. (I'm convinced the main reason why people don't use Haskell more is because it's so hard to explain what a monad is, and why you want to know.)
This is overkill because in most cases, programmers can accomplish their work without so much formality. The ParaSail blog post alludes to the reason why: "... when ad hoc matching is used, there is a presumption that the module formal (or target type) is very stable and new operations are not expected to be added to it." The typical use case of C++ templates (for example) is for things that "look like double" or "look like int." That suggests a second approach:
2. "Archetype by prototype" or "looks like type T"
If T is something simple like "double" or "int", then you save a lot of syntax and / or library writing. If T is complicated, this forces you to write an archetype class, which is probably a good thing if T is complicated! For the numerical linear algebra algorithms i write and maintain, this would help remind me whether an algorithm (that claims to be generic on a Scalar type) can handle complex numbers. (Complex arithmetic changes linear algebra in subtle ways that don't only have to do with syntax.)
This approach does not require special syntax for defining archetypes. However, it might still be nice to have such a syntax, and also to have a small library of archetype classes. I could see this being useful for iterators, for example.
"Archetype by prototype" would impose requirements "lazily," much like users of C++ template metaprogramming expect. The "actual archetype" would consist only of those operations with the concrete datatype that are actually written in code. For example, if you never ask for the cosine of a Scalar, you don't need to have an implementation of cosine that takes a Scalar argument. It would be nice to have some development environment support for "extracting" the "actual archetype." C++ compilers do some of this already in their error messages, if you have the patience to read them (after about a page of compiler backtrace, you get something like "no implementation of cos(Scalar) exists for Scalar = YourType").
In summary, i propose "archetype by prototype" as an alternative to "concepts" for making compile-time polymorphism more usable.
Subscribe to:
Posts (Atom)
