From f339778a59a3dcf786c96ea13c4041ee7a0d9027 Mon Sep 17 00:00:00 2001
From: Christian Zimmermann <christian.zimmermann@ur.de>
Date: Tue, 12 Feb 2019 20:26:38 +0100
Subject: [PATCH 1/9] generalize single range (static size parameter, -> -1 if
 non-static size)

---
 CMakeLists.txt                                |   2 +-
 src/include/ranges/anonymous_range.h          |  26 +-
 .../ranges/range_types/classic_range.h        |  32 +-
 src/include/ranges/range_types/null_range.h   |  34 +--
 src/include/ranges/range_types/space_range.h  |  32 +-
 src/include/ranges/range_types/spin_range.h   |  32 +-
 src/include/ranges/rbase_def.h                |  12 +-
 src/include/ranges/single_range.h             | 273 ++++++++++--------
 src/lib/ranges/anonymous_range.cc             |   2 +-
 src/lib/ranges/range_types/classic_range.cc   |  36 +--
 src/lib/ranges/range_types/null_range.cc      |  42 +--
 src/lib/ranges/range_types/space_range.cc     |  36 +--
 src/lib/ranges/range_types/spin_range.cc      |  38 +--
 src/tests/op_perf_test.cc                     |  18 +-
 14 files changed, 320 insertions(+), 295 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index ba65885..a66ce5c 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 2.8)
 
 project(multi_array)
 
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Werror -std=c++11 -Wpedantic -O3 -g -march=native -fopenmp")
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Werror -std=c++11 -Wpedantic -O3 -g -march=native -fopenmp -funroll-loops -ftree-vectorize")
 
 enable_testing()
 
diff --git a/src/include/ranges/anonymous_range.h b/src/include/ranges/anonymous_range.h
index 4cfb980..0099258 100644
--- a/src/include/ranges/anonymous_range.h
+++ b/src/include/ranges/anonymous_range.h
@@ -14,12 +14,12 @@
 namespace MultiArrayTools
 {
 
-    typedef SingleIndex<size_t,SpaceType::ANON> AnonymousIndex;
+    typedef GenSingleIndex<size_t,SpaceType::ANON,-1> AnonymousIndex;
 
     //template <class R>
     //using SIZET = size_t;
 
-    typedef SingleRange<size_t,SpaceType::ANON> AnonymousRange;
+    typedef GenSingleRange<size_t,SpaceType::ANON,-1> AnonymousRange;
     
     // NOT THREAD SAVE!!
     class AnonymousRangeFactory : public RangeFactoryBase
@@ -53,7 +53,7 @@ namespace MultiArrayTools
     };
 
     template <>
-    class SingleRange<size_t,SpaceType::ANON> : public RangeInterface<AnonymousIndex>
+    class GenSingleRange<size_t,SpaceType::ANON,-1> : public RangeInterface<AnonymousIndex>
     {
     public:
 
@@ -64,7 +64,7 @@ namespace MultiArrayTools
 	
 	typedef RangeBase RB;
 	typedef typename RangeInterface<AnonymousIndex>::IndexType IndexType;
-	typedef SingleRange<size_t,SpaceType::ANON> RangeType;
+	typedef GenSingleRange<size_t,SpaceType::ANON,-1> RangeType;
 	typedef size_t MetaType;
 	
 	virtual size_t size() const final;
@@ -105,16 +105,16 @@ namespace MultiArrayTools
 	
     protected:
 
-	SingleRange() = default;
-	SingleRange(const AnonymousRange& in) = default;
+	GenSingleRange() = default;
+	GenSingleRange(const AnonymousRange& in) = default;
 
 	template <class... RangeTypes>
-	SingleRange(const std::tuple<std::shared_ptr<RangeTypes>...>& origs);
+	GenSingleRange(const std::tuple<std::shared_ptr<RangeTypes>...>& origs);
 	
 	template <class... RangeTypes>
-	SingleRange(std::shared_ptr<RangeTypes>... origs);
+	GenSingleRange(std::shared_ptr<RangeTypes>... origs);
 
-	SingleRange(const std::vector<std::shared_ptr<RangeBase>>& origs);
+	GenSingleRange(const std::vector<std::shared_ptr<RangeBase>>& origs);
 	
 	size_t mSize = 1;
 	bool mEmpty = true;
@@ -203,7 +203,7 @@ namespace MultiArrayTools
      ***********************/
 
     template <class... RangeTypes>
-    SingleRange<size_t,SpaceType::ANON>::SingleRange(const std::tuple<std::shared_ptr<RangeTypes>...>& origs) :
+    GenSingleRange<size_t,SpaceType::ANON,-1>::GenSingleRange(const std::tuple<std::shared_ptr<RangeTypes>...>& origs) :
 	RangeInterface<AnonymousIndex>()
     {
 	RPackNum<sizeof...(RangeTypes)-1>::RangesToVec( origs, mOrig );
@@ -214,7 +214,7 @@ namespace MultiArrayTools
     }
     
     template <class... RangeTypes>
-    SingleRange<size_t,SpaceType::ANON>::SingleRange(std::shared_ptr<RangeTypes>... origs) :
+    GenSingleRange<size_t,SpaceType::ANON,-1>::GenSingleRange(std::shared_ptr<RangeTypes>... origs) :
 	RangeInterface<AnonymousIndex>()
     {
 	auto rst = std::make_tuple(origs...);
@@ -226,13 +226,13 @@ namespace MultiArrayTools
     }
    
     template <class Range>
-    std::shared_ptr<Range> SingleRange<size_t,SpaceType::ANON>::fullsub(size_t num) const
+    std::shared_ptr<Range> GenSingleRange<size_t,SpaceType::ANON,-1>::fullsub(size_t num) const
     {
 	return std::dynamic_pointer_cast<Range>( mOrig.at(num) );
     }
 
     template <class... Ranges>
-    std::shared_ptr<MultiRange<Ranges...> > SingleRange<size_t,SpaceType::ANON>::scast(SIZET<Ranges>... sizes) const
+    std::shared_ptr<MultiRange<Ranges...> > GenSingleRange<size_t,SpaceType::ANON,-1>::scast(SIZET<Ranges>... sizes) const
     {
 	std::tuple<std::shared_ptr<Ranges>...> rtp;
 	RPackNum<sizeof...(Ranges)-1>::resolveRangeType(mOrig, rtp, 0, sizes...);
diff --git a/src/include/ranges/range_types/classic_range.h b/src/include/ranges/range_types/classic_range.h
index 3a65fea..f9b4ed6 100644
--- a/src/include/ranges/range_types/classic_range.h
+++ b/src/include/ranges/range_types/classic_range.h
@@ -12,29 +12,29 @@
 
 namespace MultiArrayTools
 {
-    typedef SingleIndex<size_t,SpaceType::NONE> ClassicIndex;
+    typedef GenSingleIndex<size_t,SpaceType::NONE,-1> ClassicIndex;
 
     template <>
-    class SingleRangeFactory<size_t,SpaceType::NONE> : public RangeFactoryBase
+    class GenSingleRangeFactory<size_t,SpaceType::NONE,-1> : public RangeFactoryBase
     {
     public:
 	
-	typedef SingleRange<size_t,SpaceType::NONE> oType;
+	typedef GenSingleRange<size_t,SpaceType::NONE,-1> oType;
 
-	SingleRangeFactory(size_t size = 0);
+	GenSingleRangeFactory(size_t size = 0);
 	std::shared_ptr<RangeBase> create();
 	
     };
 
     template <>
-    class SingleRange<size_t,SpaceType::NONE> : public RangeInterface<ClassicIndex>
+    class GenSingleRange<size_t,SpaceType::NONE,-1> : public RangeInterface<ClassicIndex>
     {
     public:
 	typedef RangeBase RB;
-	typedef typename RangeInterface<SingleIndex<size_t,SpaceType::NONE> >::IndexType IndexType;
-	typedef SingleRange<size_t,SpaceType::NONE> RangeType;
+	typedef typename RangeInterface<GenSingleIndex<size_t,SpaceType::NONE,-1> >::IndexType IndexType;
+	typedef GenSingleRange<size_t,SpaceType::NONE,-1> RangeType;
 	typedef size_t MetaType;
-        typedef SingleRangeFactory<size_t,SpaceType::NONE> FType; 
+        typedef GenSingleRangeFactory<size_t,SpaceType::NONE,-1> FType; 
         
 	virtual size_t size() const final;
 	virtual size_t dim() const final;
@@ -52,28 +52,28 @@ namespace MultiArrayTools
 	virtual IndexType end() const final;
 	//virtual std::shared_ptr<VIWB> index() const final;
 	
-	friend SingleRangeFactory<size_t,SpaceType::NONE>;
+	friend GenSingleRangeFactory<size_t,SpaceType::NONE,-1>;
 
 	static constexpr bool defaultable = true;
 	static constexpr size_t ISSTATIC = 0;
 	static constexpr size_t SIZE = -1;
 	static constexpr bool HASMETACONT = false;
 	
-	static SingleRangeFactory<size_t, SpaceType::NONE> factory(size_t size = 0)
-	{ return SingleRangeFactory<size_t, SpaceType::NONE>(size); }
+	static GenSingleRangeFactory<size_t,SpaceType::NONE,-1> factory(size_t size = 0)
+	{ return GenSingleRangeFactory<size_t,SpaceType::NONE,-1>(size); }
 	
     protected:
 
 	size_t mSize = 0; 
 	
-	SingleRange() = default;
-	SingleRange(const SingleRange& in) = delete;
+	GenSingleRange() = default;
+	GenSingleRange(const GenSingleRange& in) = delete;
 
-	SingleRange(size_t size);
+	GenSingleRange(size_t size);
     };
 
-    typedef SingleRange<size_t,SpaceType::NONE> ClassicRange;
-    typedef SingleRangeFactory<size_t,SpaceType::NONE> ClassicRF;
+    typedef GenSingleRange<size_t,SpaceType::NONE,-1> ClassicRange;
+    typedef GenSingleRangeFactory<size_t,SpaceType::NONE,-1> ClassicRF;
 }
 
 
diff --git a/src/include/ranges/range_types/null_range.h b/src/include/ranges/range_types/null_range.h
index 968c0c0..994bb29 100644
--- a/src/include/ranges/range_types/null_range.h
+++ b/src/include/ranges/range_types/null_range.h
@@ -12,19 +12,19 @@ include_range_type(NUL,-2)
 
 namespace MultiArrayTools
 {   
-    typedef SingleIndex<size_t,SpaceType::NUL> NullIndex;
+    typedef GenSingleIndex<size_t,SpaceType::NUL,0> NullIndex;
 
-    std::shared_ptr<SingleRange<size_t,SpaceType::NUL> > nullr();
+    std::shared_ptr<GenSingleRange<size_t,SpaceType::NUL,0> > nullr();
     std::shared_ptr<NullIndex> nulli();
 
     template <>
-    class SingleRangeFactory<size_t,SpaceType::NUL> : public RangeFactoryBase
+    class GenSingleRangeFactory<size_t,SpaceType::NUL,0> : public RangeFactoryBase
     {
     public:
 	
-	typedef SingleRange<size_t,SpaceType::NUL> oType;
+	typedef GenSingleRange<size_t,SpaceType::NUL,0> oType;
 
-	SingleRangeFactory();
+	GenSingleRangeFactory();
 	std::shared_ptr<RangeBase> create();
 
 	friend std::shared_ptr<oType> nullr();
@@ -34,14 +34,14 @@ namespace MultiArrayTools
     };
 
     template <>
-    class SingleRange<size_t,SpaceType::NUL> : public RangeInterface<NullIndex>
+    class GenSingleRange<size_t,SpaceType::NUL,0> : public RangeInterface<NullIndex>
     {
     public:
 	typedef RangeBase RB;
-	typedef typename RangeInterface<SingleIndex<size_t,SpaceType::NUL> >::IndexType IndexType;
-	typedef SingleRange<size_t,SpaceType::NUL> RangeType;
+	typedef typename RangeInterface<GenSingleIndex<size_t,SpaceType::NUL,0> >::IndexType IndexType;
+	typedef GenSingleRange<size_t,SpaceType::NUL,0> RangeType;
 	typedef size_t MetaType;
-        typedef SingleRangeFactory<size_t,SpaceType::NUL> FType; 
+        typedef GenSingleRangeFactory<size_t,SpaceType::NUL,0> FType; 
 	
 	virtual size_t size() const final;
 	virtual size_t dim() const final;
@@ -59,7 +59,7 @@ namespace MultiArrayTools
 	virtual IndexType end() const final;
 	//virtual std::shared_ptr<VIWB> index() const final;
 	
-	friend SingleRangeFactory<size_t,SpaceType::NUL>;
+	friend GenSingleRangeFactory<size_t,SpaceType::NUL,0>;
 
 	static constexpr bool defaultable = true;
 
@@ -67,19 +67,19 @@ namespace MultiArrayTools
 	static constexpr size_t SIZE = 1;
 	static constexpr bool HASMETACONT = false;
 	
-	static SingleRangeFactory<size_t,SpaceType::NUL> factory()
-	{ return SingleRangeFactory<size_t,SpaceType::NUL>(); }
+	static GenSingleRangeFactory<size_t,SpaceType::NUL,0> factory()
+	{ return GenSingleRangeFactory<size_t,SpaceType::NUL,0>(); }
 	
     protected:
 
-	SingleRange() = default;
-	SingleRange(const SingleRange& in) = delete;
+	GenSingleRange() = default;
+	GenSingleRange(const GenSingleRange& in) = delete;
 
-	//SingleRange(size_t spinNum);
+	//GenSingleRange(size_t spinNum);
     };
 
-    typedef SingleRange<size_t,SpaceType::NUL> NullRange;
-    typedef SingleRangeFactory<size_t,SpaceType::NUL> NullRF;
+    typedef GenSingleRange<size_t,SpaceType::NUL,0> NullRange;
+    typedef GenSingleRangeFactory<size_t,SpaceType::NUL,0> NullRF;
     
     std::shared_ptr<NullRF> mkNUL(const char* dp, size_t size);
 	
diff --git a/src/include/ranges/range_types/space_range.h b/src/include/ranges/range_types/space_range.h
index 238afb5..46d660a 100644
--- a/src/include/ranges/range_types/space_range.h
+++ b/src/include/ranges/range_types/space_range.h
@@ -11,29 +11,29 @@ include_range_type(PSPACE,3) // Periodic 1dim space
 namespace MultiArrayTools
 {
      // Periodic 1dim space
-    typedef SingleIndex<int,SpaceType::PSPACE> XSpaceIndex;
+    typedef GenSingleIndex<int,SpaceType::PSPACE,-1> XSpaceIndex;
 
     template <>
-    class SingleRangeFactory<int,SpaceType::PSPACE> : public RangeFactoryBase
+    class GenSingleRangeFactory<int,SpaceType::PSPACE,-1> : public RangeFactoryBase
     {
     public:
 	
-	typedef SingleRange<int,SpaceType::PSPACE> oType;
+	typedef GenSingleRange<int,SpaceType::PSPACE,-1> oType;
 
-	SingleRangeFactory(size_t size = 0);
+	GenSingleRangeFactory(size_t size = 0);
 	std::shared_ptr<RangeBase> create();
 	
     };
     
     template <>
-    class SingleRange<int,SpaceType::PSPACE> : public RangeInterface<XSpaceIndex>
+    class GenSingleRange<int,SpaceType::PSPACE,-1> : public RangeInterface<XSpaceIndex>
     {
     public:
 	typedef RangeBase RB;
-	typedef typename RangeInterface<SingleIndex<int,SpaceType::PSPACE> >::IndexType IndexType;
-	typedef SingleRange<int,SpaceType::PSPACE> RangeType;
+	typedef typename RangeInterface<GenSingleIndex<int,SpaceType::PSPACE,-1> >::IndexType IndexType;
+	typedef GenSingleRange<int,SpaceType::PSPACE,-1> RangeType;
 	typedef int MetaType;
-        typedef SingleRangeFactory<int,SpaceType::PSPACE> FType; 
+        typedef GenSingleRangeFactory<int,SpaceType::PSPACE,-1> FType; 
 	
 	virtual size_t size() const final;
 	virtual size_t dim() const final;
@@ -50,28 +50,28 @@ namespace MultiArrayTools
 	virtual IndexType begin() const final;
 	virtual IndexType end() const final;
 	
-	friend SingleRangeFactory<int,SpaceType::PSPACE>;
+	friend GenSingleRangeFactory<int,SpaceType::PSPACE,-1>;
 
 	static constexpr bool defaultable = true;
 	static constexpr size_t ISSTATIC = 0;
 	static constexpr size_t SIZE = -1;
 	static constexpr bool HASMETACONT = false;
 	
-	static SingleRangeFactory<int, SpaceType::PSPACE> factory(size_t size = 0)
-	{ return SingleRangeFactory<int, SpaceType::PSPACE>(size); }
+	static GenSingleRangeFactory<int,SpaceType::PSPACE,-1> factory(size_t size = 0)
+	{ return GenSingleRangeFactory<int,SpaceType::PSPACE,-1>(size); }
 	
     protected:
 
 	size_t mSize = 0;
 	
-	SingleRange() = default;
-	SingleRange(const SingleRange& in) = delete;
+	GenSingleRange() = default;
+	GenSingleRange(const GenSingleRange& in) = delete;
 	
-	SingleRange(size_t size);
+	GenSingleRange(size_t size);
     };
 
-    typedef SingleRange<int,SpaceType::PSPACE> PSpaceRange;
-    typedef SingleRangeFactory<int,SpaceType::PSPACE> PSpaceRF;
+    typedef GenSingleRange<int,SpaceType::PSPACE,-1> PSpaceRange;
+    typedef GenSingleRangeFactory<int,SpaceType::PSPACE,-1> PSpaceRF;
 
     std::shared_ptr<PSpaceRF> mkPSPACE(const char* dp, size_t size);
     
diff --git a/src/include/ranges/range_types/spin_range.h b/src/include/ranges/range_types/spin_range.h
index 32c8dff..97992b5 100644
--- a/src/include/ranges/range_types/spin_range.h
+++ b/src/include/ranges/range_types/spin_range.h
@@ -12,29 +12,29 @@ include_range_type(SPIN,2)
 
 namespace MultiArrayTools
 {
-    typedef SingleIndex<size_t,SpaceType::SPIN> SpinIndex;
+    typedef GenSingleIndex<size_t,SpaceType::SPIN,4> SpinIndex;
 
     template <>
-    class SingleRangeFactory<size_t,SpaceType::SPIN> : public RangeFactoryBase
+    class GenSingleRangeFactory<size_t,SpaceType::SPIN,4> : public RangeFactoryBase
     {
     public:
 	
-	typedef SingleRange<size_t,SpaceType::SPIN> oType;
+	typedef GenSingleRange<size_t,SpaceType::SPIN,4> oType;
 
-	SingleRangeFactory();
+	GenSingleRangeFactory();
 	std::shared_ptr<RangeBase> create();
 	
     };
 
     template <>
-    class SingleRange<size_t,SpaceType::SPIN> : public RangeInterface<SpinIndex>
+    class GenSingleRange<size_t,SpaceType::SPIN,4> : public RangeInterface<SpinIndex>
     {
     public:
 	typedef RangeBase RB;
-	typedef typename RangeInterface<SingleIndex<size_t,SpaceType::SPIN> >::IndexType IndexType;
-	typedef SingleRange<size_t,SpaceType::SPIN> RangeType;
+	typedef typename RangeInterface<GenSingleIndex<size_t,SpaceType::SPIN,4> >::IndexType IndexType;
+	typedef GenSingleRange<size_t,SpaceType::SPIN,4> RangeType;
 	typedef size_t MetaType;
-        typedef SingleRangeFactory<size_t,SpaceType::SPIN> FType; 
+        typedef GenSingleRangeFactory<size_t,SpaceType::SPIN,4> FType; 
         
 	virtual size_t size() const final;
 	virtual size_t dim() const final;
@@ -52,7 +52,7 @@ namespace MultiArrayTools
 	virtual IndexType end() const final;
 	//virtual std::shared_ptr<VIWB> index() const final;
 	
-	friend SingleRangeFactory<size_t,SpaceType::SPIN>;
+	friend GenSingleRangeFactory<size_t,SpaceType::SPIN,4>;
 
 	static constexpr bool defaultable = true;
 	static constexpr size_t mSpinNum = 4;
@@ -61,19 +61,19 @@ namespace MultiArrayTools
 	static constexpr size_t SIZE = mSpinNum;
 	static constexpr bool HASMETACONT = false;
 	
-	static SingleRangeFactory<size_t, SpaceType::SPIN> factory()
-	{ return SingleRangeFactory<size_t, SpaceType::SPIN>(); }
+	static GenSingleRangeFactory<size_t,SpaceType::SPIN,4> factory()
+	{ return GenSingleRangeFactory<size_t,SpaceType::SPIN,4>(); }
 	
     protected:
 
-	SingleRange() = default;
-	SingleRange(const SingleRange& in) = delete;
+	GenSingleRange() = default;
+	GenSingleRange(const GenSingleRange& in) = delete;
 
-	//SingleRange(size_t spinNum);
+	//GenSingleRange(size_t spinNum);
     };
 
-    typedef SingleRange<size_t,SpaceType::SPIN> SpinRange;
-    typedef SingleRangeFactory<size_t,SpaceType::SPIN> SpinRF;
+    typedef GenSingleRange<size_t,SpaceType::SPIN,4> SpinRange;
+    typedef GenSingleRangeFactory<size_t,SpaceType::SPIN,4> SpinRF;
 
     std::shared_ptr<SpinRF> mkSPIN(const char* dp, size_t size);
 }
diff --git a/src/include/ranges/rbase_def.h b/src/include/ranges/rbase_def.h
index 3cffe08..64ccc45 100644
--- a/src/include/ranges/rbase_def.h
+++ b/src/include/ranges/rbase_def.h
@@ -44,16 +44,16 @@ namespace MultiArrayTools
     class IndexInterface;
 
     // single_range.h    
-    template <typename U, SpaceType TYPE>
-    class SingleRange;
+    template <typename U, SpaceType TYPE, size_t S>
+    class GenSingleRange;
 
     // single_range.h    
-    template <typename U, SpaceType TYPE>
-    class SingleRangeFactory;
+    template <typename U, SpaceType TYPE, size_t S>
+    class GenSingleRangeFactory;
 
     // single_range.h
-    template <typename U, SpaceType TYPE>
-    class SingleIndex;
+    template <typename U, SpaceType TYPE, size_t S>
+    class GenSingleIndex;
 
     // subrange.h
     template <class Index>
diff --git a/src/include/ranges/single_range.h b/src/include/ranges/single_range.h
index 1affcaf..57d76f7 100644
--- a/src/include/ranges/single_range.h
+++ b/src/include/ranges/single_range.h
@@ -28,19 +28,19 @@ namespace MultiArrayTools
     }
 
     
-    template <typename U, SpaceType TYPE>
-    class SingleIndex : public IndexInterface<SingleIndex<U,TYPE>,U>
+    template <typename U, SpaceType TYPE, size_t S>
+    class GenSingleIndex : public IndexInterface<GenSingleIndex<U,TYPE,S>,U>
     {
     public:
 
-	typedef IndexInterface<SingleIndex<U,TYPE>,U> IB;
+	typedef IndexInterface<GenSingleIndex<U,TYPE,S>,U> IB;
 	typedef U MetaType;
-	typedef SingleRange<U,TYPE> RangeType;
-	typedef SingleIndex IType;
+	typedef GenSingleRange<U,TYPE,S> RangeType;
+	typedef GenSingleIndex IType;
 
-	//DEFAULT_MEMBERS_X(SingleIndex);
+	//DEFAULT_MEMBERS_X(GenSingleIndex);
 	
-	SingleIndex(const std::shared_ptr<SingleRange<U,TYPE> >& range);
+	GenSingleIndex(const std::shared_ptr<GenSingleRange<U,TYPE,S> >& range);
 
 	static constexpr IndexType sType() { return IndexType::SINGLE; }
 	static constexpr size_t totalDim() { return 1; }
@@ -53,9 +53,9 @@ namespace MultiArrayTools
 	
 	IndexType type() const;
 	
-	SingleIndex& operator=(size_t pos);
-	SingleIndex& operator++();
-	SingleIndex& operator--();
+	GenSingleIndex& operator=(size_t pos);
+	GenSingleIndex& operator++();
+	GenSingleIndex& operator--();
 
 	int pp(std::intptr_t idxPtrNum);
 	int mm(std::intptr_t idxPtrNum);
@@ -63,7 +63,7 @@ namespace MultiArrayTools
 	std::string stringMeta() const;
 	U meta() const;
 	const U* metaPtr() const;
-	SingleIndex& at(const U& metaPos);
+	GenSingleIndex& at(const U& metaPos);
 	size_t posAt(const U& metaPos) const;
 
 	bool isMeta(const U& metaPos) const;
@@ -84,30 +84,30 @@ namespace MultiArrayTools
 
 	template <class Expr>
 	auto ifor(size_t step, Expr ex) const
-	    -> For<SingleIndex<U,TYPE>,Expr>;
+	    -> For<GenSingleIndex<U,TYPE,S>,Expr>;
 
 	template <class Expr>
 	auto iforh(size_t step, Expr ex) const
-	    -> For<SingleIndex<U,TYPE>,Expr,ForType::HIDDEN>;
+	    -> For<GenSingleIndex<U,TYPE,S>,Expr,ForType::HIDDEN>;
 
         template <class Expr>
 	auto pifor(size_t step, Expr ex) const
-	    -> PFor<SingleIndex<U,TYPE>,Expr>;
+	    -> PFor<GenSingleIndex<U,TYPE,S>,Expr>;
         
     private:
 	std::shared_ptr<RangeType> mExplicitRangePtr;
 	const U* mMetaPtr;
     };
 
-    template <typename U, SpaceType TYPE>
-    class SingleRangeFactory : public RangeFactoryBase
+    template <typename U, SpaceType TYPE, size_t S>
+    class GenSingleRangeFactory : public RangeFactoryBase
     {
     public:
 	
-	typedef SingleRange<U,TYPE> oType;
+	typedef GenSingleRange<U,TYPE,S> oType;
 
-	SingleRangeFactory() = delete;
-	SingleRangeFactory(const std::vector<U>& space);
+	GenSingleRangeFactory() = delete;
+	GenSingleRangeFactory(const std::vector<U>& space);
 	std::shared_ptr<RangeBase> create();
 	
     };
@@ -186,17 +186,31 @@ namespace MultiArrayTools
 	}
 	
     };
-    
-    template <typename U, SpaceType TYPE>
-    class SingleRange : public RangeInterface<SingleIndex<U,TYPE> >
+
+    template <size_t S>
+    struct CheckStatic
+    {
+        static constexpr size_t ISSTATIC = true;
+        static constexpr size_t SIZE = S;
+    };
+
+    template <>
+    struct CheckStatic<-1>
+    {
+        static constexpr size_t ISSTATIC = false;
+        static constexpr size_t SIZE = -1;
+    };
+
+    template <typename U, SpaceType TYPE, size_t S>
+    class GenSingleRange : public RangeInterface<GenSingleIndex<U,TYPE,S> >
     {
     public:
 	typedef RangeBase RB;
-	typedef SingleIndex<U,TYPE> IndexType;
-	typedef SingleRange RangeType;
+	typedef GenSingleIndex<U,TYPE,S> IndexType;
+	typedef GenSingleRange RangeType;
 	typedef U MetaType;
-	typedef SingleRangeFactory<U,TYPE> FType; 
-	//typedef typename RangeInterface<SingleIndex<U,TYPE> >::IndexType IndexType;
+	typedef GenSingleRangeFactory<U,TYPE,S> FType; 
+	//typedef typename RangeInterface<GenSingleIndex<U,TYPE,S> >::IndexType IndexType;
 	
 	virtual size_t size() const final;
 	virtual size_t dim() const final;
@@ -215,25 +229,34 @@ namespace MultiArrayTools
 	virtual IndexType begin() const final;
 	virtual IndexType end() const final;
 		
-	friend SingleRangeFactory<U,TYPE>;
+	friend GenSingleRangeFactory<U,TYPE,S>;
 
 	static constexpr bool defaultable = false;
-	static constexpr size_t ISSTATIC = 0;
-	static constexpr size_t SIZE = -1;
+	static constexpr size_t ISSTATIC = CheckStatic<S>::ISSTATIC;
+	static constexpr size_t SIZE = CheckStatic<S>::SIZE;
 	static constexpr bool HASMETACONT = true;
 	
     protected:
 
-	SingleRange() = delete;
-	SingleRange(const SingleRange& in) = delete;
+	GenSingleRange() = delete;
+	GenSingleRange(const GenSingleRange& in) = delete;
 	
-	SingleRange(const std::vector<U>& space);
+	GenSingleRange(const std::vector<U>& space);
 
 	std::vector<U> mSpace;
 	//std::map<U,size_t> mMSpace;
 	MetaMap<U> mMSpace;
     };
 
+    template <typename U, SpaceType TYPE>
+    using SingleRange = GenSingleRange<U,TYPE,-1>;
+
+    template <typename U, SpaceType TYPE>
+    using SingleIndex = GenSingleIndex<U,TYPE,-1>;
+
+    template <typename U, SpaceType TYPE>
+    using SingleRangeFactory = GenSingleRangeFactory<U,TYPE,-1>;
+
 }
 
 /* ========================= *
@@ -243,7 +266,7 @@ namespace MultiArrayTools
 namespace MultiArrayTools
 {
     /******************
-     *  SingleIndex   *	     
+     *  GenSingleIndex   *	     
      ******************/
 
     template <bool HASMETACONT>
@@ -278,134 +301,134 @@ namespace MultiArrayTools
 	}
     };
 
-    template <typename U, SpaceType TYPE>
-    SingleIndex<U,TYPE>::SingleIndex(const std::shared_ptr<SingleRange<U,TYPE> >& range) :
-	IndexInterface<SingleIndex<U,TYPE>,U>(range, 0),
+    template <typename U, SpaceType TYPE, size_t S>
+    GenSingleIndex<U,TYPE,S>::GenSingleIndex(const std::shared_ptr<GenSingleRange<U,TYPE,S> >& range) :
+	IndexInterface<GenSingleIndex<U,TYPE,S>,U>(range, 0),
 	mExplicitRangePtr(std::dynamic_pointer_cast<RangeType>(IB::mRangePtr)),
