Blitz++ Benchmarks: Acoustic

This benchmark solves the 2-D acoustic wave equation using a finite difference approximation. The grid size is 650x650.

Snapshots from the simulation: 1. A pressure pulse emerges from a rectangular channel; 2. The pulse meets a bar of dense material; 3. The pulse travels faster through the dense bar.

Timing results

The Mflops rates are low because each iteration requires 6.5 Mb of data to be shipped through the CPU; this is far larger than the cache size. It's theoretically possible to get better performance through an iteration-space tiling, but this is beyond the scope of Blitz++ for now. Note that a conventional C++ array class library which used pairwise evaluation of expressions would require 39 Mb of data to go through the CPU each iteration; it would run roughly 6 times slower than the Fortran and Blitz++ versions!

Optimizations

The tuned versions do two important optimizations:

• Array relabelling: each time step requires copying the arrays: P1 <- P2, P2 <- P3. These copy operations almost double the amount of data which has to go through the CPU (depending on whether the loops are fused or not). Performance on this problem is limited by memory bandwidth, so copying arrays effectively halves the performance.
  
  Since Blitz++ arrays are really lightweight handles, arrays can be copied in constant time by swapping the handles rather than the actual array data (this is done by the cycleArrays(..) method). In Fortran, it's not possible to do this. You can get a similar effect by putting the stencilling code in a subroutine and doing three iterations at a time (see the source code for details).
  
  
• Tiling: A naive implementation of the benchmark only achieves a 30% L1 cache hit rate. This is because three scan lines of the array is 1024x3x4x4 = 48 kb, much larger than the L1 cache. To get a higher hit rate, it's necessary to used a tiled traversal, in which you divide the array into tiles of length (say) 128. Each tile is completely traversed before going on to the next. Three scan lines of the tile do fit comfortably into the L1 cache: 128x3x4x4 = 6 kb.
  
  At compile time, Blitz++ recognizes the expression as a 2D stencil, and automatically generates a tiled traversal. I wasn't able to get the XL Fortran 77 compiler to transform the naive implementation into a tiled traversal, so it had to be done manually. The XL Fortran 90 compiler also isn't able to generate a tiled algorithm automatically. Doing this optimization explicitly would ruin the high-level, abstract code and turn it into the Fortran 77 version. Better Fortran 90 compilers (such as Cray F90) are able to generate tilings automatically; on the T3E, the Fortran 90 version is roughly as fast as the Fortran 77 one. (I'm unable to provide Blitz++ results for the T3E because I don't yet have access to KAI C++ on the T3E).
  
  

The tuned Blitz++ and Fortran 77 versions do exactly the same optimizations. The tuned Fortran 90 version does as much optimization as possible while staying at the same level of abstraction as Blitz++. If you compare the tuned and untuned versions, you can see that the Blitz++ tuned version is a lot cleaner. The tuned Fortran versions lose their readability somewhat -- cycling arrays in Fortran isn't pretty.

XL Fortran is a state-of-the-art compiler. It's able to do many optimizations which KAI C++ doesn't, such as high-order loop transformations (interchanging, fusion, etc.), relaxation of IEEE arithmetic rules, profile-directed optimizations, inter-procedural optimization, and optimizing for specific cache structures. Blitz++ is able to compete with it by doing similar optimizations at the language level: it uses template techniques to generate optimized code.

Initial conditions setup