marray.hxx

Go to the documentation of this file.

#pragma once
#ifndef MARRAY_HXX
#define MARRAY_HXX

#include <cassert>
#include <cstddef>
#include <stdexcept> // runtime_error
#include <limits>
#include <string>
#include <sstream>
#include <cstring> // memcpy
#include <iterator> // reverse_iterator, distance
#include <vector>
#include <set>
#include <iostream> // cout
#include <memory> // allocator
#include <numeric> // accumulate

namespace marray {

enum StringStyle {TableStyle, MatrixStyle}; 
enum CoordinateOrder {FirstMajorOrder, LastMajorOrder}; 
struct InitializationSkipping { }; 

static const bool Const = true; 
static const bool Mutable = false; 
static const CoordinateOrder defaultOrder = LastMajorOrder; 
static const InitializationSkipping SkipInitialization = InitializationSkipping(); 

template<class E, class T>
    class ViewExpression;
// \cond suppress doxygen
template<class E, class T, class UnaryFunctor>
    class UnaryViewExpression;
template<class E1, class T1, class E2, class T2, class BinaryFunctor>
    class BinaryViewExpression;
template<class E, class T, class S, class BinaryFunctor>
    class BinaryViewExpressionScalarFirst;
template<class E, class T, class S, class BinaryFunctor>
    class BinaryViewExpressionScalarSecond;
// \endcond suppress doxygen
template<class T, bool isConst = false, class A = std::allocator<size_t> >
    class View;
#ifdef HAVE_CPP0X_TEMPLATE_ALIASES
    template<class T, class A> using ConstView = View<T, true, A>;
#endif
template<class T, bool isConst, class A = std::allocator<size_t> >
    class Iterator;
template<class T, class A = std::allocator<size_t> > class Vector;
template<class T, class A = std::allocator<size_t> > class Matrix;
template<class T, class A = std::allocator<size_t> > class Marray;

// assertion testing
#ifdef NDEBUG
    const bool MARRAY_NO_DEBUG = true; 
    const bool MARRAY_NO_ARG_TEST = true; 
#else
    const bool MARRAY_NO_DEBUG = false; 
    const bool MARRAY_NO_ARG_TEST = false; 
#endif

// \cond suppress doxygen
namespace marray_detail {
    // meta-programming
    template <bool PREDICATE, class TRUECASE, class FALSECASE>
        struct IfBool;
    template <class TRUECASE, class FALSECASE>
        struct IfBool<true, TRUECASE, FALSECASE>
        { typedef TRUECASE type; };
    template <class TRUECASE, class FALSECASE>
        struct IfBool<false, TRUECASE, FALSECASE>
        { typedef FALSECASE type; };

    template <class T1, class T2>
        struct IsEqual
        { static const bool type = false; };
    template <class T>
        struct IsEqual<T, T>
        { static const bool type = true; };

    template<class T> struct TypeTraits
        { static const unsigned char position = 255; };
    template<> struct TypeTraits<char>
        { static const unsigned char position = 0; };
    template<> struct TypeTraits<unsigned char>
        { static const unsigned char position = 1; };
    template<> struct TypeTraits<short>
        { static const unsigned char position = 2; };
    template<> struct TypeTraits<unsigned short>
        { static const unsigned char position = 3; };
    template<> struct TypeTraits<int>
        { static const unsigned char position = 4; };
    template<> struct TypeTraits<unsigned int>
        { static const unsigned char position = 5; };
    template<> struct TypeTraits<long>
        { static const unsigned char position = 6; };
    template<> struct TypeTraits<unsigned long>
        { static const unsigned char position = 7; };
    template<> struct TypeTraits<float>
        { static const unsigned char position = 8; };
    template<> struct TypeTraits<double>
        { static const unsigned char position = 9; };
    template<> struct TypeTraits<long double>
        { static const unsigned char position = 10; };
    template<class A, class B> struct PromoteType
        { typedef typename IfBool<TypeTraits<A>::position >= TypeTraits<B>::position, A, B>::type type; };

    // assertion testing
    template<class A> inline void Assert(A assertion) {
#ifdef NDEBUG
#else
        if(!assertion) throw std::runtime_error("Assertion failed.");
#endif
    }

    // geometry of views
    template<class A = std::allocator<size_t> > class Geometry;
    template<class ShapeIterator, class StridesIterator>
        inline void stridesFromShape(ShapeIterator, ShapeIterator,
            StridesIterator, const CoordinateOrder& = defaultOrder);

    // operations on entries of views
    template<class Functor, class T, class A>
        inline void operate(View<T, false, A>&, Functor);
    template<class Functor, class T, class A>
        inline void operate(View<T, false, A>&, const T&, Functor);
    template<class Functor, class T1, class T2, bool isConst, class A>
        inline void operate(View<T1, false, A>&, const View<T2, isConst, A>&, Functor);
    template<class Functor, class T1, class A, class E, class T2>
        inline void operate(View<T1, false, A>& v, const ViewExpression<E, T2>& expression, Functor f);

    // helper classes
    template<unsigned short N, class Functor, class T, class A>
        struct OperateHelperUnary;
    template<unsigned short N, class Functor, class T1, class T2, class A>
        struct OperateHelperBinaryScalar;
    template<unsigned short N, class Functor, class T1, class T2,
             bool isConst, class A1, class A2>
        struct OperateHelperBinary;
    template<bool isConstTo, class TFrom, class TTo, class AFrom, class ATo>
        struct AssignmentOperatorHelper;
    template<bool isIntegral>
        struct AccessOperatorHelper;

    // unary in-place functors
    template<class T>
        struct Negative { void operator()(T& x) { x = -x; } };
    template<class T>
        struct PrefixIncrement { void operator()(T& x) { ++x; } };
    template<class T>
        struct PostfixIncrement { void operator()(T& x) { x++; } };
    template<class T>
        struct PrefixDecrement { void operator()(T& x) { --x; } };
    template<class T>
        struct PostfixDecrement { void operator()(T& x) { x--; } };

    // binary in-place functors
    template<class T1, class T2>
        struct Assign { void operator()(T1& x, const T2& y) { x = y; } };
    template<class T1, class T2>
        struct PlusEqual { void operator()(T1& x, const T2& y) { x += y; } };
    template<class T1, class T2>
        struct MinusEqual { void operator()(T1& x, const T2& y) { x -= y; } };
    template<class T1, class T2>
        struct TimesEqual { void operator()(T1& x, const T2& y) { x *= y; } };
    template<class T1, class T2>
        struct DividedByEqual { void operator()(T1& x, const T2& y) { x /= y; } };

    // unary functors
    template<class T>
        struct Negate { T operator()(const T& x) const { return -x; } };

    // binary functors
    template<class T1, class T2, class U>
        struct Plus { U operator()(const T1& x, const T2& y) const { return x + y; } };
    template<class T1, class T2, class U>
        struct Minus { U operator()(const T1& x, const T2& y) const { return x - y; } };
    template<class T1, class T2, class U>
        struct Times { U operator()(const T1& x, const T2& y) const { return x * y; } };
    template<class T1, class T2, class U>
        struct DividedBy { U operator()(const T1& x, const T2& y) const { return x / y; } };
}
// \endcond suppress doxygen

template<class T, bool isConst, class A>
class View
: public ViewExpression<View<T, isConst, A>, T>
{
public:
    typedef T value_type;
    typedef T ValueType;
    typedef typename marray_detail::IfBool<isConst, const T*, T*>::type pointer;
    typedef const T* const_pointer;
    typedef typename marray_detail::IfBool<isConst, const T&, T&>::type reference;
    typedef const T& const_reference;
    typedef Iterator<T, isConst, A> iterator;
    typedef Iterator<T, true, A> const_iterator;
    typedef std::reverse_iterator<iterator> reverse_iterator;
    typedef std::reverse_iterator<const_iterator> const_reverse_iterator;
    typedef ViewExpression<View<T, isConst, A>, T> base;
    typedef typename A::template rebind<value_type>::other allocator_type;

    // construction
    View(const allocator_type& = allocator_type());
    View(pointer, const allocator_type& = allocator_type());
    View(const View<T, false, A>&);
    template<class ShapeIterator>
        View(ShapeIterator, ShapeIterator, pointer,
            const CoordinateOrder& = defaultOrder,
            const CoordinateOrder& = defaultOrder,
            const allocator_type& = allocator_type());
    template<class ShapeIterator, class StrideIterator>
        View(ShapeIterator, ShapeIterator, StrideIterator,
            pointer, const CoordinateOrder&,
            const allocator_type& = allocator_type());
    #ifdef HAVE_CPP0X_INITIALIZER_LISTS
        View(std::initializer_list<size_t>, pointer,
            const CoordinateOrder& = defaultOrder,
            const CoordinateOrder& = defaultOrder,
            const allocator_type& = allocator_type());
        View(std::initializer_list<size_t>, std::initializer_list<size_t>,
            pointer, const CoordinateOrder&,
            const allocator_type& = allocator_type());
    #endif

    // assignment
    View<T, isConst, A>& operator=(const T&);
    View<T, isConst, A>& operator=(const View<T, true, A>&); // over-write default
    View<T, isConst, A>& operator=(const View<T, false, A>&); // over-write default
    template<class TLocal, bool isConstLocal, class ALocal>
        View<T, isConst, A>& operator=(const View<TLocal, isConstLocal, ALocal>&);
    template<class E, class Te>
        View<T, isConst, A>&
        operator=(const ViewExpression<E, Te>&);

    void assign(const allocator_type& = allocator_type());
    template<class ShapeIterator>
        void assign(ShapeIterator, ShapeIterator, pointer,
            const CoordinateOrder& = defaultOrder,
            const CoordinateOrder& = defaultOrder,
            const allocator_type& = allocator_type());
    template<class ShapeIterator, class StrideIterator>
        void assign(ShapeIterator, ShapeIterator, StrideIterator,
            pointer, const CoordinateOrder&,
            const allocator_type& = allocator_type());
    #ifdef HAVE_CPP0X_INITIALIZER_LISTS
        void assign(std::initializer_list<size_t>, pointer,
            const CoordinateOrder& = defaultOrder,
            const CoordinateOrder& = defaultOrder,
            const allocator_type& = allocator_type());
        void assign(std::initializer_list<size_t>,
            std::initializer_list<size_t>, pointer,
            const CoordinateOrder&,
            const allocator_type& = allocator_type());
    #endif

    // query
    const size_t dimension() const;
    const size_t size() const;
    const size_t shape(const size_t) const;
    const size_t* shapeBegin() const;
    const size_t* shapeEnd() const;
    const size_t strides(const size_t) const;
    const size_t* stridesBegin() const;
    const size_t* stridesEnd() const;
    const CoordinateOrder& coordinateOrder() const;
    const bool isSimple() const;
    template<class TLocal, bool isConstLocal, class ALocal>
        bool overlaps(const View<TLocal, isConstLocal, ALocal>&) const;

    // element access
    template<class U> reference operator()(U);
    template<class U> reference operator()(U) const;
    #ifndef HAVE_CPP0X_VARIADIC_TEMPLATES
        reference operator()(const size_t, const size_t);
        reference operator()(const size_t, const size_t) const;
        reference operator()(const size_t, const size_t, const size_t);
        reference operator()(const size_t, const size_t, const size_t) const;
        reference operator()(const size_t, const size_t, const size_t,
            const size_t);
        reference operator()(const size_t, const size_t, const size_t,
            const size_t) const;
        reference operator()(const size_t, const size_t, const size_t,
             const size_t, const size_t);
        reference operator()(const size_t, const size_t, const size_t,
            const size_t, const size_t) const;
        reference operator()(const size_t, const size_t, const size_t,
            const size_t, const size_t, const size_t, const size_t,
            const size_t, const size_t, const size_t);
        reference operator()(const size_t, const size_t, const size_t,
            const size_t, const size_t, const size_t, const size_t,
            const size_t, const size_t, const size_t) const;
    #else
        reference operator()(const size_t);
        reference operator()(const size_t) const;
        template<typename... Args>
            reference operator()(const size_t, const Args...);
        template<typename... Args>
            reference operator()(const size_t, const Args...) const;
        private:
            size_t elementAccessHelper(const size_t, const size_t);
            size_t elementAccessHelper(const size_t, const size_t) const;
            template<typename... Args>
                size_t elementAccessHelper(const size_t, const size_t,
                    const Args...);
            template<typename... Args>
                size_t elementAccessHelper(const size_t, const size_t,
                    const Args...) const;
        public:
    #endif

    // sub-views
    template<class BaseIterator, class ShapeIterator>
        void view(BaseIterator, ShapeIterator, View<T, isConst, A>&) const;
    template<class BaseIterator, class ShapeIterator>
        void view(BaseIterator, ShapeIterator, const CoordinateOrder&,
            View<T, isConst, A>&) const;
    template<class BaseIterator, class ShapeIterator>
        View<T, isConst, A> view(BaseIterator, ShapeIterator) const;
    template<class BaseIterator, class ShapeIterator>
        View<T, isConst, A> view(BaseIterator, ShapeIterator,
            const CoordinateOrder&) const;
    template<class BaseIterator, class ShapeIterator>
        void constView(BaseIterator, ShapeIterator, View<T, true, A>&) const;
    template<class BaseIterator, class ShapeIterator>
        void constView(BaseIterator, ShapeIterator, const CoordinateOrder&,
            View<T, true, A>&) const;
    template<class BaseIterator, class ShapeIterator>
        View<T, true, A> constView(BaseIterator, ShapeIterator) const;
    template<class BaseIterator, class ShapeIterator>
        View<T, true, A> constView(BaseIterator, ShapeIterator,
            const CoordinateOrder&) const;
    #ifdef HAVE_CPP0X_INITIALIZER_LISTS
        void view(std::initializer_list<size_t>,
            std::initializer_list<size_t>, View<T, isConst, A>&) const;
        void view(std::initializer_list<size_t>,
            std::initializer_list<size_t>, const CoordinateOrder&,
            View<T, isConst, A>&) const;
        void constView(std::initializer_list<size_t>,
            std::initializer_list<size_t>, View<T, true, A>&) const;
        void constView(std::initializer_list<size_t>,
            std::initializer_list<size_t>, const CoordinateOrder&,
            View<T, true, A>&) const;
    #endif

    // iterator access
    iterator begin();
    iterator end();
    const_iterator begin() const;
    const_iterator end() const;
    reverse_iterator rbegin();
    reverse_iterator rend();
    const_reverse_iterator rbegin() const;
    const_reverse_iterator rend() const;

    // coordinate transformation
    template<class ShapeIterator>
        void reshape(ShapeIterator, ShapeIterator);
    template<class CoordinateIterator>
        void permute(CoordinateIterator);
    void transpose(const size_t, const size_t);
    void transpose();
    void shift(const int);
    void squeeze();

    template<class ShapeIterator>
        View<T, isConst, A> reshapedView(ShapeIterator, ShapeIterator) const;
    template<class CoordinateIterator>
        View<T, isConst, A> permutedView(CoordinateIterator) const;
    View<T, isConst, A> transposedView(const size_t, const size_t) const;
    View<T, isConst, A> transposedView() const;
    View<T, isConst, A> shiftedView(const int) const;
    View<T, isConst, A> boundView(const size_t, const size_t = 0) const;
    View<T, isConst, A> squeezedView() const;

    #ifdef HAVE_CPP0X_INITIALIZER_LISTS
        void reshape(std::initializer_list<size_t>);
        void permute(std::initializer_list<size_t>);

        View<T, isConst, A> reshapedView(std::initializer_list<size_t>) const;
        View<T, isConst, A> permutedView(std::initializer_list<size_t>) const;
    #endif

    // conversion between coordinates, index and offset
    template<class CoordinateIterator>
        void coordinatesToIndex(CoordinateIterator, size_t&) const;
    template<class CoordinateIterator>
        void coordinatesToOffset(CoordinateIterator, size_t&) const;
    template<class CoordinateIterator>
        void indexToCoordinates(size_t, CoordinateIterator) const;
    void indexToOffset(size_t, size_t&) const;
    #ifdef HAVE_CPP0X_INITIALIZER_LISTS
        void coordinatesToIndex(std::initializer_list<size_t>,
            size_t&) const;
        void coordinatesToOffset(std::initializer_list<size_t>,
            size_t&) const;
    #endif

    // output as string
    std::string asString(const StringStyle& = MatrixStyle) const;

private:
    typedef typename marray_detail::Geometry<A> geometry_type;

    View(pointer, const geometry_type&);

    void updateSimplicity();
    void testInvariant() const;

    // unsafe direct memory access
    const T& operator[](const size_t) const;
    T& operator[](const size_t);

    // data and memory address
    pointer data_;
    geometry_type geometry_;

template<class TLocal, bool isConstLocal, class ALocal>
    friend class View;
template<class TLocal, class ALocal>
    friend class Marray;
template<class TLocal, class ALocal>
    friend class Vector;
template<class TLocal, class ALocal>
    friend class Matrix;
// \cond suppress doxygen
template<bool isConstTo, class TFrom, class TTo, class AFrom, class ATo>
    friend struct marray_detail::AssignmentOperatorHelper;
friend struct marray_detail::AccessOperatorHelper<true>;
friend struct marray_detail::AccessOperatorHelper<false>;

template<class E, class U>
    friend class ViewExpression;
template<class E, class U, class UnaryFunctor>
    friend class UnaryViewExpression;
template<class E1, class T1, class E2, class T2, class BinaryFunctor>
    friend class BinaryViewExpression;
template<class E, class U, class S, class BinaryFunctor>
    friend class BinaryViewExpressionScalarFirst;
template<class E, class U, class S, class BinaryFunctor>
    friend class BinaryViewExpressionScalarSecond;

template<class Functor, class T1, class Alocal, class E, class T2>
    friend void marray_detail::operate(View<T1, false, Alocal>& v, const ViewExpression<E, T2>& expression, Functor f);
// \endcond end suppress doxygen
};

template<class T, bool isConst, class A>
class Iterator
{
public:
    // STL random access iterator typedefs
    typedef typename std::random_access_iterator_tag iterator_category;
    typedef T value_type;
    typedef std::ptrdiff_t difference_type;
    typedef typename marray_detail::IfBool<isConst, const T*, T*>::type pointer;
    typedef typename marray_detail::IfBool<isConst, const T&, T&>::type reference;

    // non-standard typedefs
    typedef typename marray_detail::IfBool<isConst, const View<T, true, A>*,
        View<T, false, A>*>::type view_pointer;
    typedef typename marray_detail::IfBool<isConst, const View<T, true, A>&,
        View<T, false, A>&>::type view_reference;

    // construction
    Iterator();
    Iterator(const View<T, false, A>&, const size_t = 0);
    Iterator(View<T, false, A>&, const size_t = 0);
    Iterator(const View<T, true, A>&, const size_t = 0);
    Iterator(const Iterator<T, false, A>&);
        // conversion from mutable to const

    // STL random access iterator operations
    reference operator*() const;
    pointer operator->() const;
    reference operator[](const size_t) const;
    Iterator<T, isConst, A>& operator+=(const difference_type&);
    Iterator<T, isConst, A>& operator-=(const difference_type&);
    Iterator<T, isConst, A>& operator++(); // prefix

    Iterator<T, isConst, A>& operator--(); // prefix
    Iterator<T, isConst, A> operator++(int); // postfix
    Iterator<T, isConst, A> operator--(int); // postfix
    Iterator<T, isConst, A> operator+(const difference_type&) const;
    Iterator<T, isConst, A> operator-(const difference_type&) const;
    template<bool isConstLocal>
        difference_type operator-(const Iterator<T, isConstLocal, A>&) const;
    template<bool isConstLocal>
        bool operator==(const Iterator<T, isConstLocal, A>&) const;
    template<bool isConstLocal>
        bool operator!=(const Iterator<T, isConstLocal, A>&) const;
    template<bool isConstLocal>
        bool operator<(const Iterator<T, isConstLocal, A>&) const;
    template<bool isConstLocal>
        bool operator>(const Iterator<T, isConstLocal, A>&) const;
    template<bool isConstLocal>
        bool operator<=(const Iterator<T, isConstLocal, A>&) const;
    template<bool isConstLocal>
        bool operator>=(const Iterator<T, isConstLocal, A>&) const;

    // operations beyond the STL standard
    bool hasMore() const;
    size_t index() const;
    template<class CoordinateIterator>
        void coordinate(CoordinateIterator) const;

private:
    void testInvariant() const;

