Next: Where expr, Previous: Math functions 2, Up: Array Expressions [Contents][Index]
There are four macros which make it easy to turn your own scalar functions into functions defined on arrays. They are:
BZ_DECLARE_FUNCTION(f) // 1 BZ_DECLARE_FUNCTION_RET(f,return_type) // 2 BZ_DECLARE_FUNCTION2(f) // 3 BZ_DECLARE_FUNCTION2_RET(f,return_type) // 4
Use version 1 when you have a function which takes one argument and returns a result of the same type. For example:
#include <blitz/array.h>
using namespace blitz;
double myFunction(double x)
{
return 1.0 / (1 + x);
}
BZ_DECLARE_FUNCTION(myFunction)
int main()
{
Array<double,2> A(4,4), B(4,4); // ...
B = myFunction(A);
}
Use version 2 when you have a one argument function whose return type is different than the argument type, such as
int g(double x);
Use version 3 for a function which takes two arguments and returns a result of the same type, such as:
double g(double x, double y);
Use version 4 for a function of two arguments which returns a different type, such as:
int g(double x, double y);
Blitz++ arrays support a tensor-like notation. Here’s an example of real-world tensor notation:
ijk ij k A = B C
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:
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;
Here’s the Blitz++ code which is equivalent to the tensor expression:
A = B(i,j) * C(k);
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:
ijk ij k jk i C = A x - A y
In Blitz++, this would be coded as:
using namespace blitz::tensor; C = A(i,j) * x(k) - A(j,k) * y(i);
This tensor expression can be visualized in the following way:
Here’s an example which computes an outer product of two one-dimensional arrays:
#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;
}
And the output:
(0,3) x (0,3) [ 1 0 0 1 2 0 0 2 3 0 0 3 4 0 0 4 ]
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:
ij ik kj C = A B
The repeated k index is interpreted as meaning
c = sum of (a * b ) over k ij ik kj
In Blitz++, repeated indices do not imply contraction. If you want
to contract (sum along) an index, you must use the sum() function:
Array<float,2> A, B, C; // ... firstIndex i; secondIndex j; thirdIndex k; C = sum(A(i,k) * B(k,j), k);
The sum() function is an example of an array reduction,
described in the next section.
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:
Array<float,2> A, B; // ... Array<float,4> C; // ... fourthIndex l; C = A(l,j) * B(k,i);
This is equivalent to the tensor notation
ijkl lj ki C = A B
Tensor-like notation can be mixed with other array notations:
Array<float,2> A, B; // ... Array<double,4> C; // ... C = cos(A(l,j)) * sin(B(k,i)) + 1./(i+j+k+l);
An important efficiency note about tensor-like notation: the right-hand side of an expression is completely evaluated for every element in the destination array. For example, in this code:
Array<float,1> x(4), y(4); Array<float,2> A(4,4): A = cos(x(i)) * sin(y(j));
The resulting implementation will look something like this:
for (int n=0; n < 4; ++n)
for (int m=0; m < 4; ++m)
A(n,m) = cos(x(n)) * sin(y(m));
The functions cos and sin will be invoked sixteen times each.
It’s possible that a good optimizing compiler could hoist the cos
evaluation out of the inner loop, but don’t hold your breath – there’s a
lot of complicated machinery behind the scenes to handle tensor notation,
and most optimizing compilers are easily confused. In a situation like the
above, you are probably best off manually creating temporaries for
cos(x) and sin(y) first.
Currently, Blitz++ arrays support two forms of reduction:
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:
sum()Summation (may be promoted to a higher-precision type)
product()Product
mean()Arithmetic mean (promoted to floating-point type if necessary)
min()Minimum value
max()Maximum value
minmax()Simultaneous minimum and maximum value (returns a value of type MinMaxValue<T>)
minIndex()Index of the minimum value (TinyVector<int,N_rank>)
maxIndex()Index of the maximum value (TinyVector<int,N_rank>)
count()Counts the number of times the expression is logical true (int)
any()True if the expression is true anywhere (bool)
all()True if the expression is true everywhere (bool)
Caution: minIndex() and maxIndex() return TinyVectors,
even when the rank of the array (or array expression) is 1.
Reductions can be combined with where expressions (Where expr)
to reduce over some part of an array. For example, sum(where(A > 0,
A, 0)) sums only the positive elements in an array.
Here’s an example which computes the sum of each row of a two-dimensional array:
Array<float,2> A; // ... Array<float,1> rs; // ... firstIndex i; secondIndex j; rs = sum(A, j);
The reduction sum() takes two arguments:
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:
Array<float,3> A; // ... Array<float,2> B; // ... secondIndex j; // Reduce over dimension 2 of a 3-D array? B = sum(A, j);
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:
B = sum(A(i,k,j), k);
Writing 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:
sum()Summation
product()Product
mean()Arithmetic mean (promoted to floating-point type if necessary)
min()Minimum value
max()Maximum value
minIndex()Index of the minimum value (int)
maxIndex()Index of the maximum value (int)
count()Counts the number of times the expression is logical true (int)
any()True if the expression is true anywhere (bool)
all()True if the expression is true everywhere (bool)
first()First index at which the expression is logical true (int); if the expression
is logical true nowhere, then tiny(int()) (INT_MIN) is returned.
last()Last index at which the expression is logical true (int); if the expression
is logical true nowhere, then huge(int()) (INT_MAX) is returned.
The reductions 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.
To illustrate, here’s an example:
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(4);
firstIndex i;
secondIndex j;
z = sum(A(j,i), j);
The array z now contains the sum of A along each column:
[ 10 5 12 7 ]
This table shows what the result stored in 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 ]
Note: the odd numbers for first() are tiny(int()) i.e. the smallest
number representable by an int. The exact value is machine-dependent.
The result of a reduction is an array expression, so reductions can be used as operands in an array expression:
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);
Note that this is not allowed:
Array<int,2> A; firstIndex i; secondIndex j; // Completely sum the array? int result = sum(sum(A, j), i);
You cannot reduce an array to zero dimensions! Instead, use one of the global functions described in the previous section.
Next: Where expr, Previous: Math functions 2, Up: Array Expressions [Contents][Index]