[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

3. Array Expressions

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

 
Array<int,1> A, B, C, D;    // ...

A = B + C + D;

will result in code similar to

 
for (int i=A.lbound(firstDim); i <= A.ubound(firstDim); ++i)
    A[i] = B[i] + C[i] + D[i];

[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

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

 
A(Range(2,10)) = A(Range(1,9))

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 A(1).

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:


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

3.2 Expression operands

An expression can contain any mix of these operands:


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

3.3 Array operands


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

Using subarrays in an expression

Subarrays may be used in an expression. For example, this code example performs a 5-point average on a two-dimensional array:

 
Array<float,2> A(64,64), B(64,64);   // ...
Range I(1,62), J(1,62);

A(I,J) = (B(I,J) + B(I+1,J) + B(I-1,J) 
                 + B(I,J+1) + B(I,J-1)) / 5;

[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

Mixing arrays with different storage formats

Arrays with different storage formats (for example, C-style and Fortran-style) can be mixed in the same expression. Blitz++ will handle the different storage formats automatically. However:


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

3.4 Expression operators

These binary operators are supported:

 
+ - * / % > < >= <= == != && || ^ & | 

Note: operator << and >> are reserved for use in input/output. If you need a bit-shift operation on arrays, you may define one yourself; see Declaring your own math functions on arrays.

These unary operators are supported:

 
- ~ !

The operators > < >= <= == != && || ! result in a bool-valued expression.

All operators are applied elementwise.

You can only use operators which are well-defined for the number type stored in the arrays. For example, bitwise XOR (^) is meaningful for integers, so this code is all right:

 
Array<int,3> A, B, C;   // ...
A = B ^ C;

Bitwise XOR is not meaningful on floating point types, so this code will generate a compiler error:

 
Array<float,1> A, B, C;   // ...
C = B ^ C;

Here's the compiler error generated by KAI C++ for the above code:

 
"../../blitz/ops.h", line 85: error: expression must have integral or enum type
  BZ_DEFINE_OP(BitwiseXor,^);
  ^
          detected during:
            instantiation of "blitz::BitwiseXor<float, float>::T_numtype
                      blitz::BitwiseXor<float, float>::apply(float, float)" at
                      line 210 of "../../blitz/arrayexpr.h"
            instantiation of ...
                     .
                     .

If you are creating arrays using a type you have created yourself, you will need to overload whatever operators you want to use on arrays. For example, if I create a class Polynomial, and want to write code such as:

 
Array<Polynomial,2> A, B, C;   // ...
C = A * B;

I would have to provide operator* for Polynomial by implementing

 
Polynomial Polynomial::operator*(Polynomial);)

or

 
Polynomial operator*(Polynomial, Polynomial);)

[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

3.5 Assignment operators

These assignment operators are supported:

 
= += -= *= /= %= ^= &= |= >>= <<=

An array object should appear on the left side of the operator. The right side can be:


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

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 firstIndex, secondIndex, thirdIndex, ..., tenthIndex, eleventhIndex. Here's an example of using an index placeholder:

 
Array<float,1> A(10);
firstIndex i;
A = i;

This generates code which is similar to:

 
for (int i=0; i < A.length(); ++i)
    A(i) = i;

Here's an example which fills an array with a sampled sine wave:

 
Array<float,1> A(16);
firstIndex i;

A = sin(2 * M_PI * i / 16.);

If your destination array has rank greater than 1, you may use multiple index placeholders:

 
// 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)));

Here's a plot of the resulting array:

sinsoid

Array filled using an index placeholder expression.

You can use index placeholder expressions in up to 11 dimensions. Here's a three dimensional example:

 
// 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)));

You can mix array operands and index placeholders:

 
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 ]

For your convenience, there is a namespace within blitz called tensor which declares all the index placeholders:

 
namespace blitz {
  namespace tensor {
    firstIndex i;
    secondIndex j;
    thirdIndex k;
     ...
    eleventhIndex t;
  }
}

So instead of declaring your own index placeholder objects, you can just say

 
namespace blitz::tensor;

when you would like to use them. Alternately, you can just preface all the index placeholders with tensor::, for example:

 
A = sin(2 * M_PI * tensor::i / 16.);

This will make your code more readable, since it is immediately clear that i is an index placeholder, rather than a scalar value.


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

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 Array<int> to an Arrray<float> will be promoted to float. Generally, the result is promoted to whichever type has greater precision.


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

Type promotion for user-defined types

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').


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

Manual casts

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++:

 
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 ]

