source: Sophya/trunk/SophyaExt/Blitz/blitz/array.h@ 2583

/***************************************************************************
 * blitz/array.h      Declaration of the Array<P_numtype, N_rank> class
 *
 * $Id: array.h,v 1.1.1.1 1999-04-09 17:59:01 ansari Exp $
 *
 * Copyright (C) 1997,1998 Todd Veldhuizen <tveldhui@seurat.uwaterloo.ca>
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * Suggestions:          blitz-suggest@cybervision.com
 * Bugs:                 blitz-bugs@cybervision.com
 *
 * For more information, please see the Blitz++ Home Page:
 *    http://seurat.uwaterloo.ca/blitz/
 *
 ***************************************************************************
 * $Log: not supported by cvs2svn $
 * Revision 1.2  1998/03/14 00:04:47  tveldhui
 * 0.2-alpha-05
 *
 * Revision 1.1  1997/07/16 14:51:20  tveldhui
 * Update: Alpha release 0.2 (Arrays)
 *
 */

/*
 * Wish list for array classes.
 *  - Arrays whose dimensions are unknown at compile time.
 *  - where()/elsewhere()/elsewhere() as in Dan Quinlan's implementation
 *  - block reduction operations
 *  - conversion to/from matrix & vector
 *  - apply(T func(T))
 *  - apply(T func(const T&))
 *  - apply<T func(T)>
 */

#ifndef BZ_ARRAY_H
#define BZ_ARRAY_H

#ifndef BZ_BLITZ_H
 #include <blitz/blitz.h>
#endif

#ifndef BZ_MEMBLOCK_H
 #include <blitz/memblock.h>
#endif

#ifndef BZ_RANGE_H
 #include <blitz/range.h>
#endif

#ifndef BZ_TINYVEC_H
 #include <blitz/tinyvec.h>
#endif

#ifndef BZ_TRAVERSAL_H
 #include <blitz/traversal.h>
#endif

#ifndef BZ_INDEXEXPR_H
 #include <blitz/indexexpr.h>
#endif

#ifndef BZ_PRETTYPRINT_H
 #include <blitz/prettyprint.h>
#endif

#include <blitz/array/slice.h>     // Subarrays and slicing
#include <blitz/array/map.h>       // Tensor index notation
#include <blitz/array/multi.h>     // Multicomponent arrays
#include <blitz/array/domain.h>    // RectDomain class

BZ_NAMESPACE(blitz)

/*
 * Forward declarations
 */

template<class T_numtype, int N_rank>
class ArrayIterator;

template<class P_expr>
class _bz_ArrayExpr;

#ifdef BZ_NEW_EXPRESSION_TEMPLATES
template<class T>
class ETBase { };
#endif

template<class T_array, class T_index>
class IndirectArray;

/*
 * Declaration of class GeneralStorage<N_rank>
 *
 * This class describes a storage format for an N-dimensional array.
 * The dimensions can be stored in an arbitrary order (for example, as
 * a C-style row major array or Fortran-style column major array, or
 * something else entirely).  Each dimension can be stored in either
 * ascending (the most common) or descending order.  Each dimension
 * can have its own base (starting index value: e.g. 0 for C-style arrays, 
 * 1 for Fortran arrays).
 *
 * GeneralArrayStorage<N> defaults to C-style arrays.  To implement
 * other storage formats, subclass and modify the constructor.  The
 * class FortranArray, below, is an example.
 *
 * Objects inheriting from GeneralArrayStorage<N> can be passed as
 * an optional constructor argument to Array objects.
 * e.g. Array<int,3> A(16,16,16, FortranArray<3>());
 * will create a 3-dimensional 16x16x16 Fortran-style array.
 */

template<int N_rank>
class GeneralArrayStorage {
public:
    class noInitializeFlag { };

    GeneralArrayStorage(noInitializeFlag)
    { }

    GeneralArrayStorage()
    {
        ordering_ = Range(N_rank - 1, 0, -1);
        ascendingFlag_ = 1;
        base_ = 0;
    }

    GeneralArrayStorage(const GeneralArrayStorage<N_rank>& x)
        : ordering_(x.ordering_), ascendingFlag_(x.ascendingFlag_),
          base_(x.base_)
    { 
    }

    ~GeneralArrayStorage()
    { }

    TinyVector<int, N_rank>& ordering()
    { return ordering_; }

    const TinyVector<int, N_rank>& ordering() const
    { return ordering_; }

    int ordering(int i) const
    { return ordering_[i]; }

    void setOrdering(int i, int order) 
    { ordering_[i] = order; }

    _bz_bool allRanksStoredAscending() const
    {
        _bz_bool result = _bz_true;
        for (int i=0; i < N_rank; ++i)
            result &= ascendingFlag_[i];
        return result;
    }

    _bz_bool isRankStoredAscending(int i) const
    { return ascendingFlag_[i]; }

    TinyVector<bool, N_rank>& ascendingFlag() 
    { return ascendingFlag_; }

    const TinyVector<bool, N_rank>& ascendingFlag() const
    { return ascendingFlag_; }

    void setAscendingFlag(int i, int ascendingFlag) 
    { ascendingFlag_[i] = ascendingFlag; }

    TinyVector<int, N_rank>& base()
    { return base_; }

    const TinyVector<int, N_rank>& base() const
    { return base_; }

    int base(int i) const
    { return base_[i]; }

    void setBase(int i, int base)
    { base_[i] = base; }

    void setBase(const TinyVector<int, N_rank>& base)
    { base_ = base; }

protected:
    /*
     * ordering_[] specifies the order in which the array is stored in
     * memory.  For a newly allocated array, ordering_(0) will give the
     * rank with unit stride, and ordering_(N_rank-1) will be the rank
     * with largest stride.  An order like [2, 1, 0] corresponds to
     * C-style array storage; an order like [0, 1, 2] corresponds to
     * Fortran array storage.
     *
     * ascendingFlag_[] indicates whether the data in a rank is stored
     * in ascending or descending order.  Most of the time these values
     * will all be true (indicating ascending order).  Some peculiar 
     * formats (e.g. MS-Windows BMP image format) store the data in 
     * descending order.
     *  
     * base_[] gives the first valid index for each rank.  For a C-style
     * array, all the base_ elements will be zero; for a Fortran-style
     * array, they will be one.  base_[] can be set arbitrarily using
     * the Array constructor which takes a Range argument, e.g.
     * Array<float,2> A(Range(30,40),Range(23,33));
     * will create an array with base_[] = { 30, 23 }.
     */
    TinyVector<int,  N_rank> ordering_;
    TinyVector<bool, N_rank> ascendingFlag_;
    TinyVector<int,  N_rank> base_;
};

/*
 * Class FortranArray specializes GeneralArrayStorage to provide Fortran
 * style arrays (column major ordering, base of 1).  The noInitializeFlag()
 * passed to the base constructor indicates that the subclass will take
 * care of initializing the ordering_, ascendingFlag_ and base_ members.
 */

template<int N_rank>
class FortranArray : public GeneralArrayStorage<N_rank> {
public:
    FortranArray()
        : GeneralArrayStorage<N_rank>(noInitializeFlag())
    {
        ordering_ = Range(0, N_rank - 1);
        ascendingFlag_ = 1;
        base_ = 1;
    }
};

/*
 * Class ColumnMajorArray specializes GeneralArrayStorage to provide column
 * major arrays (column major ordering, base of 0).
 */

template<int N_rank>
class ColumnMajorArray : public GeneralArrayStorage<N_rank> {
public:
    ColumnMajorArray()
        : GeneralArrayStorage<N_rank>(noInitializeFlag())
    {
        ordering_ = Range(0, N_rank - 1);
        ascendingFlag_ = 1;
        base_ = 0;
    }
};

/*
 * Declaration of class Array
 */

// NEEDS_WORK: Array should inherit protected from MemoryBlockReference.
// To make this work, need to expose MemoryBlockReference::numReferences()
// and make Array<P,N2> a friend of Array<P,N> for slicing.

template<class P_numtype, int N_rank>
class Array : public MemoryBlockReference<P_numtype> 
#ifdef BZ_NEW_EXPRESSION_TEMPLATES
    , public ETBase<Array<P_numtype,N_rank> >
