Blitz++: Fast, Accurate Numerical Computing in C++

The C++ programming language offers many useful features for tackling complex scientific computing problems. Inheritance, polymorphism, generic programming, and operator overloading are of primary importance. Unfortunately, most of these advanced features carry a hefty performance cost. Until recently, C++ performance lagged behind Fortran's performance by anywhere from 20 percent to orders of magnitude. As a result, the adoption of C++ for scientific computing has been slower than it should be, given the platform's popularity.

So, how can you turbo charge C++ so that you can have your advanced language features but lose the poor performance? This is precisely the goal of the Blitz++ project: to develop techniques that will enable C++ to rival and even surpass the speed of Fortran for numerical computing, while preserving an object-oriented interface. The Blitz++ Numerical Library is being constructed as the solution under dual licenses—GNU GPL and Blitz++ Artistic License.

Blitz++ provides a single array class called Array<T_numtype, N_rank>. This class provides a dynamically allocated N-dimensional array with reference counting, arbitrary storage ordering, subarrays and slicing, flexible expression handling, and many other useful features.

Can C++ Compete with Fortran?

Recent benchmarks show C++ encroaching steadily on Fortran's high-performance monopoly. For some benchmarks, C++ is even faster than Fortran! No, this is not a result of better optimizing compilers, preprocessors, or even language extensions. It's achieved simply through the use of template techniques. Through the magic of template metaprogramming, optimizations such as loop fusion, unrolling, tiling, and algorithm specialization can be performed automatically at compile time. Another successful package that uses this technique is, of course, the Boost C++ Libraries.

Another goal of Blitz++ is to extend the conventional dense array model (in other words, not sparse) with new and useful features. Some examples of such extensions are flexible storage formats, tensor notation, and index placeholders. Blitz++ is more than just cool arrays; you also get matrix math, random number generation, and a variety of supporting functionality.

Blitz++ Platform Support

As is true of any template metaprogramming tool, Blitz++ is only as good as the ISO/ANSI C++ compliance of the compiler you choose. Specifically, your compiler of choice must have rock-solid member templates, partial specialization, and partial ordering of templates. The good news is that Visual Studio .NET 2003 and Visual Studio 2005 are supported directly by way of project files. Additionally, a whole raft of compilers are supported, including Intel C++ for Windows, gcc, IBM XL C/C++ for AIX, and Sun Studio 10 for Solaris. Older versions of Borland C++ and Visual Studio C++ (6.0) definitely will not work. Of course, Blitz++ comes with a complete self-test suite so you can be confident that your OS and compiler are 100 percent compliant.

Getting Started

To use the Array class, you must include the header <blitz/array.h> and use the namespace blitz:

#include <blitz/array.h>
using namespace blitz;

The array class, called Array<T_numtype, N_rank>, has two templated parameters:

  • T_numtype is the numeric type to be stored in the array. T_numtype can be an integral type (scalars, including int and long double), complex type, or any user-defined type with appropriate numeric semantics.
  • N_rank is the dimensionality of the array (in other words, three for a three-dimensional array).
Array<int,1> x;                 // A one-dimensional array of int
Array<double,2> y;              // A two-dimensional array of double
Array<complex<float>, 12> z;    // A twelve-dimensional array of
// complex<float>

When no constructor arguments are provided, the array is empty and no memory is allocated. To create an array that contains some data, simply provide the size of the array as constructor arguments:

Array<double,2> y(4,4);    // A 4x4 array of double

The contents of a newly created array are garbage. To initialize the array to all zeroes, you can write the following:

y = 0;

Putting it all together, you can write a "hello world"-type program to add two 2x2 matrices together:

#include <blitz/array.h>
using namespace blitz;

int main() {
   Array<float,2> A(2,2), B(2,2), C(2,2);
   A = 1, 0, 2, 2;
   B = 0, 0, 7, 0;

   C = A + B;

   cout  << "C = " << C << endl;
   return 0;
}

and the output:

C = 2 x 2
       [        1            0]
                9             ]

The Array constructor can be manipulated in a lot of ways to achieve a variety of results: computation during initialization, setting the extent (length), setting the base of an array to be something other than zero, and building an array as the reference of another array. Here's a cool example of creating an array "B" that contains the square roots of elements in "A":

Array<float,2> A(N,N);
Array<float,2> B(sqrt(A));

Blitz++: Fast, Accurate Numerical Computing in C++

A Complete Example: Integration Techniques

Now, take a look at three different integration methods: naïve summation, two-point trapezoid, and a three-point Simpson's rule. The following is the expression for the example:

[Blitz1.png]

 1 #include <blitz/vector.h>
 2 using namespace blitz;
 3
 4 int main()
 5 {
 6    // Integrate the expression
 7
 8    // to estimate the error function erf(x) on the interval [0,1].
 9    //
10    // Three methods are compared:
11    //    1. Naive summation
12    //    2. Two-point trapezoid
13    //    3. Three-point Simpson's
14
15    const int numSamples = 1024;        // Number of samples
16    double dt = 1.0 / numSamples;       // Distance between samples
17    Range R(0, numSamples - 1);         // Index set for sample
                                          // vectors
18    Range displayRange(0, 900, 100);    // Indices to output to
                                          // screen
19
20    // For comparison purposes, use the built-in erf(x) function.
21    Vector<double> z = erf(R * dt);
22    cout << "  Built-in: " << z(displayRange+1) << endl;
23
24
25    // Naive summation
26    Vector<double> z1 = accumulate(M_2_SQRTPI *
                          exp(-sqr(R * dt)) * dt);
27    cout << "   Naive: " << z1(displayRange+1) << endl;
28
29    // Calculate the rms error
30    double error1 = sqrt(mean(sqr(z1 - z)));
31    cout << "RMS Error: " << error1 << endl;
32
33
34    // For the following rules, it's easier to work with a
35    // sampled integrand.
36    Vector<double> a = M_2_SQRTPI * exp(-sqr(R * dt));
37
38
39    // Two-point trapezoid
40    Range J(1, numSamples-1);
41    Vector<double> z2 = accumulate(dt / 2.0 * (a(J) + a(J-1)));
42
43    cout << " Trapezoid: " << z2(displayRange) << endl;
44    cout << "RMS Error: "  << sqrt(mean(sqr(z2-z(J)))) << endl;
45
46
47    // Three-point Simpson's (parabolic)
48    Range I(1, numSamples-2);
49    Vector<double> z3 = accumulate(dt / 6.0 * (a(I-1) + 4 *
                          a(I) + a(I+1)));

50
51    cout << " Parabolic: " << z3(displayRange) << endl;
52    cout << " RMS Error: " << sqrt(mean(sqr(z3 - z(I)))) << endl;
53
54
55    return 0;
56 }
57

Lines 26, 41, and 49 use the accumulate functor. A functor describes an action that should be taken on each member of a list. The action usually will vary according to the list member's type.

Parallel Computation

With multi-core CPUs on the horizon, everyone is thinking about how to exploit the extra cores by adding computation threads. Although Blitz++ can be used for parallel computing, it was not designed primarily for that. For this reason, you may want to investigate some other available libraries such as POOMA before choosing to implement a parallel code using Blitz++. Using threads on Blitz++ may require a library recompile, depending on platform.

In threadsafe mode, Blitz++ array reference counts are safeguarded by a mutex. By default, pthread mutexes are used. Blitz++ does not do locking for every array element access because it would result in an unacceptable performance hit. Responsibility for appropriate synchronization is squarely on your shoulders.

Tough Numeric Problems? No Problem

You've just barely scratched the surface of what Blitz++ can do for you. Some of the topics this article couldn't cover include sophisticated array expressions, stencils, slicing, and indirection. Give it a try because Blitz++ is ready, willing, and able to handle your toughest numeric problems.

Book Recommendation

Numerical Recipes in C++: The Art of Scientific Computing, 2nd Ed.

by Press, Tukolsky, Vetterling, and Flannery

ISBN: 0521750334

The Numerical Recipes series is a classic among numerical algorithm books. Thoroughly self-contained, it proceeds from mathematical and theoretical considerations to actual, practical computer routines. In addition to C++, Numerical Recipes books are also available for C, Fortran 77, Fortran 90, and Pascal.

About the Author

Victor Volkman has been writing for C/C++ Users Journal and other programming journals since the late 1980s. He is a graduate of Michigan Tech and a faculty advisor board member for Washtenaw Community College CIS department. Volkman is the editor of numerous books, including C/C++ Treasure Chest and is the owner of Loving Healing Press. He can help you in your quest for open source tools and libraries, just drop an e-mail to sysop@HAL9K.com.



About the Author

Victor Volkman

Victor Volkman has been writing for C/C++ Users Journal and other programming journals since the late 1980s. He is a graduate of Michigan Tech and a faculty advisor board member for Washtenaw Community College CIS department. Volkman is the editor of numerous books, including C/C++ Treasure Chest and is the owner of Loving Healing Press. He can help you in your quest for open source tools and libraries, just drop an e-mail to sysop@HAL9K.com.

Comments

  • There are no comments yet. Be the first to comment!

Leave a Comment
  • Your email address will not be published. All fields are required.

Top White Papers and Webcasts

  • The explosion in mobile devices and applications has generated a great deal of interest in APIs. Today's businesses are under increased pressure to make it easy to build apps, supply tools to help developers work more quickly, and deploy operational analytics so they can track users, developers, application performance, and more. Apigee Edge provides comprehensive API delivery tools and both operational and business-level analytics in an integrated platform. It is available as on-premise software or through …

  • Live Event Date: September 17, 2014 @ 1:00 p.m. ET / 10:00 a.m. PT Another day, another end-of-support deadline. You've heard enough about the hazards of not migrating to Windows Server 2008 or 2012. What you may not know is that there's plenty in it for you and your business, like increased automation and performance, time-saving technical features, and a lower total cost of ownership. Check out this upcoming eSeminar and join Rich Holmes, Pomeroy's practice director of virtualization, as he discusses the …

Most Popular Programming Stories

More for Developers

Latest Developer Headlines

RSS Feeds