    // attributes
    view_pointer view_;
    pointer pointer_;
    size_t index_;
    std::vector<size_t> coordinates_;

friend class Marray<T, A>;
friend class Iterator<T, !isConst, A>; // for comparison operators
};

template<class T, class A>
class Marray
: public View<T, false, A>
{
public:
    typedef View<T, false, A> base;
    typedef typename base::value_type value_type;
    typedef typename base::pointer pointer;
    typedef typename base::const_pointer const_pointer;
    typedef typename base::reference reference;
    typedef typename base::const_reference const_reference;
    typedef typename base::iterator iterator;
    typedef typename base::reverse_iterator reverse_iterator;
    typedef typename base::const_iterator const_iterator;
    typedef typename base::const_reverse_iterator const_reverse_iterator;
    typedef typename A::template rebind<value_type>::other allocator_type;

    // constructors and destructor
    Marray(const allocator_type& = allocator_type());
    Marray(const T&, const CoordinateOrder& = defaultOrder,
        const allocator_type& = allocator_type());
    template<class ShapeIterator>
        Marray(ShapeIterator, ShapeIterator, const T& = T(),
            const CoordinateOrder& = defaultOrder,
            const allocator_type& = allocator_type());
    template<class ShapeIterator>
        Marray(const InitializationSkipping&, ShapeIterator, ShapeIterator,
            const CoordinateOrder& = defaultOrder,
            const allocator_type& = allocator_type());
    #ifdef HAVE_CPP0X_INITIALIZER_LISTS
        Marray(std::initializer_list<size_t>, const T& = T(),
            const CoordinateOrder& = defaultOrder,
            const allocator_type& = allocator_type());
    #endif
    Marray(const Marray<T, A>&);
    template<class E, class Te>
        Marray(const ViewExpression<E, Te>&,
            const allocator_type& = allocator_type());
    template<class TLocal, bool isConstLocal, class ALocal>
        Marray(const View<TLocal, isConstLocal, ALocal>&);
    ~Marray();

    // assignment
    Marray<T, A>& operator=(const T&);
    Marray<T, A>& operator=(const Marray<T, A>&); // over-write default
    template<class TLocal, bool isConstLocal, class ALocal>
        Marray<T, A>& operator=(const View<TLocal, isConstLocal, ALocal>&);
    template<class E, class Te>
        Marray<T, A>& operator=(const ViewExpression<E, Te>&);
    void assign(const allocator_type& = allocator_type());

    // resize
    template<class ShapeIterator>
        void resize(ShapeIterator, ShapeIterator, const T& = T());
    template<class ShapeIterator>
        void resize(const InitializationSkipping&, ShapeIterator, ShapeIterator);
    #ifdef HAVE_CPP0X_INITIALIZER_LISTS
        void resize(std::initializer_list<size_t>, const T& = T());
        void resize(const InitializationSkipping&, std::initializer_list<size_t>);
    #endif

private:
    typedef typename base::geometry_type geometry_type;

    // hiding inherited functions
    template<class CoordinateIterator>
      void permute(CoordinateIterator) {}
    void transpose(const size_t, const size_t) {}
    void transpose() {}
    void shift(const int) {}
    void squeeze() {}

    void testInvariant() const;
    template<bool SKIP_INITIALIZATION, class ShapeIterator>
        void resizeHelper(ShapeIterator, ShapeIterator, const T& = T());

    allocator_type dataAllocator_;

friend class Vector<T, A>;
friend class Matrix<T, A>;
};

template<class T, class A>
class Vector
: public Marray<T, A>
{
public:
    typedef Marray<T, A> base;
    typedef typename base::value_type value_type;
    typedef typename base::pointer pointer;
    typedef typename base::const_pointer const_pointer;
    typedef typename base::reference reference;
    typedef typename base::const_reference const_reference;
    typedef typename base::iterator iterator;
    typedef typename base::reverse_iterator reverse_iterator;
    typedef typename base::const_iterator const_iterator;
    typedef typename base::const_reverse_iterator const_reverse_iterator;
    typedef typename base::allocator_type allocator_type;

    // constructors and destructor
    Vector(const allocator_type& = allocator_type());
    template<class TLocal, bool isConstLocal, class ALocal>
        Vector(const View<TLocal, isConstLocal, ALocal>&);
    Vector(const size_t, const T& = T(),
        const allocator_type& = allocator_type());
    Vector(const InitializationSkipping&, const size_t,
        const allocator_type& = allocator_type());
    template<class E, class Te>
        Vector(const ViewExpression<E, Te>&,
            const allocator_type& = allocator_type());
    #ifdef HAVE_CPP0X_INITIALIZER_LISTS
        Vector(std::initializer_list<T> list,
            const allocator_type& = allocator_type());
    #endif

    // assignment operator
    Vector<T, A>& operator=(const T&);
    Vector<T, A>& operator=(const Vector<T, A>&); // over-write default
    template<class TLocal, bool isConstLocal, class ALocal>
        Vector<T, A>& operator=(const View<TLocal, isConstLocal, ALocal>&);
    template<class E, class Te>
        Vector<T, A>& operator=(const ViewExpression<E, Te>&);

    // element access
    T& operator[](const size_t);
    const T& operator[](const size_t) const;

    // reshape, resize
    void reshape(const size_t);
    void resize(const size_t, const T& = T());
    void resize(const InitializationSkipping&, const size_t);

private:
    void testInvariant() const;
};

template<class T, class A>
class Matrix
: public Marray<T, A>
{
public:
    typedef Marray<T, A> base;

    typedef typename base::value_type value_type;
    typedef typename base::pointer pointer;

    typedef typename base::const_pointer const_pointer;
    typedef typename base::reference reference;
    typedef typename base::const_reference const_reference;
    typedef typename base::iterator iterator;
    typedef typename base::reverse_iterator reverse_iterator;
    typedef typename base::const_iterator const_iterator;
    typedef typename base::const_reverse_iterator const_reverse_iterator;
    typedef typename base::allocator_type allocator_type;

    // constructors and destructor
    Matrix(const allocator_type& = allocator_type());
    template<class TLocal, bool isConstLocal, class ALocal>
        Matrix(const View<TLocal, isConstLocal, ALocal>&);
    Matrix(const size_t, const size_t, const T& = T(),
        const CoordinateOrder& = defaultOrder,
        const allocator_type& = allocator_type());
    Matrix(const InitializationSkipping&,
        const size_t, const size_t,
        const CoordinateOrder& = defaultOrder,
        const allocator_type& = allocator_type());
    template<class E, class Te>
        Matrix(const ViewExpression<E, Te>&,
            const allocator_type& = allocator_type());

    // assignment operator
    Matrix<T, A>& operator=(const T&);
    Matrix<T, A>& operator=(const Matrix<T, A>&); // over-write default
        // overwrite standard operator=
    template<class TLocal, bool isConstLocal, class ALocal>
        Matrix<T, A>& operator=(const View<TLocal, isConstLocal, ALocal>&);
    template<class E, class Te>
        Matrix<T, A>& operator=(const ViewExpression<E, Te>&);

    // resize and reshape
    void reshape(const size_t, const size_t);
    void resize(const size_t, const size_t, const T& = T());
    void resize(const InitializationSkipping&, const size_t, const size_t);

private:
    void testInvariant() const;
};

// implementation of View

#ifdef HAVE_CPP0X_INITIALIZER_LISTS
template<class T, bool isConst, class A>
inline void
View<T, isConst, A>::coordinatesToIndex
(
    std::initializer_list<size_t> coordinate,
    size_t& out
) const
{
    coordinatesToIndex(coordinate.begin(), out);
}
#endif

template<class T, bool isConst, class A>
template<class CoordinateIterator>
inline void
View<T, isConst, A>::coordinatesToIndex
(
    CoordinateIterator it,
    size_t& out
) const
{
    testInvariant();
    out = 0;
    for(size_t j=0; j<this->dimension(); ++j, ++it) {
        marray_detail::Assert(MARRAY_NO_ARG_TEST || static_cast<size_t>(*it) < shape(j));
        out += static_cast<size_t>(*it) * geometry_.shapeStrides(j);
    }
}

#ifdef HAVE_CPP0X_INITIALIZER_LISTS
template<class T, bool isConst, class A>
inline void
View<T, isConst, A>::coordinatesToOffset
(
    std::initializer_list<size_t> coordinate,
    size_t& out
) const
{
    coordinatesToOffset(coordinate.begin(), out);
}
#endif

template<class T, bool isConst, class A>
template<class CoordinateIterator>
inline void
View<T, isConst, A>::coordinatesToOffset
(
    CoordinateIterator it,
    size_t& out
) const
{
    testInvariant();
    out = 0;
    for(size_t j=0; j<this->dimension(); ++j, ++it) {
        marray_detail::Assert(MARRAY_NO_ARG_TEST || static_cast<size_t>(*it) < shape(j));
        out += static_cast<size_t>(*it) * strides(j);
    }
}

template<class T, bool isConst, class A>
template<class CoordinateIterator>
inline void
View<T, isConst, A>::indexToCoordinates
(
    size_t index, // copy to work on
    CoordinateIterator outit
) const
{
    testInvariant();
    marray_detail::Assert(MARRAY_NO_DEBUG || this->dimension() > 0);
    marray_detail::Assert(MARRAY_NO_ARG_TEST || index < this->size());
    if(coordinateOrder() == FirstMajorOrder) {
        for(size_t j=0; j<this->dimension(); ++j, ++outit) {
            *outit = size_t(index / geometry_.shapeStrides(j));
            index = index % geometry_.shapeStrides(j);
        }
    }
    else { // last major order
        size_t j = this->dimension()-1;
        outit += j;
        for(;;) {
            *outit = size_t(index / geometry_.shapeStrides(j));
            index = index % geometry_.shapeStrides(j);
            if(j == 0) {
                break;
            }
            else {
                --outit;
                --j;
            }
        }
    }
}

template<class T, bool isConst, class A>
inline void
View<T, isConst, A>::indexToOffset
(
    size_t index,
    size_t& out
) const
{
    testInvariant();
    marray_detail::Assert(MARRAY_NO_ARG_TEST || index < this->size());
    if(isSimple()) {
        out = index;
    }
    else {
        out = 0;
        if(coordinateOrder() == FirstMajorOrder) {
            for(size_t j=0; j<this->dimension(); ++j) {
                out += geometry_.strides(j) * (index / geometry_.shapeStrides(j));
                index = index % geometry_.shapeStrides(j);
            }
        }
        else { // last major order
            if(this->dimension() == 0) {
                marray_detail::Assert(MARRAY_NO_ARG_TEST || index == 0);
                return;
            }
            else {
                size_t j = this->dimension()-1;
                for(;;) {
                    out += geometry_.strides(j) * (index / geometry_.shapeStrides(j));
                    index = index % geometry_.shapeStrides(j);
                    if(j == 0) {
                        break;
                    }
                    else {
                        --j;
                    }
                }
            }
        }
    }
}

template<class T, bool isConst, class A>
inline
View<T, isConst, A>::View
(
    const allocator_type& allocator
)
: data_(0),
  geometry_(geometry_type(allocator))
{
    testInvariant();
}

// private constructor
template<class T, bool isConst, class A>
inline
View<T, isConst, A>::View
(
    pointer data,
    const geometry_type& geometry
)
: data_(data),
  geometry_(geometry)
{
    testInvariant();
}

template<class T, bool isConst, class A>
inline
View<T, isConst, A>::View
(
    pointer data,
    const allocator_type& allocator
)
: data_(data),
  geometry_(geometry_type(0, defaultOrder, 1, true, allocator))
{
    testInvariant();
}

template<class T, bool isConst, class A>
inline
View<T, isConst, A>::View
(
    const View<T, false, A>& in
)
: data_(in.data_),
  geometry_(in.geometry_)
{
    testInvariant();
}


template<class T, bool isConst, class A>
template<class ShapeIterator>
inline
View<T, isConst, A>::View
(
    ShapeIterator begin,
    ShapeIterator end,
    pointer data,
    const CoordinateOrder& externalCoordinateOrder,
    const CoordinateOrder& internalCoordinateOrder,
    const allocator_type& allocator
)
:   data_(data),
    geometry_(begin, end, externalCoordinateOrder,
        internalCoordinateOrder, allocator)
{
    testInvariant();
}

template<class T, bool isConst, class A>
template<class ShapeIterator, class StrideIterator>
inline
View<T, isConst, A>::View
(
    ShapeIterator begin,
    ShapeIterator end,
    StrideIterator it,
    pointer data,
    const CoordinateOrder& internalCoordinateOrder,
    const allocator_type& allocator
)
: data_(data),
  geometry_(begin, end, it, internalCoordinateOrder, allocator)
{
    testInvariant();
}

#ifdef HAVE_CPP0X_INITIALIZER_LISTS
template<class T, bool isConst, class A>
inline
View<T, isConst, A>::View
(
    std::initializer_list<size_t> shape,
    pointer data,
    const CoordinateOrder& externalCoordinateOrder,
    const CoordinateOrder& internalCoordinateOrder,
    const allocator_type& allocator
)
:   data_(data),
    geometry_(shape.begin(), shape.end(), externalCoordinateOrder,
              internalCoordinateOrder, allocator)
{
    testInvariant();
}

template<class T, bool isConst, class A>
inline
View<T, isConst, A>::View
(
    std::initializer_list<size_t> shape,
    std::initializer_list<size_t> strides,
    pointer data,
    const CoordinateOrder& internalCoordinateOrder,
    const allocator_type& allocator
)
:   data_(data),
    geometry_(shape.begin(), shape.end(), strides.begin(),
              internalCoordinateOrder, allocator)
{
    testInvariant();
}
#endif

template<class T, bool isConst, class A>
inline void
View<T, isConst, A>::assign
(
    const allocator_type& allocator
)
{
    geometry_ = geometry_type(allocator);
    data_ = 0;
    testInvariant();
}

template<class T, bool isConst, class A>
template<class ShapeIterator>
inline void
View<T, isConst, A>::assign
(
    ShapeIterator begin,
    ShapeIterator end,
    pointer data,
    const CoordinateOrder& externalCoordinateOrder,
    const CoordinateOrder& internalCoordinateOrder,
    const allocator_type& allocator
)
{
    // the invariant is not tested as a pre-condition of this
    // function to allow for unsafe manipulations prior to its
    // call
    geometry_ = typename marray_detail::Geometry<A>(begin, end,
        externalCoordinateOrder, internalCoordinateOrder, allocator);
    data_ = data;
    testInvariant();
}

template<class T, bool isConst, class A>
template<class ShapeIterator, class StrideIterator>
inline void
View<T, isConst, A>::assign
(
    ShapeIterator begin,
    ShapeIterator end,
    StrideIterator it,
    pointer data,
    const CoordinateOrder& internalCoordinateOrder,
    const allocator_type& allocator
)
{
    // the invariant is not tested as a pre-condition of this
    // function to allow for unsafe manipulations prior to its
    // call
    geometry_ = typename marray_detail::Geometry<A>(begin, end,
        it, internalCoordinateOrder, allocator);
    data_ = data;
    testInvariant();
}

#ifdef HAVE_CPP0X_INITIALIZER_LISTS
template<class T, bool isConst, class A>
inline void
View<T, isConst, A>::assign
(
    std::initializer_list<size_t> shape,
    pointer data,
    const CoordinateOrder& externalCoordinateOrder,
    const CoordinateOrder& internalCoordinateOrder,
    const allocator_type& allocator
)
{
    assign(shape.begin(), shape.end(), data, externalCoordinateOrder,
        internalCoordinateOrder, allocator);
}

template<class T, bool isConst, class A>
inline void
View<T, isConst, A>::assign
(
    std::initializer_list<size_t> shape,
    std::initializer_list<size_t> strides,
    pointer data,
    const CoordinateOrder& internalCoordinateOrder,
    const allocator_type& allocator
)
{
    assign(shape.begin(), shape.end(), strides.begin(), data,
        internalCoordinateOrder, allocator);
}
#endif

template<class T, bool isConst, class A>
template<class U>
inline typename View<T, isConst, A>::reference
View<T, isConst, A>::operator()
(
    U u
)
{
    return marray_detail::AccessOperatorHelper<std::numeric_limits<U>::is_integer>::execute(*this, u);
}

template<class T, bool isConst, class A>
template<class U>
inline typename View<T, isConst, A>::reference
View<T, isConst, A>::operator()
(
    U u
) const
{
    return marray_detail::AccessOperatorHelper<std::numeric_limits<U>::is_integer>::execute(*this, u);
}

#ifndef HAVE_CPP0X_VARIADIC_TEMPLATES

template<class T, bool isConst, class A>
inline typename View<T, isConst, A>::reference
View<T, isConst, A>::operator()
(
    const size_t c0,
    const size_t c1
)
{
    testInvariant();
    marray_detail::Assert(MARRAY_NO_DEBUG || (data_ != 0 && dimension() == 2));
    marray_detail::Assert(MARRAY_NO_ARG_TEST || (c0 < shape(0) && c1 < shape(1)));
    return data_[c0*strides(0) + c1*strides(1)];
}

template<class T, bool isConst, class A>
inline typename View<T, isConst, A>::reference
View<T, isConst, A>::operator()
(
    const size_t c0,
    const size_t c1
) const
{
    testInvariant();
    marray_detail::Assert(MARRAY_NO_DEBUG || (data_ != 0 && dimension() == 2));
    marray_detail::Assert(MARRAY_NO_ARG_TEST || (c0 < shape(0) && c1 < shape(1)));
    return data_[c0*strides(0) + c1*strides(1)];
}

template<class T, bool isConst, class A>
inline typename View<T, isConst, A>::reference
View<T, isConst, A>::operator()
(
    const size_t c0,
    const size_t c1,
    const size_t c2
)
{
    testInvariant();
    marray_detail::Assert(MARRAY_NO_DEBUG || (data_ != 0 && dimension() == 3));
    marray_detail::Assert(MARRAY_NO_ARG_TEST || (c0 < shape(0) && c1 < shape(1)
        && c2 < shape(2)));
    return data_[c0*strides(0) + c1*strides(1) + c2*strides(2)];
}

template<class T, bool isConst, class A>
inline typename View<T, isConst, A>::reference
View<T, isConst, A>::operator()
(
    const size_t c0,
    const size_t c1,
    const size_t c2
) const
{
    testInvariant();
    marray_detail::Assert(MARRAY_NO_DEBUG || (data_ != 0 && dimension() == 3));
    marray_detail::Assert(MARRAY_NO_ARG_TEST || (c0 < shape(0) && c1 < shape(1)
        && c2 < shape(2)));
    return data_[c0*strides(0) + c1*strides(1) + c2*strides(2)];
}

template<class T, bool isConst, class A>
inline typename View<T, isConst, A>::reference
View<T, isConst, A>::operator()
(
    const size_t c0,
    const size_t c1,
    const size_t c2,
    const size_t c3
)
{
    testInvariant();
    marray_detail::Assert(MARRAY_NO_DEBUG || (data_ != 0 && dimension() == 4));
    marray_detail::Assert(MARRAY_NO_ARG_TEST || (c0 < shape(0) && c1 < shape(1)
        && c2 < shape(2) && c3 < shape(3)));
    return data_[c0*strides(0) + c1*strides(1) + c2*strides(2)
        + c3*strides(3)];
}

template<class T, bool isConst, class A>
inline typename View<T, isConst, A>::reference
View<T, isConst, A>::operator()
(
    const size_t c0,
    const size_t c1,
    const size_t c2,
    const size_t c3
) const
{
    testInvariant();
    marray_detail::Assert(MARRAY_NO_DEBUG || (data_ != 0 && dimension() == 4));
    marray_detail::Assert(MARRAY_NO_ARG_TEST || (c0 < shape(0) && c1 < shape(1)
        && c2 < shape(2) && c3 < shape(3)));
    return data_[c0*strides(0) + c1*strides(1) + c2*strides(2)
        + c3*strides(3)];
}

template<class T, bool isConst, class A>
inline typename View<T, isConst, A>::reference
View<T, isConst, A>::operator()
(
    const size_t c0,
    const size_t c1,
    const size_t c2,
    const size_t c3,
    const size_t c4
)
{
    testInvariant();
    marray_detail::Assert(MARRAY_NO_DEBUG || (data_ != 0 && dimension() == 5));
    marray_detail::Assert(MARRAY_NO_ARG_TEST || (c0 < shape(0) && c1 < shape(1)
        && c2 < shape(2) && c3 < shape(3) && c4 < shape(4)));
    return data_[c0*strides(0) + c1*strides(1) + c2*strides(2)
        + c3*strides(3) + c4*strides(4)];
}

template<class T, bool isConst, class A>
inline typename View<T, isConst, A>::reference
View<T, isConst, A>::operator()
(
    const size_t c0,
    const size_t c1,
    const size_t c2,
    const size_t c3,
    const size_t c4
) const
{
    testInvariant();
    marray_detail::Assert(MARRAY_NO_DEBUG || (data_ != 0 && dimension() == 5));
    marray_detail::Assert(MARRAY_NO_ARG_TEST || (c0 < shape(0) && c1 < shape(1)
        && c2 < shape(2) && c3 < shape(3) && c4 < shape(4)));
    return data_[c0*strides(0) + c1*strides(1) + c2*strides(2)
        + c3*strides(3) + c4*strides(4)];
}


template<class T, bool isConst, class A>
inline typename View<T, isConst, A>::reference
View<T, isConst, A>::operator()
(
    const size_t c0,
    const size_t c1,
    const size_t c2,
    const size_t c3,
    const size_t c4,
    const size_t c5,
    const size_t c6,
    const size_t c7,
    const size_t c8,
    const size_t c9
)
{
    testInvariant();
    marray_detail::Assert(MARRAY_NO_DEBUG || (data_ != 0 && dimension() == 10));
    marray_detail::Assert(MARRAY_NO_ARG_TEST || (c0 < shape(0) && c1 < shape(1)
        && c2 < shape(2) && c3 < shape(3) && c4 < shape(4)
        && c5 < shape(5) && c6 < shape(6) && c7 < shape(7)
        && c8 < shape(8) && c9 < shape(9)));
    return data_[c0*strides(0) + c1*strides(1) + c2*strides(2)
        + c3*strides(3) + c4*strides(4) + c5*strides(5) + c6*strides(6)
        + c7*strides(7) + c8*strides(8) + c9*strides(9)];
}

template<class T, bool isConst, class A>
inline typename View<T, isConst, A>::reference
View<T, isConst, A>::operator()
(
    const size_t c0,
    const size_t c1,
    const size_t c2,
    const size_t c3,
    const size_t c4,
    const size_t c5,
    const size_t c6,
    const size_t c7,
    const size_t c8,
    const size_t c9
) const
{
    testInvariant();
    marray_detail::Assert(MARRAY_NO_DEBUG || (data_ != 0 && dimension() == 10));
    marray_detail::Assert(MARRAY_NO_ARG_TEST || (c0 < shape(0) && c1 < shape(1)
        && c2 < shape(2) && c3 < shape(3) && c4 < shape(4)
        && c5 < shape(5) && c6 < shape(6) && c7 < shape(7)
        && c8 < shape(8) && c9 < shape(9)));
    return data_[c0*strides(0) + c1*strides(1) + c2*strides(2)
        + c3*strides(3) + c4*strides(4) + c5*strides(5) + c6*strides(6)
        + c7*strides(7) + c8*strides(8) + c9*strides(9)];
}

#else

template<class T, bool isConst, class A>
inline size_t
View<T, isConst, A>::elementAccessHelper
(
    const size_t Dim,
    const size_t value
)
{
    marray_detail::Assert(MARRAY_NO_ARG_TEST || (value < shape(Dim-1) ) );
    return strides(Dim-1) * value;
}

template<class T, bool isConst, class A>
inline size_t
View<T, isConst, A>::elementAccessHelper
(
    const size_t Dim,
    const size_t value
) const
{
    marray_detail::Assert(MARRAY_NO_ARG_TEST || (value < shape(Dim-1) ) );
    return strides(Dim-1) * value;
}

template<class T, bool isConst, class A>
template<typename... Args>
inline size_t
View<T, isConst, A>::elementAccessHelper
(
    const size_t Dim,
    const size_t value,
    const Args... args
)
{
    marray_detail::Assert(MARRAY_NO_ARG_TEST || (value < shape(Dim-1-sizeof...(args)) ) );
    return value * strides(Dim-1-sizeof...(args)) + elementAccessHelper(Dim, args...);
}

template<class T, bool isConst, class A>
template<typename... Args>
inline size_t
View<T, isConst, A>::elementAccessHelper
(
    const size_t Dim,
    const size_t value,
    const Args... args
) const
{
    marray_detail::Assert(MARRAY_NO_ARG_TEST || (value < shape(Dim-1-sizeof...(args)) ) );
    return value * strides(Dim-1-sizeof...(args)) + elementAccessHelper(Dim, args...);
}

template<class T, bool isConst, class A>
inline typename View<T, isConst, A>::reference
View<T, isConst, A>::operator()
(
    const size_t value
)
{
    testInvariant();
    if(dimension() == 0) {
        marray_detail::Assert(MARRAY_NO_ARG_TEST || value == 0);
        return data_[0];
    }
    else {
        size_t offset;
        indexToOffset(value, offset);
        return data_[offset];
    }
}

template<class T, bool isConst, class A>
inline typename View<T, isConst, A>::reference
View<T, isConst, A>::operator()
(
    const size_t value
) const
{
    testInvariant();
    if(dimension() == 0) {
        marray_detail::Assert(MARRAY_NO_ARG_TEST || value == 0);
        return data_[0];
    }
    else {
        size_t offset;
        indexToOffset(value, offset);
        return data_[offset];
    }
}

template<class T, bool isConst, class A>
template<typename... Args>
inline typename View<T, isConst, A>::reference
View<T, isConst, A>::operator()
(
    const size_t value,
    const Args... args
)
{
    testInvariant();
    marray_detail::Assert( MARRAY_NO_DEBUG || ( data_ != 0 && sizeof...(args)+1 == dimension() ) );
    return data_[strides(0)*value + elementAccessHelper(sizeof...(args)+1, args...) ];
}

template<class T, bool isConst, class A>
template<typename... Args>
inline typename View<T, isConst, A>::reference
View<T, isConst, A>::operator()
(
    const size_t value,
    const Args... args
) const
{
    testInvariant();
    marray_detail::Assert( MARRAY_NO_DEBUG || ( data_ != 0 && sizeof...(args)+1 == dimension() ) );
    return data_[ strides(0) * static_cast<size_t>(value)
        + static_cast<size_t>(elementAccessHelper(sizeof...(args)+1, args...)) ];
}

#endif // #ifndef HAVE_CPP0X_VARIADIC_TEMPLATES

template<class T, bool isConst, class A>
inline const size_t
View<T, isConst, A>::size() const
{
    return geometry_.size();
}

template<class T, bool isConst, class A>
inline const size_t
View<T, isConst, A>::dimension() const
{
    marray_detail::Assert(MARRAY_NO_DEBUG || this->data_ != 0);
    return geometry_.dimension();
}

template<class T, bool isConst, class A>
inline const size_t
View<T, isConst, A>::shape
(
    const size_t dimension
) const
{
    testInvariant();
    marray_detail::Assert(MARRAY_NO_DEBUG || data_ != 0);
    marray_detail::Assert(MARRAY_NO_ARG_TEST || dimension < this->dimension());
    return geometry_.shape(dimension);
}

template<class T, bool isConst, class A>
inline const size_t*
View<T, isConst, A>::shapeBegin() const
{
    testInvariant();
    marray_detail::Assert(MARRAY_NO_DEBUG || data_ != 0);
    return geometry_.shapeBegin();
}

template<class T, bool isConst, class A>
inline const size_t*
View<T, isConst, A>::shapeEnd() const
{
    testInvariant();
    marray_detail::Assert(MARRAY_NO_DEBUG || data_ != 0);
    return geometry_.shapeEnd();
}

template<class T, bool isConst, class A>
inline const size_t
View<T, isConst, A>::strides
(
    const size_t dimension
) const
{
    testInvariant();
    marray_detail::Assert(MARRAY_NO_DEBUG || data_ != 0);
    marray_detail::Assert(MARRAY_NO_ARG_TEST || dimension < this->dimension());
    return geometry_.strides(dimension);
}

template<class T, bool isConst, class A>
inline const size_t*
View<T, isConst, A>::stridesBegin() const
{
    testInvariant();
    marray_detail::Assert(MARRAY_NO_DEBUG || data_ != 0);
    return geometry_.stridesBegin();
}

template<class T, bool isConst, class A>
inline const size_t*
View<T, isConst, A>::stridesEnd() const
{
    testInvariant();
    marray_detail::Assert(MARRAY_NO_DEBUG || data_ != 0);
    return geometry_.stridesEnd();
}

template<class T, bool isConst, class A>
inline const CoordinateOrder&
View<T, isConst, A>::coordinateOrder() const
{
    testInvariant();
    return geometry_.coordinateOrder();
}

template<class T, bool isConst, class A>
inline const bool
View<T, isConst, A>::isSimple() const
{
    testInvariant();
    return geometry_.isSimple();
}

template<class T, bool isConst, class A>
inline View<T, isConst, A>&
View<T, isConst, A>::operator=
(
    const View<T, true, A>& in
)
{
    testInvariant();
    marray_detail::AssignmentOperatorHelper<isConst, T, T, A, A>::execute(in, *this);
    testInvariant();
    return *this;
}

template<class T, bool isConst, class A>
inline View<T, isConst, A>&
View<T, isConst, A>::operator=
(
    const View<T, false, A>& in
)
{
    testInvariant();
    marray_detail::AssignmentOperatorHelper<isConst, T, T, A, A>::execute(in, *this);
    testInvariant();
    return *this;
}

template<class T, bool isConst, class A>
template<class TLocal, bool isConstLocal, class ALocal>
inline View<T, isConst, A>&
View<T, isConst, A>::operator=
(
    const View<TLocal, isConstLocal, ALocal>& in
)
{
    testInvariant();
    marray_detail::AssignmentOperatorHelper<isConst, TLocal, T, ALocal, A>::execute(in, *this);
    testInvariant();
    return *this;
}

template<class T, bool isConst, class A>
inline View<T, isConst, A>&
View<T, isConst, A>::operator=
(
    const T& value
)
{
    marray_detail::Assert(MARRAY_NO_DEBUG || data_ != 0);
    if(isSimple()) {
        for(size_t j=0; j<geometry_.size(); ++j) {
            data_[j] = value;
        }
    }
    else if(dimension() == 1)
        marray_detail::OperateHelperBinaryScalar<1, marray_detail::Assign<T, T>, T, T, A>::operate(*this, value, marray_detail::Assign<T, T>(), data_);
    else if(dimension() == 2)
        marray_detail::OperateHelperBinaryScalar<2, marray_detail::Assign<T, T>, T, T, A>::operate(*this, value, marray_detail::Assign<T, T>(), data_);
    else if(dimension() == 3)
        marray_detail::OperateHelperBinaryScalar<3, marray_detail::Assign<T, T>, T, T, A>::operate(*this, value, marray_detail::Assign<T, T>(), data_);
    else if(dimension() == 4)
        marray_detail::OperateHelperBinaryScalar<4, marray_detail::Assign<T, T>, T, T, A>::operate(*this, value, marray_detail::Assign<T, T>(), data_);
    else if(dimension() == 5)
        marray_detail::OperateHelperBinaryScalar<5, marray_detail::Assign<T, T>, T, T, A>::operate(*this, value, marray_detail::Assign<T, T>(), data_);
    else if(dimension() == 6)
        marray_detail::OperateHelperBinaryScalar<6, marray_detail::Assign<T, T>, T, T, A>::operate(*this, value, marray_detail::Assign<T, T>(), data_);
    else if(dimension() == 7)
        marray_detail::OperateHelperBinaryScalar<7, marray_detail::Assign<T, T>, T, T, A>::operate(*this, value, marray_detail::Assign<T, T>(), data_);
    else if(dimension() == 8)
        marray_detail::OperateHelperBinaryScalar<8, marray_detail::Assign<T, T>, T, T, A>::operate(*this, value, marray_detail::Assign<T, T>(), data_);
    else if(dimension() == 9)
        marray_detail::OperateHelperBinaryScalar<9, marray_detail::Assign<T, T>, T, T, A>::operate(*this, value, marray_detail::Assign<T, T>(), data_);
    else if(dimension() == 10)
        marray_detail::OperateHelperBinaryScalar<10, marray_detail::Assign<T, T>, T, T, A>::operate(*this, value, marray_detail::Assign<T, T>(), data_);
    else {
        for(iterator it = begin(); it.hasMore(); ++it) {
            *it = value;
        }
    }
    return *this;
}

template<class T, bool isConst, class A>
template<class E, class Te>
inline View<T, isConst, A>&
View<T, isConst, A>::operator=
(
    const ViewExpression<E, Te>& expression
)
{
    marray_detail::operate(*this, expression, marray_detail::Assign<T, Te>());
    return *this;
}

template<class T, bool isConst, class A>
template<class BaseIterator, class ShapeIterator>
inline void
View<T, isConst, A>::view
(
    BaseIterator bit,
    ShapeIterator sit,
    View<T, isConst, A>& out
) const
{
    view(bit, sit, coordinateOrder(), out);
}

template<class T, bool isConst, class A>
template<class BaseIterator, class ShapeIterator>
inline void
View<T, isConst, A>::view
(
    BaseIterator bit,
    ShapeIterator sit,
    const CoordinateOrder& internalCoordinateOrder,
    View<T, isConst, A>& out
) const
{
    testInvariant();
    size_t offset = 0;
    coordinatesToOffset(bit, offset);
    out.assign(sit, sit+dimension(), geometry_.stridesBegin(),
        data_+offset, internalCoordinateOrder);
}

template<class T, bool isConst, class A>
template<class BaseIterator, class ShapeIterator>
inline View<T, isConst, A>
View<T, isConst, A>::view
(
    BaseIterator bit,
    ShapeIterator sit
) const
{
    View<T, isConst, A> v;
    this->view(bit, sit, v);
    return v;
}

template<class T, bool isConst, class A>
template<class BaseIterator, class ShapeIterator>
inline View<T, isConst, A>
View<T, isConst, A>::view
(
    BaseIterator bit,
    ShapeIterator sit,
    const CoordinateOrder& internalCoordinateOrder
) const
{
    View<T, isConst, A> v;
    this->view(bit, sit, internalCoordinateOrder, v);
    return v;
}

template<class T, bool isConst, class A>
template<class BaseIterator, class ShapeIterator>
inline void
View<T, isConst, A>::constView
(
    BaseIterator bit,
    ShapeIterator sit,
    View<T, true, A>& out
) const
{
    constView(bit, sit, coordinateOrder(), out);
}

template<class T, bool isConst, class A>
template<class BaseIterator, class ShapeIterator>
inline void
View<T, isConst, A>::constView
(
    BaseIterator bit,
    ShapeIterator sit,
    const CoordinateOrder& internalCoordinateOrder,
    View<T, true, A>& out
) const
{
    testInvariant();
    size_t offset = 0;
    coordinatesToOffset(bit, offset);
    out.assign(sit, sit+dimension(),
        geometry_.stridesBegin(),
        static_cast<const T*>(data_) + offset,
        internalCoordinateOrder);
}

template<class T, bool isConst, class A>
template<class BaseIterator, class ShapeIterator>
inline View<T, true, A>
View<T, isConst, A>::constView
(
    BaseIterator bit,
    ShapeIterator sit
) const
{
    View<T, true, A> v;
    this->constView(bit, sit, v);
    return v;
}

template<class T, bool isConst, class A>
template<class BaseIterator, class ShapeIterator>
inline View<T, true, A>
View<T, isConst, A>::constView
(
    BaseIterator bit,
    ShapeIterator sit,
    const CoordinateOrder& internalCoordinateOrder
) const
{
    View<T, true, A> v;
    this->constView(bit, sit, internalCoordinateOrder, v);
    return v;
}

#ifdef HAVE_CPP0X_INITIALIZER_LISTS
template<class T, bool isConst, class A>
inline void
View<T, isConst, A>::view
(
    std::initializer_list<size_t> b,
    std::initializer_list<size_t> s,
    const CoordinateOrder& internalCoordinateOrder,
    View<T, isConst, A>& out
) const
{
    view(b.begin(), s.begin(), internalCoordinateOrder, out);
}

template<class T, bool isConst, class A>
inline void
View<T, isConst, A>::view
(
    std::initializer_list<size_t> b,
    std::initializer_list<size_t> s,
    View<T, isConst, A>& out
) const
{
    view(b.begin(), s.begin(), coordinateOrder(), out);
}

template<class T, bool isConst, class A>
inline void
View<T, isConst, A>::constView
(
    std::initializer_list<size_t> b,
    std::initializer_list<size_t> s,
    const CoordinateOrder& internalCoordinateOrder,
    View<T, true, A>& out
) const
{
    constView(b.begin(), s.begin(), internalCoordinateOrder, out);
}

template<class T, bool isConst, class A>
inline void
View<T, isConst, A>::constView
(
    std::initializer_list<size_t> b,
    std::initializer_list<size_t> s,
    View<T, true, A>& out
) const
{
    constView(b.begin(), s.begin(), coordinateOrder(), out);
}
#endif

template<class T, bool isConst, class A>
template<class ShapeIterator>
inline void
View<T, isConst, A>::reshape
(
    ShapeIterator begin,
    ShapeIterator end
)
{
    testInvariant();
    marray_detail::Assert(MARRAY_NO_DEBUG || isSimple());
    if(!MARRAY_NO_ARG_TEST) {
        size_t size = std::accumulate(begin, end, 1, std::multiplies<size_t>());
        marray_detail::Assert(size == this->size());
    }
    assign(begin, end, data_, coordinateOrder(), coordinateOrder());
    testInvariant();
}

template<class T, bool isConst, class A>
template<class ShapeIterator>
inline View<T, isConst, A>
View<T, isConst, A>::reshapedView
(
    ShapeIterator begin,
    ShapeIterator end
) const
{
    View<T, isConst, A> out = *this;
    out.reshape(begin, end);
    return out;
}

#ifdef HAVE_CPP0X_INITIALIZER_LISTS
template<class T, bool isConst, class A>
inline void
View<T, isConst, A>::reshape
(
    std::initializer_list<size_t> shape
)
{
    reshape(shape.begin(), shape.end());
}

template<class T, bool isConst, class A>
inline View<T, isConst, A>
View<T, isConst, A>::reshapedView
(
    std::initializer_list<size_t> shape
) const
{
    return reshapedView(shape.begin(), shape.end());
}
#endif

template<class T, bool isConst, class A>
View<T, isConst, A>
View<T, isConst, A>::boundView
(
    const size_t dimension,
    const size_t value
) const
{
    testInvariant();
    marray_detail::Assert(MARRAY_NO_ARG_TEST || (dimension < this->dimension()
        && value < shape(dimension)));
    if(this->dimension() == 1) {
        View v(&((*this)(value)));
        v.geometry_.coordinateOrder() = coordinateOrder();
        return v;
    }
    else {
        View v;
        v.geometry_.resize(this->dimension()-1);
        v.geometry_.coordinateOrder() = coordinateOrder();
        v.geometry_.size() = size() / shape(dimension);
        for(size_t j=0, k=0; j<this->dimension(); ++j) {
            if(j != dimension) {
                v.geometry_.shape(k) = shape(j);
                v.geometry_.strides(k) = strides(j);
                ++k;
            }
        }
        marray_detail::stridesFromShape(v.geometry_.shapeBegin(), v.geometry_.shapeEnd(),
            v.geometry_.shapeStridesBegin(), v.geometry_.coordinateOrder());
        v.data_ = data_ + strides(dimension) * value;
        v.updateSimplicity();
        v.testInvariant();
        return v;
    }
}

template<class T, bool isConst, class A>
void
View<T, isConst, A>::squeeze()
{
    testInvariant();
    if(dimension() != 0) {
        size_t newDimension = dimension();
        for(size_t j=0; j<dimension(); ++j) {
            if(shape(j) == 1) {
                --newDimension;
            }
        }
        if(newDimension != dimension()) {
            if(newDimension == 0) {
                geometry_.resize(0);
                geometry_.size() = 1;
            }
            else {
                for(size_t j=0, k=0; j<geometry_.dimension(); ++j) {
                    if(geometry_.shape(j) != 1) {
                        geometry_.shape(k) = geometry_.shape(j);
                        geometry_.strides(k) = geometry_.strides(j);
                        ++k;
                    }
                }
                geometry_.resize(newDimension);
                marray_detail::stridesFromShape(geometry_.shapeBegin(), geometry_.shapeEnd(),
                    geometry_.shapeStridesBegin(), geometry_.coordinateOrder());
                updateSimplicity();
            }
        }
    }
    testInvariant();
}

template<class T, bool isConst, class A>
inline View<T, isConst, A>
View<T, isConst, A>::squeezedView() const
{
    View<T, isConst, A> v = *this;
    v.squeeze();
    return v;
}

#ifdef HAVE_CPP0X_INITIALIZER_LISTS
template<class T, bool isConst, class A>
void
View<T, isConst, A>::permute
(
    std::initializer_list<size_t> permutation
)
{
    permute(permutation.begin());
}
#endif

template<class T, bool isConst, class A>
template<class CoordinateIterator>
void
View<T, isConst, A>::permute
(
    CoordinateIterator begin
)
{
    testInvariant();
    if(!MARRAY_NO_ARG_TEST) {
        marray_detail::Assert(dimension() != 0);
        std::set<size_t> s1, s2;
        CoordinateIterator it = begin;
        for(size_t j=0; j<dimension(); ++j) {
            s1.insert(j);
            s2.insert(*it);
            ++it;
        }
        marray_detail::Assert(s1 == s2);
    }
    // update shape, shape strides, strides, and simplicity
    std::vector<size_t> newShape = std::vector<size_t>(dimension());
    std::vector<size_t> newStrides = std::vector<size_t>(dimension());
    for(size_t j=0; j<dimension(); ++j) {
        newShape[j] = geometry_.shape(static_cast<size_t>(*begin));
        newStrides[j] = geometry_.strides(static_cast<size_t>(*begin));
        ++begin;
    }
    for(size_t j=0; j<dimension(); ++j) {
        geometry_.shape(j) = newShape[j];
        geometry_.strides(j) = newStrides[j];
    }
    marray_detail::stridesFromShape(geometry_.shapeBegin(), geometry_.shapeEnd(),
        geometry_.shapeStridesBegin(), geometry_.coordinateOrder());
    updateSimplicity();
    testInvariant();
}

template<class T, bool isConst, class A>
template<class CoordinateIterator>
inline View<T, isConst, A>
View<T, isConst, A>::permutedView
(
    CoordinateIterator begin
) const
{
    View<T, isConst, A> out = *this;
    out.permute(begin);
    return out;
}

template<class T, bool isConst, class A>
void
View<T, isConst, A>::transpose
(
    const size_t c1,
    const size_t c2
)
{
    testInvariant();
    marray_detail::Assert(MARRAY_NO_ARG_TEST ||
        (dimension() != 0 && c1 < dimension() && c2 < dimension()));

    size_t j1 = c1;
    size_t j2 = c2;
    size_t c;
    size_t d;

    // transpose shape
    c = geometry_.shape(j2);
    geometry_.shape(j2) = geometry_.shape(j1);
    geometry_.shape(j1) = c;

    // transpose strides
    d = geometry_.strides(j2);
    geometry_.strides(j2) = geometry_.strides(j1);
    geometry_.strides(j1) = d;

    // update shape strides
    marray_detail::stridesFromShape(geometry_.shapeBegin(), geometry_.shapeEnd(),
        geometry_.shapeStridesBegin(), geometry_.coordinateOrder());

    updateSimplicity();
    testInvariant();
}

template<class T, bool isConst, class A>
void
View<T, isConst, A>::transpose()
{
    testInvariant();
    for(size_t j=0; j<dimension()/2; ++j) {
        size_t k = dimension()-j-1;

        // transpose shape
        size_t tmp = geometry_.shape(j);
        geometry_.shape(j) = geometry_.shape(k);
        geometry_.shape(k) = tmp;

        // transpose strides
        tmp = geometry_.strides(j);
        geometry_.strides(j) = geometry_.strides(k);
        geometry_.strides(k) = tmp;
    }
    marray_detail::stridesFromShape(geometry_.shapeBegin(), geometry_.shapeEnd(),
        geometry_.shapeStridesBegin(), geometry_.coordinateOrder());
    updateSimplicity();
    testInvariant();

}

template<class T, bool isConst, class A>
inline View<T, isConst, A>
View<T, isConst, A>::transposedView
(
    const size_t c1,
    const size_t c2
) const
{
    View<T, isConst, A> out = *this;
    out.transpose(c1, c2);
    return out;
}

template<class T, bool isConst, class A>
inline View<T, isConst, A>
View<T, isConst, A>::transposedView() const
{
    View<T, isConst, A> out = *this;
    out.transpose();
    return out;
}

template<class T, bool isConst, class A>
inline void
View<T, isConst, A>::shift
(
    const int n
)
{
    testInvariant();
    marray_detail::Assert(MARRAY_NO_DEBUG || dimension() != 0);
    if(n <= -static_cast<int>(dimension()) || n >= static_cast<int>(dimension())) {
        shift(n % static_cast<int>(dimension()));
    }
    else {
        if(n > 0) {
            shift(n - static_cast<int>(dimension()));
        }
        else {
            std::vector<size_t> p(dimension());
            for(size_t j=0; j<dimension(); ++j) {
                p[j] = static_cast<size_t>((static_cast<int>(j) - n)) % dimension();
            }
            permute(p.begin());
        }
    }
    testInvariant();
}

template<class T, bool isConst, class A>
inline View<T, isConst, A>
View<T, isConst, A>::shiftedView
(
    const int n
) const
{
    View<T, isConst, A> out = *this;
    out.shift(n);
    return out;
}


template<class T, bool isConst, class A>
inline typename View<T, isConst, A>::iterator
View<T, isConst, A>::begin()
{
    testInvariant();
    return Iterator<T, isConst, A>(*this, 0);
}

template<class T, bool isConst, class A>
inline typename View<T, isConst, A>::iterator
View<T, isConst, A>::end()
{
    testInvariant();
    return Iterator<T, isConst, A>(*this, geometry_.size());
}

template<class T, bool isConst, class A>
inline typename View<T, isConst, A>::const_iterator
View<T, isConst, A>::begin() const
{
    testInvariant();
    return Iterator<T, true>(*this, 0);
}

template<class T, bool isConst, class A>
inline typename View<T, isConst, A>::const_iterator

View<T, isConst, A>::end() const
{
    testInvariant();
    return Iterator<T, true>(*this, geometry_.size());
}

template<class T, bool isConst, class A>
inline typename View<T, isConst, A>::reverse_iterator
View<T, isConst, A>::rbegin()
{
    return reverse_iterator(end());
}

template<class T, bool isConst, class A>
inline typename View<T, isConst, A>::reverse_iterator
View<T, isConst, A>::rend()
{
    return reverse_iterator(begin());
}

template<class T, bool isConst, class A>
inline typename View<T, isConst, A>::const_reverse_iterator
View<T, isConst, A>::rbegin() const
{
    return const_reverse_iterator(end());
}

template<class T, bool isConst, class A>
inline typename View<T, isConst, A>::const_reverse_iterator
View<T, isConst, A>::rend() const
{
    return const_reverse_iterator(begin());
}

template<class T, bool isConst, class A>
inline void
View<T, isConst, A>::updateSimplicity()
{
    // no invariant test here because this function
    // is called during unsafe updates of a view
    geometry_.updateSimplicity();
}

template<class T, bool isConst, class A>
inline const T&
View<T, isConst, A>::operator[]
(
    const size_t offset
)
const
{
    return data_[offset];
}

template<class T, bool isConst, class A>
inline T&
View<T, isConst, A>::operator[]
(
    const size_t offset
)
{
    return data_[offset];
}

template<class T, bool isConst, class A>
inline void
View<T, isConst, A>::testInvariant() const
{
    if(!MARRAY_NO_DEBUG) {
        if(geometry_.dimension() == 0) {
            marray_detail::Assert(geometry_.isSimple() == true);
            if(data_ != 0) { // scalar
                marray_detail::Assert(geometry_.size() == 1);
            }
        }
        else {
            marray_detail::Assert(data_ != 0);

            // test size_ to be consistent with shape_
            size_t testSize = 1;
            for(size_t j=0; j<geometry_.dimension(); ++j) {
                testSize *= geometry_.shape(j);
            }
            marray_detail::Assert(geometry_.size() == testSize);

            // test shapeStrides_ to be consistent with shape_
            if(geometry_.coordinateOrder() == FirstMajorOrder) {
                size_t tmp = 1;
                for(size_t j=0; j<geometry_.dimension(); ++j) {
                    marray_detail::Assert(geometry_.shapeStrides(geometry_.dimension()-j-1) == tmp);
                    tmp *= geometry_.shape(geometry_.dimension()-j-1);
                }
            }
            else {
                size_t tmp = 1;
                for(size_t j=0; j<geometry_.dimension(); ++j) {
                    marray_detail::Assert(geometry_.shapeStrides(j) == tmp);
                    tmp *= geometry_.shape(j);
                }
            }

            // test the simplicity condition
            if(geometry_.isSimple()) {
                for(size_t j=0; j<geometry_.dimension(); ++j) {
                    marray_detail::Assert(geometry_.strides(j) == geometry_.shapeStrides(j));
                }
            }
        }
    }
}

template<class T, bool isConst, class A>
template<class TLocal, bool isConstLocal, class ALocal>
inline bool View<T, isConst, A>::overlaps
(
    const View<TLocal, isConstLocal, ALocal>& v
) const
{
    testInvariant();
    if(!MARRAY_NO_ARG_TEST) {
        v.testInvariant();
    }
    if(data_ == 0 || v.data_ == 0) {
        return false;
    }
    else {
        const void* dataPointer_ = data_;
        const void* vDataPointer_ = v.data_;
        const void* maxPointer = & (*this)(this->size()-1);
        const void* maxPointerV = & v(v.size()-1);
        if(    (dataPointer_   <= vDataPointer_ && vDataPointer_ <= maxPointer)
            || (vDataPointer_ <= dataPointer_   && dataPointer_   <= maxPointerV) )
        {
            return true;
        }
    }
    return false;
}

template<class T, bool isConst, class A>
std::string
View<T, isConst, A>::asString
(
    const StringStyle& style
) const
{
    testInvariant();
    std::ostringstream out(std::ostringstream::out);
    if(style == MatrixStyle) {
        if(dimension() == 0) {
            // scalar
            out << "A = " << (*this)(0) << std::endl;
        }
        else if(dimension() == 1) {
            // vector
            out << "A = (";
            for(size_t j=0; j<this->size(); ++j) {
                out << (*this)(j) << ", ";
            }
            out << "\b\b)" << std::endl;
        }
        else if(dimension() == 2) {
            // matrix
            if(coordinateOrder() == FirstMajorOrder) {
                out << "A(r,c) =" << std::endl;
                for(size_t y=0; y<this->shape(0); ++y) {
                    for(size_t x=0; x<this->shape(1); ++x) {
                        out << (*this)(y, x) << ' ';
                    }
                    out << std::endl;
                }
            }
            else {
                out << "A(c,r) =" << std::endl;
                for(size_t y=0; y<this->shape(1); ++y) {
                    for(size_t x=0; x<this->shape(0); ++x) {
                        out << (*this)(x, y) << ' ';
                    }
                    out << std::endl;
                }
            }
        }
        else {
            // higher dimensional
            std::vector<size_t> c1(dimension());
            std::vector<size_t> c2(dimension());
            unsigned short q = 2;
            if(coordinateOrder() == FirstMajorOrder) {
                q = dimension()-3;
            }
            for(const_iterator it = this->begin(); it.hasMore(); ++it) {
                it.coordinate(c2.begin());
                if(it.index() == 0 || c2[q] != c1[q]) {
                    if(it.index() != 0) {
                        out << std::endl << std::endl;
                    }
                    if(coordinateOrder() == FirstMajorOrder) {
                        out << "A(";
                        for(size_t j=0; j<dimension()-2; ++j) {
                            out << c2[j] << ",";
                        }
                    }
                    else {
                        out << "A(c,r,";
                        for(size_t j=2; j<dimension(); ++j) {
                            out << c2[j] << ",";
                        }
                    }
                    out << '\b';
                    if(coordinateOrder() == FirstMajorOrder) {
                        out << ",r,c";
                    }
                    out << ") =" << std::endl;
                }
                else if(c2[1] != c1[1]) {
                    out << std::endl;
                }
                out << *it << " ";
                c1 = c2;
            }
            out << std::endl;
        }
        out << std::endl;
    }
    else if(style == TableStyle) {
        if(dimension() == 0) {
            // scalar
            out << "A = " << (*this)(0) << std::endl;
        }
        else {
            // non-scalar
            std::vector<size_t> c(dimension());
            for(const_iterator it = this->begin(); it.hasMore(); ++it) {
                out << "A(";
                it.coordinate(c.begin());
                for(size_t j=0; j<c.size(); ++j) {
                    out << c[j] << ',';
                }
                out << "\b) = " << *it << std::endl;
            }
        }
        out << std::endl;
    }
    return out.str();
}

// implementation of arithmetic operators of View

template<class T1, class T2, bool isConst, class A>
inline View<T1, false, A>&
operator+=
(
    View<T1, false, A>& v,
    const View<T2, isConst, A>& w
)
{
    marray_detail::operate(v, w, marray_detail::PlusEqual<T1, T2>());
    return v;
}

// prefix
template<class T, class A>
inline View<T, false, A>&
operator++
(
    View<T, false, A>& v
)
{
    marray_detail::operate(v, marray_detail::PrefixIncrement<T>());
    return v;
}

// postfix
template<class T, class A>
inline Marray<T, A>
operator++
(
    Marray<T, A>& in,
    int dummy
)
{
    Marray<T, A> out = in;
    marray_detail::operate(in, marray_detail::PostfixIncrement<T>());
    return out;
}

template<class T1, class T2, bool isConst, class A>
inline View<T1, false, A>&
operator-=
(
    View<T1, false, A>& v,
    const View<T2, isConst, A>& w
)
{
    marray_detail::operate(v, w, marray_detail::MinusEqual<T1, T2>());
    return v;
}

// prefix
template<class T, class A>
inline View<T, false, A>&
operator--
(
    View<T, false, A>& v
)
{
    marray_detail::operate(v, marray_detail::PrefixDecrement<T>());
    return v;
}

// postfix
template<class T, class A>
inline Marray<T, A>
operator--
(
    Marray<T, A>& in,
    int dummy
)
{
    Marray<T, A> out = in;
    marray_detail::operate(in, marray_detail::PostfixDecrement<T>());
    return out;
}

template<class T1, class T2, bool isConst, class A>
inline View<T1, false, A>&
operator*=
(
    View<T1, false, A>& v,
    const View<T2, isConst, A>& w
)
{
    marray_detail::operate(v, w, marray_detail::TimesEqual<T1, T2>());
    return v;
}

template<class T1, class T2, bool isConst, class A>
inline View<T1, false, A>&
operator/=
(
    View<T1, false, A>& v,
    const View<T2, isConst, A>& w
)
{
    marray_detail::operate(v, w, marray_detail::DividedByEqual<T1, T2>());
    return v;
}

template<class E1, class T1, class E2, class T2>
inline const BinaryViewExpression<E1, T1, E2, T2, marray_detail::Plus<T1, T2, typename marray_detail::PromoteType<T1, T2>::type> >
operator+
(
    const ViewExpression<E1, T1>& expression1,
    const ViewExpression<E2, T2>& expression2
)
{
    typedef typename marray_detail::PromoteType<T1, T2>::type promoted_type;
    typedef marray_detail::Plus<T1, T2, promoted_type> Functor;
    typedef BinaryViewExpression<E1, T1, E2, T2, Functor> return_type;
    return return_type(expression1, expression2);
}

template<class E, class T>
inline const ViewExpression<E,T>&
operator+
(
    const ViewExpression<E,T>& expression
) // unary
{
    return expression;
}

#define MARRAY_UNARY_OPERATOR(datatype, operation, functorname) \
template<class T, class A> \
inline View<T, false, A>& \
operator operation \
( \
    View<T, false, A>& v, \
    const datatype& x \
) \
{ \
    marray_detail::operate(v, static_cast<T>(x), marray_detail:: functorname <T, T>()); \
    return v; \
} \

#define MARRAY_UNARY_OPERATOR_ALL_TYPES(op, functorname) \
    MARRAY_UNARY_OPERATOR(char, op, functorname) \
    MARRAY_UNARY_OPERATOR(unsigned char, op, functorname) \
    MARRAY_UNARY_OPERATOR(short, op, functorname) \
    MARRAY_UNARY_OPERATOR(unsigned short, op, functorname) \
    MARRAY_UNARY_OPERATOR(int, op, functorname) \
    MARRAY_UNARY_OPERATOR(unsigned int, op, functorname) \
    MARRAY_UNARY_OPERATOR(long, op, functorname) \
    MARRAY_UNARY_OPERATOR(unsigned long, op, functorname) \
    MARRAY_UNARY_OPERATOR(float, op, functorname) \
    MARRAY_UNARY_OPERATOR(double, op, functorname) \
    MARRAY_UNARY_OPERATOR(long double, op, functorname) \

MARRAY_UNARY_OPERATOR_ALL_TYPES(+=, PlusEqual)
MARRAY_UNARY_OPERATOR_ALL_TYPES(-=, MinusEqual)
MARRAY_UNARY_OPERATOR_ALL_TYPES(*=, TimesEqual)
MARRAY_UNARY_OPERATOR_ALL_TYPES(/=, DividedByEqual)

template<class E1, class T1, class E2, class T2>
inline const BinaryViewExpression<E1, T1, E2, T2,
    marray_detail::Minus<T1, T2, typename marray_detail::PromoteType<T1, T2>::type> >
operator-
(
    const ViewExpression<E1, T1>& expression1,
    const ViewExpression<E2, T2>& expression2
)
{
    return BinaryViewExpression<E1, T1, E2, T2,
        marray_detail::Minus<T1, T2,
            typename marray_detail::PromoteType<T1, T2>::type> >(
            expression1, expression2);
}

template<class E, class T>
inline const UnaryViewExpression<E,T,marray_detail::Negate<T> >
operator-
(
    const ViewExpression<E,T>& expression
) // unary
{
    return UnaryViewExpression<E,T,marray_detail::Negate<T> >(
        expression);
}

template<class E1, class T1, class E2, class T2>
inline const BinaryViewExpression<E1, T1, E2, T2,
    marray_detail::Times<T1, T2, typename marray_detail::PromoteType<T1, T2>::type> >
operator*
(
    const ViewExpression<E1, T1>& expression1,
    const ViewExpression<E2, T2>& expression2
)
{
    return BinaryViewExpression<E1, T1, E2, T2,
        marray_detail::Times<T1, T2,
            typename marray_detail::PromoteType<T1, T2>::type > >(
            expression1, expression2);
}

template<class E1, class T1, class E2, class T2>
inline const BinaryViewExpression<E1, T1, E2, T2,
    marray_detail::DividedBy<T1, T2, typename marray_detail::PromoteType<T1, T2>::type> >
operator/
(
    const ViewExpression<E1, T1>& expression1,
    const ViewExpression<E2, T2>& expression2
)
{
    return BinaryViewExpression<E1, T1, E2, T2,
        marray_detail::DividedBy<T1, T2,
            typename marray_detail::PromoteType<T1, T2>::type > >(
            expression1, expression2);
}

#define MARRAY_BINARY_OPERATOR(datatype, operation, functorname) \
template<class E, class T> \
inline const BinaryViewExpressionScalarSecond< \
    E, T, datatype, marray_detail:: functorname < \
        T, datatype, typename marray_detail::PromoteType<T, datatype>::type \
    > \
> \
operator operation \
( \
    const ViewExpression<E, T>& expression, \
    const datatype& scalar \
) \
{ \
    typedef typename marray_detail::PromoteType<T, datatype>::type \
        promoted_type; \
    typedef marray_detail:: functorname <T, datatype, promoted_type> Functor; \
    typedef BinaryViewExpressionScalarSecond<E, T, datatype, Functor> \
        expression_type; \
    return expression_type(expression, scalar); \
} \
\
template<class E, class T> \
inline const BinaryViewExpressionScalarFirst \
< \
    E, T, datatype, marray_detail:: functorname < \
        datatype, T, typename marray_detail::PromoteType<datatype, T>::type \
    > \
> \
operator operation \
( \
    const datatype& scalar, \
    const ViewExpression<E, T>& expression \
) \
{ \
    typedef typename marray_detail::PromoteType<T, datatype>::type \
        promoted_type; \
    typedef marray_detail:: functorname <datatype, T, promoted_type> Functor; \
    typedef BinaryViewExpressionScalarFirst<E, T, datatype, Functor> \
        expression_type; \
    return expression_type(expression, scalar); \
}

#define MARRAY_BINARY_OPERATOR_ALL_TYPES(op, functorname) \
    MARRAY_BINARY_OPERATOR(char, op, functorname) \
    MARRAY_BINARY_OPERATOR(unsigned char, op, functorname) \
    MARRAY_BINARY_OPERATOR(short, op, functorname) \
    MARRAY_BINARY_OPERATOR(unsigned short, op, functorname) \
    MARRAY_BINARY_OPERATOR(int, op, functorname) \
    MARRAY_BINARY_OPERATOR(unsigned int, op, functorname) \
    MARRAY_BINARY_OPERATOR(long, op, functorname) \
    MARRAY_BINARY_OPERATOR(unsigned long, op, functorname) \
    MARRAY_BINARY_OPERATOR(float, op, functorname) \
    MARRAY_BINARY_OPERATOR(double, op, functorname) \
    MARRAY_BINARY_OPERATOR(long double, op, functorname) \

MARRAY_BINARY_OPERATOR_ALL_TYPES(+, Plus)
MARRAY_BINARY_OPERATOR_ALL_TYPES(-, Minus)
MARRAY_BINARY_OPERATOR_ALL_TYPES(*, Times)
MARRAY_BINARY_OPERATOR_ALL_TYPES(/, DividedBy)

// implementation of Marray

template<class T, class A>
inline void
Marray<T, A>::assign
(
    const allocator_type& allocator
)
{
    if(this->data_ != 0) {
        dataAllocator_.deallocate(this->data_, this->size());
        this->data_ = 0;
    }
    dataAllocator_ = allocator;
    base::assign();
}

template<class T, class A>
inline
Marray<T, A>::Marray
(
    const allocator_type& allocator
)
: base(allocator),
  dataAllocator_(allocator)
{
    testInvariant();
}

template<class T, class A>
inline
Marray<T, A>::Marray
(
    const T& value,
    const CoordinateOrder& coordinateOrder,
    const allocator_type& allocator
)
:   dataAllocator_(allocator)
{
    this->data_ = dataAllocator_.allocate(1);
    this->data_[0] = value;
    this->geometry_ = geometry_type(0, coordinateOrder, 1, true, allocator);
    testInvariant();
}

template<class T, class A>
inline
Marray<T, A>::Marray
(
    const Marray<T, A>& in
)
:   dataAllocator_(in.dataAllocator_)
{
    if(!MARRAY_NO_ARG_TEST) {
        in.testInvariant();
    }
    if(in.data_ == 0) {
        this->data_ = 0;
    }
    else {
        this->data_ = dataAllocator_.allocate(in.size());
        memcpy(this->data_, in.data_, (in.size())*sizeof(T));
    }
    this->geometry_ = in.geometry_;
    testInvariant();
}

template<class T, class A>
template<class TLocal, bool isConstLocal, class ALocal>
inline
Marray<T, A>::Marray
(
    const View<TLocal, isConstLocal, ALocal>& in
)
: dataAllocator_()
{
    if(!MARRAY_NO_ARG_TEST) {
        in.testInvariant();
    }

    // adapt geometry
    this->geometry_ = in.geometry_;
    for(size_t j=0; j<in.dimension(); ++j) {
        this->geometry_.strides(j) = in.geometry_.shapeStrides(j); // !
    }
    this->geometry_.isSimple() = true;

    // copy data
    if(in.size() == 0) {
        this->data_ = 0;
    }
    else {
        this->data_ = dataAllocator_.allocate(in.size());
    }
    if(in.isSimple() && marray_detail::IsEqual<T, TLocal>::type) {
        memcpy(this->data_, in.data_, (in.size())*sizeof(T));
    }
    else {
        typename View<TLocal, isConstLocal, ALocal>::const_iterator it = in.begin();
        for(size_t j=0; j<this->size(); ++j, ++it)  {
            this->data_[j] = static_cast<T>(*it);
        }
    }

    testInvariant();
}

template<class T, class A>
template<class E, class Te>
inline
Marray<T, A>::Marray
(
    const ViewExpression<E, Te>& expression,
    const allocator_type& allocator
)
:   dataAllocator_(allocator)
{
    this->data_ = dataAllocator_.allocate(expression.size());
    if(expression.dimension() == 0) {
        this->geometry_ = geometry_type(0,
            static_cast<const E&>(expression).coordinateOrder(),
            1, true, dataAllocator_);
    }
    else {
        this->geometry_ = geometry_type(
            static_cast<const E&>(expression).shapeBegin(),
            static_cast<const E&>(expression).shapeEnd(),
            static_cast<const E&>(expression).coordinateOrder(),
            static_cast<const E&>(expression).coordinateOrder(),
            dataAllocator_);

    }
    const E& e = static_cast<const E&>(expression);
    if(e.dimension() == 0) {
        marray_detail::Assert(MARRAY_NO_ARG_TEST || e.size() < 2);
        this->data_[0] = static_cast<T>(e(0));
    }
    else {
        marray_detail::Assert(MARRAY_NO_ARG_TEST || e.size() != 0);
        marray_detail::operate(*this, e, marray_detail::Assign<T, Te>());
    }
    testInvariant();
}

template<class T, class A>
template<class ShapeIterator>
inline
Marray<T, A>::Marray
(
    ShapeIterator begin,
    ShapeIterator end,
    const T& value,
    const CoordinateOrder& coordinateOrder,
    const allocator_type& allocator
)
: dataAllocator_(allocator)
{
    size_t size = std::accumulate(begin, end, static_cast<size_t>(1), std::multiplies<size_t>());
    marray_detail::Assert(MARRAY_NO_ARG_TEST || size != 0);
    base::assign(begin, end, dataAllocator_.allocate(size), coordinateOrder,
        coordinateOrder, allocator);
    for(size_t j=0; j<size; ++j) {
        this->data_[j] = value;
    }
    testInvariant();
}

template<class T, class A>
template<class ShapeIterator>
inline
Marray<T, A>::Marray
(
    const InitializationSkipping& is,
    ShapeIterator begin,
    ShapeIterator end,
    const CoordinateOrder& coordinateOrder,
    const allocator_type& allocator
)
: dataAllocator_(allocator)
{
    size_t size = std::accumulate(begin, end, 1, std::multiplies<size_t>());
    marray_detail::Assert(MARRAY_NO_ARG_TEST || size != 0);
    base::assign(begin, end, dataAllocator_.allocate(size), coordinateOrder,
        coordinateOrder, allocator);
    testInvariant();
}

#ifdef HAVE_CPP0X_INITIALIZER_LISTS
template<class T, class A>
inline
Marray<T, A>::Marray
(
    std::initializer_list<size_t> shape,
    const T& value,
    const CoordinateOrder& coordinateOrder,
    const allocator_type& allocator

)
: dataAllocator_(allocator)
{
    size_t size = std::accumulate(shape.begin(), shape.end(),
        1, std::multiplies<size_t>());
    marray_detail::Assert(MARRAY_NO_ARG_TEST || size != 0);
    base::assign(shape.begin(), shape.end(), dataAllocator_.allocate(size),
                 coordinateOrder, coordinateOrder, allocator);
    for(size_t j=0; j<size; ++j) {
        this->data_[j] = value;
    }
    testInvariant();
}
#endif

template<class T, class A>
inline
Marray<T, A>::~Marray()
{
    dataAllocator_.deallocate(this->data_, this->size());
}

template<class T, class A>

Marray<T, A>&
Marray<T, A>::operator=
(
    const Marray<T, A>& in
)
{
    testInvariant();
    if(!MARRAY_NO_ARG_TEST) {
        in.testInvariant();
    }
    if(this != &in) { // no self-assignment
        // copy data
        if(in.data_ == 0) { // un-initialized
            // free
            dataAllocator_.deallocate(this->data_, this->size());
            this->data_ = 0;
        }
        else {
            if(this->size() != in.size()) {
                // re-alloc
                dataAllocator_.deallocate(this->data_, this->size());
                this->data_ = dataAllocator_.allocate(in.size());
            }
            // copy data
            memcpy(this->data_, in.data_, in.size() * sizeof(T));
        }
        this->geometry_ = in.geometry_;
    }
    testInvariant();
    return *this;
}

template<class T, class A>
template<class TLocal, bool isConstLocal, class ALocal>
Marray<T, A>&
Marray<T, A>::operator=
(
    const View<TLocal, isConstLocal, ALocal>& in
)
{
    if(!MARRAY_NO_ARG_TEST) {
        in.testInvariant();
    }
    if( (void*)(this) != (void*)(&in) ) { // no self-assignment
        if(in.data_ == 0) {
            dataAllocator_.deallocate(this->data_, this->size());
            this->data_ = 0;
            this->geometry_ = in.geometry_;
        }
        else if(this->overlaps(in)) {
            Marray<T, A> m = in; // temporary copy
            (*this) = m;
        }
        else {
            // re-alloc memory if necessary
            if(this->size() != in.size()) {
                dataAllocator_.deallocate(this->data_, this->size());
                this->data_ = dataAllocator_.allocate(in.size());
            }

            // copy geometry
            this->geometry_.resize(in.dimension());
            for(size_t j=0; j<in.dimension(); ++j) {
                this->geometry_.shape(j) = in.geometry_.shape(j);
                this->geometry_.shapeStrides(j) = in.geometry_.shapeStrides(j);
                this->geometry_.strides(j) = in.geometry_.shapeStrides(j); // !
            }
            this->geometry_.size() = in.size();
            this->geometry_.isSimple() = true;
            this->geometry_.coordinateOrder() = in.coordinateOrder();

            // copy data
            if(in.isSimple() && marray_detail::IsEqual<T, TLocal>::type) {
                memcpy(this->data_, in.data_, (in.size())*sizeof(T));
            }
            else if(in.dimension() == 1)
                marray_detail::OperateHelperBinary<1, marray_detail::Assign<T, TLocal>, T, TLocal, isConstLocal, A, ALocal>::operate(*this, in, marray_detail::Assign<T, TLocal>(), this->data_, &in(0));
            else if(in.dimension() == 2)
                marray_detail::OperateHelperBinary<2, marray_detail::Assign<T, TLocal>, T, TLocal, isConstLocal, A, ALocal>::operate(*this, in, marray_detail::Assign<T, TLocal>(), this->data_, &in(0));
            else if(in.dimension() == 3)
                marray_detail::OperateHelperBinary<3, marray_detail::Assign<T, TLocal>, T, TLocal, isConstLocal, A, ALocal>::operate(*this, in, marray_detail::Assign<T, TLocal>(), this->data_, &in(0));
            else if(in.dimension() == 4)
                marray_detail::OperateHelperBinary<4, marray_detail::Assign<T, TLocal>, T, TLocal, isConstLocal, A, ALocal>::operate(*this, in, marray_detail::Assign<T, TLocal>(), this->data_, &in(0));
            else if(in.dimension() == 5)
                marray_detail::OperateHelperBinary<5, marray_detail::Assign<T, TLocal>, T, TLocal, isConstLocal, A, ALocal>::operate(*this, in, marray_detail::Assign<T, TLocal>(), this->data_, &in(0));
            else if(in.dimension() == 6)
                marray_detail::OperateHelperBinary<6, marray_detail::Assign<T, TLocal>, T, TLocal, isConstLocal, A, ALocal>::operate(*this, in, marray_detail::Assign<T, TLocal>(), this->data_, &in(0));
            else if(in.dimension() == 7)
                marray_detail::OperateHelperBinary<7, marray_detail::Assign<T, TLocal>, T, TLocal, isConstLocal, A, ALocal>::operate(*this, in, marray_detail::Assign<T, TLocal>(), this->data_, &in(0));
            else if(in.dimension() == 8)
                marray_detail::OperateHelperBinary<8, marray_detail::Assign<T, TLocal>, T, TLocal, isConstLocal, A, ALocal>::operate(*this, in, marray_detail::Assign<T, TLocal>(), this->data_, &in(0));
            else if(in.dimension() == 9)
                marray_detail::OperateHelperBinary<9, marray_detail::Assign<T, TLocal>, T, TLocal, isConstLocal, A, ALocal>::operate(*this, in, marray_detail::Assign<T, TLocal>(), this->data_, &in(0));
            else if(in.dimension() == 10)
                marray_detail::OperateHelperBinary<10, marray_detail::Assign<T, TLocal>, T, TLocal, isConstLocal, A, ALocal>::operate(*this, in, marray_detail::Assign<T, TLocal>(), this->data_, &in(0));
            else {
                typename View<TLocal, isConstLocal, ALocal>::const_iterator it = in.begin();
                for(size_t j=0; j<this->size(); ++j, ++it) {
                    this->data_[j] = static_cast<T>(*it);
                }
            }
        }
    }
    testInvariant();
    return *this;
}

template<class T, class A>
inline Marray<T, A>&
Marray<T, A>::operator=
(
    const T& value
)
{
    marray_detail::Assert(MARRAY_NO_DEBUG || this->data_ != 0);
    for(size_t j=0; j<this->size(); ++j) {
        this->data_[j] = value;
    }
    return *this;
}

template<class T, class A>
template<class E, class Te>
inline Marray<T, A>&
Marray<T, A>::operator=
(
    const ViewExpression<E, Te>& expression
)
{
    if(expression.overlaps(*this)) {
        Marray<T, A> m(expression); // temporary copy
        (*this) = m; // recursive call
    }
    else {
        // re-allocate memory (if necessary)
        if(this->size() != expression.size()) {
            dataAllocator_.deallocate(this->data_, this->size());
            this->data_ = dataAllocator_.allocate(expression.size());
        }

        // copy geometry
        this->geometry_.resize(expression.dimension());
        for(size_t j=0; j<expression.dimension(); ++j) {
            this->geometry_.shape(j) = expression.shape(j);
        }
        this->geometry_.size() = expression.size();
        this->geometry_.isSimple() = true;
        this->geometry_.coordinateOrder() = expression.coordinateOrder();
        if(this->geometry_.dimension() != 0) {
            marray_detail::stridesFromShape(this->geometry_.shapeBegin(), this->geometry_.shapeEnd(),
                this->geometry_.shapeStridesBegin(), this->geometry_.coordinateOrder());
            marray_detail::stridesFromShape(this->geometry_.shapeBegin(), this->geometry_.shapeEnd(),
                this->geometry_.stridesBegin(), this->geometry_.coordinateOrder());
        }

        // copy data
        marray_detail::operate(*this, expression, marray_detail::Assign<T, Te>());
    }
    return *this;
}

template<class T, class A>
template<bool SKIP_INITIALIZATION, class ShapeIterator>
inline void
Marray<T, A>::resizeHelper
(
    ShapeIterator begin,
    ShapeIterator end,
    const T& value
)
{
    testInvariant();
    // compute size
    std::vector<size_t> newShape;
    size_t newSize = 1;
    for(ShapeIterator it = begin; it != end; ++it) {
        size_t x = static_cast<size_t>(*it);
        marray_detail::Assert(MARRAY_NO_ARG_TEST || x > 0);
        newShape.push_back(x);
        newSize *= x;
    }
    // allocate new
    value_type* newData = dataAllocator_.allocate(newSize);
    if(!SKIP_INITIALIZATION) {
        for(size_t j=0; j<newSize; ++j) {
            newData[j] = value;
        }
    }
    // copy old data in region of overlap
    if(this->data_ != 0) {
        if(newSize == 1 || this->dimension() == 0) {
            newData[0] = this->data_[0];
        }
        else {
            std::vector<size_t> base1(this->dimension());
            std::vector<size_t> base2(newShape.size());
            std::vector<size_t> shape1(this->dimension(), 1);
            std::vector<size_t> shape2(newShape.size(), 1);
            for(size_t j=0; j<std::min(this->dimension(), newShape.size()); ++j) {
                shape1[j] = std::min(this->shape(j), newShape[j]);
                shape2[j] = shape1[j];
            }
            View<T, true, A> view1;
            this->constView(base1.begin(), shape1.begin(), view1);
            View<T, false, A> viewT(newShape.begin(), newShape.end(),
                newData, this->coordinateOrder(),
                this->coordinateOrder());
            View<T, false, A> view2;
            viewT.view(base2.begin(), shape2.begin(), view2);
            view1.squeeze();
            view2.squeeze();
            view2 = view1; // copy
        }
        dataAllocator_.deallocate(this->data_, this->size());
        this->data_ = 0;
    }
    base::assign(begin, end, newData, this->geometry_.coordinateOrder(),
        this->geometry_.coordinateOrder());
    testInvariant();
}

template<class T, class A>
template<class ShapeIterator>
void
Marray<T, A>::resize
(
    ShapeIterator begin,
    ShapeIterator end,
    const T& value
)
{
    if(std::distance(begin,end)==0) {
       if(this->size()>0) {
       *this=Marray<T,A>((*this).operator()(0));
       }
       else{
          *this=Marray<T,A>(T());
       }

    }
    else
    {
      resizeHelper<false>(begin, end, value);
    }
}

template<class T, class A>
template<class ShapeIterator>
void
Marray<T, A>::resize
(
    const InitializationSkipping& is,
    ShapeIterator begin,
    ShapeIterator end
)
{
    resizeHelper<true>(begin, end);
}

#ifdef HAVE_CPP0X_INITIALIZER_LISTS
template<class T, class A>
inline void
Marray<T, A>::resize
(
    std::initializer_list<size_t> shape,
    const T& value
)
{
    resizeHelper<false>(shape.begin(), shape.end(), value);
}

template<class T, class A>
inline void
Marray<T, A>::resize
(
    const InitializationSkipping& si,
    std::initializer_list<size_t> shape
)
{
    resizeHelper<true>(shape.begin(), shape.end());
}
#endif

template<class T, class A>
inline void
Marray<T, A>::testInvariant() const
{
    View<T, false, A>::testInvariant();
    marray_detail::Assert(MARRAY_NO_DEBUG || this->geometry_.isSimple());
}

// iterator implementation

template<class T, bool isConst, class A>
inline void
Iterator<T, isConst, A>::testInvariant() const
{
    if(!MARRAY_NO_DEBUG) {
        if(view_ == 0) {
            marray_detail::Assert(coordinates_.size() == 0
                && index_ == 0
                && pointer_ == 0);
        }
        else {
            if(view_->size() == 0) { // un-initialized view
                marray_detail::Assert(coordinates_.size() == 0
                    && index_ == 0
                    && pointer_ == 0);
            }
            else { // initialized view
                marray_detail::Assert(index_ <= view_->size());
                if(index_ == view_->size()) { // end iterator
                    marray_detail::Assert(pointer_ == &((*view_)(view_->size()-1)) + 1);
                }
                else {
                    marray_detail::Assert(pointer_ == &((*view_)(index_)));
                }
                if(!view_->isSimple()) {
                    marray_detail::Assert(coordinates_.size() == view_->dimension());
                    if(index_ == view_->size()) { // end iterator
                        if(view_->coordinateOrder() == LastMajorOrder) {
                            marray_detail::Assert(coordinates_[0] == view_->shape(0));
                            for(size_t j=1; j<coordinates_.size(); ++j) {
                                marray_detail::Assert(coordinates_[j] == view_->shape(j)-1);
                            }
                        }
                        else { // FirstMajorOrder
                            size_t d = view_->dimension() - 1;
                            marray_detail::Assert(coordinates_[d] == view_->shape(d));
                            for(size_t j=0; j<d; ++j) {
                                marray_detail::Assert(coordinates_[j] == view_->shape(j)-1);
                            }
                        }
                    }
                    else {
                        std::vector<size_t> testCoord(coordinates_.size());
                        view_->indexToCoordinates(index_, testCoord.begin());
                        for(size_t j=0; j<coordinates_.size(); ++j) {
                            marray_detail::Assert(coordinates_[j] == testCoord[j]);
                        }
                    }
                }
            }
        }
    }
}

template<class T, bool isConst, class A>
inline Iterator<T, isConst, A>::Iterator()
:   view_(0),
    pointer_(0),
    index_(0),
    coordinates_(std::vector<size_t>())
{
    testInvariant();
}

template<class T, bool isConst, class A>
inline Iterator<T, isConst, A>::Iterator
(
    const View<T, true, A>& view,
    const size_t index
)
:   view_(&view),
    pointer_(0),
    index_(index),
    coordinates_(std::vector<size_t>(view.dimension()))
    // Note for developers: If isConst==false, the construction view_(&view)
    // fails due to incompatible types. This is intended because it should
    // not be possible to construct a mutable iterator on constant data.
{
    if(view.size() == 0) { // un-initialized view
        marray_detail::Assert(MARRAY_NO_ARG_TEST || index == 0);
    }
    else {
        if(view.isSimple()) {
            marray_detail::Assert(MARRAY_NO_ARG_TEST || index <= view.size());
            pointer_ = &view(0) + index;
        }
        else {
            if(index >= view.size()) { // end iterator
                if(view_->coordinateOrder() == LastMajorOrder) {
                    coordinates_[0] = view.shape(0);
                    for(size_t j=1; j<view.dimension(); ++j) {
                        coordinates_[j] = view.shape(j)-1;
                    }
                }
                else { // FirstMajorOrder
                    size_t d = view_->dimension() - 1;
                    coordinates_[d] = view.shape(d);
                    for(size_t j=0; j<d; ++j) {
                        coordinates_[j] = view.shape(j)-1;
                    }
                }
                pointer_ = &view(view.size()-1) + 1;
            }
            else {
                view.indexToCoordinates(index, coordinates_.begin());
                pointer_ = &view(index);
            }
        }
    }
    testInvariant();
}

template<class T, bool isConst, class A>
inline Iterator<T, isConst, A>::Iterator
(
    const View<T, false, A>& view,
    const size_t index
)
:   view_(reinterpret_cast<view_pointer>(&view)),
    pointer_(0),
    index_(index),
    coordinates_(std::vector<size_t>(view.dimension()))
    // Note for developers: If isConst==true, the construction
    // view_(reinterpret_cast<view_pointer>(&view)) works as well.
    // This is intended because it should be possible to construct
    // a constant iterator on mutable data.
{
    if(view.size() == 0) { // un-initialized view
        marray_detail::Assert(MARRAY_NO_ARG_TEST || index == 0);
    }
    else {
        if(view.isSimple()) {
            marray_detail::Assert(MARRAY_NO_ARG_TEST || index <= view.size());
            pointer_ = &view(0) + index;
        }
        else {
            if(index >= view.size()) { // end iterator
                if(view_->coordinateOrder() == LastMajorOrder) {
                    coordinates_[0] = view.shape(0);
                    for(size_t j=1; j<view.dimension(); ++j) {
                        coordinates_[j] = view.shape(j)-1;
                    }
                }
                else { // FirstMajorOrder
                    size_t d = view_->dimension() - 1;
                    coordinates_[d] = view.shape(d);
                    for(size_t j=0; j<d; ++j) {
                        coordinates_[j] = view.shape(j)-1;
                    }
                }
                pointer_ = &view(view.size()-1) + 1;
            }
            else {
                view.indexToCoordinates(index, coordinates_.begin());
                pointer_ = &view(index);
            }
        }
    }
    testInvariant();
}

template<class T, bool isConst, class A>
inline Iterator<T, isConst, A>::Iterator
(
    View<T, false, A>& view,
    const size_t index
)
:   view_(reinterpret_cast<view_pointer>(&view)),
    pointer_(0),
    index_(index),
    coordinates_(std::vector<size_t>(view.dimension()))
    // Note for developers: If isConst==true, the construction
    // view_(reinterpret_cast<view_pointer>(&view)) works as well.
    // This is intended because it should be possible to construct
    // a constant iterator on mutable data.
{
    if(view.size() == 0) { // un-initialized view
        marray_detail::Assert(MARRAY_NO_ARG_TEST || index == 0);
    }
    else {
        if(view.isSimple()) {
            marray_detail::Assert(MARRAY_NO_ARG_TEST || index <= view.size());
            pointer_ = &view(0) + index;
        }
        else {
            if(index >= view.size()) { // end iterator
                if(view_->coordinateOrder() == LastMajorOrder) {
                    coordinates_[0] = view.shape(0);
                    for(size_t j=1; j<view.dimension(); ++j) {
                        coordinates_[j] = view.shape(j)-1;
                    }
                }
                else { // FirstMajorOrder
                    size_t d = view_->dimension() - 1;
                    coordinates_[d] = view.shape(d);
                    for(size_t j=0; j<d; ++j) {
                        coordinates_[j] = view.shape(j)-1;
                    }
                }
                pointer_ = &view(view.size()-1) + 1;
            }
            else {
                view.indexToCoordinates(index, coordinates_.begin());
                pointer_ = &view(index);
            }
        }
    }
    testInvariant();
}

template<class T, bool isConst, class A>
inline Iterator<T, isConst, A>::Iterator
(
    const Iterator<T, false, A>& in
)
:   view_(view_pointer(in.view_)),
    pointer_(pointer(in.pointer_)),
    index_(in.index_),
    coordinates_(in.coordinates_)
{
    testInvariant();
}

template<class T, bool isConst, class A>
inline typename Iterator<T, isConst, A>::reference
Iterator<T, isConst, A>::operator*() const
{
    marray_detail::Assert(MARRAY_NO_DEBUG || (view_ != 0 && index_ < view_->size()));
    return *pointer_;
}

template<class T, bool isConst, class A>
inline typename Iterator<T, isConst, A>::pointer
Iterator<T, isConst, A>::operator->() const
{
    marray_detail::Assert(MARRAY_NO_DEBUG || (view_ != 0 && index_ < view_->size()));
    return pointer_;
}

template<class T, bool isConst, class A>
inline typename Iterator<T, isConst, A>::reference
Iterator<T, isConst, A>::operator[]
(
    const size_t x
) const
{
    marray_detail::Assert(MARRAY_NO_DEBUG || (view_ != 0 && x+index_ < view_->size()));
    return (*view_)(x+index_);
}

template<class T, bool isConst, class A>
inline Iterator<T, isConst, A>&
Iterator<T, isConst, A>::operator+=
(
    const difference_type& x
)
{
    marray_detail::Assert(MARRAY_NO_DEBUG || view_ != 0);
    if(index_ < view_->size()) { // view initialized and iterator not at the end
        if(index_ + x < view_->size()) {
            index_ += x;
            if(view_->isSimple()) {
                pointer_ += x;
            }
            else {
                pointer_ = &((*view_)(index_));
                view_->indexToCoordinates(index_, coordinates_.begin());
            }
        }
        else {
            // set to end iterator
            index_ = view_->size();
            if(view_->isSimple()) {
                pointer_ = &(*view_)(0) + view_->size();
            }
            else {
                pointer_ = (&(*view_)(view_->size()-1)) + 1;
                view_->indexToCoordinates(view_->size()-1, coordinates_.begin());
                if(view_->coordinateOrder() == LastMajorOrder) {
                    ++coordinates_[0];
                }
                else { // FirstMajorOrder
                    ++coordinates_[view_->dimension()-1];
                }
            }
        }
    }
    testInvariant();
    return *this;
}

template<class T, bool isConst, class A>
inline Iterator<T, isConst, A>&
Iterator<T, isConst, A>::operator-=
(
    const difference_type& x
)
{
    marray_detail::Assert(MARRAY_NO_DEBUG || view_ != 0);
    marray_detail::Assert(MARRAY_NO_ARG_TEST || static_cast<difference_type>(index_) >= x);
    index_ -= x;
    if(view_->isSimple()) {
        pointer_ -= x;
    }
    else {
        pointer_ = &((*view_)(index_));
        view_->indexToCoordinates(index_, coordinates_.begin());
    }
    testInvariant();
    return *this;
}

template<class T, bool isConst, class A>
inline Iterator<T, isConst, A>&
Iterator<T, isConst, A>::operator++()
{
    marray_detail::Assert(MARRAY_NO_DEBUG || view_ != 0);
    if(index_ < view_->size()) { // view initialized and iterator not at the end
        ++index_;
        if(view_->isSimple()) {
            ++pointer_;
        }
        else {
            if(index_ < view_->size()) {
                if(view_->coordinateOrder() == LastMajorOrder) {
                    for(size_t j=0; j<coordinates_.size(); ++j) {
                        if(coordinates_[j] == view_->shape(j)-1) {
                            pointer_ -= view_->strides(j) * coordinates_[j];
                            coordinates_[j] = 0;
                        }
                        else {
                            pointer_ += view_->strides(j);
                            ++coordinates_[j];
                            break;
                        }
                    }
                }
                else { // FirstMajorOrder
                    size_t j = coordinates_.size() - 1;
                    for(;;) {
                        if(coordinates_[j] == view_->shape(j)-1) {
                            pointer_ -= view_->strides(j) * coordinates_[j];
                            coordinates_[j] = 0;
                        }
                        else {
                            pointer_ += view_->strides(j);
                            ++coordinates_[j];
                            break;
                        }
                        if(j == 0) {
                            break;
                        }
                        else {
                            --j;
                        }
                    }
                }
            }
            else {
                // set to end iterator
                pointer_ = &((*view_)(view_->size()-1)) + 1;
                if(view_->coordinateOrder() == LastMajorOrder) {
                    ++coordinates_[0];
                }
                else { // FirstMajorOrder
                    ++coordinates_[view_->dimension()-1];
                }
            }
        }
    }
    testInvariant();
    return *this;
}

template<class T, bool isConst, class A>
inline Iterator<T, isConst, A>&
Iterator<T, isConst, A>::operator--()
{
    marray_detail::Assert(MARRAY_NO_DEBUG || (view_ != 0 && index_ > 0));
    --index_;
    if(view_->isSimple()) {
        --pointer_;
    }
    else {
        if(index_ == view_->size()) {
            // decrement from end iterator
            --pointer_;
            if(view_->coordinateOrder() == LastMajorOrder) {
                --coordinates_[0];
            }
            else { // FirstMajorOrder
                --coordinates_[view_->dimension() - 1];
            }
        }
        else {
            if(view_->coordinateOrder() == LastMajorOrder) {
                for(size_t j=0; j<coordinates_.size(); ++j) {
                    if(coordinates_[j] == 0) {
                        coordinates_[j] = view_->shape(j)-1;
                        pointer_ += view_->strides(j) * coordinates_[j];
                    }
                    else {
                        pointer_ -= view_->strides(j);
                        --coordinates_[j];
                        break;
                    }
                }
            }
            else { // FirstMajorOrder
                size_t j = view_->dimension() - 1;
                for(;;) {
                    if(coordinates_[j] == 0) {
                        coordinates_[j] = view_->shape(j)-1;
                        pointer_ += view_->strides(j) * coordinates_[j];
                    }
                    else {
                        pointer_ -= view_->strides(j);
                        --coordinates_[j];
                        break;
                    }
                    if(j == 0) {
                        break;
                    }
                    else {
                        --j;
                    }
                }
            }
        }
    }
    testInvariant();
    return *this;
}

template<class T, bool isConst, class A>
inline Iterator<T, isConst, A>
Iterator<T, isConst, A>::operator++(int)
{
    marray_detail::Assert(MARRAY_NO_DEBUG || view_ != 0);
    Iterator<T, isConst, A> copy = *this;
    ++(*this);
    return copy;
}

template<class T, bool isConst, class A>
inline Iterator<T, isConst, A>
Iterator<T, isConst, A>::operator--(int)
{
    marray_detail::Assert(MARRAY_NO_DEBUG || (view_ != 0 && index_ > 0));
    Iterator<T, isConst, A> copy = *this;
    --(*this);
    return copy;
}

template<class T, bool isConst, class A>
inline Iterator<T, isConst, A>
Iterator<T, isConst, A>::operator+
(
    const difference_type& x
) const
{
    marray_detail::Assert(MARRAY_NO_DEBUG || view_ != 0);
    Iterator<T, isConst, A> tmp = *this;
    tmp += x;
    return tmp;
}

template<class T, bool isConst, class A>
inline Iterator<T, isConst, A>
Iterator<T, isConst, A>::operator-
(
    const difference_type& x
) const
{
    marray_detail::Assert(MARRAY_NO_DEBUG || view_ != 0);
    Iterator<T, isConst, A> tmp = *this;
    tmp -= x;
    return tmp;
}

template<class T, bool isConst, class A>
template<bool isConstLocal>
inline typename Iterator<T, isConst, A>::difference_type
Iterator<T, isConst, A>::operator-
(
    const Iterator<T, isConstLocal, A>& it
) const
{
    marray_detail::Assert(MARRAY_NO_DEBUG || view_ != 0);
    marray_detail::Assert(MARRAY_NO_ARG_TEST || it.view_ != 0);
    return difference_type(index_)-difference_type(it.index_);
}

template<class T, bool isConst, class A>
template<bool isConstLocal>
inline bool
Iterator<T, isConst, A>::operator==
(
    const Iterator<T, isConstLocal, A>& it
) const
{
    marray_detail::Assert(MARRAY_NO_DEBUG || view_ != 0);
    marray_detail::Assert(MARRAY_NO_ARG_TEST || (it.view_ != 0 && (void*)it.view_ == (void*)view_));
    return index_ == it.index_;
}

template<class T, bool isConst, class A>
template<bool isConstLocal>
inline bool
Iterator<T, isConst, A>::operator!=
(
    const Iterator<T, isConstLocal, A>& it
) const
{
    marray_detail::Assert(MARRAY_NO_DEBUG || view_ != 0);
    marray_detail::Assert(MARRAY_NO_ARG_TEST || it.view_ != 0);
    marray_detail::Assert(MARRAY_NO_ARG_TEST ||
        static_cast<const void*>(it.view_) == static_cast<const void*>(view_));
    return index_ != it.index_;
}

template<class T, bool isConst, class A>
template<bool isConstLocal>
inline bool
Iterator<T, isConst, A>::operator<
(
    const Iterator<T, isConstLocal, A>& it
) const
{
    marray_detail::Assert(MARRAY_NO_DEBUG || view_ != 0);
    marray_detail::Assert(MARRAY_NO_ARG_TEST || (it.view_ != 0 && it.view_ == view_));
    return(index_ < it.index_);
}

template<class T, bool isConst, class A>
template<bool isConstLocal>
inline bool
Iterator<T, isConst, A>::operator>
(
    const Iterator<T, isConstLocal, A>& it
) const
{
    marray_detail::Assert(MARRAY_NO_DEBUG || view_ != 0);
    marray_detail::Assert(MARRAY_NO_ARG_TEST || (it.view_ != 0 && it.view_ == view_));
    return(index_ > it.index_);
}

template<class T, bool isConst, class A>
template<bool isConstLocal>
inline bool
Iterator<T, isConst, A>::operator<=
(
    const Iterator<T, isConstLocal, A>& it
) const
{
    marray_detail::Assert(MARRAY_NO_DEBUG || view_ != 0);
    marray_detail::Assert(MARRAY_NO_ARG_TEST || (it.view_ != 0 && it.view_ == view_));
    return(index_ <= it.index_);
}

template<class T, bool isConst, class A>
template<bool isConstLocal>
inline bool
Iterator<T, isConst, A>::operator>=
(
    const Iterator<T, isConstLocal, A>& it
) const
{
    marray_detail::Assert(MARRAY_NO_DEBUG || view_ != 0);
    marray_detail::Assert(MARRAY_NO_ARG_TEST || (it.view_ != 0 && it.view_ == view_));
    return(index_ >= it.index_);
}

template<class T, bool isConst, class A>
inline bool
Iterator<T, isConst, A>::hasMore() const
{
    marray_detail::Assert(MARRAY_NO_DEBUG || view_ != 0);
    return index_ < view_->size();
}

template<class T, bool isConst, class A>
inline size_t
Iterator<T, isConst, A>::index() const
{
    return index_;
}

template<class T, bool isConst, class A>
template<class CoordinateIterator>
inline void
Iterator<T, isConst, A>::coordinate
(
    CoordinateIterator it
) const
{
    marray_detail::Assert(MARRAY_NO_DEBUG || view_ != 0);
    marray_detail::Assert(MARRAY_NO_ARG_TEST || index_ < view_->size());
    if(view_->isSimple()) {
        view_->indexToCoordinates(index_, it);
    }
    else {
        for(size_t j=0; j<coordinates_.size(); ++j, ++it) {
            *it = coordinates_[j];
        }
    }
}

// Vector implementation

template<class T, class A>
inline
Vector<T, A>::Vector
(
    const allocator_type& allocator
)
: base(allocator)
{
    testInvariant();
}

template<class T, class A>
template<class TLocal, bool isConstLocal, class ALocal>
inline
Vector<T, A>::Vector
(
    const View<TLocal, isConstLocal, ALocal>& in
)
{
    in.testInvariant();
    marray_detail::Assert(MARRAY_NO_ARG_TEST ||
        in.data_ == 0 // un-initialized
        || (in.dimension() == 0 && in.size() == 1) // scalar
        || in.dimension() == 1 // vector
    );
    this->geometry_.size() = in.size();
    this->geometry_.coordinateOrder() = in.coordinateOrder();
    if(in.data_ != 0) { // in is initialized
        this->geometry_.resize(1);
        this->geometry_.shape(0) = in.size();
        this->geometry_.shapeStrides(0) = 1;
        this->geometry_.strides(0) = 1;
        this->data_ = this->dataAllocator_.allocate(this->size());
        if(in.dimension() == 0) { // in is a scalar
            this->data_[0] = static_cast<T>(in(0));
        }
        else {
            if(in.isSimple() && marray_detail::IsEqual<T, TLocal>::type) {
                memcpy(this->data_, in.data_, (this->size())*sizeof(T));
            }
            else {
                for(size_t j=0; j<in.size(); ++j) {
                    this->data_[j] = static_cast<T>(in(j));
                }
            }
        }
    }
    testInvariant();
}

template<class T, class A>
inline
Vector<T, A>::Vector
(
    const size_t size,
    const T& value,
    const allocator_type& allocator
)
: base(allocator)
{
    if(size != 0) {
        size_t shape[1] = {size};
        this->data_ = this->dataAllocator_.allocate(size);
        base::base::assign(&shape[0], &shape[1], this->data_);
        for(size_t j=0; j<size; ++j) {
            this->data_[j] = value;
        }
    }
    testInvariant();
}

template<class T, class A>
inline
Vector<T, A>::Vector
(
    const InitializationSkipping& is,
    const size_t size,
    const allocator_type& allocator
)
: base(allocator)
{
    if(size != 0) {
        size_t shape[1] = {size};
        this->data_ = this->dataAllocator_.allocate(size);
        base::base::assign(&shape[0], &shape[1], this->data_);
    }
    testInvariant();
}

template<class T, class A>
template<class E, class Te>
inline
Vector<T, A>::Vector
(
    const ViewExpression<E, Te>& expression,
    const allocator_type& allocator
)
: base(expression, allocator)
{
    marray_detail::Assert(MARRAY_NO_ARG_TEST || expression.dimension() == 1);
    marray_detail::Assert(MARRAY_NO_ARG_TEST || expression.size() > 0);
    testInvariant();
}

#ifdef HAVE_CPP0X_INITIALIZER_LISTS
template<class T, class A>
inline
Vector<T, A>::Vector
(
    std::initializer_list<T> list,
    const allocator_type& allocator
)
: base(allocator)
{
    if(list.size() != 0) {
        size_t shape[1] = {list.size()};
        this->data_ = this->dataAllocator_.allocate(list.size());
        base::base::assign(&shape[0], &shape[1], this->data_);
        int j=0;
        for(const T *p = list.begin(); p != list.end(); ++p, ++j) {
            this->data_[j] = *p;
        }
    }
    testInvariant();
}
#endif

template<class T, class A>
inline Vector<T, A>&
Vector<T, A>::operator=
(
    const T& value
)
{
    marray_detail::Assert(MARRAY_NO_DEBUG || this->data_ != 0);
    for(size_t j=0; j<this->size(); ++j) {
        this->data_[j] = value;
    }
    return *this;
}

template<class T, class A>
inline Vector<T, A>&
Vector<T, A>::operator=
(
    const Vector<T, A>& in
)
{
    in.testInvariant();
    base::operator=(in);
    testInvariant();
    return *this;
}

template<class T, class A>
template<class TLocal, bool isConstLocal, class ALocal>
inline Vector<T, A>&
Vector<T, A>::operator=
(
    const View<TLocal, isConstLocal, ALocal>& in
)
{
    in.testInvariant();
    marray_detail::Assert(MARRAY_NO_ARG_TEST ||
        in.data_ == 0 || // empty
        (in.dimension() == 0 && in.size() == 1) || // scalar
        in.dimension() == 1); // vector
    if(in.geometry_.dimension() == 0 && in.geometry_.size() == 1) { // in is a View to a scalar
        // copy data
        if(this->size() != 1) {
            this->dataAllocator_.deallocate(this->data_, this->size());
            this->data_ = this->dataAllocator_.allocate(1);
        }
        this->data_[0] = static_cast<T>(in(0));
        this->geometry_.resize(1);
        this->geometry_.shape(0) = 1;
        this->geometry_.shapeStrides(0) = 1;
        this->geometry_.strides(0) = 1;
        this->geometry_.size() = 1;
        this->geometry_.isSimple() = true;
        this->geometry_.coordinateOrder() = in.coordinateOrder();
    }
    else {
        base::operator=(in);
    }
    testInvariant();
    return *this;
}

template<class T, class A>
template<class E, class Te>
inline Vector<T, A>&
Vector<T, A>::operator=
(
    const ViewExpression<E, Te>& expression
)
{
    marray_detail::Assert(MARRAY_NO_ARG_TEST || expression.dimension() == 1);
    marray_detail::Assert(MARRAY_NO_ARG_TEST || expression.size() > 0);
    base::operator=(expression);
    return *this;
}

template<class T, class A>
inline void
Vector<T, A>::reshape
(
    const size_t size
)
{
    marray_detail::Assert(size == this->size());
}

template<class T, class A>
inline void
Vector<T, A>::resize
(
    const size_t size,
    const T& value
)
{
    if(size == 0) {
        base::assign();
    }
    else {
        size_t shape[1] = {size};
        base::resize(&shape[0], &shape[1], value);
    }
    testInvariant();
}

template<class T, class A>
inline void
Vector<T, A>::resize
(
    const InitializationSkipping& is,
    const size_t size
)
{
    if(size == 0) {
        base::assign();
    }
    else {
        size_t shape[1] = {size};
        base::resize(is, &shape[0], &shape[1]);
    }
    testInvariant();
}

template<class T, class A>
inline T&
Vector<T, A>::operator[]
(
    const size_t index
)
{
    testInvariant();
    return this->operator()(index);
}

template<class T, class A>
inline const T&
Vector<T, A>::operator[]
(
    const size_t index
) const
{
    testInvariant();
    return this->operator()(index);
}

template<class T, class A>
inline void
Vector<T, A>::testInvariant() const
{
    View<T, false, A>::testInvariant();
    marray_detail::Assert(MARRAY_NO_DEBUG || this->data_ == 0 ||
        (this->geometry_.isSimple() && this->geometry_.dimension() == 1));
}

// Matrix implementation

template<class T, class A>
inline
Matrix<T, A>::Matrix
(
    const allocator_type& allocator
)
: base(allocator)
{
    testInvariant();
}

template<class T, class A>
    template<class TLocal, bool isConstLocal, class ALocal>
inline
Matrix<T, A>::Matrix
(
    const View<TLocal, isConstLocal, ALocal>& in
)
{
    in.testInvariant();
    marray_detail::Assert(MARRAY_NO_ARG_TEST ||
        in.data_ == 0 || // not initialized
        (in.dimension() == 0 && in.size() == 1) || // scalar
        in.dimension() == 2); // matrix
    this->geometry_.size() = in.size();
    this->geometry_.coordinateOrder() = in.coordinateOrder();
    if(in.data_ != 0) { // 'in' is uninitialized
        this->geometry_.resize(2);
        if(in.dimension() == 0) { // in is a scalar
            this->geometry_.shape(0) = 1;
            this->geometry_.shape(1) = 1;
            this->geometry_.shapeStrides(0) = 1;
            this->geometry_.shapeStrides(1) = 1;
            this->geometry_.strides(0) = 1;
            this->geometry_.strides(1) = 1;
            this->data_ = this->dataAllocator_.allocate(1);
            this->data_[0] = static_cast<T>(in(0));
        }
        else {
            this->geometry_.shape(0) = in.shape(0);
            this->geometry_.shape(1) = in.shape(1);
            this->geometry_.shapeStrides(0) = in.geometry_.shapeStrides(0);
            this->geometry_.shapeStrides(1) = in.geometry_.shapeStrides(1);
            this->geometry_.strides(0) = in.geometry_.shapeStrides(0); // !
            this->geometry_.strides(1) = in.geometry_.shapeStrides(1); // !
            this->data_ = this->dataAllocator_.allocate(this->size());
            if(in.isSimple() && marray_detail::IsEqual<T, TLocal>::type) {
                memcpy(this->data_, in.data_, (this->size())*sizeof(T));
            }
            else {
                for(size_t j=0; j<in.size(); ++j) {
                    this->data_[j] = static_cast<T>(in(j));
                }
            }
        }
    }
    testInvariant();
}

template<class T, class A>
inline
Matrix<T, A>::Matrix
(
    const size_t n1,
    const size_t n2,
    const T& value,
    const CoordinateOrder& coordinateOrder,
    const allocator_type& allocator
)
: base(allocator)
{
    if(n1 > 0 && n2 > 0) {
        size_t shape[2] = {n1, n2};
        T* data = this->dataAllocator_.allocate(n1 * n2);
        base::base::assign(&shape[0], &shape[2], data,
            coordinateOrder, coordinateOrder);
        for(size_t j=0; j<this->size(); ++j) {
            this->data_[j] = value;
        }
    }
    testInvariant();
}


template<class T, class A>
inline
Matrix<T, A>::Matrix
(
    const InitializationSkipping& is,
    const size_t n1,
    const size_t n2,
    const CoordinateOrder& coordinateOrder,
    const allocator_type& allocator
)
: base(allocator)
{
    if(n1 > 0 && n2 > 0) {
        size_t shape[2] = {n1, n2};
        this->data_ = this->dataAllocator_.allocate(n1 * n2);
        base::base::assign(&shape[0], &shape[2], this->data_,
            coordinateOrder, coordinateOrder);
    }
    testInvariant();
}

template<class T, class A>
template<class E, class Te>
inline
Matrix<T, A>::Matrix
(
    const ViewExpression<E, Te>& expression,
    const allocator_type& allocator
)
: base(expression, allocator)
{
    marray_detail::Assert(MARRAY_NO_ARG_TEST || expression.dimension() == 2);
    marray_detail::Assert(MARRAY_NO_ARG_TEST || expression.size() > 0);
    testInvariant();
}

template<class T, class A>
inline Matrix<T, A>&
Matrix<T, A>::operator=
(
    const T& value
)
{
    marray_detail::Assert(MARRAY_NO_DEBUG || this->data_ != 0);
    for(size_t j=0; j<this->size(); ++j) {
        this->data_[j] = value;
    }
    return *this;
}

template<class T, class A>
inline Matrix<T, A>&
Matrix<T, A>::operator=
(
    const Matrix<T, A>& in
)
{
    in.testInvariant();
    base::operator=(in);
    testInvariant();
    return *this;
}

template<class T, class A>
template<class TLocal, bool isConstLocal, class ALocal>
inline Matrix<T, A>&
Matrix<T, A>::operator=
(
    const View<TLocal, isConstLocal, ALocal>& in
)
{
    marray_detail::Assert(MARRAY_NO_ARG_TEST ||
        in.data_ != 0 || // empty
        (in.dimension() == 0 && in.size() == 1) || // scalar
        in.dimension() == 2);
    if(in.dimension() == 0 && in.size() == 1) { // in is a View to a scalar
        // copy data
        if(this->size() != 1) {
            this->dataAllocator_.deallocate(this->data_, this->size());
            this->data_ = this->dataAllocator_.allocate(1);
        }
        this->data_[0] = static_cast<T>(in(0));
        this->geometry_.resize(2);
        this->geometry_.shape(0) = 1;
        this->geometry_.shape(1) = 1;
        this->geometry_.shapeStrides(0) = 1;
        this->geometry_.shapeStrides(1) = 1;
        this->geometry_.strides(0) = 1;
        this->geometry_.strides(1) = 1;
        this->geometry_.size() = 1;
        this->geometry_.isSimple() = true;
        this->geometry_.coordinateOrder() = in.coordinateOrder();
    }
    else {
        base::operator=(in);
    }
    testInvariant();
    return *this;
}

template<class T, class A>
template<class E, class Te>
inline Matrix<T, A>&
Matrix<T, A>::operator=
(
    const ViewExpression<E, Te>& expression
)
{
    marray_detail::Assert(MARRAY_NO_ARG_TEST || expression.dimension() == 2);
    marray_detail::Assert(MARRAY_NO_ARG_TEST || expression.size() > 0);
    base::operator=(expression);
    return *this;
}

template<class T, class A>
inline void
Matrix<T, A>::resize
(
    const size_t n1,
    const size_t n2,
    const T& value
)
{
    if(n1 == 0 || n2 == 0) {
        base::assign();
    }
    else {
        size_t shape[2] = {n1, n2};
        base::resize(&shape[0], &shape[2], value);
    }
    testInvariant();
}

template<class T, class A>
inline void
Matrix<T, A>::resize
(
    const InitializationSkipping& is,
    const size_t n1,
    const size_t n2
)
{
    if(n1 == 0 || n2 == 0) {
        base::assign();
    }
    else {
        size_t shape[2] = {n1, n2};
        base::resize(is, &shape[0], &shape[2]);
    }
    testInvariant();
}

template<class T, class A>
inline void
Matrix<T, A>::reshape
(
    const size_t n1,
    const size_t n2
)
{
    marray_detail::Assert(MARRAY_NO_ARG_TEST || (n2 > 0 && n1 > 0));
    size_t shape[2] = {n1, n2};
    base::reshape(&shape[0], &shape[2]);
    testInvariant();
}

template<class T, class A>
inline void
Matrix<T, A>::testInvariant() const
{
    View<T, false, A>::testInvariant();
    marray_detail::Assert(MARRAY_NO_DEBUG || this->data_ == 0 ||
        (this->geometry_.isSimple() && this->geometry_.dimension() == 2));
}

// expression templates

template<class E, class T>
class ViewExpression {
public:
    typedef E expression_type;
    typedef T value_type;

