arrays-expr.texi

c++经典教材 Blitz++ v0.8

@cindex expression templates
@cindex Array expressions
@cindex Array no temporaries
@cindex temporaries
@cindex Array temporaries

Array expressions in Blitz++ are implemented using the @emph{expression
templates} technique.  Unless otherwise noted, expression evaluation will
never generate temporaries or multiple loops; an expression such as

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

A = B + C + D;
@end example

will result in code similar to

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

@node Array expressions
@section Expression evaluation order
@cindex Array expression evaluation order
@cindex expression evaluation order
@cindex order of expression evaluation
@cindex traversal order

A commonly asked question about Blitz++ is what order it uses to evaluate
array expressions.  For example, in code such as

@example
A(Range(2,10)) = A(Range(1,9))
@end example

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 @code{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:

@itemize @bullet

@item  try to go through the destination array in the order it is laid out
in memory (i.e.@:  row-major for row-major arrays, column-major for
column-major arrays).

@item  if the expression is a stencil, Blitz will do tiling to improve cache
use.  Under some circumstances blitz will even use a traversal based on a
hilbert curve (a fractal) for 3D arrays.  

@end itemize

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:

@itemize @bullet

@item  if your array is stored in reverse order, i.e.@: because of a
A.reverse(firstDim) or funny storage order, blitz will go through the array
from end to beginning in array coordinates, but from beginning to end in
memory locations.

@item  many compilers/architecture combinations are equally fast at reverse
order.  But blitz has a specialized version for stride = +1, and it would be
wasteful to also specialize for the case stride = -1.  So 1D arrays are
traversed from beginning to end (in memory storage order).

@end itemize

@section Expression operands
@cindex Array expression operands

An expression can contain any mix of these operands:

@itemize @bullet
@item     An array of any type, so long as it is of the same rank.
Expressions which contain a mixture of array types are handled through the
type promotion mechanism described below.

@item     Scalars of type @code{int}, @code{float}, @code{double},
@code{long double}, or @code{complex<T>}

@item     Index placeholders, described below

@item     Other expressions (e.g. @code{A+(B+C)})
@end itemize

@section Array operands

@unnumberedsubsec Using subarrays in an expression

@cindex Array using subarrays in expressions

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

@example
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;
@end example

@unnumberedsubsec Mixing arrays with different storage formats

@cindex Array expressions which mix arrays of 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:

@itemize @bullet

@item     Evaluation may be slower, since a different traversal order may be
used.

@item     If you are using index placeholders (see below) or reductions in
the expression, you may @strong{not} mix array objects with different
starting bases.  

@end itemize

@section Expression operators
@cindex operators, array expressions
@cindex Array operators
@cindex Array expression operators

These binary operators are supported:

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

Note: operator @code{<<} and @code{>>} are reserved for use in input/output.
If you need a bit-shift operation on arrays, you may define one yourself;
see @ref{User et}.

These unary operators are supported:

@example
- ~ !
@end example

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

@cindex Array operators applied elementwise
All operators are applied @emph{elementwise}.

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

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

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

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

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

@example
"../../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 ...
                     .
                     .
@end example

@cindex Array arrays of user type
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 @code{Polynomial}, and want to write code such as:

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

I would have to provide @code{operator*} for @code{Polynomial} by
implementing

@example
Polynomial Polynomial::operator*(Polynomial);)
@end example

or

@example
Polynomial operator*(Polynomial, Polynomial);)
@end example

@section Assignment operators

@cindex Array assignment operators
These assignment operators are supported:

@example
= += -= *= /= %= ^= &= |= >>= <<=
@end example

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

@itemize @bullet

@item    A constant (or literal) of type @code{T_numtype}

@item    An array of appropriate rank, possibly of a different numeric type

@item    An array expression, with appropriate rank and shape

@end itemize

@node Index placeholders
@section Index placeholders
@cindex Array index placeholders
@cindex index placeholders

Blitz++ provides objects called @emph{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 @code{firstIndex}, @code{secondIndex},
@code{thirdIndex}, ..., @code{tenthIndex}, @code{eleventhIndex}.
@findex firstIndex 
@findex secondIndex 
@findex thirdIndex
@findex fourthIndex
Here's an example of using an index placeholder:

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

This generates code which is similar to:

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

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

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

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

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

@cindex index placeholders multiple

@example
// 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)));
@end example

Here's a plot of the resulting array:

@center @image{sinsoid}
@center Array filled using an index placeholder expression.

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

@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)));
@end example

You can mix array operands and index placeholders:

@example
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 ]
@end example

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

@cindex tensor namespace
@cindex @code{i} (index placeholder)
@cindex @code{j} (index placeholder)
@cindex @code{k} (index placeholder)
@cindex @code{l} (index placeholder)
@cindex @code{m} (index placeholder)
@cindex @code{n} (index placeholder)

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

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

@findex blitz::tensor namespace

@example
namespace blitz::tensor;
@end example
  
when you would like to use them.  Alternately, you can just preface all the
index placeholders with @code{tensor::}, for example:

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

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

@section Type promotion
@cindex type promotion
@cindex Array 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 @code{Array<int>} to an
@code{Arrray<float>} will be promoted to @code{float}.  Generally, the
result is promoted to whichever type has greater precision.

@unnumberedsubsec Type promotion for user-defined types

@cindex type promotion for user-defined types
@cindex Array 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:

@itemize @bullet

@item     If both types are intrinsic (e.g. bool, int, float) then type
promotion follows the standard C rules.  This generally means that the
result will be promoted to whichever type has greater precision.  In
Blitz++, these rules have been extended to incorporate
@code{complex<float>}, @code{complex<double>}, and @code{complex<long
double>}.

@item     If one of the types is intrinsic (or complex), and the other is a
user-defined type, then the result is promoted to the user-defined type.

@item     If both types are user-defined, then the result is promoted to
whichever type requires more storage space (as determined by
@code{sizeof()}).  The rationale is that more storage space probably
indicates more precision.

@end itemize

If you wish to alter the default type promotion rules above, you have two
choices:

@itemize @bullet

@findex promote_trait 

@item     If the type promotion behaviour isn't dependent on the type of
operation performed, then you can provide appropriate specializations for
the class @code{promote_trait<A,B>} which is declared in
@code{<blitz/promote.h>}.

@item     If type promotion does depend on the type of operation, then you
will need to specialize the appropriate function objects in
@code{<blitz/ops.h>}.

@end itemize

Note that you can do these specializations in your own header files (you
don't have to edit @file{promote.h} or @file{ops.h}).

@unnumberedsubsec Manual casts

@cindex casts
@cindex Array 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++:

@example
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 ]
@end example

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

@findex cast()

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

The first argument to @code{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:

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

But this is not yet supported.

@node Math functions 1
@section Single-argument math functions

All of the functions described in this section are @emph{element-wise}.  For
example, this code--

@example
Array<float,2> A, B;   //
A = sin(B);
@end example

results in @code{A(i,j) = sin(B(i,j))} for all (i,j).

@unnumberedsubsec ANSI C++ math functions

These math functions are available on all platforms:

@cindex math functions
@cindex complex math functions