B.1 Blitz++

The Blitz++ project brings high-performance numerical computing to C++. In some respects, it is what valarray<> should have been. Blitz++ has powerful array and matrix classes, operators, and functions, with strides, subarrays, and so on. The package is written to minimize the number of unnecessary temporary objects and take advantage of compile-time computation (via template metaprogramming) whenever possible.

One of the key optimizations is that arithmetic operators and mathematical functions involving Blitz++ arrays do not compute values immediately. Instead, they return expression objects. When an expression object is assigned to an array, the expression is computed, storing the results directly in the target of the assignment, without the need for large temporary arrays.

Example B-1 shows a program that demonstrates a few of the features of Blitz++.

Example B-1. Working with matrices

#include <iostream>
#include <ostream>
#include "blitz/array.h"

// Blitz formats output in a manner that is best suited for subsequent input by a
// program, not for reading by a human. The print(  ) function uses a format that
// is slightly better for human readers. (Further improvement is left as an
// exercise.)
template<typename T>
void print(const blitz::Array<T, 3>& a)
{
  std::cout << a.extent(0) << " x " << a.extent(1) <<
              " x " << a.extent(2) << ":\n[";
  for (size_t i = a.lbound(0); i <= a.ubound(0); ++i)
  {
    // Blitz can print a 2-D array well enough, so extract each 2-D plane and
    // print it. Note that plane shares storage with a, so even if a is large,
    // memory is not wasted by needlessly copying data.
    blitz::Array<T, 2>
      plane(a(i, blitz::Range::all(  ), blitz::Range::all(  )));
    std::cout << plane;
  }
  std::cout << "]\n";
}

int main(  )
{
  // Math with TinyVector and TinyMatrix uses template metaprogramming to produce
  // fast code, even for complicated operations, such as matrix multiplication.
  blitz::TinyMatrix<double,2,4> a;
  blitz::TinyMatrix<double,4,3> b;
  blitz::TinyMatrix<double,2,3> c;
  a = 1, 2, 3,  4,          // Set elements of a.
      0, 1, -2, -1;
  b = 3.14159;              // Set all elements of b.
  c = blitz::product(a, b); // Matrix multiplication
  std::cout << a << b << c;

  // Arrays can have more than two dimensions and offer more computing power.
  blitz::Array<double,3> d(2,3,4), e(2,3,4), f(2,3,4);
  // Set elements of d to values that depend on the indices.
  d = blitz::firstIndex(  ) + blitz::secondIndex(  ) +
      blitz::thirdIndex(  );
  // Set a subarray of d to 42.
  d(blitz::Range::all(  ),blitz::Range(1,2),blitz::Range(2,3))
    = 42;
  print(d);

  // Call sin for each element of d.
  e = sin(d);
  print(e);

  // If an element of e is negative, set corresponding element of f to 1; set
  // other elements to -1.
  f = blitz::where(e < 0, 1.0, -1.0);
  print(f);

  // Elementwise multiplication
  e *= f;
  print(e);
  // Add all elements of f.
  std::cout << blitz::sum(f) << '\n';
}