Array expressions in Blitz++ are implemented using the expression
templates technique. Unless otherwise noted, expression evaluation
will never generate temporaries or multiple loops; an expression such
as
will result in code similar to
Array<int,1> A, B, C, D; // ...
A = B + C + D;
for (int i=A.lbound(firstDim); i <= A.ubound(firstDim); ++i)
A[i] = B[i] + C[i] + D[i];
3.1: Expression evaluation order |
A commonly asked question about Blitz++ is what order it uses
to evaluate array expressions. For example, in code such
as
does the expression get evaluated at indices 1, 2, ..., 9 or
at 9, 8, ..., 1? This makes a big difference to the result:
in one case, the array will be shifted to the right by one
element; in the other case, most of the array elements will
be set to the value in Blitz always selects the traversal order it thinks will be fastest.
For 1D arrays, this means it will go from beginning to the end
of the array in memory (see notes below). For multidimensional arrays, it
will do one of two things:
Because the traversal order is not always predictable, it is
safest to put the result in a new array if you are doing a
stencil-style expression. Blitz guarantees this
will always work correctly. If you try to put the result in one
of the operands, you have to guess correctly which traversal order
blitz will choose. This is easy for the 1D case, but hard for
the multidimensional case.
Some special notes about 1D array traversals:
A(Range(2,10)) = A(Range(1,9))
A(1)
.
3.2: Expression operands |
3.3: Array operands |
3.4: Expression operators |
These binary operators are supported:
Note: operator These unary operators are supported:
The operators
All operators are applied elementwise.
+ - * / % > < >= <= == != && || ^ & |
<<
and >>
are reserved for use in input/output.
If you need a bit-shift operation on arrays, you may define one
yourself; see 3.10.
- ~ !
> < >= <= == != && || !
result in a bool-valued
expression.
3.5: Assignment operators |
3.6: Index placeholders |
Blitz++ provides objects called index placeholders which represent
array indices. They can be used directly in expressions.
There is a distinct index placeholder type associated with each
dimension of an array. The types are called This generates code which is similar to:
Here's an example which fills an array with a sampled sine wave:
If your destination array has rank greater than 1, you may use
multiple index placeholders:
Here's a plot of the resulting array:
You can use index placeholder expressions in up to 11 dimensions.
Here's a three dimensional example:
You can mix array operands and index placeholders:
For your convenience, there is a namespace within blitz
called So instead of declaring your own index placeholder objects,
you can just say
firstIndex
,
secondIndex
, thirdIndex
, ..., tenthIndex
, eleventhIndex
.
Here's an example of using an index placeholder:
Array<float,1> A(10);
firstIndex i;
A = i;
for (int i=0; i < A.length(); ++i)
A(i) = i;
Array<float,1> A(16);
firstIndex i;
A = sin(2 * M_PI * i / 16.);
// Fill a two-dimensional array with a radially
// symmetric, decaying sinusoid
// Create the array
int N = 64;
Array<float,2> F(N,N);
// Some parameters
float midpoint = (N-1)/2.;
int cycles = 3;
float omega = 2.0 * M_PI * cycles / double(N);
float tau = - 10.0 / N;
// Index placeholders
firstIndex i;
secondIndex j;
// Fill the array
F = cos(omega * sqrt(pow2(i-midpoint) + pow2(j-midpoint)))
* exp(tau * sqrt(pow2(i-midpoint) + pow2(j-midpoint)));
Figure 3: Array filled using an index placeholder expression
// Fill a three-dimensional array with a Gaussian function
Array<float,3> A(16,16,16);
firstIndex i;
secondIndex j;
thirdIndex k;
float midpoint = 15/2.;
float c = - 1/3.0;
A = exp(c * (sqr(i-midpoint) + sqr(j-midpoint)
+ sqr(k-midpoint)));
Array<int,1> A(5), B(5);
firstIndex i;
A = 0, 1, 1, 0, 2;
B = i * A; // Results in [ 0, 1, 2, 0, 8 ]
tensor
which declares all the index placeholders:
namespace blitz {
namespace tensor {
firstIndex i;
secondIndex j;
thirdIndex k;
...
eleventhIndex t;
}
}
3.7: Type promotion |
When operands of different numeric types are used in
an expression, the result gets promoted according to the
usual C-style type promotion. For example, the result of
adding an
The rules for type promotion of user-defined types (or
types from another library) are a bit complicated.
Here's how a pair of operand types are promoted:
If you wish to alter the default type promotion rules
above, you have two choices:
Note that you can do these specializations in your own header
files (you don't have to edit promote.h or ops.h).
There are some inconvenient aspects of C-style type promotion. For
example, when you divide two integers in C, the result gets truncated.
The same problem occurs when dividing two integer arrays in Blitz++:
The usual solution to this problem is to cast one of the operands to a floating
type. For this purpose, Blitz++ provides a function Array<int>
to an Arrray<float>
will be
promoted to float
. Generally, the result is promoted
to whichever type has greater precision.
Type promotion for user-defined types
complex<float>
, complex<double>
, and
complex<long double>
.
sizeof()
). The rationale is that
more storage space probably indicates more precision.
Manual casts
Array<int,1> A(4), B(4);
Array<float,1> C(4);
A = 1, 2, 3, 5;
B = 2, 2, 2, 7;
C = A / B; // Result: [ 0 1 1 0 ]
cast(expr,type)
which will cast the result of expr as type:
3.8: Single-argument math functions |
All of the functions described in this section are element-wise. For example, this code--
Array<float,2> A, B; // A = sin(B);
results in A(i,j) = sin(B(i,j))
for all (i,j).
These math functions are available on all platforms:
These functions are only available on platforms which provide the
IEEE Math library (libm.a) and/or System V Math Library (libmsaa.a).
Apparently not all platforms provide all of these functions, so
what you can use on your platform may be a subset of these.
If you choose to use one of these functions, be aware that you
may be limiting the portability of your code.
On some platforms, the preprocessor symbols
In the current version, Blitz++ divides these functions into two
groups: IEEE and System V. This distinction is probably artificial.
If one of the functions in a group is missing,
Blitz++ won't allow you to use any of them. You can see the
division of these functions in the files You may have to link with None of these functions are available for There may be better descriptions of these functions in your
system man pages.
atan2()
in section
3.9.
complex<T>
.
complex<T>
.
complex<T>
.
complex<T>
.
complex<T>
.
complex<T>
.
complex<T>
.
complex<T>
.
complex<T>
.
complex<T>
.
complex<T>
.
IEEE/System V math functions
_XOPEN_SOURCE
and/or _XOPEN_SOURCE_EXTENDED
need to be defined
to use these functions. These symbols can be enabled
by compiling with -DBZ_ENABLE_XOPEN_SOURCE
.
(In previous version of Blitz++, _XOPEN_SOURCE
and
_XOPEN_SOURCE_EXTENDED
were declared by default.
This was found to cause too many problems, so users
must manually enable them with -DBZ_ENABLE_XOPEN_SOURCE
.).
Blitz++/compiler/ieeemath.cpp
and Blitz++/compiler/sysvmath.cpp
. This arrangement is
unsatisfactory and will probably change in a future version.
-lm
and/or -lmsaa
to use these
functions.
complex<T>
.
rint()
rounds up or down or to the nearest
integer depends on the current floating-point rounding mode.
If you haven't altered the rounding mode, rint()
should be
equivalent to nearest()
. If rounding mode is set to round
towards +INF, rint()
is equivalent to ceil()
. If the
mode is round toward -INF, rint()
is equivalent to
floor()
. If the mode is round toward zero, rint()
is equivalent to trunc()
.
)
3.9: Two-argument math functions |
The math functions described in this section take two arguments. Most combinations of these types may be used as arguments:
float
, double
, long double
,
or complex<T>
These math functions are available on all platforms, and work for complex numbers.
complex<T>
.
blitz::
scope qualifier is needed to
disambiguate the ANSI C++ function template
polar(T,T)
. This qualifier will hopefully
disappear in a future version.
complex<T>
.
See the notes about IEEE/System V math functions in the previous section. None of these functions work for complex numbers. They will all cast their arguments to double precision.
nearest(x/y)
(the nearest integer to
x/y). The return value will lie in the range
[ -y/2, +y/2 ]. If y is zero or x is +INF or -INF,
NaNQ is returned.
3.10: Declaring your own math functions on arrays |
There are four macros which make it easy to turn
your own scalar functions into functions defined
on arrays. They are:
3.11: Tensor notation |
Blitz++ arrays support a tensor-like notation. Here's an example of
real-world tensor notation:
A is a rank 3 tensor (a three dimensional array), B is a rank 2 tensor
(a two dimensional array), and C is a rank 1 tensor (a one dimensional
array). The above expression sets
A(i,j,k) = B(i,j) * C(k).
To implement this product using Blitz++, we'll need the arrays and some
index placeholders:
Here's the Blitz++ code which is equivalent to the tensor expression:
The index placeholder arguments tell an array how to map its dimensions
onto the dimensions of the destination array.
For example, here's some real-world tensor notation:
In Blitz++, this would be coded as:
This tensor expression can be visualized in the following way:
Here's an example which computes an outer product of two one-dimensional
arrays:
And the output:
Index placeholders can not be used on the left-hand side of an
expression. If you need to reorder the indices, you must do this
on the right-hand side.
In real-world tensor notation, repeated indices imply a contraction
(or summation). For example, this tensor expression computes
a matrix-matrix product:
The repeated k index is interpreted as meaning
In Blitz++, repeated indices do not imply contraction. If you
want to contract (sum along) an index, you must use the The Index placeholders can be used in any order in an expression.
This example computes a kronecker product of a pair of
two-dimensional arrays, and permutes the indices along the
way:
This is equivalent to the tensor notation
Tensor-like notation can be mixed with other array notations:
ijk ij k
A = B C
Array<float,3> A(4,4,4);
Array<float,2> B(4,4);
Array<float,1> C(4);
firstIndex i; // Alternately, could just say
secondIndex j; // using namespace blitz::tensor;
thirdIndex k;
A = B(i,j) * C(k);
ijk ij k jk i
C = A x - A y
using namespace blitz::tensor;
C = A(i,j) * x(k) - A(j,k) * y(i);
Figure 4: Examples of array indexing, subarrays, and slicing.
#include <blitz/array.h>
using namespace blitz;
int main()
{
Array<float,1> x(4), y(4);
Array<float,2> A(4,4);
x = 1, 2, 3, 4;
y = 1, 0, 0, 1;
firstIndex i;
secondIndex j;
A = x(i) * y(j);
cout << A << endl;
return 0;
}
4 x 4
1 0 0 1
2 0 0 2
3 0 0 3
4 0 0 4
ij ik kj
C = A B
c = sum of {a * b } over k
ij ik kj
sum()
function:
Array<float,2> A, B, C; // ...
firstIndex i;
secondIndex j;
thirdIndex k;
C = sum(A(i,k) * B(k,j), k);
sum()
function is an example of an array reduction,
described in the next section.
Array<float,2> A, B; // ...
Array<float,4> C; // ...
fourthIndex l;
C = A(l,j) * B(k,i);
ijkl lj ki
C = A B
Array<float,2> A, B; // ...
Array<double,4> C; // ...
C = cos(A(l,j)) * sin(B(k,i)) + 1./(i+j+k+l);
3.12: Array reductions |
Currently, Blitz++ arrays support two forms of reduction: (A
future version will likely support block reductions.)
3.13: Complete reductions |
Complete reductions transform an array (or array expression) into a scalar. Here are some examples:
Array<float,2> A(3,3); A = 0, 1, 2, 3, 4, 5, 6, 7, 8; cout << sum(A) << endl // 36 << min(A) << endl // 0 << count(A >= 4) << endl; // 5
Here are the available complete reductions:
Note: minIndex() and maxIndex() return TinyVectors, even when the rank of the array (or array expression) is 1.
Reductions can be combined with where
expressions (3.15)
to reduce over some part of an array. For example,
sum(where(A > 0, A, 0))
sums only the positive elements
in an array.
3.14: Partial Reductions |
Here's an example which computes the sum of each row of a two-dimensional
array:
The reduction Reductions have an important restriction: It is currently only
possible to reduce over the last dimension of an array or array
expression. Reducing a dimension other than the last would require
Blitz++ to reorder the dimensions to fill the hole left behind.
For example, in order for this reduction to work:
Blitz++ would have to remap the dimensions so that the third dimension
became the second. It's not currently smart enough to do this.
However, there is a simple workaround which solves some of the
problems created by this limitation: you can do the reordering
manually, prior to the reduction:
Writing The reductions To illustrate, here's an example:
The array This table shows what the result stored in Note: the odd numbers for first() are
The result of a reduction is an array expression, so reductions
can be used as operands in an array expression:
Note that this is not allowed:
You cannot reduce an array to zero dimensions!
Instead, use one of the global functions described in the
previous section.
Array<float,2> A; // ...
Array<float,1> rs; // ...
firstIndex i;
secondIndex j;
rs = sum(A, j);
sum()
takes two arguments:
Array<float,3> A; // ...
Array<float,2> B; // ...
secondIndex j;
// Reduce over dimension 2 of a 3-D array?
B = sum(A, j);
B = sum(A(i,k,j), k);
A(i,k,j)
interchanges the second and third dimensions,
permitting you to reduce over the second dimension.
Here's a list of the reduction operations currently supported:
tiny(int())
(INT_MIN) is returned.
huge(int())
(INT_MAX) is returned.
any()
, all()
, and first()
have short-circuit
semantics: the reduction will halt as soon as the answer is known. For
example, if you use any()
, scanning of the expression will stop as
soon as the first true value is encountered.
Array<int, 2> A(4,4);
A = 3, 8, 0, 1,
1, -1, 9, 3,
2, -5, -1, 1,
4, 3, 4, 2;
Array<float, 1> z;
firstIndex i;
secondIndex j;
z = sum(A(j,i), j);
z
now contains the sum of A
along each column:
[ 10 5 12 7 ]
z
would be if
sum()
were replaced with other reductions:
sum [ 10 5 12 7 ]
mean [ 2.5 1.25 3 1.75 ]
min [ 1 -5 -1 1 ]
minIndex [ 1 2 2 0 ]
max [ 4 8 9 3 ]
maxIndex [ 3 0 1 1 ]
first((A < 0), j) [ -2147483648 1 2 -2147483648 ]
product [ 24 120 0 6 ]
count((A(j,i) > 0), j) [ 4 2 2 4 ]
any(abs(A(j,i)) > 4, j) [ 0 1 1 0 ]
all(A(j,i) > 0, j) [ 1 0 0 1 ]
tiny(int())
i.e. the
smallest number representable by an int. The exact value is
machine-dependent.
Array<int,3> A;
Array<int,2> B;
Array<int,1> C; // ...
secondIndex j;
thirdIndex k;
B = sqrt(sum(sqr(A), k));
// Do two reductions in a row
C = sum(sum(A, k), j);
Array<int,2> A;
firstIndex i;
secondIndex j;
// Completely sum the array?
int result = sum(sum(A, j), i);
3.15: where statements |
Blitz++ provides the "where" function as an array expression
version of the "?:" operator. The syntax is:
Wherever
where(array-expr1, array-expr2, array-expr3)
array-expr1
is true, array-expr2
is returned.
Where array-expr1
is false, array-expr3
is returned.
For example, suppose we wanted to sum the squares of only
the positive elements of an array. This can be implemented
using a where function:
double posSquareSum = sum(where(A > 0, pow2(A), 0));