2.4 Indexing, subarrays, and slicing

This section describes how to access the elements of an array. There are three main ways:

• Indexing obtains a single element
• Creating a subarray which refers to a smaller portion of an array
• Slicing to produce a smaller-dimensional view of a portion of an array

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 RectDomain and StridedDomain

The classes RectDomain and StridedDomain, defined in blitz/domain.h, offer a dimension-independent notation for subarrays.

RectDomain and StridedDomain can be thought of as a TinyVector<Range,N>. Both have a vector of lower- and upper-bounds; StridedDomain has a stride vector. For example, the subarray:

Array<int,2> B = A(Range(4,7), Range(8,11));  // 4..7, 8..11

could be obtained using RectDomain this way:

TinyVector<int,2> lowerBounds(4, 8);
TinyVector<int,2> upperBounds(7, 11);
RectDomain<2> subdomain(lowerBounds, upperBounds);

Array<int,2> B = A(subdomain);

Here are the prototypes of RectDomain and StridedDomain.

template<int N_rank>
class RectDomain {

public:
    RectDomain(const TinyVector<int,N_rank>& lbound,
        const TinyVector<int,N_rank>& ubound);

    const TinyVector<int,N_rank>& lbound() const;
    int lbound(int i) const;
    const TinyVector<int,N_rank>& ubound() const;
    int ubound(int i) const;
    Range operator[](int rank) const;
    void shrink(int amount);
    void shrink(int dim, int amount);
    void expand(int amount);
    void expand(int dim, int amount);
};

template<int N_rank>
class StridedDomain {

public:
    StridedDomain(const TinyVector<int,N_rank>& lbound,
        const TinyVector<int,N_rank>& ubound,
        const TinyVector<int,N_rank>& stride);

    const TinyVector<int,N_rank>& lbound() const;
    int lbound(int i) const;
    const TinyVector<int,N_rank>& ubound() const;
    int ubound(int i) const;
    const TinyVector<int,N_rank>& stride() const;
    int stride(int i) const;
    Range operator[](int rank) const;
    void shrink(int amount);
    void shrink(int dim, int amount);
    void expand(int amount);
    void expand(int dim, int amount);
};

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

Caution: 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.5 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:

#include <blitz/array.h>

using namespace blitz;

int main()
{
    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 ]

    return 0;
}

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

Note that if you use the same Range frequently, you can just construct one object and use it multiple times. For example:

Range all = Range::all();
A(0,all,all) = A(N-1,all,all);
A(all,0,all) = A(all,N-1,all);
A(all,all,0) = A(all,all,N-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:

A = (0,7) x (0,7)
[ 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.6 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 Array members).

Very important: whenever you have an assignment operator (=, +=, -=, etc.) the lhs must have the same shape as the rhs. If you want the array on the left hand side to be resized to the proper shape, you must do so by calling the resize method, for example:

A.resize(B.shape());    // Make A the same size as B
A = B;

2.4.7 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 = (0,5) x (0,5)
[ 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 ]