The usual solution to this problem is to cast one of the operands to a floating type. For this purpose, Blitz++ provides a function cast(expr,type) which will cast the result of expr as type:

 
C = A / cast(B, float());   // Result: [ 0.5  1  1.5  0.714 ]

The first argument to cast() is an array or expression. The second argument is a dummy object of the type to which you want to cast. Once compilers support templates more thoroughly, it will be possible to use this cast syntax:

 
C = A / cast<float>(B);

But this is not yet supported.


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

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).


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

ANSI C++ math functions

These math functions are available on all platforms:

abs()

Absolute value

acos()

Inverse cosine. For real arguments, the return value is in the range [0, \pi].

arg()

Argument of a complex number (atan2(Im,Re)).

asin()

Inverse sine. For real arguments, the return value is in the range [-\pi/2, \pi/2].

atan()

Inverse tangent. For real arguments, the return value is in the range [-\pi/2, \pi/2]. See also atan2() in section Two-argument math functions.

ceil()

Ceiling function: smallest floating-point integer value not less than the argument.

cexp()

Complex exponential; same as exp().

conj()

Conjugate of a complex number.

cos()

Cosine. Works for complex<T>.

cosh()

Hyperbolic cosine. Works for complex<T>.

csqrt()

Complex square root; same as sqrt().

exp()

Exponential. Works for complex<T>.

fabs()

Same as abs().

floor()

Floor function: largest floating-point integer value not greater than the argument.

log()

Natural logarithm. Works for complex<T>.

log10()

Base 10 logarithm. Works for complex<T>.

pow2(), pow3(), pow4(), pow5(), pow6(), pow7(), pow8()

These functions compute an integer power. They expand to a series of multiplications, so they can be used on any type for which multiplication is well-defined.

sin()

Sine. Works for complex<T>.

sinh()

Hyperbolic sine. Works for complex<T>.

sqr()

Same as pow2(). Computes x*x. Works for complex<T>.

sqrt()

Square root. Works for complex<T>.

tan()

Tangent. Works for complex<T>.

tanh()

Hyperbolic tangent. Works for complex<T>.


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

IEEE/System V math functions

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 _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.).

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 `Blitz++/compiler/ieeemath.cpp' and `Blitz++/compiler/sysvmath.cpp'. This arrangement is unsatisfactory and will probably change in a future version.

You may have to link with -lm and/or -lmsaa to use these functions.

None of these functions are available for complex<T>.

acosh()

Inverse hyperbolic cosine

asinh()

Inverse hyperbolic sine

atanh()

Inverse hyperbolic tangent

_class()

Classification of floating point values. The return type is integer and will be one of:

FP_PLUS_NORM

Positive normalized, nonzero

FP_MINUS_NORM

Negative normalized, nonzero

FP_PLUS_DENORM

Positive denormalized, nonzero

FP_MINUS_DENORM

Negative denormalized, nonzero

FP_PLUS_ZERO

+0.0

FP_MINUS_ZERO

-0.0

FP_PLUS_INF

Positive infinity

FP_MINUS_INF

Negative infinity

FP_NANS

Signalling Not a Number (NaNS)

FP_NANQ

Quiet Not a Number (NaNQ)

cbrt()

Cubic root

expm1()

Computes exp(x)-1

erf()