#endif
{

public:
    //////////////////////////////////////////////
    // Public Types
    //////////////////////////////////////////////

    /*
     * T_numtype  is the numeric type stored in the array.
     * T_index    is a vector type which can be used to access elements
     *            of many-dimensional arrays.
     * T_array    is the array type itself -- Array<T_numtype, N_rank>
     * T_iterator is an iterator for the array.
     */

    typedef P_numtype                T_numtype;
    typedef TinyVector<int, N_rank>  T_index;
    typedef Array<T_numtype, N_rank> T_array;
    typedef ArrayIterator<T_numtype, N_rank> T_iterator;

    enum { _bz_rank = N_rank };

    //////////////////////////////////////////////
    // Constructors                             //
    //////////////////////////////////////////////

    
    /*
     * Construct an array from an array expression.
     */

    template<class T_expr>
    Array(_bz_ArrayExpr<T_expr> expr);

    /*
     * Any missing length arguments will have their value taken from the
     * last argument.  For example,
     *   Array<int,3> A(32,64);
     * will create a 32x64x64 array.  This is handled by setupStorage().
     */

    Array(GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
        : storage_(storage)
    {
        length_ = 0;
        stride_ = 0;
        zeroOffset_ = 0;
        storageContiguous_ = false;
    }

    _bz_explicit Array(int length0, 
        GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
        : storage_(storage)
    {
        length_[0] = length0;
        setupStorage(0);
    }

    Array(int length0, int length1,
        GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
        : storage_(storage)
    {
        BZPRECONDITION(N_rank >= 2);
        TAU_TYPE_STRING(p1, "Array<T,N>::Array() [T="
            + CT(T_numtype) + ",N=" + CT(N_rank) + "]");
        TAU_PROFILE(p1, "void (int,int)", TAU_BLITZ);

        length_[0] = length0;
        length_[1] = length1;
        setupStorage(1);
    }

    Array(int length0, int length1, int length2,
        GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
        : storage_(storage)
    {
        BZPRECONDITION(N_rank >= 3);
        length_[0] = length0;
        length_[1] = length1;
        length_[2] = length2;
        setupStorage(2);
    }

    Array(int length0, int length1, int length2, int length3,
        GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
        : storage_(storage)
    {
        BZPRECONDITION(N_rank >= 4);
        length_[0] = length0;
        length_[1] = length1;
        length_[2] = length2;
        length_[3] = length3;
        setupStorage(3);
    }

    Array(int length0, int length1, int length2, int length3, int length4,
        GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
        : storage_(storage)
    {
        BZPRECONDITION(N_rank >= 5);
        length_[0] = length0;
        length_[1] = length1;
        length_[2] = length2;
        length_[3] = length3;
        length_[4] = length4;
        setupStorage(4);
    }

    Array(int length0, int length1, int length2, int length3, int length4,
        int length5,
        GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
        : storage_(storage)
    {
        BZPRECONDITION(N_rank >= 6);
        length_[0] = length0;
        length_[1] = length1;
        length_[2] = length2;
        length_[3] = length3;
        length_[4] = length4;
        length_[5] = length5;
        setupStorage(5);
    }

    Array(int length0, int length1, int length2, int length3, int length4,
        int length5, int length6,
        GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
        : storage_(storage)
    {
        BZPRECONDITION(N_rank >= 7);
        length_[0] = length0;
        length_[1] = length1;
        length_[2] = length2;
        length_[3] = length3;
        length_[4] = length4;
        length_[5] = length5;
        length_[6] = length6;
        setupStorage(6);
    }

    Array(int length0, int length1, int length2, int length3, int length4,
        int length5, int length6, int length7,
        GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
        : storage_(storage)
    {
        BZPRECONDITION(N_rank >= 8);
        length_[0] = length0;
        length_[1] = length1;
        length_[2] = length2;
        length_[3] = length3;
        length_[4] = length4;
        length_[5] = length5;
        length_[6] = length6;
        length_[7] = length7;
        setupStorage(7);
    }

    Array(int length0, int length1, int length2, int length3, int length4,
        int length5, int length6, int length7, int length8,
        GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
        : storage_(storage)
    {
        BZPRECONDITION(N_rank >= 9);
        length_[0] = length0;
        length_[1] = length1;
        length_[2] = length2;
        length_[3] = length3;
        length_[4] = length4;
        length_[5] = length5;
        length_[6] = length6;
        length_[7] = length7;
        length_[8] = length8;
        setupStorage(8);
    }

    Array(int length0, int length1, int length2, int length3, int length4,
        int length5, int length6, int length7, int length8, int length9,
        GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
        : storage_(storage)
    {
        BZPRECONDITION(N_rank >= 10);
        length_[0] = length0;
        length_[1] = length1;
        length_[2] = length2;
        length_[3] = length3;
        length_[4] = length4;
        length_[5] = length5;
        length_[6] = length6;
        length_[7] = length7;
        length_[8] = length8;
        length_[9] = length9;
        setupStorage(9);
    }

    Array(int length0, int length1, int length2, int length3, int length4,
        int length5, int length6, int length7, int length8, int length9,
        int length10,
        GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
        : storage_(storage)
    {
        BZPRECONDITION(N_rank >= 11);
        length_[0] = length0;
        length_[1] = length1;
        length_[2] = length2;
        length_[3] = length3;
        length_[4] = length4;
        length_[5] = length5;
        length_[6] = length6;
        length_[7] = length7;
        length_[8] = length8;
        length_[9] = length9;
        length_[10] = length10;
        setupStorage(10);
    }

    /*
     * Construct an array from an existing block of memory.  Ownership
     * is not acquired (i.e. the memory block will not be freed by Blitz++).
     */
    Array(T_numtype* _bz_restrict dataFirst, TinyVector<int, N_rank> shape,
        GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
      : MemoryBlockReference<T_numtype>(product(shape), dataFirst),
        storage_(storage)
    {
        BZPRECONDITION(dataFirst != 0);

        length_ = shape;
        computeStrides();
        data_ += zeroOffset_;
        storageContiguous_ = _bz_true;
    }

    /*
     * Construct an array from an existing block of memory, with a
     * given set of strides.  Ownership is not acquired (i.e. the memory
     * block will not be freed by Blitz++).
     */
    Array(T_numtype* _bz_restrict dataFirst, TinyVector<int, N_rank> shape,
        TinyVector<int, N_rank> stride, 
        GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
      : MemoryBlockReference<T_numtype>(product(shape), dataFirst),
        storage_(storage)
    {
        BZPRECONDITION(dataFirst != 0);

        length_ = shape;
        stride_ = stride;
        calculateZeroOffset();
        data_ += zeroOffset_;

        // NEEDS_WORK: the setting of storageContiguous should be determined
        // by inspecting the strides.
        storageContiguous_ = _bz_false;   
    }

    /*
     * This constructor takes an extent (length) vector and storage format.
     */

    Array(const TinyVector<int, N_rank>& extent, 
        GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
        : storage_(storage)
    {
        length_ = extent;
        setupStorage(N_rank - 1);
    }

    /*
     * This construct takes a vector of bases (lbounds) and a vector of
     * extents.
     */

    Array(const TinyVector<int, N_rank>& lbounds,
        const TinyVector<int, N_rank>& extent,
        const GeneralArrayStorage<N_rank>& storage);

    /*
     * These constructors allow arbitrary bases (starting indices) to be set.
     * e.g. Array<int,2> A(Range(10,20), Range(20,30))
     * will create an 11x11 array whose indices are 10..20 and 20..30
     */
    Array(Range r0, 
        GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
        : storage_(storage)
    {
        BZPRECONDITION(r0.isAscendingContiguous());

        length_[0] = r0.length();
        storage_.setBase(0, r0.first());
        setupStorage(0);
    }

    Array(Range r0, Range r1,
        GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
        : storage_(storage)
    {
        BZPRECONDITION(r0.isAscendingContiguous() && 
            r1.isAscendingContiguous());

        length_[0] = r0.length();
        storage_.setBase(0, r0.first());
        length_[1] = r1.length();
        storage_.setBase(1, r1.first());

        setupStorage(1);
    }

    Array(Range r0, Range r1, Range r2,
        GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
        : storage_(storage)
    {
        BZPRECONDITION(r0.isAscendingContiguous() &&
            r1.isAscendingContiguous() && r2.isAscendingContiguous());

        length_[0] = r0.length();
        storage_.setBase(0, r0.first());
        length_[1] = r1.length();
        storage_.setBase(1, r1.first());
        length_[2] = r2.length();
        storage_.setBase(2, r2.first());

        setupStorage(2);
    }

    Array(Range r0, Range r1, Range r2, Range r3,
        GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
        : storage_(storage)
    {
        BZPRECONDITION(r0.isAscendingContiguous() &&
            r1.isAscendingContiguous() && r2.isAscendingContiguous()
            && r3.isAscendingContiguous());

        length_[0] = r0.length();
        storage_.setBase(0, r0.first());
        length_[1] = r1.length();
        storage_.setBase(1, r1.first());
        length_[2] = r2.length();
        storage_.setBase(2, r2.first());
        length_[3] = r3.length();
        storage_.setBase(3, r3.first());

        setupStorage(3);
    }

    Array(Range r0, Range r1, Range r2, Range r3, Range r4,
        GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
        : storage_(storage)
    {
        BZPRECONDITION(r0.isAscendingContiguous() &&
            r1.isAscendingContiguous() && r2.isAscendingContiguous()
            && r3.isAscendingContiguous() && r4.isAscendingContiguous());

        length_[0] = r0.length();
        storage_.setBase(0, r0.first());
        length_[1] = r1.length();
        storage_.setBase(1, r1.first());
        length_[2] = r2.length();
        storage_.setBase(2, r2.first());
        length_[3] = r3.length();
        storage_.setBase(3, r3.first());
        length_[4] = r4.length();
        storage_.setBase(4, r4.first());

        setupStorage(4);
    }

    Array(Range r0, Range r1, Range r2, Range r3, Range r4, Range r5,
        GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
        : storage_(storage)
    {
        BZPRECONDITION(r0.isAscendingContiguous() &&
            r1.isAscendingContiguous() && r2.isAscendingContiguous()
            && r3.isAscendingContiguous() && r4.isAscendingContiguous()
            && r5.isAscendingContiguous());

        length_[0] = r0.length();
        storage_.setBase(0, r0.first());
        length_[1] = r1.length();
        storage_.setBase(1, r1.first());
        length_[2] = r2.length();
        storage_.setBase(2, r2.first());
        length_[3] = r3.length();
        storage_.setBase(3, r3.first());
        length_[4] = r4.length();
        storage_.setBase(4, r4.first());
        length_[5] = r5.length();
        storage_.setBase(5, r5.first());

        setupStorage(5);
    }

    Array(Range r0, Range r1, Range r2, Range r3, Range r4, Range r5,
        Range r6,
        GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
        : storage_(storage)
    {
        BZPRECONDITION(r0.isAscendingContiguous() &&
            r1.isAscendingContiguous() && r2.isAscendingContiguous()
            && r3.isAscendingContiguous() && r4.isAscendingContiguous()
            && r5.isAscendingContiguous() && r6.isAscendingContiguous());

        length_[0] = r0.length();
        storage_.setBase(0, r0.first());
        length_[1] = r1.length();
        storage_.setBase(1, r1.first());
        length_[2] = r2.length();
        storage_.setBase(2, r2.first());
        length_[3] = r3.length();
        storage_.setBase(3, r3.first());
        length_[4] = r4.length();
        storage_.setBase(4, r4.first());
        length_[5] = r5.length();
        storage_.setBase(5, r5.first());
        length_[6] = r6.length();
        storage_.setBase(6, r6.first());

        setupStorage(6);
    }

    Array(Range r0, Range r1, Range r2, Range r3, Range r4, Range r5,
        Range r6, Range r7,
        GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
        : storage_(storage)
    {
        BZPRECONDITION(r0.isAscendingContiguous() &&
            r1.isAscendingContiguous() && r2.isAscendingContiguous()
            && r3.isAscendingContiguous() && r4.isAscendingContiguous()
            && r5.isAscendingContiguous() && r6.isAscendingContiguous()
            && r7.isAscendingContiguous());

        length_[0] = r0.length();
        storage_.setBase(0, r0.first());
        length_[1] = r1.length();
        storage_.setBase(1, r1.first());
        length_[2] = r2.length();
        storage_.setBase(2, r2.first());
        length_[3] = r3.length();
        storage_.setBase(3, r3.first());
        length_[4] = r4.length();
        storage_.setBase(4, r4.first());
        length_[5] = r5.length();
        storage_.setBase(5, r5.first());
        length_[6] = r6.length();
        storage_.setBase(6, r6.first());
        length_[7] = r7.length();
        storage_.setBase(7, r7.first());

        setupStorage(7);
    }

    Array(Range r0, Range r1, Range r2, Range r3, Range r4, Range r5,
        Range r6, Range r7, Range r8,
        GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
        : storage_(storage)
    {
        BZPRECONDITION(r0.isAscendingContiguous() &&
            r1.isAscendingContiguous() && r2.isAscendingContiguous()
            && r3.isAscendingContiguous() && r4.isAscendingContiguous()
            && r5.isAscendingContiguous() && r6.isAscendingContiguous()
            && r7.isAscendingContiguous() && r8.isAscendingContiguous());

        length_[0] = r0.length();
        storage_.setBase(0, r0.first());
        length_[1] = r1.length();
        storage_.setBase(1, r1.first());
        length_[2] = r2.length();
        storage_.setBase(2, r2.first());
        length_[3] = r3.length();
        storage_.setBase(3, r3.first());
        length_[4] = r4.length();
        storage_.setBase(4, r4.first());
        length_[5] = r5.length();
        storage_.setBase(5, r5.first());
        length_[6] = r6.length();
        storage_.setBase(6, r6.first());
        length_[7] = r7.length();
        storage_.setBase(7, r7.first());
        length_[8] = r8.length();
        storage_.setBase(8, r8.first());

        setupStorage(8);
    }

    Array(Range r0, Range r1, Range r2, Range r3, Range r4, Range r5,
        Range r6, Range r7, Range r8, Range r9,
        GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
        : storage_(storage)
    {
        BZPRECONDITION(r0.isAscendingContiguous() &&
            r1.isAscendingContiguous() && r2.isAscendingContiguous()
            && r3.isAscendingContiguous() && r4.isAscendingContiguous()
            && r5.isAscendingContiguous() && r6.isAscendingContiguous()
            && r7.isAscendingContiguous() && r8.isAscendingContiguous()
            && r9.isAscendingContiguous());

        length_[0] = r0.length();
        storage_.setBase(0, r0.first());
        length_[1] = r1.length();
        storage_.setBase(1, r1.first());
        length_[2] = r2.length();
        storage_.setBase(2, r2.first());
        length_[3] = r3.length();
        storage_.setBase(3, r3.first());
        length_[4] = r4.length();
        storage_.setBase(4, r4.first());
        length_[5] = r5.length();
        storage_.setBase(5, r5.first());
        length_[6] = r6.length();
        storage_.setBase(6, r6.first());
        length_[7] = r7.length();
        storage_.setBase(7, r7.first());
        length_[8] = r8.length();
        storage_.setBase(8, r8.first());
        length_[9] = r9.length();
        storage_.setBase(9, r9.first());

        setupStorage(9);
    }

    Array(Range r0, Range r1, Range r2, Range r3, Range r4, Range r5,
        Range r6, Range r7, Range r8, Range r9, Range r10,
        GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
        : storage_(storage)
    {
        BZPRECONDITION(r0.isAscendingContiguous() &&
            r1.isAscendingContiguous() && r2.isAscendingContiguous()
            && r3.isAscendingContiguous() && r4.isAscendingContiguous()
            && r5.isAscendingContiguous() && r6.isAscendingContiguous()
            && r7.isAscendingContiguous() && r8.isAscendingContiguous()
            && r9.isAscendingContiguous() && r10.isAscendingContiguous());

        length_[0] = r0.length();
        storage_.setBase(0, r0.first());
        length_[1] = r1.length();
        storage_.setBase(1, r1.first());
        length_[2] = r2.length();
        storage_.setBase(2, r2.first());
        length_[3] = r3.length();
        storage_.setBase(3, r3.first());
        length_[4] = r4.length();
        storage_.setBase(4, r4.first());
        length_[5] = r5.length();
        storage_.setBase(5, r5.first());
        length_[6] = r6.length();
        storage_.setBase(6, r6.first());
        length_[7] = r7.length();
        storage_.setBase(7, r7.first());
        length_[8] = r8.length();
        storage_.setBase(8, r8.first());
        length_[9] = r9.length();
        storage_.setBase(9, r9.first());
        length_[10] = r10.length();
        storage_.setBase(10, r10.first());

        setupStorage(10);
    }

    /*
     * Create a reference of another array
     */
    Array(const Array<T_numtype, N_rank>& array)
    {
        // NEEDS_WORK: this const_cast is a tad ugly.
        reference(const_cast<T_array&>(array));
    }

    /*
     * These constructors are used for creating interlaced arrays (see
     * <blitz/arrayshape.h>
     */
    Array(const TinyVector<int,N_rank-1>& shape,
        int lastExtent, const GeneralArrayStorage<N_rank>& storage);
    //Array(const TinyVector<Range,N_rank-1>& shape,
    //    int lastExtent, const GeneralArrayStorage<N_rank>& storage);

    /*
     * These constructors make the array a view of a subportion of another
     * array.  If there fewer than N_rank Range arguments provided, no
     * slicing is performed in the unspecified ranks.
     * e.g. Array<int,3> A(20,20,20);
     *      Array<int,3> B(A, Range(5,15));
     * is equivalent to:
     *      Array<int,3> B(A, Range(5,15), Range::all(), Range::all());
     */
    Array(Array<T_numtype, N_rank>& array, Range r0)
    {
        constructSubarray(array, r0);
    }

    Array(Array<T_numtype, N_rank>& array, Range r0, Range r1)
    {
        constructSubarray(array, r0, r1);
    }

    Array(Array<T_numtype, N_rank>& array, Range r0, Range r1, Range r2)
    {
        constructSubarray(array, r0, r1, r2);
    }

    Array(Array<T_numtype, N_rank>& array, Range r0, Range r1, Range r2,
        Range r3)
    {
        constructSubarray(array, r0, r1, r2, r3);
    }

    Array(Array<T_numtype, N_rank>& array, Range r0, Range r1, Range r2,
        Range r3, Range r4)
    {
        constructSubarray(array, r0, r1, r2, r3, r4);
    }

    Array(Array<T_numtype, N_rank>& array, Range r0, Range r1, Range r2,
        Range r3, Range r4, Range r5)
    {
        constructSubarray(array, r0, r1, r2, r3, r4, r5);
    }

    Array(Array<T_numtype, N_rank>& array, Range r0, Range r1, Range r2,
        Range r3, Range r4, Range r5, Range r6)
    {
        constructSubarray(array, r0, r1, r2, r3, r4, r5, r6);
    }

    Array(Array<T_numtype, N_rank>& array, Range r0, Range r1, Range r2,
        Range r3, Range r4, Range r5, Range r6, Range r7)
    {
        constructSubarray(array, r0, r1, r2, r3, r4, r5, r6, r7);
    }

    Array(Array<T_numtype, N_rank>& array, Range r0, Range r1, Range r2,
        Range r3, Range r4, Range r5, Range r6, Range r7, Range r8)
    {
        constructSubarray(array, r0, r1, r2, r3, r4, r5, r6, r7, r8);
    }

    Array(Array<T_numtype, N_rank>& array, Range r0, Range r1, Range r2,
        Range r3, Range r4, Range r5, Range r6, Range r7, Range r8, Range r9)
    {
        constructSubarray(array, r0, r1, r2, r3, r4, r5, r6, r7, r8, r9);
    }

    Array(Array<T_numtype, N_rank>& array, Range r0, Range r1, Range r2,
        Range r3, Range r4, Range r5, Range r6, Range r7, Range r8, Range r9,
        Range r10)
    {
        constructSubarray(array, r0, r1, r2, r3, r4, r5, r6, r7, r8, r9, r10);
    }

    Array(Array<T_numtype, N_rank>& array, 
        const RectDomain<N_rank>& subdomain)
    {
        constructSubarray(array, subdomain);
    }

    /*
     * This constructor is invoked by the operator()'s which take
     * a combination of integer and Range arguments.  It's not intended
     * for end-user use.
     */
    template<int N_rank2, class R0, class R1, class R2, class R3, class R4,
        class R5, class R6, class R7, class R8, class R9, class R10>
    Array(Array<T_numtype,N_rank2>& array, R0 r0, R1 r1, R2 r2,
        R3 r3, R4 r4, R5 r5, R6 r6, R7 r7, R8 r8, R9 r9, R10 r10)
    {
        constructSlice(array, r0, r1, r2, r3, r4, r5, r6, r7, r8, r9, r10);
    }

    //////////////////////////////////////////////
    // Member functions
    //////////////////////////////////////////////

    const TinyVector<int, N_rank>&    base() const
    { return storage_.base(); }

    int                               base(int rank) const
    { return storage_.base(rank); }

    T_iterator                        begin() const
    { return T_iterator(*this); }

    // Deprecated: now extractComponent(...)
    template<class T_numtype2>
    Array<T_numtype2,N_rank>          chopComponent(T_numtype2 a, int compNum,
                                          int numComponents)
    { return extractComponent(a, compNum, numComponents); }

    int                               cols() const
    { return length_[1]; }

    int                               columns() const
    { return length_[1]; }

    T_array                           copy() const;

    // data_ always refers to the point (0,0,...,0) which may
    // not be in the array if the base is not zero in each rank.
    // These data() routines return a pointer to the first
    // element in the array (but note that it may not be
    // stored first in memory if some ranks are stored descending).

    const T_numtype* _bz_restrict     data() const
    { return data_ + dot(storage_.base(), stride_); }

    T_numtype* _bz_restrict           data() 
    { return data_ + dot(storage_.base(), stride_); }

    // These dataZero() routines refer to the point (0,0,...,0)
    // which may not be in the array if the bases are nonzero.
    
    const T_numtype* _bz_restrict     dataZero() const
    { return data_; }

    T_numtype* _bz_restrict           dataZero()
    { return data_; }

    // These dataFirst() routines refer to the element in the
    // array which falls first in memory.
    
    const T_numtype* _bz_restrict     dataFirst() const
    {
        return data_ + dot(storage_.base() 
            + (1 - storage_.ascendingFlag()) * (length_ - 1), stride_);
    }

    T_numtype* _bz_restrict           dataFirst()
    {
        return data_ + dot(storage_.base()
            + (1 - storage_.ascendingFlag()) * (length_ - 1), stride_);
    }

    int                               depth() const
    { return length_[2]; }

    int                               dimensions() const
    { return N_rank; }

    RectDomain<N_rank>                domain() const
    {
        return RectDomain<N_rank>(lbound(), ubound());
    }

    void dumpStructureInformation(ostream& os = cout) const;

    // NEEDS_WORK -- end()

    int                               extent(int rank) const
    { return length_[rank]; }

    const TinyVector<int,N_rank>&     extent() const
    { return length_; }

    template<class T_numtype2>
    Array<T_numtype2,N_rank>          extractComponent(T_numtype2, int compNum,
                                          int numComponents);

    void                              free() 
    {
        changeToNullBlock();
        length_ = 0;
    }
 
    _bz_bool                          isMajorRank(int rank) const
    { return storage_.ordering(rank) == 0; }

    _bz_bool                          isMinorRank(int rank) const
    { return storage_.ordering(rank) != 0; }

    _bz_bool                          isRankStoredAscending(int rank) const
    { return storage_.isRankStoredAscending(rank); }

    _bz_bool                          isStorageContiguous() const
    { return storageContiguous_; }

    int                               lbound(int rank) const
    { return base(rank); }

    TinyVector<int,N_rank>            lbound() const
    {
        return base();
    }

    int                               length(int rank) const
    { return length_[rank]; }

    const TinyVector<int, N_rank>&    length() const
    { return length_; }

    void                              makeUnique();

    int                               numElements() const
    { return product(length_); }

    // NEEDS_WORK -- Expose the numReferences() method
    // MemoryBlockReference<T_numtype>::numReferences;

    // The storage_.ordering_ array is a list of dimensions from
    // the most minor (stride 1) to major dimension.  Generally,
    // ordering(0) will return the dimension which has the smallest
    // stride, and ordering(N_rank-1) will return the dimension with
    // the largest stride.
    int                               ordering(int storageRankIndex) const
    { return storage_.ordering(storageRankIndex); }

    const TinyVector<int, N_rank>&    ordering() const
    { return storage_.ordering(); }

    void                              transposeSelf(int r0, int r1, int r2=0, 
        int r3=0, int r4=0, int r5=0, int r6=0, int r7=0, int r8=0, int 
        r9=0, int r10=0);
    T_array                           transpose(int r0, int r1, int r2=0,
        int r3=0, int r4=0, int r5=0, int r6=0, int r7=0, int r8=0, int
        r9=0, int r10=0);

    int                               rank() const
    { return N_rank; }

    void                              reference(T_array&);

    void                              resize(int extent);
    void                              resize(int extent1, int extent2);
    void                              resize(int extent1, int extent2,
                                        int extent3);
    void                              resize(int extent1, int extent2,
                                        int extent3, int extent4);
    void                              resize(int extent1, int extent2,
                                        int extent3, int extent4, int extent5);
    void                              resize(int extent1, int extent2,
                                        int extent3, int extent4, int extent5,
                                        int extent6);
    void                              resize(int extent1, int extent2,
                                        int extent3, int extent4, int extent5,
                                        int extent6, int extent7);
    void                              resize(int extent1, int extent2,
                                        int extent3, int extent4, int extent5,
                                        int extent6, int extent7, int extent8);
    void                              resize(int extent1, int extent2,
                                        int extent3, int extent4, int extent5,
                                        int extent6, int extent7, int extent8,
                                        int extent9);
    void                              resize(int extent1, int extent2,
                                        int extent3, int extent4, int extent5,
                                        int extent6, int extent7, int extent8,
                                        int extent9, int extent10);
    void                              resize(int extent1, int extent2,
                                        int extent3, int extent4, int extent5,
                                        int extent6, int extent7, int extent8,
                                        int extent9, int extent10, 
                                        int extent11);


    void                              resize(Range r1);
    void                              resize(Range r1, Range r2);
    void                              resize(Range r1, Range r2, Range r3);
    void                              resize(Range r1, Range r2, Range r3,
                                        Range r4);
    void                              resize(Range r1, Range r2, Range r3,
                                        Range r4, Range r5);
    void                              resize(Range r1, Range r2, Range r3,
                                        Range r4, Range r5, Range r6);
    void                              resize(Range r1, Range r2, Range r3,
                                        Range r4, Range r5, Range r6,
                                        Range r7);
    void                              resize(Range r1, Range r2, Range r3,
                                        Range r4, Range r5, Range r6,
                                        Range r7, Range r8);
    void                              resize(Range r1, Range r2, Range r3,
                                        Range r4, Range r5, Range r6,
                                        Range r7, Range r8, Range r9);
    void                              resize(Range r1, Range r2, Range r3,
                                        Range r4, Range r5, Range r6,
                                        Range r7, Range r8, Range r9, 
                                        Range r10);
    void                              resize(Range r1, Range r2, Range r3,
                                        Range r4, Range r5, Range r6,
                                        Range r7, Range r8, Range r9, 
                                        Range r10, Range r11);

    void                              resize(const TinyVector<int,N_rank>&);
 

    void                              resizeAndPreserve(int extent);
    void                              resizeAndPreserve(int extent1, 
                                        int extent2);
    void                              resizeAndPreserve(int extent1, 
                                        int extent2, int extent3);
    void                              resizeAndPreserve(int extent1,
                                        int extent2, int extent3, int extent4);
    void                              resizeAndPreserve(int extent1,
                                        int extent2, int extent3, int extent4,
                                        int extent5);
    void                              resizeAndPreserve(int extent1,
                                        int extent2, int extent3, int extent4,
                                        int extent5, int extent6);
    void                              resizeAndPreserve(int extent1,
                                        int extent2, int extent3, int extent4,
                                        int extent5, int extent6, int extent7);
    void                              resizeAndPreserve(int extent1,
                                        int extent2, int extent3, int extent4,
                                        int extent5, int extent6, int extent7,
                                        int extent8);
    void                              resizeAndPreserve(int extent1,
                                        int extent2, int extent3, int extent4,
                                        int extent5, int extent6, int extent7,
                                        int extent8, int extent9);
    void                              resizeAndPreserve(int extent1,
                                        int extent2, int extent3, int extent4,
                                        int extent5, int extent6, int extent7,
                                        int extent8, int extent9, 
                                        int extent10);
    void                              resizeAndPreserve(int extent1,
                                        int extent2, int extent3, int extent4,
                                        int extent5, int extent6, int extent7,
                                        int extent8, int extent9, int extent10,
                                        int extent11);

    // NEEDS_WORK -- resizeAndPreserve(Range,...)
    // NEEDS_WORK -- resizeAndPreserve(const TinyVector<int,N_rank>&);
    // NEEDS_WORK -- resizeAndPreserve(const Domain<N_rank>&);

    T_array                           reverse(int rank);
    void                              reverseSelf(int rank);

    int                               rows() const
    { return length_[0]; }

    void                              slice(int rank, Range r);

    const TinyVector<int, N_rank>&    shape() const
    { return length_; }

    int                               size() const
    { return numElements(); }

    const TinyVector<int, N_rank>&    stride() const
    { return stride_; }

    int                               stride(int rank) const
    { return stride_[rank]; }

    int                               ubound(int rank) const
    { return base(rank) + length_(rank) - 1; }

    TinyVector<int, N_rank>           ubound() const
    { 
        TinyVector<int, N_rank> ub;
        ub = base() + extent() - 1;
        return ub;
    }

    int                               zeroOffset() const
    { return zeroOffset_; }

    //////////////////////////////////////////////
    // Debugging routines
    //////////////////////////////////////////////

    _bz_bool isInRange(int i0) const
    {
        return unsigned(i0 - base(0)) < length_[0];
    }

    _bz_bool isInRange(int i0, int i1) const
    {
        return unsigned(i0 - base(0)) < length_[0]
            && unsigned(i1 - base(1)) < length_[1];
    }

    _bz_bool isInRange(int i0, int i1, int i2) const
    {
        return unsigned(i0 - base(0)) < length_[0]
            && unsigned(i1 - base(1)) < length_[1]
            && unsigned(i2 - base(2)) < length_[2];
    }

    _bz_bool isInRange(int i0, int i1, int i2, int i3) const
    {
        return unsigned(i0 - base(0)) < length_[0]
            && unsigned(i1 - base(1)) < length_[1]
            && unsigned(i2 - base(2)) < length_[2]
            && unsigned(i3 - base(3)) < length_[3];
    }

    _bz_bool isInRange(int i0, int i1, int i2, int i3, int i4) const
    {
        return unsigned(i0 - base(0)) < length_[0]
            && unsigned(i1 - base(1)) < length_[1]
            && unsigned(i2 - base(2)) < length_[2]
            && unsigned(i3 - base(3)) < length_[3]
            && unsigned(i4 - base(4)) < length_[4];
    }

    _bz_bool isInRange(int i0, int i1, int i2, int i3, int i4,
        int i5) const
    {
        return unsigned(i0 - base(0)) < length_[0]
            && unsigned(i1 - base(1)) < length_[1]
            && unsigned(i2 - base(2)) < length_[2]
            && unsigned(i3 - base(3)) < length_[3]
            && unsigned(i4 - base(4)) < length_[4]
            && unsigned(i5 - base(5)) < length_[5];
    }

    _bz_bool isInRange(int i0, int i1, int i2, int i3, int i4,
        int i5, int i6) const
    {
        return unsigned(i0 - base(0)) < length_[0]
            && unsigned(i1 - base(1)) < length_[1]
            && unsigned(i2 - base(2)) < length_[2]
            && unsigned(i3 - base(3)) < length_[3]
            && unsigned(i4 - base(4)) < length_[4]
            && unsigned(i5 - base(5)) < length_[5]
            && unsigned(i6 - base(6)) < length_[6];
    }

    _bz_bool isInRange(int i0, int i1, int i2, int i3, int i4,
        int i5, int i6, int i7) const
    {
        return unsigned(i0 - base(0)) < length_[0]
            && unsigned(i1 - base(1)) < length_[1]
            && unsigned(i2 - base(2)) < length_[2]
            && unsigned(i3 - base(3)) < length_[3]
            && unsigned(i4 - base(4)) < length_[4]
            && unsigned(i5 - base(5)) < length_[5]
            && unsigned(i6 - base(6)) < length_[6]
            && unsigned(i7 - base(7)) < length_[7];
    }

    _bz_bool isInRange(int i0, int i1, int i2, int i3, int i4,
        int i5, int i6, int i7, int i8) const
    {
        return unsigned(i0 - base(0)) < length_[0]
            && unsigned(i1 - base(1)) < length_[1]
            && unsigned(i2 - base(2)) < length_[2]
            && unsigned(i3 - base(3)) < length_[3]
            && unsigned(i4 - base(4)) < length_[4]
            && unsigned(i5 - base(5)) < length_[5]
            && unsigned(i6 - base(6)) < length_[6]
            && unsigned(i7 - base(7)) < length_[7]
            && unsigned(i8 - base(8)) < length_[8];
    }

    _bz_bool isInRange(int i0, int i1, int i2, int i3, int i4,
        int i5, int i6, int i7, int i8, int i9) const
    {
        return unsigned(i0 - base(0)) < length_[0]
            && unsigned(i1 - base(1)) < length_[1]
            && unsigned(i2 - base(2)) < length_[2]
            && unsigned(i3 - base(3)) < length_[3]
            && unsigned(i4 - base(4)) < length_[4]
            && unsigned(i5 - base(5)) < length_[5]
            && unsigned(i6 - base(6)) < length_[6]
            && unsigned(i7 - base(7)) < length_[7]
            && unsigned(i8 - base(8)) < length_[8]
            && unsigned(i9 - base(9)) < length_[9];
    }

    _bz_bool isInRange(int i0, int i1, int i2, int i3, int i4,
        int i5, int i6, int i7, int i8, int i9, int i10) const
    {
        return unsigned(i0 - base(0)) < length_[0]
            && unsigned(i1 - base(1)) < length_[1]
            && unsigned(i2 - base(2)) < length_[2]
            && unsigned(i3 - base(3)) < length_[3]
            && unsigned(i4 - base(4)) < length_[4]
            && unsigned(i5 - base(5)) < length_[5]
            && unsigned(i6 - base(6)) < length_[6]
            && unsigned(i7 - base(7)) < length_[7]
            && unsigned(i8 - base(8)) < length_[8]
            && unsigned(i9 - base(9)) < length_[9]
            && unsigned(i10 - base(10)) < length_[10];
    }

    _bz_bool isInRange(const T_index& index) const
    {
        for (int i=0; i < N_rank; ++i)
            if (unsigned(index[i] - base(i)) >= length_[i])
                return _bz_false;

        return _bz_true;
    }

    _bz_bool assertInRange(const T_index& index) const
    {
        BZPRECHECK(isInRange(index), "Array index out of range: " << index
            << endl << "Lower bounds: " << storage_.base() << endl
            << "Upper bounds: " << (storage_.base() + length_ - 1) << endl);
        return _bz_true;
    }

    _bz_bool assertInRange(int i0) const
    {
        BZPRECHECK(isInRange(i0), "Array index out of range: " << i0
            << endl << "Lower bounds: " << storage_.base() << endl
            << "Upper bounds: " << (storage_.base() + length_ - 1) << endl);
        return _bz_true;
    }

    _bz_bool assertInRange(int i0, int i1) const
    {
        BZPRECHECK(isInRange(i0,i1), "Array index out of range: (" 
            << i0 << ", " << i1 << ")"
            << endl << "Lower bounds: " << storage_.base() << endl
            << "Upper bounds: " << (storage_.base() + length_ - 1) << endl);
        return _bz_true;
    }

    _bz_bool assertInRange(int i0, int i1, int i2) const
    {
        BZPRECHECK(isInRange(i0,i1,i2), "Array index out of range: ("
            << i0 << ", " << i1 << ", " << i2 << ")"
            << endl << "Lower bounds: " << storage_.base() << endl
            << "Upper bounds: " << (storage_.base() + length_ - 1) << endl);
        return _bz_true;
    }

    _bz_bool assertInRange(int i0, int i1, int i2, int i3) const
    {
        BZPRECHECK(isInRange(i0,i1,i2,i3), "Array index out of range: ("
            << i0 << ", " << i1 << ", " << i2 << ", " << i3 << ")"
            << endl << "Lower bounds: " << storage_.base() << endl
            << "Upper bounds: " << (storage_.base() + length_ - 1) << endl);
        return _bz_true;
    }

    _bz_bool assertInRange(int i0, int i1, int i2, int i3, int i4) const
    {
        BZPRECHECK(isInRange(i0,i1,i2,i3,i4), "Array index out of range: ("
            << i0 << ", " << i1 << ", " << i2 << ", " << i3 
            << ", " << i4 << ")"
            << endl << "Lower bounds: " << storage_.base() << endl
            << "Upper bounds: " << (storage_.base() + length_ - 1) << endl);
        return _bz_true;
    }

    _bz_bool assertInRange(int i0, int i1, int i2, int i3, int i4, 
        int i5) const
    {
        BZPRECHECK(isInRange(i0,i1,i2,i3,i4,i5), "Array index out of range: ("
            << i0 << ", " << i1 << ", " << i2 << ", " << i3
            << ", " << i4 << ", " << i5 << ")"
            << endl << "Lower bounds: " << storage_.base() << endl
            << "Upper bounds: " << (storage_.base() + length_ - 1) << endl);
        return _bz_true;
    }

    _bz_bool assertInRange(int i0, int i1, int i2, int i3, int i4,
        int i5, int i6) const
    {
        BZPRECHECK(isInRange(i0,i1,i2,i3,i4,i5,i6), 
            "Array index out of range: ("
            << i0 << ", " << i1 << ", " << i2 << ", " << i3
            << ", " << i4 << ", " << i5 << ", " << i6 << ")"
            << endl << "Lower bounds: " << storage_.base() << endl
            << "Upper bounds: " << (storage_.base() + length_ - 1) << endl);
        return _bz_true;
    }

    _bz_bool assertInRange(int i0, int i1, int i2, int i3, int i4,
        int i5, int i6, int i7) const
    {
        BZPRECHECK(isInRange(i0,i1,i2,i3,i4,i5,i6,i7),
            "Array index out of range: ("
            << i0 << ", " << i1 << ", " << i2 << ", " << i3
            << ", " << i4 << ", " << i5 << ", " << i6 << ", " << i7 << ")"
            << endl << "Lower bounds: " << storage_.base() << endl
            << "Upper bounds: " << (storage_.base() + length_ - 1) << endl);
        return _bz_true;
    }

    _bz_bool assertInRange(int i0, int i1, int i2, int i3, int i4,
        int i5, int i6, int i7, int i8) const
    {
        BZPRECHECK(isInRange(i0,i1,i2,i3,i4,i5,i6,i7,i8),
            "Array index out of range: ("
            << i0 << ", " << i1 << ", " << i2 << ", " << i3
            << ", " << i4 << ", " << i5 << ", " << i6 << ", " << i7 
            << ", " << i8 << ")"
            << endl << "Lower bounds: " << storage_.base() << endl
            << "Upper bounds: " << (storage_.base() + length_ - 1) << endl);
        return _bz_true;
    }

    _bz_bool assertInRange(int i0, int i1, int i2, int i3, int i4,
        int i5, int i6, int i7, int i8, int i9) const
    {
        BZPRECHECK(isInRange(i0,i1,i2,i3,i4,i5,i6,i7,i8,i9),
            "Array index out of range: ("
            << i0 << ", " << i1 << ", " << i2 << ", " << i3
            << ", " << i4 << ", " << i5 << ", " << i6 << ", " << i7
            << ", " << i8 << ", " << i9 << ")"
            << endl << "Lower bounds: " << storage_.base() << endl
            << "Upper bounds: " << (storage_.base() + length_ - 1) << endl);
        return _bz_true;
    }

    _bz_bool assertInRange(int i0, int i1, int i2, int i3, int i4,
        int i5, int i6, int i7, int i8, int i9, int i10) const
    {
        BZPRECHECK(isInRange(i0,i1,i2,i3,i4,i5,i6,i7,i8,i9,i10),
            "Array index out of range: ("
            << i0 << ", " << i1 << ", " << i2 << ", " << i3
            << ", " << i4 << ", " << i5 << ", " << i6 << ", " << i7
            << ", " << i8 << ", " << i9 << ", " << i10 << ")"
            << endl << "Lower bounds: " << storage_.base() << endl
            << "Upper bounds: " << (storage_.base() + length_ - 1) << endl);
        return _bz_true;
    }

    //////////////////////////////////////////////
    // Subscripting operators
    //////////////////////////////////////////////

    template<int N_rank2>
    T_numtype operator()(const TinyVector<int,N_rank2>& index) const
    {
        assertInRange(index);
        return data_[dot(index, stride_)];
    }

    template<int N_rank2>
    T_numtype& _bz_restrict operator()(const TinyVector<int,N_rank2>& index) 
    {
        assertInRange(index);
        return data_[dot(index, stride_)];
    }

    T_numtype operator()(TinyVector<int,1> index) const
    {
        assertInRange(index[0]);
        return data_[index[0] * stride_[0]];
    }

    T_numtype& operator()(TinyVector<int,1> index)
    {
        assertInRange(index[0]);
        return data_[index[0] * stride_[0]];
    }

    T_numtype operator()(TinyVector<int,2> index) const
    {
        assertInRange(index[0], index[1]);
        return data_[index[0] * stride_[0] + index[1] * stride_[1]];
    }

    T_numtype& operator()(TinyVector<int,2> index)
    {
        assertInRange(index[0], index[1]);
        return data_[index[0] * stride_[0] + index[1] * stride_[1]];
    }

    T_numtype operator()(TinyVector<int,3> index) const
    {
        assertInRange(index[0], index[1], index[2]);
        return data_[index[0] * stride_[0] + index[1] * stride_[1]
            + index[2] * stride_[2]];
    }

    T_numtype& operator()(TinyVector<int,3> index)
    {
        assertInRange(index[0], index[1], index[2]);
        return data_[index[0] * stride_[0] + index[1] * stride_[1]
            + index[2] * stride_[2]];
    }

    T_numtype operator()(const TinyVector<int,4>& index) const
    {
        assertInRange(index[0], index[1], index[2], index[3]);
        return data_[index[0] * stride_[0] + index[1] * stride_[1]
            + index[2] * stride_[2] + index[3] * stride_[3]];
    }

    T_numtype& operator()(const TinyVector<int,4>& index)
    {
        assertInRange(index[0], index[1], index[2], index[3]);
        return data_[index[0] * stride_[0] + index[1] * stride_[1]
            + index[2] * stride_[2] + index[3] * stride_[3]];
    }

    T_numtype operator()(const TinyVector<int,5>& index) const
    {
        assertInRange(index[0], index[1], index[2], index[3],
            index[4]);
        return data_[index[0] * stride_[0] + index[1] * stride_[1]
            + index[2] * stride_[2] + index[3] * stride_[3]
            + index[4] * stride_[4]];
    }

    T_numtype& operator()(const TinyVector<int,5>& index)
    {
        assertInRange(index[0], index[1], index[2], index[3],
            index[4]);
        return data_[index[0] * stride_[0] + index[1] * stride_[1]
            + index[2] * stride_[2] + index[3] * stride_[3]
            + index[4] * stride_[4]];
    }

    T_numtype operator()(const TinyVector<int,6>& index) const
    {
        assertInRange(index[0], index[1], index[2], index[3],
            index[4], index[5]);
        return data_[index[0] * stride_[0] + index[1] * stride_[1]
            + index[2] * stride_[2] + index[3] * stride_[3]
            + index[4] * stride_[4] + index[5] * stride_[5]];
    }

    T_numtype& operator()(const TinyVector<int,6>& index)
    {
        assertInRange(index[0], index[1], index[2], index[3],
            index[4], index[5]);
        return data_[index[0] * stride_[0] + index[1] * stride_[1]
            + index[2] * stride_[2] + index[3] * stride_[3]
            + index[4] * stride_[4] + index[5] * stride_[5]];
    }

    T_numtype operator()(const TinyVector<int,7>& index) const
    {
        assertInRange(index[0], index[1], index[2], index[3],
            index[4], index[5], index[6]);
        return data_[index[0] * stride_[0] + index[1] * stride_[1]
            + index[2] * stride_[2] + index[3] * stride_[3]
            + index[4] * stride_[4] + index[5] * stride_[5]
            + index[6] * stride_[6]];
    }

    T_numtype& operator()(const TinyVector<int,7>& index)
    {
        assertInRange(index[0], index[1], index[2], index[3],
            index[4], index[5], index[6]);
        return data_[index[0] * stride_[0] + index[1] * stride_[1]
            + index[2] * stride_[2] + index[3] * stride_[3]
            + index[4] * stride_[4] + index[5] * stride_[5]
            + index[6] * stride_[6]];
    }

    T_numtype operator()(const TinyVector<int,8>& index) const
    {
        assertInRange(index[0], index[1], index[2], index[3],
            index[4], index[5], index[6], index[7]);
        return data_[index[0] * stride_[0] + index[1] * stride_[1]
            + index[2] * stride_[2] + index[3] * stride_[3]
            + index[4] * stride_[4] + index[5] * stride_[5]
            + index[6] * stride_[6] + index[7] * stride_[7]];
    }

    T_numtype& operator()(const TinyVector<int,8>& index)
    {
        assertInRange(index[0], index[1], index[2], index[3],
            index[4], index[5], index[6], index[7]);
        return data_[index[0] * stride_[0] + index[1] * stride_[1]
            + index[2] * stride_[2] + index[3] * stride_[3]
            + index[4] * stride_[4] + index[5] * stride_[5]
            + index[6] * stride_[6] + index[7] * stride_[7]];
    }

    T_numtype operator()(const TinyVector<int,9>& index) const
    {
        assertInRange(index[0], index[1], index[2], index[3],
            index[4], index[5], index[6], index[7], index[8]);
        return data_[index[0] * stride_[0] + index[1] * stride_[1]
            + index[2] * stride_[2] + index[3] * stride_[3]
            + index[4] * stride_[4] + index[5] * stride_[5]
            + index[6] * stride_[6] + index[7] * stride_[7]
            + index[8] * stride_[8]];
    }

    T_numtype& operator()(const TinyVector<int,9>& index)
    {
        assertInRange(index[0], index[1], index[2], index[3],
            index[4], index[5], index[6], index[7], index[8]);
        return data_[index[0] * stride_[0] + index[1] * stride_[1]
            + index[2] * stride_[2] + index[3] * stride_[3]
            + index[4] * stride_[4] + index[5] * stride_[5]
            + index[6] * stride_[6] + index[7] * stride_[7]
            + index[8] * stride_[8]];
    }

    T_numtype operator()(const TinyVector<int,10>& index) const
    {
        assertInRange(index[0], index[1], index[2], index[3],
            index[4], index[5], index[6], index[7], index[8], index[9]);
        return data_[index[0] * stride_[0] + index[1] * stride_[1]
            + index[2] * stride_[2] + index[3] * stride_[3]
            + index[4] * stride_[4] + index[5] * stride_[5]
            + index[6] * stride_[6] + index[7] * stride_[7]
            + index[8] * stride_[8] + index[9] * stride_[9]];
    }

    T_numtype& operator()(const TinyVector<int,10>& index)
    {
        assertInRange(index[0], index[1], index[2], index[3],
            index[4], index[5], index[6], index[7], index[8], index[9]);
        return data_[index[0] * stride_[0] + index[1] * stride_[1]
            + index[2] * stride_[2] + index[3] * stride_[3]
            + index[4] * stride_[4] + index[5] * stride_[5]
            + index[6] * stride_[6] + index[7] * stride_[7]
            + index[8] * stride_[8] + index[9] * stride_[9]];
    }

    T_numtype operator()(const TinyVector<int,11>& index) const
    {
        assertInRange(index[0], index[1], index[2], index[3],
            index[4], index[5], index[6], index[7], index[8], index[9],
            index[10]);
        return data_[index[0] * stride_[0] + index[1] * stride_[1]
            + index[2] * stride_[2] + index[3] * stride_[3]
            + index[4] * stride_[4] + index[5] * stride_[5]
            + index[6] * stride_[6] + index[7] * stride_[7]
            + index[8] * stride_[8] + index[9] * stride_[9]
            + index[10] * stride_[10]];
    }

    T_numtype& operator()(const TinyVector<int,11>& index)
    {
        assertInRange(index[0], index[1], index[2], index[3],
            index[4], index[5], index[6], index[7], index[8], index[9],
            index[10]);
        return data_[index[0] * stride_[0] + index[1] * stride_[1]
            + index[2] * stride_[2] + index[3] * stride_[3]
            + index[4] * stride_[4] + index[5] * stride_[5]
            + index[6] * stride_[6] + index[7] * stride_[7]
            + index[8] * stride_[8] + index[9] * stride_[9]
            + index[10] * stride_[10]];
    }

    T_numtype operator()(int i0) const
    { 
        assertInRange(i0);
        return data_[i0 * stride_[0]]; 
    }

    T_numtype& _bz_restrict operator()(int i0) 
    {
        assertInRange(i0);
        return data_[i0 * stride_[0]];
    }

    T_numtype operator()(int i0, int i1) const
    { 
        assertInRange(i0, i1);
        return data_[i0 * stride_[0] + i1 * stride_[1]];
    }

    T_numtype& _bz_restrict operator()(int i0, int i1)
    {
        assertInRange(i0, i1);
        return data_[i0 * stride_[0] + i1 * stride_[1]];
    }

    T_numtype operator()(int i0, int i1, int i2) const
    {
        assertInRange(i0, i1, i2);
        return data_[i0 * stride_[0] + i1 * stride_[1]
            + i2 * stride_[2]];
    }

    T_numtype& _bz_restrict operator()(int i0, int i1, int i2) 
    {
        assertInRange(i0, i1, i2);
        return data_[i0 * stride_[0] + i1 * stride_[1]
            + i2 * stride_[2]];
    }

    T_numtype operator()(int i0, int i1, int i2, int i3) const
    {
        assertInRange(i0, i1, i2, i3);
        return data_[i0 * stride_[0] + i1 * stride_[1]
            + i2 * stride_[2] + i3 * stride_[3]];
    }

    T_numtype& _bz_restrict operator()(int i0, int i1, int i2, int i3)
    {
        assertInRange(i0, i1, i2, i3);
        return data_[i0 * stride_[0] + i1 * stride_[1]
            + i2 * stride_[2] + i3 * stride_[3]];
    }

    T_numtype operator()(int i0, int i1, int i2, int i3,
        int i4) const
    {
        assertInRange(i0, i1, i2, i3, i4);
        return data_[i0 * stride_[0] + i1 * stride_[1]
            + i2 * stride_[2] + i3 * stride_[3] + i4 * stride_[4]];
    }

    T_numtype& _bz_restrict operator()(int i0, int i1, int i2, int i3,
        int i4)
    {
        assertInRange(i0, i1, i2, i3, i4);
        return data_[i0 * stride_[0] + i1 * stride_[1]
            + i2 * stride_[2] + i3 * stride_[3] + i4 * stride_[4]];
    }

    T_numtype operator()(int i0, int i1, int i2, int i3,
        int i4, int i5) const
    {
        assertInRange(i0, i1, i2, i3, i4, i5);
        return data_[i0 * stride_[0] + i1 * stride_[1]
            + i2 * stride_[2] + i3 * stride_[3] + i4 * stride_[4]
            + i5 * stride_[5]];
    }

    T_numtype& _bz_restrict operator()(int i0, int i1, int i2, int i3,
        int i4, int i5)
    {
        assertInRange(i0, i1, i2, i3, i4, i5);
        return data_[i0 * stride_[0] + i1 * stride_[1]
            + i2 * stride_[2] + i3 * stride_[3] + i4 * stride_[4]
            + i5 * stride_[5]];
    }

    T_numtype operator()(int i0, int i1, int i2, int i3,
        int i4, int i5, int i6) const
    {
        assertInRange(i0, i1, i2, i3, i4, i5, i6);
        return data_[i0 * stride_[0] + i1 * stride_[1]
            + i2 * stride_[2] + i3 * stride_[3] + i4 * stride_[4]
            + i5 * stride_[5] + i6 * stride_[6]];
    }

    T_numtype& _bz_restrict operator()(int i0, int i1, int i2, int i3,
        int i4, int i5, int i6)
    {
        assertInRange(i0, i1, i2, i3, i4, i5, i6);
        return data_[i0 * stride_[0] + i1 * stride_[1]
            + i2 * stride_[2] + i3 * stride_[3] + i4 * stride_[4]
            + i5 * stride_[5] + i6 * stride_[6]];
    }

    T_numtype operator()(int i0, int i1, int i2, int i3,
        int i4, int i5, int i6, int i7) const
    {
        assertInRange(i0, i1, i2, i3, i4, i5, i6, i7);
        return data_[i0 * stride_[0] + i1 * stride_[1]
            + i2 * stride_[2] + i3 * stride_[3] + i4 * stride_[4]
            + i5 * stride_[5] + i6 * stride_[6] + i7 * stride_[7]];
    }

    T_numtype& _bz_restrict operator()(int i0, int i1, int i2, int i3,
        int i4, int i5, int i6, int i7)
    {
        assertInRange(i0, i1, i2, i3, i4, i5, i6, i7);
        return data_[i0 * stride_[0] + i1 * stride_[1]
            + i2 * stride_[2] + i3 * stride_[3] + i4 * stride_[4]
            + i5 * stride_[5] + i6 * stride_[6] + i7 * stride_[7]];
    }

    T_numtype operator()(int i0, int i1, int i2, int i3,
        int i4, int i5, int i6, int i7, int i8) const
    {
        assertInRange(i0, i1, i2, i3, i4, i5, i6, i7, i8);
        return data_[i0 * stride_[0] + i1 * stride_[1]
            + i2 * stride_[2] + i3 * stride_[3] + i4 * stride_[4]
            + i5 * stride_[5] + i6 * stride_[6] + i7 * stride_[7]
            + i8 * stride_[8]];
    }

    T_numtype& _bz_restrict operator()(int i0, int i1, int i2, int i3,
        int i4, int i5, int i6, int i7, int i8)
    {
        assertInRange(i0, i1, i2, i3, i4, i5, i6, i7, i8);
        return data_[i0 * stride_[0] + i1 * stride_[1]
            + i2 * stride_[2] + i3 * stride_[3] + i4 * stride_[4]
            + i5 * stride_[5] + i6 * stride_[6] + i7 * stride_[7]
            + i8 * stride_[8]];
    }

    T_numtype operator()(int i0, int i1, int i2, int i3,
        int i4, int i5, int i6, int i7, int i8, int i9) const
    {
        assertInRange(i0, i1, i2, i3, i4, i5, i6, i7, i8, i9);
        return data_[i0 * stride_[0] + i1 * stride_[1]
            + i2 * stride_[2] + i3 * stride_[3] + i4 * stride_[4]
            + i5 * stride_[5] + i6 * stride_[6] + i7 * stride_[7]
            + i8 * stride_[8] + i9 * stride_[9]];
    }

    T_numtype& _bz_restrict operator()(int i0, int i1, int i2, int i3,
        int i4, int i5, int i6, int i7, int i8, int i9)
    {
        assertInRange(i0, i1, i2, i3, i4, i5, i6, i7, i8, i9);
        return data_[i0 * stride_[0] + i1 * stride_[1]
            + i2 * stride_[2] + i3 * stride_[3] + i4 * stride_[4]
            + i5 * stride_[5] + i6 * stride_[6] + i7 * stride_[7]
            + i8 * stride_[8] + i9 * stride_[9]];
    }

    T_numtype operator()(int i0, int i1, int i2, int i3,
        int i4, int i5, int i6, int i7, int i8, int i9, int i10) const
    {
        assertInRange(i0, i1, i2, i3, i4, i5, i6, i7, i8, 
            i9, i10);
        return data_[i0 * stride_[0] + i1 * stride_[1]
            + i2 * stride_[2] + i3 * stride_[3] + i4 * stride_[4]
            + i5 * stride_[5] + i6 * stride_[6] + i7 * stride_[7]
            + i8 * stride_[8] + i9 * stride_[9] + i10 * stride_[10]];
    }

    T_numtype& _bz_restrict operator()(int i0, int i1, int i2, int i3,
        int i4, int i5, int i6, int i7, int i8, int i9, int i10)
    {
        assertInRange(i0, i1, i2, i3, i4, i5, i6, i7, i8, 
            i9, i10);
        return data_[i0 * stride_[0] + i1 * stride_[1]
            + i2 * stride_[2] + i3 * stride_[3] + i4 * stride_[4]
            + i5 * stride_[5] + i6 * stride_[6] + i7 * stride_[7]
            + i8 * stride_[8] + i9 * stride_[9] + i10 * stride_[10]];
    }

    /*
     * Slicing to produce subarrays.  If the number of Range arguments is
     * fewer than N_rank, then missing arguments are treated like Range::all().
     */

    T_array operator()(const RectDomain<N_rank>& subdomain)
    {
        return T_array(*this, subdomain);
    }

    T_array operator()(Range r0)
    {
        return T_array(*this, r0);
    }

    T_array operator()(Range r0, Range r1)
    {
        return T_array(*this, r0, r1);
    }

    T_array operator()(Range r0, Range r1, Range r2)
    {
        return T_array(*this, r0, r1, r2);
    }

    T_array operator()(Range r0, Range r1, Range r2, Range r3)
    {
        return T_array(*this, r0, r1, r2, r3);
    }

    T_array operator()(Range r0, Range r1, Range r2, Range r3, Range r4)
    {
        return T_array(*this, r0, r1, r2, r3, r4);
    }

    T_array operator()(Range r0, Range r1, Range r2, Range r3, Range r4,
        Range r5)
    {
        return T_array(*this, r0, r1, r2, r3, r4, r5);
    }

    T_array operator()(Range r0, Range r1, Range r2, Range r3, Range r4,
        Range r5, Range r6)
    {
        return T_array(*this, r0, r1, r2, r3, r4, r5, r6);
    }

    T_array operator()(Range r0, Range r1, Range r2, Range r3, Range r4,
        Range r5, Range r6, Range r7)
    {
        return T_array(*this, r0, r1, r2, r3, r4, r5, r6, r7);
    }

    T_array operator()(Range r0, Range r1, Range r2, Range r3, Range r4,
        Range r5, Range r6, Range r7, Range r8)
    {
        return T_array(*this, r0, r1, r2, r3, r4, r5, r6, r7, r8);
    }

    T_array operator()(Range r0, Range r1, Range r2, Range r3, Range r4,
        Range r5, Range r6, Range r7, Range r8, Range r9)
    {
        return T_array(*this, r0, r1, r2, r3, r4, r5, r6, r7, r8, r9);
    }

    T_array operator()(Range r0, Range r1, Range r2, Range r3, Range r4,
        Range r5, Range r6, Range r7, Range r8, Range r9, Range r10)
    {
        return T_array(*this, r0, r1, r2, r3, r4, r5, r6, r7, r8, r9, r10);
    }

    // Allow any mixture of Range, int and Vector<int> objects as
    // operands for operator():   A(Range(3,7), 5, Range(2,4))

    /*
     * These versions of operator() allow any combination of int
     * and Range operands to be used.  Each int operand reduces
     * the rank of the resulting array by one.  
     *
     * e.g.  Array<int,4> A(20,20,20,20);
     *       Array<int,2> B = A(Range(5,15), 3, 5, Range(8,9));
     *
     * SliceInfo is a helper class defined in <blitz/arrayslice.h>.
     * It counts the number of Range vs. int arguments and does some
     * other helpful things.
     *
     * Once partial specialization becomes widely implemented, these
     * operators may be expanded to accept Vector<int> arguments
     * and produce ArrayPick<T,N> objects.
     *
     * This operator() is not provided with a single argument because
     * the appropriate cases exist above.
     */

#ifdef BZ_PARTIAL_ORDERING

    template<class T1, class T2>
    _bz_typename SliceInfo<T_numtype,T1,T2>::T_slice
    operator()(T1 r1, T2 r2)
    {
        return SliceInfo<T_numtype,T1,T2>::T_slice(*this, r1, r2,
            nilArraySection(), nilArraySection(), nilArraySection(),
            nilArraySection(), nilArraySection(), nilArraySection(),
            nilArraySection(), nilArraySection(), nilArraySection());
    }

    template<class T1, class T2, class T3>
    _bz_typename SliceInfo<T_numtype,T1,T2,T3>::T_slice 
    operator()(T1 r1, T2 r2, T3 r3)
    {
        return SliceInfo<T_numtype,T1,T2,T3>::T_slice(*this, r1, r2, r3,
            nilArraySection(), nilArraySection(), nilArraySection(),
            nilArraySection(), nilArraySection(), nilArraySection(),
            nilArraySection(), nilArraySection());
    }

    template<class T1, class T2, class T3, class T4>
    _bz_typename SliceInfo<T_numtype,T1,T2,T3,T4>::T_slice
    operator()(T1 r1, T2 r2, T3 r3, T4 r4)
    {
        return SliceInfo<T_numtype,T1,T2,T3,T4>::T_slice(*this, r1, r2, r3,
            r4, nilArraySection(), nilArraySection(),
            nilArraySection(), nilArraySection(), nilArraySection(),
            nilArraySection(), nilArraySection());
    }

    template<class T1, class T2, class T3, class T4, class T5>
    _bz_typename SliceInfo<T_numtype,T1,T2,T3,T4,T5>::T_slice
    operator()(T1 r1, T2 r2, T3 r3, T4 r4, T5 r5)
    {
        return SliceInfo<T_numtype,T1,T2,T3,T4,T5>::T_slice(*this, r1, r2, r3,
            r4, r5, nilArraySection(),
            nilArraySection(), nilArraySection(), nilArraySection(),
            nilArraySection(), nilArraySection());
    }

    template<class T1, class T2, class T3, class T4, class T5, class T6>
    _bz_typename SliceInfo<T_numtype,T1,T2,T3,T4,T5,T6>::T_slice
    operator()(T1 r1, T2 r2, T3 r3, T4 r4, T5 r5, T6 r6)
    {
        return SliceInfo<T_numtype,T1,T2,T3,T4,T5,T6>::T_slice(*this, r1, r2, 
            r3, r4, r5, r6, 
            nilArraySection(), nilArraySection(), nilArraySection(),
            nilArraySection(), nilArraySection());
    }

    template<class T1, class T2, class T3, class T4, class T5, class T6,
        class T7>
    _bz_typename SliceInfo<T_numtype,T1,T2,T3,T4,T5,T6,T7>::T_slice
    operator()(T1 r1, T2 r2, T3 r3, T4 r4, T5 r5, T6 r6, T7 r7)
    {
        return SliceInfo<T_numtype,T1,T2,T3,T4,T5,T6,T7>::T_slice(*this, r1, 
            r2, r3, r4, r5, r6, r7,
            nilArraySection(), nilArraySection(),
            nilArraySection(), nilArraySection());
    }

    template<class T1, class T2, class T3, class T4, class T5, class T6,
        class T7, class T8>
    _bz_typename SliceInfo<T_numtype,T1,T2,T3,T4,T5,T6,T7,T8>::T_slice
    operator()(T1 r1, T2 r2, T3 r3, T4 r4, T5 r5, T6 r6, T7 r7, T8 r8)
    {
        return SliceInfo<T_numtype,T1,T2,T3,T4,T5,T6,T7,T8>::T_slice(*this, 
            r1, r2, r3, r4, r5, r6, r7, r8,
            nilArraySection(), nilArraySection(), nilArraySection());
    }

    template<class T1, class T2, class T3, class T4, class T5, class T6,
        class T7, class T8, class T9>
    _bz_typename SliceInfo<T_numtype,T1,T2,T3,T4,T5,T6,T7,T8,T9>::T_slice
    operator()(T1 r1, T2 r2, T3 r3, T4 r4, T5 r5, T6 r6, T7 r7, T8 r8, 
        T9 r9)
    {
        return SliceInfo<T_numtype,T1,T2,T3,T4,T5,T6,T7,T8,T9>::T_slice(*this, 
            r1, r2, r3, r4, r5, r6, r7, r8, r9, 
            nilArraySection(), nilArraySection());
    }

    template<class T1, class T2, class T3, class T4, class T5, class T6,
        class T7, class T8, class T9, class T10>
    _bz_typename SliceInfo<T_numtype,T1,T2,T3,T4,T5,T6,T7,T8,T9,T10>::T_slice
    operator()(T1 r1, T2 r2, T3 r3, T4 r4, T5 r5, T6 r6, T7 r7, T8 r8,
        T9 r9, T10 r10)
    {
        return SliceInfo<T_numtype,T1,T2,T3,T4,T5,T6,T7,T8,T9,
            T10>::T_slice(*this, r1, r2, r3, r4, r5, r6, r7, r8, 
            r9, r10, nilArraySection());
    }

    template<class T1, class T2, class T3, class T4, class T5, class T6,
        class T7, class T8, class T9, class T10, class T11>
    _bz_typename SliceInfo<T_numtype,T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,
        T11>::T_slice
    operator()(T1 r1, T2 r2, T3 r3, T4 r4, T5 r5, T6 r6, T7 r7, T8 r8,
        T9 r9, T10 r10, T11 r11)
    {
        return SliceInfo<T_numtype,T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,
            T11>::T_slice(*this, r1, r2, r3, r4, r5, r6, r7, r8, r9, r10, r11);
    }

#endif // BZ_PARTIAL_ORDERING

    /*
     * These versions of operator() are provided to support tensor-style
     * array notation, e.g.
     *
     * Array<float, 2> A, B;
     * firstIndex i;
     * secondIndex j;
     * thirdIndex k;
     * Array<float, 3> C = A(i,j) * B(j,k);
     */

    template<int N0>
    _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0> >
    operator()(IndexPlaceholder<N0>)
    { 
        return _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0> >(*this);
    }

    template<int N0, int N1>
    _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0, N1> >
    operator()(IndexPlaceholder<N0>, IndexPlaceholder<N1>)
    {
        return _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0,
            N1> >(*this);
    } 

    template<int N0, int N1, int N2>
    _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0, N1, N2> >
    operator()(IndexPlaceholder<N0>, IndexPlaceholder<N1>,
        IndexPlaceholder<N2>)
    {
        return _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0,
            N1, N2> >(*this);
    }

    template<int N0, int N1, int N2, int N3>
    _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0, N1, N2, N3> >
    operator()(IndexPlaceholder<N0>, IndexPlaceholder<N1>,
        IndexPlaceholder<N2>, IndexPlaceholder<N3>)
    {
        return _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0,
            N1, N2, N3> >(*this);
    }

    template<int N0, int N1, int N2, int N3, int N4>
    _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0, N1, N2, N3, N4> >
    operator()(IndexPlaceholder<N0>, IndexPlaceholder<N1>,
        IndexPlaceholder<N2>, IndexPlaceholder<N3>, IndexPlaceholder<N4>)
    {
        return _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0,
            N1, N2, N3, N4> >(*this);
    }

    template<int N0, int N1, int N2, int N3, int N4, int N5>
    _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0, N1, N2, N3, 
        N4, N5> >
    operator()(IndexPlaceholder<N0>, IndexPlaceholder<N1>,
        IndexPlaceholder<N2>, IndexPlaceholder<N3>, IndexPlaceholder<N4>,
        IndexPlaceholder<N5>)
    {
        return _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0,
            N1, N2, N3, N4, N5> >(*this);
    }

    template<int N0, int N1, int N2, int N3, int N4, int N5, int N6>
    _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0, N1, N2, N3,
        N4, N5, N6> >
    operator()(IndexPlaceholder<N0>, IndexPlaceholder<N1>,
        IndexPlaceholder<N2>, IndexPlaceholder<N3>, IndexPlaceholder<N4>,
        IndexPlaceholder<N5>, IndexPlaceholder<N6>)
    {
        return _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0,
            N1, N2, N3, N4, N5, N6> >(*this);
    }

    template<int N0, int N1, int N2, int N3, int N4, int N5, int N6,
        int N7>
    _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0, N1, N2, N3,
        N4, N5, N6, N7> >
    operator()(IndexPlaceholder<N0>, IndexPlaceholder<N1>,
        IndexPlaceholder<N2>, IndexPlaceholder<N3>, IndexPlaceholder<N4>,
        IndexPlaceholder<N5>, IndexPlaceholder<N6>, IndexPlaceholder<N7>)
    {
        return _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0,
            N1, N2, N3, N4, N5, N6, N7> >(*this);
    }

    template<int N0, int N1, int N2, int N3, int N4, int N5, int N6,
        int N7, int N8>
    _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0, N1, N2, N3,
        N4, N5, N6, N7, N8> >
    operator()(IndexPlaceholder<N0>, IndexPlaceholder<N1>,
        IndexPlaceholder<N2>, IndexPlaceholder<N3>, IndexPlaceholder<N4>,
        IndexPlaceholder<N5>, IndexPlaceholder<N6>, IndexPlaceholder<N7>,
        IndexPlaceholder<N8>)
    {
        return _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0,
            N1, N2, N3, N4, N5, N6, N7, N8> >(*this);
    }

    template<int N0, int N1, int N2, int N3, int N4, int N5, int N6,
        int N7, int N8, int N9>
    _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0, N1, N2, N3,
        N4, N5, N6, N7, N8, N9> >
    operator()(IndexPlaceholder<N0>, IndexPlaceholder<N1>,
        IndexPlaceholder<N2>, IndexPlaceholder<N3>, IndexPlaceholder<N4>,
        IndexPlaceholder<N5>, IndexPlaceholder<N6>, IndexPlaceholder<N7>,
        IndexPlaceholder<N8>, IndexPlaceholder<N9>)
    {
        return _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0,
            N1, N2, N3, N4, N5, N6, N7, N8, N9> >(*this);
    }

    template<int N0, int N1, int N2, int N3, int N4, int N5, int N6,
        int N7, int N8, int N9, int N10>
    _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0, N1, N2, N3,
        N4, N5, N6, N7, N8, N9, N10> >
    operator()(IndexPlaceholder<N0>, IndexPlaceholder<N1>,
        IndexPlaceholder<N2>, IndexPlaceholder<N3>, IndexPlaceholder<N4>,
        IndexPlaceholder<N5>, IndexPlaceholder<N6>, IndexPlaceholder<N7>,
        IndexPlaceholder<N8>, IndexPlaceholder<N9>, IndexPlaceholder<N10>)
    {
        return _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0,
            N1, N2, N3, N4, N5, N6, N7, N8, N9, N10> >(*this);
    }

    //////////////////////////////////////////////
    // Support for multicomponent arrays
    //////////////////////////////////////////////

    /*
     * See <blitz/array/multi.h> for an explanation of the traits class
     * multicomponent_traits.
     */
    Array<_bz_typename multicomponent_traits<T_numtype>::T_element,N_rank>
    operator[](int component)
    {
        typedef _bz_typename multicomponent_traits<T_numtype>::T_element 
            T_compType;

        return extractComponent(T_compType(), 
            component, multicomponent_traits<T_numtype>::numComponents);
    }

    //////////////////////////////////////////////
    // Indirection
    //////////////////////////////////////////////
  
    template<class T_indexContainer>
    IndirectArray<T_array, T_indexContainer>
    operator[](const T_indexContainer& index)
    {
        return IndirectArray<T_array, T_indexContainer>(*this,
            const_cast<T_indexContainer&>(index));
    }
 
    //////////////////////////////////////////////
    // Assignment Operators
    //////////////////////////////////////////////

    // Scalar operand
    // NEEDS_WORK : need a precondition check on
    // isStorageContiguous when operator, is used.
    ListInitializationSwitch<T_array> operator=(T_numtype x)
    {
        return ListInitializationSwitch<T_array>(*this, x);
    }

    T_array& initialize(T_numtype);

    // Was:
    // T_array& operator=(T_numtype);

    T_array& operator+=(T_numtype);
    T_array& operator-=(T_numtype);
    T_array& operator*=(T_numtype);
    T_array& operator/=(T_numtype);
    T_array& operator%=(T_numtype);
    T_array& operator^=(T_numtype);
    T_array& operator&=(T_numtype);
    T_array& operator|=(T_numtype);
    T_array& operator>>=(T_numtype);
    T_array& operator<<=(T_numtype);

    // Array operands
    T_array& operator=(const Array<T_numtype,N_rank>&);

    template<class T_numtype2> 
    T_array& operator=(const Array<T_numtype2,N_rank>&);
    template<class T_numtype2>
    T_array& operator+=(const Array<T_numtype2,N_rank>&);
    template<class T_numtype2>
    T_array& operator-=(const Array<T_numtype2,N_rank>&);
    template<class T_numtype2>
    T_array& operator*=(const Array<T_numtype2,N_rank>&);
    template<class T_numtype2>
    T_array& operator/=(const Array<T_numtype2,N_rank>&);
    template<class T_numtype2>
    T_array& operator%=(const Array<T_numtype2,N_rank>&);
    template<class T_numtype2>
    T_array& operator^=(const Array<T_numtype2,N_rank>&);
    template<class T_numtype2>
    T_array& operator&=(const Array<T_numtype2,N_rank>&);
    template<class T_numtype2>
    T_array& operator|=(const Array<T_numtype2,N_rank>&);
    template<class T_numtype2>
    T_array& operator>>=(const Array<T_numtype2,N_rank>&);
    template<class T_numtype2>
    T_array& operator<<=(const Array<T_numtype2,N_rank>&);

    // Array expression operands
    template<class T_expr>
    inline T_array& operator=(_bz_ArrayExpr<T_expr> expr);
    template<class T_expr>
    inline T_array& operator+=(_bz_ArrayExpr<T_expr> expr);
    template<class T_expr>
    inline T_array& operator-=(_bz_ArrayExpr<T_expr> expr);
    template<class T_expr>
    inline T_array& operator*=(_bz_ArrayExpr<T_expr> expr);
    template<class T_expr>
    inline T_array& operator/=(_bz_ArrayExpr<T_expr> expr);
    template<class T_expr>
    inline T_array& operator%=(_bz_ArrayExpr<T_expr> expr);
    template<class T_expr>
    inline T_array& operator^=(_bz_ArrayExpr<T_expr> expr);
    template<class T_expr>
    inline T_array& operator&=(_bz_ArrayExpr<T_expr> expr);
    template<class T_expr>
    inline T_array& operator|=(_bz_ArrayExpr<T_expr> expr);
    template<class T_expr>
    inline T_array& operator>>=(_bz_ArrayExpr<T_expr> expr);
    template<class T_expr>
    inline T_array& operator<<=(_bz_ArrayExpr<T_expr> expr);

    // NEEDS_WORK -- Index placeholder operand

    // NEEDS_WORK -- Random operand


public:
    // Undocumented implementation routines

    template<class T_expr, class T_update>
    inline T_array& evaluate(_bz_ArrayExpr<T_expr> expr,
        T_update);

#ifdef BZ_HAVE_STD
    template<class T_expr, class T_update>
    inline T_array& evaluateWithFastTraversal(
        const TraversalOrder<N_rank - 1>& order, _bz_ArrayExpr<T_expr> expr,
        T_update);
#endif

#ifdef BZ_ARRAY_2D_STENCIL_TILING
    template<class T_expr, class T_update>
    inline T_array& evaluateWithTiled2DTraversal(_bz_ArrayExpr<T_expr> expr,
        T_update);
#endif

    template<class T_expr, class T_update>
    inline T_array& evaluateWithIndexTraversal1(_bz_ArrayExpr<T_expr> expr,
        T_update);

    template<class T_expr, class T_update>
    inline T_array& evaluateWithIndexTraversalN(_bz_ArrayExpr<T_expr> expr,
        T_update);

    template<class T_expr, class T_update>
    inline T_array& evaluateWithStackTraversal1(_bz_ArrayExpr<T_expr> expr,
        T_update);

    template<class T_expr, class T_update>
    inline T_array& evaluateWithStackTraversalN(_bz_ArrayExpr<T_expr> expr,
        T_update);

    T_numtype* _bz_restrict getInitializationIterator()
    { return dataFirst(); }

    _bz_bool canCollapse(int outerRank, int innerRank) const
    { 
#ifdef BZ_DEBUG_TRAVERSE
        BZ_DEBUG_MESSAGE("stride(" << innerRank << ")=" << stride(innerRank)
          << ", extent()=" << extent(innerRank) << ", stride(outerRank)="
          << stride(outerRank));
#endif
        return (stride(innerRank) * extent(innerRank) == stride(outerRank)); 
    }

protected:
    //////////////////////////////////////////////
    // Implementation routines
    //////////////////////////////////////////////

    _bz_inline2 void computeStrides();
    _bz_inline2 void setupStorage(int rank);
    void constructSubarray(Array<T_numtype, N_rank>& array, 
        const RectDomain<N_rank>&);
    void constructSubarray(Array<T_numtype, N_rank>& array, Range r0);
    void constructSubarray(Array<T_numtype, N_rank>& array, Range r0, Range r1);
    void constructSubarray(Array<T_numtype, N_rank>& array, Range r0,
        Range r1, Range r2);
    void constructSubarray(Array<T_numtype, N_rank>& array, Range r0,
        Range r1, Range r2, Range r3);
    void constructSubarray(Array<T_numtype, N_rank>& array, Range r0,
        Range r1, Range r2, Range r3, Range r4);
    void constructSubarray(Array<T_numtype, N_rank>& array, Range r0,
        Range r1, Range r2, Range r3, Range r4, Range r5);
    void constructSubarray(Array<T_numtype, N_rank>& array, Range r0,
        Range r1, Range r2, Range r3, Range r4, Range r5, Range r6);
    void constructSubarray(Array<T_numtype, N_rank>& array, Range r0,
        Range r1, Range r2, Range r3, Range r4, Range r5, Range r6,
        Range r7);
    void constructSubarray(Array<T_numtype, N_rank>& array, Range r0,
        Range r1, Range r2, Range r3, Range r4, Range r5, Range r6,
        Range r7, Range r8);
    void constructSubarray(Array<T_numtype, N_rank>& array, Range r0,
        Range r1, Range r2, Range r3, Range r4, Range r5, Range r6,
        Range r7, Range r8, Range r9);
    void constructSubarray(Array<T_numtype, N_rank>& array, Range r0,
        Range r1, Range r2, Range r3, Range r4, Range r5, Range r6,
        Range r7, Range r8, Range r9, Range r10);

    void calculateZeroOffset();

    template<int N_rank2, class R0, class R1, class R2, class R3, class R4, 
        class R5, class R6, class R7, class R8, class R9, class R10>
    void constructSlice(Array<T_numtype, N_rank2>& array, R0 r0, R1 r1, R2 r2, 
        R3 r3, R4 r4, R5 r5, R6 r6, R7 r7, R8 r8, R9 r9, R10 r10);

    template<int N_rank2>
    void slice(int& setRank, Range r, Array<T_numtype,N_rank2>& array,
        TinyVector<int,N_rank2>& rankMap, int sourceRank);

    template<int N_rank2>
    void slice(int& setRank, int i, Array<T_numtype,N_rank2>& array,
        TinyVector<int,N_rank2>& rankMap, int sourceRank);

    template<int N_rank2>
    void slice(int& setRank, nilArraySection, Array<T_numtype,N_rank2>& array,
        TinyVector<int,N_rank2>& rankMap, int sourceRank)
    { }

    void doTranspose(int destRank, int sourceRank, T_array& array);

protected:
    //////////////////////////////////////////////
    // Data members
    //////////////////////////////////////////////

    // NB: adding new data members may require changes to ctors, reference()

    /*
     * For a description of the storage_ members, see the comments for class 
     * GeneralArrayStorage<N_rank> above.
     *
     * length_[] contains the extent of each rank.  E.g. a 10x20x30 array
     *           would have length_ = { 10, 20, 30}.
     * stride_[] contains the stride to move to the next element along each
     *           rank.
     * zeroOffset_ is the distance from the first element in the array 
     *           to the point (0,0,...,0).  If base_ is zero and all ranks are 
     *           stored ascending, then zeroOffset_ is zero.  This value
     *           is needed because to speed up indexing, the data_ member
     *           (inherited from MemoryBlockReference) always refers to
     *           (0,0,...,0).
     * storageContiguous_ is true if all the elements in the array are
     *           stored contiguously in memory.  If this is true, then
     *           the array can be converted to a one-dimensional array.
     *           If the array is sliced, then this becomes false.
     */
    GeneralArrayStorage<N_rank> storage_;
    TinyVector<int, N_rank> length_;
    TinyVector<int, N_rank> stride_;
    int zeroOffset_;
    _bz_bool storageContiguous_;
};

/*
 * Rank numbers start with zero, which may be confusing to users coming
 * from Fortran.  To make code more readable, the following constants
 * may help.  Example: instead of
 * 
 * int firstRankExtent = A.extent(0);
 *
 * One can write:
 *
 * int firstRankExtent = A.extent(firstRank);
 */

const int firstRank    = 0;
const int secondRank   = 1;
const int thirdRank    = 2;
const int fourthRank   = 3;
const int fifthRank    = 4;
const int sixthRank    = 5;
const int seventhRank  = 6;
const int eighthRank   = 7;
const int ninthRank    = 8;
const int tenthRank    = 9;
const int eleventhRank = 10;

const int firstDim    = 0;
const int secondDim   = 1;
const int thirdDim    = 2;
const int fourthDim   = 3;
const int fifthDim    = 4;
const int sixthDim    = 5;
const int seventhDim  = 6;
const int eighthDim   = 7;
const int ninthDim    = 8;
const int tenthDim    = 9;
const int eleventhDim = 10;

/*
 * Global Functions
 */

template<class T_numtype>
ostream& operator<<(ostream&, const Array<T_numtype,1>&);

template<class T_numtype>
ostream& operator<<(ostream&, const Array<T_numtype,2>&);

// NEEDS_WORK -- output formatting for N > 2?

BZ_NAMESPACE_END

/*
 * Include implementations of the member functions and some additional
 * global functions.
 */

#include <blitz/array/iter.h>       // Array iterators
#include <blitz/array/expr.h>       // Array expression objects
#include <blitz/array/methods.cc>   // Member functions
#include <blitz/array/eval.cc>      // Array expression evaluation
#include <blitz/array/ops.cc>       // Assignment operators
#include <blitz/array/io.cc>        // Output formatting
#include <blitz/array/bops.cc>      // Expression templates, two operands
#include <blitz/array/uops.cc>      // Expression templates, math functions
#include <blitz/array/misc.cc>      // Expression templates, miscellaneous
#include <blitz/array/reduce.h>     // Array reduction expression templates
#include <blitz/array/interlace.cc> // Allocation of interlaced arrays
#include <blitz/array/resize.cc>    // Array resize, resizeAndPreserve
#include <blitz/array/slicing.cc>   // Slicing and subarrays
#include <blitz/array/cycle.cc>     // Cycling arrays
#include <blitz/array/complex.cc>   // Special support for complex arrays
#include <blitz/array/zip.h>        // Zipping multicomponent types
#include <blitz/array/where.h>      // where(X,Y,Z)
#include <blitz/array/stencil.h>    // Stencil objects
#include <blitz/array/indirect.h>   // Indirection

#endif // BZ_ARRAY_H