2.1: Getting started |
Array<T_numtype,N_rank>
. (In the future, a
TinyArray<T_numtype,...>
class may be implemented for
small arrays whose dimensions are known at compile time.
Other possible additions: sparse arrays and arrays with
symmetries.)
This array class provides a
dynamically allocated N-dimensional array, with reference
counting, arbitrary storage ordering, subarrays and slicing,
flexible expression handling, and many other useful
features.
2.1.1: Template parameters |
The Array
class takes two template parameters:
T_numtype
is the numeric type to be stored in the array.
T_numtype
can be an integral type (bool, char, unsigned char,
short int, short unsigned int, int, unsigned int, long,
unsigned long), floating point type (float, double, long double)
complex type (complex<float>, complex<double>, complex<long double>)
or any user-defined type with appropriate numeric semantics.
N_rank
is the rank (or dimensionality) of the array.
This should be a positive integer.
To use the Array
class, include the header
<blitz/array.h>
and use the namespace blitz
:
#include <blitz/array.h> using namespace blitz; Array<int,1> x; // A one-dimensional array of int Array<double,2> y; // A two-dimensional array of double . . Array<complex<float>, 12> z; // A twelve-dimensional array of complex<float>
When no constructor arguments are provided, the array is empty, and no memory is allocated. To create an array which contains some data, provide the size of the array as constructor arguments:
Array<double,2> y(4,4); // A 4x4 array of double
The contents of a newly-created array are garbage. To initialize the array, you can write:
y = 0;
and all the elements of the array will be set to zero. If the
contents of the array are known, you can initialize it using
a comma-delimited list of values. For example, this code
excerpt sets y
equal to a 4x4 identity matrix:
y = 1, 0, 0, 0 0, 1, 0, 0 0, 0, 1, 0 0, 0, 0, 1;
2.1.2: Array types |
The Array<T,N>
class supports a variety of arrays:
Array<int,1>
and Array<float,3>
Array<complex<float>,2>
Polynomial
, then Array<Polynomial,2>
is an array of
Polynomial
objects.
TinyVector
and TinyMatrix
,
in which each element is a fixed-size vector or array.
For example, Array<TinyVector<float,3>,3>
is a three-dimensional
vector field.
Array<Array<int,1>,1>
, in
which each element is a variable-length array.
2.1.3: A simple example |
Here's an example program which creates two 3x3 arrays, initializes them, and adds them:
#include <blitz/array.h> using namespace blitz; int main() { Array<float,2> A(3,3), B(3,3), C(3,3); A = 1, 0, 0, 2, 2, 2, 1, 0, 0; B = 0, 0, 7, 0, 8, 0, 9, 9, 9; C = A + B; cout << "A = " << A << endl << "B = " << B << endl << "C = " << C << endl; return 0; } |
and the output:
A = 3 x 3 1 0 0 2 2 2 1 0 0 B = 3 x 3 0 0 7 0 8 0 9 9 9 C = 3 x 3 1 0 7 2 10 2 10 9 9 |
2.1.4: Storage orders |
Blitz++ is very flexible about the way arrays are stored in memory.
The default storage format is row-major, C-style arrays whose indices start at zero.
Fortran-style arrays can also be created. Fortran arrays are stored in column-major order, and have indices which start at one. To create a Fortran-style array, use this syntax:
Array<int,2> A(3, 3, FortranArray<2>());
The last parameter, FortranArray<2>()
, tells the Array
constructor
to use a 2-dimensional Fortran-style array format. The template parameter
of FortranArray
must be the rank of the array being created.
FortranArray<N>
is a subclass of the type GeneralArrayStorage<N>
,
which encapsulates information about how an array is laid out in memory.
By altering the contents of a GeneralArrayStorage<N>
object, you
can lay out your arrays any way you want: the dimensions can
be ordered arbitrarily and stored in ascending or descending order,
and the starting indices can be arbitrary.
Creating custom array storage formats is described in a later section (2.14).
2.2: Public types |
Array
class declares these public types:
T_numtype
is the element type stored in the array. For example,
the type Array<double,2>::T_numtype
would be double
.
T_index
is a vector index into the array. The class
TinyVector
is used for this purpose.
T_array
is the array type itself (Array<T_numtype,N_rank>
)
T_iterator
is an iterator type. NB: this iterator is
not yet fully implemented, and is NOT STL compatible at the
present time.
2.3: Constructors |
2.3.1: Default constructor |
Array(); Array(GeneralArrayStorage<N_rank> storage)
The default constructor creates a C-style array of zero size. Any attempt to access data in the array may result in a run-time error, because there isn't any data to access!
An optional argument specifies a storage order for the array.
Arrays created using the default constructor can subsequently be given
data by the resize()
, resizeAndPreserve()
, or reference()
member functions.
2.3.2: Constructors which take extent parameters |
Array(int extent1); Array(int extent1, int extent2); Array(int extent1, int extent2, int extent3); ... Array(int extent1, int extent2, int extent3, ..., int extent11)
These constructors take arguments which specify the size of the array to be constructed. You should provide as many arguments as there are dimensions in the array. (If you provide fewer than N_rank arguments, the missing arguments will be filled in using the last provided argument. However, for code clarity, it makes sense to provide all N_rank parameters.)
An optional last parameter specifies a storage format:
Array(int extent1, GeneralArrayStorage<N_rank> storage); Array(int extent1, int extent2, GeneralArrayStorage<N_rank> storage); ...
For high-rank arrays, it may be convenient to use this constructor:
Array(const TinyVector<int, N_rank>& extent); Array(const TinyVector<int, N_rank>& extent, GeneralArrayStorage<N_rank> storage);
The argument extent
is a vector containing the extent (length) of
the array in each dimension. The optional second parameter indicates
a storage format. Note that you can construct TinyVector<int,N>
objects on the fly with the shape(i1,i2,...)
global function. For
example, Array<int,2> A(shape(3,5))
will create a 3x5 array.
A similar constructor lets you provide both a vector of base index values (lbounds) and extents:
Array(const TinyVector<int, N_rank>& lbound, const TinyVector<int, N_rank>& extent); Array(const TinyVector<int, N_rank>& lbound, const TinyVector<int, N_rank>& extent, GeneralArrayStorage<N_rank> storage);
The argument lbound
is a vector containing the base index value
(or lbound) of the array in each dimension.
The argument extent
is a vector containing the extent (length) of
the array in each dimension. The optional third parameter indicates
a storage format. As with the above constructor, you can use the
shape(i1,i2,...)
global function to create the lbound
and extent
parameters.
2.3.3: Constructors with Range arguments |
These constructors allow arbitrary bases (starting indices) to be set:
Array(Range r1); Array(Range r1, Range r2); Array(Range r1, Range r2, Range r3); ... Array(Range r1, Range r2, Range r3, ..., Range r11);
For example, this code:
Array<int,2> A(Range(10,20), Range(20,30));
will create an 11x11 array whose indices are 10..20 and 20..30. An optional last parameter provides a storage order:
Array(Range r1, GeneralArrayStorage<N_rank> storage); Array(Range r1, Range r2, GeneralArrayStorage<N_rank> storage); ...
2.3.4: Referencing another array |
This constructor makes a shared view of another array's data:
Array(Array<T_numtype, N_rank>& array);
After this constructor is used, both Array
objects refer to the
same data. Any changes made to one array will appear in the
other array. If you want to make a duplicate copy of an array,
use the copy()
member function.
2.3.5: Creating an array from pre-existing data |
These constructor create array objects from pre-existing data:
Array(T_numtype* dataFirst, TinyVector<int, N_rank> shape); Array(T_numtype* dataFirst, TinyVector<int, N_rank> shape, GeneralArrayStorage<N_rank> storage);
The first argument is a pointer to the array data. It should point to
the element of the array which is stored first in memory. The second
argument indicates the shape of the array. You can create this argument
using the shape()
function. For example:
double data[] = { 1, 2, 3, 4 }; Array<double,2> A(data, shape(2,2)); // Make a 2x2 array
The shape()
function takes N integer arguments and returns a
TinyVector<int,N>
.
By default, Blitz++ arrays are row-major. If you want to work with data which is stored in column-major order (e.g. a Fortran array), use the second version of the constructor:
Array<double,2> B(data, shape(2,2), FortranArray<2>());
Note that ownership of the array data is not acquired by Blitz++. This means that it is up to you to free the memory if necessary after you're done with it.
Another version of this constructor allows you to pass an arbitrary vector of strides:
Array(T_numtype* _bz_restrict dataFirst, TinyVector<int, N_rank> shape, TinyVector<int, N_rank> stride, GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
2.3.6: Interlacing arrays |
For some platforms, it can be advantageous to store a set of arrays
interlaced together in memory. Blitz++ provides support for
this through the routines interlaceArrays()
and allocateArrays()
.
An example:
Array<int,2> A, B; interlaceArrays(shape(10,10), A, B);
The first parameter of interlaceArrays()
is the shape for the
arrays (10x10). The subsequent arguments are the set of arrays to
be interlaced together. Up to 11 arrays may be interlaced.
All arrays must store the same data type and be of the same
rank.
In the above example, storage is allocated so that
A(0,0)
is followed immediately by B(0,0)
in memory,
which is folloed by A(0,1)
and B(0,1)
, and so on.
A related routine is allocateArrays()
, which has identical syntax:
Array<int,2> A, B; allocateArrays(shape(10,10), A, B);
Unlike interlaceArrays()
, which always interlaces the arrays,
the routine allocateArrays()
may or may not interlace them,
depending on whether interlacing is considered advantageous for your
platform. If the tuning flag BZ_INTERLACE_ARRAYS
is
defined in <blitz/tuning.h>
, then the arrays are
interlaced.
Note that the performance effects of interlacing are
unpredictable: in some situations it can be a benefit, and in
most others it can slow your code down substantially. You should
only use interlaceArrays()
after
running some benchmarks to determine whether interlacing
is beneficial for your particular algorithm and architecture.
2.3.7: A note about reference counting |
Blitz++ arrays use reference counting. When you create a new array,
a memory block is allocated. The Array
object acts like a handle
for this memory block. A memory block can be shared among multiple
Array
objects -- for example, when you take subarrays and slices.
The memory block keeps track of how many Array
objects are
referring to it. When a memory block is orphaned -- when no
Array
objects are referring to it -- it automatically
deletes itself and frees the allocated memory.
2.4: Indexing, subarrays, and slicing |
Indexing, subarrays and slicing all use the overloaded parenthesis operator().
As a running example, we'll consider the three dimensional array pictured below, which has index ranges (0..7, 0..7, 0..7). Shaded portions of the array show regions which have been obtained by indexing, creating a subarray, and slicing.
2.4.1: Indexing |
There are two ways to get a single element from
an array. The simplest is to provide a set of integer operands to
operator()
:
A(7,0,0) = 5; cout << "A(7,0,0) = " << A(7,0,0) << endl;This version of indexing is available for arrays of rank one through eleven. If the array object isn't const, the return type of
operator()
is a reference; if the array object is const, the
return type is a value.
You can also get an element by providing an operand
of type TinyVector<int,N_rank> where N_rank
is the
rank of the array object:
TinyVector<int,3> index; index = 7, 0, 0; A(index) = 5; cout << "A(7,0,0) = " << A(index) << endl;
This version of operator()
is also available in a const-overloaded
version.
It's possible to use fewer than N_rank
indices. However, missing
indices are assumed to be zero, which will cause bounds errors
if the valid index range does not include zero (e.g. Fortran arrays).
For this reason, and for code clarity, it's a bad idea to omit indices.
2.4.2: Subarrays |
You can obtain a subarray by providing Range
operands to operator()
. A Range
object represents a
set of regularly spaced index values. For example,
Array<int,3> B = A(Range(5,7), Range(5,7), Range(0,2));
The object B now refers to elements (5..7,5..7,0..2) of the array A.
The returned subarray is of type Array<T_numtype,N_rank>
. This
means that subarrays can be used wherever arrays can be: in expressions,
as lvalues, etc. Some examples:
// A three-dimensional stencil (used in solving PDEs) Range I(1,6), J(1,6), K(1,6); B = (A(I,J,K) + A(I+1,J,K) + A(I-1,J,K) + A(I,J+1,K) + A(I,J-1,K) + A(I,J+1,K) + A(I,J,K+1) + A(I,J,K-1)) / 7.0; // Set a subarray of A to zero A(Range(5,7), Range(5,7), Range(5,7)) = 0.;
The bases of the subarray are equal to the bases of the original array:
Array<int,2> D(Range(1,5), Range(1,5)); // 1..5, 1..5 Array<int,2> E = D(Range(2,3), Range(2,3)); // 1..2, 1..2
An array can be used on both sides of an expression only if the subarrays don't overlap. If the arrays overlap, the result may depend on the order in which the array is traversed.
2.4.3: Slicing |
A combination of integer and Range operands produces a slice. Each integer operand reduces the rank of the array by one. For example:
Array<int,2> F = A(Range::all(), 2, Range::all()); Array<int,1> G = A(2, 7, Range::all());
Range and integer operands can be used in any combination, for arrays up to rank 11.
Note: Using a combination of integer and Range operands requires
a newer language feature (partial ordering of member templates) which
not all compilers support. If your compiler does provide this
feature, BZ_PARTIAL_ORDERING
will be defined in <blitz/config.h>
.
If not, you can use this workaround:
Array<int,3> F = A(Range::all(), Range(2,2), Range::all()); Array<int,3> G = A(Range(2,2), Range(7,7), Range::all());
2.4.4: More about Range objects |
A Range
object represents an ordered set of uniformly spaced
integers. Here are some examples of using Range objects to obtain
subarrays:
Array<int,1> A(7); A = 0, 1, 2, 3, 4, 5, 6; cout << A(Range::all()) << endl // [ 0 1 2 3 4 5 6 ] << A(Range(3,5)) << endl // [ 3 4 5 ] << A(Range(3,toEnd)) << endl // [ 3 4 5 6 ] << A(Range(fromStart,3)) << endl // [ 0 1 2 3 ] << A(Range(1,5,2)) << endl // [ 1 3 5 ] << A(Range(5,1,-2)) << endl // [ 5 3 1 ] << A(Range(fromStart,toEnd,2)) << endl; // [ 0 2 4 6 ]
The optional third constructor argument specifies
a stride. For example, Range(1,5,2)
refers to elements
[1 3 5]. Strides can also be negative: Range(5,1,-2)
refers to elements [5 3 1].
Here's an example of using strides with a two-dimensional array:
#include <blitz/array.h> using namespace blitz; int main() { Array<int,2> A(8,8); A = 0; Array<int,2> B = A(Range(1,7,3), Range(1,5,2)); B = 1; cout << "A = " << A << endl; return 0; } |
Here's an illustration of the B
subarray:
And the program output:
0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 0 |
2.4.5: A note about assignment |
The assignment operator (=) always results in the expression on the right-hand side (rhs) being copied to the lhs (i.e. the data on the lhs is overwritten with the result from the rhs). This is different from some array packages in which the assignment operator makes the lhs a reference (or alias) to the rhs. To further confuse the issue, the copy constructor for arrays does have reference semantics. Here's an example which should clarify things:
Array<int,1> A(5), B(10); A = B(Range(0,4)); // Statement 1 Array<int,1> C = B(Range(0,4)); // Statement 2
Statement 1 results in a portion of B
's data being copied into
A
. After Statement 1, both A
and B
have their own
(nonoverlapping) blocks of data. Contrast this behaviour with that
of Statement 2, which is not an assignment (it uses the copy
constructor). After Statement 2 is executed, the array C
is a
reference (or alias) to B's data.
So to summarize: If you want to copy the rhs, use an assignment operator. If you want to reference (or alias) the rhs, use the copy constructor (or alternately, the reference() member function in 2.6.2).
2.4.6: An example |
#include <blitz/array.h> using namespace blitz; int main() { Array<int,2> A(6,6), B(3,3); // Set the upper left quadrant of A to 5 A(Range(0,2), Range(0,2)) = 5; // Set the upper right quadrant of A to an identity matrix B = 1, 0, 0, 0, 1, 0, 0, 0, 1; A(Range(0,2), Range(3,5)) = B; // Set the fourth row to 1 A(3, Range::all()) = 1; // Set the last two rows to 0 A(Range(4, toEnd), Range::all()) = 0; // Set the bottom right element to 8 A(5,5) = 8; cout << "A = " << A << endl; return 0; } |
The output:
A = 6 x 6 5 5 5 1 0 0 5 5 5 0 1 0 5 5 5 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 8 |
2.5: Debug mode |
BZ_DEBUG
. For most compilers, the command
line argument -DBZ_DEBUG
should work.
In debugging mode, your programs will run very slowly. This is because Blitz++ is doing lots of precondition checking and bounds checking. When it detects something fishy, it will likely halt your program and display an error message.
For example, this program attempts to access an element of a 4x4 array which doesn't exist:
#include <blitz/array.h> using namespace blitz; int main() { Array<complex<float>, 2> Z(4,4); // Since this is a C-style array, the valid index // ranges are 0..3 and 0..3 Z(4,4) = complex<float>(1.0, 0.0); return 0; } |
When compiled with -DBZ_DEBUG
, the out of bounds indices are
detected and an error message results:
[Blitz++] Precondition failure: Module ../../blitz/array.h line 1070 Array index out of range: (4, 4) Lower bounds: [ 0 0 ] Upper bounds: [ 3 3 ] Assertion failed: __EX, file ../../blitz/array.h, line 1070 IOT/Abort trap (core dumped) |
Precondition failures send their error messages to the standard
error stream (cerr
). After displaying the error message,
assert(0)
is invoked.
2.6: Member functions |
2.6.1: A note about dimension parameters |
Several of the member functions take a dimension parameter which
is an integer in the range 0 .. N_rank - 1. For example, the method
extent(int n)
returns the extent (or length) of the
array in dimension n
.
These parameters are problematic:
reverse()
member function won't stand a chance of
understanding what A.reverse(2)
does.
As a solution to this problem, Blitz++ provides a series of symbolic constants which you can use to refer to dimensions:
const int firstDim = 0; const int secondDim = 1; const int thirdDim = 2; . . const int eleventhDim = 10;
These symbols should be used in place of the numerals 0, 1, ... N_rank - 1. For example:
A.reverse(thirdDim);
This code is clearer: you can see that the parameter refers to a dimension, and it isn't much of a leap to realize that it's reversing the element ordering in the third dimension.
If you find firstDim
, secondDim
, ... aesthetically unpleasing,
there are equivalent symbols firstRank
, secondRank
,
thirdRank
, ..., eleventhRank
.
Why stop at eleven?
The symbols had to stop somewhere, and eleven seemed an appropriate place to stop. Besides, if you're working in more than eleven dimensions your code is going to be confusing no matter what help Blitz++ provides.
2.6.2: Member function descriptions |
const TinyVector<int, N_rank>& base() const;
int base(int dimension) const;
The base of a dimension is the first valid index value. A typical
C-style array will have base of zero; a Fortran-style array will have
base of one. The base can be different for each dimension, but only
if you deliberately use a Range-argument constructor
or design a custom storage ordering.
The first version returns a reference to the vector of base values.
The second version returns the base for just one dimension; it's
equivalent to the lbound()
member function. See the
note on dimension parameters such as firstDim
above.
int cols() const;
int columns() const;
Both of these functions return the extent of the array in the
second dimension. Equivalent to extent(secondDim)
.
See also rows()
and depth()
.
Array<T_numtype, N_rank> copy() const;
This method creates a copy of the array's data, using the same
storage ordering as the current array. The returned array is
guaranteed to be stored contiguously in memory, and to be
the only object referring to its memory block (i.e. the data
isn't shared with any other array object).
const T_numtype* [restrict] data() const;
T_numtype* [restrict] data();
const T_numtype* [restrict] dataZero() const;
T_numtype* [restrict] dataZero();
const T_numtype* [restrict] dataFirst() const;
T_numtype* [restrict] dataFirst();
These member functions all return pointers to the array data.
The NCEG restrict
qualifier is used only if your compiler
supports it.
If you're working with the default storage order (C-style arrays
with base zero), you'll only need to use data()
.
Otherwise, things get complicated:
data()
returns a pointer to the element whose indices
are equal to the array base. With a C-style array, this means
the element (0,0,...,0); with a Fortran-style array, this means
the element (1,1,...,1). If A
is an array object,
A.data()
is equivalent to (&A(A.base(firstDim), A.base(secondDim), ...)).
If any of the dimensions are stored in reverse order,
data()
will not refer to the element which comes first in
memory.
dataZero()
returns a pointer to the element (0,0,...,0),
even if such an element does not exist in the array. What's
the point of having such a pointer? Say you want to access
the element (i,j,k). If you add to the pointer the dot product
of (i,j,k) with the stride vector (A.stride()
), you get
a pointer to the element (i,j,k).
dataFirst()
returns a pointer to the element of the array
which comes first in memory. Note however, that under some
circumstances (e.g. subarrays), the data will
not be stored contiguously in memory. You have to be very
careful when meddling directly with an array's data.
Other relevant functions are: isStorageContiguous()
and
zeroOffset()
.
int depth() const;
Returns the extent of the array in the third dimension. This
function is equivalent to extent(thirdDim)
.
See also rows()
and columns()
.
int dimensions() const;
Returns the number of dimensions (rank) of the array. The return
value is the second template parameter (N_rank) of the
Array
object. Same as rank()
.
RectDomain<N_rank> domain() const;
Returns the domain of the array. The domain consists of
a vector of lower bounds and a vector of upper bounds for
the indices. NEEDS_WORK-- need a section to explain
methods of RectDomain<N>
.
int extent(int dimension) const;
The first version the extent (length) of the array in the specified dimension.
See the note about dimension
parameters such as firstDim
in the previous section.
Array<T_numtype2,N_rank> extractComponent(T_numtype2,
int componentNumber, int numComponents);
This method returns an array view of a single component of a
multicomponent array. In a multicomponent array, each element
is a tuple of fixed size. The components are numbered 0, 1, ...,
numComponents-1
. Example:
Array<TinyVector<int,3>,2> A(128,128); // A 128x128 array of int[3] Array<int,2> B = A.extractComponent(int(), 1, 3);
Now the B array refers to the 2nd component of every element in A.
Note: for complex arrays, special global functions real(A)
and
imag(A)
are provided to obtain real and imaginary components of
an array. See the Global Functions section.
void free();
This method resizes an array to zero size. If the array data is
not being shared with another array object, then it is freed.
bool isMajorRank(int dimension) const;
Returns true if the dimension has the largest stride. For
C-style arrays (the default), the first dimension always has
the largest stride. For Fortran-style arrays, the last dimension
has the largest stride. See also isMinorRank()
below and
the note about dimension parameters such as firstDim
in the
previous section.
bool isMinorRank(int dimension) const;
Returns true if the dimension does not have the largest stride.
See also isMajorRank()
.
bool isRankStoredAscending(int dimension) const;
Returns true if the dimension is stored in ascending order in memory.
This is the default. It will only return false if you have reversed
a dimension using reverse()
or have created a custom storage order
with a descending dimension.
bool isStorageContiguous() const;
Returns true if the array data is stored contiguously in memory.
If you slice the array or work on subarrays, there can be
skips -- the array data is interspersed with other data not
part of the array. See also the various data..()
functions.
If you need to ensure that the storage is contiguous, try
reference(copy())
.
int lbound(int dimension) const;
TinyVector<int,N_rank> lbound() const;
The first version returns the lower bound of the valid index range
for a dimension. The second version returns a vector of lower bounds
for all dimensions.
The
lower bound is the first valid index value. If you're
using a C-style array (the default), the lbound will be zero;
Fortran-style arrays have lbound equal to one. The lbound can
be different for each dimension, but only if you deliberately
set them that way using a Range constructor or a custom storage
ordering. This function is equivalent to base(dimension)
.
See the note about dimension parameters such as firstDim
in the previous section.
void makeUnique();
If the array's data is being shared with another Blitz++ array
object, this member function creates a copy so the array object has
a unique view of the data. Note: if the array was created
from pre-existing data (by passing a data pointer to
the array constructor), this method will create a copy.
int numElements() const;
Returns the total number of elements in the array, calculated
by taking the product of the extent in each dimension. Same
as size()
.
const TinyVector<int, N_rank>& ordering() const;
int ordering(int storageRankIndex) const;
These member functions return information about how the data is
ordered in memory. The first version returns the complete ordering
vector; the second version returns a single element from the
ordering vector. The argument for the second version must be
in the range 0 .. N_rank-1.
The ordering vector is a list of dimensions in increasing order
of stride; ordering(0)
will return the dimension number
with the smallest stride, and ordering(N_rank-1)
will return
the dimension number with largest stride. For a C-style
array, the ordering vector contains the elements
(N_rank-1, N_rank-2, ..., 0). For a Fortran-style array,
the ordering vector is (0, 1, ..., N_rank-1). See also the
description of custom storage orders in section 2.14.
int rank() const;
Returns the rank (number of dimensions) of the array. The return
value is equal to N_rank. Equivalent to dimensions()
.
void reference(Array<T_numtype,N_rank>& A);
This causes the array to adopt another array's data as its own.
After this member function is used, the array object and the array
A
are indistinguishable -- they have identical sizes, index
ranges, and data. The data is shared between the two arrays.
void resize(int extent1, ...);
void resize(const TinyVector<int,N_rank>&);
These functions resize an array to the specified size. If
the array is already the size specified, then no memory is
allocated. After resizing, the contents of the array are
garbage. See also resizeAndPreserve()
.
void resizeAndPreserve(int extent1, ...);
These functions resize an array to the specified size. If
the array is already the size specified, then no change
occurs (the array is not reallocated and copied).
The contents of the array are preserved whenever possible;
if the new array size is smaller, then some data will be
lost. Any new elements created by resizing the array
are left uninitialized.
Array<T,N> reverse(int dimension);
void reverseSelf(int dimension);
This method reverses the array in the specified dimension.
For example, if reverse(firstDim)
is invoked on a
2-dimensional array, then the ordering of rows in the
array will be reversed; reverse(secondDim)
would
reverse the order of the columns. Note that this is
implemented by twiddling the strides of the array, and
doesn't cause any data copying. The first version
returns a reversed "view" of the array data; the second
version applies the reversal to the array itself.
int rows() const;
Returns the extent (length) of the array in the first dimension.
This function is equivalent to extent(firstDim)
.
See also columns()
, and depth()
.
int size() const;
Returns the total number of elements in the array, calculated
by taking the product of the extent in each dimension. Same
as numElements()
.
const TinyVector<int, N_rank>& shape() const;
Returns the vector of extents (lengths) of the array.
const TinyVector<int, N_rank>& stride() const;
int stride(int dimension) const;
The first version returns the stride vector; the second version returns
the stride associated with a dimension. A stride is the distance
between pointers to two array elements which are adjacent in a
dimension. For example, A.stride(firstDim)
is equal to
&A(1,0,0) - &A(0,0,0)
. The stride for the second dimension,
A.stride(secondDim)
, is equal to &A(0,1,0) - &A(0,0,0)
,
and so on. For more information about strides, see the description
of custom storage formats in Section 2.14.
See also the description of parameters like firstDim
and
secondDim
in the previous section.
Array<T,N> transpose(int dimension1, int dimension2, ...);
void transposeSelf(int dimension1,
int dimension2, ...);
These methods permute the dimensions of the array. The
dimensions of the array are reordered so that the first dimension is
dimension1
, the second is dimension2
, and so on. The
arguments should be a permutation of the symbols firstDim,
secondDim, ...
. Note that this is implemented by twiddling the
strides of the array, and doesn't cause any data copying.
The first version returns a transposed "view"
of the array data; the second version transposes the array itself.
int ubound(int dimension) const;
TinyVector<int,N_rank> ubound() const;
The first version returns the upper bound of the valid index range for a dimension.
The second version returns a vector of upper bounds for all dimensions.
The upper bound is the last valid index value. If you're using
a C-style array (the default), the ubound will be equal to the
extent(dimension)-1
. Fortran-style arrays will have
ubound equal to extent(dimension)
. The ubound can be
different for each dimension. The return value of
ubound(dimension)
will always be equal to
lbound(dimension) + extent(dimension) - 1
.
See the note about dimension parameters such as firstDim
in the previous section.
int zeroOffset() const;
This function has to do with the storage of arrays in memory. You
may want to refer to the description of the data..()
member
functions and of custom storage orders in Section 2.14
for clarification. The return value of zeroOffset()
is the
distance from the first element in the array to the (possibly nonexistant)
element (0,0,...,0)
. In this context, "first element" returns to the
element (base(firstDim),base(secondDim),...)
.
2.7: Global functions |
void allocateArrays(TinyVector<int,N>& shape, Array<T,N>& A, Array<T,N>& B, ...);This function will allocate interlaced arrays, but only if interlacing is desirable for your architecture. This is controlled by the
BZ_INTERLACE_ARRAYS
flag in <blitz/tuning.h>. You can provide up
to 11 arrays as parameters. Any views currently associated with the
array objects are lost. Here is a typical use:
Array<int,2> A, B, C; allocateArrays(shape(64,64),A,B,C);
If array interlacing is enabled, then the arrays are stored in memory like this: A(0,0), B(0,0), C(0,0), A(0,1), B(0,1), ... If interlacing is disabled, then the arrays are allocated in the normal fashion: each array has its own block of memory. Once interlaced arrays are allocated, they can be used just like regular arrays.
void cycleArrays(Array<T,N>& A, Array<T,N>& B);
void cycleArrays(Array<T,N>& A, Array<T,N>& B,
Array<T,N>& C);
void cycleArrays(Array<T,N>& A, Array<T,N>& B,
Array<T,N>& C, Array<T,N>& D);
void cycleArrays(Array<T,N>& A, Array<T,N>& B,
Array<T,N>& C, Array<T,N>& D,
Array<T,N>& E);
These routines are useful for time-stepping PDEs. They take a set of
arrays such as [A,B,C,D] and cyclically rotate them to [B,C,D,A]; i.e.
the A array then refers to what was B's data, the B array refers to what
was C's data, and the D array refers to what was A's data. These functions
operate in constant time, since only the handles change (i.e. no data
is copied; only pointers change).
Array<T,N> imag(Array<complex<T>,N>&);
This method returns a view of the imaginary portion of the array.
void interlaceArrays(TinyVector<int,N>& shape,
Array<T,N>& A, Array<T,N>& B, ...);
This function is similar to allocateArrays()
above, except that
the arrays are always interlaced, regardless of the setting of
the BZ_INTERLACE_ARRAYS
flag.
Array<T,N> real(Array<complex<T>,N>&);
This method returns a view of the real portion of the array.
TinyVector<int,1> shape(int L);
TinyVector<int,2> shape(int L, int M);
TinyVector<int,3> shape(int L, int M, int N);
TinyVector<int,4> shape(int L, int M, int N, int O);
... [up to 11 dimensions]
These functions may be used to create shape parameters.
They package the set of integer arguments as a TinyVector
of
appropriate length. For an example use, see allocateArrays()
above.
2.8: Expressions |
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];
2.8.1: Expression operands |
An expression can contain any mix of these operands:
int
, float
, double
, long double
,
or complex<T>
A+(B+C)
)
2.8.2: Array operands |
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;
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:
2.8.3: Expression operators |
These binary operators are supported:
+ - * / % > < >= <= == != && || ^ & | >> <<
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);
2.8.4: Assignment operators |
These assignment operators are supported:
= += -= *= /= %= ^= &= |= >>= <<=
An array object should appear on the left side of the operator. The right side can be:
T_numtype
In the near future, random number generators and index placeholders will be supported.
2.8.5: 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:
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
verb(using namespace blitz::tensor;
)
when you would like to use them. Alternately, you
can just preface all the index placeholders with
tensor::
, for example:
verb(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.
2.8.6: 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.
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:
complex<float>
, complex<double>
, and
complex<long double>
.
sizeof()
). The rationale is that
more storage space probably indicates more precision.
If you wish to alter the default type promotion rules above, you have two choices:
promote_trait<A,B>
which
is declared in <blitz/promote.h>
.
<blitz/ops.h>
.
Note that you can do these specializations in your own header files (you don't have to edit promote.h or ops.h).
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.
2.8.7: 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).
ANSI C++ math functions
These math functions are available on all platforms:
complex<T>
.
complex<T>
.
atan2()
in section
2.8.8. Works for 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
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. Blitz++ takes care of this in
<blitz/blitz.h>
. If these definitions
cause trouble for you, you can compile with -DBZ_DISABLE_XOPEN_SOURCE
which will prevent these symbols from being defined.
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>
.
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()
.
There may be better descriptions of these functions in your system man pages.
2.8.8: Two-argument math functions |
The math functions described in this subsection take two arguments. Most combinations of these types may be used as arguments:
float
, double
, long double
,
or complex<T>
ANSI C++ math functions
These math functions are available on all platforms, and work for complex numbers.
complex<T>
.
complex<T>
.
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.
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.
2.8.9: 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:
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.
2.8.10: Array reductions |
Currently, Blitz++ arrays support two forms of reduction: (A future version will likely support block reductions.)
2.8.11: 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:
tiny(int())
(INT_MIN)
is returned.
2.8.12: 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:
tiny(int())
(INT_MIN) 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 ]
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.
2.8.13: 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));
2.9: Stencil objects |
2.9.1: Motivation: a nicer notation for stencils |
Suppose we wanted to implement the 3-D acoustic wave equation using finite differencing. Here is how a single iteration would look using subarray syntax:
Range I(1,N-2), J(1,N-2), K(1,N-2); P3(I,J,K) = (2-6*c(I,J,K)) * P2(I,J,K) + c(I,J,K)*(P2(I-1,J,K) + P2(I+1,J,K) + P2(I,J-1,K) + P2(I,J+1,K) + P2(I,J,K-1) + P2(I,J,K+1)) - P1(I,J,K);
This syntax is a bit klunky. With stencil objects, the implementation becomes:
BZ_DECLARE_STENCIL4(acoustic3D_stencil,P1,P2,P3,c) P3 = 2 * P2 + c * Laplacian3D(P2) - P1; BZ_END_STENCIL . . applyStencil(acoustic3D_stencil(), P1, P2, P3, c);
2.9.2: Declaring stencil objects |
A stencil declaration may not be inside a function. It can appear inside a class declaration (in which case the stencil object is a nested type).
Stencil objects are declared using the macros
BZ_DECLARE_STENCIL1
, BZ_DECLARE_STENCIL2
, etc.
The number suffix is how many arrays are involved
in the stencil (in the above example, 4 arrays-- P1,P2,P3,c --
are used, so the macro BZ_DECLARE_STENCIL4 is invoked).
The first argument is a name for the stencil object. Subsequent arguments are names for the arrays on which the stencil operates.
After the stencil declaration, the macro BZ_END_STENCIL
must appear.
In between the two macros, you can have multiple assignment statements, if/else/elseif constructs, function calls, loops, etc.
Here are some simple examples:
BZ_DECLARE_STENCIL2(smooth2D,A,B) A = (B(0,0) + B(0,1) + B(0,-1) + B(1,0) + B(-1,0)) / 5.0; BZ_END_STENCIL BZ_DECLARE_STENCIL4(acoustic2D,P1,P2,P3,c) A = 2 * P2 + c * (-4 * P2(0,0) + P2(0,1) + P2(0,-1) + P2(1,0) + P2(-1,0)) - P1; BZ_END_STENCIL BZ_DECLARE_STENCIL8(prop2D,E1,E2,E3,M1,M2,M3,cE,cM) E3 = 2 * E2 + cE * Laplacian2D(E2) - E1; M3 = 2 * M2 + cM * Laplacian2D(M2) - M1; BZ_END_STENCIL BZ_DECLARE_STENCIL3(smooth2Db,A,B,c) if ((c > 0.0) && (c < 1.0)) A = c * (B(0,0) + B(0,1) + B(0,-1) + B(1,0) + B(-1,0)) / 5.0 + (1-c)*B; else A = 0; BZ_END_STENCIL
Currently, a stencil can take up to 11 array parameters.
You can use the notation A(i,j,k)
to read the element at
an offset (i,j,k)
from the current element. If you omit
the parentheses (i.e. as in "A"), then the current element
is read.
You can invoke stencil operators which calculate finite differences and laplacians.
2.9.3: Stencil operators |
This section lists all the stencil operators provided by Blitz++.
They assume that an array represents evenly spaced data points
separated by a distance of h
. A 2nd-order accurate operator
has error term O(h^2); a 4th-order accurate operator has
error term O(h^4).
All of the stencils have factors associated with them. For
example, the central12
operator is a discrete first derivative
which is 2nd-order accurate. Its factor is 2h; this means that
to get the first derivative of an array A, you need to use
central12(A,firstDim)/(2*h)
. Typically when designing
stencils, one factors out all of the h
terms for efficiency.
The factor terms always consist of an integer multiplier (often 1)
and a power of h
. For ease of use, all of the operators
listed below are provided in a second ``normalized'' version in which
the integer multiplier is 1. The normalized versions have
an n
appended to the name, for example central12n
is
the normalized version of central12
, and has factor h
instead of 2h
.
These operators are defined in blitz/array/stencilops.h
if
you wish to see the implementation.
Note that the above are available in normalized versions
central12n
, central22n
, ..., central44n
which
have factors of h
, h^2
, h^3
, or h^4
as
appropriate.
Note that the above are available in normalized versions
forward11n
, forward21n
, ..., forward42n
which
have factors of h
, h^2
, h^3
, or h^4
as
appropriate.
Note that the above are available in normalized versions
backward11n
, backward21n
, ..., backward42n
which
have factors of h
, h^2
, h^3
, or h^4
as
appropriate.
Note that the above are available in normalized versions
Laplacian2D4n
, Laplacian3D4n
which have factors
h^2.
These return TinyVector
s of the appropriate numeric type and length:
2.9.3.6: Grad-squared operators
There are also grad-squared operators, which return TinyVector
s of
second derivatives:
Note that the above are available in normalized versions
gradSqr2Dn
, gradSqr2D4n
, gradSqr3Dn
, gradSqr3D4n
which have factors h^2.
The curl operators return three-dimensional TinyVector
s of the
appropriate numeric type:
Note that the above are available in normalized versions curln
and
curl4n
, which have factors of h
.
The divergence operators return a scalar value.
Note that the above are available in normalized versions
divn
and div4n
, which have factors of h
.
2.9.3.9: Mixed partial derivatives
There are also normalized versions of the above, mixed22n
and
mixed24n
which have factors h^2
.
2.9.4: Declaring your own stencil operators |
You can declare your own stencil operators using the macro
BZ_DECLARE_STENCIL_OPERATOR1
. For example, here is the
declaration of Laplacian2D
:
BZ_DECLARE_STENCIL_OPERATOR1(Laplacian2D, A) return -4*A(0,0) + A(-1,0) + A(1,0) + A(0,-1) + A(0,1); BZ_END_STENCIL_OPERATOR
To declare a difference operator, use this syntax:
BZ_DECLARE_DIFF(central12,A) { return A.shift(1,dim) - A.shift(-1,dim); }
The method shift(offset,dim)
retrieves the element at
offset
in dimension dim
.
Stencil operator declarations cannot occur inside a function or class.
2.9.5: Composing stencils |
Watch this space...
2.9.6: Applying a stencil |
The syntax for applying a stencil is:
applyStencil(stencilname(),A,B,C...,F);
Where stencilname
is the name of the stencil, and
A,B,C,...,F
are the arrays on which the stencil
operates.
For examples, see examples/stencil.cpp
and examples/stencil2.cpp
.
Blitz++ interrogates the stencil object to find out how large its footprint is. It only applies the stencil over the region of the arrays where it won't overrun the boundaries.
2.9.7: What still needs work |
2.10: Multicomponent and complex arrays |
Here are some examples of multicomponent arrays:
// A 3-dimensional array; each element is a length 3 vector of float Array<TinyVector<float,3>,3> A; // A complex 2-dimensional array Array<complex<double>,2> B; // A 2-dimensional image containing RGB tuples struct RGB24 { unsigned char r, g, b; }; Array<RGB24,2> C;
2.10.1: Extracting components |
Blitz++ provides some special support for such arrays. The most important is the ability to extract a single component. For example:
Array<TinyVector<float,3>,2> A(128,128); Array<float,2> B = A.extractComponent(float(), 1, 3); B = 0;
The call to extractComponent
returns an array of floats; this
array is a view of the second component of each element of A.
The arguments of extractComponent
are: (1) the type of the
component (in this example, float); (2) the component number
to extract (numbered 0, 1, ... N-1); and (3) the number of
components in the array.
This is a little bit messy, so Blitz++ provides a handy shortcut
using operator[]
:
Array<TinyVector<float,3>,2> A(128,128); A[1] = 0;
The number inside the square brackets is the component number.
However, for this operation to work, Blitz++ has to already know
how many components there are, and what type they are. It
knows this already for TinyVector
and complex<T>
.
If you use your own type, though, you will have to tell
Blitz++ this information using the macro BZ_DECLARE_MULTICOMPONENT_TYPE()
.
This macro has three arguments:
BZ_DECLARE_MULTICOMPONENT_TYPE(T_element, T_componentType, numComponents)
T_element
is the element type of the array. T_componentType
is the type of the components of that element. numComponents
is
the number of components in each element.
An example will clarify this. Suppose we wanted to make a
colour image, stored in 24-bit HSV (hue-saturation-value) format.
We can make a class HSV24
which represents a single pixel:
#include <blitz/array.h> using namespace blitz; class HSV24 { public: // These constants will makes the code below cleaner; we can // refer to the components by name, rather than number. static const int hue=0, saturation=1, value=2; HSV24() { } HSV24(int hue, int saturation, int value) : h_(hue), s_(saturation), v_(value) { } // Some other stuff here, obviously private: unsigned char h_, s_, v_; };
Right after the class declaration, we will invoke the
macro BZ_DECLARE_MULTICOMPONENT_TYPE
to tell Blitz++
about HSV24:
// HSV24 has 3 components of type unsigned char BZ_DECLARE_MULTICOMPONENT_TYPE(HSV24, unsigned char, 3);
Now we can create HSV images and modify the individual components:
int main() { Array<HSV24,2> A(128,128); // A 128x128 HSV image ... // Extract a greyscale version of the image Array<unsigned char,2> A_greyscale = A[HSV24::value]; // Bump up the saturation component to get a // pastel effect A[HSV24::saturation] *= 1.3; // Brighten up the middle of the image Range middle(32,96); A[HSV24::value](middle,middle) *= 1.2; }
2.10.2: Special support for complex arrays |
Since complex arrays are used frequently, Blitz++ provides two special methods for getting the real and imaginary components:
Array<complex<float>,2> A(32,32); real(A) = 1.0; imag(A) = 0.0;
The function real(A)
returns an array view of the
real component; imag(A)
returns a view of the imaginary
component.
2.10.3: Zipping together expressions |
Blitz++ provides a function zip()
which lets you
combine two or more expressions into a single component.
For example, you can combine two real expressions into
a complex expression, or three integer expressions into
an HSV24 expression. The function has this syntax:
resultexpr zip(expr1, expr2, T_element) resultexpr zip(expr1, expr2, expr3, T_element) ** not available yet resultexpr zip(expr1, expr2, expr3, expr4, T_element) ** not available yet
The types resultexpr
, expr1
and expr2
are
array expressions. The third argument is the type you
want to create. For example:
int N = 16; Array<complex<float>,1> A(N); Array<float,1> theta(N); ... A = zip(cos(theta), sin(theta), complex<float>());
The above line is equivalent to:
for (int i=0; i < N; ++i) A[i] = complex<float>(cos(theta[i]), sin(theta[i]));
2.11: Indirection |
I
of
rows and a list J
of columns, and you want to modify
the array at all (i,j) positions where i is in I
and
j is in J
. This is a cartesian product of
the index sets I
and J
.
In all cases, Blitz++ expects a Standard Template Library
container. Some useful STL containers are
list<>
, vector<>
, deque<>
and set<>
.
Documentation of these classes is often provided
with your compiler, or see also the good documentation
at http://www.objectspace.com/toolkits/documentation/stdref/index.html.
STL containers are used because they are widely available
and provide easier manipulation of "sets" than Blitz++
arrays. For example, you can easily expand and merge
sets which are stored in STL containers; doing this is
not so easy with Blitz++ arrays, which are designed for
numerical work.
STL containers are generally included by writing
#include <list> // for list<> #include <vector> // for vector<> #include <deque> // for deque<> #include <set> // for set<>
The []
operator is overloaded on arrays so that
the syntax array[container]
provides an indirect
view of the array. So far, this indirect view may
only be used as an lvalue (i.e. on the left-hand side
of an assignment statement).
The examples in the next sections are available
in the Blitz++ distribution in <examples/indirect.cpp>
.
2.11.1: Indirection using lists of array positions |
The simplest kind of indirection uses a list of points. For one-dimensional arrays, you can just use an STL container of integers. Example:
Array<int,1> A(5), B(5); A = 0; B = 1, 2, 3, 4, 5; vector<int> I; I.push_back(2); I.push_back(4); I.push_back(1); A[I] = B;
After this code, the array A contains
[ 0 2 3 0 5 ]
.
Note that arrays on the right-hand-side of the assignment must have the same shape as the array on the left-hand-side (before indirection). In the statement "A[I]=B", A and B must have the same shape, not I and B.
For multidimensional arrays, you can use an
STL container of TinyVector<int,N_rank>
objects. Example:
Array<int,2> A(4,4), B(4,4); A = 0; B = 10*tensor::i + tensor::j; typedef TinyVector<int,2> coord; list<coord> I; I.push_back(coord(1,1)); I.push_back(coord(2,2)); A[I] = B;
After this code, the array A contains:
0 0 0 0 0 11 0 0 0 0 22 0 0 0 0 0
(The tensor::i
notation is explained in the section on
index placeholders 2.8.5.
2.11.2: Cartesian-product indirection |
The Cartesian product of the sets I, J and K is the set of (i,j,k) tuples for which i is in I, j is in J, and k is in K.
Blitz++ implements cartesian-product indirection using an adaptor which takes a set of STL containers and iterates through their Cartesian product. Note that the cartesian product is never explicitly created. You create the Cartesian-product adaptor by calling the function:
template<class T_container>
indexSet(T_container& c1, T_container& c2, ...)
The returned adaptor can then be used in
the []
operator of an array object.
Here is a two-dimensional example:
Array<int,2> A(6,6), B(6,6); A = 0; B = 10*tensor::i + tensor::j; vector<int> I, J; I.push_back(1); I.push_back(2); I.push_back(4); J.push_back(0); J.push_back(2); J.push_back(5); A[indexSet(I,J)] = B;
After this code, the A array contains:
0 0 0 0 0 0 10 0 12 0 0 15 20 0 22 0 0 25 0 0 0 0 0 0 40 0 42 0 0 45 0 0 0 0 0 0
All the containers used in a cartesian product
must be the same type (e.g. all vector<int>
or
all set<TinyVector<int,2> >
, but they may
be different sizes. Singleton containers
(containers containing a single value) are fine.
2.11.3: Indirection with lists of strips |
You can also do indirection with a container of one-dimensional strips. This is useful when you want to manipulate some arbitrarily-shaped, well-connected subdomain of an array. By representing the subdomain as a list of strips, you allow Blitz++ to operate on vectors, rather than scattered points; this is much more efficient.
Strips are represented by objects of type
RectDomain<N>
, where N
is the
dimensionality of the array. The RectDomain<N>
class can be used to represent any rectangular
subdomain, but for indirection it is only
used to represent strips.
You create a strip by using this function:
RectDomain<N> strip(TinyVector<int,N> start,
int stripDimension, int ubound);
The start
parameter is where the strip
starts; stripDimension
is the dimension
in which the strip runs; ubound
is the
last index value for the strip. For
example, to create a 2-dimensional strip
from (2,5) to (2,9), one would write:
TinyVector<int,2> start(2,5);
RectDomain<2> myStrip = strip(start,secondDim,9);
Here is a more substantial example which creates a list of strips representing a circle subset of an array:
const int N = 7; Array<int,2> A(N,N), B(N,N); typedef TinyVector<int,2> coord; A = 0; B = 1; double centre_i = (N-1)/2.0; double centre_j = (N-1)/2.0; double radius = 0.8 * N/2.0; // circle will contain a list of strips which represent a circular // subdomain. list<RectDomain<2> > circle; for (int i=0; i < N; ++i) { double jdist2 = pow2(radius) - pow2(i-centre_i); if (jdist2 < 0.0) continue; int jdist = int(sqrt(jdist2)); coord startPos(i, int(centre_j - jdist)); circle.push_back(strip(startPos, secondDim, int(centre_j + jdist))); } // Set only those points in the circle subdomain to 1 A[circle] = B;
After this code, the A array contains:
0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 1 1 1 1 1 0 0 1 1 1 1 1 0 0 1 1 1 1 1 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0
2.12: Creating arrays of a user type |
Array
class with types you have created yourself,
or types from another library. If you want to do arithmetic on the
array, whatever operators you use on the arrays have to be defined
on the underlying type.
For example, here's a simple class for doing fixed point computations in the interval [0,1]:
#include <blitz/numinquire.h> // for huge() class FixedPoint { public: // The type to use for the mantissa typedef unsigned int T_mantissa; FixedPoint() { } explicit FixedPoint(T_mantissa mantissa) { mantissa_ = mantissa; } FixedPoint(double value) { assert((value >= 0.0) && (value <= 1.0)); mantissa_ = value * huge(T_mantissa()); } FixedPoint operator+(FixedPoint x) { return FixedPoint(mantissa_ + x.mantissa_); } double value() const { return mantissa_ / double(huge(T_mantissa())); } private: T_mantissa mantissa_; }; ostream& operator<<(ostream& os, const FixedPoint& a) { os << a.value(); return os; } |
The function huge(T)
returns the largest representable value
for type T; in the example above, it's equal to UINT_MAX
.
The FixedPoint
class declares three useful operations: conversion
from double
, addition, and outputing to an ostream
. We can
use all of these operations on an Array<FixedPoint>
object:
#include <blitz/array.h> using namespace blitz; int main() { // Create an array using the FixedPoint class: Array<FixedPoint, 2> A(4,4), B(4,4); A = 0.5, 0.3, 0.8, 0.2, 0.1, 0.3, 0.2, 0.9, 0.0, 1.0, 0.7, 0.4, 0.2, 0.3, 0.8, 0.4; B = A + 0.05; cout << "B = " << B << endl; return 0; } |
Note that the array A
is initialized using a comma-delimited list
of double
; this makes use of the constructor FixedPoint(double)
.
The assignment B = A + 0.05
uses
FixedPoint::operator+(FixedPoint)
, with an implicit conversion
from double
to FixedPoint
. Formatting the array B
onto
the standard output stream is done using the output operator
defined for FixedPoint
.
Here's the program output:
B = 4 x 4 0.55 0.35 0.85 0.25 0.15 0.35 0.25 0.95 0.05 0.05 0.75 0.45 0.25 0.35 0.85 0.45 |
2.13: Output formatting |
The current version of Blitz++ includes some rudimentary output formatting for one- and two-dimensional arrays (and array expressions). Here's an example:
#include <blitz/array.h> using namespace blitz; int main() { Array<int,2> A(4,5,FortranArray<2>()); firstIndex i; secondIndex j; A = 10*i + j; cout << "A = " << A << endl; Array<float,1> B(20); B = exp(-i/100.); cout << "B = " << endl << B << endl; return 0; } |
And the output:
A = 4 x 5 11 12 13 14 15 21 22 23 24 25 31 32 33 34 35 41 42 43 44 45 B = [ 1 0.99005 0.980199 0.970446 0.960789 0.951229 0.941765 0.932394 0.923116 0.913931 0.904837 0.895834 0.88692 0.878095 0.869358 0.860708 0.852144 0.843665 0.83527 0.826959 ] |
2.14: Array storage orders |
Blitz++ is very flexible about the way arrays are stored in memory. Starting indices can be 0, 1, or arbitrary numbers; arrays can be stored in row major, column major or an order based on any permutation of the dimensions; each dimension can be stored in either ascending or descending order. An N dimensional array can be stored in N!2N possible ways.
Before getting into the messy details, a review of array storage formats is useful. If you're already familiar with strides and bases, you might want to skip on to the next section.
2.14.1: Fortran and C-style arrays |
Suppose we want to store this two-dimensional array in memory:
[ 1 2 3 ] [ 4 5 6 ] [ 7 8 9 ]
Row major vs. column major
To lay the array out in memory, it's necessary to map the indices (i,j) into a one-dimensional block. Here are two ways the array might appear in memory:
[ 1 2 3 4 5 6 7 8 9 ] [ 1 4 7 2 5 8 3 6 9 ]
The first order corresponds to a C or C++ style array, and is called row-major ordering: the data is stored first by row, and then by column. The second order corresponds to a Fortran style array, and is called column-major ordering: the data is stored first by column, and then by row.
The simplest way of mapping the indices (i,j) into one-dimensional
memory is to take a linear combination. (Taking a linear
combination is sufficient for dense, asymmetric arrays, such as
are provided by the Blitz++ Array
class.) Here's the appropriate
linear combination for row major ordering:
memory offset = 3*i + 1*j
And for column major ordering:
memory offset = 1*i + 3*j
The coefficients of the (i,j) indices are called strides. For a row major storage of this array, the row stride is 3 -- you have to skip three memory locations to move down a row. The column stride is 1 -- you move one memory location to move to the next column. This is also known as unit stride. For column major ordering, the row and column strides are 1 and 3, respectively.
Bases
To throw another complication into this scheme, C-style arrays have indices which start at zero, and Fortran-style arrays have indices which start at one. The first valid index value is called the base. To account for a non-zero base, it's necessary to include an offset term in addition to the linear combination. Here's the mapping for a C-style array with i=0..3 and j=0..3:
memory offset = 0 + 3*i + 1*j
No offset is necessary since the indices start at zero for C-style arrays. For a Fortran-style array with i=1..4 and j=1..4, the mapping would be:
memory offset = -4 + 3*i + 1*j
By default, Blitz++ creates arrays in the C-style storage format (base zero, row major ordering). To create a Fortran-style array, you can use this syntax:
Array<int,2> A(3, 3, FortranArray<2>());
The third parameter, FortranArray<2>()
, tells the Array
constructor
to use a storage format appropriate for two-dimensional Fortran arrays
(base one, column major ordering).
A similar object, ColumnMajor<N>
, tells the Array
constructor
to use column major ordering, with base zero:
Array<int,2> B(3, 3, ColumnMajor<2>());
This creates a 3x3 array with indices i=0..2 and j=0..2.
In addition to supporting the 0 and 1 conventions for C and Fortran-style arrays, Blitz++ allows you to choose arbitrary bases, possibly different for each dimension. For example, this declaration creates an array whose indices have ranges i=5..8 and j=2..5:
Array<int,2> A(Range(5,8), Range(2,5));
2.14.2: Creating custom storage orders |
All Array
constructors take an optional parameter of type
GeneralArrayStorage<N_rank>
. This parameter encapsulates
a complete description of the storage format. If you want a
storage format other than C or Fortran-style, you have two
choices:
GeneralArrayStorage<N_rank>
,
customize the storage format, and use the object as a argument
for the Array
constructor.
GeneralArrayStorage<N_rank>
. This is useful if you
will be using the storage format many times. This approach
(inheriting from GeneralArrayStorage<N_rank>
) was used
to create the FortranArray<N_rank>
objects. If you want
to take this approach, you can use the declaration of
FortranArray<N_rank>
in <blitz/array.h>
as a guide.
The next sections describe how to modify a GeneralArrayStorage<N_rank>
object to suit your needs.
In higher dimensions
In more than two dimensions, the choice of storage order becomes more complicated. Suppose we had a 3x3x3 array. To map the indices (i,j,k) into memory, we might choose one of these mappings:
memory offset = 9*i + 3*j + 1*k memory offset = 1*i + 3*j + 9*k
The first corresponds to a C-style array, and the second to a Fortran-style array. But there are other choices; we can permute the strides (1,3,9) any which way:
memory offset = 1*i + 9*j + 3*k memory offset = 3*i + 1*j + 9*k memory offset = 3*i + 9*j + 1*k memory offset = 9*i + 1*j + 3*k
For an N dimensional array, there are N! such permutations. Blitz++
allows you to select any permutation of the dimensions as a storage
order. First you need to create an object of type
GeneralArrayStorage<N_rank>
:
GeneralArrayStorage<3> storage;
GeneralArrayStorage<N_rank>
contains a vector called ordering
which controls the order in which dimensions are stored in memory.
The ordering
vector will contain a permutation of the numbers
0, 1, ..., N_rank-1. Since some people are used to the first
dimension being 1 rather than 0, a set of symbols (firstDim, secondDim,
..., eleventhDim) are provided which make the code more legible.
The ordering
vector lists the dimensions in increasing order of
stride. You can access this vector using the member function
ordering()
. A C-style array, the default, would have:
storage.ordering() = thirdDim, secondDim, firstDim;
meaning that the third index (k) is associated with the smallest stride, and the first index (i) is associated with the largest stride. A Fortran-style array would have:
storage.ordering() = firstDim, secondDim, thirdDim;
Reversed dimensions
To add yet another wrinkle, there are some applications where the rows or columns need to be stored in reverse order. (For example, certain bitmap formats store image rows from bottom to top rather than top to bottom.)
Blitz++ allows you to store each dimension in either ascending or
descending order. By default, arrays are always stored in ascending
order. The GeneralArrayStorage<N_rank>
object contains a
vector called ascendingFlag
which indicates whether each
dimension is stored ascending (true
) or descending (false
).
To alter the contents of this vector, use the ascendingFlag()
method:
// Store the third dimension in descending order storage.ascendingFlag() = true, true, false; // Store all the dimensions in descending order storage.ascendingFlag() = false, false, false;
Setting the base vector
GeneralArrayStorage<N_rank>
also has a base
vector which
contains the base index value for each dimension. By default,
the base vector is set to zero. FortranArray<N_rank>
sets
the base vector to one.
To set your own set of bases, you have two choices:
base
vector inside the
GeneralArrayStorage<N_rank>
object. The method
base()
returns a mutable reference to the base
vector which you can use to set the bases.
Range
arguments to the
Array
constructor.
Here are some examples of the first approach:
// Set all bases equal to 5 storage.base() = 5; // Set the bases to [ 1 0 1 ] storage.base() = 1, 0, 1;
And of the second approach:
// Have bases of 5, but otherwise C-style storage Array<int,3> A(Range(5,7), Range(5,7), Range(5,7)); // Have bases of [ 1 0 1 ] and use a custom storage Array<int,3> B(Range(1,4), Range(0,3), Range(1,4), storage);
Working simultaneously with different storage orders
Once you have created an array object, you will probably never have to worry about its storage order. Blitz++ should handle arrays of different storage orders transparently. It's possible to mix arrays of different storage orders in one expression, and still get the correct result.
Note however, that mixing different storage orders in an expression may incur a performance penalty, since Blitz++ will have to pay more attention to differences in indexing than it normally would.
You may not mix arrays with different domains in the same expression. For example, adding a base zero to a base one array is a no-no. The reason for this restriction is that certain expressions become ambiguous, for example:
Array<int,1> A(Range(0,5)), B(Range(1,6)); A=0; B=0; using namespace blitz::tensor; int result = sum(A+B+i);
Should the index i
take its domain from array A
or array B
? To avoid such ambiguities, users are
forbidden from mixing arrays with different domains
in an expression.
Debug dumps of storage order information
In debug mode (-DBZ_DEBUG
), class Array
provides a member
function dumpStructureInformation()
which displays information
about the array storage:
Array<float,4> A(3,7,8,2,FortranArray<4>()); A.dumpStructureInformation(cerr);
The optional argument is an ostream
to dump information
to. It defaults to cout
. Here's the output:
Dump of Array<float, 4>: ordering_ = [ 0 1 2 3 ] ascendingFlag_ = [ 1111 ] base_ = [ 1 1 1 1 ] length_ = [ 3 7 8 2 ] stride_ = [ 1 3 21 168 ] zeroOffset_ = -193 numElements() = 336 storageContiguous = 1
A note about storage orders and initialization
When initializing arrays with comma delimited lists, note that the array is filled in storage order: from the first memory location to the last memory location. This won't cause any problems if you stick with C-style arrays, but it can be confusing for Fortran-style arrays:
Array<int,2> A(3, 3, FortranArray<2>()); A = 1, 2, 3, 4, 5, 6, 7, 8, 9; cout << A << endl;
The output from this code excerpt will be:
A = 3 x 3 1 4 7 2 5 8 3 6 9
This is because Fortran-style arrays are stored in column major order.
2.14.3: Storage orders example |
/***************************************************************************** * storage.cpp Blitz++ Array custom storage orders example * * $Id: storage.cpp,v 1.1 1997/07/16 19:38:23 tveldhui Exp $ * * $Log: storage.cpp,v $ * Revision 1.1 1997/07/16 19:38:23 tveldhui * Update: Alpha release 0.2 (Arrays) * ***************************************************************************** */ #include <blitz/array.h> BZ_USING_NAMESPACE(blitz) int main() { // 3x3 C-style row major storage, base zero Array<int,2> A(3, 3); // 3x3 column major storage, base zero Array<int,2> B(3, 3, ColumnMajorArray<2>()); // A custom storage format: // Indices have range 0..3, 0..3 // Column major ordering // Rows are stored ascending, columns stored descending GeneralArrayStorage<2> storage; storage.ordering() = firstRank, secondRank; storage.base() = 0, 0; storage.ascendingFlag() = true, false; Array<int,2> C(3, 3, storage); // Set each array equal to // [ 1 2 3 ] // [ 4 5 6 ] // [ 7 8 9 ] A = 1, 2, 3, 4, 5, 6, 7, 8, 9; cout << "A = " << A << endl; // Comma-delimited lists initialize in memory-storage order only. // Hence we list the values in column-major order to initialize B: B = 1, 4, 7, 2, 5, 8, 3, 6, 9; cout << "B = " << B << endl; // Array C is stored in column major, plus the columns are stored // in descending order! C = 3, 6, 9, 2, 5, 8, 1, 4, 7; cout << "C = " << C << endl; Array<int,2> D(3,3); D = A + B + C; #ifdef BZ_DEBUG A.dumpStructureInformation(); B.dumpStructureInformation(); C.dumpStructureInformation(); D.dumpStructureInformation(); #endif cout << "D = " << D << endl; return 0; } |
And the output:
A = 3 x 3 1 2 3 4 5 6 7 8 9 B = 3 x 3 1 2 3 4 5 6 7 8 9 C = 3 x 3 1 2 3 4 5 6 7 8 9 D = 3 x 3 3 6 9 12 15 18 21 24 27 |