Computes the error function: @erf(x) = 2/sqrt(Pi) * integral(exp(-t^2), t=0..x)

Note that for large values of the parameter, calculating can result in extreme loss of accuracy. Instead, use erfc().

erfc()

Computes the complementary error function erfc(x) = 1 - erf(x).

finite()

Returns a nonzero integer if the parameter is a finite number (i.e. not +INF, -INF, NaNQ or NaNS).

ilogb()

Returns an integer which is equal to the unbiased exponent of the parameter.

blitz_isnan()

Returns a nonzero integer if the parameter is NaNQ or NaNS (quiet or signalling Not a Number).

itrunc()

Round a floating-point number to a signed integer. Returns the nearest signed integer to the parameter in the direction of 0.

j0()

Bessel function of the first kind, order 0.

j1()

Bessel function of the first kind, order 1.

lgamma()

Natural logarithm of the gamma function. The gamma function is defined as: Gamma(x) = integral(e^(-t) * t^(x-1), t=0..infinity))

logb()

Returns a floating-point double that is equal to the unbiased exponent of the parameter.

log1p()

Calculates log(1+x), where x is the parameter.

nearest()

Returns the nearest floating-point integer value to the parameter. If the parameter is exactly halfway between two integer values, an even value is returned.

rint()

Rounds the parameter and returns a floating-point integer value. Whether 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().

rsqrt()

Reciprocal square root.

uitrunc()

Returns the nearest unsigned integer to the parameter in the direction of zero.

y0()

Bessel function of the second kind, order 0.

y1()

Bessel function of the second kind, order 1.

There may be better descriptions of these functions in your system man pages.


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

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:


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

ANSI C++ math functions

These math functions are available on all platforms, and work for complex numbers.

atan2(x,y)

Inverse tangent of (y/x). The signs of both parameters are used to determine the quadrant of the return value, which is in the range [-\pi, \pi]. Works for complex<T>.

blitz::polar(r,t)

Computes ; i.e. converts polar-form to Cartesian form complex numbers. The blitz:: scope qualifier is needed to disambiguate the ANSI C++ function template polar(T,T). This qualifier will hopefully disappear in a future version.

pow(x,y)

Computes x to the exponent y. Works for complex<T>.


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

IEEE/System V math functions

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.

copysign(x,y)

Returns the x parameter with the same sign as the y parameter.

drem(x,y)

Computes a floating point remainder. The return value r is equal to r = x - n * y, where n is equal to 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.

fmod(x,y)

Computes a floating point modulo remainder. The return value r is equal to r = x - n * y, where n is selected so that r has the same sign as x and magnitude less than abs(y). In order words, if x > 0, r is in the range [0, |y|], and if x < 0, r is in the range [-|y|, 0].

hypot(x,y)

Computes so that underflow does not occur and overflow occurs only if the final result warrants it.

nextafter(x,y)

Returns the next representable number after x in the direction of y.

remainder(x,y)

Equivalent to drem(x,y).

scalb(x,y)

Calculates.

unordered(x,y)

Returns a nonzero value if a floating-point comparison between x and y would be unordered. Otherwise, it returns zero.


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

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:

 
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);

[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

3.11 Tensor notation

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:

tensor1

Examples of array indexing, subarrays, and slicing.

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:

 
4 x 4
[         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.


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

3.12 Array reductions

Currently, Blitz++ arrays support two forms of reduction:


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

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:

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

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)

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 (where statements) to reduce over some part of an array. For example, sum(where(A > 0, A, 0)) sums only the positive elements in an array.


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

3.14 Partial Reductions

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;
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.


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

3.15 where statements

Blitz++ provides the where function as an array expression version of the ( ? : ) operator. The syntax is:

 
where(array-expr1, array-expr2, array-expr3)

Wherever 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));

[ << ] [ >> ]           [Top] [Contents] [Index] [ ? ]

This document was generated by Charlie & on April, 4 2005 using texi2html 1.76.