    const size_t dimension() const
        { return static_cast<const E&>(*this).dimension(); }
    const size_t size() const
        { return static_cast<const E&>(*this).size(); }
    const size_t shape(const size_t j) const
        { return static_cast<const E&>(*this).shape(j); }
    const size_t* shapeBegin() const
        { return static_cast<const E&>(*this).shapeBegin(); }
    const size_t* shapeEnd() const
        { return static_cast<const E&>(*this).shapeEnd(); }
    template<class Tv, bool isConst, class A>
        bool overlaps(const View<Tv, isConst, A>& v) const
            { return static_cast<const E&>(*this).overlaps(v); }
    const CoordinateOrder& coordinateOrder() const
        { return static_cast<const E&>(*this).coordinateOrder(); }
    const bool isSimple() const
        { return static_cast<const E&>(*this).isSimple(); }
    template<class Accessor>
        const T& operator()(Accessor it) const
            { return static_cast<const E&>(*this)(it); }
    const T& operator()(const size_t c0, const size_t c1) const
        { return static_cast<const E&>(*this)(c0, c1); }
    const T& operator()(const size_t c0, const size_t c1, const size_t c2) const
        { return static_cast<const E&>(*this)(c0, c1, c2); }
    const T& operator()(const size_t c0, const size_t c1, const size_t c2, const size_t c3) const
        { return static_cast<const E&>(*this)(c0, c1, c2, c3); }
    const T& operator()(const size_t c0, const size_t c1, const size_t c2, const size_t c3, const size_t c4) const
        { return static_cast<const E&>(*this)(c0, c1, c2, c3, c4); }
    const T& operator[](const size_t offset) const
        { return static_cast<const E&>(*this)[offset]; }
    operator E&()
        { return static_cast<E&>(*this); }
    operator E const&() const
        { return static_cast<const E&>(*this); }