-	mMetaPtr(MetaPtrHandle<SingleIndex<U,TYPE>::RangeType::HASMETACONT>::set
+	mMetaPtr(MetaPtrHandle<GenSingleIndex<U,TYPE,S>::RangeType::HASMETACONT>::set
 		 ( dynamic_cast<RangeType*>(IB::mRangePtr.get() ) ) ) {}
 
-    template <typename U, SpaceType TYPE>
-    IndexType SingleIndex<U,TYPE>::type() const
+    template <typename U, SpaceType TYPE, size_t S>
+    IndexType GenSingleIndex<U,TYPE,S>::type() const
     {
 	return IndexType::SINGLE;
     }
 
-    template <typename U, SpaceType TYPE>
-    SingleIndex<U,TYPE>& SingleIndex<U,TYPE>::operator=(size_t pos)
+    template <typename U, SpaceType TYPE, size_t S>
+    GenSingleIndex<U,TYPE,S>& GenSingleIndex<U,TYPE,S>::operator=(size_t pos)
     {
 	IB::mPos = pos;
 	return *this;
     }
 
-    template <typename U, SpaceType TYPE>
-    SingleIndex<U,TYPE>& SingleIndex<U,TYPE>::operator++()
+    template <typename U, SpaceType TYPE, size_t S>
+    GenSingleIndex<U,TYPE,S>& GenSingleIndex<U,TYPE,S>::operator++()
     {
 	++IB::mPos;
 	return *this;
     }
 
-    template <typename U, SpaceType TYPE>
-    SingleIndex<U,TYPE>& SingleIndex<U,TYPE>::operator--()
+    template <typename U, SpaceType TYPE, size_t S>
+    GenSingleIndex<U,TYPE,S>& GenSingleIndex<U,TYPE,S>::operator--()
     {
 	--IB::mPos;
 	return *this;
     }
 
-    template <typename U, SpaceType TYPE>
-    int SingleIndex<U,TYPE>::pp(std::intptr_t idxPtrNum)
+    template <typename U, SpaceType TYPE, size_t S>
+    int GenSingleIndex<U,TYPE,S>::pp(std::intptr_t idxPtrNum)
     {
 	++(*this);
 	return 1;
     }
 
-    template <typename U, SpaceType TYPE>
-    int SingleIndex<U,TYPE>::mm(std::intptr_t idxPtrNum)
+    template <typename U, SpaceType TYPE, size_t S>
+    int GenSingleIndex<U,TYPE,S>::mm(std::intptr_t idxPtrNum)
     {
 	--(*this);
 	return 1;
     }
 
-    template <typename U, SpaceType TYPE>
-    std::string SingleIndex<U,TYPE>::stringMeta() const
+    template <typename U, SpaceType TYPE, size_t S>
+    std::string GenSingleIndex<U,TYPE,S>::stringMeta() const
     {
-	return std::dynamic_pointer_cast<SingleRange<U,TYPE> const>( IB::mRangePtr )->stringMeta(IB::mPos);
+	return std::dynamic_pointer_cast<GenSingleRange<U,TYPE,S> const>( IB::mRangePtr )->stringMeta(IB::mPos);
     }
     
-    template <typename U, SpaceType TYPE>
-    U SingleIndex<U,TYPE>::meta() const
+    template <typename U, SpaceType TYPE, size_t S>
+    U GenSingleIndex<U,TYPE,S>::meta() const
     {
-	return MetaPtrHandle<SingleIndex<U,TYPE>::RangeType::HASMETACONT>::getMeta
+	return MetaPtrHandle<GenSingleIndex<U,TYPE,S>::RangeType::HASMETACONT>::getMeta
 	    ( mMetaPtr, IB::mPos, mExplicitRangePtr );
     }
 
-    template <typename U, SpaceType TYPE>
-    const U* SingleIndex<U,TYPE>::metaPtr() const
+    template <typename U, SpaceType TYPE, size_t S>
+    const U* GenSingleIndex<U,TYPE,S>::metaPtr() const
     {
 	return mMetaPtr;
     }
 
-    template <typename U, SpaceType TYPE>
-    bool SingleIndex<U,TYPE>::isMeta(const U& metaPos) const
+    template <typename U, SpaceType TYPE, size_t S>
+    bool GenSingleIndex<U,TYPE,S>::isMeta(const U& metaPos) const
     {
 	return mExplicitRangePtr->isMeta(metaPos);
     }
     
-    template <typename U, SpaceType TYPE>
-    SingleIndex<U,TYPE>& SingleIndex<U,TYPE>::at(const U& metaPos)
+    template <typename U, SpaceType TYPE, size_t S>
+    GenSingleIndex<U,TYPE,S>& GenSingleIndex<U,TYPE,S>::at(const U& metaPos)
     {
-	(*this) = std::dynamic_pointer_cast<SingleRange<U,TYPE> const>( IB::mRangePtr )->getMeta( metaPos );
+	(*this) = std::dynamic_pointer_cast<GenSingleRange<U,TYPE,S> const>( IB::mRangePtr )->getMeta( metaPos );
 	return *this;
     }
 
-    template <typename U, SpaceType TYPE>
-    size_t SingleIndex<U,TYPE>::posAt(const U& metaPos) const
+    template <typename U, SpaceType TYPE, size_t S>
+    size_t GenSingleIndex<U,TYPE,S>::posAt(const U& metaPos) const
     {
-	return std::dynamic_pointer_cast<SingleRange<U,TYPE> const>( IB::mRangePtr )->getMeta( metaPos );
+	return std::dynamic_pointer_cast<GenSingleRange<U,TYPE,S> const>( IB::mRangePtr )->getMeta( metaPos );
     }
 
-    template <typename U, SpaceType TYPE>
-    size_t SingleIndex<U,TYPE>::dim() // = 1
+    template <typename U, SpaceType TYPE, size_t S>
+    size_t GenSingleIndex<U,TYPE,S>::dim() // = 1
     {
 	return 1;
     }
 
-    template <typename U, SpaceType TYPE>
-    bool SingleIndex<U,TYPE>::last()
+    template <typename U, SpaceType TYPE, size_t S>
+    bool GenSingleIndex<U,TYPE,S>::last()
     {
 	return IB::mPos == IB::mMax - 1;
     }
 
-    template <typename U, SpaceType TYPE>
-    bool SingleIndex<U,TYPE>::first()
+    template <typename U, SpaceType TYPE, size_t S>
+    bool GenSingleIndex<U,TYPE,S>::first()
     {
 	return IB::mPos == 0;
     }
 
-    template <typename U, SpaceType TYPE>
-    std::shared_ptr<typename SingleIndex<U,TYPE>::RangeType> SingleIndex<U,TYPE>::range()
+    template <typename U, SpaceType TYPE, size_t S>
+    std::shared_ptr<typename GenSingleIndex<U,TYPE,S>::RangeType> GenSingleIndex<U,TYPE,S>::range()
     {
 	return mExplicitRangePtr;
     }
 
-    template <typename U, SpaceType TYPE>
+    template <typename U, SpaceType TYPE, size_t S>
     template <size_t N>
-    void SingleIndex<U,TYPE>::getPtr() {}
+    void GenSingleIndex<U,TYPE,S>::getPtr() {}
     
-    template <typename U, SpaceType TYPE>
-    size_t SingleIndex<U,TYPE>::getStepSize(size_t n)
+    template <typename U, SpaceType TYPE, size_t S>
+    size_t GenSingleIndex<U,TYPE,S>::getStepSize(size_t n)
     {
 	return 1;
     }
 
-    template <typename U, SpaceType TYPE>
-    std::string SingleIndex<U,TYPE>::id() const
+    template <typename U, SpaceType TYPE, size_t S>
+    std::string GenSingleIndex<U,TYPE,S>::id() const
     {
 	return std::string("sin") + std::to_string(IB::mId);
     }
 
-    template <typename U, SpaceType TYPE>
-    void SingleIndex<U,TYPE>::print(size_t offset)
+    template <typename U, SpaceType TYPE, size_t S>
+    void GenSingleIndex<U,TYPE,S>::print(size_t offset)
     {
 	if(offset == 0){
 	    std::cout << " === " << std::endl;
@@ -415,31 +438,31 @@ namespace MultiArrayTools
 		  << "](" << IB::mRangePtr << "): " << meta() << std::endl;
     }
 
-    template <typename U, SpaceType TYPE>
+    template <typename U, SpaceType TYPE, size_t S>
     template <class Expr>
-    auto SingleIndex<U,TYPE>::ifor(size_t step, Expr ex) const
-	-> For<SingleIndex<U,TYPE>,Expr>
+    auto GenSingleIndex<U,TYPE,S>::ifor(size_t step, Expr ex) const
+	-> For<GenSingleIndex<U,TYPE,S>,Expr>
     {
 	//static const size_t LAYER = typename Expr::LAYER; 
-	return For<SingleIndex<U,TYPE>,Expr>(this, step, ex);
+	return For<GenSingleIndex<U,TYPE,S>,Expr>(this, step, ex);
     }
 
-    template <typename U, SpaceType TYPE>
+    template <typename U, SpaceType TYPE, size_t S>
     template <class Expr>
-    auto SingleIndex<U,TYPE>::iforh(size_t step, Expr ex) const
-	-> For<SingleIndex<U,TYPE>,Expr,ForType::HIDDEN>
+    auto GenSingleIndex<U,TYPE,S>::iforh(size_t step, Expr ex) const
+	-> For<GenSingleIndex<U,TYPE,S>,Expr,ForType::HIDDEN>
     {
 	//static const size_t LAYER = typename Expr::LAYER; 
-	return For<SingleIndex<U,TYPE>,Expr,ForType::HIDDEN>(this, step, ex);
+	return For<GenSingleIndex<U,TYPE,S>,Expr,ForType::HIDDEN>(this, step, ex);
     }
 
-    template <typename U, SpaceType TYPE>
+    template <typename U, SpaceType TYPE, size_t S>
     template <class Expr>
-    auto SingleIndex<U,TYPE>::pifor(size_t step, Expr ex) const
-	-> PFor<SingleIndex<U,TYPE>,Expr>
+    auto GenSingleIndex<U,TYPE,S>::pifor(size_t step, Expr ex) const
+	-> PFor<GenSingleIndex<U,TYPE,S>,Expr>
     {
 	//static const size_t LAYER = typename Expr::LAYER; 
-	return PFor<SingleIndex<U,TYPE>,Expr>(this, step, ex);
+	return PFor<GenSingleIndex<U,TYPE,S>,Expr>(this, step, ex);
     }
 
     
@@ -447,14 +470,14 @@ namespace MultiArrayTools
      *   SingleRange    *
      ********************/
 
-    template <typename U, SpaceType TYPE>
-    SingleRangeFactory<U,TYPE>::SingleRangeFactory(const std::vector<U>& space)
+    template <typename U, SpaceType TYPE, size_t S>
+    GenSingleRangeFactory<U,TYPE,S>::GenSingleRangeFactory(const std::vector<U>& space)
     {
-	mProd = std::shared_ptr<oType>( new SingleRange<U,TYPE>( space ) );
+	mProd = std::shared_ptr<oType>( new GenSingleRange<U,TYPE,S>( space ) );
     }
 
-    template <typename U, SpaceType TYPE>
-    std::shared_ptr<RangeBase> SingleRangeFactory<U,TYPE>::create()
+    template <typename U, SpaceType TYPE, size_t S>
+    std::shared_ptr<RangeBase> GenSingleRangeFactory<U,TYPE,S>::create()
     {
 	setSelf();
 	return mProd;
@@ -464,9 +487,9 @@ namespace MultiArrayTools
      *   SingleRange    *
      ********************/
     
-    template <typename U, SpaceType TYPE>
-    SingleRange<U,TYPE>::SingleRange(const std::vector<U>& space) :
-	RangeInterface<SingleIndex<U,TYPE> >(),
+    template <typename U, SpaceType TYPE, size_t S>
+    GenSingleRange<U,TYPE,S>::GenSingleRange(const std::vector<U>& space) :
+	RangeInterface<GenSingleIndex<U,TYPE,S> >(),
 	mSpace(space), mMSpace(mSpace)
     {
 	//for(size_t i = 0; i != mSpace.size(); ++i){
@@ -474,50 +497,50 @@ namespace MultiArrayTools
 	//}
     }
     
-    template <typename U, SpaceType TYPE>
-    const U& SingleRange<U,TYPE>::get(size_t pos) const
+    template <typename U, SpaceType TYPE, size_t S>
+    const U& GenSingleRange<U,TYPE,S>::get(size_t pos) const
     {
 	return mSpace[pos];
     }
 
-    template <typename U, SpaceType TYPE>
-    size_t SingleRange<U,TYPE>::getMeta(const U& metaPos) const
+    template <typename U, SpaceType TYPE, size_t S>
+    size_t GenSingleRange<U,TYPE,S>::getMeta(const U& metaPos) const
     {
 	return mMSpace.at(metaPos);
     }
 
-    template <typename U, SpaceType TYPE>
-    size_t SingleRange<U,TYPE>::size() const
+    template <typename U, SpaceType TYPE, size_t S>
+    size_t GenSingleRange<U,TYPE,S>::size() const
     {
 	return mSpace.size();
     }
 
-    template <typename U, SpaceType TYPE>
-    size_t SingleRange<U,TYPE>::dim() const
+    template <typename U, SpaceType TYPE, size_t S>
+    size_t GenSingleRange<U,TYPE,S>::dim() const
     {
 	return 1;
     }
 
-    template <typename U, SpaceType TYPE>
-    bool SingleRange<U,TYPE>::isMeta(const U& metaPos) const
+    template <typename U, SpaceType TYPE, size_t S>
+    bool GenSingleRange<U,TYPE,S>::isMeta(const U& metaPos) const
     {
 	return mMSpace.count(metaPos) != 0;
     }
     
-    template <typename U, SpaceType TYPE>
-    SpaceType SingleRange<U,TYPE>::spaceType() const
+    template <typename U, SpaceType TYPE, size_t S>
+    SpaceType GenSingleRange<U,TYPE,S>::spaceType() const
     {
 	return TYPE;
     }
     
-    template <typename U, SpaceType TYPE>
-    std::string SingleRange<U,TYPE>::stringMeta(size_t pos) const
+    template <typename U, SpaceType TYPE, size_t S>
+    std::string GenSingleRange<U,TYPE,S>::stringMeta(size_t pos) const
     {
 	return xToString(get(pos));
     }
 
-    template <typename U, SpaceType TYPE>
-    std::vector<char> SingleRange<U,TYPE>::data() const
+    template <typename U, SpaceType TYPE, size_t S>
+    std::vector<char> GenSingleRange<U,TYPE,S>::data() const
     {
         DataHeader h = dataHeader();
 	std::vector<char> out;
@@ -530,8 +553,8 @@ namespace MultiArrayTools
 	return out;
     }
     
-    template <typename U, SpaceType TYPE>
-    DataHeader SingleRange<U,TYPE>::dataHeader() const
+    template <typename U, SpaceType TYPE, size_t S>
+    DataHeader GenSingleRange<U,TYPE,S>::dataHeader() const
     {
 	DataHeader h;
 	h.spaceType = static_cast<int>( TYPE );
@@ -541,19 +564,19 @@ namespace MultiArrayTools
         return h;
     }
     
-    template <typename U, SpaceType TYPE>
-    typename SingleRange<U,TYPE>::IndexType SingleRange<U,TYPE>::begin() const
+    template <typename U, SpaceType TYPE, size_t S>
+    typename GenSingleRange<U,TYPE,S>::IndexType GenSingleRange<U,TYPE,S>::begin() const
     {
-	SingleIndex<U,TYPE> i( std::dynamic_pointer_cast<SingleRange<U,TYPE> >
+	GenSingleIndex<U,TYPE,S> i( std::dynamic_pointer_cast<GenSingleRange<U,TYPE,S> >
 			       ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
 	i = 0;
 	return i;
     }
     
-    template <typename U, SpaceType TYPE>
-    typename SingleRange<U,TYPE>::IndexType SingleRange<U,TYPE>::end() const
+    template <typename U, SpaceType TYPE, size_t S>
+    typename GenSingleRange<U,TYPE,S>::IndexType GenSingleRange<U,TYPE,S>::end() const
     {
-	SingleIndex<U,TYPE> i( std::dynamic_pointer_cast<SingleRange<U,TYPE> >
+	GenSingleIndex<U,TYPE,S> i( std::dynamic_pointer_cast<GenSingleRange<U,TYPE,S> >
 			       ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
 	i = size();
 	return i;
diff --git a/src/lib/ranges/anonymous_range.cc b/src/lib/ranges/anonymous_range.cc
index 3849270..4b5a337 100644
--- a/src/lib/ranges/anonymous_range.cc
+++ b/src/lib/ranges/anonymous_range.cc
@@ -196,7 +196,7 @@ namespace MultiArrayTools
         return mOrig;
     }
     
-    SingleRange<size_t,SpaceType::ANON>::SingleRange(const std::vector<std::shared_ptr<RangeBase>>& origs) :
+    GenSingleRange<size_t,SpaceType::ANON,-1>::GenSingleRange(const std::vector<std::shared_ptr<RangeBase>>& origs) :
         RangeInterface<AnonymousIndex>(),
         mOrig(origs)
     {
diff --git a/src/lib/ranges/range_types/classic_range.cc b/src/lib/ranges/range_types/classic_range.cc
index 9cdba7d..db147d9 100644
--- a/src/lib/ranges/range_types/classic_range.cc
+++ b/src/lib/ranges/range_types/classic_range.cc
@@ -5,60 +5,60 @@
 namespace MultiArrayTools
 {
     /********************
-     *   SingleRange    *
+     *   GenSingleRange    *
      ********************/
     
-    SingleRangeFactory<size_t,SpaceType::NONE>::SingleRangeFactory(size_t size)
+    GenSingleRangeFactory<size_t,SpaceType::NONE,-1>::GenSingleRangeFactory(size_t size)
     {
 	// Quasi Singleton
 	if(not mProd){
-	    mProd = std::shared_ptr<oType>( new SingleRange<size_t,SpaceType::NONE>(size) );
+	    mProd = std::shared_ptr<oType>( new GenSingleRange<size_t,SpaceType::NONE,-1>(size) );
 	    setSelf();
 	}
     }
 
-    std::shared_ptr<RangeBase> SingleRangeFactory<size_t,SpaceType::NONE>::create()
+    std::shared_ptr<RangeBase> GenSingleRangeFactory<size_t,SpaceType::NONE,-1>::create()
     {
 	return mProd;
     }
     
     /********************
-     *   SingleRange    *
+     *   GenSingleRange    *
      ********************/
 
-    SingleRange<size_t,SpaceType::NONE>::SingleRange(size_t size) : mSize(size) {}
+    GenSingleRange<size_t,SpaceType::NONE,-1>::GenSingleRange(size_t size) : mSize(size) {}
     
-    size_t SingleRange<size_t,SpaceType::NONE>::get(size_t pos) const
+    size_t GenSingleRange<size_t,SpaceType::NONE,-1>::get(size_t pos) const
     {
 	return pos;
     }
     
-    size_t SingleRange<size_t,SpaceType::NONE>::getMeta(size_t metaPos) const
+    size_t GenSingleRange<size_t,SpaceType::NONE,-1>::getMeta(size_t metaPos) const
     {
 	return metaPos;
     }
     
-    size_t SingleRange<size_t,SpaceType::NONE>::size() const
+    size_t GenSingleRange<size_t,SpaceType::NONE,-1>::size() const
     {
 	return mSize;
     }
     
-    size_t SingleRange<size_t,SpaceType::NONE>::dim() const
+    size_t GenSingleRange<size_t,SpaceType::NONE,-1>::dim() const
     {
 	return 1;
     }
 
-    SpaceType SingleRange<size_t,SpaceType::NONE>::spaceType() const
+    SpaceType GenSingleRange<size_t,SpaceType::NONE,-1>::spaceType() const
     {
 	return SpaceType::NONE;
     }
     
-    std::string SingleRange<size_t,SpaceType::NONE>::stringMeta(size_t pos) const
+    std::string GenSingleRange<size_t,SpaceType::NONE,-1>::stringMeta(size_t pos) const
     {
 	return std::to_string(get(pos));
     }
 
-    std::vector<char> SingleRange<size_t,SpaceType::NONE>::data() const
+    std::vector<char> GenSingleRange<size_t,SpaceType::NONE,-1>::data() const
     {
 	DataHeader h = dataHeader();
 	std::vector<char> out;
@@ -70,7 +70,7 @@ namespace MultiArrayTools
 	return out;
     }
 
-    DataHeader SingleRange<size_t,SpaceType::NONE>::dataHeader() const
+    DataHeader GenSingleRange<size_t,SpaceType::NONE,-1>::dataHeader() const
     {
 	DataHeader h;
 	h.spaceType = static_cast<int>( SpaceType::NONE );
@@ -80,17 +80,17 @@ namespace MultiArrayTools
     }
 
     
-    typename SingleRange<size_t,SpaceType::NONE>::IndexType SingleRange<size_t,SpaceType::NONE>::begin() const
+    typename GenSingleRange<size_t,SpaceType::NONE,-1>::IndexType GenSingleRange<size_t,SpaceType::NONE,-1>::begin() const
     {
-	SingleIndex<size_t,SpaceType::NONE> i( std::dynamic_pointer_cast<SingleRange<size_t,SpaceType::NONE> >
+	GenSingleIndex<size_t,SpaceType::NONE,-1> i( std::dynamic_pointer_cast<GenSingleRange<size_t,SpaceType::NONE,-1> >
 			       ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
 	i = 0;
 	return i;
     }
     
-    typename SingleRange<size_t,SpaceType::NONE>::IndexType SingleRange<size_t,SpaceType::NONE>::end() const
+    typename GenSingleRange<size_t,SpaceType::NONE,-1>::IndexType GenSingleRange<size_t,SpaceType::NONE,-1>::end() const
     {
-	SingleIndex<size_t,SpaceType::NONE> i( std::dynamic_pointer_cast<SingleRange<size_t,SpaceType::NONE> >
+	GenSingleIndex<size_t,SpaceType::NONE,-1> i( std::dynamic_pointer_cast<GenSingleRange<size_t,SpaceType::NONE,-1> >
 			       ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
 	i = size();
 	return i;
diff --git a/src/lib/ranges/range_types/null_range.cc b/src/lib/ranges/range_types/null_range.cc
index d42ebb5..46ac26b 100644
--- a/src/lib/ranges/range_types/null_range.cc
+++ b/src/lib/ranges/range_types/null_range.cc
@@ -10,11 +10,11 @@ namespace MultiArrayTools
     }
 	
     /********************
-     *   SingleRange    *
+     *   GenSingleRange    *
      ********************/
 
-    std::shared_ptr<SingleRange<size_t,SpaceType::NUL>>
-    SingleRangeFactory<size_t,SpaceType::NUL>::mRInstance = nullptr;
+    std::shared_ptr<GenSingleRange<size_t,SpaceType::NUL,0>>
+    GenSingleRangeFactory<size_t,SpaceType::NUL,0>::mRInstance = nullptr;
     
     std::shared_ptr<NullRange> nullr()
     {
@@ -28,12 +28,12 @@ namespace MultiArrayTools
 	return std::make_shared<NullIndex>(nullr());
     }
 
-    SingleRangeFactory<size_t,SpaceType::NUL>::SingleRangeFactory()
+    GenSingleRangeFactory<size_t,SpaceType::NUL,0>::GenSingleRangeFactory()
     {
 	// Singleton
 	if(not mRInstance){
 	    if(not mProd){
-		mProd = std::shared_ptr<oType>( new SingleRange<size_t,SpaceType::NUL>() );
+		mProd = std::shared_ptr<oType>( new GenSingleRange<size_t,SpaceType::NUL,0>() );
 		setSelf();
 	    }
 	    mRInstance = std::dynamic_pointer_cast<NullRange>( mProd );
@@ -42,46 +42,46 @@ namespace MultiArrayTools
 	}
     }
 
-    std::shared_ptr<RangeBase> SingleRangeFactory<size_t,SpaceType::NUL>::create()
+    std::shared_ptr<RangeBase> GenSingleRangeFactory<size_t,SpaceType::NUL,0>::create()
     {
 	return mProd;
     }
     
     /********************
-     *   SingleRange    *
+     *   GenSingleRange    *
      ********************/
     
-    size_t SingleRange<size_t,SpaceType::NUL>::get(size_t pos) const
+    size_t GenSingleRange<size_t,SpaceType::NUL,0>::get(size_t pos) const
     {
 	return 0;
     }
     
-    size_t SingleRange<size_t,SpaceType::NUL>::getMeta(size_t metapos) const
+    size_t GenSingleRange<size_t,SpaceType::NUL,0>::getMeta(size_t metapos) const
     {
 	return 0;
     }
     
-    size_t SingleRange<size_t,SpaceType::NUL>::size() const
+    size_t GenSingleRange<size_t,SpaceType::NUL,0>::size() const
     {
 	return 1;
     }
     
-    size_t SingleRange<size_t,SpaceType::NUL>::dim() const
+    size_t GenSingleRange<size_t,SpaceType::NUL,0>::dim() const
     {
 	return 1;
     }
 
-    SpaceType SingleRange<size_t,SpaceType::NUL>::spaceType() const
+    SpaceType GenSingleRange<size_t,SpaceType::NUL,0>::spaceType() const
     {
 	return SpaceType::NUL;
     }
     
-    std::string SingleRange<size_t,SpaceType::NUL>::stringMeta(size_t pos) const
+    std::string GenSingleRange<size_t,SpaceType::NUL,0>::stringMeta(size_t pos) const
     {
 	return std::to_string(get(pos));
     }
 
-    std::vector<char> SingleRange<size_t,SpaceType::NUL>::data() const
+    std::vector<char> GenSingleRange<size_t,SpaceType::NUL,0>::data() const
     {
 	DataHeader h = dataHeader(); 
 	std::vector<char> out;
@@ -91,7 +91,7 @@ namespace MultiArrayTools
 	return out;
     }
 
-    DataHeader SingleRange<size_t,SpaceType::NUL>::dataHeader() const
+    DataHeader GenSingleRange<size_t,SpaceType::NUL,0>::dataHeader() const
     {
 	DataHeader h;
 	h.spaceType = static_cast<int>( SpaceType::NUL );
@@ -100,17 +100,17 @@ namespace MultiArrayTools
         return h;
     }
    
-    typename SingleRange<size_t,SpaceType::NUL>::IndexType SingleRange<size_t,SpaceType::NUL>::begin() const
+    typename GenSingleRange<size_t,SpaceType::NUL,0>::IndexType GenSingleRange<size_t,SpaceType::NUL,0>::begin() const
     {
-	SingleIndex<size_t,SpaceType::NUL> i( std::dynamic_pointer_cast<SingleRange<size_t,SpaceType::NUL> >
+	GenSingleIndex<size_t,SpaceType::NUL,0> i( std::dynamic_pointer_cast<GenSingleRange<size_t,SpaceType::NUL,0> >
 			       ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
 	i = 0;
 	return i;
     }
     
-    typename SingleRange<size_t,SpaceType::NUL>::IndexType SingleRange<size_t,SpaceType::NUL>::end() const
+    typename GenSingleRange<size_t,SpaceType::NUL,0>::IndexType GenSingleRange<size_t,SpaceType::NUL,0>::end() const
     {
-	SingleIndex<size_t,SpaceType::NUL> i( std::dynamic_pointer_cast<SingleRange<size_t,SpaceType::NUL> >
+	GenSingleIndex<size_t,SpaceType::NUL,0> i( std::dynamic_pointer_cast<GenSingleRange<size_t,SpaceType::NUL,0> >
 			       ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
 	i = size();
 	return i;
@@ -118,12 +118,12 @@ namespace MultiArrayTools
 
     // put this in the interface class !!!
     /*
-    std::shared_ptr<VIWB> SingleRange<size_t,SpaceType::NUL>::index() const
+    std::shared_ptr<VIWB> GenSingleRange<size_t,SpaceType::NUL,0>::index() const
     {
 	typedef IndexWrapper<IndexType> IW;
 	return std::make_shared<IW>
             ( std::make_shared<IndexType>
-	      ( std::dynamic_pointer_cast<SingleRange<size_t,SpaceType::NUL> >
+	      ( std::dynamic_pointer_cast<GenSingleRange<size_t,SpaceType::NUL,0> >
 		( std::shared_ptr<RangeBase>( RB::mThis ) ) ) );
     }
     */
diff --git a/src/lib/ranges/range_types/space_range.cc b/src/lib/ranges/range_types/space_range.cc
index e4ca7f3..ae4758e 100644
--- a/src/lib/ranges/range_types/space_range.cc
+++ b/src/lib/ranges/range_types/space_range.cc
@@ -11,60 +11,60 @@ namespace MultiArrayTools
     }
     
     /********************
-     *   SingleRange    *
+     *   GenSingleRange    *
      ********************/
     
-    SingleRangeFactory<int,SpaceType::PSPACE>::SingleRangeFactory(size_t size)
+    GenSingleRangeFactory<int,SpaceType::PSPACE,-1>::GenSingleRangeFactory(size_t size)
     {
 	// Quasi Singleton
 	if(not mProd){
-	    mProd = std::shared_ptr<oType>( new SingleRange<int,SpaceType::PSPACE>(size) );
+	    mProd = std::shared_ptr<oType>( new GenSingleRange<int,SpaceType::PSPACE,-1>(size) );
 	    setSelf();
 	}
     }
 
-    std::shared_ptr<RangeBase> SingleRangeFactory<int,SpaceType::PSPACE>::create()
+    std::shared_ptr<RangeBase> GenSingleRangeFactory<int,SpaceType::PSPACE,-1>::create()
     {
 	return mProd;
     }
     
     /********************
-     *   SingleRange    *
+     *   GenSingleRange    *
      ********************/
 
-    SingleRange<int,SpaceType::PSPACE>::SingleRange(size_t size) : mSize(size) { }
+    GenSingleRange<int,SpaceType::PSPACE,-1>::GenSingleRange(size_t size) : mSize(size) { }
     
-    int SingleRange<int,SpaceType::PSPACE>::get(size_t pos) const
+    int GenSingleRange<int,SpaceType::PSPACE,-1>::get(size_t pos) const
     {
 	return pos > mSize / 2 ? pos - mSize : pos;
     }
     
-    size_t SingleRange<int,SpaceType::PSPACE>::getMeta(int metaPos) const
+    size_t GenSingleRange<int,SpaceType::PSPACE,-1>::getMeta(int metaPos) const
     {
 	return metaPos < 0 ? metaPos + mSize : metaPos;
     }
     
-    size_t SingleRange<int,SpaceType::PSPACE>::size() const
+    size_t GenSingleRange<int,SpaceType::PSPACE,-1>::size() const
     {
 	return mSize;
     }
     
-    size_t SingleRange<int,SpaceType::PSPACE>::dim() const
+    size_t GenSingleRange<int,SpaceType::PSPACE,-1>::dim() const
     {
 	return 1;
     }
 
-    SpaceType SingleRange<int,SpaceType::PSPACE>::spaceType() const
+    SpaceType GenSingleRange<int,SpaceType::PSPACE,-1>::spaceType() const
     {
 	return SpaceType::PSPACE;
     }
     
-    std::string SingleRange<int,SpaceType::PSPACE>::stringMeta(size_t pos) const
+    std::string GenSingleRange<int,SpaceType::PSPACE,-1>::stringMeta(size_t pos) const
     {
 	return std::to_string(get(pos));
     }
 
-    std::vector<char> SingleRange<int,SpaceType::PSPACE>::data() const
+    std::vector<char> GenSingleRange<int,SpaceType::PSPACE,-1>::data() const
     {
 	DataHeader h = dataHeader();
 	std::vector<char> out;
@@ -76,7 +76,7 @@ namespace MultiArrayTools
 	return out;
     }
 
-    DataHeader SingleRange<int,SpaceType::PSPACE>::dataHeader() const
+    DataHeader GenSingleRange<int,SpaceType::PSPACE,-1>::dataHeader() const
     {
 	DataHeader h;
 	h.spaceType = static_cast<int>( SpaceType::PSPACE );
@@ -85,17 +85,17 @@ namespace MultiArrayTools
         return h;
     }
     
-    typename SingleRange<int,SpaceType::PSPACE>::IndexType SingleRange<int,SpaceType::PSPACE>::begin() const
+    typename GenSingleRange<int,SpaceType::PSPACE,-1>::IndexType GenSingleRange<int,SpaceType::PSPACE,-1>::begin() const
     {
-	SingleIndex<int,SpaceType::PSPACE> i( std::dynamic_pointer_cast<SingleRange<int,SpaceType::PSPACE> >
+	GenSingleIndex<int,SpaceType::PSPACE,-1> i( std::dynamic_pointer_cast<GenSingleRange<int,SpaceType::PSPACE,-1> >
 			       ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
 	i = 0;
 	return i;
     }
     
-    typename SingleRange<int,SpaceType::PSPACE>::IndexType SingleRange<int,SpaceType::PSPACE>::end() const
+    typename GenSingleRange<int,SpaceType::PSPACE,-1>::IndexType GenSingleRange<int,SpaceType::PSPACE,-1>::end() const
     {
-	SingleIndex<int,SpaceType::PSPACE> i( std::dynamic_pointer_cast<SingleRange<int,SpaceType::PSPACE> >
+	GenSingleIndex<int,SpaceType::PSPACE,-1> i( std::dynamic_pointer_cast<GenSingleRange<int,SpaceType::PSPACE,-1> >
 			       ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
 	i = size();
 	return i;
diff --git a/src/lib/ranges/range_types/spin_range.cc b/src/lib/ranges/range_types/spin_range.cc
index b391481..8ba73f4 100644
--- a/src/lib/ranges/range_types/spin_range.cc
+++ b/src/lib/ranges/range_types/spin_range.cc
@@ -10,58 +10,58 @@ namespace MultiArrayTools
     }
     
     /********************
-     *   SingleRange    *
+     *   GenSingleRange    *
      ********************/
     
-    SingleRangeFactory<size_t,SpaceType::SPIN>::SingleRangeFactory()
+    GenSingleRangeFactory<size_t,SpaceType::SPIN,4>::GenSingleRangeFactory()
     {
 	// Quasi Singleton
 	if(not mProd){
-	    mProd = std::shared_ptr<oType>( new SingleRange<size_t,SpaceType::SPIN>() );
+	    mProd = std::shared_ptr<oType>( new GenSingleRange<size_t,SpaceType::SPIN,4>() );
 	    setSelf();
 	}
     }
 
-    std::shared_ptr<RangeBase> SingleRangeFactory<size_t,SpaceType::SPIN>::create()
+    std::shared_ptr<RangeBase> GenSingleRangeFactory<size_t,SpaceType::SPIN,4>::create()
     {
 	return mProd;
     }
     
     /********************
-     *   SingleRange    *
+     *   GenSingleRange    *
      ********************/
     
-    size_t SingleRange<size_t,SpaceType::SPIN>::get(size_t pos) const
+    size_t GenSingleRange<size_t,SpaceType::SPIN,4>::get(size_t pos) const
     {
 	return pos;
     }
     
-    size_t SingleRange<size_t,SpaceType::SPIN>::getMeta(size_t metaPos) const
+    size_t GenSingleRange<size_t,SpaceType::SPIN,4>::getMeta(size_t metaPos) const
     {
 	return metaPos;
     }
     
-    size_t SingleRange<size_t,SpaceType::SPIN>::size() const
+    size_t GenSingleRange<size_t,SpaceType::SPIN,4>::size() const
     {
 	return mSpinNum;
     }
     
-    size_t SingleRange<size_t,SpaceType::SPIN>::dim() const
+    size_t GenSingleRange<size_t,SpaceType::SPIN,4>::dim() const
     {
 	return 1;
     }
 
-    SpaceType SingleRange<size_t,SpaceType::SPIN>::spaceType() const
+    SpaceType GenSingleRange<size_t,SpaceType::SPIN,4>::spaceType() const
     {
 	return SpaceType::SPIN;
     }
     
-    std::string SingleRange<size_t,SpaceType::SPIN>::stringMeta(size_t pos) const
+    std::string GenSingleRange<size_t,SpaceType::SPIN,4>::stringMeta(size_t pos) const
     {
 	return std::to_string(get(pos));
     }
 
-    std::vector<char> SingleRange<size_t,SpaceType::SPIN>::data() const
+    std::vector<char> GenSingleRange<size_t,SpaceType::SPIN,4>::data() const
     {
 	DataHeader h = dataHeader();
 	std::vector<char> out;
@@ -71,7 +71,7 @@ namespace MultiArrayTools
 	return out;
     }
     
-    DataHeader SingleRange<size_t,SpaceType::SPIN>::dataHeader() const
+    DataHeader GenSingleRange<size_t,SpaceType::SPIN,4>::dataHeader() const
     {
 	DataHeader h;
 	h.spaceType = static_cast<int>( SpaceType::SPIN );
@@ -80,17 +80,17 @@ namespace MultiArrayTools
         return h;
     }    
     
-    typename SingleRange<size_t,SpaceType::SPIN>::IndexType SingleRange<size_t,SpaceType::SPIN>::begin() const
+    typename GenSingleRange<size_t,SpaceType::SPIN,4>::IndexType GenSingleRange<size_t,SpaceType::SPIN,4>::begin() const
     {
-	SingleIndex<size_t,SpaceType::SPIN> i( std::dynamic_pointer_cast<SingleRange<size_t,SpaceType::SPIN> >
+	GenSingleIndex<size_t,SpaceType::SPIN,4> i( std::dynamic_pointer_cast<GenSingleRange<size_t,SpaceType::SPIN,4> >
 			       ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
 	i = 0;
 	return i;
     }
     
-    typename SingleRange<size_t,SpaceType::SPIN>::IndexType SingleRange<size_t,SpaceType::SPIN>::end() const
+    typename GenSingleRange<size_t,SpaceType::SPIN,4>::IndexType GenSingleRange<size_t,SpaceType::SPIN,4>::end() const
     {
-	SingleIndex<size_t,SpaceType::SPIN> i( std::dynamic_pointer_cast<SingleRange<size_t,SpaceType::SPIN> >
+	GenSingleIndex<size_t,SpaceType::SPIN,4> i( std::dynamic_pointer_cast<GenSingleRange<size_t,SpaceType::SPIN,4> >
 			       ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
 	i = size();
 	return i;
@@ -98,12 +98,12 @@ namespace MultiArrayTools
 
     // put this in the interface class !!!
     /*
-    std::shared_ptr<VIWB> SingleRange<size_t,SpaceType::SPIN>::index() const
+    std::shared_ptr<VIWB> GenSingleRange<size_t,SpaceType::SPIN,4>::index() const
     {
 	typedef IndexWrapper<IndexType> IW;
 	return std::make_shared<IW>
             ( std::make_shared<IndexType>
-	      ( std::dynamic_pointer_cast<SingleRange<size_t,SpaceType::SPIN> >
+	      ( std::dynamic_pointer_cast<GenSingleRange<size_t,SpaceType::SPIN,4> >
 		( std::shared_ptr<RangeBase>( RB::mThis ) ) ) );
     }
     */
diff --git a/src/tests/op_perf_test.cc b/src/tests/op_perf_test.cc
index 4c623ec..62accee 100644
--- a/src/tests/op_perf_test.cc
+++ b/src/tests/op_perf_test.cc
@@ -134,13 +134,14 @@ namespace {
 	typedef MultiRangeFactory<SR,SR,SR,SR,SR,SR,SR,SR> SR8F;
 	typedef SR8F::oType SR8;
 
-        static const size_t os = 3000;
-	static const size_t s = 65536*os;
-	
+        static const size_t os = 30;
+        static const size_t is = 65536;
+        static const size_t s = is*os;
+        
 	OpTest_Spin()
 	{
-	    data.resize(s);
-	    for(size_t i = 0; i != s; ++i){
+	    data.resize(is);
+	    for(size_t i = 0; i != is; ++i){
 		double arg = static_cast<double>( i - s ) - 0.1; 
 		data[i] = sin(arg);
 		//VCHECK(data[i]);
@@ -162,7 +163,7 @@ namespace {
 
     void OpTest_Spin::contract()
     {
-	MultiArray<double,CR,SR,SR,SR,SR,SR,SR,SR,SR> ma( cr, sr, sr, sr, sr, sr, sr, sr, sr, data);
+	MultiArray<double,SR,SR,SR,SR,SR,SR,SR,SR> ma( sr, sr, sr, sr, sr, sr, sr, sr, data);
 	MultiArray<double,CR,SR,SR> res1( cr, sr, sr );
 
         auto ii = MAT::getIndex<CR>(cr);
@@ -177,7 +178,8 @@ namespace {
 
 	std::clock_t begin = std::clock();
 	//for(size_t i = 0; i != os; ++i){
-        res1(ii ,delta, deltap).par() += ma(ii, delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix);
+        //res1(ii ,delta, deltap).par() += ma(ii, delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix);
+        res1(ii ,delta, deltap) += ma(delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix);
             //}
 	std::clock_t end = std::clock();
 	std::cout << "MultiArray time: " << static_cast<double>( end - begin ) / CLOCKS_PER_SEC
@@ -199,7 +201,7 @@ namespace {
                             for(size_t d = 0; d != 4; ++d){
                                 for(size_t p = 0; p != 4; ++p){
                                     const size_t tidx = i*4*4 + d*4 + p;
-                                    const size_t sidx = i*65536 + d*4*4*4*4*4*4*4 + a*5*4*4*4*4*4 + b*5*4*4*4 + c*5*4 + p;
+                                    const size_t sidx = /*i*65536 +*/ d*4*4*4*4*4*4*4 + a*5*4*4*4*4*4 + b*5*4*4*4 + c*5*4 + p;
                                     vres[tidx] += data[sidx];
                                 }
                             }

From 0326ebc5050dd8a22dc9e4146353bd966a0a42d1 Mon Sep 17 00:00:00 2001
From: Christian Zimmermann <christian.zimmermann@ur.de>
Date: Tue, 12 Feb 2019 22:34:09 +0100
Subject: [PATCH 2/9] conversions

---
 src/include/conversions.h     |  64 +++++++++++++++++++
 src/include/type_operations.h | 112 ++++++++++++++++++++++++++++------
 2 files changed, 156 insertions(+), 20 deletions(-)
 create mode 100644 src/include/conversions.h

diff --git a/src/include/conversions.h b/src/include/conversions.h
new file mode 100644
index 0000000..87d1a1e
--- /dev/null
+++ b/src/include/conversions.h
@@ -0,0 +1,64 @@
+
+#ifndef __ma_conversions_h__
+#define __ma_conversions_h__
+
+#include "multi_array.h"
+#include "slice.h"
+
+namespace MultiArrayTools
+{
+
+    namespace ConversionSizes
+    {
+	template <size_t N>
+	struct OrigSize
+	{
+	    template <typename C, typename T>
+	    struct FromTo
+	    {
+		static void check() { static_assert( not N % (sizeof(T) / sizeof(C)), "conversion does not fit" ); }
+		static constexpr size_t SIZE = N * sizeof(T) / sizeof(C);
+	    };
+	};
+
+	template <>
+	struct OrigSize<-1>
+	{
+	    template <typename C, typename T>
+	    struct FromTo
+	    {
+		static void check() {}
+		static constexpr size_t SIZE = -1;
+	    };
+	};
+    }
+
+    namespace
+    {
+	template <typename C, typename T, size_t N>
+	using SC = ConversionSizes::OrigSize<N>::template FromTo<C,T>;
+
+	template <typename C, typename T, class Range>
+	using SCR = SC<C,T,typename Range::SIZE>;
+
+	template <typename C, typename T, class Range>
+	using SCRR = GenSingleRange<typename Range::value_type,SpaceType::NONE,SCR<C,T,typename Range::SIZE>;
+    }
+    
+    template <typename C, typename T, class Range, class... Ranges>
+    Slice<C,Ranges...,SCRR<C,T,Range>> tcast(MultiArray<T,Ranges...,Range>& ma)
+    {
+	return Slice<C,Ranges...,SCRR<C,T,Range>>
+	    ( ma.range()->space(), reinterpret_cast<C*>( ma.data() ) );
+    }
+    
+    template <typename C, typename T, class Range, class... Ranges>
+    ConstSlice<C,Ranges...,SCRR<C,T,Range>> tcast(const MultiArray<T,Ranges...,Range>& ma)
+    {
+	return ConstSlice<C,Ranges...,SCRR<C,T,Range>>
+	    ( ma.range()->space(), reinterpret_cast<const C*>( ma.data() ) );
+    }
+
+}
+
+#endif
diff --git a/src/include/type_operations.h b/src/include/type_operations.h
index 6870dac..efa35c6 100644
--- a/src/include/type_operations.h
+++ b/src/include/type_operations.h
@@ -82,7 +82,60 @@ namespace MultiArrayTools
 	std::transform(a.begin(), a.end(), b.begin(), a.begin(), std::plus<T>());
 	return a;
     }
-    
+
+    template <typename T>
+    std::vector<T>& operator-=(std::vector<T>& a, const std::vector<T>& b)
+    {
+	std::transform(a.begin(), a.end(), b.begin(), a.begin(), std::minus<T>());
+	return a;
+    }
+
+    template <typename T>
+    std::vector<T>& operator*=(std::vector<T>& a, const std::vector<T>& b)
+    {
+	std::transform(a.begin(), a.end(), b.begin(), a.begin(), std::multiplies<T>());
+	return a;
+    }
+
+    template <typename T>
+    std::vector<T>& operator/=(std::vector<T>& a, const std::vector<T>& b)
+    {
+	std::transform(a.begin(), a.end(), b.begin(), a.begin(), std::divides<T>());
+	return a;
+    }
+
+    template <typename T>
+    std::vector<T> operator+(std::vector<T>& a, const std::vector<T>& b)
+    {
+	std::vector<T> out(a.size());
+	std::transform(a.begin(), a.end(), b.begin(), out.begin(), std::plus<T>());
+	return out;
+    }
+
+    template <typename T>
+    std::vector<T> operator-(std::vector<T>& a, const std::vector<T>& b)
+    {
+	std::vector<T> out(a.size());
+	std::transform(a.begin(), a.end(), b.begin(), out.begin(), std::minus<T>());
+	return out;
+    }
+
+    template <typename T>
+    std::vector<T> operator*(std::vector<T>& a, const std::vector<T>& b)
+    {
+	std::vector<T> out(a.size());
+	std::transform(a.begin(), a.end(), b.begin(), out.begin(), std::multiplies<T>());
+	return out;
+    }
+
+    template <typename T>
+    std::vector<T> operator/(std::vector<T>& a, const std::vector<T>& b)
+    {
+	std::vector<T> out(a.size());
+	std::transform(a.begin(), a.end(), b.begin(), out.begin(), std::divides<T>());
+	return out;
+    }
+
     template <class OperationClass, typename T>
     class OperationTemplate<std::vector<T>,OperationClass> : public OperationBase<std::vector<T>,OperationClass>
     {
@@ -101,28 +154,47 @@ namespace MultiArrayTools
 	friend OperationClass;
     };
 
-    inline std::array<int,2>& operator+=(std::array<int,2>& a, const std::array<int,2>& b)
+    template <typename T, size_t N>
+    inline std::array<T,N> operator+(std::array<T,N>& a, const std::array<T,N>& b)
     {
-	std::get<0>(a) += std::get<0>(b);
-	std::get<1>(a) += std::get<1>(b);
-	return a;
+	std::array<T,N> out;
+	for(size_t i = 0; i != N; ++i){
+	    out[i] = a[i] + b[i];
+	}
+	return out;
     }
+
+    template <typename T, size_t N>
+    inline std::array<T,N> operator-(std::array<T,N>& a, const std::array<T,N>& b)
+    {
+	std::array<T,N> out;
+	for(size_t i = 0; i != N; ++i){
+	    out[i] = a[i] - b[i];
+	}
+	return out;
+    }
+
+    template <typename T, size_t N>
+    inline std::array<T,N> operator*(std::array<T,N>& a, const std::array<T,N>& b)
+    {
+	std::array<T,N> out;
+	for(size_t i = 0; i != N; ++i){
+	    out[i] = a[i] * b[i];
+	}
+	return out;
+    }
+
+    template <typename T, size_t N>
+    inline std::array<T,N> operator/(std::array<T,N>& a, const std::array<T,N>& b)
+    {
+	std::array<T,N> out;
+	for(size_t i = 0; i != N; ++i){
+	    out[i] = a[i] / b[i];
+	}
+	return out;
+    }
+
     
-    inline std::array<int,3>& operator+=(std::array<int,3>& a, const std::array<int,3>& b)
-    {
-	std::get<0>(a) += std::get<0>(b);
-	std::get<1>(a) += std::get<1>(b);
-	std::get<2>(a) += std::get<2>(b);
-	return a;
-    }
-    
-    inline std::tuple<int,int,int>& operator+=(std::tuple<int,int,int>& a, const std::tuple<int,int,int>& b)
-    {
-	std::get<0>(a) += std::get<0>(b);
-	std::get<1>(a) += std::get<1>(b);
-	std::get<2>(a) += std::get<2>(b);
-	return a;
-    }
 
 } // namespace MultiArrayTools
 

From 6d776f853dd9e22951e4686a3616af6b4d85381a Mon Sep 17 00:00:00 2001
From: Christian Zimmermann <christian.zimmermann@ur.de>
Date: Wed, 13 Feb 2019 21:59:13 +0100
Subject: [PATCH 3/9] im com

---
 CMakeLists.txt                                |  2 +-
 src/include/allocator.h                       | 74 +++++++++++++++
 src/include/allocator.h~                      | 62 ++++++++++++
 src/include/conversions.h                     | 94 +++++++++++++++----
 src/include/functional_multi_array.h          |  2 +-
 src/include/helper_tools.cc.h                 |  2 +-
 src/include/helper_tools.h                    |  2 +-
 src/include/map_range.cc.h                    | 14 +--
 src/include/map_range.h                       |  2 +-
 src/include/map_range_factory_product_map.h   |  2 +-
 src/include/mbase_def.h                       |  1 +
 src/include/multi_array.cc.h                  |  8 +-
 src/include/multi_array.h                     | 19 ++--
 src/include/ranges/anonymous_range.h          | 20 ++--
 src/include/ranges/dynamic_meta.h             |  8 +-
 src/include/ranges/dynamic_range.cc.h         | 20 ++--
 src/include/ranges/dynamic_range.h            | 22 ++---
 src/include/ranges/multi_range.h              |  8 +-
 .../ranges/multi_range_factory_product_map.h  |  2 +-
 src/include/ranges/range_base.h               |  2 +-
 .../ranges/range_types/classic_range.h        |  2 +-
 src/include/ranges/range_types/null_range.h   |  2 +-
 src/include/ranges/range_types/space_range.h  |  2 +-
 src/include/ranges/range_types/spin_range.h   |  2 +-
 src/include/ranges/rbase_def.h                |  2 +-
 src/include/ranges/rpack_num.h                | 36 +++----
 src/include/ranges/single_range.h             | 62 +++++++++---
 src/include/ranges/subrange.h                 | 24 ++---
 src/include/ranges/type_map.h                 | 22 +++--
 src/include/ranges/value_range.h              |  6 +-
 src/include/ranges/x_to_string.h              |  4 +-
 src/include/slice.cc.h                        | 11 ++-
 src/include/slice.h                           |  2 +
 src/include/type_operations.h                 | 93 +++++++++---------
 src/include/xfor/xfor.h                       | 13 +--
 src/lib/map_range_factory_product_map.cc      |  2 +-
 src/lib/ranges/anonymous_range.cc             | 18 ++--
 .../ranges/multi_range_factory_product_map.cc |  2 +-
 src/lib/ranges/range_base.cc                  |  4 +-
 src/lib/ranges/range_types/classic_range.cc   |  4 +-
 src/lib/ranges/range_types/null_range.cc      |  4 +-
 src/lib/ranges/range_types/space_range.cc     |  4 +-
 src/lib/ranges/range_types/spin_range.cc      |  4 +-
 src/tests/ma_unit_test.cc                     |  8 +-
 src/tests/op_perf_test.cc                     | 52 +++++-----
 src/tests/op_unit_test.cc                     | 46 ++++-----
 src/tests/ranges/anonymous_unit_test.cc       | 14 +--
 src/tests/ranges/index_unit_test.cc           |  2 +-
 48 files changed, 533 insertions(+), 280 deletions(-)
 create mode 100644 src/include/allocator.h
 create mode 100644 src/include/allocator.h~

diff --git a/CMakeLists.txt b/CMakeLists.txt
index a66ce5c..77f8ffb 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 2.8)
 
 project(multi_array)
 
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Werror -std=c++11 -Wpedantic -O3 -g -march=native -fopenmp -funroll-loops -ftree-vectorize")
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Werror -std=c++11 -g -Wpedantic -Ofast -mavx2 -funroll-loops -fopenmp")
 
 enable_testing()
 
diff --git a/src/include/allocator.h b/src/include/allocator.h
new file mode 100644
index 0000000..2cfdd2c
--- /dev/null
+++ b/src/include/allocator.h
@@ -0,0 +1,74 @@
+
+#ifndef __ma_allocator__
+#define __ma_allocator__
+
+#include <cstdlib>
+#include <new>
+#include <vector>
+#include <cstdint>
+#include <cassert>
+#include <iostream>
+
+namespace MultiArrayHelper
+{
+
+    template <typename T>
+    struct Allocator
+    {
+	typedef T value_type;
+	static constexpr size_t type_size = sizeof(value_type);	
+	static constexpr size_t N = 32;
+	
+	Allocator() = default;
+
+	template <typename U>
+	Allocator(const Allocator<U>& x) {}
+
+	T* allocate(size_t n)
+	{
+	    //constexpr size_t NN = N*type_size;
+	    const size_t nn = n*type_size;
+	    if(n > static_cast<size_t>(-1-N)) throw std::bad_alloc();
+	    auto p = static_cast<char*>(std::malloc(nn + N));
+	    if(not p) throw std::bad_alloc();
+	    auto ip = reinterpret_cast<std::intptr_t>(p);
+	    auto diff = N - (ip % N);
+	    p += diff;
+	    //auto ipx = reinterpret_cast<std::intptr_t>(p);
+	    //std::cout << "ixp: " << ipx << std::endl;
+	    //std::cout << "ixp %: " << ipx % N << std::endl;
+	    auto ps = reinterpret_cast<std::intptr_t*>(p);
+	    ps[-1] = diff;
+	    return reinterpret_cast<T*>(p);
+	}
+
+	void deallocate(T* p, size_t n)
+	{
+	    auto ps = reinterpret_cast<std::intptr_t*>(p);
+	    std::free(reinterpret_cast<char*>(p)-ps[-1]);
+	}
+    };
+
+    template <class T, class U>
+    bool operator==(const Allocator<T>& a, const Allocator<U>& b)
+    {
+	return true;
+    }
+
+    template <class T, class U>
+    bool operator!=(const Allocator<T>& a, const Allocator<U>& b)
+    {
+	return false;
+    }
+
+    
+} // namespace MultiArrayHelper
+
+namespace MultiArrayTools
+{
+    template <typename T>
+    using vector = std::vector<T,MultiArrayHelper::Allocator<T>>;
+
+} // namespace MultiArrayTools
+
+#endif
diff --git a/src/include/allocator.h~ b/src/include/allocator.h~
new file mode 100644
index 0000000..9636d7a
--- /dev/null
+++ b/src/include/allocator.h~
@@ -0,0 +1,62 @@
+
+#ifndef __ma_allocator__
+#define __ma_allocator__
+
+#include <cstdlib>
+#include <new>
+#include <vector>
+
+namespace MultiArrayHelper
+{
+
+    template <typename T, size_t N>
+    struct Allocator
+    {
+	typedef T value_type;
+	static constexpr type_size = sizeof(value_type);	
+	
+	Allocator() = default;
+
+	template <typename U>
+	Allocator(const Allocator<U>& x) {}
+
+	T* allocate(size_t n)
+	{
+	    //constexpr size_t NN = N*type_size;
+	    const size_t nn = n*type_size;
+	    if(n > static_cast<size_t>(-1-NN)) throw std::bad_alloc();
+	    auto p = static_cast<char*>(std::malloc(nn + N));
+	    if(not p) throw std::bad_alloc();
+	    auto ip = static_cast<std::intptr_t>(p);
+	    mDiff = N - (ip % N);
+	    p += mDiff;
+	    return static_cast<T*>(p);
+	}
+
+	void deacllocate(T* p, size_t n)
+	{
+	    std::free(p-mDiff);
+	}
+	
+    private:
+	std::intptr_t mDiff;
+    };
+
+    template <class T, class U>
+    bool operator==(const Allocator<T>& a, const Allocator<T>& b)
+    {
+	return true;
+    }
+
+    template <class T, class U>
+    bool operator!=(const Allocator<T>& a, const Allocator<T>& b)
+    {
+	return false;
+    }
+
+    template <typename T>
+    using vector = vector<T,Allocator>;
+    
+} // namespace MultiArrayHelper
+
+#endif
diff --git a/src/include/conversions.h b/src/include/conversions.h
index 87d1a1e..60a3b53 100644
--- a/src/include/conversions.h
+++ b/src/include/conversions.h
@@ -36,29 +36,87 @@ namespace MultiArrayTools
     namespace
     {
 	template <typename C, typename T, size_t N>
-	using SC = ConversionSizes::OrigSize<N>::template FromTo<C,T>;
+	using SC = typename ConversionSizes::OrigSize<N>::template FromTo<C,T>;
 
 	template <typename C, typename T, class Range>
-	using SCR = SC<C,T,typename Range::SIZE>;
+	using SCR = SC<C,T,Range::SIZE>;
 
 	template <typename C, typename T, class Range>
-	using SCRR = GenSingleRange<typename Range::value_type,SpaceType::NONE,SCR<C,T,typename Range::SIZE>;
-    }
-    
-    template <typename C, typename T, class Range, class... Ranges>
-    Slice<C,Ranges...,SCRR<C,T,Range>> tcast(MultiArray<T,Ranges...,Range>& ma)
-    {
-	return Slice<C,Ranges...,SCRR<C,T,Range>>
-	    ( ma.range()->space(), reinterpret_cast<C*>( ma.data() ) );
-    }
-    
-    template <typename C, typename T, class Range, class... Ranges>
-    ConstSlice<C,Ranges...,SCRR<C,T,Range>> tcast(const MultiArray<T,Ranges...,Range>& ma)
-    {
-	return ConstSlice<C,Ranges...,SCRR<C,T,Range>>
-	    ( ma.range()->space(), reinterpret_cast<const C*>( ma.data() ) );
+	using SCRR = GenSingleRange<typename Range::MetaType,SpaceType::NONE,SCR<C,T,Range>::SIZE>;
     }
 
+    template <size_t N>
+    struct SubTuple
+    {
+	template <class RTP, class... Ranges>
+	static inline auto mk(const RTP& rtp, const Ranges&... rs)
+	    -> decltype(SubTuple<N-1>::mk(rtp, std::get<N>(rtp), rs...))
+	{
+	    return SubTuple<N-1>::mk(rtp, std::get<N>(rtp), rs...);
+	}
+    };
+
+    template <>
+    struct SubTuple<0>
+    {
+	template <class RTP, class... Ranges>
+	static inline auto mk(const RTP& rtp, const Ranges&... rs)
+	    -> decltype(std::make_tuple(std::get<0>(rtp), rs...))
+	{
+	    return std::make_tuple(std::get<0>(rtp), rs...);
+	}
+    };
+
+    template <class... Ranges>
+    using LastR = typename std::tuple_element<sizeof...(Ranges)-1,std::tuple<Ranges...>>::type;
+    
+    template <typename C, typename T, class... Ranges>
+    auto rtcast(const std::tuple<std::shared_ptr<Ranges>...>& rtp)
+	-> decltype(std::tuple_cat(SubTuple<sizeof...(Ranges)-2>::mk(rtp),
+				   std::make_tuple( std::dynamic_pointer_cast<SCRR<C,T,LastR<Ranges...>>>
+						    ( SCRR<C,T,LastR<Ranges...>>::factory().create() ) ) ))
+    {
+	return std::tuple_cat(SubTuple<sizeof...(Ranges)-2>::mk(rtp),
+				   std::make_tuple( std::dynamic_pointer_cast<SCRR<C,T,LastR<Ranges...>>>
+						    ( SCRR<C,T,LastR<Ranges...>>::factory().create() ) ) );
+    }
+
+    template <typename T, class... Ranges>
+    inline Slice<T,Ranges...> rangeTpToSlice( const std::tuple<std::shared_ptr<Ranges>...>& rtp, T* data )
+    {
+	return Slice<T,Ranges...>(rtp, data);
+    }
+
+    template <typename T, class... Ranges>
+    inline ConstSlice<T,Ranges...> rangeTpToSlice( const std::tuple<std::shared_ptr<Ranges>...>& rtp, const T* data )
+    {
+	return ConstSlice<T,Ranges...>(rtp, data);
+    }
+
+    template <typename C, typename T, class... Ranges>
+    auto tcast(MultiArray<T,Ranges...>& ma)
+	-> decltype(rangeTpToSlice
+	    ( rtcast<C,T>( ma.range()->space() ),
+	      reinterpret_cast<C*>( ma.data() ) ))
+    {   
+	// VCHECK(reinterpret_cast<std::intptr_t>(ma.data()) % 32);
+	//VCHECK(reinterpret_cast<std::intptr_t>(reinterpret_cast<C*>(ma.data())) % 32);
+	return rangeTpToSlice
+	    ( rtcast<C,T>( ma.range()->space() ),
+	      reinterpret_cast<C*>( ma.data() ) );
+    }
+    
+    template <typename C, typename T, class... Ranges>
+    auto tcast(const MultiArray<T,Ranges...>& ma)
+	-> decltype(rangeTpToSlice
+	    ( rtcast<C,T>( ma.range()->space() ),
+	      reinterpret_cast<const C*>( ma.data() ) ))
+    {
+	//VCHECK(reinterpret_cast<std::intptr_t>(ma.data()) % 32);
+	//VCHECK(reinterpret_cast<std::intptr_t>(reinterpret_cast<C*>(ma.data())) % 32);
+	return rangeTpToSlice
+	    ( rtcast<C,T>( ma.range()->space() ),
+	      reinterpret_cast<const C*>( ma.data() ) );
+    }
 }
-
 #endif
diff --git a/src/include/functional_multi_array.h b/src/include/functional_multi_array.h
index 7f4c50b..2e10629 100644
--- a/src/include/functional_multi_array.h
+++ b/src/include/functional_multi_array.h
@@ -15,7 +15,7 @@ namespace MultiArrayTools
 	static auto mk(const std::shared_ptr<Index>& i)
 	    -> MultiArray<typename Index::MetaType, typename Index::RangeType>
 	{
-	    std::vector<typename Index::MetaType> vv(i->range()->size());
+	    vector<typename Index::MetaType> vv(i->range()->size());
 	    for(Index j = (*i); j.pos() != j.max(); ++j){
 		vv[j.pos()] = j.meta();
 	    }	    
diff --git a/src/include/helper_tools.cc.h b/src/include/helper_tools.cc.h
index fc09571..28750a3 100644
--- a/src/include/helper_tools.cc.h
+++ b/src/include/helper_tools.cc.h
@@ -109,7 +109,7 @@ namespace MultiArrayTools
     }
 
     template <class Range>
-    auto createRange(const std::vector<char>& cvec)
+    auto createRange(const vector<char>& cvec)
 	-> std::shared_ptr<Range>
     {
 	const char* dp = cvec.data();
diff --git a/src/include/helper_tools.h b/src/include/helper_tools.h
index 1eaa9a6..aa35dd1 100644
--- a/src/include/helper_tools.h
+++ b/src/include/helper_tools.h
@@ -61,7 +61,7 @@ namespace MultiArrayTools
 	-> std::shared_ptr<typename RangeFactory::oType>;
 
     template <class Range>
-    auto createRange(const std::vector<char>& cvec)
+    auto createRange(const vector<char>& cvec)
 	-> std::shared_ptr<Range>;
 
     template <size_t N, class MArray>
diff --git a/src/include/map_range.cc.h b/src/include/map_range.cc.h
index c4841d3..2f661a6 100644
--- a/src/include/map_range.cc.h
+++ b/src/include/map_range.cc.h
@@ -347,7 +347,7 @@ namespace MultiArrayTools
 	    }
 	}
 	if(not check){
-	    std::vector<std::intptr_t> pv(sizeof...(Ranges));
+	    vector<std::intptr_t> pv(sizeof...(Ranges));
 	    RPackNum<sizeof...(Ranges)-1>::RangesToVec(ptp, pv);
 	    pv.push_back( reinterpret_cast<std::intptr_t>
 			  ( &std::dynamic_pointer_cast<oType>( mProd )->mMapf ) );
@@ -376,8 +376,8 @@ namespace MultiArrayTools
                 mult[mapf[ii]]++;
             }
 	
-            std::vector<typename MapF::value_type> outmeta(mult.size());
-            std::vector<size_t> outmult(mult.size());
+            vector<typename MapF::value_type> outmeta(mult.size());
+            vector<size_t> outmult(mult.size());
 
             size_t cnt = 0;
             for(auto& x: mult){
@@ -404,12 +404,12 @@ namespace MultiArrayTools
             for(auto ii = mapf.begin(); ii.max() != ii.pos(); ++ii) {
                 max = mapf[ii]+1 > max ? mapf[ii]+1 : max;
             }
-            std::vector<size_t> mult(max,0);
+            vector<size_t> mult(max,0);
             for(auto ii = mapf.begin(); ii.max() != ii.pos(); ++ii) {
                 mult[mapf[ii]]++;
             }
 	
-            std::vector<size_t> outmult(mult.size());
+            vector<size_t> outmult(mult.size());
 
             size_t cnt = 0;
             for(auto& x: mult){
@@ -505,10 +505,10 @@ namespace MultiArrayTools
     }
 
     template <class MapF, SpaceType XSTYPE, class... Ranges>
-    std::vector<char> GenMapRange<MapF,XSTYPE,Ranges...>::data() const
+    vector<char> GenMapRange<MapF,XSTYPE,Ranges...>::data() const
     {
 	DataHeader h = dataHeader();
-	std::vector<char> out;
+	vector<char> out;
 	//out.reserve(h.metaSize + sizeof(DataHeader));
 	char* hcp = reinterpret_cast<char*>(&h);
 	out.insert(out.end(), hcp, hcp + sizeof(DataHeader));
diff --git a/src/include/map_range.h b/src/include/map_range.h
index 61c7d43..0406c7c 100644
--- a/src/include/map_range.h
+++ b/src/include/map_range.h
@@ -259,7 +259,7 @@ namespace MultiArrayTools
         virtual DataHeader dataHeader() const final;
         
 	virtual std::string stringMeta(size_t pos) const final;
-	virtual std::vector<char> data() const final;
+	virtual vector<char> data() const final;
 	
 	const Space& space() const;
 	
diff --git a/src/include/map_range_factory_product_map.h b/src/include/map_range_factory_product_map.h
index 992a822..d02d7c0 100644
--- a/src/include/map_range_factory_product_map.h
+++ b/src/include/map_range_factory_product_map.h
@@ -18,7 +18,7 @@ namespace MultiArrayTools
 	friend class GenMapRangeFactory;
 	
     private:
-	static std::map<std::shared_ptr<RangeBase>,std::vector<std::intptr_t> > mAleadyCreated;
+	static std::map<std::shared_ptr<RangeBase>,vector<std::intptr_t> > mAleadyCreated;
     };
     
 }
diff --git a/src/include/mbase_def.h b/src/include/mbase_def.h
index 72a8b7f..2c53cb2 100644
--- a/src/include/mbase_def.h
+++ b/src/include/mbase_def.h
@@ -3,6 +3,7 @@
 #define __mbase_def_h__
 
 #include "ranges/rbase_def.h"
+#include "allocator.h"
 
 namespace MultiArrayTools
 {
diff --git a/src/include/multi_array.cc.h b/src/include/multi_array.cc.h
index 55af2c7..1395490 100644
--- a/src/include/multi_array.cc.h
+++ b/src/include/multi_array.cc.h
@@ -7,7 +7,7 @@ namespace MultiArrayTools
     Scalar<T> scalar(const T& in)
     {
 	NullRF nrf;
-	return Scalar<T>( std::dynamic_pointer_cast<NullRange>( nrf.create() ), std::vector<T>( { in } ) );
+	return Scalar<T>( std::dynamic_pointer_cast<NullRange>( nrf.create() ), vector<T>( { in } ) );
     }
     
     /*******************
@@ -24,7 +24,7 @@ namespace MultiArrayTools
 
     template <typename T, class... SRanges>
     MultiArray<T,SRanges...>::MultiArray(const typename CRange::Space& space,
-					 std::vector<T>&& vec) :
+					 vector<T>&& vec) :
 	MutableMultiArrayBase<T,SRanges...>(space),
 	mCont(vec)
     {
@@ -52,7 +52,7 @@ namespace MultiArrayTools
     }
     
     template <typename T, class... SRanges>
-    MultiArray<T,SRanges...>::MultiArray(const std::shared_ptr<SRanges>&... ranges, const std::vector<T>& vec) :
+    MultiArray<T,SRanges...>::MultiArray(const std::shared_ptr<SRanges>&... ranges, const vector<T>& vec) :
 	MutableMultiArrayBase<T,SRanges...>(ranges...),
 	mCont(vec)
     {
@@ -63,7 +63,7 @@ namespace MultiArrayTools
     }
     
     template <typename T, class... SRanges>
-    MultiArray<T,SRanges...>::MultiArray(const std::shared_ptr<SRanges>&... ranges, std::vector<T>&& vec) :
+    MultiArray<T,SRanges...>::MultiArray(const std::shared_ptr<SRanges>&... ranges, vector<T>&& vec) :
 	MutableMultiArrayBase<T,SRanges...>(ranges...),
 	mCont(vec)
     {
diff --git a/src/include/multi_array.h b/src/include/multi_array.h
index 18ec445..aaa54ff 100644
--- a/src/include/multi_array.h
+++ b/src/include/multi_array.h
@@ -24,7 +24,6 @@ namespace MultiArrayTools
 	    return ma;
 	}
     };
-
     
     template <typename T, class... SRanges>
     class MultiArray : public MutableMultiArrayBase<T,SRanges...>
@@ -41,15 +40,15 @@ namespace MultiArrayTools
 	DEFAULT_MEMBERS(MultiArray);
 	MultiArray(const std::shared_ptr<SRanges>&... ranges);
 	MultiArray(const std::shared_ptr<SRanges>&... ranges, const T& val);
-	MultiArray(const std::shared_ptr<SRanges>&... ranges, const std::vector<T>& vec);
-	MultiArray(const std::shared_ptr<SRanges>&... ranges, std::vector<T>&& vec);
+	MultiArray(const std::shared_ptr<SRanges>&... ranges, const vector<T>& vec);
+	MultiArray(const std::shared_ptr<SRanges>&... ranges, vector<T>&& vec);
 	MultiArray(const typename CRange::Space& space);
-	MultiArray(const typename CRange::Space& space, std::vector<T>&& vec);
+	MultiArray(const typename CRange::Space& space, vector<T>&& vec);
 	MultiArray(MultiArray<T,AnonymousRange>&& ama, SIZET<SRanges>... sizes);
 	
 	// Only if ALL ranges have default extensions:
-	//MultiArray(const std::vector<T>& vec);
-	//MultiArray(std::vector<T>&& vec);
+	//MultiArray(const vector<T>& vec);
+	//MultiArray(vector<T>&& vec);
 	
 	// template <class Range2, class Range3>
 	// MultiArray(const MultiArray<MultiArray<T,Range2>,Range3> in);
@@ -76,8 +75,8 @@ namespace MultiArrayTools
 	
 	virtual const T* data() const override;
 	virtual T* data() override;
-	virtual std::vector<T>& vdata() { return mCont; }
-	virtual const std::vector<T>& vdata() const { return mCont; }
+	virtual vector<T>& vdata() { return mCont; }
+	virtual const vector<T>& vdata() const { return mCont; }
 
 	virtual std::shared_ptr<MultiArrayBase<T,AnonymousRange> > anonymous(bool slice = false) const override;
 	virtual std::shared_ptr<MultiArrayBase<T,AnonymousRange> > anonymousMove() override;
@@ -99,7 +98,7 @@ namespace MultiArrayTools
 	
     private:
 	
-	std::vector<T> mCont;
+	vector<T> mCont;
     };
 
     template <typename T>
@@ -119,7 +118,7 @@ namespace MultiArrayTools
 	    const size_t smas = sma.size();
 	    const size_t mas = ma.size();
 	    auto cr = ma.range()->cat(sma.range());
-	    std::vector<T> ov;
+	    vector<T> ov;
 	    ov.reserve(mas * smas);
 
 	    for(auto& x: ma){
diff --git a/src/include/ranges/anonymous_range.h b/src/include/ranges/anonymous_range.h
index 0099258..78ad36c 100644
--- a/src/include/ranges/anonymous_range.h
+++ b/src/include/ranges/anonymous_range.h
@@ -36,7 +36,7 @@ namespace MultiArrayTools
 	template <class... RangeTypes>
 	AnonymousRangeFactory(std::shared_ptr<RangeTypes>... origs);
 
-	AnonymousRangeFactory(const std::vector<std::shared_ptr<RangeBase>>& origs);
+	AnonymousRangeFactory(const vector<std::shared_ptr<RangeBase>>& origs);
 	
 	template <class Range>
 	void append(std::shared_ptr<Range> r);
@@ -45,9 +45,9 @@ namespace MultiArrayTools
 
     private:
 	
-	std::shared_ptr<RangeBase> checkIfCreated(const std::vector<std::shared_ptr<RangeBase> >& pvec);
+	std::shared_ptr<RangeBase> checkIfCreated(const vector<std::shared_ptr<RangeBase> >& pvec);
 	
-	static std::map<std::shared_ptr<RangeBase>,std::vector<std::intptr_t> > mAleadyCreated;
+	static std::map<std::shared_ptr<RangeBase>,vector<std::intptr_t> > mAleadyCreated;
 
 	bool mProductCreated = false;
     };
@@ -81,7 +81,7 @@ namespace MultiArrayTools
         virtual DataHeader dataHeader() const final;
         
 	virtual std::string stringMeta(size_t pos) const final;
-	virtual std::vector<char> data() const final;
+	virtual vector<char> data() const final;
 	
 	std::shared_ptr<RangeBase> sub(size_t num) const;
 
@@ -91,10 +91,10 @@ namespace MultiArrayTools
 	template <class... Ranges>
 	std::shared_ptr<MultiRange<Ranges...> > scast(SIZET<Ranges>... sizes) const; // save cast
 
-        const std::vector<std::shared_ptr<RangeBase> >& orig() const;
+        const vector<std::shared_ptr<RangeBase> >& orig() const;
         
         std::shared_ptr<AnonymousRange> sreplace(const std::shared_ptr<RangeBase> in, size_t num) const;
-        std::shared_ptr<AnonymousRange> sreplace(const std::vector<std::shared_ptr<RangeBase>>& in, size_t num) const;
+        std::shared_ptr<AnonymousRange> sreplace(const vector<std::shared_ptr<RangeBase>>& in, size_t num) const;
 	
 	bool isEmpty() const;
 	
@@ -114,12 +114,12 @@ namespace MultiArrayTools
 	template <class... RangeTypes>
 	GenSingleRange(std::shared_ptr<RangeTypes>... origs);
 
-	GenSingleRange(const std::vector<std::shared_ptr<RangeBase>>& origs);
+	GenSingleRange(const vector<std::shared_ptr<RangeBase>>& origs);
 	
 	size_t mSize = 1;
 	bool mEmpty = true;
 	
-	std::vector<std::shared_ptr<RangeBase> > mOrig;
+	vector<std::shared_ptr<RangeBase> > mOrig;
     };
 }
 
@@ -172,7 +172,7 @@ namespace MultiArrayHelper
 
     template <>
     inline void resolveSetRange<AnonymousRange>(std::shared_ptr<AnonymousRange>& rp,
-						const std::vector<std::shared_ptr<RangeBase> >& orig,
+						const vector<std::shared_ptr<RangeBase> >& orig,
 						size_t origpos, size_t size)
     {
     	AnonymousRangeFactory arf;
@@ -184,7 +184,7 @@ namespace MultiArrayHelper
     }
 
     template <>
-    inline void setRangeToVec<AnonymousRange>(std::vector<std::shared_ptr<RangeBase> >& v,
+    inline void setRangeToVec<AnonymousRange>(vector<std::shared_ptr<RangeBase> >& v,
 					      std::shared_ptr<AnonymousRange> r)
     {
 	if(not r->isEmpty()){
diff --git a/src/include/ranges/dynamic_meta.h b/src/include/ranges/dynamic_meta.h
index 5f054e8..29346f7 100644
--- a/src/include/ranges/dynamic_meta.h
+++ b/src/include/ranges/dynamic_meta.h
@@ -5,6 +5,8 @@
 #include <cstdlib>
 #include <utility>
 
+#include "allocator.h"
+
 namespace MultiArrayTools
 {
 
@@ -14,7 +16,7 @@ namespace MultiArrayTools
     struct DynamicMetaSetter
     {
 	template <typename... Us>
-	static inline void set(const std::tuple<Us...>& in, std::vector<DynamicMetaElem>& meta)
+	static inline void set(const std::tuple<Us...>& in, vector<DynamicMetaElem>& meta)
 	{
 	    typedef typename std::tuple_element<N,std::tuple<Us...>>::type elemtype;
 	    meta[N].first = reinterpret_cast<const char*>( &std::get<N>(in) );
@@ -27,7 +29,7 @@ namespace MultiArrayTools
     struct DynamicMetaSetter<0>
     {
 	template <typename... Us>
-	static inline void set(const std::tuple<Us...>& in, std::vector<DynamicMetaElem>& meta)
+	static inline void set(const std::tuple<Us...>& in, vector<DynamicMetaElem>& meta)
 	{
 	    typedef typename std::tuple_element<0,std::tuple<Us...>>::type elemtype;
 	    meta[0].first = reinterpret_cast<const char*>( &std::get<0>(in) );
@@ -39,7 +41,7 @@ namespace MultiArrayTools
     class DynamicMetaT
     {
     private:
-        std::vector<DynamicMetaElem> mMeta;
+        vector<DynamicMetaElem> mMeta;
 
     public:
         DynamicMetaT() = default;
diff --git a/src/include/ranges/dynamic_range.cc.h b/src/include/ranges/dynamic_range.cc.h
index 254b5c5..69a7980 100644
--- a/src/include/ranges/dynamic_range.cc.h
+++ b/src/include/ranges/dynamic_range.cc.h
@@ -72,7 +72,7 @@ namespace MultiArrayTools
     }
 
     template <class EC>
-    DynamicRangeFactory<EC>::DynamicRangeFactory(const std::vector<std::shared_ptr<RangeBase>>& origs)
+    DynamicRangeFactory<EC>::DynamicRangeFactory(const vector<std::shared_ptr<RangeBase>>& origs)
     {
         mProd = std::shared_ptr<oType>( new DynamicRange<EC>( origs ) );
     }
@@ -99,10 +99,10 @@ namespace MultiArrayTools
 
     // INSTANCIATE IF NEEDED!!
     template <class EC>
-    std::map<std::shared_ptr<RangeBase>,std::vector<std::intptr_t> > DynamicRangeFactory<EC>::mAleadyCreated;
+    std::map<std::shared_ptr<RangeBase>,vector<std::intptr_t> > DynamicRangeFactory<EC>::mAleadyCreated;
 
     template <class EC>
-    std::shared_ptr<RangeBase> DynamicRangeFactory<EC>::checkIfCreated(const std::vector<std::shared_ptr<RangeBase> >& pvec)
+    std::shared_ptr<RangeBase> DynamicRangeFactory<EC>::checkIfCreated(const vector<std::shared_ptr<RangeBase> >& pvec)
     {
 	std::shared_ptr<RangeBase> out;
 	bool check = false;
@@ -122,7 +122,7 @@ namespace MultiArrayTools
 	    }
 	}
 	if(not check){
-	    std::vector<std::intptr_t> app(pvec.size());
+	    vector<std::intptr_t> app(pvec.size());
 	    for(size_t i = 0; i != app.size(); ++i){
 		app[i] = reinterpret_cast<std::intptr_t>( pvec[i].get() );
 	    }
@@ -223,7 +223,7 @@ namespace MultiArrayTools
     }
 
     template <class EC>
-    DynamicIndex<EC>& DynamicIndex<EC>::operator()(const std::vector<std::shared_ptr<IndexW<EC>>>& ivec)
+    DynamicIndex<EC>& DynamicIndex<EC>::operator()(const vector<std::shared_ptr<IndexW<EC>>>& ivec)
     {
         assert(mIVec.size() == ivec.size());
 	for(size_t i = 0; i != mIVec.size(); ++i){
@@ -237,7 +237,7 @@ namespace MultiArrayTools
     template <class... Indices>
     DynamicIndex<EC>& DynamicIndex<EC>::operator()(const std::shared_ptr<Indices>&... is)
     {
-	std::vector<std::shared_ptr<IndexW<EC>>> tmp =
+	vector<std::shared_ptr<IndexW<EC>>> tmp =
 	    { std::make_shared<IndexWrapper<Indices,EC>>(is)... };
 	
 	assert(mIVec.size() == tmp.size());
@@ -489,10 +489,10 @@ namespace MultiArrayTools
     }
 
     template <class EC>
-    std::vector<char> DynamicRange<EC>::data() const
+    vector<char> DynamicRange<EC>::data() const
     {
 	DataHeader h = dataHeader();
-	std::vector<char> out;
+	vector<char> out;
 	char* hcp = reinterpret_cast<char*>(&h);
 	out.insert(out.end(), hcp, hcp + sizeof(DataHeader));
 	for(auto& x: mOrig){
@@ -571,7 +571,7 @@ namespace MultiArrayTools
     }
 
     template <class EC>
-    DynamicRange<EC>::DynamicRange(const std::vector<std::shared_ptr<RangeBase>>& origs) :
+    DynamicRange<EC>::DynamicRange(const vector<std::shared_ptr<RangeBase>>& origs) :
         RangeInterface<DynamicIndex<EC>>(),
         mOrig(origs)
     {
@@ -602,7 +602,7 @@ namespace MultiArrayTools
     }
 
     template <class EC>
-    const std::vector<std::shared_ptr<RangeBase> >& DynamicRange<EC>::orig() const
+    const vector<std::shared_ptr<RangeBase> >& DynamicRange<EC>::orig() const
     {
         return mOrig;
     }
diff --git a/src/include/ranges/dynamic_range.h b/src/include/ranges/dynamic_range.h
index 4388eb0..61e3525 100644
--- a/src/include/ranges/dynamic_range.h
+++ b/src/include/ranges/dynamic_range.h
@@ -172,7 +172,7 @@ namespace MultiArrayTools
     class DynamicIndex : public IndexInterface<DynamicIndex<EC>,DynamicMetaT>
     {
     private:
-        typedef std::vector<std::pair<std::shared_ptr<IndexW<EC>>,size_t>> IVecT;
+        typedef vector<std::pair<std::shared_ptr<IndexW<EC>>,size_t>> IVecT;
         
         IVecT mIVec;
 
@@ -197,7 +197,7 @@ namespace MultiArrayTools
 	DynamicIndex& operator--();
 
         DynamicIndex& operator()(const IVecT& ivec);
-        DynamicIndex& operator()(const std::vector<std::shared_ptr<IndexW<EC>>>& ivec);
+        DynamicIndex& operator()(const vector<std::shared_ptr<IndexW<EC>>>& ivec);
 		
 	template <class... Indices>
 	DynamicIndex& operator()(const std::shared_ptr<Indices>&... is);
@@ -259,7 +259,7 @@ namespace MultiArrayTools
 	template <class... RangeTypes>
 	DynamicRangeFactory(std::shared_ptr<RangeTypes>... origs);
 
-	DynamicRangeFactory(const std::vector<std::shared_ptr<RangeBase>>& origs);
+	DynamicRangeFactory(const vector<std::shared_ptr<RangeBase>>& origs);
 
 	template <class Range>
 	void append(std::shared_ptr<Range> r);
@@ -268,9 +268,9 @@ namespace MultiArrayTools
 
     private:
 	
-	std::shared_ptr<RangeBase> checkIfCreated(const std::vector<std::shared_ptr<RangeBase> >& pvec);
+	std::shared_ptr<RangeBase> checkIfCreated(const vector<std::shared_ptr<RangeBase> >& pvec);
 	
-	static std::map<std::shared_ptr<RangeBase>,std::vector<std::intptr_t> > mAleadyCreated;
+	static std::map<std::shared_ptr<RangeBase>,vector<std::intptr_t> > mAleadyCreated;
 
 	bool mProductCreated = false;
     };
@@ -299,12 +299,12 @@ namespace MultiArrayTools
 	template <class... RangeTypes>
 	DynamicRange(std::shared_ptr<RangeTypes>... origs);
 
-        DynamicRange(const std::vector<std::shared_ptr<RangeBase>>& origs);
+        DynamicRange(const vector<std::shared_ptr<RangeBase>>& origs);
 
 	size_t mSize = 1;
 	bool mEmpty = true;
 	
-	std::vector<std::shared_ptr<RangeBase> > mOrig;
+	vector<std::shared_ptr<RangeBase> > mOrig;
 
     public:
 	virtual size_t size() const final;
@@ -320,7 +320,7 @@ namespace MultiArrayTools
         virtual DataHeader dataHeader() const final;
         
 	virtual std::string stringMeta(size_t pos) const final;
-	virtual std::vector<char> data() const final;
+	virtual vector<char> data() const final;
 	
 	std::shared_ptr<RangeBase> sub(size_t num) const;
 
@@ -330,7 +330,7 @@ namespace MultiArrayTools
 	template <class... Ranges>
 	std::shared_ptr<MultiRange<Ranges...> > scast(SIZET<Ranges>... sizes) const; // save cast
 
-        const std::vector<std::shared_ptr<RangeBase> >& orig() const;
+        const vector<std::shared_ptr<RangeBase> >& orig() const;
         
 	void sreplace(const std::shared_ptr<RangeBase> in, size_t num);
 	
@@ -352,7 +352,7 @@ namespace MultiArrayHelper
 
     template <class EC>
     inline void resolveSetRange(std::shared_ptr<DynamicRange<EC>>& rp,
-                                const std::vector<std::shared_ptr<RangeBase> >& orig,
+                                const vector<std::shared_ptr<RangeBase> >& orig,
                                 size_t origpos, size_t size)
     {
     	DynamicRangeFactory<EC> arf;
@@ -364,7 +364,7 @@ namespace MultiArrayHelper
     }
 
     template <class EC>
-    inline void setRangeToVec(std::vector<std::shared_ptr<RangeBase> >& v,
+    inline void setRangeToVec(vector<std::shared_ptr<RangeBase> >& v,
                               std::shared_ptr<DynamicRange<EC>> r)
     {
 	if(not r->isEmpty()){
diff --git a/src/include/ranges/multi_range.h b/src/include/ranges/multi_range.h
index 9f6ae4e..dcf0eb5 100644
--- a/src/include/ranges/multi_range.h
+++ b/src/include/ranges/multi_range.h
@@ -201,7 +201,7 @@ namespace MultiArrayTools
         virtual DataHeader dataHeader() const final;
         
 	virtual std::string stringMeta(size_t pos) const final;
-	virtual std::vector<char> data() const final;
+	virtual vector<char> data() const final;
 	
 	const Space& space() const;
 	
@@ -519,7 +519,7 @@ namespace MultiArrayTools
 	    }
 	}
 	if(not check){
-	    std::vector<std::intptr_t> pv(sizeof...(Ranges));
+	    vector<std::intptr_t> pv(sizeof...(Ranges));
 	    RPackNum<sizeof...(Ranges)-1>::RangesToVec(ptp, pv);
 	    MultiRangeFactoryProductMap::mAleadyCreated[mProd] = pv;
 	    out = mProd;
@@ -590,10 +590,10 @@ namespace MultiArrayTools
     }
 
     template <class... Ranges>
-    std::vector<char> MultiRange<Ranges...>::data() const
+    vector<char> MultiRange<Ranges...>::data() const
     {
 	DataHeader h = dataHeader();
-	std::vector<char> out;
+	vector<char> out;
 	//out.reserve(h.metaSize + sizeof(DataHeader));
 	char* hcp = reinterpret_cast<char*>(&h);
 	out.insert(out.end(), hcp, hcp + sizeof(DataHeader));
diff --git a/src/include/ranges/multi_range_factory_product_map.h b/src/include/ranges/multi_range_factory_product_map.h
index c6b10df..e9dcb38 100644
--- a/src/include/ranges/multi_range_factory_product_map.h
+++ b/src/include/ranges/multi_range_factory_product_map.h
@@ -17,7 +17,7 @@ namespace MultiArrayTools
 	friend class MultiRangeFactory;
 	
     private:
-	static std::map<std::shared_ptr<RangeBase>,std::vector<std::intptr_t> > mAleadyCreated;
+	static std::map<std::shared_ptr<RangeBase>,vector<std::intptr_t> > mAleadyCreated;
     };
     
 }
diff --git a/src/include/ranges/range_base.h b/src/include/ranges/range_base.h
index fb5605b..846070e 100644
--- a/src/include/ranges/range_base.h
+++ b/src/include/ranges/range_base.h
@@ -81,7 +81,7 @@ namespace MultiArrayTools
         std::intptr_t id() const;
         
 	virtual std::string stringMeta(size_t pos) const = 0;
-	virtual std::vector<char> data() const = 0; // usefull when writing to files, etc...
+	virtual vector<char> data() const = 0; // usefull when writing to files, etc...
 
 	virtual SpaceType spaceType() const = 0;
         virtual DataHeader dataHeader() const = 0;
diff --git a/src/include/ranges/range_types/classic_range.h b/src/include/ranges/range_types/classic_range.h
index f9b4ed6..ee4ae3c 100644
--- a/src/include/ranges/range_types/classic_range.h
+++ b/src/include/ranges/range_types/classic_range.h
@@ -43,7 +43,7 @@ namespace MultiArrayTools
         virtual DataHeader dataHeader() const final;
         
 	virtual std::string stringMeta(size_t pos) const final;
-	virtual std::vector<char> data() const final;
+	virtual vector<char> data() const final;
 	
 	size_t get(size_t pos) const;
 	size_t getMeta(size_t metaPos) const;
diff --git a/src/include/ranges/range_types/null_range.h b/src/include/ranges/range_types/null_range.h
index 994bb29..a7db938 100644
--- a/src/include/ranges/range_types/null_range.h
+++ b/src/include/ranges/range_types/null_range.h
@@ -47,7 +47,7 @@ namespace MultiArrayTools
 	virtual size_t dim() const final;
 
 	virtual std::string stringMeta(size_t pos) const final;
-	virtual std::vector<char> data() const final;
+	virtual vector<char> data() const final;
 
 	virtual SpaceType spaceType() const final;
         virtual DataHeader dataHeader() const final;
diff --git a/src/include/ranges/range_types/space_range.h b/src/include/ranges/range_types/space_range.h
index 46d660a..8235ff4 100644
--- a/src/include/ranges/range_types/space_range.h
+++ b/src/include/ranges/range_types/space_range.h
@@ -39,7 +39,7 @@ namespace MultiArrayTools
 	virtual size_t dim() const final;
 
 	virtual std::string stringMeta(size_t pos) const final;
-	virtual std::vector<char> data() const final;
+	virtual vector<char> data() const final;
 
 	virtual SpaceType spaceType() const final;
         virtual DataHeader dataHeader() const final;
diff --git a/src/include/ranges/range_types/spin_range.h b/src/include/ranges/range_types/spin_range.h
index 97992b5..d720f94 100644
--- a/src/include/ranges/range_types/spin_range.h
+++ b/src/include/ranges/range_types/spin_range.h
@@ -40,7 +40,7 @@ namespace MultiArrayTools
 	virtual size_t dim() const final;
 
 	virtual std::string stringMeta(size_t pos) const final;
-	virtual std::vector<char> data() const final;
+	virtual vector<char> data() const final;
 
 	virtual SpaceType spaceType() const final;
         virtual DataHeader dataHeader() const final;
diff --git a/src/include/ranges/rbase_def.h b/src/include/ranges/rbase_def.h
index 64ccc45..41c0a28 100644
--- a/src/include/ranges/rbase_def.h
+++ b/src/include/ranges/rbase_def.h
@@ -14,7 +14,7 @@
   << " in " << __func__ << ": " << #a << " = " << a << std::endl;
 #endif
 
-
+#include "allocator.h"
 
 namespace MultiArrayTools
 {
diff --git a/src/include/ranges/rpack_num.h b/src/include/ranges/rpack_num.h
index a583455..7cabbe0 100644
--- a/src/include/ranges/rpack_num.h
+++ b/src/include/ranges/rpack_num.h
@@ -46,7 +46,7 @@ namespace MultiArrayHelper
 
     
     template <class Range>
-    inline void resolveSetRange(std::shared_ptr<Range>& rp, const std::vector<std::shared_ptr<RangeBase> >& orig,
+    inline void resolveSetRange(std::shared_ptr<Range>& rp, const vector<std::shared_ptr<RangeBase> >& orig,
 				size_t origpos, size_t size)
     {
 	assert(size == 1);
@@ -54,7 +54,7 @@ namespace MultiArrayHelper
     }
 
     template <class Range>
-    inline void setRangeToVec(std::vector<std::shared_ptr<RangeBase> >& v,
+    inline void setRangeToVec(vector<std::shared_ptr<RangeBase> >& v,
 			      std::shared_ptr<Range> r)
     {
 	v.insert(v.begin(), r);
@@ -251,7 +251,7 @@ namespace MultiArrayHelper
 
 	template <class... Ranges>
 	static inline void RangesToVec(const std::tuple<std::shared_ptr<Ranges>...>& rst,
-				       std::vector<std::shared_ptr<RangeBase> >& v)
+				       vector<std::shared_ptr<RangeBase> >& v)
 	{
 	    setRangeToVec(v, std::get<N>(rst));
 	    RPackNum<N-1>::RangesToVec(rst, v);
@@ -259,7 +259,7 @@ namespace MultiArrayHelper
 
 	template <class... Ranges>
 	static inline void RangesToVec(const std::tuple<std::shared_ptr<Ranges>...>& rst,
-				       std::vector<std::intptr_t>& v)
+				       vector<std::intptr_t>& v)
 	{
 	    v[N] = reinterpret_cast<std::intptr_t>( std::get<N>(rst).get() );
 	    RPackNum<N-1>::RangesToVec(rst, v);
@@ -324,7 +324,7 @@ namespace MultiArrayHelper
 	}
 
 	template <class RangeTuple, typename... SIZET>
-	static inline void resolveRangeType(const std::vector<std::shared_ptr<RangeBase> >& orig,
+	static inline void resolveRangeType(const vector<std::shared_ptr<RangeBase> >& orig,
 					    RangeTuple& rtp, size_t off, size_t size, SIZET... sizes)
 	{
 	    constexpr size_t tps = std::tuple_size<RangeTuple>::value;
@@ -334,7 +334,7 @@ namespace MultiArrayHelper
 
 	template <class... Ranges>
 	static inline bool checkIfCreated(const std::tuple<std::shared_ptr<Ranges>...>& p,
-					  const std::vector<std::intptr_t>& a)
+					  const vector<std::intptr_t>& a)
 	{
 	    return reinterpret_cast<std::intptr_t>( std::get<N>(p).get() ) == a[N] and
 		RPackNum<N-1>::checkIfCreated(p,a);
@@ -347,22 +347,22 @@ namespace MultiArrayHelper
 	}
 
 	template <class... Ranges>
-	static inline void fillRangeDataVec(std::vector<char>& out,
+	static inline void fillRangeDataVec(vector<char>& out,
 					    const std::tuple<std::shared_ptr<Ranges>...>& tp)
 	{
-	    std::vector<char> part = std::get<sizeof...(Ranges)-N-1>(tp)->data();
+	    vector<char> part = std::get<sizeof...(Ranges)-N-1>(tp)->data();
 	    out.insert(out.end(), part.begin(), part.end());
 	    RPackNum<N-1>::fillRangeDataVec(out, tp);
 	}
 
 	template <size_t SIZE, class Range, class... Ranges>
-	static inline bool compareSpaceTypes(const std::vector<std::shared_ptr<RangeBase> >& rbvec)
+	static inline bool compareSpaceTypes(const vector<std::shared_ptr<RangeBase> >& rbvec)
 	{
 	    return rbvec[SIZE-N-1]->spaceType() == Range::STYPE and RPackNum<N-1>::template compareSpaceTypes<SIZE,Ranges...>(rbvec);
 	}
 
 	template <class... Ranges>
-	static inline void setSpace(const std::vector<std::shared_ptr<RangeBase> >& rbvec,
+	static inline void setSpace(const vector<std::shared_ptr<RangeBase> >& rbvec,
 				    std::tuple<std::shared_ptr<Ranges>...>& stp)
 	{
 	    typedef typename std::remove_reference<decltype(*std::get<N>( stp ))>::type RType;
@@ -529,14 +529,14 @@ namespace MultiArrayHelper
 
 	template <class... Ranges>
 	static inline void RangesToVec(const std::tuple<std::shared_ptr<Ranges>...>& rst,
-				       std::vector<std::intptr_t>& v)
+				       vector<std::intptr_t>& v)
 	{
 	    v[0] = reinterpret_cast<std::intptr_t>( std::get<0>(rst).get() );;
 	}
 	
 	template <class... Ranges>
 	static inline void RangesToVec(const std::tuple<std::shared_ptr<Ranges>...>& rst,
-				       std::vector<std::shared_ptr<RangeBase> >& v)
+				       vector<std::shared_ptr<RangeBase> >& v)
 	{
 	    setRangeToVec(v, std::get<0>(rst));
 	}
@@ -596,7 +596,7 @@ namespace MultiArrayHelper
 	}
 
 	template <class RangeTuple, typename... SIZET>
-	static inline void resolveRangeType(const std::vector<std::shared_ptr<RangeBase> >& orig,
+	static inline void resolveRangeType(const vector<std::shared_ptr<RangeBase> >& orig,
 					    RangeTuple& rtp, size_t off, size_t size)
 	{
 	    constexpr size_t tps = std::tuple_size<RangeTuple>::value;
@@ -605,7 +605,7 @@ namespace MultiArrayHelper
 
 	template <class... Ranges>
 	static inline bool checkIfCreated(const std::tuple<std::shared_ptr<Ranges>...>& p,
-					  const std::vector<std::intptr_t>& a)
+					  const vector<std::intptr_t>& a)
 	{
 	    return reinterpret_cast<std::intptr_t>( std::get<0>(p).get() ) == a[0];
 	}
@@ -617,21 +617,21 @@ namespace MultiArrayHelper
 	}
 
 	template <class... Ranges>
-	static inline void fillRangeDataVec(std::vector<char>& out,
+	static inline void fillRangeDataVec(vector<char>& out,
 					    const std::tuple<std::shared_ptr<Ranges>...>& tp)
 	{
-	    std::vector<char> part = std::get<sizeof...(Ranges)-1>(tp)->data();
+	    vector<char> part = std::get<sizeof...(Ranges)-1>(tp)->data();
 	    out.insert(out.end(), part.begin(), part.end());
 	}
 
 	template <size_t SIZE, class Range>
-	static inline bool compareSpaceTypes(const std::vector<std::shared_ptr<RangeBase> >& rbvec)
+	static inline bool compareSpaceTypes(const vector<std::shared_ptr<RangeBase> >& rbvec)
 	{
 	    return rbvec[SIZE-1]->spaceType() == Range::STYPE;
 	}
 
 	template <class... Ranges>
-	static inline void setSpace(const std::vector<std::shared_ptr<RangeBase> >& rbvec,
+	static inline void setSpace(const vector<std::shared_ptr<RangeBase> >& rbvec,
 				    std::tuple<std::shared_ptr<Ranges>...>& stp)
 	{
 	    typedef typename std::remove_reference<decltype(*std::get<0>( stp ))>::type RType;
diff --git a/src/include/ranges/single_range.h b/src/include/ranges/single_range.h
index 57d76f7..dd88c27 100644
--- a/src/include/ranges/single_range.h
+++ b/src/include/ranges/single_range.h
@@ -107,9 +107,8 @@ namespace MultiArrayTools
 	typedef GenSingleRange<U,TYPE,S> oType;
 
 	GenSingleRangeFactory() = delete;
-	GenSingleRangeFactory(const std::vector<U>& space);
+	GenSingleRangeFactory(const vector<U>& space);
 	std::shared_ptr<RangeBase> create();
-	
     };
 
     template <typename U>
@@ -125,7 +124,7 @@ namespace MultiArrayTools
 	MetaMap& operator=(MetaMap&& in) = default;
 
 	//MetaMap(const std::map<U,size_t>& in) : mMap(in) {}
-	MetaMap(const std::vector<U>& in)
+	MetaMap(const vector<U>& in)
 	{
 	    for(size_t i = 0; i != in.size(); ++i){
 		mMap[in[i]] = i;
@@ -140,7 +139,7 @@ namespace MultiArrayTools
     class MetaMap<std::array<int,2> >
     {
     private:
-	std::vector<size_t> mMap;
+	vector<size_t> mMap;
 	int min1;
 	int min2;
 	int max1;
@@ -156,7 +155,7 @@ namespace MultiArrayTools
 	MetaMap& operator=(const MetaMap& in) = default;
 	MetaMap& operator=(MetaMap&& in) = default;
 
-	MetaMap(const std::vector<U>& in) : min1(in[0][0]),
+	MetaMap(const vector<U>& in) : min1(in[0][0]),
 					    min2(in[0][1]),
 					    max1(in[0][0]),
 					    max2(in[0][1])
@@ -201,6 +200,20 @@ namespace MultiArrayTools
         static constexpr size_t SIZE = -1;
     };
 
+    template <SpaceType TYPE>
+    struct CheckDefault
+    {
+	static constexpr size_t ISDEFAULT = false;
+	static constexpr size_t HASMETACONT = true; 
+    };
+
+    template <>
+    struct CheckDefault<SpaceType::NONE>
+    {
+	static constexpr size_t ISDEFAULT = true;
+	static constexpr size_t HASMETACONT = false; 
+    };
+    
     template <typename U, SpaceType TYPE, size_t S>
     class GenSingleRange : public RangeInterface<GenSingleIndex<U,TYPE,S> >
     {
@@ -219,7 +232,7 @@ namespace MultiArrayTools
         virtual DataHeader dataHeader() const final;
         
 	virtual std::string stringMeta(size_t pos) const final;
-	virtual std::vector<char> data() const final;
+	virtual vector<char> data() const final;
 	
 	bool isMeta(const U& metaPos) const;
 	
@@ -231,19 +244,26 @@ namespace MultiArrayTools
 		
 	friend GenSingleRangeFactory<U,TYPE,S>;
 
-	static constexpr bool defaultable = false;
+	static GenSingleRangeFactory<U,TYPE,S> factory()
+	{
+	    static_assert( not CheckDefault<TYPE>::HASMETACONT,
+			   "asked for default factory for meta data containing range" );
+	    return GenSingleRangeFactory<U,TYPE,S>(vector<U>());
+	}
+	
+	static constexpr bool defaultable = CheckDefault<TYPE>::ISDEFAULT;
 	static constexpr size_t ISSTATIC = CheckStatic<S>::ISSTATIC;
 	static constexpr size_t SIZE = CheckStatic<S>::SIZE;
-	static constexpr bool HASMETACONT = true;
+	static constexpr bool HASMETACONT = CheckDefault<TYPE>::HASMETACONT;
 	
     protected:
 
 	GenSingleRange() = delete;
 	GenSingleRange(const GenSingleRange& in) = delete;
 	
-	GenSingleRange(const std::vector<U>& space);
+	GenSingleRange(const vector<U>& space);
 
-	std::vector<U> mSpace;
+	vector<U> mSpace;
 	//std::map<U,size_t> mMSpace;
 	MetaMap<U> mMSpace;
     };
@@ -471,7 +491,7 @@ namespace MultiArrayTools
      ********************/
 
     template <typename U, SpaceType TYPE, size_t S>
-    GenSingleRangeFactory<U,TYPE,S>::GenSingleRangeFactory(const std::vector<U>& space)
+    GenSingleRangeFactory<U,TYPE,S>::GenSingleRangeFactory(const vector<U>& space)
     {
 	mProd = std::shared_ptr<oType>( new GenSingleRange<U,TYPE,S>( space ) );
     }
@@ -488,7 +508,7 @@ namespace MultiArrayTools
      ********************/
     
     template <typename U, SpaceType TYPE, size_t S>
-    GenSingleRange<U,TYPE,S>::GenSingleRange(const std::vector<U>& space) :
+    GenSingleRange<U,TYPE,S>::GenSingleRange(const vector<U>& space) :
 	RangeInterface<GenSingleIndex<U,TYPE,S> >(),
 	mSpace(space), mMSpace(mSpace)
     {
@@ -509,10 +529,22 @@ namespace MultiArrayTools
 	return mMSpace.at(metaPos);
     }
 
+    template <size_t SIZE>
+    inline size_t getStatSizeOrDyn(size_t size)
+    {
+	return SIZE;
+    }
+
+    template <>
+    inline size_t getStatSizeOrDyn<-1>(size_t size)
+    {
+	return size;
+    }
+
     template <typename U, SpaceType TYPE, size_t S>
     size_t GenSingleRange<U,TYPE,S>::size() const
     {
-	return mSpace.size();
+	return getStatSizeOrDyn<S>(mSpace.size());
     }
 
     template <typename U, SpaceType TYPE, size_t S>
@@ -540,10 +572,10 @@ namespace MultiArrayTools
     }
 
     template <typename U, SpaceType TYPE, size_t S>
-    std::vector<char> GenSingleRange<U,TYPE,S>::data() const
+    vector<char> GenSingleRange<U,TYPE,S>::data() const
     {
         DataHeader h = dataHeader();
-	std::vector<char> out;
+	vector<char> out;
 	out.reserve(h.metaSize + sizeof(DataHeader));
 	char* hcp = reinterpret_cast<char*>(&h);
 	out.insert(out.end(), hcp, hcp + sizeof(DataHeader));
diff --git a/src/include/ranges/subrange.h b/src/include/ranges/subrange.h
index 76e5fa2..d963549 100644
--- a/src/include/ranges/subrange.h
+++ b/src/include/ranges/subrange.h
@@ -100,7 +100,7 @@ namespace MultiArrayTools
 
         SubRangeFactory() = delete;
         SubRangeFactory(const std::shared_ptr<Range>& fullRange,
-                        const std::vector<size_t>& subset);
+                        const vector<size_t>& subset);
         std::shared_ptr<RangeBase> create();
     };
     
@@ -121,7 +121,7 @@ namespace MultiArrayTools
         virtual DataHeader dataHeader() const final;
         
 	virtual std::string stringMeta(size_t pos) const final;
-	virtual std::vector<char> data() const final;
+	virtual vector<char> data() const final;
 	
 	bool isMeta(const MetaType& metaPos) const;
 	
@@ -132,7 +132,7 @@ namespace MultiArrayTools
 	virtual IndexType end() const final;
 
         std::shared_ptr<Range> fullRange() const;
-        const std::vector<size_t>& subset() const;
+        const vector<size_t>& subset() const;
         std::shared_ptr<SingleRange<MetaType,SpaceType::ANY>> outRange() const;
         
         friend SubRangeFactory<Range>;
@@ -147,10 +147,10 @@ namespace MultiArrayTools
         SubRange() = delete;
         SubRange(const SubRange& in) = delete;
 
-        SubRange(const std::shared_ptr<Range>& fullRange, const std::vector<size_t>& subset);
+        SubRange(const std::shared_ptr<Range>& fullRange, const vector<size_t>& subset);
         
         std::shared_ptr<Range> mFullRange;
-        std::vector<size_t> mSubSet;
+        vector<size_t> mSubSet;
     };
     
 } // namespace MultiArrayTools
@@ -350,7 +350,7 @@ namespace MultiArrayTools
 
     template <class Range>
     SubRangeFactory<Range>::SubRangeFactory(const std::shared_ptr<Range>& fullRange,
-                                            const std::vector<size_t>& subset)
+                                            const vector<size_t>& subset)
     {
         mProd = std::shared_ptr<oType>( new SubRange<Range>( fullRange, subset ) );
     }
@@ -368,7 +368,7 @@ namespace MultiArrayTools
 
     template <class Range>
     SubRange<Range>::SubRange(const std::shared_ptr<Range>& fullRange,
-                              const std::vector<size_t>& subset) :
+                              const vector<size_t>& subset) :
         RangeInterface<SubIndex<typename Range::IndexType>>(),
         mFullRange(fullRange), mSubSet(subset) {}
 
@@ -408,14 +408,14 @@ namespace MultiArrayTools
     }
 
     template <class Range>
-    std::vector<char> SubRange<Range>::data() const
+    vector<char> SubRange<Range>::data() const
     {
         DataHeader h = dataHeader();
-        std::vector<char> out;
+        vector<char> out;
         out.reserve(h.metaSize + sizeof(DataHeader));
         char* hcp = reinterpret_cast<char*>(&h);
         out.insert(out.end(), hcp, hcp + sizeof(DataHeader));
-        std::vector<MetaType> subvec(mSubSet.size());
+        vector<MetaType> subvec(mSubSet.size());
         size_t i = 0;
         for(auto& x: mSubSet){
             subvec[i++] = mFullRange->get(x);
@@ -477,7 +477,7 @@ namespace MultiArrayTools
     }
 
     template <class Range>
-    const std::vector<size_t>& SubRange<Range>::subset() const
+    const vector<size_t>& SubRange<Range>::subset() const
     {
         return mSubSet;
     }
@@ -485,7 +485,7 @@ namespace MultiArrayTools
     template <class Range>
     std::shared_ptr<SingleRange<typename SubRange<Range>::MetaType,SpaceType::ANY>> SubRange<Range>::outRange() const
     {
-        std::vector<MetaType> ometa(mSubSet.size());
+        vector<MetaType> ometa(mSubSet.size());
         size_t i = 0;
         for(auto& x: mSubSet){
             ometa[i++] = mFullRange->get(x);
diff --git a/src/include/ranges/type_map.h b/src/include/ranges/type_map.h
index fb3c608..1088afd 100644
--- a/src/include/ranges/type_map.h
+++ b/src/include/ranges/type_map.h
@@ -7,6 +7,8 @@
 #include <array>
 #include <iostream>
 
+#include "allocator.h"
+
 namespace MultiArrayTools
 {
 
@@ -23,7 +25,7 @@ namespace MultiArrayTools
     };
 
     template <typename T>
-    inline void stringCat(std::vector<char>& out, const std::vector<T>& in)
+    inline void stringCat(vector<char>& out, const vector<T>& in)
     {
 	//for(auto& x: in) { std::cout << x << std::endl; }
 	const char* scp = reinterpret_cast<const char*>(in.data());
@@ -31,20 +33,20 @@ namespace MultiArrayTools
     }
 
     template <typename T>
-    inline size_t metaSize(const std::vector<T>& in)
+    inline size_t metaSize(const vector<T>& in)
     {
 	return in.size() * sizeof(T);
     }
 	
     template <typename T>
-    inline void metaCat(std::vector<T>& vec, const char* begin, size_t size)
+    inline void metaCat(vector<T>& vec, const char* begin, size_t size)
     {
 	const T* tp = reinterpret_cast<const T*>( begin );
 	vec.insert(vec.end(), tp, tp + size / sizeof(T));
     }
 
     template <>
-    inline void stringCat<std::string>(std::vector<char>& out, const std::vector<std::string>& in)
+    inline void stringCat<std::string>(vector<char>& out, const vector<std::string>& in)
     {
 	//for(auto& x: in) { std::cout << x << std::endl; }
 	std::string tmp = "";
@@ -56,7 +58,7 @@ namespace MultiArrayTools
     }
 
     template <>
-    inline size_t metaSize<std::string>(const std::vector<std::string>& in)
+    inline size_t metaSize<std::string>(const vector<std::string>& in)
     {
 	size_t out = 0;
 	for(auto& x: in){
@@ -66,7 +68,7 @@ namespace MultiArrayTools
     }
     
     template <>
-    inline void metaCat<std::string>(std::vector<std::string>& vec, const char* begin, size_t size)
+    inline void metaCat<std::string>(vector<std::string>& vec, const char* begin, size_t size)
     {
 
 	std::string tmp(begin, size);
@@ -92,10 +94,10 @@ namespace MultiArrayTools
     include_type(float,4)
     include_type(double,5)
     include_type(std::string,6)
-    include_type(std::vector<size_t>,101)
-    include_type(std::vector<int>,102)
-    include_type(std::vector<double>,105)
-    include_type(std::vector<std::string>,106)
+    include_type(vector<size_t>,101)
+    include_type(vector<int>,102)
+    include_type(vector<double>,105)
+    include_type(vector<std::string>,106)
     include_type(std::array<size_t XCOMMAX() 2>,201)
     include_type(std::array<int XCOMMAX() 2>,202)
     include_type(std::array<size_t XCOMMAX() 3>,301)
diff --git a/src/include/ranges/value_range.h b/src/include/ranges/value_range.h
index 146b010..93b9a23 100644
--- a/src/include/ranges/value_range.h
+++ b/src/include/ranges/value_range.h
@@ -112,7 +112,7 @@ namespace MultiArrayTools
         virtual DataHeader dataHeader() const final;
         
 	virtual std::string stringMeta(size_t pos) const final;
-	virtual std::vector<char> data() const final;
+	virtual vector<char> data() const final;
 
         U get(size_t pos) const;
         
@@ -361,11 +361,11 @@ namespace MultiArrayTools
     }
     
     template <typename U>
-    std::vector<char> ValueRange<U>::data() const
+    vector<char> ValueRange<U>::data() const
     {
 	assert(0);
 	DataHeader h = dataHeader();
-	std::vector<char> out;
+	vector<char> out;
 	out.reserve(h.metaSize + sizeof(DataHeader));
 	char* hcp = reinterpret_cast<char*>(&h);
 	out.insert(out.end(), hcp, hcp + sizeof(DataHeader));
diff --git a/src/include/ranges/x_to_string.h b/src/include/ranges/x_to_string.h
index 8d00dd4..7dfbc6e 100644
--- a/src/include/ranges/x_to_string.h
+++ b/src/include/ranges/x_to_string.h
@@ -27,7 +27,7 @@ namespace MultiArrayHelper
     inline std::string xToString<DynamicMetaT>(const DynamicMetaT& x);
     
     template <typename T>
-    inline std::string xToString(const std::vector<T>& x);
+    inline std::string xToString(const vector<T>& x);
 
     template <typename T, size_t N>
     inline std::string xToString(const std::array<T,N>& x);
@@ -90,7 +90,7 @@ namespace MultiArrayHelper
     }
     
     template <typename T>
-    inline std::string xToString(const std::vector<T>& x)
+    inline std::string xToString(const vector<T>& x)
     {
 	std::string out = "[ ";
 	for(auto& y: x){
diff --git a/src/include/slice.cc.h b/src/include/slice.cc.h
index b777fca..919161a 100644
--- a/src/include/slice.cc.h
+++ b/src/include/slice.cc.h
@@ -115,7 +115,16 @@ namespace MultiArrayTools
     Slice<T,SRanges...>::Slice(const std::shared_ptr<SRanges>&... ranges, T* data) :
 	MutableMultiArrayBase<T,SRanges...>(ranges...),
 	mData(data) {}
-    
+
+    template <typename T, class... SRanges>
+    Slice<T,SRanges...>::Slice(const std::tuple<std::shared_ptr<SRanges>...>& ranges,
+			       T* data) :
+	MutableMultiArrayBase<T,SRanges...>(ranges),
+	mData(data)
+    {
+	MAB::mInit = true;
+    }
+
     template <typename T, class... SRanges>
     const T& Slice<T,SRanges...>::operator[](const IType& i) const
     {
diff --git a/src/include/slice.h b/src/include/slice.h
index e727a46..32daf22 100644
--- a/src/include/slice.h
+++ b/src/include/slice.h
@@ -66,6 +66,8 @@ namespace MultiArrayTools
 
 	DEFAULT_MEMBERS(Slice);
 
+	Slice(const std::tuple<std::shared_ptr<SRanges>...>& ranges,
+	      T* data = nullptr);
 	Slice(const std::shared_ptr<SRanges>&... ranges, T* data = nullptr);
 	    
 	virtual const T& operator[](const IType& i) const final;
diff --git a/src/include/type_operations.h b/src/include/type_operations.h
index efa35c6..564b83b 100644
--- a/src/include/type_operations.h
+++ b/src/include/type_operations.h
@@ -70,77 +70,77 @@ namespace MultiArrayTools
 
 	getter(size_t i) : mPos(i) {}
 
-	inline T operator()(const std::vector<T>& in)
+	inline T operator()(const vector<T>& in)
 	{
 	    return in[mPos];
 	}
     };
 
     template <typename T>
-    std::vector<T>& operator+=(std::vector<T>& a, const std::vector<T>& b)
+    vector<T>& operator+=(vector<T>& a, const vector<T>& b)
     {
 	std::transform(a.begin(), a.end(), b.begin(), a.begin(), std::plus<T>());
 	return a;
     }
 
     template <typename T>
-    std::vector<T>& operator-=(std::vector<T>& a, const std::vector<T>& b)
+    vector<T>& operator-=(vector<T>& a, const vector<T>& b)
     {
 	std::transform(a.begin(), a.end(), b.begin(), a.begin(), std::minus<T>());
 	return a;
     }
 
     template <typename T>
-    std::vector<T>& operator*=(std::vector<T>& a, const std::vector<T>& b)
+    vector<T>& operator*=(vector<T>& a, const vector<T>& b)
     {
 	std::transform(a.begin(), a.end(), b.begin(), a.begin(), std::multiplies<T>());
 	return a;
     }
 
     template <typename T>
-    std::vector<T>& operator/=(std::vector<T>& a, const std::vector<T>& b)
+    vector<T>& operator/=(vector<T>& a, const vector<T>& b)
     {
 	std::transform(a.begin(), a.end(), b.begin(), a.begin(), std::divides<T>());
 	return a;
     }
 
     template <typename T>
-    std::vector<T> operator+(std::vector<T>& a, const std::vector<T>& b)
+    vector<T> operator+(vector<T>& a, const vector<T>& b)
     {
-	std::vector<T> out(a.size());
+	vector<T> out(a.size());
 	std::transform(a.begin(), a.end(), b.begin(), out.begin(), std::plus<T>());
 	return out;
     }
 
     template <typename T>
-    std::vector<T> operator-(std::vector<T>& a, const std::vector<T>& b)
+    vector<T> operator-(vector<T>& a, const vector<T>& b)
     {
-	std::vector<T> out(a.size());
+	vector<T> out(a.size());
 	std::transform(a.begin(), a.end(), b.begin(), out.begin(), std::minus<T>());
 	return out;
     }
 
     template <typename T>
-    std::vector<T> operator*(std::vector<T>& a, const std::vector<T>& b)
+    vector<T> operator*(vector<T>& a, const vector<T>& b)
     {
-	std::vector<T> out(a.size());
+	vector<T> out(a.size());
 	std::transform(a.begin(), a.end(), b.begin(), out.begin(), std::multiplies<T>());
 	return out;
     }
 
     template <typename T>
-    std::vector<T> operator/(std::vector<T>& a, const std::vector<T>& b)
+    vector<T> operator/(vector<T>& a, const vector<T>& b)
     {
-	std::vector<T> out(a.size());
+	vector<T> out(a.size());
 	std::transform(a.begin(), a.end(), b.begin(), out.begin(), std::divides<T>());
 	return out;
     }
 
     template <class OperationClass, typename T>
-    class OperationTemplate<std::vector<T>,OperationClass> : public OperationBase<std::vector<T>,OperationClass>
+    class OperationTemplate<vector<T>,OperationClass> : public OperationBase<vector<T>,OperationClass>
     {
     public:
-	typedef OperationBase<std::vector<T>,OperationClass> OB;
+	typedef OperationBase<vector<T>,OperationClass> OB;
 
 	auto operator[](size_t i)
 	    -> Operation<T,getter<T>,OperationClass>
@@ -154,46 +154,55 @@ namespace MultiArrayTools
 	friend OperationClass;
     };
 
-    template <typename T, size_t N>
-    inline std::array<T,N> operator+(std::array<T,N>& a, const std::array<T,N>& b)
+    typedef struct v256 { double _x[4]; } v256;
+
+    template <int N>
+    inline void xadd(double* o, const double* a, const double* b)
     {
-	std::array<T,N> out;
-	for(size_t i = 0; i != N; ++i){
-	    out[i] = a[i] + b[i];
+	//#pragma omp simd aligned(o, a, b: 16)
+	for(int i = 0; i < N; i++) {
+	    o[i] = a[i] + b[i];
+	}  
+    }
+    
+    inline v256 operator+(const v256& a, const v256& b)
+    {
+	alignas(32) v256 out;
+	xadd<4>( reinterpret_cast<double*>(&out), reinterpret_cast<const double*>(&a),
+		 reinterpret_cast<const double*>(&b) );
+	return out;
+    }
+    /*
+    inline v256 operator-(const v256& a, const v256& b)
+    {
+        	alignas(32) v256 out;
+#pragma omp simd aligned(outp, ap, bp: 32)
+	for(int i = 0; i < IN; ++i){
+	    outp[i] = ap[i] - bp[i];
 	}
 	return out;
     }
 
-    template <typename T, size_t N>
-    inline std::array<T,N> operator-(std::array<T,N>& a, const std::array<T,N>& b)
+    inline v256 operator*(const v256& a, const v256& b)
     {
-	std::array<T,N> out;
-	for(size_t i = 0; i != N; ++i){
-	    out[i] = a[i] - b[i];
+	alignas(32) v256 out;
+#pragma omp simd aligned(outp, ap, bp: 32)
+	for(int i = 0; i < IN; ++i){
+	    outp[i] = ap[i] * bp[i];
 	}
 	return out;
     }
-
-    template <typename T, size_t N>
-    inline std::array<T,N> operator*(std::array<T,N>& a, const std::array<T,N>& b)
+    
+    inline v256 operator/(const v256& a, const v256& b)
     {
-	std::array<T,N> out;
-	for(size_t i = 0; i != N; ++i){
-	    out[i] = a[i] * b[i];
+	alignas(32) v256 out;
+#pragma omp simd aligned(outp, ap, bp: 32)
+	for(int i = 0; i < IN; ++i){
+	    outp[i] = ap[i] / bp[i];
 	}
 	return out;
     }
-
-    template <typename T, size_t N>
-    inline std::array<T,N> operator/(std::array<T,N>& a, const std::array<T,N>& b)
-    {
-	std::array<T,N> out;
-	for(size_t i = 0; i != N; ++i){
-	    out[i] = a[i] / b[i];
-	}
-	return out;
-    }
-
+    */
     
 
 } // namespace MultiArrayTools
diff --git a/src/include/xfor/xfor.h b/src/include/xfor/xfor.h
index c854af3..3cda525 100644
--- a/src/include/xfor/xfor.h
+++ b/src/include/xfor/xfor.h
@@ -9,6 +9,7 @@
 #include "xfor/for_utils.h"
 #include "xfor/exttype.h"
 
+#include "allocator.h"
 #include <omp.h>
 
 #define VCHECK(a) std::cout << __FILE__ << ": @" << __LINE__ \
@@ -16,7 +17,7 @@
 
 namespace MultiArrayHelper
 {
-
+    using namespace MultiArrayTools;
     // 'HIDDEN FOR' CLASS for nested for loops in contractions a.s.o.
     // (NO COUNTING OF MASTER POSITION !!!!!)
 
@@ -152,7 +153,7 @@ namespace MultiArrayHelper
 	typedef decltype(mkExt(0).extend(mExpr.rootSteps())) ExtType;
 	ExtType mExt;
 
-        const std::vector<size_t>* mSubSet;
+        const vector<size_t>* mSubSet;
         
         mutable ExtType mRootSteps;
 
@@ -169,10 +170,10 @@ namespace MultiArrayHelper
 	
 	SubExpr(const std::shared_ptr<IndexClass>& indPtr,
 		std::intptr_t siptr,
-                const std::vector<size_t>* subset, Expr expr);
+                const vector<size_t>* subset, Expr expr);
 
 	SubExpr(const IndexClass* indPtr, std::intptr_t siptr,
-                const std::vector<size_t>* subset, Expr expr);
+                const vector<size_t>* subset, Expr expr);
 
 	
 	inline void operator()(size_t mlast, DExt last) const override final;
@@ -650,7 +651,7 @@ namespace MultiArrayHelper
     template <class IndexClass, class Expr>
     SubExpr<IndexClass,Expr>::SubExpr(const std::shared_ptr<IndexClass>& indPtr,
 				      std::intptr_t siptr,
-                                      const std::vector<size_t>* subset, Expr expr) :
+                                      const vector<size_t>* subset, Expr expr) :
 	mIndPtr(indPtr.get()), mSIPtr(siptr), mSPos(mIndPtr->pos()), mMax(mIndPtr->max()),
         mExpr(expr),
 	mExt( mkExt(0).extend( mExpr.rootSteps( reinterpret_cast<std::intptr_t>( mIndPtr )) ) ),
@@ -661,7 +662,7 @@ namespace MultiArrayHelper
 
     template <class IndexClass, class Expr>
     SubExpr<IndexClass,Expr>::SubExpr(const IndexClass* indPtr, std::intptr_t siptr,
-                                      const std::vector<size_t>* subset, Expr expr) :
+                                      const vector<size_t>* subset, Expr expr) :
 	mIndPtr(indPtr), mSIPtr(siptr), mSPos(mIndPtr->pos()), mMax(mIndPtr->max()),
         mExpr(expr),
 	mExt( mkExt(0).extend( mExpr.rootSteps( reinterpret_cast<std::intptr_t>( mIndPtr )) ) ),
diff --git a/src/lib/map_range_factory_product_map.cc b/src/lib/map_range_factory_product_map.cc
index 563ec74..28fbe39 100644
--- a/src/lib/map_range_factory_product_map.cc
+++ b/src/lib/map_range_factory_product_map.cc
@@ -3,5 +3,5 @@
 
 namespace MultiArrayTools
 {
-    std::map<std::shared_ptr<RangeBase>,std::vector<std::intptr_t> > MapRangeFactoryProductMap::mAleadyCreated;
+    std::map<std::shared_ptr<RangeBase>,vector<std::intptr_t> > MapRangeFactoryProductMap::mAleadyCreated;
 }
diff --git a/src/lib/ranges/anonymous_range.cc b/src/lib/ranges/anonymous_range.cc
index 4b5a337..ae3ecb7 100644
--- a/src/lib/ranges/anonymous_range.cc
+++ b/src/lib/ranges/anonymous_range.cc
@@ -13,14 +13,14 @@ namespace MultiArrayTools
     }
 
     
-    std::map<std::shared_ptr<RangeBase>,std::vector<std::intptr_t> > AnonymousRangeFactory::mAleadyCreated;
+    std::map<std::shared_ptr<RangeBase>,vector<std::intptr_t> > AnonymousRangeFactory::mAleadyCreated;
 
-    AnonymousRangeFactory::AnonymousRangeFactory(const std::vector<std::shared_ptr<RangeBase>>& origs)
+    AnonymousRangeFactory::AnonymousRangeFactory(const vector<std::shared_ptr<RangeBase>>& origs)
     {
         mProd = std::shared_ptr<oType>( new AnonymousRange( origs ) );
     }
     
-    std::shared_ptr<RangeBase> AnonymousRangeFactory::checkIfCreated(const std::vector<std::shared_ptr<RangeBase> >& pvec)
+    std::shared_ptr<RangeBase> AnonymousRangeFactory::checkIfCreated(const vector<std::shared_ptr<RangeBase> >& pvec)
     {
 	std::shared_ptr<RangeBase> out;
 	bool check = false;
@@ -40,7 +40,7 @@ namespace MultiArrayTools
 	    }
 	}
 	if(not check){
-	    std::vector<std::intptr_t> app(pvec.size());
+	    vector<std::intptr_t> app(pvec.size());
 	    for(size_t i = 0; i != app.size(); ++i){
 		app[i] = reinterpret_cast<std::intptr_t>( pvec[i].get() );
 	    }
@@ -108,10 +108,10 @@ namespace MultiArrayTools
 	return out;
     }
 
-    std::vector<char> AnonymousRange::data() const
+    vector<char> AnonymousRange::data() const
     {
 	DataHeader h = dataHeader();
-	std::vector<char> out;
+	vector<char> out;
 	char* hcp = reinterpret_cast<char*>(&h);
 	out.insert(out.end(), hcp, hcp + sizeof(DataHeader));
 	for(auto& x: mOrig){
@@ -168,7 +168,7 @@ namespace MultiArrayTools
         return std::dynamic_pointer_cast<AnonymousRange>(arf.create());
     }
 
-    std::shared_ptr<AnonymousRange> AnonymousRange::sreplace(const std::vector<std::shared_ptr<RangeBase>>& in, size_t num) const
+    std::shared_ptr<AnonymousRange> AnonymousRange::sreplace(const vector<std::shared_ptr<RangeBase>>& in, size_t num) const
     {
         size_t nsize = 1;
         for(auto& x: in){
@@ -191,12 +191,12 @@ namespace MultiArrayTools
         return std::dynamic_pointer_cast<AnonymousRange>(arf.create());
     }
 
-    const std::vector<std::shared_ptr<RangeBase> >& AnonymousRange::orig() const
+    const vector<std::shared_ptr<RangeBase> >& AnonymousRange::orig() const
     {
         return mOrig;
     }
     
-    GenSingleRange<size_t,SpaceType::ANON,-1>::GenSingleRange(const std::vector<std::shared_ptr<RangeBase>>& origs) :
+    GenSingleRange<size_t,SpaceType::ANON,-1>::GenSingleRange(const vector<std::shared_ptr<RangeBase>>& origs) :
         RangeInterface<AnonymousIndex>(),
         mOrig(origs)
     {
diff --git a/src/lib/ranges/multi_range_factory_product_map.cc b/src/lib/ranges/multi_range_factory_product_map.cc
index 30da2dc..226ce7b 100644
--- a/src/lib/ranges/multi_range_factory_product_map.cc
+++ b/src/lib/ranges/multi_range_factory_product_map.cc
@@ -3,5 +3,5 @@
 
 namespace MultiArrayTools
 {
-    std::map<std::shared_ptr<RangeBase>,std::vector<std::intptr_t> > MultiRangeFactoryProductMap::mAleadyCreated;
+    std::map<std::shared_ptr<RangeBase>,vector<std::intptr_t> > MultiRangeFactoryProductMap::mAleadyCreated;
 }
diff --git a/src/lib/ranges/range_base.cc b/src/lib/ranges/range_base.cc
index bfe422e..6f1d5eb 100644
--- a/src/lib/ranges/range_base.cc
+++ b/src/lib/ranges/range_base.cc
@@ -15,7 +15,7 @@ namespace MultiArrayTools
     template <class... Ranges>
     using STP = std::tuple<std::shared_ptr<Ranges>...>;
 
-    typedef std::vector<std::shared_ptr<RangeBase> > RVEC;
+    typedef vector<std::shared_ptr<RangeBase> > RVEC;
 
     template <class... Ranges>
     inline bool compareSpaceTypes(const RVEC& rvec)
@@ -127,7 +127,7 @@ namespace MultiArrayTools
 		    assert(0);
 		}
 #define register_type(x) else if(x == h.metaType) {	\
-		    std::vector<TypeMap<x>::type> vd;\
+		    vector<TypeMap<x>::type> vd;\
 		    metaCat(vd, *dp, h.metaSize); \
 		    out = std::make_shared<SingleRangeFactory<TypeMap<x>::type, \
 							      SpaceType::ANY> >(vd); }
diff --git a/src/lib/ranges/range_types/classic_range.cc b/src/lib/ranges/range_types/classic_range.cc
index db147d9..7fbd227 100644
--- a/src/lib/ranges/range_types/classic_range.cc
+++ b/src/lib/ranges/range_types/classic_range.cc
@@ -58,10 +58,10 @@ namespace MultiArrayTools
 	return std::to_string(get(pos));
     }
 
-    std::vector<char> GenSingleRange<size_t,SpaceType::NONE,-1>::data() const
+    vector<char> GenSingleRange<size_t,SpaceType::NONE,-1>::data() const
     {
 	DataHeader h = dataHeader();
-	std::vector<char> out;
+	vector<char> out;
 	out.reserve(h.metaSize + sizeof(DataHeader));
 	char* hcp = reinterpret_cast<char*>(&h);
 	out.insert(out.end(), hcp, hcp + sizeof(DataHeader));
diff --git a/src/lib/ranges/range_types/null_range.cc b/src/lib/ranges/range_types/null_range.cc
index 46ac26b..107d4b0 100644
--- a/src/lib/ranges/range_types/null_range.cc
+++ b/src/lib/ranges/range_types/null_range.cc
@@ -81,10 +81,10 @@ namespace MultiArrayTools
 	return std::to_string(get(pos));
     }
 
-    std::vector<char> GenSingleRange<size_t,SpaceType::NUL,0>::data() const
+    vector<char> GenSingleRange<size_t,SpaceType::NUL,0>::data() const
     {
 	DataHeader h = dataHeader(); 
-	std::vector<char> out;
+	vector<char> out;
 	out.reserve(h.metaSize + sizeof(DataHeader));
 	char* hcp = reinterpret_cast<char*>(&h);
 	out.insert(out.end(), hcp, hcp + sizeof(DataHeader));
diff --git a/src/lib/ranges/range_types/space_range.cc b/src/lib/ranges/range_types/space_range.cc
index ae4758e..b867f2a 100644
--- a/src/lib/ranges/range_types/space_range.cc
+++ b/src/lib/ranges/range_types/space_range.cc
@@ -64,10 +64,10 @@ namespace MultiArrayTools
 	return std::to_string(get(pos));
     }
 
-    std::vector<char> GenSingleRange<int,SpaceType::PSPACE,-1>::data() const
+    vector<char> GenSingleRange<int,SpaceType::PSPACE,-1>::data() const
     {
 	DataHeader h = dataHeader();
-	std::vector<char> out;
+	vector<char> out;
 	out.reserve(h.metaSize + sizeof(DataHeader));
 	char* hcp = reinterpret_cast<char*>(&h);
 	out.insert(out.end(), hcp, hcp + sizeof(DataHeader));
diff --git a/src/lib/ranges/range_types/spin_range.cc b/src/lib/ranges/range_types/spin_range.cc
index 8ba73f4..191ae2f 100644
--- a/src/lib/ranges/range_types/spin_range.cc
+++ b/src/lib/ranges/range_types/spin_range.cc
@@ -61,10 +61,10 @@ namespace MultiArrayTools
 	return std::to_string(get(pos));
     }
 
-    std::vector<char> GenSingleRange<size_t,SpaceType::SPIN,4>::data() const
+    vector<char> GenSingleRange<size_t,SpaceType::SPIN,4>::data() const
     {
 	DataHeader h = dataHeader();
-	std::vector<char> out;
+	vector<char> out;
 	out.reserve(h.metaSize + sizeof(DataHeader));
 	char* hcp = reinterpret_cast<char*>(&h);
 	out.insert(out.end(), hcp, hcp + sizeof(DataHeader));
diff --git a/src/tests/ma_unit_test.cc b/src/tests/ma_unit_test.cc
index 705c748..20847ee 100644
--- a/src/tests/ma_unit_test.cc
+++ b/src/tests/ma_unit_test.cc
@@ -15,7 +15,7 @@ namespace {
     template <class Factory, typename T>
     void swapFactory(std::shared_ptr<RangeFactoryBase>& fptr, std::initializer_list<T> ilist)
     {
-	std::vector<T> tmp = ilist;
+	vector<T> tmp = ilist;
 	auto nptr = std::make_shared<Factory>( tmp );
 	fptr = nptr;
     }
@@ -50,7 +50,7 @@ namespace {
 
 	std::shared_ptr<RangeFactoryBase> rfbptr;
 	std::shared_ptr<SRange> srptr;
-	std::vector<double> vv = { 3.141, 2.718, 1.618, 0.693, 0.577 };
+	vector<double> vv = { 3.141, 2.718, 1.618, 0.693, 0.577 };
     };
 
     
@@ -91,7 +91,7 @@ namespace {
 	std::shared_ptr<SRange> sr3ptr;
 	std::shared_ptr<SRange> sr4ptr;
 	std::shared_ptr<MRange> mrptr;
-	std::vector<double> vv = { 2.917, 9.436, 0.373, 7.192, 7.315, 1.536, 4.892, 0.280,
+	vector<double> vv = { 2.917, 9.436, 0.373, 7.192, 7.315, 1.536, 4.892, 0.280,
 				   8.870, 4.790, 8.215, 5.063, 1.530, 3.084, 1.609, 4.847,
 				   8.175, 0.112, 6.712, 6.408, 1.959, 0.331, 4.209, 2.951 };
     };
@@ -114,7 +114,7 @@ namespace {
     
     TEST_F(MATest_1Dim, ForLoop)
     {
-	std::vector<double> v2 = { 0.693 , 2.718, 3.141, 1.618, 9.98 };
+	vector<double> v2 = { 0.693 , 2.718, 3.141, 1.618, 9.98 };
 	MultiArray<double,MATest_1Dim::SRange> ma(srptr, std::move( v2 ) );
 	size_t cnt = 0;
 	for(auto el: ma){
diff --git a/src/tests/op_perf_test.cc b/src/tests/op_perf_test.cc
index 62accee..90af492 100644
--- a/src/tests/op_perf_test.cc
+++ b/src/tests/op_perf_test.cc
@@ -5,6 +5,7 @@
 #include <cassert>
 
 #include "multi_array_header.h"
+#include "conversions.h"
 
 #include <ctime>
 #include <cmath>
@@ -32,15 +33,15 @@ namespace {
     template <class Factory, typename T>
     void swapFactory(std::shared_ptr<RangeFactoryBase>& fptr, std::initializer_list<T> ilist)
     {
-	std::vector<T> tmp = ilist;
+	vector<T> tmp = ilist;
 	auto nptr = std::make_shared<Factory>( tmp );
 	fptr = nptr;
     }
 
     template <class Factory, typename T>
-    void swapFactory(std::shared_ptr<RangeFactoryBase>& fptr, std::vector<T>& ilist)
+    void swapFactory(std::shared_ptr<RangeFactoryBase>& fptr, vector<T>& ilist)
     {
-	std::vector<T> tmp = ilist;
+	vector<T> tmp = ilist;
 	auto nptr = std::make_shared<Factory>( tmp );
 	fptr = nptr;
     }
@@ -78,14 +79,14 @@ namespace {
 	OpTest_Performance()
 	{
 	
-	    std::vector<size_t> initvec1(vs1);
-	    cv1.resize(vs1);
+	    vector<size_t> initvec1(vs1);
+ 	    cv1.resize(vs1);
 	    for(size_t i = 0; i != vs1; ++i){
 		initvec1[i] = i;
 		cv1[i] = sqrt( static_cast<double>(i)*0.53 );
 	    }
 
-	    std::vector<size_t> initvec2(vs2);
+	    vector<size_t> initvec2(vs2);
 	    cv2.resize(vs2*vs1);
 	    for(size_t i = 0; i != vs2; ++i){
 		initvec2[i] = i;
@@ -112,14 +113,14 @@ namespace {
 	//const size_t vs2 = 1000;
 	const size_t vs1 = 4000;
 	const size_t vs2 = 2500;
-
+ 
 	std::shared_ptr<RangeFactoryBase> rfbptr;
 	std::shared_ptr<SRange> sr1ptr;
 	std::shared_ptr<SRange> sr2ptr;
 	std::shared_ptr<MRange> mrptr;
 
-	std::vector<double> cv1;
-	std::vector<double> cv2;
+	vector<double> cv1;
+	vector<double> cv2;
     };
 #endif
     class OpTest_Spin
@@ -132,12 +133,12 @@ namespace {
 	typedef SpinRF SRF;
 	typedef SpinRange SR;
 	typedef MultiRangeFactory<SR,SR,SR,SR,SR,SR,SR,SR> SR8F;
-	typedef SR8F::oType SR8;
-
-        static const size_t os = 30;
-        static const size_t is = 65536;
-        static const size_t s = is*os;
-        
+ 	typedef SR8F::oType SR8;
+ 
+          static const size_t os = 300;
+          static const size_t is = 65536;
+          static const size_t s = is*os;
+	
 	OpTest_Spin()
 	{
 	    data.resize(is);
@@ -152,11 +153,11 @@ namespace {
             cr = std::dynamic_pointer_cast<CR>(cf.create());
 	}
 
-	void contract();
+     	void contract();
 	
     private:
 	
-	std::vector<double> data;
+	vector<double> data;
 	std::shared_ptr<SR> sr;
         std::shared_ptr<CR> cr;
     };
@@ -172,26 +173,27 @@ namespace {
 	auto beta = MAT::getIndex<SR>();
 	auto gamma = MAT::getIndex<SR>();
 	auto delta = MAT::getIndex<SR>();
-	auto deltap = MAT::getIndex<SR>();
+	//auto deltap = MAT::getIndex<SR>();
+	auto deltap = MAT::getIndex<GenSingleRange<size_t,SpaceType::NONE,1>>();
 	
 	auto mix = MAT::mkMIndex( alpha, beta, gamma, jj );
 
 	std::clock_t begin = std::clock();
 	//for(size_t i = 0; i != os; ++i){
         //res1(ii ,delta, deltap).par() += ma(ii, delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix);
-        res1(ii ,delta, deltap) += ma(delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix);
-            //}
+        tcast<v256>(res1)(ii ,delta, deltap) += tcast<v256>(ma)(delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix);
+	//}
 	std::clock_t end = std::clock();
 	std::cout << "MultiArray time: " << static_cast<double>( end - begin ) / CLOCKS_PER_SEC
 		  << std::endl;
 	
-	std::vector<double> vres(4*4*os);
+	vector<double> vres(4*4*os);
 	for(size_t d = 0; d != 4; ++d){
 	    for(size_t p = 0; p != 4; ++p){
 		const size_t tidx = d*4 + p;
 		vres[tidx] = 0.;
 	    }
-	}
+	}                  
 	std::clock_t begin2 = std::clock();
         for(size_t j = 0; j != os; ++j) {
             for(size_t i = 0; i != os; ++i){
@@ -232,7 +234,7 @@ namespace {
 	assert( xround(res1.at(mkts(0,3,2))) == xround(vres[14]) );
 	assert( xround(res1.at(mkts(0,3,3))) == xround(vres[15]) );
 
-	std::cout << "std::vector - for loop time: " << static_cast<double>( end2 - begin2 ) / CLOCKS_PER_SEC
+	std::cout << "vector - for loop time: " << static_cast<double>( end2 - begin2 ) / CLOCKS_PER_SEC
 		  << std::endl;
 	std::cout << "ratio: " << static_cast<double>( end - begin ) / static_cast<double>( end2 - begin2 ) << std::endl;
     }    
@@ -256,7 +258,7 @@ namespace {
 	std::cout << "MultiArray time: " << static_cast<double>( end - begin ) / CLOCKS_PER_SEC
 		  << std::endl;
 
-	std::vector<double> res2(vs1*vs2);
+	vector<double> res2(vs1*vs2);
 	std::clock_t begin2 = std::clock();
 	
 	for(size_t i = 0; i != vs2; ++i){
@@ -266,7 +268,7 @@ namespace {
 	}
 
 	std::clock_t end2 = std::clock();
-	std::cout << "std::vector - for loop time: " << static_cast<double>( end2 - begin2 ) / CLOCKS_PER_SEC
+	std::cout << "vector - for loop time: " << static_cast<double>( end2 - begin2 ) / CLOCKS_PER_SEC
 		  << std::endl;
 
 	std::cout << "ratio: " << static_cast<double>( end - begin ) / static_cast<double>( end2 - begin2 ) << std::endl;
diff --git a/src/tests/op_unit_test.cc b/src/tests/op_unit_test.cc
index 4bd8bfb..350f964 100644
--- a/src/tests/op_unit_test.cc
+++ b/src/tests/op_unit_test.cc
@@ -31,15 +31,15 @@ namespace {
     template <class Factory, typename T>
     void swapFactory(std::shared_ptr<RangeFactoryBase>& fptr, std::initializer_list<T> ilist)
     {
-	std::vector<T> tmp = ilist;
+	vector<T> tmp = ilist;
 	auto nptr = std::make_shared<Factory>( tmp );
 	fptr = nptr;
     }
 
     template <class Factory, typename T>
-    void swapFactory(std::shared_ptr<RangeFactoryBase>& fptr, std::vector<T>& ilist)
+    void swapFactory(std::shared_ptr<RangeFactoryBase>& fptr, vector<T>& ilist)
     {
-	std::vector<T> tmp = ilist;
+	vector<T> tmp = ilist;
 	auto nptr = std::make_shared<Factory>( tmp );
 	fptr = nptr;
     }
@@ -106,8 +106,8 @@ namespace {
 
 	std::shared_ptr<RangeFactoryBase> rfbptr;
 	std::shared_ptr<SRange> srptr;
-	std::vector<double> v1 = { 2.917, 9.436, 0.373, 7.192 };
-	std::vector<double> v2 = { 8.870, 4.790, 8.215, 5.063 };
+	vector<double> v1 = { 2.917, 9.436, 0.373, 7.192 };
+	vector<double> v2 = { 8.870, 4.790, 8.215, 5.063 };
     };
     
     class OpTest_MDim : public ::testing::Test
@@ -151,11 +151,11 @@ namespace {
 	std::shared_ptr<SRange> sr4ptr;
 	std::shared_ptr<MRange> mr1ptr;
 	std::shared_ptr<MRange> mr2ptr;
-	std::vector<double> v1 = { 2.917, 9.436, 0.373 };
-	std::vector<double> v2 = { 8.870, 4.790 };
-	std::vector<double> v3 = { 0.353, 4.005, 1.070, 2.310, 9.243, 2.911 };
-	std::vector<double> v4 = { 1.470, 2.210 };
-	std::vector<double> v5 = { 30.932, -33.693, -26.205, -15.504, 21.227, 17.829,
+	vector<double> v1 = { 2.917, 9.436, 0.373 };
+	vector<double> v2 = { 8.870, 4.790 };
+	vector<double> v3 = { 0.353, 4.005, 1.070, 2.310, 9.243, 2.911 };
+	vector<double> v4 = { 1.470, 2.210 };
+	vector<double> v5 = { 30.932, -33.693, -26.205, -15.504, 21.227, 17.829,
 				   -14.364, -1.868, -25.703, 13.836, 23.563, 41.339 };
     };
 
@@ -184,10 +184,10 @@ namespace {
 	std::shared_ptr<SRange> sr2ptr;
 	std::shared_ptr<SRange> sr3ptr;
 
-        std::vector<double> v1 = { 2.917, 9.436, 0.373, 0.353, 4.005, 1.070,
+        vector<double> v1 = { 2.917, 9.436, 0.373, 0.353, 4.005, 1.070,
                                    -14.364, -1.868, -25.703, 13.836, 23.563, 41.339 };
 
-        std::vector<double> v2 = { 0.353, 4.005, 1.070, 2.310, 9.243, 2.911 };
+        vector<double> v2 = { 0.353, 4.005, 1.070, 2.310, 9.243, 2.911 };
     };
 
     
@@ -222,7 +222,7 @@ namespace {
 	std::shared_ptr<SRange> sr1ptr;
 	std::shared_ptr<SRange> sr2ptr;
 	std::shared_ptr<MpRange> mpr1ptr;
-	std::vector<double> v1 = { -31.71, -77.16, -18.81,
+	vector<double> v1 = { -31.71, -77.16, -18.81,
 				   -67.06, 72.31, -54.48,
 				   -50.91, -11.62, -59.57,
 				   -42.53, 80.41, 6.35 };
@@ -241,14 +241,14 @@ namespace {
 	OpTest_Performance()
 	{
 	
-	    std::vector<size_t> initvec1(vs1);
+	    vector<size_t> initvec1(vs1);
 	    cv1.resize(vs1);
 	    for(size_t i = 0; i != vs1; ++i){
 		initvec1[i] = i;
 		cv1[i] = sqrt( static_cast<double>(i)*0.53 );
 	    }
 
-	    std::vector<size_t> initvec2(vs2);
+	    vector<size_t> initvec2(vs2);
 	    cv2.resize(vs2*vs1);
 	    for(size_t i = 0; i != vs2; ++i){
 		initvec2[i] = i;
@@ -277,8 +277,8 @@ namespace {
 	std::shared_ptr<SRange> sr2ptr;
 	std::shared_ptr<MRange> mrptr;
 
-	std::vector<double> cv1;
-	std::vector<double> cv2;
+	vector<double> cv1;
+	vector<double> cv2;
     };
 
     class OpTest_Spin : public ::testing::Test
@@ -303,7 +303,7 @@ namespace {
 	    sr = std::dynamic_pointer_cast<SR>(f.create());
 	}
 
-	std::vector<double> data;
+	vector<double> data;
 	std::shared_ptr<SR> sr;
     };
 
@@ -337,8 +337,8 @@ namespace {
 	std::shared_ptr<SR> sr1ptr;
 	std::shared_ptr<SR> sr2ptr;
 
-	std::vector<double> mv1;
-	std::vector<double> mv2;
+	vector<double> mv1;
+	vector<double> mv2;
 
 	double vcontract(size_t i)
 	{
@@ -458,7 +458,7 @@ namespace {
 
 	res2(delta, deltap).par() += ma(delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(alpha).c(beta).c(gamma);
 	
-	std::vector<double> vres(4*4);
+	vector<double> vres(4*4);
 
 	std::clock_t begin2 = std::clock();
 	for(size_t d = 0; d != 4; ++d){
@@ -517,7 +517,7 @@ namespace {
 	EXPECT_EQ( xround(res2.at(mkts(3,2))), xround(vres[14]) );
 	EXPECT_EQ( xround(res2.at(mkts(3,3))), xround(vres[15]) );
 
-	std::cout << "std::vector - for loop time: " << static_cast<double>( end2 - begin2 ) / CLOCKS_PER_SEC
+	std::cout << "vector - for loop time: " << static_cast<double>( end2 - begin2 ) / CLOCKS_PER_SEC
 		  << std::endl;
 	std::cout << "ratio: " << static_cast<double>( end - begin ) / static_cast<double>( end2 - begin2 ) << std::endl;
     }    
@@ -543,7 +543,7 @@ namespace {
         MultiArray<double,SRange,SRange,SRange> ma1(sr1ptr, sr2ptr, sr3ptr, v1);
         MultiArray<double,SRange,SRange> ma2(sr3ptr, sr2ptr, v2);
 
-        SubRangeFactory<SRange> subf(sr2ptr, std::vector<size_t>({0,2}));
+        SubRangeFactory<SRange> subf(sr2ptr, vector<size_t>({0,2}));
         auto subptr = MAT::createExplicit(subf);
         
         MultiArray<double,SRange,SRange> res(sr3ptr,sr1ptr,0.);
diff --git a/src/tests/ranges/anonymous_unit_test.cc b/src/tests/ranges/anonymous_unit_test.cc
index b12065e..cba81cb 100644
--- a/src/tests/ranges/anonymous_unit_test.cc
+++ b/src/tests/ranges/anonymous_unit_test.cc
@@ -18,7 +18,7 @@ namespace {
     template <class Factory, typename T>
     void swapFactory(std::shared_ptr<RangeFactoryBase>& fptr, std::initializer_list<T> ilist)
     {
-	std::vector<T> tmp = ilist;
+	vector<T> tmp = ilist;
 	auto nptr = std::make_shared<Factory>( tmp );
 	fptr = nptr;
     }
@@ -160,7 +160,7 @@ namespace {
 
     TEST_F(AnonymousTest, ToString1)
     {
-	std::vector<char> vv = sr1ptr->data();
+	vector<char> vv = sr1ptr->data();
 	const char* dp = vv.data();
 	auto ff = createRangeFactory(&dp);
 	auto ss = std::dynamic_pointer_cast<SRange>( ff->create() );
@@ -183,9 +183,9 @@ namespace {
 	auto sr = std::dynamic_pointer_cast<SpinRange>( srf.create() );
 	auto nr = std::dynamic_pointer_cast<NameRange>( nrf.create() );
 	
-	std::vector<char> cv = cr->data();
-	std::vector<char> sv = sr->data();
-	std::vector<char> nv = nr->data(); 
+	vector<char> cv = cr->data();
+	vector<char> sv = sr->data();
+	vector<char> nv = nr->data(); 
 	const char* cp = cv.data();
 	const char* sp = sv.data();
 	const char* np = nv.data();
@@ -212,7 +212,7 @@ namespace {
 
     TEST_F(AnonymousTest, ToStringMulti)
     {
-	std::vector<char> vv = m3rptr->data();
+	vector<char> vv = m3rptr->data();
 	const char* dp = vv.data();
 	
 	auto ff2 = std::dynamic_pointer_cast<M3RF>( createRangeFactory(&dp) );
@@ -235,7 +235,7 @@ namespace {
     {
 	AnonymousRangeFactory arf2(sr1ptr,m3rptr,sr2ptr);
 	auto ar = std::dynamic_pointer_cast<AnonymousRange>( arf2.create() );
-	std::vector<char> vv = ar->data();
+	vector<char> vv = ar->data();
 	const char* dp = vv.data();
 	
 	auto ff2 = std::dynamic_pointer_cast<AnonymousRangeFactory>( createRangeFactory(&dp) );
diff --git a/src/tests/ranges/index_unit_test.cc b/src/tests/ranges/index_unit_test.cc
index d4289f6..5da0265 100644
--- a/src/tests/ranges/index_unit_test.cc
+++ b/src/tests/ranges/index_unit_test.cc
@@ -17,7 +17,7 @@ namespace {
     template <class Factory, typename T>
     void swapFactory(std::shared_ptr<RangeFactoryBase>& fptr, std::initializer_list<T> ilist)
     {
-	std::vector<T> tmp = ilist;
+	vector<T> tmp = ilist;
 	auto nptr = std::make_shared<Factory>( tmp );
 	fptr = nptr;
     }

From 53c433f8dfc9fdea49dca2879d5fc331a007f2d8 Mon Sep 17 00:00:00 2001
From: Christian Zimmermann <christian.zimmermann@ur.de>
Date: Thu, 14 Feb 2019 14:39:59 +0100
Subject: [PATCH 4/9] vectorizes, but not that optimal...

---
 src/include/arith.h                 |  6 ++++++
 src/include/multi_array_operation.h |  6 +++---
 src/include/type_operations.h       | 27 +++++++++++++++++++++------
 3 files changed, 30 insertions(+), 9 deletions(-)

diff --git a/src/include/arith.h b/src/include/arith.h
index c3d94b7..5bf1a09 100644
--- a/src/include/arith.h
+++ b/src/include/arith.h
@@ -86,6 +86,12 @@ namespace MultiArrayTools
 	{
 	    return a1 + a2;
 	}
+
+	static inline T& sapply(T& a1, T a2)
+	{
+	    return a1 += a2;
+	}
+
     };
 
     template <typename T>
diff --git a/src/include/multi_array_operation.h b/src/include/multi_array_operation.h
index a5a3dfe..31e997f 100644
--- a/src/include/multi_array_operation.h
+++ b/src/include/multi_array_operation.h
@@ -91,9 +91,9 @@ namespace MultiArrayTools
     template <typename T>
     struct SelfIdentity
     {
-        static inline T apply(const T& a, const T& b)
+        static inline T& sapply(T& a, const T& b)
         {
-            return b;
+            return a = b;
         }
     };
 
@@ -140,7 +140,7 @@ namespace MultiArrayTools
 	OperationMaster(T* data, const OpClass& second,
 			IndexType& index, bool doParallel = false);
 
-	inline void set(size_t pos, T val) { mDataPtr[pos] = AOp::apply(mDataPtr[pos],val); }
+	inline void set(size_t pos, T val) { AOp::sapply(mDataPtr[pos],val); }
         
 	//inline void add(size_t pos, T val) { mDataPtr[pos] += val; }
 	inline T get(size_t pos) const;
diff --git a/src/include/type_operations.h b/src/include/type_operations.h
index 564b83b..8508fb5 100644
--- a/src/include/type_operations.h
+++ b/src/include/type_operations.h
@@ -154,25 +154,40 @@ namespace MultiArrayTools
 	friend OperationClass;
     };
 
-    typedef struct v256 { double _x[4]; } v256;
+    typedef struct v256 { alignas(32) double _x[4]; } v256;
 
     template <int N>
     inline void xadd(double* o, const double* a, const double* b)
     {
-	//#pragma omp simd aligned(o, a, b: 16)
+#pragma omp simd aligned(o, a, b: 32)
 	for(int i = 0; i < N; i++) {
 	    o[i] = a[i] + b[i];
 	}  
     }
+
+    template <int N>
+    inline void xsadd(double* o, const double* a)
+    {
+#pragma omp simd aligned(o, a: 32)
+	for(int i = 0; i < N; i++) {
+	    o[i] += a[i];
+	}  
+    }
     
     inline v256 operator+(const v256& a, const v256& b)
     {
-	alignas(32) v256 out;
-	xadd<4>( reinterpret_cast<double*>(&out), reinterpret_cast<const double*>(&a),
+	alignas(32) v256 o;
+	xadd<4>( reinterpret_cast<double*>(&o), reinterpret_cast<const double*>(&a),
 		 reinterpret_cast<const double*>(&b) );
-	return out;
+	return o;
     }
-    /*
+
+    inline v256& operator+=(v256& o, const v256& a)
+    {
+	xsadd<4>( reinterpret_cast<double*>(&o), reinterpret_cast<const double*>(&a) );
+	return o;
+    }
+/*
     inline v256 operator-(const v256& a, const v256& b)
     {
         	alignas(32) v256 out;

From 0559372b0fca251089759a69ec7c284914c85e80 Mon Sep 17 00:00:00 2001
From: Christian Zimmermann <christian.zimmermann@ur.de>
Date: Thu, 14 Feb 2019 16:20:15 +0100
Subject: [PATCH 5/9] march=native

---
 CMakeLists.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 77f8ffb..805c8e3 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 2.8)
 
 project(multi_array)
 
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Werror -std=c++11 -g -Wpedantic -Ofast -mavx2 -funroll-loops -fopenmp")
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Werror -std=c++11 -g -Wpedantic -Ofast -march=native -funroll-loops -fopenmp")
 
 enable_testing()
 

From 59ffd5adcc0dec6d11ebbd7cff5eaa9c30b798e3 Mon Sep 17 00:00:00 2001
From: Christian Zimmermann <christian.zimmermann@ur.de>
Date: Thu, 14 Feb 2019 22:04:33 +0100
Subject: [PATCH 6/9] use omp simd in performance comparision loop

---
 src/tests/op_perf_test.cc | 14 ++++++++++----
 1 file changed, 10 insertions(+), 4 deletions(-)

diff --git a/src/tests/op_perf_test.cc b/src/tests/op_perf_test.cc
index 90af492..5f5de8c 100644
--- a/src/tests/op_perf_test.cc
+++ b/src/tests/op_perf_test.cc
@@ -195,16 +195,22 @@ namespace {
 	    }
 	}                  
 	std::clock_t begin2 = std::clock();
+	double* vrptr = vres.data();
+	double* dptr = data.data();
+	
         for(size_t j = 0; j != os; ++j) {
             for(size_t i = 0; i != os; ++i){
                 for(size_t a = 0; a != 4; ++a){
                     for(size_t b = 0; b != 4; ++b){
                         for(size_t c = 0; c != 4; ++c){
                             for(size_t d = 0; d != 4; ++d){
-                                for(size_t p = 0; p != 4; ++p){
-                                    const size_t tidx = i*4*4 + d*4 + p;
-                                    const size_t sidx = /*i*65536 +*/ d*4*4*4*4*4*4*4 + a*5*4*4*4*4*4 + b*5*4*4*4 + c*5*4 + p;
-                                    vres[tidx] += data[sidx];
+				const size_t tidx = i*4*4 + d*4;
+				const size_t sidx = /*i*65536 +*/ d*4*4*4*4*4*4*4 + a*5*4*4*4*4*4 + b*5*4*4*4 + c*5*4;
+				double* xvrptr = vrptr + tidx;
+				double* xdptr = dptr + sidx;
+#pragma omp simd aligned(xvrptr, xdptr: 32)
+                                for(int p = 0; p < 4; p++){
+                                    xvrptr[p] += xdptr[p];
                                 }
                             }
                         }

From 6f320fe9e96abf10e190e667928d6430d3811745 Mon Sep 17 00:00:00 2001
From: Christian Zimmermann <christian.zimmermann@ur.de>
Date: Tue, 26 Feb 2019 18:56:57 +0100
Subject: [PATCH 7/9] still not better...

---
 src/include/map_range.cc.h             |   4 +-
 src/include/map_range.h                |   4 +-
 src/include/multi_array_operation.cc.h | 123 +++++++++++++++++++------
 src/include/multi_array_operation.h    |  89 +++++++++++++-----
 src/include/pack_num.h                 |   6 +-
 src/include/type_operations.h          |  19 ++--
 src/include/xfor/xfor.h                |  72 +++++++--------
 7 files changed, 212 insertions(+), 105 deletions(-)

diff --git a/src/include/map_range.cc.h b/src/include/map_range.cc.h
index 2f661a6..939657c 100644
--- a/src/include/map_range.cc.h
+++ b/src/include/map_range.cc.h
@@ -29,7 +29,7 @@ namespace MultiArrayTools
 
     template <class MapF, class IndexPack, class Expr, SpaceType STYPE>
     inline void OpExpr<MapF,IndexPack,Expr,STYPE>::operator()(size_t mlast,
-							 ExtType last) const
+                                                              ExtType last)
     {
 	constexpr size_t NEXT = Op::SIZE;
 	const ExtType nxpos = last;
@@ -42,7 +42,7 @@ namespace MultiArrayTools
     }
 
     template <class MapF, class IndexPack, class Expr, SpaceType STYPE>
-    inline void OpExpr<MapF,IndexPack,Expr,STYPE>::operator()(size_t mlast) const
+    inline void OpExpr<MapF,IndexPack,Expr,STYPE>::operator()(size_t mlast)
     {
 	const ExtType last;
 	constexpr size_t NEXT = Op::SIZE;
diff --git a/src/include/map_range.h b/src/include/map_range.h
index 0406c7c..129dd10 100644
--- a/src/include/map_range.h
+++ b/src/include/map_range.h
@@ -66,8 +66,8 @@ namespace MultiArrayTools
 	
 	OpExpr(const MapF& mapf, const IndexPack& ipack, const std::shared_ptr<OIType>& oind, size_t step, Expr ex);
 
-	inline void operator()(size_t mlast, ExtType last) const;
-	inline void operator()(size_t mlast = 0) const;
+	inline void operator()(size_t mlast, ExtType last);
+	inline void operator()(size_t mlast = 0);
 
 	auto rootSteps(std::intptr_t iPtrNum = 0) const -> ExtType;
 
diff --git a/src/include/multi_array_operation.cc.h b/src/include/multi_array_operation.cc.h
index 4b1abed..83e1cc0 100644
--- a/src/include/multi_array_operation.cc.h
+++ b/src/include/multi_array_operation.cc.h
@@ -109,18 +109,25 @@ namespace MultiArrayTools
     /*****************************************
      *   OperationMaster::AssignmentExpr     *
      *****************************************/
-
+    /*
     template <typename T, class AOp, class OpClass, class... Ranges>
     OperationMaster<T,AOp,OpClass,Ranges...>::AssignmentExpr::
     AssignmentExpr(OperationMaster& m, const OpClass& sec) :
 	mM(m), mSec(sec) {}
-
+    */
+    template <typename T, class AOp, class OpClass, class... Ranges>
+    OperationMaster<T,AOp,OpClass,Ranges...>::AssignmentExpr::
+    AssignmentExpr(T* dataPtr, const OpClass& sec) :
+        mSec(sec), mDataPtr(dataPtr) {}
+    
     template <typename T, class AOp, class OpClass, class... Ranges>
     inline void OperationMaster<T,AOp,OpClass,Ranges...>::AssignmentExpr::
-    operator()(size_t start, ExtType last) const
+    operator()(size_t start, ExtType last)
     {
 	//VCHECK(mSec.template get<ExtType>(last));
-	mM.set(start, mSec.template get<ExtType>(last) );
+	//mM.set(start, mSec.template get<ExtType>(last) );
+        mDataPtr[start] = mSec.template get<ExtType>(last);
+        //AOp::sapply(mDataPtr[start], mSec.template get<ExtType>(last));
     }
     
     template <typename T, class AOp, class OpClass, class... Ranges>
@@ -130,7 +137,30 @@ namespace MultiArrayTools
     {
 	return mSec.rootSteps(iPtrNum);
     }
+
+    template <typename T, class AOp, class OpClass, class... Ranges>
+    OperationMaster<T,AOp,OpClass,Ranges...>::AddExpr::
+    AddExpr(T* dataPtr, const OpClass& sec) :
+        mSec(sec), mDataPtr(dataPtr) {}
     
+    template <typename T, class AOp, class OpClass, class... Ranges>
+    inline void OperationMaster<T,AOp,OpClass,Ranges...>::AddExpr::
+    operator()(size_t start, ExtType last)
+    {
+	//VCHECK(mSec.template get<ExtType>(last));
+	//mM.set(start, mSec.template get<ExtType>(last) );
+        mDataPtr[start] += mSec.template get<ExtType>(last);
+        //AOp::sapply(mDataPtr[start], mSec.template get<ExtType>(last));
+    }
+    
+    template <typename T, class AOp, class OpClass, class... Ranges>
+    typename OperationMaster<T,AOp,OpClass,Ranges...>::AddExpr::ExtType
+    OperationMaster<T,AOp,OpClass,Ranges...>::AddExpr::
+    rootSteps(std::intptr_t iPtrNum) const
+    {
+	return mSec.rootSteps(iPtrNum);
+    }
+
     
     /*************************
      *   OperationMaster     *
@@ -159,15 +189,17 @@ namespace MultiArrayTools
     template <typename T, class AOp, class OpClass, class... Ranges>
     void OperationMaster<T,AOp,OpClass,Ranges...>::performAssignment(std::intptr_t blockIndexNum)
     {
+        /*
 	AssignmentExpr ae(*this, mSecond); // Expression to be executed within loop
         if(mDoParallel){
-            const auto ploop = mIndex.pifor( 1, mSecond.loop(ae) );
+            auto ploop = mIndex.pifor( 1, mSecond.loop(ae) );
             ploop(); // execute overall loop(s) and so internal hidden loops and so the inherited expressions
         }
         else {
-            const auto loop = mIndex.ifor( 1, mSecond.loop(ae) );
+            auto loop = mIndex.ifor( 1, mSecond.loop(ae) );
             loop(); // execute overall loop(s) and so internal hidden loops and so the inherited expressions
         }
+        */
     }
 
     template <typename T, class AOp, class OpClass, class... Ranges>
@@ -186,11 +218,12 @@ namespace MultiArrayTools
     ConstOperationRoot(const MultiArrayBase<T,Ranges...>& ma,
 		       const std::shared_ptr<typename Ranges::IndexType>&... indices) :
 	mDataPtr(ma.data()),
+        mOrigDataPtr(ma.data()),
 	mIndex( ma.begin() )
     {
 	//VCHECK(ma.data());
 	mIndex(indices...);
-	mOff = mIndex.pos();
+	//mOff = mIndex.pos();
     }
 
     template <typename T, class... Ranges>
@@ -198,35 +231,37 @@ namespace MultiArrayTools
     ConstOperationRoot(std::shared_ptr<MultiArrayBase<T,Ranges...> > maptr,
 		       const std::shared_ptr<typename Ranges::IndexType>&... indices) :
 	mDataPtr(maptr->data()),
+        mOrigDataPtr(maptr->data()),
 	mIndex(maptr->begin()),
 	mMaPtr(maptr)
     {
 	mIndex(indices...);
-	mOff = mIndex.pos();
+	//mOff = mIndex.pos();
     }
 
     template <typename T, class... Ranges>
     ConstOperationRoot<T,Ranges...>::
     ConstOperationRoot(const T* data, const IndexType& ind) :
 	mDataPtr(data),
+        mOrigDataPtr(data),
 	mIndex( ind )
     {
-	mOff = mIndex.pos();
+	//mOff = mIndex.pos();
     }
 
     template <typename T, class... Ranges>
     template <class ET>
     inline T ConstOperationRoot<T,Ranges...>::get(ET pos) const
     {
-	return mDataPtr[pos.val()+mOff];
+	return mDataPtr[pos.val()/*+mOff*/];
     }
 
     template <typename T, class... Ranges>
     template <class ET>
-    inline const ConstOperationRoot<T,Ranges...>& ConstOperationRoot<T,Ranges...>::set(ET pos) const
+    inline ConstOperationRoot<T,Ranges...>& ConstOperationRoot<T,Ranges...>::set(ET pos)
     {
 	mIndex = pos.val();
-	mOff = mIndex.pos();
+	mDataPtr = mOrigDataPtr + mIndex.pos();
 	return *this;
     }
 
@@ -267,7 +302,7 @@ namespace MultiArrayTools
 
     template <typename T, class Op>
     template <class ET>
-    inline const StaticCast<T,Op>& StaticCast<T,Op>::set(ET pos) const
+    inline StaticCast<T,Op>& StaticCast<T,Op>::set(ET pos)
     {
 	mOp.set(pos);
 	return *this;
@@ -311,7 +346,7 @@ namespace MultiArrayTools
 
     template <class Range>
     template <class ET>
-    inline const MetaOperationRoot<Range>& MetaOperationRoot<Range>::set(ET pos) const
+    inline MetaOperationRoot<Range>& MetaOperationRoot<Range>::set(ET pos)
     {
         assert(0);
 	//(*mIndex) = pos.val();
@@ -342,38 +377,65 @@ namespace MultiArrayTools
     OperationRoot(MutableMultiArrayBase<T,Ranges...>& ma,
 		  const std::shared_ptr<typename Ranges::IndexType>&... indices) :
 	mDataPtr(ma.data()),
+        mOrigDataPtr(ma.data()),
 	mIndex( ma.begin() )
     {
 	mIndex(indices...);
-	mOff = mIndex.pos();
+	//mOff = mIndex.pos();
     }
 
     template <typename T, class... Ranges>
     OperationRoot<T,Ranges...>::
     OperationRoot(T* data, const IndexType& ind) :
 	mDataPtr(data),
+        mOrigDataPtr(data),
 	mIndex( ind )
     {
-	mOff = mIndex.pos();
+	//mOff = mIndex.pos();
     }
 
     template <typename T, class... Ranges>
     template <class OpClass>
-    OperationMaster<T,SelfIdentity<T>,OpClass,Ranges...> OperationRoot<T,Ranges...>::operator=(const OpClass& in)
+    OperationRoot<T,Ranges...>& OperationRoot<T,Ranges...>::operator=(const OpClass& in)
+    //OperationMaster<T,SelfIdentity<T>,OpClass,Ranges...> OperationRoot<T,Ranges...>::operator=(const OpClass& in)
     {
-	return OperationMaster<T,SelfIdentity<T>,OpClass,Ranges...>(mDataPtr, in, mIndex, mDoParallel);
+	typename OperationMaster<T,SelfIdentity<T>,OpClass,Ranges...>::AssignmentExpr ae(mDataPtr, in);
+        // Expression to be executed within loop
+        if(mDoParallel){
+            auto ploop = mIndex.pifor( 1, in.loop(ae) );
+            ploop(); // execute overall loop(s) and so internal hidden loops and so the inherited expressions
+        }
+        else {
+            auto loop = mIndex.ifor( 1, in.loop(ae) );
+            loop(); // execute overall loop(s) and so internal hidden loops and so the inherited expressions
+        }
+        return *this;
+	//return OperationMaster<T,SelfIdentity<T>,OpClass,Ranges...>(mDataPtr, in, mIndex, mDoParallel);
     }
 
     template <typename T, class... Ranges>
     template <class OpClass>
-    OperationMaster<T,plus<T>,OpClass,Ranges...> OperationRoot<T,Ranges...>::operator+=(const OpClass& in)
+    OperationRoot<T,Ranges...>& OperationRoot<T,Ranges...>::operator+=(const OpClass& in)
+    //OperationMaster<T,plus<T>,OpClass,Ranges...> OperationRoot<T,Ranges...>::operator+=(const OpClass& in)
     {
-	return OperationMaster<T,plus<T>,OpClass,Ranges...>(mDataPtr, in, mIndex, mDoParallel);
+        typename OperationMaster<T,plus<T>,OpClass,Ranges...>::AddExpr ae(mDataPtr, in);
+        // Expression to be executed within loop
+        if(mDoParallel){
+            auto ploop = mIndex.pifor( 1, in.loop(ae) );
+            ploop(); // execute overall loop(s) and so internal hidden loops and so the inherited expressions
+        }
+        else {
+            auto loop = mIndex.ifor( 1, in.loop(ae) );
+            loop(); // execute overall loop(s) and so internal hidden loops and so the inherited expressions
+        }
+        return *this;
+	//return OperationMaster<T,plus<T>,OpClass,Ranges...>(mDataPtr, in, mIndex, mDoParallel);
     }
 
     template <typename T, class... Ranges>
-    OperationMaster<T,SelfIdentity<T>,OperationRoot<T,Ranges...>,Ranges...>
-    OperationRoot<T,Ranges...>::operator=(const OperationRoot<T,Ranges...>& in)
+    OperationRoot<T,Ranges...>& OperationRoot<T,Ranges...>::operator=(const OperationRoot<T,Ranges...>& in)
+    //OperationMaster<T,SelfIdentity<T>,OperationRoot<T,Ranges...>,Ranges...>
+    //OperationRoot<T,Ranges...>::operator=(const OperationRoot<T,Ranges...>& in)
     {
 	return operator=<OperationRoot<T,Ranges...> >(in);
     }
@@ -389,15 +451,15 @@ namespace MultiArrayTools
     template <class ET>
     inline T OperationRoot<T,Ranges...>::get(ET pos) const
     {
-	return mDataPtr[pos.val()+mOff];
+	return mDataPtr[pos.val()/*+mOff*/];
     }
 
     template <typename T, class... Ranges>
     template <class ET>
-    inline const OperationRoot<T,Ranges...>& OperationRoot<T,Ranges...>::set(ET pos) const
+    inline OperationRoot<T,Ranges...>& OperationRoot<T,Ranges...>::set(ET pos)
     {
 	mIndex = pos.val();
-	mOff = mIndex.pos();
+	mDataPtr = mOrigDataPtr + mIndex.pos();
 	return *this;
     }
 
@@ -418,7 +480,8 @@ namespace MultiArrayTools
     template <typename T, class... Ranges>
     T* OperationRoot<T,Ranges...>::data() const
     {
-	return mDataPtr + mIndex().pos();
+        auto i = mIndex;
+	return mDataPtr + i().pos();
     }
 
     template <typename T, class... Ranges>
@@ -448,7 +511,7 @@ namespace MultiArrayTools
 
     template <typename T>
     template <class ET>
-    inline const OperationValue<T>& OperationValue<T>::set(ET pos) const
+    inline OperationValue<T>& OperationValue<T>::set(ET pos)
     {
 	return *this;
     }
@@ -498,7 +561,7 @@ namespace MultiArrayTools
 
     template <typename T, class OpFunction, class... Ops>
     template <class ET>
-    inline const Operation<T,OpFunction,Ops...>& Operation<T,OpFunction,Ops...>::set(ET pos) const
+    inline Operation<T,OpFunction,Ops...>& Operation<T,OpFunction,Ops...>::set(ET pos)
     {
 	PackNum<sizeof...(Ops)-1>::setOpPos(mOps,pos);
 	return *this;
@@ -539,7 +602,7 @@ namespace MultiArrayTools
 
     template <typename T, class Op, class IndexType>
     template <class ET>
-    inline const Contraction<T,Op,IndexType>& Contraction<T,Op,IndexType>::set(ET pos) const
+    inline Contraction<T,Op,IndexType>& Contraction<T,Op,IndexType>::set(ET pos)
     {
 	mOp.set(pos);
 	return *this;
@@ -586,7 +649,7 @@ namespace MultiArrayTools
 
     template <typename T, class Op, class... Indices>
     template <class ET>
-    inline const SliceContraction<T,Op,Indices...>& SliceContraction<T,Op,Indices...>::set(ET pos) const
+    inline SliceContraction<T,Op,Indices...>& SliceContraction<T,Op,Indices...>::set(ET pos)
     {
 	mOp.set(pos);
 	return *this;
diff --git a/src/include/multi_array_operation.h b/src/include/multi_array_operation.h
index 31e997f..89846d5 100644
--- a/src/include/multi_array_operation.h
+++ b/src/include/multi_array_operation.h
@@ -91,7 +91,7 @@ namespace MultiArrayTools
     template <typename T>
     struct SelfIdentity
     {
-        static inline T& sapply(T& a, const T& b)
+        static inline T& sapply(T& a, T b)
         {
             return a = b;
         }
@@ -107,27 +107,57 @@ namespace MultiArrayTools
 	private:
 	    AssignmentExpr() = default;
 	    	    
-	    OperationMaster& mM;
-	    const OpClass& mSec;
-
+	    //OperationMaster mM;
+	    OpClass mSec;
+            T* mDataPtr;
+            
 	public:
 
 	    static constexpr size_t LAYER = 0;
 	    static constexpr size_t SIZE = OpClass::SIZE;
 	    typedef decltype(mSec.rootSteps()) ExtType;
 	    
-	    AssignmentExpr(OperationMaster& m, const OpClass& sec);
+	    //AssignmentExpr(OperationMaster& m, const OpClass& sec);
+            AssignmentExpr(T* dataPtr, const OpClass& sec);
 
 	    AssignmentExpr(const AssignmentExpr& in) = default;
 	    AssignmentExpr(AssignmentExpr&& in) = default;
 	    
-	    inline void operator()(size_t start = 0) const; 
-	    inline void operator()(size_t start, ExtType last) const;
+	    inline void operator()(size_t start = 0); 
+	    inline void operator()(size_t start, ExtType last);
 
 	    auto rootSteps(std::intptr_t iPtrNum = 0) const -> ExtType;
 	    
 	};
-	
+
+        class AddExpr
+	{
+	private:
+	    AddExpr() = default;
+	    	    
+	    //OperationMaster mM;
+	    OpClass mSec;
+            T* mDataPtr;
+            
+	public:
+
+	    static constexpr size_t LAYER = 0;
+	    static constexpr size_t SIZE = OpClass::SIZE;
+	    typedef decltype(mSec.rootSteps()) ExtType;
+	    
+	    //AssignmentExpr(OperationMaster& m, const OpClass& sec);
+            AddExpr(T* dataPtr, const OpClass& sec);
+
+	    AddExpr(const AddExpr& in) = default;
+	    AddExpr(AddExpr&& in) = default;
+	    
+	    inline void operator()(size_t start = 0); 
+	    inline void operator()(size_t start, ExtType last);
+
+	    auto rootSteps(std::intptr_t iPtrNum = 0) const -> ExtType;
+	    
+	};
+
 	typedef T value_type;
 	//typedef OperationBase<T> OB;
 	typedef ContainerRange<T,Ranges...> CRange;
@@ -180,7 +210,7 @@ namespace MultiArrayTools
 	inline T get(ET pos) const;
 
 	template <class ET>
-	inline const ConstOperationRoot& set(ET pos) const;
+	inline ConstOperationRoot& set(ET pos);
 
 	MExt<void> rootSteps(std::intptr_t iPtrNum = 0) const; // nullptr for simple usage with decltype
 
@@ -193,8 +223,9 @@ namespace MultiArrayTools
 
 	//MultiArrayBase<T,Ranges...> const& mArrayRef;
 	const T* mDataPtr;
-	mutable IndexType mIndex;
-	mutable size_t mOff = 0;
+        const T* mOrigDataPtr;
+	IndexType mIndex;
+	//size_t mOff = 0;
 	std::shared_ptr<MultiArrayBase<T,Ranges...> > mMaPtr; // never remove this ptr, otherwise we lose temporary container instances!
     };
 
@@ -202,7 +233,7 @@ namespace MultiArrayTools
     class StaticCast : public OperationTemplate<T,StaticCast<T,Op> >
     {
     private:
-	const Op& mOp;
+	Op mOp;
 
     public:
 	
@@ -219,7 +250,7 @@ namespace MultiArrayTools
 	inline T get(ET pos) const;
 
 	template <class ET>
-	inline const StaticCast& set(ET pos) const;
+	inline StaticCast& set(ET pos);
 
 	auto rootSteps(std::intptr_t iPtrNum = 0) const
 	    -> decltype(mOp.rootSteps(iPtrNum));
@@ -253,7 +284,7 @@ namespace MultiArrayTools
 	inline value_type get(ET pos) const;
 
 	template <class ET>
-	inline const MetaOperationRoot& set(ET pos) const;
+	inline MetaOperationRoot& set(ET pos);
 
 	MExt<void> rootSteps(std::intptr_t iPtrNum = 0) const; // nullptr for simple usage with decltype
 
@@ -284,7 +315,7 @@ namespace MultiArrayTools
 		      const std::shared_ptr<typename Ranges::IndexType>&... indices);
 
 	OperationRoot(T* data, const IndexType& ind);
-
+        /*
 	template <class OpClass>
 	OperationMaster<T,SelfIdentity<T>,OpClass,Ranges...> operator=(const OpClass& in);
 
@@ -292,14 +323,23 @@ namespace MultiArrayTools
 	OperationMaster<T,plus<T>,OpClass,Ranges...> operator+=(const OpClass& in);
 
 	OperationMaster<T,SelfIdentity<T>,OperationRoot,Ranges...> operator=(const OperationRoot& in);
+        */
 
+        template <class OpClass>
+        OperationRoot& operator=(const OpClass& in);
+
+        template <class OpClass>
+        OperationRoot& operator+=(const OpClass& in);
+
+        OperationRoot& operator=(const OperationRoot& in);
+        
         OperationRoot& par();
         
 	template <class ET>
 	inline T get(ET pos) const;
 
 	template <class ET>
-	inline const OperationRoot& set(ET pos) const;
+	inline OperationRoot& set(ET pos);
 
 	MExt<void> rootSteps(std::intptr_t iPtrNum = 0) const; // nullptr for simple usage with decltype
 
@@ -316,8 +356,9 @@ namespace MultiArrayTools
 
 	//MutableMultiArrayBase<T,Ranges...>& mArrayRef;
 	T* mDataPtr;
-	mutable IndexType mIndex;
- 	mutable size_t mOff = 0;
+        T* mOrigDataPtr;
+	IndexType mIndex;
+ 	//size_t mOff = 0;
         bool mDoParallel = false;
     };
 
@@ -338,7 +379,7 @@ namespace MultiArrayTools
 	inline T get(ET pos) const;
 
 	template <class ET>
-	inline const OperationValue& set(ET pos) const;
+	inline OperationValue& set(ET pos);
 
 	MExt<void> rootSteps(std::intptr_t iPtrNum = 0) const; // nullptr for simple usage with decltype
 
@@ -415,7 +456,7 @@ namespace MultiArrayTools
 	inline T get(ET pos) const;
 
 	template <class ET>
-	inline const Operation& set(ET pos) const;
+	inline Operation& set(ET pos);
 
 	auto rootSteps(std::intptr_t iPtrNum = 0) const // nullptr for simple usage with decltype
 	    -> decltype(PackNum<sizeof...(Ops)-1>::mkSteps(iPtrNum, mOps));
@@ -473,7 +514,7 @@ namespace MultiArrayTools
 	
     private:
 
-	const Op& mOp;
+	Op mOp;
 	std::shared_ptr<IndexType> mInd;
 
     public:
@@ -485,7 +526,7 @@ namespace MultiArrayTools
 	inline T get(ET pos) const;
 
 	template <class ET>
-	inline const Contraction& set(ET pos) const;
+	inline Contraction& set(ET pos);
 
 	auto rootSteps(std::intptr_t iPtrNum = 0) const  // nullptr for simple usage with decltype
 	    -> decltype(mOp.rootSteps(iPtrNum));
@@ -509,7 +550,7 @@ namespace MultiArrayTools
 
     private:
 
-	const Op& mOp;
+	mutable Op mOp;
 	mutable std::shared_ptr<MultiArray<T,typename Indices::RangeType...> > mCont;
 	mutable OperationRoot<T,typename Indices::RangeType...> mTarOp;
 		
@@ -522,7 +563,7 @@ namespace MultiArrayTools
 	inline const value_type& get(ET pos) const;
 
 	template <class ET>
-	inline const SliceContraction& set(ET pos) const;
+	inline SliceContraction& set(ET pos);
 
 	auto rootSteps(std::intptr_t iPtrNum = 0) const  // nullptr for simple usage with decltype
 	    -> decltype(mOp.rootSteps(iPtrNum));
diff --git a/src/include/pack_num.h b/src/include/pack_num.h
index d9bf56e..d369b17 100644
--- a/src/include/pack_num.h
+++ b/src/include/pack_num.h
@@ -71,7 +71,7 @@ namespace MultiArrayHelper
 	}
 
 	template <size_t LAST, typename T, class ETuple, class OpTuple, class OpFunction, typename... Args>
-	static inline T mkOpExpr(std::shared_ptr<OpFunction> f, const ETuple& pos, const OpTuple& ops, const Args&... args)
+	static inline T mkOpExpr(std::shared_ptr<OpFunction> f, const ETuple& pos, const OpTuple& ops, Args... args)
 	{
 	    typedef typename std::remove_reference<decltype(std::get<N>(ops))>::type NextOpType;
 	    static_assert(LAST > NextOpType::SIZE, "inconsistent array positions");
@@ -116,7 +116,7 @@ namespace MultiArrayHelper
 	}
 
 	template <size_t LAST,class OpTuple, class ETuple>
-	static inline void setOpPos(const OpTuple& ot, const ETuple& et)
+	static inline void setOpPos(OpTuple& ot, const ETuple& et)
 	{
 	    typedef typename std::remove_reference<decltype(std::get<N>(ot))>::type NextOpType;
 	    static_assert(LAST > NextOpType::SIZE, "inconsistent array positions");
@@ -203,7 +203,7 @@ namespace MultiArrayHelper
 	}
 
 	template <size_t LAST,class OpTuple, class ETuple>
-	static inline void setOpPos(const OpTuple& ot, const ETuple& et)
+	static inline void setOpPos(OpTuple& ot, const ETuple& et)
 	{
 	    typedef typename std::remove_reference<decltype(std::get<0>(et))>::type NextOpType;
 	    static_assert(LAST > NextOpType::SIZE, "inconsistent array positions");
diff --git a/src/include/type_operations.h b/src/include/type_operations.h
index 8508fb5..8f970d9 100644
--- a/src/include/type_operations.h
+++ b/src/include/type_operations.h
@@ -154,7 +154,10 @@ namespace MultiArrayTools
 	friend OperationClass;
     };
 
-    typedef struct v256 { alignas(32) double _x[4]; } v256;
+    struct v256
+    {
+        alignas(32) double _x[4];
+    };
 
     template <int N>
     inline void xadd(double* o, const double* a, const double* b)
@@ -173,16 +176,16 @@ namespace MultiArrayTools
 	    o[i] += a[i];
 	}  
     }
-    
-    inline v256 operator+(const v256& a, const v256& b)
+
+    inline v256 operator+(v256 a, v256 b)
     {
-	alignas(32) v256 o;
+	v256 o;
 	xadd<4>( reinterpret_cast<double*>(&o), reinterpret_cast<const double*>(&a),
 		 reinterpret_cast<const double*>(&b) );
 	return o;
     }
 
-    inline v256& operator+=(v256& o, const v256& a)
+    inline v256& operator+=(v256& o, v256 a)
     {
 	xsadd<4>( reinterpret_cast<double*>(&o), reinterpret_cast<const double*>(&a) );
 	return o;
@@ -190,7 +193,7 @@ namespace MultiArrayTools
 /*
     inline v256 operator-(const v256& a, const v256& b)
     {
-        	alignas(32) v256 out;
+        	v256 out;
 #pragma omp simd aligned(outp, ap, bp: 32)
 	for(int i = 0; i < IN; ++i){
 	    outp[i] = ap[i] - bp[i];
@@ -200,7 +203,7 @@ namespace MultiArrayTools
 
     inline v256 operator*(const v256& a, const v256& b)
     {
-	alignas(32) v256 out;
+	v256 out;
 #pragma omp simd aligned(outp, ap, bp: 32)
 	for(int i = 0; i < IN; ++i){
 	    outp[i] = ap[i] * bp[i];
@@ -210,7 +213,7 @@ namespace MultiArrayTools
     
     inline v256 operator/(const v256& a, const v256& b)
     {
-	alignas(32) v256 out;
+	v256 out;
 #pragma omp simd aligned(outp, ap, bp: 32)
 	for(int i = 0; i < IN; ++i){
 	    outp[i] = ap[i] / bp[i];
diff --git a/src/include/xfor/xfor.h b/src/include/xfor/xfor.h
index 3cda525..7af5458 100644
--- a/src/include/xfor/xfor.h
+++ b/src/include/xfor/xfor.h
@@ -35,8 +35,8 @@ namespace MultiArrayHelper
 	ExpressionBase& operator=(const ExpressionBase& in) = default;
 	ExpressionBase& operator=(ExpressionBase&& in) = default;
 
-	virtual void operator()(size_t mlast, DExt last) const = 0;
-	virtual void operator()(size_t mlast = 0) const = 0;
+	virtual void operator()(size_t mlast, DExt last) = 0;
+	virtual void operator()(size_t mlast = 0) = 0;
 
         virtual DExt dRootSteps(std::intptr_t iPtrNum = 0) const = 0;
         virtual DExt dExtension() const = 0;
@@ -127,9 +127,9 @@ namespace MultiArrayHelper
 			 Expr expr);
 
 	
-	inline void operator()(size_t mlast, DExt last) const override final;
-	inline void operator()(size_t mlast, ExtType last) const;
-	inline void operator()(size_t mlast = 0) const override final;
+	inline void operator()(size_t mlast, DExt last) override final;
+	inline void operator()(size_t mlast, ExtType last);
+	inline void operator()(size_t mlast = 0) override final;
 
         DExt dRootSteps(std::intptr_t iPtrNum = 0) const override final;
         DExt dExtension() const override final;
@@ -176,9 +176,9 @@ namespace MultiArrayHelper
                 const vector<size_t>* subset, Expr expr);
 
 	
-	inline void operator()(size_t mlast, DExt last) const override final;
-	inline void operator()(size_t mlast, ExtType last) const;
-	inline void operator()(size_t mlast = 0) const override final;
+	inline void operator()(size_t mlast, DExt last) override final;
+	inline void operator()(size_t mlast, ExtType last) ;
+	inline void operator()(size_t mlast = 0) override final;
 
         DExt dRootSteps(std::intptr_t iPtrNum = 0) const override final;
         DExt dExtension() const override final;
@@ -224,9 +224,9 @@ namespace MultiArrayHelper
 	For(const IndexClass* indPtr,
 	    size_t step, Expr expr);
 
-	inline void operator()(size_t mlast, DExt last) const override final;
-	inline void operator()(size_t mlast, ExtType last) const;
-	inline void operator()(size_t mlast = 0) const override final;
+	inline void operator()(size_t mlast, DExt last) override final;
+	inline void operator()(size_t mlast, ExtType last) ;
+	inline void operator()(size_t mlast = 0) override final;
 
         PFor<IndexClass,Expr> parallel() const;
         
@@ -272,9 +272,9 @@ namespace MultiArrayHelper
 	PFor(const IndexClass* indPtr,
 	    size_t step, Expr expr);
 
-	inline void operator()(size_t mlast, DExt last) const override final;
-	inline void operator()(size_t mlast, ExtType last) const;
-	inline void operator()(size_t mlast = 0) const override final;
+	inline void operator()(size_t mlast, DExt last) override final;
+	inline void operator()(size_t mlast, ExtType last) ;
+	inline void operator()(size_t mlast = 0) override final;
 
         DExt dRootSteps(std::intptr_t iPtrNum = 0) const override final;
         DExt dExtension() const override final;
@@ -311,8 +311,8 @@ namespace MultiArrayHelper
         template <class Expr>
         DynamicExpression(Expr ex) : mNext( std::make_shared<Expr>(ex) ) {}
         
-	inline void operator()(size_t mlast, DExt last) const override final;
-	inline void operator()(size_t mlast = 0) const override final;
+	inline void operator()(size_t mlast, DExt last) override final;
+	inline void operator()(size_t mlast = 0) override final;
 
         inline DExt dRootSteps(std::intptr_t iPtrNum = 0) const override final;
         inline DExt dExtension() const override final;
@@ -345,9 +345,9 @@ namespace MultiArrayHelper
 
 	ExpressionHolder(DynamicExpression expr);
 	
-	inline void operator()(size_t mlast, DExt last) const override final;
-	inline void operator()(size_t mlast, ExtType last) const;
-	inline void operator()(size_t mlast = 0) const override final;
+	inline void operator()(size_t mlast, DExt last) override final;
+	inline void operator()(size_t mlast, ExtType last) ;
+	inline void operator()(size_t mlast = 0) override final;
 
         DExt dRootSteps(std::intptr_t iPtrNum = 0) const override final;
         DExt dExtension() const override final;
@@ -391,14 +391,14 @@ namespace MultiArrayHelper
     }
 
     template <class IndexClass, class Expr, ForType FT>
-    inline void For<IndexClass,Expr,FT>::operator()(size_t mlast, DExt last) const
+    inline void For<IndexClass,Expr,FT>::operator()(size_t mlast, DExt last)
     {
         operator()(mlast, *reinterpret_cast<ExtType const*>(last.first));
     }
     
     template <class IndexClass, class Expr, ForType FT>
     inline void For<IndexClass,Expr,FT>::operator()(size_t mlast,
-						    ExtType last) const
+						    ExtType last)
     {
 	typedef typename IndexClass::RangeType RangeType;
 	for(size_t pos = 0u; pos != ForBound<RangeType::ISSTATIC>::template bound<RangeType::SIZE>(mMax); ++pos){
@@ -411,7 +411,7 @@ namespace MultiArrayHelper
     }
 
     template <class IndexClass, class Expr, ForType FT>
-    inline void For<IndexClass,Expr,FT>::operator()(size_t mlast) const
+    inline void For<IndexClass,Expr,FT>::operator()(size_t mlast)
     {
 	typedef typename IndexClass::RangeType RangeType;
 	const ExtType last;
@@ -484,14 +484,14 @@ namespace MultiArrayHelper
     }
 
     template <class IndexClass, class Expr>
-    inline void PFor<IndexClass,Expr>::operator()(size_t mlast, DExt last) const
+    inline void PFor<IndexClass,Expr>::operator()(size_t mlast, DExt last)
     {
         operator()(mlast, *reinterpret_cast<ExtType const*>(last.first));
     }
     
     template <class IndexClass, class Expr>
     inline void PFor<IndexClass,Expr>::operator()(size_t mlast,
-						    ExtType last) const
+						    ExtType last)
     {
         CHECK;
 	typedef typename IndexClass::RangeType RangeType;
@@ -511,7 +511,7 @@ namespace MultiArrayHelper
     }
 
     template <class IndexClass, class Expr>
-    inline void PFor<IndexClass,Expr>::operator()(size_t mlast) const
+    inline void PFor<IndexClass,Expr>::operator()(size_t mlast)
     {
         CHECK;
 	typedef typename IndexClass::RangeType RangeType;
@@ -588,14 +588,14 @@ namespace MultiArrayHelper
     }
 
     template <class IndexClass, class Expr>
-    inline void SingleExpression<IndexClass,Expr>::operator()(size_t mlast, DExt last) const
+    inline void SingleExpression<IndexClass,Expr>::operator()(size_t mlast, DExt last)
     {
         operator()(mlast, *reinterpret_cast<ExtType const*>(last.first));
     }
 
     template <class IndexClass, class Expr>
     inline void SingleExpression<IndexClass,Expr>::operator()(size_t mlast,
-							      ExtType last) const
+							      ExtType last)
     {
 	//typedef typename IndexClass::RangeType RangeType;
 	const size_t pos = mIndPtr->pos();
@@ -605,7 +605,7 @@ namespace MultiArrayHelper
     }
 
     template <class IndexClass, class Expr>
-    inline void SingleExpression<IndexClass,Expr>::operator()(size_t mlast) const
+    inline void SingleExpression<IndexClass,Expr>::operator()(size_t mlast)
     {
 	//typedef typename IndexClass::RangeType RangeType;
 	const ExtType last;
@@ -672,14 +672,14 @@ namespace MultiArrayHelper
     }
 
     template <class IndexClass, class Expr>
-    inline void SubExpr<IndexClass,Expr>::operator()(size_t mlast, DExt last) const
+    inline void SubExpr<IndexClass,Expr>::operator()(size_t mlast, DExt last)
     {
         operator()(mlast, *reinterpret_cast<ExtType const*>(last.first));
     }
     
     template <class IndexClass, class Expr>
     inline void SubExpr<IndexClass,Expr>::operator()(size_t mlast,
-                                                     ExtType last) const
+                                                     ExtType last)
     {
         const size_t pos = (*mSubSet)[last.val()];
         const size_t mnpos = mlast;
@@ -688,7 +688,7 @@ namespace MultiArrayHelper
     }
 
     template <class IndexClass, class Expr>
-    inline void SubExpr<IndexClass,Expr>::operator()(size_t mlast) const
+    inline void SubExpr<IndexClass,Expr>::operator()(size_t mlast)
     {
 	const ExtType last;
         const size_t pos = (*mSubSet)[last.val()];
@@ -731,12 +731,12 @@ namespace MultiArrayHelper
      *   DynamicExpression   *
      ***************************/
 
-    inline void DynamicExpression::operator()(size_t mlast, DExt last) const
+    inline void DynamicExpression::operator()(size_t mlast, DExt last)
     {
 	(*mNext)(mlast,last);
     }
 
-    inline void DynamicExpression::operator()(size_t mlast) const
+    inline void DynamicExpression::operator()(size_t mlast)
     {
 	(*mNext)(mlast);
     }
@@ -759,13 +759,13 @@ namespace MultiArrayHelper
     ExpressionHolder<Expr>::ExpressionHolder(DynamicExpression expr) : mExpr(expr) {}
 
     template <class Expr>
-    inline void ExpressionHolder<Expr>::operator()(size_t mlast, DExt last) const
+    inline void ExpressionHolder<Expr>::operator()(size_t mlast, DExt last)
     {
 	mExpr(mlast,last);
     }
 
     template <class Expr>
-    inline void ExpressionHolder<Expr>::operator()(size_t mlast, ExtType last) const
+    inline void ExpressionHolder<Expr>::operator()(size_t mlast, ExtType last)
     {
 	mExpr(mlast,
               std::make_pair<size_t const*,size_t>
@@ -774,7 +774,7 @@ namespace MultiArrayHelper
     }
 
     template <class Expr>
-    inline void ExpressionHolder<Expr>::operator()(size_t mlast) const
+    inline void ExpressionHolder<Expr>::operator()(size_t mlast)
     {
 	mExpr(mlast);
     }

From 09f80749104645afc921a8c84026f3e8eb2a35cd Mon Sep 17 00:00:00 2001
From: Christian Zimmermann <christian.zimmermann@ur.de>
Date: Wed, 27 Feb 2019 13:37:53 +0100
Subject: [PATCH 8/9] NOW ALMOST

---
 src/include/multi_array_operation.cc.h |  9 +++++----
 src/include/multi_array_operation.h    | 24 ++++++++++++++++--------
 src/include/type_operations.h          | 10 +++++-----
 src/include/xfor/exttype.h             |  8 ++++----
 4 files changed, 30 insertions(+), 21 deletions(-)

diff --git a/src/include/multi_array_operation.cc.h b/src/include/multi_array_operation.cc.h
index 83e1cc0..cd9bb24 100644
--- a/src/include/multi_array_operation.cc.h
+++ b/src/include/multi_array_operation.cc.h
@@ -251,7 +251,7 @@ namespace MultiArrayTools
 
     template <typename T, class... Ranges>
     template <class ET>
-    inline T ConstOperationRoot<T,Ranges...>::get(ET pos) const
+    inline const T& ConstOperationRoot<T,Ranges...>::get(ET pos) const
     {
 	return mDataPtr[pos.val()/*+mOff*/];
     }
@@ -449,7 +449,7 @@ namespace MultiArrayTools
     
     template <typename T, class... Ranges>
     template <class ET>
-    inline T OperationRoot<T,Ranges...>::get(ET pos) const
+    inline T& OperationRoot<T,Ranges...>::get(ET pos) const
     {
 	return mDataPtr[pos.val()/*+mOff*/];
     }
@@ -504,7 +504,7 @@ namespace MultiArrayTools
 
     template <typename T>
     template <class ET>
-    inline T OperationValue<T>::get(ET pos) const
+    inline const T& OperationValue<T>::get(ET pos) const
     {
 	return mVal;
     }
@@ -595,7 +595,8 @@ namespace MultiArrayTools
     // forward loop !!!!
     template <typename T, class Op, class IndexType>
     template <class ET>
-    inline T Contraction<T,Op,IndexType>::get(ET pos) const
+    inline auto Contraction<T,Op,IndexType>::get(ET pos) const
+        -> decltype(mOp.template get<ET>(pos))
     {
 	return mOp.template get<ET>(pos);
     }
diff --git a/src/include/multi_array_operation.h b/src/include/multi_array_operation.h
index 89846d5..450e183 100644
--- a/src/include/multi_array_operation.h
+++ b/src/include/multi_array_operation.h
@@ -197,6 +197,7 @@ namespace MultiArrayTools
 	typedef ContainerIndex<T,typename Ranges::IndexType...> IndexType;
 
 	static constexpr size_t SIZE = 1;
+        static constexpr bool CONT = true;
 	
 	ConstOperationRoot(const MultiArrayBase<T,Ranges...>& ma,
 			   const std::shared_ptr<typename Ranges::IndexType>&... indices);
@@ -207,7 +208,7 @@ namespace MultiArrayTools
 	ConstOperationRoot(const T* data, const IndexType& ind);
 
 	template <class ET>
-	inline T get(ET pos) const;
+	inline const T& get(ET pos) const;
 
 	template <class ET>
 	inline ConstOperationRoot& set(ET pos);
@@ -243,6 +244,7 @@ namespace MultiArrayTools
 	typedef typename Op::IndexType IndexType;
 
 	static constexpr size_t SIZE = Op::SIZE;
+        static constexpr bool CONT = false;
 	
 	StaticCast(const Op& op);
 
@@ -277,7 +279,8 @@ namespace MultiArrayTools
 	typedef OperationBase<value_type,MetaOperationRoot<Range> > OT;
 	
 	static constexpr size_t SIZE = 1;
-	
+        static constexpr bool CONT = false;
+        
 	MetaOperationRoot(const std::shared_ptr<IndexType>& ind);
 
 	template <class ET>
@@ -310,6 +313,7 @@ namespace MultiArrayTools
 	typedef ContainerIndex<T,typename Ranges::IndexType...> IndexType;
 
 	static constexpr size_t SIZE = 1;
+        static constexpr bool CONT = true;
 	
 	OperationRoot(MutableMultiArrayBase<T,Ranges...>& ma,
 		      const std::shared_ptr<typename Ranges::IndexType>&... indices);
@@ -336,7 +340,7 @@ namespace MultiArrayTools
         OperationRoot& par();
         
 	template <class ET>
-	inline T get(ET pos) const;
+	inline T& get(ET pos) const;
 
 	template <class ET>
 	inline OperationRoot& set(ET pos);
@@ -372,11 +376,12 @@ namespace MultiArrayTools
 	typedef ContainerIndex<T,NullIndex> IndexType;
 
 	static constexpr size_t SIZE = 1;
+        static constexpr bool CONT = true;
 
 	OperationValue(const T& val);
 
 	template <class ET>
-	inline T get(ET pos) const;
+	inline const T& get(ET pos) const;
 
 	template <class ET>
 	inline OperationValue& set(ET pos);
@@ -441,6 +446,7 @@ namespace MultiArrayTools
 	
 	static constexpr size_t SIZE = RootSum<Ops...>::SIZE;
 	static constexpr bool FISSTATIC = OpFunction::FISSTATIC;
+        static constexpr bool CONT = false;
 
     private:
 	std::tuple<Ops...> mOps;
@@ -500,7 +506,6 @@ namespace MultiArrayTools
 	return OpMaker<OpFunction::FISSTATIC>::mkOperation(f, ops...);
     }
 
-
     
     template <typename T, class Op, class IndexType>
     class Contraction : public OperationTemplate<T,Contraction<T,Op,IndexType> >
@@ -511,7 +516,8 @@ namespace MultiArrayTools
 	typedef OperationBase<T,Contraction<T,Op,IndexType> > OT;
 
 	static constexpr size_t SIZE = Op::SIZE;
-	
+        static constexpr bool CONT = Op::CONT;
+        
     private:
 
 	Op mOp;
@@ -523,7 +529,8 @@ namespace MultiArrayTools
 	Contraction(const Op& op, std::shared_ptr<IndexType> ind);
 
 	template <class ET>
-	inline T get(ET pos) const;
+	inline auto get(ET pos) const
+            -> decltype(mOp.template get<ET>(pos));
 
 	template <class ET>
 	inline Contraction& set(ET pos);
@@ -547,7 +554,8 @@ namespace MultiArrayTools
 	typedef OperationTemplate<T,SliceContraction<T,Op,Indices...> > OT;
 
 	static constexpr size_t SIZE = Op::SIZE;
-
+        static constexpr bool CONT = false;
+        
     private:
 
 	mutable Op mOp;
diff --git a/src/include/type_operations.h b/src/include/type_operations.h
index 8f970d9..d40d68e 100644
--- a/src/include/type_operations.h
+++ b/src/include/type_operations.h
@@ -177,17 +177,17 @@ namespace MultiArrayTools
 	}  
     }
 
-    inline v256 operator+(v256 a, v256 b)
+    inline v256 operator+(const v256& a, const v256& b)
     {
 	v256 o;
-	xadd<4>( reinterpret_cast<double*>(&o), reinterpret_cast<const double*>(&a),
-		 reinterpret_cast<const double*>(&b) );
+	xadd<4>( o._x, a._x, b._x );
 	return o;
     }
 
-    inline v256& operator+=(v256& o, v256 a)
+    inline v256& operator+=(v256& o, const v256& a)
     {
-	xsadd<4>( reinterpret_cast<double*>(&o), reinterpret_cast<const double*>(&a) );
+	//xsadd<4>( reinterpret_cast<double*>(&o), reinterpret_cast<const double*>(&a) );
+        xsadd<4>( o._x, a._x );
 	return o;
     }
 /*
diff --git a/src/include/xfor/exttype.h b/src/include/xfor/exttype.h
index 65718b3..ffaccf2 100644
--- a/src/include/xfor/exttype.h
+++ b/src/include/xfor/exttype.h
@@ -36,7 +36,7 @@ namespace MultiArrayHelper
 	template <size_t N>
 	inline MExt(const std::array<size_t,N>& arr);
 
-	inline size_t val() const;
+	inline const size_t& val() const;
 	inline const X& next() const;
 	
 	inline MExt operator+(const MExt& in) const;
@@ -76,7 +76,7 @@ namespace MultiArrayHelper
 	template <size_t N>
 	inline MExt(const std::array<size_t,N>& arr);
 
-	inline size_t val() const;
+	inline const size_t& val() const;
 	inline size_t next() const { return 0; }
 	
 	inline MExt operator+(const MExt& in) const;
@@ -151,7 +151,7 @@ namespace MultiArrayHelper
 	mExt(y.val()), mNext(y.next(), z) {}
     
     template <class X>
-    inline size_t MExt<X>::val() const
+    inline const size_t& MExt<X>::val() const
     {
 	return mExt;
     }
@@ -196,7 +196,7 @@ namespace MultiArrayHelper
 	mExt(std::get<NUM>(arr)) {}
 
     //template <>
-    inline size_t MExt<void>::val() const
+    inline const size_t& MExt<void>::val() const
     {
 	return mExt;
     }

From 25fadae9741bcc2d6e5c815c923996fbac775db0 Mon Sep 17 00:00:00 2001
From: Christian Zimmermann <christian.zimmermann@ur.de>
Date: Wed, 27 Feb 2019 19:02:06 +0100
Subject: [PATCH 9/9] make comparison loop absolutely equivalent to ma test
 loop -> there is no notable difference in the performance anymore

---
 src/tests/op_perf_test.cc | 50 +++++++++++++++++++++------------------
 1 file changed, 27 insertions(+), 23 deletions(-)

diff --git a/src/tests/op_perf_test.cc b/src/tests/op_perf_test.cc
index 5f5de8c..131db7e 100644
--- a/src/tests/op_perf_test.cc
+++ b/src/tests/op_perf_test.cc
@@ -9,6 +9,7 @@
 
 #include <ctime>
 #include <cmath>
+#include <chrono>
 
 #define ONLY_SPIN
 
@@ -135,7 +136,7 @@ namespace {
 	typedef MultiRangeFactory<SR,SR,SR,SR,SR,SR,SR,SR> SR8F;
  	typedef SR8F::oType SR8;
  
-          static const size_t os = 300;
+          static const size_t os = 3000;
           static const size_t is = 65536;
           static const size_t s = is*os;
 	
@@ -176,17 +177,9 @@ namespace {
 	//auto deltap = MAT::getIndex<SR>();
 	auto deltap = MAT::getIndex<GenSingleRange<size_t,SpaceType::NONE,1>>();
 	
-	auto mix = MAT::mkMIndex( alpha, beta, gamma, jj );
+	auto mix = MAT::mkMIndex( jj, alpha, beta, gamma );
 
-	std::clock_t begin = std::clock();
-	//for(size_t i = 0; i != os; ++i){
-        //res1(ii ,delta, deltap).par() += ma(ii, delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix);
-        tcast<v256>(res1)(ii ,delta, deltap) += tcast<v256>(ma)(delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix);
-	//}
-	std::clock_t end = std::clock();
-	std::cout << "MultiArray time: " << static_cast<double>( end - begin ) / CLOCKS_PER_SEC
-		  << std::endl;
-	
+        
 	vector<double> vres(4*4*os);
 	for(size_t d = 0; d != 4; ++d){
 	    for(size_t p = 0; p != 4; ++p){
@@ -194,16 +187,15 @@ namespace {
 		vres[tidx] = 0.;
 	    }
 	}                  
-	std::clock_t begin2 = std::clock();
+        auto begin2 = std::chrono::system_clock::now();
 	double* vrptr = vres.data();
 	double* dptr = data.data();
-	
-        for(size_t j = 0; j != os; ++j) {
-            for(size_t i = 0; i != os; ++i){
-                for(size_t a = 0; a != 4; ++a){
-                    for(size_t b = 0; b != 4; ++b){
-                        for(size_t c = 0; c != 4; ++c){
-                            for(size_t d = 0; d != 4; ++d){
+        for(size_t i = 0; i != os; ++i){
+            for(size_t d = 0; d != 4; ++d){
+                for(size_t j = 0; j != os; ++j) {
+                    for(size_t a = 0; a != 4; ++a){
+                        for(size_t b = 0; b != 4; ++b){
+                            for(size_t c = 0; c != 4; ++c){
 				const size_t tidx = i*4*4 + d*4;
 				const size_t sidx = /*i*65536 +*/ d*4*4*4*4*4*4*4 + a*5*4*4*4*4*4 + b*5*4*4*4 + c*5*4;
 				double* xvrptr = vrptr + tidx;
@@ -218,8 +210,20 @@ namespace {
                 }
             }
 	}
-	std::clock_t end2 = std::clock();
+        auto end2 = std::chrono::system_clock::now();
+	std::cout << "vector - for loop time: " << std::chrono::duration<double>(end2-begin2).count()
+		  << std::endl;
 
+        
+        auto begin = std::chrono::system_clock::now();
+	//for(size_t i = 0; i != os; ++i){
+        //res1(ii ,delta, deltap).par() += ma(ii, delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix);
+        tcast<v256>(res1)(ii ,delta, deltap).par() += tcast<v256>(ma)(delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix);
+	//}
+        auto end = std::chrono::system_clock::now();
+	std::cout << "MultiArray time: " << std::chrono::duration<double>(end-begin).count()
+		  << std::endl;
+	
 	assert( xround(res1.at(mkts(0,0,0))) == xround(vres[0]) );
 	assert( xround(res1.at(mkts(0,0,1))) == xround(vres[1]) );
 	assert( xround(res1.at(mkts(0,0,2))) == xround(vres[2]) );
@@ -240,9 +244,9 @@ namespace {
 	assert( xround(res1.at(mkts(0,3,2))) == xround(vres[14]) );
 	assert( xround(res1.at(mkts(0,3,3))) == xround(vres[15]) );
 
-	std::cout << "vector - for loop time: " << static_cast<double>( end2 - begin2 ) / CLOCKS_PER_SEC
-		  << std::endl;
-	std::cout << "ratio: " << static_cast<double>( end - begin ) / static_cast<double>( end2 - begin2 ) << std::endl;
+	std::cout << "ratio: "
+                  << std::chrono::duration<double>(end-begin).count() / std::chrono::duration<double>(end2-begin2).count()
+                  << std::endl;
     }    
 #ifndef ONLY_SPIN    
     void OpTest_Performance::PCheck()