array-impl.h

A C++ class library for scientific computing

// -*- C++ -*-
/***************************************************************************
 * blitz/array-impl.h    Definition of the Array<P_numtype, N_rank> class
 *
 * $Id: array-impl.h,v 1.25 2005/10/13 23:46:43 julianc Exp $
 *
 * Copyright (C) 1997-2001 Todd Veldhuizen <[e-mail]>
 *
 * 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:          [e-mail]
 * Bugs:                 [e-mail]
 *
 * For more information, please see the Blitz++ Home Page:
 *    http://oonumerics.org/blitz/
 *
 ***************************************************************************/

/*
 * 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

#include <blitz/blitz.h>
#include <blitz/memblock.h>
#include <blitz/range.h>
#include <blitz/tinyvec.h>

#ifdef BZ_ARRAY_SPACE_FILLING_TRAVERSAL
#include <blitz/traversal.h>
#endif

#include <blitz/indexexpr.h>
#include <blitz/prettyprint.h>

#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
#include <blitz/array/storage.h>   // GeneralArrayStorage


BZ_NAMESPACE(blitz)

/*
 * Forward declarations
 */

template<typename T_numtype, int N_rank>
class ArrayIterator;

template<typename T_numtype, int N_rank>
class ConstArrayIterator;

template<typename T_numtype, int N_rank>
class FastArrayIterator;

template<typename P_expr>
class _bz_ArrayExpr;

template<typename T_array, typename T_index>
class IndirectArray;

template <typename P_numtype,int N_rank>
void swap(Array<P_numtype,N_rank>&,Array<P_numtype,N_rank>&);

template <typename P_numtype, int N_rank>
void find(Array<TinyVector<int,N_rank>,1>&,const Array<P_numtype,N_rank>&);

/*
 * 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<typename P_numtype, int N_rank>
class Array : public MemoryBlockReference<P_numtype> 
#ifdef BZ_NEW_EXPRESSION_TEMPLATES
    , public ETBase<Array<P_numtype,N_rank> >
#endif
{

private:
    typedef MemoryBlockReference<P_numtype> T_base;
    using T_base::data_;
    using T_base::changeToNullBlock;
    using T_base::numReferences;

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 a a fast iterator for the array, used for expression
     *            templates
     * iterator   is a STL-style iterator
     * const_iterator is an STL-style const iterator
     */

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

    typedef ArrayIterator<T_numtype,N_rank> iterator;
    typedef ConstArrayIterator<T_numtype,N_rank> const_iterator;

    static const int _bz_rank = N_rank;

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

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

    template<typename T_expr>
    explicit 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;
    }

    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 (this is provided for backwards compatibility).
     */
    Array(T_numtype* restrict dataFirst, TinyVector<int, N_rank> shape,
        GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
      : MemoryBlockReference<T_numtype>(product(shape), dataFirst, 
          neverDeleteData),
        storage_(storage)
    {
        BZPRECONDITION(dataFirst != 0);

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

    /*
     * 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* 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, 
          neverDeleteData),