    // \cond suppress doxygen
    class ExpressionIterator {
    public:
        ExpressionIterator(const ViewExpression<E, T>& expression)
        : expression_(expression), // cast!
          data_(&expression_(0)),
          offset_(0)
            {}
        void incrementCoordinate(const size_t coordinateIndex)
            { offset_ += expression_.strides(coordinateIndex); }
        void resetCoordinate(const size_t coordinateIndex)
            { offset_ -= expression_.strides(coordinateIndex)

                         * (expression_.shape(coordinateIndex) - 1); }
        const T& operator*() const
            { // return expression_[offset_];
              // would require making this nested class a friend of View
              // which in turn would require a forward declaration of
              // this class. work around:
              return data_[offset_]; }
    private:
        const E& expression_;
        const T* data_;
        size_t offset_;
    };
    // \endcond suppress doxygen
};

// \cond suppress doxygen
template<class E, class T, class UnaryFunctor>
class UnaryViewExpression
: public ViewExpression<UnaryViewExpression<E, T, UnaryFunctor>, T>
{
public:
    typedef E expression_type;
    typedef T value_type;

    UnaryViewExpression(const ViewExpression<E, T>& e)
        : e_(e), // cast!
          unaryFunctor_(UnaryFunctor())
        {}
    const size_t dimension() const
        { return e_.dimension(); }
    const size_t size() const
        { return e_.size(); }
    const size_t shape(const size_t j) const
        { return e_.shape(j); }
    const size_t* shapeBegin() const
        { return e_.shapeBegin(); }
    const size_t* shapeEnd() const
        { return e_.shapeEnd(); }
    template<class Tv, bool isConst, class A>
        bool overlaps(const View<Tv, isConst, A>& v) const
            { return e_.overlaps(v); }
    const CoordinateOrder& coordinateOrder() const
        { return e_.coordinateOrder(); }
    const bool isSimple() const
        { return e_.isSimple(); }
    template<class Accessor>
        const T operator()(Accessor it) const
            { return unaryFunctor_(e_(it)); }
    const T operator()(const size_t c0, const size_t c1) const
        { return unaryFunctor_(e_(c0, c1)); }
    const T operator()(const size_t c0, const size_t c1,const size_t c2) const
        { return unaryFunctor_(e_(c0, c1, c2)); }
    const T operator()(const size_t c0, const size_t c1,const size_t c2, const size_t c3) const
        { return unaryFunctor_(e_(c0, c1, c2, c3)); }
    const T operator()(const size_t c0, const size_t c1,const size_t c2, const size_t c3, const size_t c4) const
        { return unaryFunctor_(e_(c0, c1, c2, c3, c4)); }
    const T operator[](const size_t offset) const
        { return unaryFunctor_(e_[offset]); }

    class ExpressionIterator {
    public:
        ExpressionIterator(const UnaryViewExpression<E, T, UnaryFunctor>& expression)
        : unaryFunctor_(expression.unaryFunctor_),
          iterator_(expression.e_)
            {}
        void incrementCoordinate(const size_t coordinateIndex)
            { iterator_.incrementCoordinate(coordinateIndex); }
        void resetCoordinate(const size_t coordinateIndex)
            { iterator_.resetCoordinate(coordinateIndex); }
        const T operator*() const
            { return unaryFunctor_(*iterator_); }
    private:
        UnaryFunctor unaryFunctor_;
        typename E::ExpressionIterator iterator_;
    };

private:
    const E& e_;
    UnaryFunctor unaryFunctor_;
};

template<class E1, class T1, class E2, class T2, class BinaryFunctor>
class BinaryViewExpression
: public ViewExpression<BinaryViewExpression<E1, T1, E2, T2, BinaryFunctor>,
                        typename marray_detail::PromoteType<T1, T2>::type>
{
public:
    typedef E1 expression_type_1;
    typedef E2 expression_type_2;
    typedef T1 value_type_1;
    typedef T2 value_type_2;
    typedef typename marray_detail::PromoteType<T1, T2>::type value_type;
    typedef BinaryFunctor functor_type;
    typedef ViewExpression<BinaryViewExpression<E1, T1, E2, T2, BinaryFunctor>,
        value_type> base;

    BinaryViewExpression(const ViewExpression<E1, T1>& e1,
        const ViewExpression<E2, T2>& e2)
        : e1_(e1), e2_(e2), // cast!
          binaryFunctor_(BinaryFunctor())
        {
            if(!MARRAY_NO_DEBUG) {
                marray_detail::Assert(e1_.size() != 0 && e2_.size() != 0);
                marray_detail::Assert(e1_.dimension() == e2_.dimension());
                for(size_t j=0; j<e1_.dimension(); ++j) {
                    marray_detail::Assert(e1_.shape(j) == e2_.shape(j));
                }
            }
        }
    const size_t dimension() const
        { return e1_.dimension() < e2_.dimension() ? e2_.dimension() : e1_.dimension(); }
    const size_t size() const
        { return e1_.size() < e2_.size() ? e2_.size() : e1_.size(); }
    const size_t shape(const size_t j) const
        { return e1_.dimension() < e2_.dimension() ? e2_.shape(j) : e1_.shape(j); }
    const size_t* shapeBegin() const
        { return e1_.dimension() < e2_.dimension() ? e2_.shapeBegin() : e1_.shapeBegin(); }
    const size_t* shapeEnd() const
        { return e1_.dimension() < e2_.dimension() ? e2_.shapeEnd() : e1_.shapeEnd(); }
    template<class Tv, bool isConst, class A>
        bool overlaps(const View<Tv, isConst, A>& v) const
            { return e1_.overlaps(v) || e2_.overlaps(v); }
    const CoordinateOrder& coordinateOrder() const
        { return e1_.coordinateOrder(); }
    const bool isSimple() const
        { return e1_.isSimple() && e2_.isSimple()
                 && e1_.coordinateOrder() == e2_.coordinateOrder(); }
    template<class Accessor>
        const value_type operator()(Accessor it) const
            { return binaryFunctor_(e1_(it), e2_(it)); }
    const value_type operator()(const size_t c0, const size_t c1) const
        { return binaryFunctor_(e1_(c0, c1), e2_(c0, c1)); }
    const value_type operator()(const size_t c0, const size_t c1, const size_t c2) const
        { return binaryFunctor_(e1_(c0, c1, c2), e2_(c0, c1, c2)); }
    const value_type operator()(const size_t c0, const size_t c1, const size_t c2, const size_t c3) const
        { return binaryFunctor_(e1_(c0, c1, c2, c3), e2_(c0, c1, c2, c3)); }
    const value_type operator()(const size_t c0, const size_t c1, const size_t c2, const size_t c3, const size_t c4) const
        { return binaryFunctor_(e1_(c0, c1, c2, c3, c4), e2_(c0, c1, c2, c3, c4)); }
    const value_type operator[](const size_t offset) const
        { return binaryFunctor_(e1_[offset], e2_[offset]); }

    class ExpressionIterator {
    public:
        ExpressionIterator(const BinaryViewExpression<E1, T1, E2, T2, BinaryFunctor>& expression)
        : binaryFunctor_(expression.binaryFunctor_),
          iterator1_(expression.e1_),
          iterator2_(expression.e2_)
            {}
        void incrementCoordinate(const size_t coordinateIndex)
            {   iterator1_.incrementCoordinate(coordinateIndex);
                iterator2_.incrementCoordinate(coordinateIndex); }
        void resetCoordinate(const size_t coordinateIndex)
            {   iterator1_.resetCoordinate(coordinateIndex);
                iterator2_.resetCoordinate(coordinateIndex); }
        const value_type operator*() const
            { return binaryFunctor_(*iterator1_, *iterator2_); }
    private:
        BinaryFunctor binaryFunctor_;
        typename E1::ExpressionIterator iterator1_;
        typename E2::ExpressionIterator iterator2_;
    };

private:
    const expression_type_1& e1_;
    const expression_type_2& e2_;
    BinaryFunctor binaryFunctor_;
};

template<class E, class T, class S, class BinaryFunctor>
class BinaryViewExpressionScalarFirst
: public ViewExpression<BinaryViewExpressionScalarFirst<E, T, S, BinaryFunctor>,
                        typename marray_detail::PromoteType<T, S>::type> {
public:
    typedef E expression_type;
    typedef T value_type_1;
    typedef S scalar_type;
    typedef typename marray_detail::PromoteType<T, S>::type value_type;
    typedef BinaryFunctor functor_type;
    typedef ViewExpression<BinaryViewExpressionScalarFirst<E, T, S, BinaryFunctor>,
        value_type> base;

    BinaryViewExpressionScalarFirst(const ViewExpression<E, T>& e,
        const scalar_type& scalar)
        : e_(e), // cast!
          scalar_(scalar), binaryFunctor_(BinaryFunctor())
        { }
    const size_t dimension() const
        { return e_.dimension(); }
    const size_t size() const
        { return e_.size(); }
    const size_t shape(const size_t j) const
        { return e_.shape(j); }
    const size_t* shapeBegin() const
        { return e_.shapeBegin(); }
    const size_t* shapeEnd() const
        { return e_.shapeEnd(); }
    template<class Tv, bool isConst, class A>
        bool overlaps(const View<Tv, isConst, A>& v) const
            { return e_.overlaps(v); }
    const CoordinateOrder& coordinateOrder() const
        { return e_.coordinateOrder(); }
    const bool isSimple() const
        { return e_.isSimple(); }
    template<class Accessor>
        const value_type operator()(Accessor it) const
            { return binaryFunctor_(scalar_, e_(it)); }
    const value_type operator()(const size_t c0, const size_t c1) const
        { return binaryFunctor_(scalar_, e_(c0, c1)); }
    const value_type operator()(const size_t c0, const size_t c1, const size_t c2) const
        { return binaryFunctor_(scalar_, e_(c0, c1, c2)); }
    const value_type operator()(const size_t c0, const size_t c1, const size_t c2, const size_t c3) const
        { return binaryFunctor_(scalar_, e_(c0, c1, c2, c3)); }
    const value_type operator()(const size_t c0, const size_t c1, const size_t c2, const size_t c3, const size_t c4) const
        { return binaryFunctor_(scalar_, e_(c0, c1, c2, c3, c4)); }
    const value_type operator[](const size_t offset) const
        { return binaryFunctor_(scalar_, e_[offset]); }

    class ExpressionIterator {
    public:
        ExpressionIterator(const BinaryViewExpressionScalarFirst<E, T, S, BinaryFunctor>& expression)
        : binaryFunctor_(expression.binaryFunctor_),
          scalar_(expression.scalar_),
          iterator_(expression.e_)
            {}
        void incrementCoordinate(const size_t coordinateIndex)
            { iterator_.incrementCoordinate(coordinateIndex); }
        void resetCoordinate(const size_t coordinateIndex)
            { iterator_.resetCoordinate(coordinateIndex); }
        const T operator*() const
            { return binaryFunctor_(scalar_, *iterator_); }
    private:
        BinaryFunctor binaryFunctor_;
        const typename BinaryViewExpressionScalarFirst<E, T, S, BinaryFunctor>::scalar_type& scalar_;
        typename E::ExpressionIterator iterator_;
    };

private:
    const expression_type& e_;
    const scalar_type scalar_;
    BinaryFunctor binaryFunctor_;
};

template<class E, class T, class S, class BinaryFunctor>
class BinaryViewExpressionScalarSecond
: public ViewExpression<BinaryViewExpressionScalarSecond<E, T, S, BinaryFunctor>,
                        typename marray_detail::PromoteType<T, S>::type> {
public:
    typedef T value_type_1;
    typedef E expression_type;
    typedef S scalar_type;
    typedef typename marray_detail::PromoteType<T, S>::type value_type;
    typedef BinaryFunctor functor_type;
    typedef ViewExpression<BinaryViewExpressionScalarSecond<E, T, S, BinaryFunctor>,
        value_type> base;

    BinaryViewExpressionScalarSecond(const ViewExpression<E, T>& e,
        const scalar_type& scalar)
        : e_(e), // cast!
          scalar_(scalar), binaryFunctor_(BinaryFunctor())
        { }
    const size_t dimension() const
        { return e_.dimension(); }
    const size_t size() const
        { return e_.size(); }
    const size_t shape(const size_t j) const
        { return e_.shape(j); }
    const size_t* shapeBegin() const
        { return e_.shapeBegin(); }
    const size_t* shapeEnd() const
        { return e_.shapeEnd(); }
    template<class Tv, bool isConst, class A>
        bool overlaps(const View<Tv, isConst, A>& v) const
            { return e_.overlaps(v); }
    const CoordinateOrder& coordinateOrder() const
        { return e_.coordinateOrder(); }
    const bool isSimple() const
        { return e_.isSimple(); }
    template<class Accessor>
        const value_type operator()(Accessor it) const
            { return binaryFunctor_(e_(it), scalar_); }
    const value_type operator()(const size_t c0, const size_t c1) const
        { return binaryFunctor_(e_(c0, c1), scalar_); }
    const value_type operator()(const size_t c0, const size_t c1, const size_t c2) const
        { return binaryFunctor_(e_(c0, c1, c2), scalar_); }
    const value_type operator()(const size_t c0, const size_t c1, const size_t c2, const size_t c3) const
        { return binaryFunctor_(e_(c0, c1, c2, c3), scalar_); }
    const value_type operator()(const size_t c0, const size_t c1, const size_t c2, const size_t c3, const size_t c4) const
        { return binaryFunctor_(e_(c0, c1, c2, c3, c4), scalar_); }
    const value_type operator[](const size_t offset) const
        { return binaryFunctor_(e_[offset], scalar_); }

    class ExpressionIterator {
    public:
        ExpressionIterator(const BinaryViewExpressionScalarSecond<E, T, S, BinaryFunctor>& expression)
        : binaryFunctor_(expression.binaryFunctor_),
          scalar_(expression.scalar_),
          iterator_(expression.e_)
            {}
        void incrementCoordinate(const size_t coordinateIndex)
            { iterator_.incrementCoordinate(coordinateIndex); }
        void resetCoordinate(const size_t coordinateIndex)
            { iterator_.resetCoordinate(coordinateIndex); }
        const T operator*() const
            { return binaryFunctor_(*iterator_, scalar_); }
    private:
        BinaryFunctor binaryFunctor_;
        const typename BinaryViewExpressionScalarSecond<E, T, S, BinaryFunctor>::scalar_type& scalar_;
        typename E::ExpressionIterator iterator_;
    };

private:
    const expression_type& e_;
    const scalar_type scalar_;
    BinaryFunctor binaryFunctor_;
};
// \endcond suppress doxygen

// implementation of marray_detail

// \cond suppress doxygen
namespace marray_detail {

template<class A>
class Geometry
{
public:
    typedef typename A::template rebind<size_t>::other allocator_type;

    Geometry(const allocator_type& = allocator_type());
    Geometry(const size_t, const CoordinateOrder& = defaultOrder,
        const size_t = 0, const bool = true,
        const allocator_type& = allocator_type());
    template<class ShapeIterator>
        Geometry(ShapeIterator, ShapeIterator,
            const CoordinateOrder& = defaultOrder,
            const CoordinateOrder& = defaultOrder,
            const allocator_type& = allocator_type());
    template<class ShapeIterator, class StrideIterator>
        Geometry(ShapeIterator, ShapeIterator, StrideIterator,
            const CoordinateOrder& = defaultOrder,
            const allocator_type& = allocator_type());
    Geometry(const Geometry<A>&);
    ~Geometry();

    Geometry<A>& operator=(const Geometry<A>&);

    void resize(const size_t dimension);
    const size_t dimension() const;
    const size_t shape(const size_t) const;
    size_t& shape(const size_t);
    const size_t shapeStrides(const size_t) const;
    size_t& shapeStrides(const size_t);
    const size_t strides(const size_t) const;
    size_t& strides(const size_t);
    const size_t* shapeBegin() const;
    size_t* shapeBegin();
    const size_t* shapeEnd() const;
    size_t* shapeEnd();
    const size_t* shapeStridesBegin() const;
    size_t* shapeStridesBegin();
    const size_t* shapeStridesEnd() const;
    size_t* shapeStridesEnd();
    const size_t* stridesBegin() const;
    size_t* stridesBegin();
    const size_t* stridesEnd() const;
    size_t* stridesEnd();
    const size_t size() const;
    size_t& size();
    const CoordinateOrder& coordinateOrder() const;
    CoordinateOrder& coordinateOrder();
    const bool isSimple() const;
    void updateSimplicity();
    bool& isSimple();

private:
    allocator_type allocator_;
    size_t* shape_;
    size_t* shapeStrides_;
        // Intended redundancy: shapeStrides_ could be
        // computed from shape_ and coordinateOrder_
    size_t* strides_;
    size_t dimension_;
    size_t size_;
        // intended redundancy: size_ could be computed from shape_
    CoordinateOrder coordinateOrder_;
    bool isSimple_;
        // simple array: an array which is unstrided (i.e. the strides
        // equal the shape strides), cf. the function testInvariant of
        // View for the formal definition.
};

template<class A>
inline
Geometry<A>::Geometry
(
    const typename Geometry<A>::allocator_type& allocator
)
: allocator_(allocator),
  shape_(0),
  shapeStrides_(0),
  strides_(0),
  dimension_(0),
  size_(0),
  coordinateOrder_(defaultOrder),
  isSimple_(true)
{
}

template<class A>
inline
Geometry<A>::Geometry
(
    const Geometry<A>& g
)
: allocator_(g.allocator_),
  shape_(g.dimension_==0 ? 0 : allocator_.allocate(g.dimension_*3)),
  shapeStrides_(shape_ + g.dimension_),
  strides_(shapeStrides_ + g.dimension_),
  dimension_(g.dimension_),
  size_(g.size_),
  coordinateOrder_(g.coordinateOrder_),
  isSimple_(g.isSimple_)
{
    /*
    for(size_t j=0; j<dimension_; ++j) {
        shape_[j] = g.shape_[j];
        shapeStrides_[j] = g.shapeStrides_[j];
        strides_[j] = g.strides_[j];
    }
    */
    memcpy(shape_, g.shape_, (dimension_*3)*sizeof(size_t));
}

template<class A>
inline
Geometry<A>::Geometry
(
    const size_t dimension,
    const CoordinateOrder& order,
    const size_t size,
    const bool isSimple,
    const typename Geometry<A>::allocator_type& allocator
)
: allocator_(allocator),
  shape_(allocator_.allocate(dimension*3)),
  shapeStrides_(shape_+dimension),
  strides_(shapeStrides_+dimension),
  dimension_(dimension),
  size_(size),
  coordinateOrder_(order),
  isSimple_(isSimple)
{
}

template<class A>
template<class ShapeIterator>
inline
Geometry<A>::Geometry
(
    ShapeIterator begin,
    ShapeIterator end,
    const CoordinateOrder& externalCoordinateOrder,
    const CoordinateOrder& internalCoordinateOrder,
    const typename Geometry<A>::allocator_type& allocator
)
: allocator_(allocator),
  shape_(allocator_.allocate(std::distance(begin, end) * 3)),
  shapeStrides_(shape_ + std::distance(begin, end)),
  strides_(shapeStrides_ + std::distance(begin, end)),
  dimension_(std::distance(begin, end)),
  size_(1),
  coordinateOrder_(internalCoordinateOrder),
  isSimple_(true)
{
    if(dimension_ != 0) { // if the array is not a scalar
        isSimple_ = (externalCoordinateOrder == internalCoordinateOrder);
        for(size_t j=0; j<dimension(); ++j, ++begin) {
            const size_t s = static_cast<size_t>(*begin);
            shape(j) = s;
            size() *= s;
        }
        stridesFromShape(shapeBegin(), shapeEnd(), stridesBegin(),
            externalCoordinateOrder);
        stridesFromShape(shapeBegin(), shapeEnd(), shapeStridesBegin(),
            internalCoordinateOrder);
    }
}

template<class A>
template<class ShapeIterator, class StrideIterator>
inline
Geometry<A>::Geometry
(
    ShapeIterator begin,
    ShapeIterator end,
    StrideIterator it,
    const CoordinateOrder& internalCoordinateOrder,
    const typename Geometry<A>::allocator_type& allocator
)
: allocator_(allocator),
  shape_(allocator_.allocate(std::distance(begin, end) * 3)),
  shapeStrides_(shape_ + std::distance(begin, end)),
  strides_(shapeStrides_ + std::distance(begin, end)),
  dimension_(std::distance(begin, end)),
  size_(1),
  coordinateOrder_(internalCoordinateOrder),
  isSimple_(true)
{
    if(dimension() != 0) {
        for(size_t j=0; j<dimension(); ++j, ++begin, ++it) {
            const size_t s = static_cast<size_t>(*begin);
            shape(j) = s;
            size() *= s;
            strides(j) = *it;
        }
        stridesFromShape(shapeBegin(), shapeEnd(), shapeStridesBegin(),
            internalCoordinateOrder);
        updateSimplicity();
    }
}

template<class A>
inline
Geometry<A>::~Geometry()
{
    allocator_.deallocate(shape_, dimension_*3);
}

template<class A>
inline Geometry<A>&
Geometry<A>::operator=
(
    const Geometry<A>& g
)
{
    if(&g != this) { // no self-assignment
        if(g.dimension_ != dimension_) {
            allocator_.deallocate(shape_, dimension_*3);
            dimension_ = g.dimension_;
            shape_ = allocator_.allocate(dimension_*3);
            shapeStrides_ = shape_+dimension_;
            strides_ = shapeStrides_+dimension_;
            dimension_ = g.dimension_;
        }
        /*
        for(size_t j=0; j<dimension_; ++j) {
            shape_[j] = g.shape_[j];
            shapeStrides_[j] = g.shapeStrides_[j];
            strides_[j] = g.strides_[j];
        }
        */
        memcpy(shape_, g.shape_, (dimension_*3)*sizeof(size_t));
        size_ = g.size_;
        coordinateOrder_ = g.coordinateOrder_;
        isSimple_ = g.isSimple_;
    }
    return *this;
}

template<class A>
inline void
Geometry<A>::resize
(
    const size_t dimension
)
{
    if(dimension != dimension_) {
        size_t* newShape = allocator_.allocate(dimension*3);
        size_t* newShapeStrides = newShape + dimension;
        size_t* newStrides = newShapeStrides + dimension;
        for(size_t j=0; j<( (dimension < dimension_) ? dimension : dimension_); ++j) {
            // save existing entries
            newShape[j] = shape(j);
            newShapeStrides[j] = shapeStrides(j);
            newStrides[j] = strides(j);
        }
        allocator_.deallocate(shape_, dimension_*3);
        shape_ = newShape;
        shapeStrides_ = newShapeStrides;
        strides_ = newStrides;
        dimension_ = dimension;
    }
}

template<class A>
inline const size_t
Geometry<A>::dimension() const
{
    return dimension_;
}

template<class A>
inline const size_t
Geometry<A>::shape(const size_t j) const
{
    Assert(MARRAY_NO_DEBUG || j<dimension_);
    return shape_[j];
}

template<class A>
inline size_t&
Geometry<A>::shape(const size_t j)
{
    Assert(MARRAY_NO_DEBUG || j<dimension_);
    return shape_[j];
}

template<class A>
inline const size_t
Geometry<A>::shapeStrides
(
    const size_t j
) const
{
    Assert(MARRAY_NO_DEBUG || j<dimension_);
    return shapeStrides_[j];
}

template<class A>
inline size_t&
Geometry<A>::shapeStrides
(
    const size_t j
)
{
    Assert(MARRAY_NO_DEBUG || j<dimension_);
    return shapeStrides_[j];
}

template<class A>
inline const size_t
Geometry<A>::strides
(
    const size_t j
) const
{
    Assert(MARRAY_NO_DEBUG || j<dimension_);
    return strides_[j];
}

template<class A>
inline size_t&
Geometry<A>::strides
(
    const size_t j
)
{
    Assert(MARRAY_NO_DEBUG || j<dimension_);
    return strides_[j];
}

template<class A>
inline const size_t*
Geometry<A>::shapeBegin() const
{
    return shape_;
}

template<class A>
inline size_t*
Geometry<A>::shapeBegin()
{
    return shape_;
}

template<class A>
inline const size_t*
Geometry<A>::shapeEnd() const
{
    return shape_ + dimension_;
}

template<class A>
inline size_t*
Geometry<A>::shapeEnd()
{
    return shape_ + dimension_;
}

template<class A>
inline const size_t*
Geometry<A>::shapeStridesBegin() const
{
    return shapeStrides_;
}

template<class A>
inline size_t*
Geometry<A>::shapeStridesBegin()
{
    return shapeStrides_;
}

template<class A>
inline const size_t*
Geometry<A>::shapeStridesEnd() const
{
    return shapeStrides_ + dimension_;
}

template<class A>
inline size_t*
Geometry<A>::shapeStridesEnd()
{
    return shapeStrides_ + dimension_;
}

template<class A>
inline const size_t*
Geometry<A>::stridesBegin() const
{
    return strides_;
}

template<class A>
inline size_t*
Geometry<A>::stridesBegin()
{
    return strides_;
}

template<class A>
inline const size_t*
Geometry<A>::stridesEnd() const
{
    return strides_ + dimension_;
}

template<class A>
inline size_t*
Geometry<A>::stridesEnd()
{
    return strides_ + dimension_;
}

template<class A>
inline const size_t
Geometry<A>::size() const
{
    return size_;
}

template<class A>
inline size_t&
Geometry<A>::size()
{
    return size_;
}

template<class A>
inline const CoordinateOrder&
Geometry<A>::coordinateOrder() const
{
    return coordinateOrder_;
}

template<class A>
inline CoordinateOrder&
Geometry<A>::coordinateOrder()
{
    return coordinateOrder_;
}

template<class A>
inline const bool
Geometry<A>::isSimple() const
{
    return isSimple_;
}

template<class A>
inline bool&
Geometry<A>::isSimple()
{
    return isSimple_;
}

template<class A>
inline void
Geometry<A>::updateSimplicity()
{
    for(size_t j=0; j<dimension(); ++j) {
        if(shapeStrides(j) != strides(j)) {
            isSimple_ = false;
            return;
        }
    }
    isSimple_ = true;
    // a 0-dimensional geometry is simple
}

template<class ShapeIterator, class StridesIterator>
inline void
stridesFromShape
(
    ShapeIterator begin,
    ShapeIterator end,
    StridesIterator strideBegin,
    const CoordinateOrder& coordinateOrder
)
{
    Assert(MARRAY_NO_ARG_TEST || std::distance(begin, end) != 0);
    size_t dimension = std::distance(begin, end);
    ShapeIterator shapeIt;
    StridesIterator strideIt;
    if(coordinateOrder == FirstMajorOrder) {
        shapeIt = begin + (dimension-1);
        strideIt = strideBegin + (dimension-1);
        *strideIt = 1;
        for(size_t j=1; j<dimension; ++j) {
            size_t tmp = *strideIt;
            --strideIt;
            (*strideIt) = tmp * (*shapeIt);
            --shapeIt;
        }
    }
    else {
        shapeIt = begin;
        strideIt = strideBegin;
        *strideIt = 1;
        for(size_t j=1; j<dimension; ++j) {
            size_t tmp = *strideIt;
            ++strideIt;
            (*strideIt) = tmp * (*shapeIt);
            ++shapeIt;
        }
    }
}

template<unsigned short N, class Functor, class T, class A>
struct OperateHelperUnary
{
    static inline void operate
    (
        View<T, false, A>& v,
        Functor f,
        T* data
    )
    {
        for(size_t j=0; j<v.shape(N-1); ++j) {
            OperateHelperUnary<N-1, Functor, T, A>::operate(v, f, data);
            data += v.strides(N-1);
        }
        data -= v.shape(N-1) * v.strides(N-1);
    }
};


template<class Functor, class T, class A>
struct OperateHelperUnary<0, Functor, T, A>
{
    static inline void operate
    (
        View<T, false, A>& v,
        Functor f,
        T* data
    )
    {
        f(*data);
    }
};

template<unsigned short N, class Functor, class T1, class T2, class A>
struct OperateHelperBinaryScalar
{
    static inline void operate
    (
        View<T1, false, A>& v,
        const T2& x,
        Functor f,
        T1* data
    )
    {
        for(size_t j=0; j<v.shape(N-1); ++j) {
            OperateHelperBinaryScalar<N-1, Functor, T1, T2, A>::operate(
                v, x, f, data);
            data += v.strides(N-1);
        }
        data -= v.shape(N-1) * v.strides(N-1);
    }
};

template<class Functor, class T1, class T2, class A>
struct OperateHelperBinaryScalar<0, Functor, T1, T2, A>
{
    static inline void operate
    (
        View<T1, false, A>& v,
        const T2& x,
        Functor f,
        T1* data
    )
    {
        f(*data, x);
    }
};

template<unsigned short N, class Functor, class T1, class T2,
         bool isConst, class A1, class A2>
struct OperateHelperBinary
{
    static inline void operate
    (
        View<T1, false, A1>& v,
        const View<T2, isConst, A2>& w,
        Functor f,
        T1* data1,
        const T2* data2
    )
    {
        for(size_t j=0; j<v.shape(N-1); ++j) {
            OperateHelperBinary<N-1, Functor, T1, T2, isConst, A1, A2>::operate(
                v, w, f, data1, data2);
            data1 += v.strides(N-1);
            data2 += w.strides(N-1);
        }
        data1 -= v.shape(N-1) * v.strides(N-1);
        data2 -= w.shape(N-1) * w.strides(N-1);
    }
};

template<class Functor, class T1, class T2, bool isConst, class A1, class A2>
struct OperateHelperBinary<0, Functor, T1, T2, isConst, A1, A2>
{
    static inline void operate
    (
        View<T1, false, A1>& v,
        const View<T2, isConst, A2>& w,
        Functor f,
        T1* data1,
        const T2* data2
    )
    {
        f(*data1, *data2);
    }
};

template<class TFrom, class TTo, class AFrom, class ATo>
struct AssignmentOperatorHelper<false, TFrom, TTo, AFrom, ATo>
{
    // from constant to mutable
    //
    // here, 'to' must be initialized (which is asserted) because
    // otherwise, the pointer to.data_ to mutable data would have
    // to be initialized with the pointer from.data_ to constant
    // data which we don't do.
    //
    static void execute
    (
        const View<TFrom, true, AFrom>& from,
        View<TTo, false, ATo>& to
    )
    {
        typedef typename View<TFrom, true, AFrom>::const_iterator FromIterator;
        typedef typename View<TTo, false, ATo>::iterator ToIterator;
        if(!MARRAY_NO_ARG_TEST) {
            Assert(from.data_ != 0 && from.dimension() == to.dimension());
            for(size_t j=0; j<from.dimension(); ++j) {
                Assert(from.shape(j) == to.shape(j));
            }
        }
        if(from.overlaps(to)) {
            Marray<TFrom, AFrom> m = from; // temporary copy
            execute(m, to);
        }
        else if(from.coordinateOrder() == to.coordinateOrder()
                && from.isSimple() && to.isSimple()
                && IsEqual<TFrom, TTo>::type) {
            memcpy(to.data_, from.data_, (from.size())*sizeof(TFrom));
        }
        else if(from.dimension() == 1)
            OperateHelperBinary<1, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
        else if(from.dimension() == 2)
            OperateHelperBinary<2, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
        else if(from.dimension() == 3)
            OperateHelperBinary<3, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
        else if(from.dimension() == 4)
            OperateHelperBinary<4, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
        else if(from.dimension() == 5)
            OperateHelperBinary<5, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
        else if(from.dimension() == 6)
            OperateHelperBinary<6, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
        else if(from.dimension() == 7)
            OperateHelperBinary<7, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
        else if(from.dimension() == 8)
            OperateHelperBinary<8, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
        else if(from.dimension() == 9)
            OperateHelperBinary<9, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
        else if(from.dimension() == 10)
            OperateHelperBinary<10, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
        else {
            FromIterator itFrom = from.begin();
            ToIterator itTo = to.begin();
            for(; itFrom.hasMore(); ++itFrom, ++itTo) {
                *itTo = static_cast<TTo>(*itFrom);
            }
        }
    }

    static void execute
    (
        const View<TFrom, false, AFrom>& from,
        View<TTo, false, ATo>& to
    )
    {
        typedef typename View<TFrom, false, AFrom>::const_iterator FromIterator;
        typedef typename View<TTo, false, ATo>::iterator ToIterator;
        if(static_cast<const void*>(&from) != static_cast<const void*>(&to)) { // no self-assignment
            if(to.data_ == 0) { // if the view 'to' is not initialized
                // initialize the view 'to' with source data
                Assert(MARRAY_NO_ARG_TEST || sizeof(TTo) == sizeof(TFrom));
                to.data_ = static_cast<TTo*>(static_cast<void*>(from.data_)); // copy pointer
                to.geometry_ = from.geometry_;
            }
            else { // if the view 'to' is initialized
                if(!MARRAY_NO_ARG_TEST) {
                    Assert(from.data_ != 0 && from.dimension() == to.dimension());
                    for(size_t j=0; j<from.dimension(); ++j) {
                        Assert(from.shape(j) == to.shape(j));
                    }
                }
                if(from.overlaps(to)) {
                    Marray<TFrom, AFrom> m = from; // temporary copy
                    execute(m, to);
                }
                else if(from.coordinateOrder() == to.coordinateOrder()
                        && from.isSimple() && to.isSimple()
                        && IsEqual<TFrom, TTo>::type) {
                    memcpy(to.data_, from.data_, (from.size())*sizeof(TFrom));
                }
                else if(from.dimension() == 1)
                    OperateHelperBinary<1, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
                else if(from.dimension() == 2)
                    OperateHelperBinary<2, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
                else if(from.dimension() == 3)
                    OperateHelperBinary<3, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
                else if(from.dimension() == 4)
                    OperateHelperBinary<4, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
                else if(from.dimension() == 5)
                    OperateHelperBinary<5, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
                else if(from.dimension() == 6)
                    OperateHelperBinary<6, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
                else if(from.dimension() == 7)
                    OperateHelperBinary<7, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
                else if(from.dimension() == 8)
                    OperateHelperBinary<8, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
                else if(from.dimension() == 9)
                    OperateHelperBinary<9, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
                else if(from.dimension() == 10)
                    OperateHelperBinary<10, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
                else {
                    FromIterator itFrom = from.begin();
                    ToIterator itTo = to.begin();
                    for(; itFrom.hasMore(); ++itFrom, ++itTo) {
                        *itTo = static_cast<TTo>(*itFrom);
                    }
                }
            }
        }
    }
};

template<class TFrom, class TTo, class AFrom, class ATo>
struct AssignmentOperatorHelper<true, TFrom, TTo, AFrom, ATo>
{
    template<bool isConstFrom>
    static void execute
    (
        const View<TFrom, isConstFrom, AFrom>& from,
        View<TTo, true, ATo>& to
    )
    {
        Assert(MARRAY_NO_ARG_TEST || sizeof(TFrom) == sizeof(TTo));
        to.data_ = static_cast<const TTo*>(
            static_cast<const void*>(from.data_)); // copy pointer
        to.geometry_ = from.geometry_;
    }
};

template<>
struct AccessOperatorHelper<true>
{
    // access by scalar index
    template<class T, class U, bool isConst, class A>
    static typename View<T, isConst, A>::reference
    execute(const View<T, isConst, A>& v, const U& index)
    {
        v.testInvariant();
        Assert(MARRAY_NO_DEBUG || v.data_ != 0);
        Assert(MARRAY_NO_DEBUG || v.dimension() != 0 || index == 0);
        size_t offset;
        v.indexToOffset(index, offset);
        return v.data_[offset];
    }
};

template<>
struct AccessOperatorHelper<false>
{
    // access by iterator
    template<class T, class U, bool isConst, class A>
    static typename View<T, isConst, A>::reference
    execute(const View<T, isConst, A>& v, U it)
    {
        v.testInvariant();
        Assert(MARRAY_NO_DEBUG ||    v.data_ != 0 );
        Assert(MARRAY_NO_DEBUG || v.dimension() != 0 || *it == 0);
        size_t offset;
        v.coordinatesToOffset(it, offset);
        return v.data_[offset];
    }
};

template<class Functor, class T, class A>
inline void
operate
(
    View<T, false, A>& v,
    Functor f
)
{
    if(v.isSimple()) {
        T* data = &v(0);
        for(size_t j=0; j<v.size(); ++j) {
            f(data[j]);
        }
    }
    else if(v.dimension() == 1)
        OperateHelperUnary<1, Functor, T, A>::operate(v, f, &v(0));
    else if(v.dimension() == 2)
        OperateHelperUnary<2, Functor, T, A>::operate(v, f, &v(0));
    else if(v.dimension() == 3)
        OperateHelperUnary<3, Functor, T, A>::operate(v, f, &v(0));
    else if(v.dimension() == 4)
        OperateHelperUnary<4, Functor, T, A>::operate(v, f, &v(0));
    else if(v.dimension() == 5)
        OperateHelperUnary<5, Functor, T, A>::operate(v, f, &v(0));
    else if(v.dimension() == 6)
        OperateHelperUnary<6, Functor, T, A>::operate(v, f, &v(0));
    else if(v.dimension() == 7)
        OperateHelperUnary<7, Functor, T, A>::operate(v, f, &v(0));
    else if(v.dimension() == 8)
        OperateHelperUnary<8, Functor, T, A>::operate(v, f, &v(0));
    else if(v.dimension() == 9)
        OperateHelperUnary<9, Functor, T, A>::operate(v, f, &v(0));
    else if(v.dimension() == 10)
        OperateHelperUnary<10, Functor, T, A>::operate(v, f, &v(0));
    else {
        for(typename View<T, false, A>::iterator it = v.begin(); it.hasMore(); ++it) {
            f(*it);
        }
    }
}

template<class Functor, class T, class A>
inline void
operate
(
    View<T, false, A>& v,
    const T& x,
    Functor f
)
{
    if(v.isSimple()) {
        T* data = &v(0);
        for(size_t j=0; j<v.size(); ++j) {
            f(data[j], x);
        }
    }
    else if(v.dimension() == 1)
        OperateHelperBinaryScalar<1, Functor, T, T, A>::operate(v, x, f, &v(0));
    else if(v.dimension() == 2)
        OperateHelperBinaryScalar<2, Functor, T, T, A>::operate(v, x, f, &v(0));
    else if(v.dimension() == 3)
        OperateHelperBinaryScalar<3, Functor, T, T, A>::operate(v, x, f, &v(0));
    else if(v.dimension() == 4)
        OperateHelperBinaryScalar<4, Functor, T, T, A>::operate(v, x, f, &v(0));
    else if(v.dimension() == 5)
        OperateHelperBinaryScalar<5, Functor, T, T, A>::operate(v, x, f, &v(0));
    else if(v.dimension() == 6)
        OperateHelperBinaryScalar<6, Functor, T, T, A>::operate(v, x, f, &v(0));
    else if(v.dimension() == 7)
        OperateHelperBinaryScalar<7, Functor, T, T, A>::operate(v, x, f, &v(0));
    else if(v.dimension() == 8)
        OperateHelperBinaryScalar<8, Functor, T, T, A>::operate(v, x, f, &v(0));
    else if(v.dimension() == 9)
        OperateHelperBinaryScalar<9, Functor, T, T, A>::operate(v, x, f, &v(0));
    else if(v.dimension() == 10)
        OperateHelperBinaryScalar<10, Functor, T, T, A>::operate(v, x, f, &v(0));
    else {
        for(typename View<T, false, A>::iterator it = v.begin(); it.hasMore(); ++it) {
            f(*it, x);
        }
    }
}

template<class Functor, class T1, class T2, bool isConst, class A>
inline void
operate
(
    View<T1, false, A>& v,
    const View<T2, isConst, A>& w,
    Functor f
)
{
    if(!MARRAY_NO_ARG_TEST) {
        Assert(v.size() != 0 && w.size() != 0);
        Assert(w.dimension() == 0 || v.dimension() == w.dimension());
        if(w.dimension() != 0) {
            for(size_t j=0; j<v.dimension(); ++j) {
                Assert(v.shape(j) == w.shape(j));
            }
        }
    }
    if(w.dimension() == 0) {
        T2 x = w(0);
        if(v.isSimple()) {
            T1* dataV = &v(0);
            for(size_t j=0; j<v.size(); ++j) {
                f(dataV[j], x);
            }
        }
        else if(v.dimension() == 1)
            OperateHelperBinaryScalar<1, Functor, T1, T2, A>::operate(v, x, f, &v(0));
        else if(v.dimension() == 2)
            OperateHelperBinaryScalar<2, Functor, T1, T2, A>::operate(v, x, f, &v(0));
        else if(v.dimension() == 3)
            OperateHelperBinaryScalar<3, Functor, T1, T2, A>::operate(v, x, f, &v(0));
        else if(v.dimension() == 4)
            OperateHelperBinaryScalar<4, Functor, T1, T2, A>::operate(v, x, f, &v(0));
        else if(v.dimension() == 5)
            OperateHelperBinaryScalar<5, Functor, T1, T2, A>::operate(v, x, f, &v(0));
        else if(v.dimension() == 6)
            OperateHelperBinaryScalar<6, Functor, T1, T2, A>::operate(v, x, f, &v(0));
        else if(v.dimension() == 7)
            OperateHelperBinaryScalar<7, Functor, T1, T2, A>::operate(v, x, f, &v(0));
        else if(v.dimension() == 8)
            OperateHelperBinaryScalar<8, Functor, T1, T2, A>::operate(v, x, f, &v(0));
        else if(v.dimension() == 9)
            OperateHelperBinaryScalar<9, Functor, T1, T2, A>::operate(v, x, f, &v(0));
        else if(v.dimension() == 10)
            OperateHelperBinaryScalar<10, Functor, T1, T2, A>::operate(v, x, f, &v(0));
        else {
            for(typename View<T1, false>::iterator it = v.begin(); it.hasMore(); ++it) {
                f(*it, x);
            }
        }
    }
    else {
        if(v.overlaps(w)) {
            Marray<T2> m = w; // temporary copy
            operate(v, m, f); // recursive call
        }
        else {
            if(v.coordinateOrder() == w.coordinateOrder()
                && v.isSimple() && w.isSimple()) {
                T1* dataV = &v(0);
                const T2* dataW = &w(0);
                for(size_t j=0; j<v.size(); ++j) {
                    f(dataV[j], dataW[j]);
                }
            }
            else if(v.dimension() == 1)
                OperateHelperBinary<1, Functor, T1, T2, isConst, A, A>::operate(v, w, f, &v(0), &w(0));
            else if(v.dimension() == 2)
                OperateHelperBinary<2, Functor, T1, T2, isConst, A, A>::operate(v, w, f, &v(0), &w(0));
            else if(v.dimension() == 3)
                OperateHelperBinary<3, Functor, T1, T2, isConst, A, A>::operate(v, w, f, &v(0), &w(0));
            else if(v.dimension() == 4)
                OperateHelperBinary<4, Functor, T1, T2, isConst, A, A>::operate(v, w, f, &v(0), &w(0));
            else if(v.dimension() == 5)
                OperateHelperBinary<5, Functor, T1, T2, isConst, A, A>::operate(v, w, f, &v(0), &w(0));
            else if(v.dimension() == 6)
                OperateHelperBinary<6, Functor, T1, T2, isConst, A, A>::operate(v, w, f, &v(0), &w(0));
            else if(v.dimension() == 7)
                OperateHelperBinary<7, Functor, T1, T2, isConst, A, A>::operate(v, w, f, &v(0), &w(0));
            else if(v.dimension() == 8)
                OperateHelperBinary<8, Functor, T1, T2, isConst, A, A>::operate(v, w, f, &v(0), &w(0));
            else if(v.dimension() == 9)
                OperateHelperBinary<9, Functor, T1, T2, isConst, A, A>::operate(v, w, f, &v(0), &w(0));
            else if(v.dimension() == 10)
                OperateHelperBinary<10, Functor, T1, T2, isConst, A, A>::operate(v, w, f, &v(0), &w(0));
            else {
                typename View<T1, false>::iterator itV = v.begin();
                typename View<T2, isConst>::const_iterator itW = w.begin();
                for(; itV.hasMore(); ++itV, ++itW) {
                    Assert(MARRAY_NO_DEBUG || itW.hasMore());
                    f(*itV, *itW);
                }
                Assert(MARRAY_NO_DEBUG || !itW.hasMore());
            }
        }
    }
}

template<class Functor, class T1, class A, class E, class T2>
inline void operate
(
    View<T1, false, A>& v,
    const ViewExpression<E, T2>& expression,
    Functor f
)
{
    const E& e = expression; // cast
    if(!MARRAY_NO_DEBUG) {
        Assert(v.size() != 0 && e.size() != 0);
        Assert(e.dimension() == v.dimension());
        if(v.dimension() == 0) {
            Assert(v.size() == 1 && e.size() == 1);
        }
        else {
            for(size_t j=0; j<v.dimension(); ++j) {
                Assert(v.shape(j) == e.shape(j));
            }
        }
    }
    if(e.overlaps(v)) {
        Marray<T1, A> m(e); // temporary copy
        operate(v, m, f);
    }
    else if(v.dimension() == 0) {
        f(v[0], e[0]);
    }
    else if(v.isSimple() && e.isSimple()
    && v.coordinateOrder() == e.coordinateOrder()) {
        for(size_t j=0; j<v.size(); ++j) {
            f(v[j], e[j]);
        }
    }
    else {
        // loop unrolling does not improve performance here
        typename E::ExpressionIterator itE(e);
        size_t offsetV = 0;
        std::vector<size_t> coordinate(v.dimension());
        size_t maxDimension = v.dimension() - 1;
        for(;;) {
            f(v[offsetV], *itE);
            for(size_t j=0; j<v.dimension(); ++j) {
                if(coordinate[j]+1 == v.shape(j)) {
                    if(j == maxDimension) {
                        return;
                    }
                    else {
                        offsetV -= coordinate[j] * v.strides(j);
                        itE.resetCoordinate(j);
                        coordinate[j] = 0;
                    }
                }
                else {
                    offsetV += v.strides(j);
                    itE.incrementCoordinate(j);
                    ++coordinate[j];
                    break;
                }
            }
        }
    }
}

} // namespace marray_detail
// \endcond suppress doxygen

} // namespace marray

#endif // #ifndef MARRAY_HXX