re-organize cmake / src / build structure (-> linkable)

CMakeLists.txt

@@ -6,6 +6,8 @@ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Werror -std=c++11 -Wpedantic -Wno
 
 enable_testing()
 
+set(INSTALL_PATH ${CMAKE_SOURCE_DIR}/install)
+
 find_package( GTest REQUIRED )
 if(GTest_FOUND)
   include_directories(${GTEST_INCLUDE_DIRS})
@@ -20,26 +22,11 @@ else()
   message(FATAL_ERROR "Threads not found")
 endif()
 
-include_directories(src)
+#include_directories(src)
 
-set(INDEX_CC_FILES "${INDEX_CC_FILES}" "src/ranges/range_base.cc" "src/ranges/index_info.cc")
+#set(INDEX_CC_FILES "${INDEX_CC_FILES}" "src/ranges/range_base.cc" "src/ranges/index_info.cc")
-set(MA_CC_FILES "${MA_CC_FILES}" "${INDEX_CC_FILES}" "src/operation_utils.cc")
+#set(MA_CC_FILES "${MA_CC_FILES}" "${INDEX_CC_FILES}" "src/operation_utils.cc")
-
-add_executable(iutest src/ranges/tests/index_unit_test.cc ${INDEX_CC_FILES})
-target_link_libraries(iutest ${GTEST_BOTH_LIBRARIES} ${CMAKE_THREAD_LIBS_INIT})
-add_test(NAME iutest COMMAND iutest)
-
-add_executable(mautest src/tests/ma_unit_test.cc ${MA_CC_FILES})
-target_link_libraries(mautest ${GTEST_BOTH_LIBRARIES} ${CMAKE_THREAD_LIBS_INIT})
-add_test(NAME mautest COMMAND mautest)
-
-add_executable(oputest src/tests/op_unit_test.cc ${MA_CC_FILES})
-target_link_libraries(oputest ${GTEST_BOTH_LIBRARIES} ${CMAKE_THREAD_LIBS_INIT})
-add_test(NAME oputest COMMAND oputest)
-
-add_executable(opptest src/tests/op_perf_test.cc ${MA_CC_FILES})
-#target_link_libraries(oputest)
-
-#install(TARGETS testm DESTINATION install)
 
 set(CMAKE_INSTALL_PREFIX ..)
+
+add_subdirectory(src)

install/include/multi_array_operation.h

@@ -13,7 +13,7 @@
 #include "mbase_def.h"
 
 //#include "block/block.h"
-#include "operation_utils.h"
+//#include "operation_utils.h"
 #include "ranges/rheader.h"
 #include "pack_num.h"
 
@@ -96,8 +96,7 @@ namespace MultiArrayTools
 	typedef OperationBase<T> OB;
 	typedef ContainerRange<Ranges...> CRange;
 	typedef typename MultiRange<Ranges...>::IndexType IndexType;
-	typedef MBlock<T> bType;	
+	
-
 	OperationMaster(MutableMultiArrayBase<T,Ranges...>& ma, const OpClass& second,
 			std::shared_ptr<typename CRange::IndexType>& index);
 
@@ -118,7 +117,6 @@ namespace MultiArrayTools
 	T* mDataPtr;
 	std::shared_ptr<IndexType> mIndex;
 	IndexInfo mIInfo;
-	mutable bType mBlock;
     };
 
     

install/include/pack_num.h

@@ -17,33 +17,6 @@ namespace MultiArrayHelper
     template <size_t N>
     struct PackNum
     {
-	
-	template <class... Ops>
-	static void makeBlockTypeVec(std::vector<std::pair<BlockType,size_t> >& btv,
-				     const std::tuple<Ops...>& ops,
-				     std::shared_ptr<VIWB> idxPtr, bool init)
-	{
-	    auto subvec = std::move( std::get<N>(ops).block(idxPtr, init) );
-	    btv.insert(btv.end(), subvec.begin(), subvec.end() );
-	    PackNum<N-1>::makeBlockTypeVec(btv, ops, idxPtr, init);
-	}
-
-	template <class... Ops>
-	static void makeBlockTypeVec(std::vector<std::pair<BlockType,size_t> >& btv,
-				     const std::tuple<Ops...>& ops,
-				     const IndexInfo* idxPtr, bool init)
-	{
-	    auto subvec = std::move( std::get<N>(ops).block(idxPtr, init) );
-	    btv.insert(btv.end(), subvec.begin(), subvec.end() );
-	    PackNum<N-1>::makeBlockTypeVec(btv, ops, idxPtr, init);
-	}
-	
-	template <typename T, class Func, class ArgTuple, class... Args>
-	static void unpackArgs(BlockResult<T>& res, const ArgTuple& tp, const Args&... args)
-	{
-	    PackNum<N-1>::template unpackArgs<T,Func>(res, tp, std::get<N>(tp).get(), args...);
-	}
-
 	template <typename... T>
 	static void printTuple(std::ostream& out, const std::tuple<T...>& tp)
 	{
@@ -88,34 +61,6 @@ namespace MultiArrayHelper
     template<>
     struct PackNum<0>
     {
-
-	template <typename T, class Func, class ArgTuple, class... Args>
-	static void unpackArgs(BlockResult<T>& res, const ArgTuple& tp, const Args&... args)
-	{
-	    static_assert(sizeof...(Args) == std::tuple_size<ArgTuple>::value-1,
-			  "inconsistent number of arguments");
-	    BlockBinaryOp<T,Func,decltype(std::get<0>(tp).get()), decltype(args)...> f(res);
-	    f(std::get<0>(tp).get(), args...);
-	}
-	
-	template <class... Ops>
-	static void makeBlockTypeVec(std::vector<std::pair<BlockType,size_t> >& btv,
-				     const std::tuple<Ops...>& ops,
-				     std::shared_ptr<VIWB> idxPtr, bool init)
-	{
-	    auto subvec = std::move( std::get<0>(ops).block(idxPtr, init) );
-	    btv.insert(btv.end(), subvec.begin(), subvec.end() );
-	}
-
-	template <class... Ops>
-	static void makeBlockTypeVec(std::vector<std::pair<BlockType,size_t> >& btv,
-				     const std::tuple<Ops...>& ops,
-				     const IndexInfo* idxPtr, bool init)
-	{
-	    auto subvec = std::move( std::get<0>(ops).block(idxPtr, init) );
-	    btv.insert(btv.end(), subvec.begin(), subvec.end() );
-	}
-
 	template <typename... T>
 	static void printTuple(std::ostream& out, const std::tuple<T...>& tp)
 	{	    

install/include/ranges/range_base.h

@@ -60,6 +60,11 @@ namespace MultiArrayTools
 
 	bool operator==(const RangeBase& in) const;
 	bool operator!=(const RangeBase& in) const;
+
+	//virtual bool regular() const = 0; // integer distance (e.g. 2,3,4,...)
+	//virtual bool linear() const = 0; // 1dim valuable (e.g. 2.45, 3.12, 3.56,...)
+	//virtual bool multi() const = 0; // mdim
+	//virtual bool maplike() const = 0; // meta type is ~ MultiArray<T,...>
 	
 	friend RangeFactoryBase;
 	

src/CMakeLists.txt

@@ -0,0 +1,6 @@
+
+include_directories(${CMAKE_SOURCE_DIR}/src/include)
+add_subdirectory(tests)
+add_subdirectory(lib)
+
+install(DIRECTORY include/ DESTINATION ${INSTALL_PATH}/include)

src/block/bbase_def.h

@@ -1,45 +0,0 @@
-
-
-#ifndef __bbase_def_h__
-#define __bbase_def_h__
-
-namespace MultiArrayHelper
-{
-
-    template <typename T, class OpFunc, class BlockClass1, class BlockClass2>
-    class BlockBinaryOp;
-    
-    
-    template <typename T, class BlockClass>
-    class BlockBase;
-
-    
-    template <typename T, class BlockClass>
-    class MutableBlockBase;
-
-    
-    template <typename T>
-    class Block;
-
-    
-    template <typename T>
-    class MBlock;
-
-    
-    template <typename T>
-    class BlockResult;
-
-    
-    enum class BlockType {
-	INDEF = 0,
-	BLOCK = 1,
-	VALUE = 2,
-	SPLIT = 3,
-	RESULT = 4,
-	ARRAY = 5
-    };
-
-
-} // end namespace MultiArrayHelper
-
-#endif

src/block/block.h

@@ -1,606 +0,0 @@
-// -*- C++ -*-
-
-#ifndef __ma_block_h__
-#define __ma_block_h__
-
-#include <cstdlib>
-#include <vector>
-
-#include "base_def.h"
-#include "block/bbase_def.h"
-
-namespace MultiArrayHelper
-{
-    
-    // manage vectorization in the future !!
-    
-    template <typename T, class OpFunc, class BlockClass1, class BlockClass2>
-    class BlockBinaryOp
-    {
-    public:
-	BlockBinaryOp(BlockResult<T>& mRes);
-	BlockResult<T>& operator()(const BlockClass1& arg1, const BlockClass2& arg2);
-    private:
-	BlockResult<T>& mRes;
-    };
-
-    template <typename T, class OpFunc, class BlockClass>
-    class BlockBinaryOpSelf
-    {
-    public:
-	BlockBinaryOpSelf(BlockResult<T>& res);
-	BlockResult<T>& operator()(const BlockClass& arg);
-    private:
-	BlockResult<T>& mRes;
-    };
-
-    template <typename T, class OpFunc, class BlockClass>
-    class BlockContraction
-    {
-    public:
-	BlockContraction(T& res);
-	T& operator()(const BlockClass& arg);
-    private:
-	T& mRes;
-    };
-    
-    // EVERYTHING IN HERE MUST  N O T  BE VITUAL !!
-    
-    template <typename T, class BlockClass>
-    class BlockBase
-    {
-    public:
-
-	const BlockClass& THIS() const { return static_cast<BlockClass const&>(*this); }
-	BlockClass& THIS() { return static_cast<BlockClass&>(*this); }
-
-	static BlockType sType() { return BlockClass::sType(); }
-	
-	size_t size() const;
-	bool init() const;
-	
-	BlockType type() const { return THIS().type(); }
-	const T& operator[](size_t pos) const { return THIS()[pos]; }
-	BlockClass& set(size_t npos) { return THIS().set(npos); }
-	size_t stepSize() const { return THIS().stepSize(); }
-	
-    private:
-
-	friend BlockClass;
-	friend MutableBlockBase<T,BlockClass>;
-	
-	DEFAULT_MEMBERS(BlockBase);
-	BlockBase(size_t size);
-	
-	size_t mSize = 0;
-	bool mInit = false;
-    };
-
-    template <typename T>
-    std::ostream& operator<<(std::ostream& out, const Block<T>& block)
-    {
-	out << block[0];
-	for(size_t i = 1; i != block.size(); ++i){
-	    out << ", " << block[i];
-	}
-	return out;
-    }
-    
-    template <typename T>
-    std::ostream& operator<<(std::ostream& out, const MBlock<T>& block)
-    {
-	out << block[0];
-	for(size_t i = 1; i != block.size(); ++i){
-	    out << ", " << block[i];
-	}
-	return out;
-    }
-
-    template <typename T>
-    std::ostream& operator<<(std::ostream& out, const BlockResult<T>& block)
-    {
-	out << block[0];
-	for(size_t i = 1; i != block.size(); ++i){
-	    out << ", " << block[i];
-	}
-	return out;
-    }
-
-    
-    template <typename T, class BlockClass>
-    class MutableBlockBase : public BlockBase<T,BlockClass>
-    {
-    public:
-
-	typedef BlockBase<T,BlockClass> BB;
-	
-	T& operator[](size_t pos) { return BB::THIS()[pos]; }
-	
-    private:
-
-	friend BlockClass;
-	
-	DEFAULT_MEMBERS(MutableBlockBase);
-	MutableBlockBase(size_t size);
-
-    };
-    
-    template <typename T>
-    class Block : public BlockBase<T,Block<T> >
-    {
-    public:
-	typedef BlockBase<T,Block<T> > BB;
-
-	static BlockType sType() { return BlockType::BLOCK; }
-	
-	DEFAULT_MEMBERS(Block);
-	Block(const T* data, size_t begPos, size_t size, size_t stepSize);
-
-	BlockType type() const;
-	const T& operator[](size_t pos) const;
-	Block& set(size_t npos);
-	size_t stepSize() const;
-	
-    protected:
-	const T* mData;
-	const T* mBegPtr;
-	size_t mStepSize;
-    };
-
-    template <class BlockClass>
-    class BlockArray : public BlockBase<BlockClass,BlockArray<BlockClass> >
-    {
-
-	typedef BlockBase<BlockClass,BlockArray<BlockClass> > BB;
-
-	static BlockType sType() { return BlockType::ARRAY; }
-	
-	DEFAULT_MEMBERS(BlockArray);
-
-	BlockArray(const BlockClass& block, size_t size, size_t stepSize);
-
-	BlockType type() const;
-	BlockClass& operator[](size_t pos);
-	BlockArray& set(size_t npos);
-	size_t stepSize() const;
-	
-    protected:
-	const BlockClass& mBlock;
-	size_t mStepSize; // total stepSize !!
-    };
-    
-    
-    template <typename T>
-    class MBlock : public MutableBlockBase<T,MBlock<T> >
-    {
-    public:
-
-	typedef BlockBase<T,MBlock<T> > BB;
-
-	static BlockType sType() { return BlockType::BLOCK; }
-	
-	DEFAULT_MEMBERS(MBlock);
-	MBlock(T* data, size_t begPos, size_t size, size_t stepSize);
-
-	template <class BlockClass>
-	MBlock& operator=(const BlockClass& in);
-	
-	BlockType type() const;
-	const T& operator[](size_t pos) const;
-	T& operator[](size_t pos);
-	MBlock& set(size_t npos);
-	size_t stepSize() const;
-	
-    protected:
-	T* mData;
-	T* mBegPtr;
-	size_t mStepSize;
-    };    
-    
-    template <typename T>
-    class BlockResult : public MutableBlockBase<T,BlockResult<T> >
-    {
-    public:
-	
-	typedef BlockBase<T,BlockResult<T> > BB;
-	using BB::init;
-
-	static BlockType sType() { return BlockType::RESULT; }
-	
-	BlockResult();
-	BlockResult(const BlockResult& in);
-	BlockResult(BlockResult&& in);
-	BlockResult& operator=(const BlockResult& in);
-	BlockResult& operator=(BlockResult&& in);
-
-	template <typename... ArgTypes>
-	BlockResult(size_t size, const ArgTypes&... args);
-	
-	~BlockResult();
-	
-	template <class BlockClass>
-	BlockResult& operator=(const BlockClass& in);
-
-	template <typename... Args>
-	BlockResult& init(size_t size, const Args&... args);
-	
-	BlockType type() const;
-	const T& operator[](size_t pos) const;
-	T& operator[](size_t i);
-	BlockResult& set(size_t npos);
-	size_t stepSize() const;
-	
-    protected:
-	T* mResPtr = nullptr;
-	T* mBegPtr = nullptr;
-    };
-    
-} // end namespace MultiArrayHelper
-
-/* ========================= *
- * ---   TEMPLATE CODE   --- *
- * ========================= */
-
-namespace MultiArrayHelper
-{
-
-    /*********************
-     *   BlockBinaryOp   *
-     *********************/
-
-    template <typename T, class OpFunc, class BlockClass1, class BlockClass2>
-    BlockBinaryOp<T,OpFunc,BlockClass1,BlockClass2>::BlockBinaryOp(BlockResult<T>& res) : mRes(res) {}
-    
-    template <typename T, class OpFunc, class BlockClass1, class BlockClass2>
-    BlockResult<T>& BlockBinaryOp<T,OpFunc,BlockClass1,BlockClass2>::operator()(const BlockClass1& arg1,
-								     const BlockClass2& arg2)
-    {
-	static OpFunc f;
-	//assert(mRes.init() and arg1.init() and arg2.init());
-	//assert(arg1.size() == arg2.size());
-	for(size_t i = 0; i != arg1.size(); ++i){
-	    mRes[i] = f(arg1[i], arg2[i]);
-	}
-	return mRes;
-    }
-
-    template <typename T, class OpFunc, class BlockClass>
-    BlockBinaryOpSelf<T,OpFunc,BlockClass>::BlockBinaryOpSelf(BlockResult<T>& res) : mRes(res) {}
-
-    template <typename T, class OpFunc, class BlockClass>
-    BlockResult<T>& BlockBinaryOpSelf<T,OpFunc,BlockClass>::operator()(const BlockClass& arg)
-    {
-	static OpFunc f;
-	//assert(mRes.init() and arg.init());
-	//assert(mRes.size() == arg.size());
-	for(size_t i = 0; i != arg.size(); ++i){
-	    mRes[i] = f(mRes[i], arg[i]);
-	}
-	return mRes;
-    }
-
-    template <typename T, class OpFunc, class BlockClass>
-    BlockContraction<T,OpFunc,BlockClass>::BlockContraction(T& res) : mRes(res) {}
-
-    template <typename T, class OpFunc, class BlockClass>
-    T& BlockContraction<T,OpFunc,BlockClass>::operator()(const BlockClass& arg)
-    {
-	static OpFunc f;
-
-	for(size_t i = 0; i != arg.size(); ++i){
-	    mRes = f(mRes, arg[i]);
-	}
-	return mRes;
-    }
-
-    
-    /*****************
-     *   BlockBase   *
-     *****************/
-    
-    template <typename T, class BlockClass>
-    BlockBase<T,BlockClass>::BlockBase(size_t size) : mSize(size), mInit(size != 0) {}
-
-    template <typename T, class BlockClass>
-    size_t BlockBase<T,BlockClass>::size() const
-    {
-	return mSize;
-    }
-
-    template <typename T, class BlockClass>
-    bool BlockBase<T,BlockClass>::init() const
-    {
-	return mInit;
-    }
-    
-    /************************
-     *   MutableBlockBase   *
-     ************************/
-    
-    template <typename T, class BlockClass>
-    MutableBlockBase<T,BlockClass>::MutableBlockBase(size_t size) : BlockBase<T,BlockClass>(size) {}
-
-    
-    /*************
-     *   Block   *
-     *************/
-
-    template <typename T>
-    Block<T>::Block(const T* data,
-		    size_t begPos, size_t size, size_t stepSize) :
-	BlockBase<T,Block>(size),
-	mData(data),
-	mBegPtr(data + begPos),
-	mStepSize(stepSize) {}
-
-    template <typename T>
-    BlockType Block<T>::type() const
-    {
-	return mStepSize == 0 ? BlockType::VALUE :
-	    ( mStepSize == 1 ? BlockType::BLOCK : BlockType::SPLIT );
-    }
-    
-    template <typename T>
-    const T& Block<T>::operator[](size_t i) const
-    {
-	
-	return *(mBegPtr + i * mStepSize);
-    }
-
-    template <typename T>
-    Block<T>& Block<T>::set(size_t npos)
-    {
-	mBegPtr = mData + npos;
-	return *this;
-    }
-
-    template <typename T>
-    size_t Block<T>::stepSize() const
-    {
-	return 1;
-    }
-
-    
-    /******************
-     *   BlockArray   *
-     ******************/
-    
-    template <class BlockClass>
-    BlockArray<BlockClass>::BlockArray(const BlockClass& block, size_t size, size_t stepSize) :
-	BlockBase<BlockClass,BlockArray<BlockClass> >(size),
-	mBlock(block), mStepSize(stepSize) {}
-
-    template <class BlockClass>
-    BlockType BlockArray<BlockClass>::type() const
-    {
-	return BlockType::ARRAY;
-    }
-
-    template <class BlockClass>
-    BlockClass& BlockArray<BlockClass>::operator[](size_t pos)
-    {
-	return mBlock.set(mStepSize * pos);
-    }
-
-    template <class BlockClass>
-    BlockArray<BlockClass>& BlockArray<BlockClass>::set(size_t npos)
-    {
-	mBlock.set(npos);
-	return *this;
-    }
-
-    template <class BlockClass>
-    size_t BlockArray<BlockClass>::stepSize() const
-    {
-	return mStepSize;
-    }
-	
-    
-    /**************
-     *   MBlock   *
-     **************/
-
-    template <typename T>
-    MBlock<T>::MBlock(T* data,
-		      size_t begPos, size_t size, size_t stepSize) :
-	MutableBlockBase<T,MBlock>(size),
-	mData(data),
-	mBegPtr(data + begPos),
-	mStepSize(stepSize) {}
-
-    template <typename T>
-    template <class BlockClass>
-    MBlock<T>& MBlock<T>::operator=(const BlockClass& in)
-    {
-	for(size_t i = 0; i != BB::mSize; ++i){
-	    (*this)[i] = in[i];
-	}
-	return *this;
-    }
-   
-    template <typename T>
-    BlockType MBlock<T>::type() const
-    {
-	return mStepSize == 0 ? BlockType::VALUE :
-	    ( mStepSize == 1 ? BlockType::BLOCK : BlockType::SPLIT );
-    }
-    
-    template <typename T>
-    const T& MBlock<T>::operator[](size_t i) const
-    {
-	
-	return *(mBegPtr + i * mStepSize);
-    }
-
-    template <typename T>
-    T& MBlock<T>::operator[](size_t i)
-    {
-	
-	return *(mBegPtr + i * mStepSize);
-    }
-    
-    template <typename T>
-    MBlock<T>& MBlock<T>::set(size_t npos)
-    {
-	mBegPtr = mData + npos;
-	return *this;
-    }   
-
-    template <typename T>
-    size_t MBlock<T>::stepSize() const
-    {
-	return 1;
-    }
-    
-    /*******************
-     *   BlockResult   *
-     *******************/
-
-    template <typename T>
-    BlockResult<T>::BlockResult() : MutableBlockBase<T,BlockResult>() {}
-
-    template <typename T>
-    BlockResult<T>::BlockResult(const BlockResult<T>& in) : MutableBlockBase<T,BlockResult>(in.size())
-    {
-	if(BB::mInit){
-	    mResPtr = new T[BB::mSize];
-	}
-	for(size_t i = 0; i != BB::mSize; ++i){
-	    mResPtr[i] = in.mResPtr[i];
-	}
-	mBegPtr = mResPtr;
-    }
-
-    template <typename T>
-    BlockResult<T>::BlockResult(BlockResult<T>&& in) : MutableBlockBase<T,BlockResult>(in.size())
-    {
-	if(BB::mInit){
-	    mResPtr = in.mResPtr;
-	    mBegPtr = mResPtr;
-	}
-	in.mSize = 0;
-	in.mInit = false;
-	in.mResPtr = nullptr;
-	in.mBegPtr = nullptr;
-    }
-
-    template <typename T>
-    BlockResult<T>& BlockResult<T>::operator=(const BlockResult<T>& in)
-    {
-	BB::mSize = in.size();
-	BB::mInit = BB::mInit and BB::mSize != 0;
-	if(BB::mInit){
-	    mResPtr = new T[BB::mSize];
-	}
-	for(size_t i = 0; i != BB::mSize; ++i){
-	    mResPtr[i] = in.mResPtr[i];
-	}
-	mBegPtr = mResPtr;
-	return *this;
-    }
-
-    template <typename T>
-    BlockResult<T>& BlockResult<T>::operator=(BlockResult<T>&& in)
-    {
-	BB::mSize = in.size();
-	BB::mInit = BB::mInit and BB::mSize != 0;
-	if(BB::mInit){
-	    mResPtr = in.mResPtr;
-	    mBegPtr = mResPtr;
-	}
-	in.mSize = 0;
-	in.mInit = false;
-	in.mResPtr = nullptr;
-	return *this;
-    }
-
-    template <typename T>
-    template <typename... ArgTypes>
-    BlockResult<T>::BlockResult(size_t size, const ArgTypes&... args) :
-	MutableBlockBase<T,BlockResult>(size)
-    {
-	if(BB::mInit){
-	    mResPtr = new T[BB::mSize](args...);
-	}
-	for(size_t i = 0; i != BB::mSize; ++i){
-	    mResPtr[i] = static_cast<T>( 0 );
-	}
-	mBegPtr = mResPtr;
-    }
-
-    template <typename T>
-    BlockResult<T>::~BlockResult()
-    {
-	delete[] mResPtr;
-	mResPtr = nullptr;
-	mBegPtr = nullptr;
-	BB::mInit = false;
-	BB::mSize = 0;
-    }
-    
-    template <typename T>
-    template <class BlockClass>
-    BlockResult<T>& BlockResult<T>::operator=(const BlockClass& in)
-    {
-	assert(BB::mInit);
-	assert(BB::mSize == in.size());
-	for(size_t i = 0; i != BB::mSize; ++i){
-	    (*this)[i] = in[i];
-	}
-	return *this;
-    }
-    
-    template <typename T>
-    BlockType BlockResult<T>::type() const
-    {
-	return BlockType::RESULT;
-    }
-    
-    template <typename T>
-    const T& BlockResult<T>::operator[](size_t i) const
-    {
-	return mBegPtr[i];
-    }
-
-    template <typename T>
-    T& BlockResult<T>::operator[](size_t i)
-    {
-	return mBegPtr[i];
-    }
-
-    template <typename T>
-    BlockResult<T>& BlockResult<T>::set(size_t npos)
-    {
-	mBegPtr = mResPtr + npos;
-	return *this;
-    }
-
-    template <typename T>
-    size_t BlockResult<T>::stepSize() const
-    {
-	return 1;
-    }
-    
-    template <typename T>
-    template <typename... Args>
-    BlockResult<T>& BlockResult<T>::init(size_t size, const Args&... args)
-    {
-	BB::mSize = size;
-	delete[] mResPtr;
-	if(BB::mSize != 0){
-	    BB::mInit = true;
-	    mResPtr = new T[BB::mSize](args...);
-	}
-	else {
-	    BB::mInit = false;
-	}
-	for(size_t i = 0; i != BB::mSize; ++i){
-	    mResPtr[i] = static_cast<T>( 0 );
-	}
-	mBegPtr = mResPtr;
-	return *this;
-    }
-    
-} // end namespace MultiArrayHelper
-
-#endif

src/include/arith.h

@@ -0,0 +1,46 @@
+
+#ifndef __arith_h__
+#define __arith_h__
+
+namespace MultiArrayHelper
+{
+    // OPERATIONS (STATIC)
+    template <typename T>
+    struct plus
+    {
+	static inline T apply(T a1, T a2)
+	{
+	    return a1 + a2;
+	}
+    };
+
+    template <typename T>
+    struct minus
+    {
+	static inline T apply(T a1, T a2)
+	{
+	    return a1 - a2;
+	}
+    };
+
+    template <typename T>
+    struct multiplies
+    {
+	static inline T apply(T a1, T a2)
+	{
+	    return a1 * a2;
+	}
+    };
+
+    template <typename T>
+    struct divides
+    {
+	static inline T apply(T a1, T a2)
+	{
+	    return a1 / a2;
+	}
+    };
+
+} // end namespace MultiArrayHelper
+    
+#endif

src/include/base_def.h

@@ -0,0 +1,35 @@
+// -*- C++ -*-
+
+#ifndef __base_def_h__
+#define __base_def_h__
+
+#include <cassert>
+
+#define DEBUG_MODE_X
+
+#ifdef DEBUG_MODE_X
+
+#include <iostream>
+#define CHECK std::cout << __FILE__ << ": @" << __LINE__ << " in " << __func__ << std::endl;
+#define VCHECK(a) std::cout << __FILE__ << ": @" << __LINE__ \
+  << " in " << __func__ << ": " << #a << " = " << a << std::endl;
+
+#else
+#define CHECK
+#define VCHECK(a)
+
+#endif
+
+#define DEFAULT_MEMBERS(__class_name__) __class_name__() = default; \
+    __class_name__(const __class_name__& in) = default; \
+    __class_name__& operator=(const __class_name__& in) = default; \
+    __class_name__(__class_name__&& in) = default; \
+    __class_name__& operator=(__class_name__&& in) = default
+
+#define DEFAULT_MEMBERS_X(__class_name__) __class_name__(const __class_name__& in) = default; \
+    __class_name__& operator=(const __class_name__& in) = default; \
+    __class_name__(__class_name__&& in) = default; \
+    __class_name__& operator=(__class_name__&& in) = default
+
+
+#endif

src/include/helper_tools.h

@@ -0,0 +1,72 @@
+
+#ifndef __helper_tools_h__
+#define __helper_tools_h__
+
+#include "base_def.h"
+
+namespace MultiArrayTools
+{
+
+    template <class RangeType>
+    auto getIndex(std::shared_ptr<RangeType> range)
+	-> std::shared_ptr<typename RangeType::IndexType>;
+
+    // only if 'RangeType' is defaultable and unique (Singleton)
+    template <class RangeType>
+    auto getIndex() -> std::shared_ptr<typename RangeType::IndexType>;
+
+    template <class... RangeTypes>
+    auto mkMulti(std::shared_ptr<RangeTypes>... ranges)
+	-> std::shared_ptr<MultiRange<RangeTypes...> >;
+
+    template <class... IndexTypes>
+    auto mkMIndex(std::shared_ptr<IndexTypes>... indices)
+	-> decltype( getIndex( mkMulti( indices.range()... ) ) );
+
+    
+}
+
+/* ========================= *
+ * ---   TEMPLATE CODE   --- *
+ * ========================= */
+
+namespace MultiArrayTools
+{
+    template <class RangeType>
+    auto getIndex(std::shared_ptr<RangeType> range)
+	-> std::shared_ptr<typename RangeType::IndexType>
+    {
+	return std::dynamic_pointer_cast<IndexWrapper<typename RangeType::IndexType> >
+	    ( range->index() )->get();
+    }
+
+    template <class RangeType>
+    auto getIndex() -> std::shared_ptr<typename RangeType::IndexType>
+    {
+	static_assert( RangeType::defaultable,
+		       /*typeid(typename RangeType).name() + */" is not defaultable" );
+	static auto f = RangeType::factory();
+	static auto r = std::dynamic_pointer_cast<RangeType>( f.create() );
+	return std::dynamic_pointer_cast<IndexWrapper<typename RangeType::IndexType> >
+	    ( r->index() )->get();
+    }
+    
+    template <class... RangeTypes>
+    auto mkMulti(std::shared_ptr<RangeTypes>... ranges)
+	-> std::shared_ptr<MultiRange<RangeTypes...> >
+    {
+	MultiRangeFactory<RangeTypes...> mrf( ranges... );
+	return std::dynamic_pointer_cast<MultiRange<RangeTypes...> >( mrf.create() );
+    }
+
+    template <class... IndexTypes>
+    auto mkMIndex(std::shared_ptr<IndexTypes>... indices)
+	-> decltype( getIndex( mkMulti( indices->range()... ) ) )
+    {
+	auto mi = getIndex( mkMulti( indices->range()... ) );	    
+	(*mi)( indices... );
+	return mi;
+    }
+}
+
+#endif

src/include/mbase_def.h

@@ -0,0 +1,56 @@
+
+#ifndef __mbase_def_h__
+#define __mbase_def_h__
+
+namespace MultiArrayTools
+{
+    /***********************
+     *  Provided Types     *
+     ***********************/
+    
+    // multi_array.h
+    template <typename T, class... SRanges>
+    class MultiArrayBase;
+
+    // multi_array.h
+    template <typename T, class... SRanges>
+    class MutableMultiArrayBase;
+    
+    // multi_array.h
+    template <typename T, class... SRanges>
+    class MultiArray;
+   
+    // multi_array_operation.h
+    template <typename T>
+    class OperationBase;
+
+    // multi_array_operation.h
+    template <typename T>
+    class MutableOperationBase;
+
+    // multi_array_operation.h
+    template <typename T, class OperationClass>
+    class OperationTemplate;
+    
+    // multi_array_operation.h
+    template <typename T, class OpClass, class... Ranges>
+    class OperationMaster;
+    
+    // multi_array_operation.h
+    template <typename T, class... Ranges>
+    class OperationRoot;
+
+    // multi_array_operation.h
+    template <typename T, class... Ranges>
+    class ConstOperationRoot;
+
+    // multi_array_operation.h
+    template <typename T, class OpFunction, class... Ops>
+    class Operation;
+
+    // multi_array_operation.h
+    template <typename T, class Op, class IndexType>
+    class Contraction;
+}
+
+#endif

src/include/multi_array.h

@@ -0,0 +1,872 @@
+// -*- C++ -*-
+
+#ifndef __multi_array_h__
+#define __multi_array_h__
+
+#include <cstdlib>
+#include <vector>
+#include <memory>
+#include <iterator>
+#include <algorithm>
+
+#include "base_def.h"
+#include "mbase_def.h"
+
+#include "ranges/rheader.h"
+
+namespace MultiArrayTools
+{
+
+    // Explicitely specify subranges in template argument !!!
+    template <typename T, class... SRanges>
+    class MultiArrayBase
+    {
+    public:
+
+	typedef T value_type;
+	typedef ContainerRange<SRanges...> CRange;
+	typedef typename CRange::IndexType IndexType;
+	
+	class const_iterator : public std::iterator<std::random_access_iterator_tag,T>
+	{
+	public:
+
+	    DEFAULT_MEMBERS(const_iterator);
+	    
+	    const_iterator(const MultiArrayBase& ma);
+	    const_iterator(const MultiArrayBase& ma, const typename CRange::IndexType& index);
+	    
+	    // Requirements:
+	    bool operator==(const const_iterator& it) const;
+	    bool operator!=(const const_iterator& it) const;
+
+	    const T& operator*() const;
+	    T const* operator->() const;
+
+	    const_iterator& operator++();
+	    const_iterator operator++(int);
+	    const_iterator& operator--();
+	    const_iterator operator--(int);
+
+	    const_iterator& operator+=(int diff);
+	    const_iterator& operator-=(int diff);
+	    const_iterator operator+(int num) const;
+	    const_iterator operator-(int num) const;
+
+	    int operator-(const const_iterator& it) const;
+
+	    const T& operator[](int num) const;
+
+	    bool operator<(const const_iterator& it) const;
+	    bool operator>(const const_iterator& it) const;
+	    bool operator<=(const const_iterator& it) const;
+	    bool operator>=(const const_iterator& it) const;
+
+	    // Multi Array specific:
+	    typename ContainerRange<SRanges...>::IndexType index() const;
+	    
+	protected:
+	    MultiArrayBase const* mMAPtr = nullptr;
+	    size_t mPos;
+	};
+
+	DEFAULT_MEMBERS(MultiArrayBase);
+	MultiArrayBase(const std::shared_ptr<SRanges>&... ranges);
+
+	virtual ~MultiArrayBase() = default;
+
+	virtual const T& operator[](const IndexType& i) const = 0;
+	//virtual const T& operator[](const typename CRange::IndexType& i) const = 0;
+	virtual const T& at(const typename CRange::IndexType::MetaType& meta) const = 0;
+
+	virtual const T* data() const = 0;
+	virtual const std::vector<T>& datav() const = 0;
+	
+	virtual size_t size() const; 
+	virtual bool isSlice() const = 0;
+
+	virtual const_iterator begin() const;
+	virtual const_iterator end() const;
+	
+	virtual IndexType beginIndex() const;
+	virtual IndexType endIndex() const;
+
+	virtual const std::shared_ptr<CRange>& range() const;
+
+	virtual bool isConst() const;
+
+	virtual ConstOperationRoot<T,SRanges...>
+	operator()(std::shared_ptr<typename SRanges::IndexType>&... inds) const;
+	
+	virtual bool isInit() const;
+
+    protected:
+	bool mInit = false;
+	std::shared_ptr<CRange> mRange;
+
+    };
+
+    template <typename T, class... SRanges>
+    class MutableMultiArrayBase : public MultiArrayBase<T,SRanges...>
+    {
+    public:
+
+	typedef ContainerRange<SRanges...> CRange;
+	typedef typename MultiArrayBase<T,SRanges...>::const_iterator const_iterator;
+	typedef MultiArrayBase<T,SRanges...> MAB;
+	typedef typename CRange::IndexType IndexType;
+	
+	class iterator : public std::iterator<std::random_access_iterator_tag,T>,
+			 public std::iterator<std::output_iterator_tag,T>
+	{
+	public:
+
+	    DEFAULT_MEMBERS(iterator);
+	    
+	    iterator(MutableMultiArrayBase& ma);
+	    iterator(MutableMultiArrayBase& ma, const IndexType& index);
+	    
+	    // Requirements:
+	    bool operator==(const iterator& it) const;
+	    bool operator!=(const iterator& it) const;
+
+	    const T& operator*() const;
+	    T const* operator->() const;
+	    T& operator*();
+	    T* operator->();
+
+	    iterator& operator++();
+	    iterator operator++(int);
+	    iterator& operator--();
+	    iterator operator--(int);
+
+	    iterator& operator+=(int diff);
+	    iterator& operator-=(int diff);
+	    iterator operator+(int num) const;
+	    iterator operator-(int num) const;
+
+	    int operator-(const iterator& it) const;
+
+	    const T& operator[](int num) const;
+	    T& operator[](int num);
+
+	    bool operator<(const iterator& it) const;
+	    bool operator>(const iterator& it) const;
+	    bool operator<=(const iterator& it) const;
+	    bool operator>=(const iterator& it) const;
+
+	    // Multi Array specific:
+	    typename CRange::IndexType index() const;
+	    	    
+	protected:
+	    MutableMultiArrayBase* mMAPtr = nullptr;
+	    size_t mPos;
+	};
+
+	using MultiArrayBase<T,SRanges...>::operator[];
+	using MultiArrayBase<T,SRanges...>::at;
+	using MultiArrayBase<T,SRanges...>::data;
+	using MultiArrayBase<T,SRanges...>::datav;
+	using MultiArrayBase<T,SRanges...>::begin;
+	using MultiArrayBase<T,SRanges...>::end;
+	
+	DEFAULT_MEMBERS(MutableMultiArrayBase);
+	MutableMultiArrayBase(const std::shared_ptr<SRanges>&... ranges);
+
+	virtual T& operator[](const IndexType& i) = 0;
+	virtual T& at(const typename CRange::IndexType::MetaType& meta) = 0;
+	
+	virtual T* data() = 0;
+	virtual std::vector<T>& datav() = 0;
+	
+	virtual iterator begin();
+	virtual iterator end();
+
+	virtual bool isConst() const override;
+
+	virtual ConstOperationRoot<T,SRanges...>
+	operator()(std::shared_ptr<typename SRanges::IndexType>&... inds) const override;
+	virtual OperationRoot<T,SRanges...> operator()(std::shared_ptr<typename SRanges::IndexType>&... inds);
+    };
+    
+    template <typename T, class... SRanges>
+    class MultiArray : public MutableMultiArrayBase<T,SRanges...>
+    {
+    public:
+
+	typedef ContainerRange<SRanges...> CRange;
+	typedef MultiArrayBase<T,SRanges...> MAB;
+	typedef typename MultiArrayBase<T,SRanges...>::const_iterator const_iterator;
+	typedef typename MutableMultiArrayBase<T,SRanges...>::iterator iterator;
+	typedef typename CRange::IndexType IndexType;
+	
+	DEFAULT_MEMBERS(MultiArray);
+	MultiArray(const std::shared_ptr<SRanges>&... ranges);
+	MultiArray(const std::shared_ptr<SRanges>&... ranges, const std::vector<T>& vec);
+	MultiArray(const std::shared_ptr<SRanges>&... ranges, std::vector<T>&& vec);
+
+	// Only if ALL ranges have default extensions:
+	//MultiArray(const std::vector<T>& vec);
+	//MultiArray(std::vector<T>&& vec);
+	
+	// template <class Range2, class Range3>
+	// MultiArray(const MultiArray<MultiArray<T,Range2>,Range3> in);
+
+	// implement contstructor using FunctionalMultiArray as Input !!!
+	
+	//template <class Range2, class Range3>
+	//MultiArray& operator=(const MultiArray<MultiArray<T,Range2>,Range3> in);
+	
+	virtual T& operator[](const IndexType& i) override;
+	virtual const T& operator[](const IndexType& i) const override;
+	virtual T& at(const typename CRange::IndexType::MetaType& meta) override;
+	virtual const T& at(const typename CRange::IndexType::MetaType& meta) const override;
+	
+	virtual bool isConst() const override;	
+	virtual bool isSlice() const override;
+
+	template <class... SRanges2>
+	MultiArray<T,SRanges2...> format(const std::shared_ptr<SRanges2>&... nrs); // reformat array using 'nr' which in
+	//                                                                 total must have the same size as mRange
+	
+	virtual const T* data() const override;
+	virtual T* data() override;
+
+	virtual const std::vector<T>& datav() const override;
+	virtual std::vector<T>& datav() override;
+	
+	//	virtual void manipulate(ManipulatorBase<T>& mb,
+	//			const typename CRange::IndexType& manBegin,
+	//				const typename CRange::IndexType& manEnd);
+	
+	template <typename U, class... SRanges2>
+	friend class MultiArray;
+	
+    private:
+	std::vector<T> mCont;
+    };
+
+    template <typename T, class Function, class... SRanges>
+    class FunctionalMultiArray : public MultiArrayBase<T,SRanges...>
+    {
+    public:
+
+	typedef ContainerRange<SRanges...> CRange;
+	typedef MultiArrayBase<T,CRange> MAB;
+	typedef typename MultiArrayBase<T,CRange>::const_iterator const_iterator;
+	typedef typename CRange::IndexType IndexType;
+	
+	DEFAULT_MEMBERS(FunctionalMultiArray);
+	//FunctionalMultiArray(const CRange& range);
+	FunctionalMultiArray(const std::shared_ptr<SRanges>&... ranges, const Function& func);
+
+	virtual const T& operator[](const IndexType& i) const override;
+
+	virtual bool isConst() const override;	
+	virtual bool isSlice() const override;
+
+    protected:
+	mutable T mVal;
+	Function mFunc;
+    };
+    
+}
+
+/* ========================= *
+ * ---   TEMPLATE CODE   --- *
+ * ========================= */
+
+namespace MultiArrayTools
+{
+    /**************************************
+     *  MultiArrayBase::const_iterator    *	     
+     **************************************/
+
+    template <typename T, class... SRanges>
+    MultiArrayBase<T,SRanges...>::const_iterator::const_iterator(const MultiArrayBase<T,SRanges...>& ma):
+	mMAPtr(&ma), mPos(0) { }
+
+    template <typename T, class... SRanges>
+    MultiArrayBase<T,SRanges...>::const_iterator::const_iterator(const MultiArrayBase<T,SRanges...>& ma,
+							     const typename CRange::IndexType& index):
+	mMAPtr(&ma), mPos(index.pos()) { }
+    
+    template <typename T, class... SRanges>
+    bool MultiArrayBase<T,SRanges...>::const_iterator::operator==(const const_iterator& it) const
+    {
+	return mMAPtr == it.mMAPtr and mPos == it.mPos;
+    }
+
+    template <typename T, class... SRanges>
+    bool MultiArrayBase<T,SRanges...>::const_iterator::operator!=(const const_iterator& it) const
+    {
+	return mMAPtr != it.mMAPtr or mPos != it.mPos;
+    }
+    
+    template <typename T, class... SRanges>
+    const T& MultiArrayBase<T,SRanges...>::const_iterator::operator*() const
+    {
+	return mMAPtr->data()[mPos];
+    }
+
+    template <typename T, class... SRanges>
+    T const* MultiArrayBase<T,SRanges...>::const_iterator::operator->() const
+    {
+	return &mMAPtr->data()[mPos];
+    }
+    
+    template <typename T, class... SRanges>
+    typename MultiArrayBase<T,SRanges...>::const_iterator& MultiArrayBase<T,SRanges...>::const_iterator::operator++()
+    {
+	++mPos;
+	return *this;
+    }
+
+    template <typename T, class... SRanges>
+    typename MultiArrayBase<T,SRanges...>::const_iterator MultiArrayBase<T,SRanges...>::const_iterator::operator++(int)
+    {
+	const_iterator tmp(*this);
+	++mPos;
+	return tmp;
+    }
+
+    template <typename T, class... SRanges>
+    typename MultiArrayBase<T,SRanges...>::const_iterator& MultiArrayBase<T,SRanges...>::const_iterator::operator--()
+    {
+	--mPos;
+	return *this;
+    }
+
+    template <typename T, class... SRanges>
+    typename MultiArrayBase<T,SRanges...>::const_iterator MultiArrayBase<T,SRanges...>::const_iterator::operator--(int)
+    {
+	const_iterator tmp(*this);
+	--mPos;
+	return tmp;
+    }
+    
+    template <typename T, class... SRanges>
+    typename MultiArrayBase<T,SRanges...>::const_iterator& MultiArrayBase<T,SRanges...>::const_iterator::operator+=(int diff)
+    {
+	mPos += diff;
+	return *this;
+    }
+
+    template <typename T, class... SRanges>
+    typename MultiArrayBase<T,SRanges...>::const_iterator& MultiArrayBase<T,SRanges...>::const_iterator::operator-=(int diff)
+    {
+	mPos -= diff;
+	return *this;
+    }
+
+    template <typename T, class... SRanges>
+    typename MultiArrayBase<T,SRanges...>::const_iterator MultiArrayBase<T,SRanges...>::const_iterator::operator+(int num) const
+    {
+	const_iterator tmp(*this);
+	tmp += num;
+	return tmp;
+    }
+
+    template <typename T, class... SRanges>
+    typename MultiArrayBase<T,SRanges...>::const_iterator MultiArrayBase<T,SRanges...>::const_iterator::operator-(int num) const
+    {
+	const_iterator tmp(*this);
+	tmp -= num;
+    }
+    
+    template <typename T, class... SRanges>
+    int MultiArrayBase<T,SRanges...>::const_iterator::operator-(const const_iterator& it) const
+    {
+	return mPos - it.mPos;
+    }
+    
+    template <typename T, class... SRanges>
+    const T& MultiArrayBase<T,SRanges...>::const_iterator::operator[](int num) const
+    {
+	return *(operator+(num));
+    }
+    
+    template <typename T, class... SRanges>
+    bool MultiArrayBase<T,SRanges...>::const_iterator::operator<(const const_iterator& it) const
+    {
+	return mPos < it.mPos;
+    }
+
+    template <typename T, class... SRanges>
+    bool MultiArrayBase<T,SRanges...>::const_iterator::operator>(const const_iterator& it) const
+    {
+	return mPos > it.mPos;
+    }
+
+    template <typename T, class... SRanges>
+    bool MultiArrayBase<T,SRanges...>::const_iterator::operator<=(const const_iterator& it) const
+    {
+	return mPos <= it.mPos;
+    }
+
+    template <typename T, class... SRanges>
+    bool MultiArrayBase<T,SRanges...>::const_iterator::operator>=(const const_iterator& it) const
+    {
+	return mPos >= it.mPos;
+    }
+
+    template <typename T, class... SRanges>
+    typename MultiArrayBase<T,SRanges...>::IndexType
+    MultiArrayBase<T,SRanges...>::const_iterator::index() const
+    {
+	auto i = mMAPtr->beginIndex();
+	i = mPos;
+	return i;
+    }
+    
+    /**********************
+     *  MultiArrayBase    *	     
+     **********************/
+
+    template <typename T, class... SRanges>
+    MultiArrayBase<T,SRanges...>::MultiArrayBase(const std::shared_ptr<SRanges>&... ranges)
+    {
+	ContainerRangeFactory<SRanges...> crf(ranges...);
+	mRange = std::dynamic_pointer_cast<ContainerRange<SRanges...> >( crf.create() );
+    }
+    
+    template <typename T, class... SRanges>
+    size_t MultiArrayBase<T,SRanges...>::size() const
+    {
+	return mRange->size();
+    }
+
+    template <typename T, class... SRanges>
+    typename MultiArrayBase<T,SRanges...>::const_iterator MultiArrayBase<T,SRanges...>::begin() const
+    {
+	return const_iterator(*this, beginIndex());
+    }
+    
+    template <typename T, class... SRanges>
+    typename MultiArrayBase<T,SRanges...>::const_iterator MultiArrayBase<T,SRanges...>::end() const
+    {
+	return const_iterator(*this, endIndex());
+    }
+    
+    template <typename T, class... SRanges>
+    typename MultiArrayBase<T,SRanges...>::IndexType
+    MultiArrayBase<T,SRanges...>::beginIndex() const
+    {
+	return mRange->begin();
+    }
+
+    template <typename T, class... SRanges>
+    typename MultiArrayBase<T,SRanges...>::IndexType
+    MultiArrayBase<T,SRanges...>::endIndex() const
+    {
+	return mRange->end();
+    }
+
+    template <typename T, class... SRanges>
+    const std::shared_ptr<typename MultiArrayBase<T,SRanges...>::CRange>&
+    MultiArrayBase<T,SRanges...>::range() const
+    {
+	return mRange;
+    }
+
+    template <typename T, class... SRanges>
+    bool MultiArrayBase<T,SRanges...>::isConst() const
+    {
+	return true;
+    }
+    
+    template <typename T, class... SRanges>
+    ConstOperationRoot<T,SRanges...>
+    MultiArrayBase<T,SRanges...>::operator()(std::shared_ptr<typename SRanges::IndexType>&... inds) const
+    {
+	return ConstOperationRoot<T,SRanges...>(*this, inds...);
+    }
+    
+    template <typename T, class... SRanges>
+    bool MultiArrayBase<T,SRanges...>::isInit() const
+    {
+	return mInit;
+    }
+
+    /****************************************
+     *  MutableMultiArrayBase::iterator     *	     
+     ****************************************/
+    
+    template <typename T, class... SRanges>
+    MutableMultiArrayBase<T,SRanges...>::iterator::iterator(MutableMultiArrayBase<T,SRanges...>& ma):
+	mMAPtr(&ma), mPos(0)
+    { }
+    
+    template <typename T, class... SRanges>
+    MutableMultiArrayBase<T,SRanges...>::iterator::iterator(MutableMultiArrayBase<T,SRanges...>& ma,
+							const typename CRange::IndexType& index):
+	mMAPtr(&ma), mPos(index.pos())
+    { }
+    
+    template <typename T, class... SRanges>
+    bool MutableMultiArrayBase<T,SRanges...>::iterator::operator==(const iterator& it) const
+    {
+	return mMAPtr == it.mMAPtr and mPos == it.mPos;
+    }
+
+    template <typename T, class... SRanges>
+    bool MutableMultiArrayBase<T,SRanges...>::iterator::operator!=(const iterator& it) const
+    {
+	return mMAPtr != it.mMAPtr or mPos != it.mPos;
+    }
+    
+    template <typename T, class... SRanges>
+    const T& MutableMultiArrayBase<T,SRanges...>::iterator::operator*() const
+    {
+	return mMAPtr->data()[mPos];
+    }
+
+    template <typename T, class... SRanges>
+    T const* MutableMultiArrayBase<T,SRanges...>::iterator::operator->() const
+    {
+	return &mMAPtr->data()[mPos];
+    }
+
+    template <typename T, class... SRanges>
+    T& MutableMultiArrayBase<T,SRanges...>::iterator::operator*()
+    {
+	return mMAPtr->data()[mPos];
+    }
+
+    template <typename T, class... SRanges>
+    T* MutableMultiArrayBase<T,SRanges...>::iterator::operator->()
+    {
+	return &mMAPtr->data()[mPos];
+    }
+
+    template <typename T, class... SRanges>
+    typename MutableMultiArrayBase<T,SRanges...>::iterator& MutableMultiArrayBase<T,SRanges...>::iterator::operator++()
+    {
+	++mPos;
+	return *this;
+    }
+
+    template <typename T, class... SRanges>
+    typename MutableMultiArrayBase<T,SRanges...>::iterator MutableMultiArrayBase<T,SRanges...>::iterator::operator++(int)
+    {
+	iterator tmp(*this);
+	++mPos;
+	return tmp;
+    }
+
+    template <typename T, class... SRanges>
+    typename MutableMultiArrayBase<T,SRanges...>::iterator& MutableMultiArrayBase<T,SRanges...>::iterator::operator--()
+    {
+	--mPos;
+	return *this;
+    }
+
+    template <typename T, class... SRanges>
+    typename MutableMultiArrayBase<T,SRanges...>::iterator MutableMultiArrayBase<T,SRanges...>::iterator::operator--(int)
+    {
+	iterator tmp(*this);
+	--mPos;
+	return tmp;
+    }
+    
+    template <typename T, class... SRanges>
+    typename MutableMultiArrayBase<T,SRanges...>::iterator& MutableMultiArrayBase<T,SRanges...>::iterator::operator+=(int diff)
+    {
+	mPos += diff;
+	return *this;
+    }
+
+    template <typename T, class... SRanges>
+    typename MutableMultiArrayBase<T,SRanges...>::iterator& MutableMultiArrayBase<T,SRanges...>::iterator::operator-=(int diff)
+    {
+	mPos -= diff;
+	return *this;
+    }
+
+    template <typename T, class... SRanges>
+    typename MutableMultiArrayBase<T,SRanges...>::iterator MutableMultiArrayBase<T,SRanges...>::iterator::operator+(int num) const
+    {
+	iterator tmp(*this);
+	tmp += num;
+	return tmp;
+    }
+
+    template <typename T, class... SRanges>
+    typename MutableMultiArrayBase<T,SRanges...>::iterator MutableMultiArrayBase<T,SRanges...>::iterator::operator-(int num) const
+    {
+	iterator tmp(*this);
+	tmp -= num;
+    }
+    
+    template <typename T, class... SRanges>
+    int MutableMultiArrayBase<T,SRanges...>::iterator::operator-(const iterator& it) const
+    {
+	return mPos - it.mPos;
+    }
+    
+    template <typename T, class... SRanges>
+    const T& MutableMultiArrayBase<T,SRanges...>::iterator::operator[](int num) const
+    {
+	return *(operator+(num));
+    }
+
+    template <typename T, class... SRanges>
+    T& MutableMultiArrayBase<T,SRanges...>::iterator::operator[](int num)
+    {
+	return *(operator+(num));
+    }
+    
+    template <typename T, class... SRanges>
+    bool MutableMultiArrayBase<T,SRanges...>::iterator::operator<(const iterator& it) const
+    {
+	return mPos < it.mPos;
+    }
+
+    template <typename T, class... SRanges>
+    bool MutableMultiArrayBase<T,SRanges...>::iterator::operator>(const iterator& it) const
+    {
+	return mPos > it.mPos;
+    }
+
+    template <typename T, class... SRanges>
+    bool MutableMultiArrayBase<T,SRanges...>::iterator::operator<=(const iterator& it) const
+    {
+	return mPos <= it.mPos;
+    }
+
+    template <typename T, class... SRanges>
+    bool MutableMultiArrayBase<T,SRanges...>::iterator::operator>=(const iterator& it) const
+    {
+	return mPos >= it.mPos;
+    }
+
+    template <typename T, class... SRanges>
+    typename MutableMultiArrayBase<T,SRanges...>::IndexType
+    MutableMultiArrayBase<T,SRanges...>::iterator::index() const
+    {
+	auto i = mMAPtr->beginIndex();
+	i = mPos;
+	return i;
+    }    
+    
+    /******************************
+     *  MutableMultiArrayBase     *	     
+     ******************************/
+
+    template <typename T, class... SRanges>
+    MutableMultiArrayBase<T,SRanges...>::MutableMultiArrayBase(const std::shared_ptr<SRanges>&... ranges) :
+	MultiArrayBase<T,SRanges...>(ranges...) {}
+
+    template <typename T, class... SRanges>
+    typename MutableMultiArrayBase<T,SRanges...>::iterator MutableMultiArrayBase<T,SRanges...>::begin()
+    {
+	return iterator(*this, MAB::beginIndex());
+    }
+
+    template <typename T, class... SRanges>
+    typename MutableMultiArrayBase<T,SRanges...>::iterator MutableMultiArrayBase<T,SRanges...>::end()
+    {
+	return iterator(*this, MAB::endIndex());
+    }
+    
+    template <typename T, class... SRanges>
+    bool MutableMultiArrayBase<T,SRanges...>::isConst() const
+    {
+	return false;
+    }
+
+    template <typename T, class... SRanges>
+    OperationRoot<T,SRanges...>
+    MutableMultiArrayBase<T,SRanges...>::operator()(std::shared_ptr<typename SRanges::IndexType>&... inds)
+    {
+	return OperationRoot<T,SRanges...>(*this, inds...);
+    }
+    
+    template <typename T, class... SRanges>
+    ConstOperationRoot<T,SRanges...>
+    MutableMultiArrayBase<T,SRanges...>::operator()(std::shared_ptr<typename SRanges::IndexType>&... inds) const
+    {
+	return ConstOperationRoot<T,SRanges...>(*this, inds...);
+    }
+    
+    
+    /*******************
+     *  MultiArray     *	     
+     *******************/
+
+    template <typename T, class... SRanges>
+    MultiArray<T,SRanges...>::MultiArray(const std::shared_ptr<SRanges>&... ranges) :
+	MutableMultiArrayBase<T,SRanges...>(ranges...),
+	mCont(MAB::mRange->size())
+    {
+	MAB::mInit = true;
+    }
+
+    template <typename T, class... SRanges>
+    MultiArray<T,SRanges...>::MultiArray(const std::shared_ptr<SRanges>&... ranges, const std::vector<T>& vec) :
+	MutableMultiArrayBase<T,SRanges...>(ranges...),
+	mCont(vec)
+    {
+	MAB::mInit = true;
+	if(mCont.size() > MAB::mRange->size()){
+	    mCont.erase(mCont.begin() + MAB::mRange->size(), mCont.end());
+	}
+    }
+    
+    template <typename T, class... SRanges>
+    MultiArray<T,SRanges...>::MultiArray(const std::shared_ptr<SRanges>&... ranges, std::vector<T>&& vec) :
+	MutableMultiArrayBase<T,SRanges...>(ranges...),
+	mCont(vec)
+    {
+	MAB::mInit = true;
+	if(mCont.size() > MAB::mRange->size()){
+	    mCont.erase(mCont.begin() + MAB::mRange->size(), mCont.end());
+	}
+    }
+
+    /*
+    template <typename T, class... SRanges>
+    template <class Range2, class Range3>
+    MultiArray<T,SRanges...>::MultiArray(const MultiArray<MultiArray<T,Range2>,Range3> in) :
+	MutableMultiArrayBase<T,SRanges...>(merge(in.range(), in[ in.beginIndex() ].range()))
+	// assert that Range2 has always same extension
+    {
+	MAB::mInit = true;
+	mCont.clear();
+	for(auto i = in.beginIndex(); i != in.endIndex(); ++i){
+	    mCont.insert(mCont.end(), in[i].mCont.begin(), in[i].mCont.end());
+	}
+	assert(mCont.size() == MAB::mRange->size());
+    }
+    */
+    /*
+    template <typename T, class... SRanges>
+    template <class Range2, class Range3>
+    MultiArray<T,SRanges...>& MultiArray<T,SRanges...>::operator=(const MultiArray<MultiArray<T,Range2>,Range3> in)
+    {
+	MAB::mRange.reset(new Range(merge(in.range(), in[ in.beginIndex() ].range())));
+	// assert that Range2 has always same extension
+	mCont.clear();
+	for(auto i = in.beginIndex(); i != in.endIndex(); ++i){
+	    mCont.insert(mCont.end(), in[i].mCont.begin(), in[i].mCont.end());
+	}
+	assert(mCont.size() == MAB::mRange->size());
+	return *this;
+	} */   
+
+    template <typename T, class... SRanges>
+    T& MultiArray<T,SRanges...>::operator[](const typename CRange::IndexType& i)
+    {
+	return mCont[ i.pos() ];
+    }
+
+    template <typename T, class... SRanges>
+    const T& MultiArray<T,SRanges...>::operator[](const typename CRange::IndexType& i) const
+    {
+	return mCont[ i.pos() ];
+    }
+
+    template <typename T, class... SRanges>
+    T& MultiArray<T,SRanges...>::at(const typename CRange::IndexType::MetaType& meta)
+    {
+	return mCont[ MAB::beginIndex().at(meta).pos() ];
+    }
+
+    template <typename T, class... SRanges>
+    const T& MultiArray<T,SRanges...>::at(const typename CRange::IndexType::MetaType& meta) const
+    {
+	return mCont[ MAB::beginIndex().at(meta).pos() ];
+    }
+	
+    template <typename T, class... SRanges>
+    bool MultiArray<T,SRanges...>::isConst() const
+    {
+	return false;
+    }
+    
+    template <typename T, class... SRanges>
+    bool MultiArray<T,SRanges...>::isSlice() const
+    {
+	return false;
+    }
+
+    template <typename T, class... SRanges>
+    template <class... SRanges2>
+    MultiArray<T,SRanges2...> MultiArray<T,SRanges...>::format(const std::shared_ptr<SRanges2>&... nrs)
+    {
+	return MultiArray<T,SRanges2...>( nrs... , std::move(mCont) );
+    }
+    
+    template <typename T, class... SRanges>
+    const T* MultiArray<T,SRanges...>::data() const
+    {
+	return mCont.data();
+    }
+
+    template <typename T, class... SRanges>
+    T* MultiArray<T,SRanges...>::data()
+    {
+	return mCont.data();
+    }
+
+    template <typename T, class... SRanges>
+    const std::vector<T>& MultiArray<T,SRanges...>::datav() const
+    {
+	return mCont;
+    }
+
+    template <typename T, class... SRanges>
+    std::vector<T>& MultiArray<T,SRanges...>::datav()
+    {
+	return mCont;
+    }
+
+    
+    /*
+    template <typename T, class... SRanges>
+    void MultiArray<T,SRanges...>::manipulate(ManipulatorBase<T>& mb,
+					 const typename Range::IndexType& manBegin,
+					 const typename Range::IndexType& manEnd)
+    {
+	mb.setup(mCont, manBegin.pos(), manEnd.pos());
+	mb.execute();
+    }
+    */
+
+
+    /****************************
+     *   FunctionalMultiArray   *	     
+     ****************************/
+
+    /*
+    template <typename T, class Range, class Function>
+    FunctionalMultiArray<T,Range,Function>::FunctionalMultiArray(const Range& range) :
+	MultiArrayBase<T,SRanges...>(range), mFunc() {}
+    */
+    template <typename T, class Function, class... SRanges>
+    FunctionalMultiArray<T,Function,SRanges...>::FunctionalMultiArray(const std::shared_ptr<SRanges>&... ranges,
+								  const Function& func) :
+	MultiArrayBase<T,SRanges...>(ranges...), mFunc(func) {}
+
+    template <typename T, class Function, class... SRanges>
+    const T& FunctionalMultiArray<T,Function,SRanges...>::operator[](const typename CRange::IndexType& i) const
+    {
+	mVal = mFunc(i);
+	return mVal;
+    }
+
+    template <typename T, class Function, class... SRanges>
+    bool FunctionalMultiArray<T,Function,SRanges...>::isConst() const
+    {
+	return true;
+    }
+
+    template <typename T, class Function, class... SRanges>
+    bool FunctionalMultiArray<T,Function,SRanges...>::isSlice() const
+    {
+	return false;
+    }
+    
+}
+
+#endif

src/include/multi_array_header.h

@@ -0,0 +1,60 @@
+// -*- C++ -*-
+
+#ifndef __multi_array_header_h__
+#define __multi_array_header_h__
+
+#include <cstdlib>
+//#include "base_def.h"
+//#include "range_base.h"
+//#include "index_base.h"
+//#include "single_range.h"
+//#include "multi_range.h"
+//#include "container_range.h"
+//#include "block.h"
+#include "multi_array_operation.h"
+#include "multi_array.h"
+#include "helper_tools.h"
+//#include "slice.h"
+//#include "manipulator.h"
+//#include "range_transformer.h"
+//#include "ma_functional.h"
+
+namespace MultiArrayTools
+{
+    
+    /*********************************
+     *   Some standard definitions   *
+     *********************************/
+    /*
+    // ===== Index Types =====
+
+    typedef SingleIndex<size_t,RangeType::ANY> GenericNI;
+    typedef SingleIndex<double,RangeType::ANY> GenericFI;
+    typedef SingleIndex<size_t,RangeType::LORENTZ> LorentzI;
+    typedef SingleIndex<int,RangeType::SPACE> Space1dNI;
+    typedef SingleIndex<int,RangeType::MOMENTUM> Mom1dNI;
+    typedef SingleIndex<size_t,RangeType::ENSEMBLE> EnsI;
+    typedef SingleIndex<std::array<int,3>, RangeType::MOMENTUM> LinMomI;
+    typedef SingleIndex<std::array<int,4>, RangeType::MOMENTUM> LinMom4dI;
+    // ...
+
+    // ===== Range Types =====
+
+    typedef SingleRange<size_t,RangeType::ANY> GenericNR;
+    typedef SingleRange<double,RangeType::ANY> GenericFR;
+    typedef SingleRange<size_t,RangeType::LORENTZ> LorentzR;
+    typedef SingleRange<int,RangeType::SPACE> Space1dNR;
+    typedef SingleRange<size_t,RangeType::DISTANCE> DistanceNR;
+    typedef SingleRange<int,RangeType::MOMENTUM> Mom1dNR;
+    typedef SingleRange<size_t, RangeType::ENSEMBLE> EnsR;
+    typedef MultiRange<Space1dNR,Space1dNR,Space1dNR> Space3dNR;
+    typedef MultiRange<Mom1dNR,Mom1dNR,Mom1dNR> Mom3dNR;
+    typedef SingleRange<VET, RangeType::VALUE_ERROR> ValErrR;
+    typedef SingleRange<std::array<int,3>, RangeType::MOMENTUM> LinMomR;
+    typedef SingleRange<std::array<int,4>, RangeType::MOMENTUM> LinMom4dR;
+    // ...
+    */
+
+}
+
+#endif

src/include/multi_array_operation.h

@@ -0,0 +1,599 @@
+// -*- C++ -*-
+
+#ifndef __multi_array_operation_h__
+#define __multi_array_operation_h__
+
+#include <cstdlib>
+#include <tuple>
+#include <cmath>
+#include <map>
+#include <utility>
+
+#include "base_def.h"
+#include "mbase_def.h"
+
+//#include "block/block.h"
+//#include "operation_utils.h"
+#include "ranges/rheader.h"
+#include "pack_num.h"
+
+#include "ranges/index_info.h"
+#include "arith.h"
+
+namespace MultiArrayTools
+{
+
+    namespace
+    {
+	using namespace MultiArrayHelper;
+    }
+    
+    template <typename T, class OperationClass>
+    class OperationTemplate
+    {
+    public:
+	
+	OperationClass& THIS() { return static_cast<OperationClass&>(*this); }
+	const OperationClass& THIS() const { return static_cast<OperationClass const&>(*this); }
+
+	template <class Second>
+	auto operator+(const Second& in) const
+	    -> Operation<T,plus<T>,OperationClass,Second>;
+	
+	template <class Second>
+	auto operator-(const Second& in) const
+	    -> Operation<T,minus<T>,OperationClass,Second>;
+	
+	template <class Second>
+	auto operator*(const Second& in) const
+	    -> Operation<T,multiplies<T>,OperationClass,Second>;
+	
+	template <class Second>
+	auto operator/(const Second& in) const
+	    -> Operation<T,divides<T>,OperationClass,Second>;
+
+	template <class IndexType>
+	auto c(std::shared_ptr<IndexType>& ind) const
+	    -> Contraction<T,OperationClass,IndexType>;
+
+    private:
+	friend OperationClass;
+	OperationTemplate() = default;
+    };
+    
+    template <typename T, class OpClass, class... Ranges>
+    class OperationMaster
+    {
+    public:
+
+	class AssignmentExpr
+	{
+	private:
+	    AssignmentExpr() = default;
+	    	    
+	    OperationMaster& mM;
+	    const OpClass& mSec;
+
+	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(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;
+
+	    auto rootSteps(std::intptr_t iPtrNum = 0) const -> ExtType;
+	    
+	};
+	
+	typedef T value_type;
+	typedef OperationBase<T> OB;
+	typedef ContainerRange<Ranges...> CRange;
+	typedef typename MultiRange<Ranges...>::IndexType IndexType;
+	
+	OperationMaster(MutableMultiArrayBase<T,Ranges...>& ma, const OpClass& second,
+			std::shared_ptr<typename CRange::IndexType>& index);
+
+	OperationMaster(MutableMultiArrayBase<T,Ranges...>& ma, const OpClass& second,
+			std::shared_ptr<typename CRange::IndexType>& index,
+			const IndexInfo* blockIndex);
+
+	inline void set(size_t pos, T val) { mDataPtr[pos] = val; }
+	inline void add(size_t pos, T val) { mDataPtr[pos] += val; }
+	inline T get(size_t pos) const;
+
+    private:
+
+	std::shared_ptr<IndexType> mkIndex(std::shared_ptr<typename CRange::IndexType>& index);
+	void performAssignment(std::intptr_t blockIndexNum);
+	OpClass const& mSecond;
+	MutableMultiArrayBase<T,Ranges...>& mArrayRef;
+	T* mDataPtr;
+	std::shared_ptr<IndexType> mIndex;
+	IndexInfo mIInfo;
+    };
+
+    
+    template <typename T, class... Ranges>
+    class ConstOperationRoot : /*public OperationBase<T>,*/
+			       public OperationTemplate<T,ConstOperationRoot<T,Ranges...> >
+    {
+    public:
+
+	typedef T value_type;
+	typedef OperationBase<T> OB;
+	typedef OperationTemplate<T,ConstOperationRoot<T,Ranges...> > OT;
+	typedef ContainerRange<Ranges...> CRange;
+	typedef typename CRange::IndexType IndexType;
+
+	static constexpr size_t SIZE = 1;
+	
+	ConstOperationRoot(const MultiArrayBase<T,Ranges...>& ma,
+			   const std::shared_ptr<typename Ranges::IndexType>&... indices);
+	
+	template <class ET>
+	inline T get(ET pos) const;
+
+	MExt<void> rootSteps(std::intptr_t iPtrNum = 0) const; // nullptr for simple usage with decltype
+
+	template <class Expr>
+	Expr loop(Expr exp) const;
+		
+    private:
+
+	std::shared_ptr<IndexType>
+	mkIndex(const MultiArrayBase<T,Ranges...>& ma,
+		const std::shared_ptr<typename Ranges::IndexType>&... indices);
+
+	MultiArrayBase<T,Ranges...> const& mArrayRef;
+	const T* mDataPtr;
+	std::shared_ptr<IndexType> mIndex;
+	IndexInfo mIInfo;
+    };
+    
+    template <typename T, class... Ranges>
+    class OperationRoot : public OperationTemplate<T,OperationRoot<T,Ranges...> >
+    {
+    public:
+
+	typedef T value_type;
+	typedef OperationBase<T> OB;
+	typedef OperationTemplate<T,OperationRoot<T,Ranges...> > OT;
+	typedef ContainerRange<Ranges...> CRange;
+	typedef typename CRange::IndexType IndexType;
+
+	static constexpr size_t SIZE = 1;
+	
+	OperationRoot(MutableMultiArrayBase<T,Ranges...>& ma,
+		      const std::shared_ptr<typename Ranges::IndexType>&... indices);
+
+	template <class OpClass>
+	OperationMaster<T,OpClass,Ranges...> operator=(const OpClass& in);
+
+	template <class ET>
+	inline T get(ET pos) const;
+
+	MExt<void> rootSteps(std::intptr_t iPtrNum = 0) const; // nullptr for simple usage with decltype
+
+	template <class Expr>
+	Expr loop(Expr exp) const;
+	
+    private:
+
+	std::shared_ptr<IndexType>
+	mkIndex(const MultiArrayBase<T,Ranges...>& ma,
+		const std::shared_ptr<typename Ranges::IndexType>&... indices);
+
+	MutableMultiArrayBase<T,Ranges...>& mArrayRef;
+	T* mDataPtr;
+	std::shared_ptr<IndexType> mIndex;
+	IndexInfo mIInfo;
+    };
+
+    template <class Op>
+    size_t sumRootNum()
+    {
+	return typename Op::rootNum();
+    }
+    
+    template <class Op1, class Op2, class... Ops>
+    size_t sumRootNum()
+    {
+	return typename Op1::rootNum() + sumRootNum<Op2,Ops...>();
+    }
+
+    template <size_t N>
+    struct RootSumN
+    {
+	template <class Op1, class... Ops>
+	struct rs
+	{
+	    static constexpr size_t SIZE = Op1::SIZE + RootSumN<N-1>::template rs<Ops...>::SIZE;
+	};
+    };
+
+    template <>
+    struct RootSumN<0>
+    {
+	template <class Op1>
+	struct rs
+	{
+	    static constexpr size_t SIZE = Op1::SIZE;
+	};
+    };
+
+    
+    template <class... Ops>
+    struct RootSum
+    {
+	static constexpr size_t SIZE = RootSumN<sizeof...(Ops)-1>::template rs<Ops...>::SIZE;
+    };
+
+    
+    template <typename T, class OpFunction, class... Ops>
+    class Operation : public OperationTemplate<T,Operation<T,OpFunction,Ops...> >
+    {
+    public:
+
+	typedef T value_type;
+	typedef OperationBase<T> OB;
+	typedef OperationTemplate<T,Operation<T,OpFunction,Ops...> > OT;
+	typedef OpFunction F;
+	
+	static constexpr size_t SIZE = RootSum<Ops...>::SIZE;
+
+    private:
+	std::tuple<Ops...> mOps;
+
+    public:
+	typedef decltype(PackNum<sizeof...(Ops)-1>::template mkSteps<Ops...>(0, mOps)) ETuple;
+		
+	Operation(const Ops&... ops);
+	
+	template <class ET>
+	inline T get(ET pos) const;
+
+	auto rootSteps(std::intptr_t iPtrNum = 0) const // nullptr for simple usage with decltype
+	    -> decltype(PackNum<sizeof...(Ops)-1>::mkSteps(iPtrNum, mOps));
+
+	template <class Expr>
+	auto loop(Expr exp) const
+	    -> decltype(PackNum<sizeof...(Ops)-1>::mkLoop( mOps, exp));
+	
+    };
+    
+    template <typename T, class Op, class IndexType>
+    class Contraction : public OperationTemplate<T,Contraction<T,Op,IndexType> >
+    {
+    public:
+
+	typedef T value_type;
+	typedef OperationTemplate<T,Contraction<T,Op,IndexType> > OT;
+
+	static constexpr size_t SIZE = Op::SIZE;
+	
+    private:
+
+	const Op& mOp;
+	std::shared_ptr<IndexType> mInd;
+
+    public:
+	typedef decltype(mOp.rootSteps(0)) ETuple;
+	
+	Contraction(const Op& op, std::shared_ptr<IndexType> ind);
+
+	template <class ET>
+	inline T get(ET pos) const;
+
+	auto rootSteps(std::intptr_t iPtrNum = 0) const  // nullptr for simple usage with decltype
+	    -> decltype(mOp.rootSteps(iPtrNum));
+
+	template <class Expr>
+	auto loop(Expr exp) const -> decltype(mInd->iforh(exp));
+    };
+    
+}
+
+/* ========================= *
+ * ---   TEMPLATE CODE   --- *
+ * ========================= */
+
+namespace MultiArrayTools
+{
+    namespace
+    {
+	using namespace MultiArrayHelper;
+    }
+    
+    /***************************
+     *   OperationTemplate     *
+     ***************************/
+    
+    template <typename T, class OperationClass>
+    template <class Second>
+    auto OperationTemplate<T,OperationClass>::operator+(const Second& in) const
+	-> Operation<T,plus<T>,OperationClass,Second>
+    {
+    	return Operation<T,plus<T>,OperationClass,Second>(THIS(), in);
+    }
+
+    template <typename T, class OperationClass>
+    template <class Second>
+    auto OperationTemplate<T,OperationClass>::operator-(const Second& in) const
+	-> Operation<T,minus<T>,OperationClass,Second>
+    {
+    	return Operation<T,minus<T>,OperationClass,Second>(THIS(), in);
+    }
+    
+    template <typename T, class OperationClass>
+    template <class Second>
+    auto OperationTemplate<T,OperationClass>::operator*(const Second& in) const
+	-> Operation<T,multiplies<T>,OperationClass,Second>
+    {
+    	return Operation<T,multiplies<T>,OperationClass,Second>(THIS(), in);
+    }
+
+    template <typename T, class OperationClass>
+    template <class Second>
+    auto OperationTemplate<T,OperationClass>::operator/(const Second& in) const
+	-> Operation<T,divides<T>,OperationClass,Second>
+    {
+    	return Operation<T,divides<T>,OperationClass,Second>(THIS(), in);
+    }
+
+    template <typename T, class OperationClass>
+    template <class IndexType>
+    auto OperationTemplate<T,OperationClass>::c(std::shared_ptr<IndexType>& ind) const
+	-> Contraction<T,OperationClass,IndexType>
+    {
+	return Contraction<T,OperationClass,IndexType>(THIS(), ind);
+    }
+
+    /*****************************************
+     *   OperationMaster::AssignmentExpr     *
+     *****************************************/
+
+    template <typename T, class OpClass, class... Ranges>
+    OperationMaster<T,OpClass,Ranges...>::AssignmentExpr::
+    AssignmentExpr(OperationMaster& m, const OpClass& sec) :
+	mM(m), mSec(sec) {}
+
+    template <typename T, class OpClass, class... Ranges>
+    inline void OperationMaster<T,OpClass,Ranges...>::AssignmentExpr::
+    operator()(size_t start, ExtType last) const
+    {
+	mM.add(start, mSec.template get<ExtType>(last) );
+    }
+    
+    template <typename T, class OpClass, class... Ranges>
+    typename OperationMaster<T,OpClass,Ranges...>::AssignmentExpr::ExtType
+    OperationMaster<T,OpClass,Ranges...>::AssignmentExpr::
+    rootSteps(std::intptr_t iPtrNum) const
+    {
+	return mSec.rootSteps(iPtrNum);
+    }
+    
+    
+    /*************************
+     *   OperationMaster     *
+     *************************/
+
+    template <typename T, class OpClass, class... Ranges>
+    OperationMaster<T,OpClass,Ranges...>::
+    OperationMaster(MutableMultiArrayBase<T,Ranges...>& ma, const OpClass& second,
+		    std::shared_ptr<typename CRange::IndexType>& index) :
+	mSecond(second), mArrayRef(ma), mDataPtr(mArrayRef.data()),
+	mIndex(mkIndex(index)), mIInfo(*mIndex)
+    {
+	performAssignment(0);
+    }
+
+    template <typename T, class OpClass, class... Ranges>
+    OperationMaster<T,OpClass,Ranges...>::
+    OperationMaster(MutableMultiArrayBase<T,Ranges...>& ma, const OpClass& second,
+		    std::shared_ptr<typename CRange::IndexType>& index,
+		    const IndexInfo* blockIndex) :
+	mSecond(second), mArrayRef(ma), mDataPtr(mArrayRef.data()),
+	mIndex(mkIndex(index)), mIInfo(*mIndex)
+    {
+	performAssignment(0);
+    }
+
+    template <typename T, class OpClass, class... Ranges>
+    std::shared_ptr<typename OperationMaster<T,OpClass,Ranges...>::IndexType>
+    OperationMaster<T,OpClass,Ranges...>::
+    mkIndex(std::shared_ptr<typename CRange::IndexType>& index)
+    {
+	MultiRangeFactory<Ranges...> mrf( index->range() );
+	std::shared_ptr<MultiRange<Ranges...> > mr =
+	    std::dynamic_pointer_cast<MultiRange<Ranges...> >( mrf.create() );
+	auto i = std::make_shared<IndexType>( mr->begin() );
+	(*i) = *index;
+	return i;
+    }
+    
+    template <typename T, class OpClass, class... Ranges>
+    void OperationMaster<T,OpClass,Ranges...>::performAssignment(std::intptr_t blockIndexNum)
+    {
+	AssignmentExpr ae(*this, mSecond); // Expression to be executed within loop
+	const auto loop = mSecond.template loop<decltype(mIndex->ifor(ae))>( mIndex->ifor(ae) );
+	// hidden Loops outside ! -> auto vectorizable
+	loop(); // execute overall loop(s) and so internal hidden loops and so the inherited expressions
+    }
+
+    template <typename T, class OpClass, class... Ranges>
+    inline T OperationMaster<T,OpClass,Ranges...>::get(size_t pos) const
+    {
+	return mDataPtr[pos];
+    }
+
+   
+    /****************************
+     *   ConstOperationRoot     *
+     ****************************/
+
+    template <typename T, class... Ranges>
+    ConstOperationRoot<T,Ranges...>::
+    ConstOperationRoot(const MultiArrayBase<T,Ranges...>& ma,
+		       const std::shared_ptr<typename Ranges::IndexType>&... indices) :
+	mArrayRef(ma), mDataPtr(mArrayRef.data()),
+	mIndex( mkIndex(ma,indices...) ), mIInfo(*mIndex)
+    {}
+
+    template <typename T, class... Ranges>
+    std::shared_ptr<typename ConstOperationRoot<T,Ranges...>::IndexType>
+    ConstOperationRoot<T,Ranges...>::
+    mkIndex(const MultiArrayBase<T,Ranges...>& ma,
+	    const std::shared_ptr<typename Ranges::IndexType>&... indices)
+    {
+	auto i = std::make_shared<IndexType>( ma.range() );
+	(*mIndex)(indices...);
+	return i;
+    }
+
+    template <typename T, class... Ranges>
+    template <class ET>
+    inline T ConstOperationRoot<T,Ranges...>::get(ET pos) const
+    {
+	return mDataPtr[pos.val()];
+    }
+
+    template <typename T, class... Ranges>
+    MExt<void> ConstOperationRoot<T,Ranges...>::rootSteps(std::intptr_t iPtrNum) const
+    {
+	return MExt<void>(getStepSize( mIndex->info(), iPtrNum ));
+	//return MExt<void>(getStepSize( getRootIndices( mIndex->info() ), iPtrNum ));
+    }
+
+
+    template <typename T, class... Ranges>
+    template <class Expr>
+    Expr ConstOperationRoot<T,Ranges...>::loop(Expr exp) const
+    {
+	return exp;
+    }
+    
+    /***********************
+     *   OperationRoot     *
+     ***********************/
+    
+    template <typename T, class... Ranges>
+    OperationRoot<T,Ranges...>::
+    OperationRoot(MutableMultiArrayBase<T,Ranges...>& ma,
+		  const std::shared_ptr<typename Ranges::IndexType>&... indices) :
+	mArrayRef(ma), mDataPtr(mArrayRef.data()),
+	mIndex( mkIndex( ma, indices... ) ), mIInfo(*mIndex)
+    {}
+
+    template <typename T, class... Ranges>
+    std::shared_ptr<typename OperationRoot<T,Ranges...>::IndexType>
+    OperationRoot<T,Ranges...>::
+    mkIndex(const MultiArrayBase<T,Ranges...>& ma,
+	    const std::shared_ptr<typename Ranges::IndexType>&... indices)
+    {
+	auto i = std::make_shared<IndexType>( ma.range() );
+	(*mIndex)(indices...);
+	return i;
+    }
+    
+    template <typename T, class... Ranges>
+    template <class OpClass>
+    OperationMaster<T,OpClass,Ranges...> OperationRoot<T,Ranges...>::operator=(const OpClass& in)
+    {
+	return OperationMaster<T,OpClass,Ranges...>(mArrayRef, in, mIndex);
+    }
+
+    template <typename T, class... Ranges>
+    template <class ET>
+    inline T OperationRoot<T,Ranges...>::get(ET pos) const
+    {
+	return mDataPtr[pos.val()];
+    }
+
+    template <typename T, class... Ranges>
+    MExt<void> OperationRoot<T,Ranges...>::rootSteps(std::intptr_t iPtrNum) const
+    {
+	return MExt<void>(getStepSize( mIndex->info(), iPtrNum ));
+	//return MExt<void>(getStepSize( getRootIndices( mIndex->info() ), iPtrNum ));
+    }
+
+    template <typename T, class... Ranges>
+    template <class Expr>
+    Expr OperationRoot<T,Ranges...>::loop(Expr exp) const
+    {
+	return exp;
+    }
+    
+    /*******************
+     *   Operation     *
+     *******************/
+    
+    template <typename T, class OpFunction, class... Ops>
+    Operation<T,OpFunction,Ops...>::Operation(const Ops&... ops) :
+	mOps(ops...) {}
+
+    template <typename T, class OpFunction, class... Ops>
+    template <class ET>
+    inline T Operation<T,OpFunction,Ops...>::get(ET pos) const
+    {
+	typedef std::tuple<Ops...> OpTuple;
+	return PackNum<sizeof...(Ops)-1>::
+	    template mkOpExpr<SIZE,T,ET,OpTuple,OpFunction>(pos, mOps);
+    }
+
+    template <typename T, class OpFunction, class... Ops>
+    auto Operation<T,OpFunction,Ops...>::rootSteps(std::intptr_t iPtrNum) const
+	-> decltype(PackNum<sizeof...(Ops)-1>::mkSteps(iPtrNum, mOps))
+    {
+	return PackNum<sizeof...(Ops)-1>::mkSteps(iPtrNum, mOps);
+    }
+
+    template <typename T, class OpFunction, class... Ops>
+    template <class Expr>
+    auto Operation<T,OpFunction,Ops...>::loop(Expr exp) const
+	-> decltype(PackNum<sizeof...(Ops)-1>::mkLoop( mOps, exp ))
+    {
+	return PackNum<sizeof...(Ops)-1>::mkLoop( mOps, exp );
+    }
+
+    
+    /*********************
+     *   Contraction     *
+     *********************/
+
+    template <typename T, class Op, class IndexType>
+    Contraction<T,Op,IndexType>::Contraction(const Op& op, std::shared_ptr<IndexType> ind) :
+	mOp(op),
+	mInd(ind) {}
+
+    // forward loop !!!!
+    template <typename T, class Op, class IndexType>
+    template <class ET>
+    inline T Contraction<T,Op,IndexType>::get(ET pos) const
+    {
+	return mOp.template get<ET>(pos);
+    }
+
+    template <typename T, class Op, class IndexType>
+    auto Contraction<T,Op,IndexType>::rootSteps(std::intptr_t iPtrNum) const
+	-> decltype(mOp.rootSteps(iPtrNum))
+    {
+	return mOp.rootSteps(iPtrNum);
+    }
+
+    template <typename T, class Op, class IndexType>
+    template <class Expr>
+    auto Contraction<T,Op,IndexType>::loop(Expr exp) const -> decltype(mInd->iforh(exp))
+    {
+	return mInd->iforh(exp);
+    }
+
+}
+
+#endif

src/include/pack_num.h

@@ -0,0 +1,108 @@
+// -*- C++ -*-
+
+#ifndef __pack_num_h__
+#define __pack_num_h__
+
+#include <cstdlib>
+//#include <type_traits>
+#include <tuple>
+#include <ostream>
+
+#include "base_def.h"
+#include "xfor/exttype.h"
+
+namespace MultiArrayHelper
+{
+	
+    template <size_t N>
+    struct PackNum
+    {
+	template <typename... T>
+	static void printTuple(std::ostream& out, const std::tuple<T...>& tp)
+	{
+	    out << std::get<sizeof...(T)-N-1>(tp) << ", ";
+	    PackNum<N-1>::printTuple(out, tp);
+	}
+
+	template <class... Ops>
+	static auto mkSteps(std::intptr_t ii, const std::tuple<Ops...>& otp)
+	    -> decltype(PackNum<N-1>::mkSteps(ii, otp).extend( std::get<N>(otp).rootSteps(ii)) )
+	{
+	    return PackNum<N-1>::mkSteps(ii, otp).extend( std::get<N>(otp).rootSteps(ii));
+	}
+	
+	template <class RootStepTuple, class IndexClass, class OpClass>
+	static void mkExt(std::array<RootStepTuple,IndexClass::totalDim()>& out,
+			  const std::array<std::intptr_t,IndexClass::totalDim()>& siar,
+			  const OpClass& second)
+	{
+	    std::get<N>(out) = second.rootSteps( std::get<N>(siar) );
+	    PackNum<N-1>::mkExt(out, siar, second);
+	}
+
+	template <size_t LAST, typename T, class ETuple, class OpTuple, class OpFunction, typename... Args>
+	static inline T mkOpExpr(const ETuple& pos, const OpTuple& ops, const Args&... args)
+	{
+	    typedef typename std::remove_reference<decltype(std::get<N>(ops))>::type NextOpType;
+	    static_assert(LAST > NextOpType::SIZE, "inconsistent array positions");
+	    static constexpr size_t NEXT = LAST - NextOpType::SIZE;
+	    return PackNum<N-1>::template mkOpExpr<NEXT,T,ETuple,OpTuple,OpFunction,T,Args...>
+		( pos, ops, std::get<N>(ops).get(Getter<NEXT>::template getX<ETuple>( pos )), args...);
+	}
+
+	template <class OpTuple, class Expr>
+	static auto mkLoop( const OpTuple& ot, Expr&& exp )
+	    -> decltype(std::get<N>(ot).loop( PackNum<N-1>::mkLoop(ot,exp) ))
+	{
+	    return std::get<N>(ot).loop( PackNum<N-1>::mkLoop(ot,exp) );
+	}
+    };
+    
+    template<>
+    struct PackNum<0>
+    {
+	template <typename... T>
+	static void printTuple(std::ostream& out, const std::tuple<T...>& tp)
+	{	    
+	    out << std::get<sizeof...(T)-1>(tp);
+	}
+
+	template <class... Ops>
+	static auto mkSteps(std::intptr_t ii, const std::tuple<Ops...>& otp)
+	    -> decltype(std::get<0>(otp).rootSteps(ii))
+	{
+	    return std::get<0>(otp).rootSteps(ii);
+	}
+
+	template <class RootStepTuple, class IndexClass, class OpClass>
+	static void mkExt(std::array<RootStepTuple,IndexClass::totalDim()>& out,
+			  const std::array<std::intptr_t,IndexClass::totalDim()>& siar,
+			  const OpClass& second)
+	{
+	    std::get<0>(out) = second.rootSteps( std::get<0>(siar) );
+	}
+
+	template <size_t LAST, typename T, class ETuple, class OpTuple, class OpFunction, typename... Args>
+	static inline T mkOpExpr(const ETuple& pos, const OpTuple& ops, const Args&... args)
+	{
+	    typedef typename std::remove_reference<decltype(std::get<0>(ops))>::type NextOpType;
+	    static constexpr size_t NEXT = LAST - NextOpType::SIZE;
+	    static_assert(NEXT == 0, "inconsistent array positions");
+	    return OpFunction::apply(std::get<0>(ops).get(Getter<0>::template getX<ETuple>( pos )), args...);
+	}
+
+	template <class OpTuple, class Expr>
+	static auto mkLoop( const OpTuple& ot, Expr&& exp )
+	    -> decltype(std::get<0>(ot).loop( exp ))
+	{
+	    return std::get<0>(ot).loop( exp );
+	}
+
+	
+    };
+
+
+    
+} // end namespace MultiArrayHelper
+
+#endif

src/include/ranges/anonymous_range.h

@@ -0,0 +1,168 @@
+// -*- C++ -*-
+
+#ifndef __anonymous_range_h__
+#define __anonymous_range_h__
+
+#include <cstdlib>
+//#include "base_def.h"
+#include "ranges/range_base.h"
+#include "ranges/single_range.h"
+
+namespace MultiArrayTools
+{
+
+    typedef SingleIndex<size_t,SpaceType::NONE> AnonymousIndex;
+
+    class AnonymousRangeFactory : public RangeFactoryBase
+    {
+    public:
+	
+	typedef AnonymousRange oType;
+
+	AnonymousRangeFactory() = delete;
+
+	template <class... RangeTypes>
+	AnonymousRangeFactory(const std::tuple<std::shared_ptr<RangeTypes>...>& origs);
+	
+	template <class... RangeTypes>
+	AnonymousRangeFactory(std::shared_ptr<RangeTypes>... origs);
+
+	std::shared_ptr<RangeBase> create();
+	
+    };
+    
+    class AnonymousRange : public RangeInterface<AnonymousIndex>
+    {
+	typedef RangeBase RB;
+	typedef typename RangeInterface<AnonymousIndex>::IndexType IndexType;
+	
+	virtual size_t size() const override;
+	virtual size_t dim() const override;
+	
+	size_t get(size_t pos) const;
+	size_t getMeta(size_t metaPos) const;
+	
+	virtual IndexType begin() const override;
+	virtual IndexType end() const override;
+	virtual std::shared_ptr<VIWB> index() const override;
+	
+	friend AnonymousRangeFactory;
+	
+    protected:
+
+	AnonymousRange() = delete;
+	AnonymousRange(const AnonymousRange& in) = delete;
+
+	template <class... RangeTypes>
+	AnonymousRange(const std::tuple<std::shared_ptr<RangeTypes>...>& origs);
+	
+	template <class... RangeTypes>
+	AnonymousRange(std::shared_ptr<RangeTypes>... origs);
+
+	size_t mSize;
+	
+	std::vector<std::shared_ptr<RangeBase> > mOrig;
+    };
+
+}
+
+/* ========================= *
+ * ---   TEMPLATE CODE   --- *
+ * ========================= */
+
+namespace MultiArrayTools
+{
+
+    /***********************
+     *   AnonymousRange    *
+     ***********************/
+
+    template <class... RangeTypes>
+    AnonymousRangeFactory::AnonymousRangeFactory(const std::tuple<std::shared_ptr<RangeTypes>...>& origs)
+    {
+	mProd = std::shared_ptr<oType>( new AnonymousRange( origs ) );
+    }
+
+    template <class... RangeTypes>
+    AnonymousRangeFactory::AnonymousRangeFactory(std::shared_ptr<RangeTypes>... origs)
+    {
+	mProd = std::shared_ptr<oType>( new AnonymousRange( origs... ) );
+    }
+        
+    std::shared_ptr<RangeBase> AnonymousRangeFactory::create()
+    {
+	setSelf();
+	return mProd;
+    }
+    
+    /***********************
+     *   AnonymousRange    *
+     ***********************/
+
+    template <class... RangeTypes>
+    AnonymousRange::AnonymousRange(const std::tuple<std::shared_ptr<RangeTypes>...>& origs) :
+	RangeInterface<AnonymousIndex>()
+    {
+	mOrig.resize(sizeof...(RangeTypes));
+	RPackNum<sizeof...(RangeTypes)-1>::RangesToVec( origs, mOrig );
+	RPackNum<sizeof...(RangeTypes)-1>::getSize( origs );
+    }
+    
+    template <class... RangeTypes>
+    AnonymousRange::AnonymousRange(std::shared_ptr<RangeTypes>... origs) :
+	RangeInterface<AnonymousIndex>()
+    {
+	auto rst = std::make_tuple(origs...);
+	mOrig.resize(sizeof...(RangeTypes));
+	RPackNum<sizeof...(RangeTypes)-1>::RangesToVec( rst, mOrig );
+	RPackNum<sizeof...(RangeTypes)-1>::getSize( rst );
+    }
+    
+    size_t AnonymousRange::get(size_t pos) const
+    {
+	return pos;
+    }
+        
+    size_t AnonymousRange::getMeta(size_t metaPos) const
+    {
+	return metaPos;
+    }
+        
+    size_t AnonymousRange::size() const
+    {
+	return mSize;
+    }
+        
+    size_t AnonymousRange::dim() const
+    {
+	return 1;
+    }
+        
+    typename AnonymousRange::IndexType AnonymousRange::begin() const
+    {
+	AnonymousIndex i( std::dynamic_pointer_cast<AnonymousRange>
+			  ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
+	i = 0;
+	return i;
+    }
+    
+    typename AnonymousRange::IndexType AnonymousRange::end() const
+    {
+	AnonymousIndex i( std::dynamic_pointer_cast<AnonymousRange>
+			  ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
+	i = size();
+	return i;
+    }
+
+    // put this in the interface class !!!
+    std::shared_ptr<VIWB> AnonymousRange::index() const
+    {
+	typedef IndexWrapper<IndexType> IW;
+	return std::make_shared<IW>
+	    (std::make_shared<IndexType>
+	     ( std::dynamic_pointer_cast<AnonymousRange>
+	       ( std::shared_ptr<RangeBase>( RB::mThis ) ) ) );
+    }
+}
+
+#endif

src/include/ranges/container_range.h

@@ -0,0 +1,497 @@
+// -*- C++ -*-
+
+#ifndef __container_range_h__
+#define __container_range_h__
+
+#include <cstdlib>
+#include <tuple>
+#include <memory>
+
+//#include "base_def.h"
+#include "ranges/range_base.h"
+#include "ranges/index_base.h"
+
+#include "rpack_num.h"
+
+namespace MultiArrayTools
+{
+    
+    template <class... Indices>
+    class ContainerIndex : public IndexInterface<ContainerIndex<Indices...>,
+						 std::tuple<typename Indices::MetaType...> >
+    {
+    public:
+
+	typedef IndexInterface<ContainerIndex<Indices...>,
+			       std::tuple<typename Indices::MetaType...> > IB;
+	typedef std::tuple<typename Indices::MetaType...>  MetaType;
+	typedef std::tuple<std::shared_ptr<Indices>...> IndexPack;
+	typedef ContainerRange<typename Indices::RangeType...> RangeType;
+
+	static IndexType sType() { return IndexType::CONT; }
+	static size_t sDim() { return sizeof...(Indices); }
+	static size_t totalDim() { return mkTotalDim<Indices...>(); }
+	
+    private:
+
+	bool mExternControl = false;
+	IndexPack mIPack;
+	std::array<size_t,sizeof...(Indices)+1> mBlockSizes; 
+	
+    public:
+
+	ContainerIndex() = delete;
+	
+	template <class MRange>
+	ContainerIndex(const std::shared_ptr<MRange>& range);
+	
+	template <size_t N>
+	auto get() const -> decltype( *std::get<N>( mIPack ) )&;
+
+	template <size_t N>
+	auto getPtr() const -> decltype( std::get<N>( mIPack ) )&;
+
+	const IndexPack& pack() const { return mIPack; }
+	
+	ContainerIndex& sync(); // recalculate 'IB::mPos' when externalControl == true
+	ContainerIndex& operator()(const std::shared_ptr<Indices>&... inds); // control via external indices
+	ContainerIndex& operator()(); // -> sync; just to shorten the code
+
+	// ==== >>>>> STATIC POLYMORPHISM <<<<< ====
+	
+	IndexType type() const;
+	
+	ContainerIndex& operator++();
+	ContainerIndex& operator--();
+	
+	ContainerIndex& operator=(size_t pos);
+
+	int pp(std::intptr_t idxPtrNum);
+	int mm(std::intptr_t idxPtrNum);
+	
+	MetaType meta();
+	ContainerIndex& at(const MetaType& metaPos);
+
+	size_t dim();
+	bool first();
+	bool last();
+
+	std::shared_ptr<RangeType> range();
+	
+	template <size_t N>
+	auto getPtr() -> decltype( std::get<N>( mIPack ) )&;
+
+	std::shared_ptr<VIWB> getVPtr(size_t n);
+	size_t getStepSize(size_t n);
+
+	std::vector<IndexInfo> infoVec() const;
+	
+	std::string id() const;
+	void print(size_t offset);
+	
+	template <class Exprs>
+	auto ifor(Exprs exs) const
+	    -> decltype(RPackNum<sizeof...(Indices)-1>::mkFor(mIPack, exs));
+
+	template <class Exprs>
+	auto iforh(Exprs exs) const
+	    -> decltype(RPackNum<sizeof...(Indices)-1>::mkForh(mIPack, exs));
+
+    };
+
+    
+    template <class... Ranges>
+    class ContainerRangeFactory : public RangeFactoryBase
+    {
+    public:
+
+	typedef ContainerRange<Ranges...> oType;
+	
+	ContainerRangeFactory();
+	ContainerRangeFactory(const std::shared_ptr<Ranges>&... rs);
+	ContainerRangeFactory(const typename ContainerRange<Ranges...>::SpaceType& space);
+	
+	virtual std::shared_ptr<RangeBase> create() override; 
+	
+    protected:
+	
+    };
+    
+    template <class... Ranges>
+    class ContainerRange : public RangeInterface<ContainerIndex<typename Ranges::IndexType...> >
+    {
+    public:
+
+	typedef RangeBase RB;
+	typedef std::tuple<std::shared_ptr<Ranges>...> SpaceType;
+	typedef ContainerIndex<typename Ranges::IndexType...> IndexType;
+	//typedef typename RangeInterface<ContainerIndex<typename Ranges::IndexType...> >::IndexType IndexType;
+	
+    protected:
+	ContainerRange() = default;
+	ContainerRange(const ContainerRange& in) = delete;
+	ContainerRange& operator=(const ContainerRange& in) = delete;
+	
+	ContainerRange(const std::shared_ptr<Ranges>&... rs);
+	ContainerRange(const SpaceType& space);
+
+	SpaceType mSpace;	
+
+    public:	
+	static const size_t sdim = sizeof...(Ranges);
+
+	virtual size_t dim() const override;
+	virtual size_t size() const override;
+
+	template <size_t N>
+	auto get() const -> decltype( *std::get<N>( mSpace ) )&;
+
+	template <size_t N>
+	auto getPtr() const -> decltype( std::get<N>( mSpace ) )&;
+
+	const SpaceType& space() const;
+	
+	virtual IndexType begin() const override;
+	virtual IndexType end() const override;
+	
+	virtual std::shared_ptr<VIWB> index() const override;
+
+	friend ContainerRangeFactory<Ranges...>;
+
+	static constexpr bool defaultable = false;
+	static constexpr size_t ISSTATIC = SubProp<Ranges...>::ISSTATIC;
+	static constexpr size_t SIZE = SubProp<Ranges...>::SIZE;
+    };
+   
+} // end namespace MultiArrayTools
+
+/* ========================= *
+ * ---   TEMPLATE CODE   --- *
+ * ========================= */
+
+namespace MultiArrayTools
+{
+    namespace
+    {
+	using namespace MultiArrayHelper;
+    }
+    
+    /**********************
+     *   ContainerIndex   *
+     **********************/
+
+    template <class... Indices>
+    template <class MRange>
+    ContainerIndex<Indices...>::ContainerIndex(const std::shared_ptr<MRange>& range) :
+	IndexInterface<ContainerIndex<Indices...>,std::tuple<typename Indices::MetaType...> >(range, 0)
+    {
+	RPackNum<sizeof...(Indices)-1>::construct(mIPack, *range);
+	IB::mPos = RPackNum<sizeof...(Indices)-1>::makePos(mIPack);
+	std::get<sizeof...(Indices)>(mBlockSizes) = 1;
+	RPackNum<sizeof...(Indices)-1>::initBlockSizes(mBlockSizes, mIPack);
+    }
+
+    template <class... Indices>
+    ContainerIndex<Indices...>& ContainerIndex<Indices...>::sync()
+    {
+	if(mExternControl){
+	    IB::mPos = RPackNum<sizeof...(Indices)-1>::makePos(mIPack);
+	    //VCHECK(id());
+	    //VCHECK(sizeof...(Indices));
+	    //assert(IB::mPos < IB::max());
+	}
+	return *this;
+    }
+    
+    template <class... Indices>
+    template <size_t N>
+    auto ContainerIndex<Indices...>::get() const -> decltype( *std::get<N>( mIPack ) )&
+    {
+	return *std::get<N>( mIPack );
+    }
+
+    template <class... Indices>
+    template <size_t N>
+    auto ContainerIndex<Indices...>::getPtr() const -> decltype( std::get<N>( mIPack ) )&
+    {
+	return std::get<N>( mIPack );
+    }
+
+    template <class... Indices>
+    ContainerIndex<Indices...>& ContainerIndex<Indices...>::operator()(const std::shared_ptr<Indices>&... inds)
+    {
+	RPackNum<sizeof...(Indices)-1>::swapIndices(mIPack, inds...);
+	mExternControl = true;
+	return sync();
+    }    
+
+    template <class... Indices>
+    ContainerIndex<Indices...>& ContainerIndex<Indices...>::operator()()
+    {
+	return sync();
+    }
+
+    template <class... Indices>
+    IndexType ContainerIndex<Indices...>::type() const { return IndexType::CONT; }
+
+    template <class... Indices>
+    ContainerIndex<Indices...>& ContainerIndex<Indices...>::operator++()
+    {
+	if(mExternControl){
+	    IB::mPos = RPackNum<sizeof...(Indices)-1>::makePos(mIPack);
+	}
+	RPackNum<sizeof...(Indices)-1>::pp( mIPack );
+	++IB::mPos;
+	return *this;
+    }
+
+    template <class... Indices>
+    ContainerIndex<Indices...>& ContainerIndex<Indices...>::operator--()
+    {
+	if(mExternControl){
+	    IB::mPos = RPackNum<sizeof...(Indices)-1>::makePos(mIPack);
+	}
+	RPackNum<sizeof...(Indices)-1>::mm( mIPack );
+	--IB::mPos;
+	return *this;
+
+    }
+
+    template <class... Indices>
+    ContainerIndex<Indices...>& ContainerIndex<Indices...>::operator=(size_t pos)
+    {
+	IB::mPos = pos;
+	RPackNum<sizeof...(Indices)-1>::setIndexPack(mIPack, pos);
+	return *this;
+    }
+
+    template <class... Indices>
+    int ContainerIndex<Indices...>::pp(std::intptr_t idxPtrNum)
+    {
+	int tmp = RPackNum<sizeof...(Indices)-1>::pp(mIPack, mBlockSizes, idxPtrNum);
+	IB::mPos += tmp;
+	return tmp;
+    }
+
+    template <class... Indices>
+    int ContainerIndex<Indices...>::mm(std::intptr_t idxPtrNum)
+    {
+	int tmp = RPackNum<sizeof...(Indices)-1>::mm(mIPack, mBlockSizes, idxPtrNum);
+	IB::mPos -= tmp;
+	return tmp;
+    }
+
+    template <class... Indices>
+    typename ContainerIndex<Indices...>::MetaType ContainerIndex<Indices...>::meta()
+    {
+	MetaType metaTuple;
+	RPackNum<sizeof...(Indices)-1>::getMetaPos(metaTuple, mIPack);
+	return metaTuple;
+    }
+
+    template <class... Indices>
+    ContainerIndex<Indices...>& ContainerIndex<Indices...>::at(const MetaType& metaPos)
+    {
+	RPackNum<sizeof...(Indices)-1>::setMeta(mIPack, metaPos);
+	IB::mPos = RPackNum<sizeof...(Indices)-1>::makePos(mIPack);
+	return *this;
+    }
+
+    template <class... Indices>
+    size_t ContainerIndex<Indices...>::dim()
+    {
+	return sizeof...(Indices);
+    }
+
+    template <class... Indices>
+    bool ContainerIndex<Indices...>::first()
+    {
+	return IB::pos() == 0;
+    }
+
+    template <class... Indices>
+    bool ContainerIndex<Indices...>::last()
+    {
+	return IB::pos() == IB::mMax - 1;
+    }
+
+    template <class... Indices>
+    std::shared_ptr<typename ContainerIndex<Indices...>::RangeType>
+    ContainerIndex<Indices...>::range()
+    {
+	return std::dynamic_pointer_cast<RangeType>( IB::mRangePtr );
+    }
+
+    template <class... Indices>
+    template <size_t N>
+    auto ContainerIndex<Indices...>::getPtr() -> decltype( std::get<N>( mIPack ) )&
+    {
+	return std::get<N>( mIPack );
+    }
+
+    template <class... Indices>
+    std::shared_ptr<VIWB> ContainerIndex<Indices...>::getVPtr(size_t n)
+    {
+	if(n >= sizeof...(Indices)){
+	    assert(0);
+	    // throw !!
+	}
+	ContainerIndex<Indices...> const* t = this;
+	return RPackNum<sizeof...(Indices)-1>::getIndexPtr(*t, n);
+    }
+
+    template <class... Indices>
+    size_t ContainerIndex<Indices...>::getStepSize(size_t n)
+    {
+	if(n >= sizeof...(Indices)){
+	    assert(0);
+	    // throw !!
+	}
+	return mBlockSizes[n+1];
+    }
+
+    template <class... Indices>
+    std::vector<IndexInfo> ContainerIndex<Indices...>::infoVec() const
+    {
+	std::vector<IndexInfo> out;
+	out.reserve(sizeof...(Indices));
+	RPackNum<sizeof...(Indices)-1>::buildInfoVec(out, mIPack, mBlockSizes);
+	return std::move( out );
+    }
+    
+    template <class... Indices>
+    std::string ContainerIndex<Indices...>::id() const
+    {
+	return std::string("con") + std::to_string(IB::mId);
+    }
+
+    template <class... Indices>
+    void ContainerIndex<Indices...>::print(size_t offset)
+    {
+	if(offset == 0){
+	    std::cout << " === " << std::endl;
+	}
+	for(size_t j = 0; j != offset; ++j) { std::cout << "\t"; }
+	std::cout << id() << "[" << reinterpret_cast<std::intptr_t>(this) << "]"
+		  << "(" << IB::mRangePtr << "): " << meta() << std::endl;
+	RPackNum<sizeof...(Indices)-1>::printIndex(mIPack, offset+1);
+    }
+
+    template <class... Indices>
+    template <class Exprs>
+    auto ContainerIndex<Indices...>::ifor(Exprs exs) const
+	-> decltype(RPackNum<sizeof...(Indices)-1>::mkFor(mIPack, exs))
+    {
+	return RPackNum<sizeof...(Indices)-1>::mkFor(mIPack, exs);
+    }
+
+    template <class... Indices>
+    template <class Exprs>
+    auto ContainerIndex<Indices...>::iforh(Exprs exs) const
+	-> decltype(RPackNum<sizeof...(Indices)-1>::mkForh(mIPack, exs))
+    {
+	return RPackNum<sizeof...(Indices)-1>::mkForh(mIPack, exs);
+    }
+
+    
+    /*****************************
+     *   ContainerRangeFactory   *
+     *****************************/
+
+    template <class... Ranges>
+    ContainerRangeFactory<Ranges...>::ContainerRangeFactory(const std::shared_ptr<Ranges>&... rs)
+    {
+	mProd = std::shared_ptr<ContainerRange<Ranges...> >( new ContainerRange<Ranges...>( rs... ) );
+    }
+   
+    template <class... Ranges>
+    ContainerRangeFactory<Ranges...>::
+    ContainerRangeFactory(const typename ContainerRange<Ranges...>::SpaceType& space)
+    {
+	mProd = std::shared_ptr<ContainerRange<Ranges...> >( new ContainerRange<Ranges...>( space ) );
+    }
+
+    template <class... Ranges>
+    std::shared_ptr<RangeBase> ContainerRangeFactory<Ranges...>::create()
+    {
+	setSelf();
+	return mProd;
+    }
+    
+    /**********************
+     *   ContainerRange   *
+     **********************/
+
+    template <class... Ranges>
+    ContainerRange<Ranges...>::ContainerRange(const std::shared_ptr<Ranges>&... rs) :
+	mSpace( std::make_tuple( rs... ) ) {}
+
+    template <class... Ranges>
+    ContainerRange<Ranges...>::ContainerRange(const SpaceType& space) : mSpace( space ) {}
+
+    template <class... Ranges>
+    size_t ContainerRange<Ranges...>::dim() const
+    {
+	return sizeof...(Ranges);
+    }
+
+    template <class... Ranges>
+    size_t ContainerRange<Ranges...>::size() const
+    {
+	return RPackNum<sizeof...(Ranges)-1>::getSize(mSpace);
+    }
+
+    template <class... Ranges>
+    template <size_t N>
+    auto ContainerRange<Ranges...>::get() const -> decltype( *std::get<N>( mSpace ) )&
+    {
+	return *std::get<N>( mSpace );
+    }
+
+    template <class... Ranges>
+    template <size_t N>
+    auto ContainerRange<Ranges...>::getPtr() const -> decltype( std::get<N>( mSpace ) )&
+    {
+	return std::get<N>( mSpace );
+    }
+
+    template <class... Ranges>
+    const typename ContainerRange<Ranges...>::SpaceType& ContainerRange<Ranges...>::space() const
+    {
+	return mSpace;
+    }
+        
+    template <class... Ranges>
+    typename ContainerRange<Ranges...>::IndexType ContainerRange<Ranges...>::begin() const
+    {
+	ContainerIndex<typename Ranges::IndexType...>
+	    i( std::dynamic_pointer_cast<ContainerRange<Ranges...> >
+	       ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
+	i = 0;
+	return i;
+    }
+
+    template <class... Ranges>
+    typename ContainerRange<Ranges...>::IndexType ContainerRange<Ranges...>::end() const
+    {
+	ContainerIndex<typename Ranges::IndexType...>
+	    i( std::dynamic_pointer_cast<ContainerRange<Ranges...> >
+	       ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
+	i = size();
+	return i;
+    }
+    
+    template <class... Ranges>
+    std::shared_ptr<VIWB> ContainerRange<Ranges...>::index() const
+    {
+	typedef IndexWrapper<IndexType> IW;
+	return std::make_shared<IW>
+	    ( std::make_shared<IndexType>
+	      ( std::dynamic_pointer_cast<ContainerRange<Ranges...> >
+		( std::shared_ptr<RangeBase>( RB::mThis ) ) ) );
+    }
+    
+} // end namespace MultiArrayTools
+
+
+#endif

src/include/ranges/index_base.h

@@ -0,0 +1,157 @@
+// -*- C++ -*-
+
+#ifndef __index_base_h__
+#define __index_base_h__
+
+#include <cstdlib>
+#include <string>
+#include <vector>
+
+#include "rbase_def.h"
+#include "range_base.h"
+#include "index_type.h"
+#include "vindex_wrapper.h"
+#include "index_info.h"
+
+#include "xfor/xfor.h"
+
+namespace MultiArrayTools
+{
+    
+    template <class I, typename MetaType>
+    class IndexInterface
+    {
+    public:
+	//typedef typename I::RangeType RangeType;
+	
+	//DEFAULT_MEMBERS(IndexInterface);
+
+	I& THIS() { return static_cast<I&>(*this); }
+	I const& THIS() const { return static_cast<I const&>(*this); }
+	
+	~IndexInterface() = default; 
+
+	IndexType type() const { return THIS().type(); }
+	
+	I& operator=(size_t pos) { return THIS() = pos; }
+	I& operator++() { return THIS()++; }
+	I& operator--() { return THIS()--;}
+
+	int pp(std::intptr_t idxPtrNum) { return THIS().pp(idxPtrNum); }
+	int mm(std::intptr_t idxPtrNum) { return THIS().mm(idxPtrNum); }
+	
+	bool operator==(const IndexInterface& in) const;
+	bool operator!=(const IndexInterface& in) const;
+	
+	size_t dim() const { return THIS().dim(); }
+	size_t pos() const;
+	size_t max() const;
+
+	bool last() const { return THIS().last(); }
+	bool first() const { return THIS().first(); }
+
+	
+	std::shared_ptr<RangeBase> vrange() const { return mRangePtr; }
+	/*auto range() const -> decltype( I::S_range(THIS()) ) { return I::S_range(THIS()); }
+	
+	template <size_t N>
+	auto getPtr() const -> decltype(I::template S_get<N>(THIS()))
+	{ return I::template S_get<N>(THIS()); }
+	*/
+	std::shared_ptr<VIWB> getVPtr(size_t n) const { return THIS().getVPtr(n); }
+
+	std::vector<IndexInfo> infoVec() const { return THIS().infoVec(); }
+	
+	size_t getStepSize(size_t n) const { return THIS().getStepSize(n); }
+	
+	operator size_t() const;
+
+	std::string id() const { return THIS().id(); }
+
+	MetaType meta() const { return THIS().meta(); }
+	I& at(const MetaType& meta) { return THIS().at(meta); }
+
+	void print(size_t offset = 0) const { THIS().print(offset); }
+
+	IndexInfo info() const { return IndexInfo(THIS()); }
+	
+	// CHECK / IMPLEMENT !!!!!!
+	template <class Expr>
+	auto ifor(const Expr ex) const -> decltype(THIS().template ifor<Expr>(ex))
+	{ return THIS().template ifor<Expr>(ex); }
+
+	template <class Expr>
+	auto iforh(const Expr ex) const -> decltype(THIS().template iforh<Expr>(ex))
+	{ return THIS().template iforh<Expr>(ex); }
+
+    private:
+
+	friend I;
+	
+	IndexInterface() { mId = indexId(); }
+	IndexInterface(const IndexInterface& in) = default;
+	IndexInterface& operator=(const IndexInterface& in) = default;
+	IndexInterface(IndexInterface&& in) = default;
+	IndexInterface& operator=(IndexInterface&& in) = default;
+	
+	IndexInterface(const std::shared_ptr<RangeBase>& range, size_t pos);
+	
+	std::shared_ptr<RangeBase> mRangePtr;
+	size_t mPos;
+	size_t mId;
+	size_t mMax;
+    };
+   
+}
+
+/* ========================= *
+ * ---   TEMPLATE CODE   --- *
+ * ========================= */
+
+namespace MultiArrayTools
+{
+    /**********************
+     *   IndexInterface   *
+     **********************/
+
+    template <class I, typename MetaType>
+    IndexInterface<I,MetaType>::IndexInterface(const std::shared_ptr<RangeBase>& range,
+					       size_t pos) : mRangePtr(range),
+							     mPos(pos),
+							     mMax(mRangePtr->size())
+    {
+	mId = indexId();
+    }
+
+    template <class I, typename MetaType>
+    bool IndexInterface<I,MetaType>::operator==(const IndexInterface& in) const
+    {
+	return in.mPos == mPos and in.mRangePtr.get() == mRangePtr.get();
+    }
+
+    template <class I, typename MetaType>
+    bool IndexInterface<I,MetaType>::operator!=(const IndexInterface& in) const
+    {
+	return in.mPos != mPos or in.mRangePtr.get() != mRangePtr.get();
+    }
+
+    template <class I, typename MetaType>
+    size_t IndexInterface<I,MetaType>::pos() const
+    {
+	return mPos;
+    }
+
+    template <class I, typename MetaType>
+    size_t IndexInterface<I,MetaType>::max() const
+    {
+	return mMax;
+    }
+
+    template <class I, typename MetaType>
+    IndexInterface<I,MetaType>::operator size_t() const
+    {
+	return pos();
+    }
+}
+
+#endif

src/include/ranges/index_info.h

@@ -0,0 +1,103 @@
+// -*- C++ -*-
+
+#ifndef __index_info_h__
+#define __index_info_h__
+
+#include <cstdlib>
+#include <cstdint>
+#include <vector>
+#include <memory>
+#include "vindex_base.h"
+#include "index_type.h"
+
+namespace MultiArrayTools
+{
+    class IndexInfo;
+    
+    class IndexInfo
+    {
+    public:
+
+	IndexInfo(IndexInfo&& in) = default;
+	IndexInfo& operator=(IndexInfo&& in) = default;
+	IndexInfo(const IndexInfo& in) = default;
+	IndexInfo& operator=(const IndexInfo& in) = default;
+	
+	template <class IndexClass>
+	IndexInfo(const IndexClass& ind, size_t stepSize = 1);
+
+	template <class IndexClass>
+	IndexInfo& reassign(const IndexClass& ind, size_t stepSize = 1);
+	
+	bool operator==(const IndexInfo& in) const;
+	bool operator!=(const IndexInfo& in) const;
+
+	bool operator<=(const IndexInfo& in) const;
+	bool operator<(const IndexInfo& in) const;
+	bool operator>(const IndexInfo& in) const;
+	bool operator>=(const IndexInfo& in) const;
+
+	
+	const IndexInfo* getPtr(size_t inum) const;
+	std::intptr_t getPtrNum() const;
+	size_t dim() const;
+	size_t max() const;
+	size_t getStepSize(size_t inum) const;
+	size_t getStepSize() const;
+	IndexType type() const;
+	
+    private:
+
+	IndexInfo() = default;
+	
+	std::vector<IndexInfo> mNext;
+	std::intptr_t mPtrNum;
+	size_t mDim;
+	size_t mMax;
+	size_t mStepSize;
+	IndexType mType;
+    };
+
+    template <class IndexClass>
+    IndexInfo::IndexInfo(const IndexClass& ind, size_t stepSize) :
+	mNext(ind.infoVec()),
+	mPtrNum( reinterpret_cast<std::intptr_t>( &ind ) ),
+	mDim(ind.vrange()->dim()),
+	mMax(ind.max()),
+	mStepSize(stepSize),
+	mType(ind.type())
+    {}
+
+    template <class IndexClass>
+    IndexInfo& IndexInfo::reassign(const IndexClass& ind, size_t stepSize)
+    {
+	IndexInfo ii(ind, stepSize);
+	(*this) = std::move(ii);
+	return *this;
+    }
+
+    std::vector<IndexInfo> getRootIndices(const IndexInfo& info);
+
+    inline size_t getStepSize(const IndexInfo& ii, std::intptr_t j)
+    {
+	if(ii.type() == IndexType::SINGLE){
+	    return ii.getPtrNum() == j ? 1 : 0;	    
+	}
+	else {
+	    size_t ss = 0;
+	    size_t sx = 1;
+	    for(size_t i = 0; i != ii.dim(); ++i){
+		const IndexInfo& itmp = *ii.getPtr(ii.dim()-i-1);
+		const size_t max = itmp.max();
+		const size_t tmp = getStepSize(itmp, j);
+		ss += tmp * sx;
+		sx *= max;
+	    }
+	    return ss;
+	}
+    }
+
+    size_t getStepSize(const std::vector<IndexInfo>& iv, std::intptr_t j);
+} // end namespace MultiArrayTools
+
+#endif

src/include/ranges/index_type.h

@@ -0,0 +1,14 @@
+
+#ifndef __index_type_h__
+#define __index_type_h__
+
+namespace MultiArrayTools
+{
+    enum class IndexType{
+	SINGLE = 0,
+	MULTI = 1,
+	CONT = 2	
+    };    
+}
+
+#endif

src/include/ranges/multi_range.h

@@ -0,0 +1,543 @@
+// -*- C++ -*-
+
+#ifndef __multi_range_h__
+#define __multi_range_h__
+
+#include <cstdlib>
+#include <tuple>
+#include <memory>
+
+//#include "base_def.h"
+#include "ranges/range_base.h"
+#include "ranges/index_base.h"
+
+#include "rpack_num.h"
+
+namespace MultiArrayTools
+{
+    namespace
+    {
+	using namespace MultiArrayHelper;
+    }
+    
+    template <class... Indices>
+    class MultiIndex : public IndexInterface<MultiIndex<Indices...>,
+					     std::tuple<typename Indices::MetaType...> >
+    {
+    public:
+	
+	typedef IndexInterface<MultiIndex<Indices...>,
+			       std::tuple<typename Indices::MetaType...> > IB;
+	typedef std::tuple<std::shared_ptr<Indices>...> IndexPack;
+	typedef std::tuple<typename Indices::MetaType...> MetaType;
+	typedef MultiRange<typename Indices::RangeType...> RangeType;
+
+	static IndexType sType() { return IndexType::MULTI; }
+	static size_t sDim() { return sizeof...(Indices); }
+	static size_t totalDim() { return mkTotalDim<Indices...>(); }
+	
+    private:
+	
+	IndexPack mIPack;
+	std::array<size_t,sizeof...(Indices)+1> mBlockSizes;
+
+    public:
+
+	const IndexPack& pack() const { return mIPack; }
+	
+	MultiIndex() = delete;
+	
+	// NO DEFAULT HERE !!!
+	// ( have to assign sub-indices (ptr!) correctly )
+	//MultiIndex(const MultiIndex& in);
+	//MultiIndex& operator=(const MultiIndex& in);
+	MultiIndex& operator=(ContainerIndex<Indices...>& ci);
+
+	template <class MRange>
+	MultiIndex(const std::shared_ptr<MRange>& range);
+
+	template <size_t DIR>
+	MultiIndex& up();
+
+	template <size_t DIR>
+	MultiIndex& down();
+
+	template <size_t N>
+	auto get() const -> decltype( *std::get<N>( mIPack ) )&;
+
+	template <size_t N>
+	auto getPtr() const -> decltype( std::get<N>( mIPack ) )&;
+	
+	// raplace instances (in contrast to its analogon in ContainerIndex
+	// MultiIndices CANNOT be influences be its subindices, so there is
+	// NO foreign/external controll)
+	// Do NOT share index instances between two or more MultiIndex instances
+	MultiIndex& operator()(std::shared_ptr<Indices>&... indices);
+
+	// ==== >>>>> STATIC POLYMORPHISM <<<<< ====
+	
+	IndexType type() const;
+
+	MultiIndex& operator=(size_t pos);
+
+	MultiIndex& operator++();
+	MultiIndex& operator--();
+
+	int pp(std::intptr_t idxPtrNum);
+	int mm(std::intptr_t idxPtrNum);
+
+	MetaType meta();
+	MultiIndex& at(const MetaType& metaPos);
+	
+	size_t dim();
+	bool first();
+	bool last();
+	std::shared_ptr<RangeType> range();
+
+	template <size_t N>
+	auto getPtr() -> decltype( std::get<N>( mIPack ) )&;
+	
+	std::shared_ptr<VIWB> getVPtr(size_t n);
+	size_t getStepSize(size_t n);
+
+	std::vector<IndexInfo> infoVec() const;
+	
+	std::string id() const;
+	void print(size_t offset);
+
+	template <class Exprs>
+	auto ifor(Exprs exs) const
+	    -> decltype(RPackNum<sizeof...(Indices)-1>::mkFor(mIPack, exs));
+
+	template <class Exprs>
+	auto iforh(Exprs exs) const
+	    -> decltype(RPackNum<sizeof...(Indices)-1>::mkForh(mIPack, exs));
+
+    };
+
+    /*************************
+     *   MultiRangeFactory   *
+     *************************/
+
+    template <class... Ranges>
+    class MultiRangeFactory : public RangeFactoryBase
+    {
+    public:
+	typedef MultiRange<Ranges...> oType;
+	
+	MultiRangeFactory() = delete;
+	MultiRangeFactory(const std::shared_ptr<Ranges>&... rs);
+	MultiRangeFactory(const typename MultiRange<Ranges...>::SpaceType& space);
+	MultiRangeFactory(const std::shared_ptr<ContainerRange<Ranges...> >& cr);
+	
+	virtual std::shared_ptr<RangeBase> create() override;
+    };
+    
+    /******************
+     *   MultiRange   *
+     ******************/
+      
+    template <class... Ranges>
+    class MultiRange : public RangeInterface<MultiIndex<typename Ranges::IndexType...> >
+    {
+    public:
+	typedef RangeBase RB;
+	typedef std::tuple<std::shared_ptr<Ranges>...> SpaceType;
+	typedef MultiIndex<typename Ranges::IndexType...> IndexType;	
+	//typedef typename RangeInterface<MultiIndex<typename Ranges::IndexType...> >::IndexType IndexType;
+
+    protected:
+	MultiRange() = delete;
+	MultiRange(const MultiRange& in) = delete;
+	MultiRange& operator=(const MultiRange& in) = delete;
+	
+	MultiRange(const std::shared_ptr<Ranges>&... rs);
+	MultiRange(const SpaceType& space);
+	
+	SpaceType mSpace;
+
+    public:
+	
+	static const size_t sdim = sizeof...(Ranges);
+
+	template <size_t N>
+	auto get() const -> decltype( *std::get<N>( mSpace ) )&;
+
+	template <size_t N>
+	auto getPtr() const -> decltype( std::get<N>( mSpace ) )&;
+
+	virtual size_t dim() const override;
+	virtual size_t size() const override;
+	
+	const SpaceType& space() const;
+	
+	virtual IndexType begin() const override;
+	virtual IndexType end() const override;
+
+	virtual std::shared_ptr<VIWB> index() const override;
+
+	friend MultiRangeFactory<Ranges...>;
+
+	static constexpr bool defaultable = false;
+	static constexpr size_t ISSTATIC = SubProp<Ranges...>::ISSTATIC;
+	static constexpr size_t SIZE = SubProp<Ranges...>::SIZE;
+    };
+    
+}
+
+/* ========================= *
+ * ---   TEMPLATE CODE   --- *
+ * ========================= */
+
+namespace MultiArrayTools
+{
+
+    namespace
+    {
+	using namespace MultiArrayHelper;
+    }
+
+    /******************
+     *   MultiIndex   *
+     ******************/
+
+    /*
+    template <class... Indices>
+    MultiIndex<Indices...>::MultiIndex(const MultiIndex<Indices...>& in) :
+	IndexInterface<std::tuple<typename Indices::MetaType...> >(in)
+    {
+	RPackNum<sizeof...(Indices)-1>::copy(mIPack, in);
+	IB::mPos = RPackNum<sizeof...(Indices)-1>::makePos(mIPack);
+    }
+	
+    template <class... Indices>
+    MultiIndex<Indices...>& MultiIndex<Indices...>::operator=(const MultiIndex<Indices...>& in)
+    {
+	IndexI::operator=(in);
+	RPackNum<sizeof...(Indices)-1>::copy(mIPack, in);
+	IB::mPos = RPackNum<sizeof...(Indices)-1>::makePos(mIPack);
+	return *this;
+    }
+    */
+    template <class... Indices>
+    MultiIndex<Indices...>& MultiIndex<Indices...>::operator=(ContainerIndex<Indices...>& ci)
+    {	
+	RPackNum<sizeof...(Indices)-1>::copyInst(mIPack, ci);
+	IB::mPos = RPackNum<sizeof...(Indices)-1>::makePos(mIPack);
+	return *this;
+    }
+    
+    template <class... Indices>
+    template <class MRange>
+    MultiIndex<Indices...>::MultiIndex(const std::shared_ptr<MRange>& range) :
+	IndexInterface<MultiIndex<Indices...>,std::tuple<typename Indices::MetaType...> >(range, 0)
+    {
+	RPackNum<sizeof...(Indices)-1>::construct(mIPack, *range);
+	IB::mPos = RPackNum<sizeof...(Indices)-1>::makePos(mIPack);
+	std::get<sizeof...(Indices)>(mBlockSizes) = 1;
+	RPackNum<sizeof...(Indices)-1>::initBlockSizes(mBlockSizes, mIPack); // has one more element!
+    }
+
+    template <class... Indices>
+    template <size_t DIR>
+    MultiIndex<Indices...>& MultiIndex<Indices...>::up()
+    {
+	static_assert(DIR < sizeof...(Indices), "DIR exceeds number of sub-indices");
+	IB::mPos += RPackNum<sizeof...(Indices)-DIR-1>::blockSize( mIPack );
+	RPackNum<DIR>::pp( mIPack );
+	return *this;
+    }
+
+    template <class... Indices>
+    template <size_t DIR>
+    MultiIndex<Indices...>& MultiIndex<Indices...>::down()
+    {
+	static_assert(DIR < sizeof...(Indices), "DIR exceeds number of sub-indices");
+	IB::mPos -= RPackNum<sizeof...(Indices)-DIR-1>::blockSize( mIPack );
+	RPackNum<DIR>::mm( mIPack );
+	return *this;
+    }
+    
+    template <class... Indices>
+    template <size_t N>
+    auto MultiIndex<Indices...>::get() const -> decltype( *std::get<N>( mIPack ) )&
+    {
+	return *std::get<N>(mIPack);
+    }
+
+    template <class... Indices>
+    template <size_t N>
+    auto MultiIndex<Indices...>::getPtr() const -> decltype( std::get<N>( mIPack ) )&
+    {
+	return std::get<N>(mIPack);
+    }
+    
+    template <class... Indices>
+    MultiIndex<Indices...>& MultiIndex<Indices...>::operator()(std::shared_ptr<Indices>&... indices)
+    {
+	RPackNum<sizeof...(Indices)-1>::swapIndices(mIPack, indices...);
+	RPackNum<sizeof...(Indices)-1>::setIndexPack(mIPack, IB::mPos);
+	return *this;
+    }
+
+    template <class... Indices>
+    IndexType MultiIndex<Indices...>::type() const
+    {
+	return IndexType::MULTI;
+    }
+
+    template <class... Indices>
+    MultiIndex<Indices...>& MultiIndex<Indices...>::operator=(size_t pos)
+    {
+	IB::mPos = pos;
+	RPackNum<sizeof...(Indices)-1>::setIndexPack(mIPack, pos);
+	return *this;
+    }
+
+    template <class... Indices>
+    MultiIndex<Indices...>& MultiIndex<Indices...>::operator++()
+    {
+	RPackNum<sizeof...(Indices)-1>::pp( mIPack );
+	++IB::mPos;
+	return *this;
+    }
+
+    template <class... Indices>
+    MultiIndex<Indices...>& MultiIndex<Indices...>::operator--()
+    {
+	RPackNum<sizeof...(Indices)-1>::mm( mIPack );
+	--IB::mPos;
+	return *this;
+    }
+
+    template <class... Indices>
+    int MultiIndex<Indices...>::pp(std::intptr_t idxPtrNum)
+    {
+	int tmp = RPackNum<sizeof...(Indices)-1>::pp(mIPack, mBlockSizes, idxPtrNum);
+	IB::mPos += tmp;
+	return tmp;
+    }
+
+    template <class... Indices>
+    int MultiIndex<Indices...>::mm(std::intptr_t idxPtrNum)
+    {
+	int tmp = RPackNum<sizeof...(Indices)-1>::mm(mIPack, mBlockSizes, idxPtrNum);
+	IB::mPos -= tmp;
+	return tmp;
+    }
+
+    template <class... Indices>
+    typename MultiIndex<Indices...>::MetaType MultiIndex<Indices...>::meta()
+    {
+	MetaType metaTuple;
+	RPackNum<sizeof...(Indices)-1>::getMetaPos(metaTuple, mIPack);
+	return metaTuple;
+    }
+
+    template <class... Indices>
+    MultiIndex<Indices...>& MultiIndex<Indices...>::at(const MetaType& metaPos)
+    {
+	RPackNum<sizeof...(Indices)-1>::setMeta(mIPack, metaPos);
+	IB::mPos = RPackNum<sizeof...(Indices)-1>::makePos(mIPack);
+	return *this;
+    }
+
+    template <class... Indices>
+    size_t MultiIndex<Indices...>::dim()
+    {
+	return sizeof...(Indices);
+    }
+
+    template <class... Indices>
+    bool MultiIndex<Indices...>::first()
+    {
+	return IB::mPos == 0;
+    }
+
+    template <class... Indices>
+    bool MultiIndex<Indices...>::last()
+    {
+	return IB::mPos == IB::mMax - 1;
+    }
+
+    template <class... Indices>
+    std::shared_ptr<typename MultiIndex<Indices...>::RangeType>
+    MultiIndex<Indices...>::range()
+    {
+	return std::dynamic_pointer_cast<RangeType>( IB::mRangePtr );
+    }
+
+    template <class... Indices>
+    template <size_t N>
+    auto MultiIndex<Indices...>::getPtr() -> decltype( std::get<N>( mIPack ) )&
+    {
+	return std::get<N>(mIPack);
+    }
+	
+    template <class... Indices>
+    std::shared_ptr<VIWB> MultiIndex<Indices...>::getVPtr(size_t n)
+    {
+	if(n >= sizeof...(Indices)){
+	    assert(0);
+	    // throw !!
+	}
+	MultiIndex<Indices...> const* t = this;
+	return RPackNum<sizeof...(Indices)-1>::getIndexPtr(*t, n);
+    }
+    
+    template <class... Indices>	
+    size_t MultiIndex<Indices...>::getStepSize(size_t n)
+    {
+	if(n >= sizeof...(Indices)){
+	    assert(0);
+	    // throw !!
+	}
+	return mBlockSizes[n+1];
+    }
+
+    template <class... Indices>
+    std::vector<IndexInfo> MultiIndex<Indices...>::infoVec() const
+    {
+	std::vector<IndexInfo> out;
+	out.reserve(sizeof...(Indices));
+	RPackNum<sizeof...(Indices)-1>::buildInfoVec(out, mIPack, mBlockSizes);
+	return std::move( out );
+    }
+    
+    template <class... Indices>
+    std::string MultiIndex<Indices...>::id() const
+    {
+	return std::string("mul") + std::to_string(IB::mId);
+    }
+
+    template <class... Indices>
+    void MultiIndex<Indices...>::print(size_t offset)
+    {
+	if(offset == 0){
+	    std::cout << " === " << std::endl;
+	}
+	for(size_t j = 0; j != offset; ++j) { std::cout << "\t"; }
+	std::cout << id() << "[" << reinterpret_cast<std::intptr_t>(this)
+		  << "]" << "(" << IB::mRangePtr << "): " << meta() << std::endl;
+	RPackNum<sizeof...(Indices)-1>::printIndex(mIPack, offset+1);
+    }
+
+    template <class... Indices>
+    template <class Exprs>
+    auto MultiIndex<Indices...>::ifor(Exprs exs) const
+	-> decltype(RPackNum<sizeof...(Indices)-1>::mkFor(mIPack, exs))
+    {
+	return RPackNum<sizeof...(Indices)-1>::mkFor(mIPack, exs);
+    }
+
+    template <class... Indices>
+    template <class Exprs>
+    auto MultiIndex<Indices...>::iforh(Exprs exs) const
+	-> decltype(RPackNum<sizeof...(Indices)-1>::mkForh(mIPack, exs))
+    {
+	return RPackNum<sizeof...(Indices)-1>::mkForh(mIPack, exs);
+    }
+
+    /*************************
+     *   MultiRangeFactory   *
+     *************************/
+
+    template <class... Ranges>
+    MultiRangeFactory<Ranges...>::MultiRangeFactory(const std::shared_ptr<Ranges>&... rs)
+    {
+	mProd = std::shared_ptr< MultiRange<Ranges...> >( new MultiRange<Ranges...>( rs... ) );
+    }
+    
+    template <class... Ranges>
+    MultiRangeFactory<Ranges...>::MultiRangeFactory(const typename MultiRange<Ranges...>::SpaceType& st)
+    {
+	mProd = std::shared_ptr< MultiRange<Ranges...> >( new MultiRange<Ranges...>( st ) );
+    }
+
+    template <class... Ranges>
+    MultiRangeFactory<Ranges...>::MultiRangeFactory(const std::shared_ptr<ContainerRange<Ranges...> >& cr)
+    {
+	mProd = std::shared_ptr< MultiRange<Ranges...> >( new MultiRange<Ranges...>( cr->space() ) );
+    }
+    
+    template <class... Ranges>
+    std::shared_ptr<RangeBase> MultiRangeFactory<Ranges...>::create()
+    {
+	setSelf();
+	return mProd;
+    }
+    
+    /******************
+     *   MultiRange   *
+     ******************/
+
+    template <class... Ranges>
+    MultiRange<Ranges...>::MultiRange(const std::shared_ptr<Ranges>&... rs) : mSpace(std::make_tuple(rs...)) {}
+
+    template <class... Ranges>
+    MultiRange<Ranges...>::MultiRange(const SpaceType& space) : mSpace( space ) {}
+
+    template <class... Ranges>
+    template <size_t N>
+    auto MultiRange<Ranges...>::get() const -> decltype( *std::get<N>( mSpace ) )&
+    {
+	return *std::get<N>(mSpace);
+    }
+
+    template <class... Ranges>
+    template <size_t N>
+    auto MultiRange<Ranges...>::getPtr() const -> decltype( std::get<N>( mSpace ) )&
+    {
+	return std::get<N>(mSpace);
+    }
+
+    template <class... Ranges>
+    size_t MultiRange<Ranges...>::dim() const
+    {
+	return sdim;
+    }
+    
+    template <class... Ranges>
+    size_t MultiRange<Ranges...>::size() const
+    {
+	return RPackNum<sizeof...(Ranges)-1>::getSize(mSpace);
+    }
+
+    template <class... Ranges>
+    const typename MultiRange<Ranges...>::SpaceType& MultiRange<Ranges...>::space() const
+    {
+	return mSpace;
+    }
+    
+    template <class... Ranges>
+    typename MultiRange<Ranges...>::IndexType MultiRange<Ranges...>::begin() const
+    {
+	MultiIndex<typename Ranges::IndexType...>
+	    i( std::dynamic_pointer_cast<MultiRange<Ranges...> >
+	       ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
+	i = 0;
+	return i;
+    }
+
+    template <class... Ranges>
+    typename MultiRange<Ranges...>::IndexType MultiRange<Ranges...>::end() const
+    {
+	MultiIndex<typename Ranges::IndexType...>
+	    i( std::dynamic_pointer_cast<MultiRange<Ranges...> >
+	       ( std::shared_ptr<RangeBase>( RB::mThis )) );
+	i = size();
+	return i;
+    }
+
+    template <class... Ranges>
+    std::shared_ptr<VIWB> MultiRange<Ranges...>::index() const
+    {
+	typedef IndexWrapper<IndexType> IW;
+	return std::make_shared<IW>
+	    ( std::make_shared<IndexType>
+	      ( std::dynamic_pointer_cast<MultiRange<Ranges...> >
+		( std::shared_ptr<RangeBase>( RB::mThis ) ) ) );
+    }
+}
+
+#endif

src/include/ranges/range_base.h

@@ -0,0 +1,93 @@
+// -*- C++ -*-
+
+#ifndef __range_base_h__
+#define __range_base_h__
+
+#include <cstdlib>
+#include <vector>
+#include <memory>
+
+#include "base_def.h"
+
+namespace MultiArrayTools
+{
+    class RangeBase;
+}
+#include "vindex_base.h"
+
+namespace MultiArrayTools
+{
+
+    size_t indexId();
+    
+    enum class SpaceType
+    {
+	NONE = 0,
+	ANY = 1,
+#define include_range_type(x,n) x = n,
+#include "range_types/header.h"
+#undef include_range_type
+    };
+
+    class RangeFactoryBase
+    {
+    public:
+
+	RangeFactoryBase() = default;
+	virtual ~RangeFactoryBase() = default;
+
+	// should return mProd !!
+	virtual std::shared_ptr<RangeBase> create() = 0;
+	
+    protected:
+
+	std::shared_ptr<RangeBase> mProd;
+	
+	// call this function before returning product !!
+	void setSelf();
+    };
+    
+    class RangeBase
+    {
+    public:
+
+	virtual ~RangeBase() = default;
+	    
+	virtual size_t size() const = 0;
+	virtual size_t dim() const = 0;
+	
+	virtual std::shared_ptr<VIWB> index() const = 0;
+
+	bool operator==(const RangeBase& in) const;
+	bool operator!=(const RangeBase& in) const;
+
+	//virtual bool regular() const = 0; // integer distance (e.g. 2,3,4,...)
+	//virtual bool linear() const = 0; // 1dim valuable (e.g. 2.45, 3.12, 3.56,...)
+	//virtual bool multi() const = 0; // mdim
+	//virtual bool maplike() const = 0; // meta type is ~ MultiArray<T,...>
+	
+	friend RangeFactoryBase;
+	
+    protected:
+
+	RangeBase() = default;
+	std::weak_ptr<RangeBase> mThis;
+    };
+
+    template <class Index>
+    class RangeInterface : public RangeBase
+    {
+    public:
+
+	typedef Index IndexType;
+
+	virtual Index begin() const = 0;
+	virtual Index end() const = 0;
+
+    protected:
+	RangeInterface() = default;
+    };
+    
+}
+
+#endif

src/include/ranges/range_types/header.h

@@ -0,0 +1,23 @@
+
+#ifdef include_range_type
+#define __incl_this__
+#endif
+#ifdef __single_range_h__
+// singel_range is template which is specialized here
+// therefore it must be defined before...
+#define __incl_this__
+#endif
+
+#ifdef __incl_this__
+
+#define __ranges_header__
+//#ifndef __ranges_header__
+//#define __ranges_header__
+
+#include "spin_range.h"
+
+#undef __ranges_header__
+//#endif
+#endif
+
+#undef __incl_this__

src/include/ranges/range_types/spin_range.h

@@ -0,0 +1,149 @@
+
+
+#ifdef include_range_type
+include_range_type(SPIN,2)
+#else
+
+#ifdef __ranges_header__
+// assert, that this is only used within range_types/header.h
+
+//#ifndef __spin_range_h__
+//#define __spin_range_h__
+
+namespace MultiArrayTools
+{
+    typedef SingleIndex<size_t,SpaceType::SPIN> SpinIndex;
+
+    template <>
+    class SingleRangeFactory<size_t,SpaceType::SPIN> : public RangeFactoryBase
+    {
+    public:
+	
+	typedef SingleRange<size_t,SpaceType::SPIN> oType;
+
+	SingleRangeFactory();
+	std::shared_ptr<RangeBase> create();
+	
+    };
+
+    template <>
+    class SingleRange<size_t,SpaceType::SPIN> : public RangeInterface<SpinIndex>
+    {
+    public:
+	typedef RangeBase RB;
+	typedef typename RangeInterface<SingleIndex<size_t,SpaceType::SPIN> >::IndexType IndexType;
+	
+	virtual size_t size() const override;
+	virtual size_t dim() const override;
+	
+	size_t get(size_t pos) const;
+	size_t getMeta(size_t metaPos) const;
+	
+	virtual IndexType begin() const override;
+	virtual IndexType end() const override;
+	virtual std::shared_ptr<VIWB> index() const override;
+	
+	friend SingleRangeFactory<size_t,SpaceType::SPIN>;
+
+	static constexpr bool defaultable = true;
+	static constexpr size_t mSpinNum = 4;
+
+	static constexpr size_t ISSTATIC = 1;
+	static constexpr size_t SIZE = mSpinNum;
+	
+	static SingleRangeFactory<size_t, SpaceType::SPIN> factory()
+	{ return SingleRangeFactory<size_t, SpaceType::SPIN>(); }
+	
+    protected:
+
+	SingleRange() = default;
+	SingleRange(const SingleRange& in) = delete;
+
+	//SingleRange(size_t spinNum);
+    };
+
+    typedef SingleRange<size_t,SpaceType::SPIN> SpinRange;
+    typedef SingleRangeFactory<size_t,SpaceType::SPIN> SpinRF;
+}
+
+/* ========================= *
+ * ---   TEMPLATE CODE   --- *
+ * ========================= */
+
+namespace MultiArrayTools
+{
+    /********************
+     *   SingleRange    *
+     ********************/
+    
+    SingleRangeFactory<size_t,SpaceType::SPIN>::SingleRangeFactory()
+    {
+	// Quasi Singleton
+	if(not mProd){
+	    mProd = std::shared_ptr<oType>( new SingleRange<size_t,SpaceType::SPIN>() );
+	    setSelf();
+	}
+    }
+
+    std::shared_ptr<RangeBase> SingleRangeFactory<size_t,SpaceType::SPIN>::create()
+    {
+	return mProd;
+    }
+    
+    /********************
+     *   SingleRange    *
+     ********************/
+    
+    size_t SingleRange<size_t,SpaceType::SPIN>::get(size_t pos) const
+    {
+	return pos;
+    }
+    
+    size_t SingleRange<size_t,SpaceType::SPIN>::getMeta(size_t metaPos) const
+    {
+	return metaPos;
+    }
+    
+    size_t SingleRange<size_t,SpaceType::SPIN>::size() const
+    {
+	return mSpinNum;
+    }
+    
+    size_t SingleRange<size_t,SpaceType::SPIN>::dim() const
+    {
+	return 1;
+    }
+    
+    typename SingleRange<size_t,SpaceType::SPIN>::IndexType SingleRange<size_t,SpaceType::SPIN>::begin() const
+    {
+	SingleIndex<size_t,SpaceType::SPIN> i( std::dynamic_pointer_cast<SingleRange<size_t,SpaceType::SPIN> >
+			       ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
+	i = 0;
+	return i;
+    }
+    
+    typename SingleRange<size_t,SpaceType::SPIN>::IndexType SingleRange<size_t,SpaceType::SPIN>::end() const
+    {
+	SingleIndex<size_t,SpaceType::SPIN> i( std::dynamic_pointer_cast<SingleRange<size_t,SpaceType::SPIN> >
+			       ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
+	i = size();
+	return i;
+    }
+
+    // put this in the interface class !!!
+    std::shared_ptr<VIWB> SingleRange<size_t,SpaceType::SPIN>::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::shared_ptr<RangeBase>( RB::mThis ) ) ) );
+    }
+
+}
+
+//#endif // #ifndef __spin_range_h__
+
+#endif // #ifdef __ranges_header__
+
+#endif // #ifdef include_range_type

src/include/ranges/rbase_def.h

@@ -0,0 +1,80 @@
+// -*- C++ -*-
+
+#ifndef __ranges_base_def_h__
+#define __ranges_base_def_h__
+
+namespace MultiArrayTools
+{
+    /***********************
+     *  Provided Types     *
+     ***********************/
+    
+    // range_base.h
+    enum class SpaceType;
+
+    // range_base.h
+    class RangeFactoryBase;
+    
+    // range_base.h
+    class RangeBase;
+    
+    // range_base.h
+    template <class Index>
+    class RangeInterface;
+
+    // index_base.h
+    class VirtualIndexWrapperBase;
+    typedef VirtualIndexWrapperBase VIWB;    
+    
+    // index_base.h
+    template <class I>
+    class IndexWrapper;
+    
+    // index_base.h
+    template <class I, typename MetaType>
+    class IndexInterface;
+
+    // single_range.h    
+    template <typename U, SpaceType TYPE>
+    class SingleRange;
+
+    // single_range.h    
+    template <typename U, SpaceType TYPE>
+    class SingleRangeFactory;
+
+    // single_range.h
+    template <typename U, SpaceType TYPE>
+    class SingleIndex;
+    
+    // multi_range.h
+    template <class... Ranges>
+    class MultiRangeFactory;
+    
+    // multi_range.h    
+    template <class... Ranges>
+    class MultiRange;
+
+    // multi_range.h
+    template <class... Indices>
+    class MultiIndex;
+
+    // container_range.h
+    template <class... Ranges>
+    class ContainerRangeFactory;
+
+    // container_range.h
+    template <class... Ranges>
+    class ContainerRange;
+
+    // container_range.h
+    template <class... Indices>
+    class ContainerIndex;
+
+    // anonymous_range.h
+    class AnonymousRangeFactory;
+
+    // anonymous_range.h
+    class AnonymousRange;
+}
+
+#endif

src/include/ranges/rheader.h

@@ -0,0 +1,5 @@
+
+#include "single_range.h"
+#include "multi_range.h"
+#include "container_range.h"
+//#include "anonymous_range.h" have to fix bug !!

src/include/ranges/rpack_num.h

@@ -0,0 +1,452 @@
+
+#ifndef __rpack_num_h__
+#define __rpack_num_h__
+
+#include <memory>
+#include "vindex_wrapper.h"
+#include "index_info.h"
+
+namespace MultiArrayHelper
+{
+    using namespace MultiArrayTools;
+
+    
+    template <class Index1>
+    size_t mkTotalDim()
+    {
+	return Index1::totalDim();
+    }
+    
+    template <class Index1, class Index2, class... Indices>
+    size_t mkTotalDim()
+    {
+	return Index1::totalDim() * mkTotalDim<Index2,Indices...>();
+    }
+
+
+    template <class RangeType, class... Ranges>
+    struct SubProp
+    {
+	static constexpr size_t ISSTATIC = RangeType::ISSTATIC & SubProp<Ranges...>::ISSTATIC;
+	static constexpr size_t SIZE = RangeType::SIZE * SubProp<Ranges...>::SIZE;
+    };
+
+    template <>
+    struct SubProp<void>
+    {
+	static constexpr size_t ISSTATIC = 1;
+	static constexpr size_t SIZE = 1;
+    };
+    
+    template <size_t N>
+    struct RPackNum
+    {
+	template <class IndexType>
+	static std::shared_ptr<VIWB> getIndex(const IndexType& in, size_t n)
+	{
+	    if(n == N){
+		return make_viwb( in.template get<N>() );
+	    }
+	    else {
+		return RPackNum<N-1>::getIndex(in, n);
+	    }
+	}
+
+	template <class IndexType>
+	static std::shared_ptr<VIWB> getIndexPtr(const IndexType& in, size_t n)
+	{
+	    if(n == N){
+		return make_viwb( in.template getPtr<N>() );
+	    }
+	    else {
+		return RPackNum<N-1>::getIndexPtr(in, n);
+	    }
+	}
+
+	template <class... Indices>
+	static void initBlockSizes(std::array<size_t,sizeof...(Indices)+1>& bs,
+				   std::tuple<std::shared_ptr<Indices>...>& ip)
+	{
+	    std::get<N>(bs) = RPackNum<sizeof...(Indices)-N-1>::blockSize(ip);
+	    RPackNum<N-1>::initBlockSizes(bs, ip);
+	}
+	
+	template <class... Indices>
+	static inline void pp(std::tuple<std::shared_ptr<Indices>...>& ip)
+	{
+	    auto& si = *std::get<N>(ip);
+	    if(si.last()){
+		si = 0;
+		RPackNum<N-1>::pp(ip);
+	    }
+	    else {
+		++si;
+	    }
+	}
+	
+	template <class... Indices>
+	static inline int pp(std::tuple<std::shared_ptr<Indices>...>& ip,
+			     std::array<size_t,sizeof...(Indices)+1>& bs,
+			     std::intptr_t idxPtrNum)
+	{
+	    auto& siPtr = std::get<N>(ip);
+	    //VCHECK(siPtr.id());
+	    if(reinterpret_cast<std::intptr_t>(siPtr.get()) == idxPtrNum){
+		return RPackNum<N-1>::pp(ip, bs, idxPtrNum);
+	    }
+	    else {
+		int tmp = siPtr->pp(idxPtrNum);
+		if(siPtr->pos() == siPtr->max()){
+		    (*siPtr) = 0;
+		    return RPackNum<N-1>::pp(ip, bs, idxPtrNum) - siPtr->max() + 1;
+		}
+		else {
+		    return tmp * std::get<N+1>(bs);
+		}
+	    }
+	}
+
+	template <class... Indices>
+	static inline void mm(std::tuple<std::shared_ptr<Indices>...>& ip)
+	{
+	    auto& si = *std::get<N>(ip);
+	    if(si.first()){
+		si = si.max() - 1;
+		RPackNum<N-1>::mm(ip);
+	    }
+	    else {
+		--si;
+	    }
+	}
+	
+	// !!!!
+	template <class... Indices>
+	static inline int mm(std::tuple<std::shared_ptr<Indices>...>& ip,
+			     std::array<size_t,sizeof...(Indices)+1>& bs,
+			     std::intptr_t idxPtrNum)
+	{
+	    auto& siPtr = std::get<N>(ip);
+	    if(reinterpret_cast<std::intptr_t>(siPtr.get()) == idxPtrNum){
+		return std::get<N>(bs) + RPackNum<N-1>::mm(ip, bs, idxPtrNum);
+	    }
+	    else {
+		if(siPtr->first()){
+		    (*siPtr) = siPtr->max() - 1;
+		    return RPackNum<N-1>::mm(ip, bs, idxPtrNum) - siPtr->max() + 1;
+		}
+		else {
+		    return siPtr->mm(idxPtrNum);
+		}
+	    }
+	}
+	
+	template <class RangeTuple>
+	static size_t getSize(const RangeTuple& rt)
+	{
+	    return std::get<N>(rt)->size() * RPackNum<N-1>::getSize(rt);
+	}
+	
+	template <class IndexPack, class MetaType>
+	static void getMetaPos(MetaType& target,
+			       const IndexPack& source)
+	{
+	    std::get<N>(target) = std::get<N>(source)->meta();
+	    RPackNum<N-1>::getMetaPos(target, source);
+	}
+	
+	template <class IndexPack, typename MetaType>
+	static void setMeta(IndexPack& target, const MetaType& source)
+	{
+	    std::get<N>(target)->at( std::get<N>(source) );
+	    RPackNum<N-1>::setMeta(target, source);
+	}
+	
+	template <class IndexPack>
+	static void setIndexPack(IndexPack& iPack, size_t pos)
+	{
+	    auto& i = *std::get<N>(iPack).get();
+	    const size_t ownPos = pos % i.max(); 
+	    i = ownPos;
+	    RPackNum<N-1>::setIndexPack(iPack, (pos - ownPos) / i.max() );
+	}
+	
+	template <class MRange, class... Indices>
+	static void construct(std::tuple<std::shared_ptr<Indices>...>& ip,
+			      const MRange& range)
+	{
+	    typedef typename std::remove_reference<decltype(range.template get<N>())>::type SubRangeType;
+	    typedef typename SubRangeType::IndexType SubIndexType;
+	    typedef typename std::remove_reference<decltype(*std::get<N>(ip).get())>::type TypeFromIndexPack;
+	    
+	    static_assert(std::is_same<SubIndexType,TypeFromIndexPack>::value,
+			  "inconsiśtent types");
+	    
+	    std::get<N>(ip) = std::shared_ptr<SubIndexType>( new SubIndexType( range.template getPtr<N>() ) );
+	    RPackNum<N-1>::construct(ip, range);
+	}
+
+	template <template<class...> class IndexType, class... Indices>
+	static void copyInst(std::tuple<std::shared_ptr<Indices>...>& ip,
+			     const IndexType<Indices...>& ind)
+	{
+	    std::get<N>(ip) = ind.template getPtr<N>() ;
+	    RPackNum<N-1>::copyInst(ip, ind);
+	}
+	
+	template <class... Indices>
+	static inline size_t makePos(const std::tuple<std::shared_ptr<Indices>...>& iPtrTup)
+	{
+	    //const auto& idx = *std::get<N>(iPtrTup);
+	    return std::get<N>(iPtrTup)->pos() + RPackNum<N-1>::makePos(iPtrTup) * std::get<N>(iPtrTup)->max();
+	}
+
+	template <class Pack, class IndexType, class... Indices>
+	static void swapIndices(Pack& ipack, const std::shared_ptr<IndexType>& nind,
+				const std::shared_ptr<Indices>&... ninds)
+	{
+	    std::get<std::tuple_size<Pack>::value-N-1>(ipack) = nind;
+	    RPackNum<N-1>::swapIndices(ipack, ninds...);
+	}
+
+	template <class... Indices>
+	static size_t blockSize(const std::tuple<std::shared_ptr<Indices>...>& pack)
+	{
+	    return std::get<sizeof...(Indices)-N-1>(pack)->max() * RPackNum<N-1>::blockSize(pack);
+	}
+
+
+	template <class... Ranges>
+	static void RangesToVec(const std::tuple<std::shared_ptr<Ranges>...>& rst,
+				std::vector<RangeBase> v)
+	{
+	    v[N] = std::get<N>(rst);
+	    RPackNum<N-1>::RangesToVec(rst, v);
+	}
+
+	template <class... Indices>
+	static void printIndex(const std::tuple<std::shared_ptr<Indices>...>& ip, size_t offset)
+	{
+	    std::get<N>(ip)->print(offset);
+	    RPackNum<N-1>::printIndex(ip, offset);
+	}
+
+	template <class Range, class... Ranges>
+	static void checkDefaultable()
+	{
+	    static_assert( Range::defaultable, "not defaultable" );
+	    RPackNum<N-1>::template checkDefaultable<Ranges...>();
+	}
+
+	template <class... Indices>
+	static void buildInfoVec(std::vector<IndexInfo>& out,
+				 const std::tuple<std::shared_ptr<Indices>...>& ip,
+				 const std::array<size_t,sizeof...(Indices)+1>& bs)
+	{
+	    static const size_t POS = sizeof...(Indices)-N-1;
+	    out.emplace_back(*std::get<POS>(ip), std::get<POS>(bs));
+	    RPackNum<N-1>::buildInfoVec(out, ip, bs);
+	}
+
+	template <class IndexPack, class Exprs>
+	static auto mkFor(const IndexPack& ipack, Exprs exs)
+	    -> decltype(std::get<std::tuple_size<IndexPack>::value-N-1>(ipack)
+			->ifor( RPackNum<N-1>::mkFor(ipack, exs) ) )
+	{
+	    return std::get<std::tuple_size<IndexPack>::value-N-1>(ipack)
+		->ifor( RPackNum<N-1>::mkFor(ipack, exs) );
+	}
+
+	template <class IndexPack, class Exprs>
+	static auto mkForh(const IndexPack& ipack, Exprs exs)
+	    -> decltype(std::get<std::tuple_size<IndexPack>::value-N-1>(ipack)
+			->iforh( RPackNum<N-1>::mkForh(ipack, exs) ) )
+	{
+	    return std::get<std::tuple_size<IndexPack>::value-N-1>(ipack)
+		->iforh( RPackNum<N-1>::mkForh(ipack, exs) );
+	}
+
+    };
+
+    
+    template<>
+    struct RPackNum<0>
+    {
+	template <class IndexType>
+	static std::shared_ptr<VIWB> getIndex(const IndexType& in, size_t n)
+	{
+	    return make_viwb( in.template get<0>() );
+	}
+
+	template <class IndexType>
+	static std::shared_ptr<VIWB> getIndexPtr(const IndexType& in, size_t n)
+	{
+	    return make_viwb( in.template getPtr<0>() );
+	}
+
+	template <class... Indices>
+	static void initBlockSizes(std::array<size_t,sizeof...(Indices)+1>& bs,
+				   std::tuple<std::shared_ptr<Indices>...>& ip)
+	{
+	    std::get<0>(bs) = RPackNum<sizeof...(Indices)-1>::blockSize(ip);
+	}
+
+	template <class... Indices>
+	static inline void pp(std::tuple<std::shared_ptr<Indices>...>& ip)
+	{
+	    auto& si = *std::get<0>(ip);
+	    ++si;
+	}
+	
+	template <class... Indices>
+	static inline int pp(std::tuple<std::shared_ptr<Indices>...>& ip,
+			     std::array<size_t,sizeof...(Indices)+1>& bs,
+			     std::intptr_t idxPtrNum)
+	{
+	    auto& siPtr = std::get<0>(ip);
+	    if(reinterpret_cast<std::intptr_t>(siPtr.get()) == idxPtrNum){
+		return std::get<0>(bs);
+	    }
+	    else {
+		int tmp = siPtr->pp(idxPtrNum);
+		return tmp * std::get<1>(bs);
+	    }
+	}
+
+	template <class... Indices>
+	static inline void mm(std::tuple<std::shared_ptr<Indices>...>& ip)
+	{
+	    auto& si = *std::get<0>(ip);
+	    --si;
+	}
+	
+	template <class... Indices>
+	static inline int mm(std::tuple<std::shared_ptr<Indices>...>& ip,
+			     std::array<size_t,sizeof...(Indices)+1>& bs,
+			     std::intptr_t idxPtrNum)
+	{
+	    auto& siPtr = std::get<0>(ip);
+	    if(reinterpret_cast<std::intptr_t>(siPtr.get()) == idxPtrNum){
+		return std::get<0>(bs);
+		//return 1;
+	    }
+	    else {
+		return siPtr->mm(idxPtrNum);
+	    }
+	}
+	
+	template <class RangeTuple>
+	static size_t getSize(const RangeTuple& rt)
+	{
+	    return std::get<0>(rt)->size();
+	}
+	
+	template <class IndexPack, class MetaType>
+	static void getMetaPos(MetaType& target,
+			       const IndexPack& source)
+	{
+	    std::get<0>(target) = std::get<0>(source)->meta();
+	}
+	
+	template <class IndexPack, typename MetaType>
+	static void setMeta(IndexPack& target, const MetaType& source)
+	{
+	    std::get<0>(target)->at( std::get<0>( source ) );
+	}
+	
+	template <class IndexPack>
+	static void setIndexPack(IndexPack& iPack, size_t pos)
+	{
+	    auto& i = *std::get<0>(iPack);
+	    const size_t ownPos = pos % i.max(); 
+	    i = ownPos;
+	}
+	
+	template <class MRange, class... Indices>
+	static void construct(std::tuple<std::shared_ptr<Indices>...>& ip,
+			      const MRange& range)
+	{
+	    typedef typename std::remove_reference<decltype(range.template get<0>())>::type SubRangeType;
+	    typedef typename SubRangeType::IndexType SubIndexType;
+	    typedef typename std::remove_reference<decltype(*std::get<0>(ip).get())>::type TypeFromIndexPack;
+	    
+	    static_assert(std::is_same<SubIndexType,TypeFromIndexPack>::value,
+			  "inconsiśtent types");
+	    
+	    std::get<0>(ip) = std::shared_ptr<SubIndexType>( new SubIndexType( range.template getPtr<0>() ) );
+	}
+
+	template <template<class...> class IndexType, class... Indices>
+	static void copyInst(std::tuple<std::shared_ptr<Indices>...>& ip,
+			     const IndexType<Indices...>& ind)
+	{
+	    std::get<0>(ip) = ind.template getPtr<0>();
+	}
+	
+	template <class... Indices>
+	static inline size_t makePos(const std::tuple<std::shared_ptr<Indices>...>& iPtrTup)
+	{
+	    return std::get<0>(iPtrTup)->pos();
+	}
+	
+	template <class Pack, class IndexType>
+	static void swapIndices(Pack& ipack, const std::shared_ptr<IndexType>& nind)
+	{
+	    std::get<std::tuple_size<Pack>::value-1>(ipack) = nind;
+	}
+
+	template <class... Indices>
+	static size_t blockSize(const std::tuple<std::shared_ptr<Indices>...>& pack)
+	{
+	    return std::get<sizeof...(Indices)-1>(pack)->max();
+	}
+
+	template <class... Ranges>
+	static void RangesToVec(const std::tuple<std::shared_ptr<Ranges>...>& rst,
+				std::vector<RangeBase> v)
+	{
+	    v[0] = std::get<0>(rst);
+	}
+
+	template <class... Indices>
+	static void printIndex(const std::tuple<std::shared_ptr<Indices>...>& ip, size_t offset)
+	{
+	    std::get<0>(ip)->print(offset);
+	}
+
+	template <class Range>
+	static void checkDefaultable()
+	{
+	    static_assert( Range::defaultable, "not defaultable" );
+	}
+
+	template <class... Indices>
+	static void buildInfoVec(std::vector<IndexInfo>& out,
+				 const std::tuple<std::shared_ptr<Indices>...>& ip,
+				 const std::array<size_t,sizeof...(Indices)+1>& bs)
+	{
+	    static const size_t POS = sizeof...(Indices)-1;
+	    out.emplace_back(*std::get<POS>(ip), std::get<POS>(bs));
+	}
+
+	template <class IndexPack, class Exprs>
+	static auto mkFor(const IndexPack& ipack, Exprs exs)
+	    -> decltype(std::get<std::tuple_size<IndexPack>::value-1>(ipack)
+			->ifor(exs) )
+	{
+	    return std::get<std::tuple_size<IndexPack>::value-1>(ipack)->ifor(exs);
+	}
+
+	template <class IndexPack, class Exprs>
+	static auto mkForh(const IndexPack& ipack, Exprs exs)
+	    -> decltype(std::get<std::tuple_size<IndexPack>::value-1>(ipack)
+			->iforh(exs) )
+	{
+	    return std::get<std::tuple_size<IndexPack>::value-1>(ipack)->iforh(exs);
+	}
+
+    };
+	
+
+} //  end namespace MultiArrayHelper
+
+#endif

src/include/ranges/single_range.h

@@ -0,0 +1,367 @@
+// -*- C++ -*-
+
+#ifndef __single_range_h__
+#define __single_range_h__
+
+#include <cstdlib>
+#include <vector>
+#include <memory>
+
+//#include "base_def.h"
+#include "ranges/index_base.h"
+#include "ranges/range_base.h"
+
+#include "xfor/xfor.h"
+
+using MultiArrayHelper::For;
+
+namespace MultiArrayTools
+{
+
+    template <typename U, SpaceType TYPE>
+    class SingleIndex : public IndexInterface<SingleIndex<U,TYPE>,U>
+    {
+    public:
+
+	typedef IndexInterface<SingleIndex<U,TYPE>,U> IB;
+	typedef U MetaType;
+	typedef SingleRange<U,TYPE> RangeType;
+
+	//DEFAULT_MEMBERS_X(SingleIndex);
+	
+	SingleIndex(const std::shared_ptr<SingleRange<U,TYPE> >& range);
+
+	static IndexType sType() { return IndexType::SINGLE; }
+	static size_t totalDim() { return 1; }
+
+	// ==== >>>>> STATIC POLYMORPHISM <<<<< ====
+	
+	IndexType type() const;
+	
+	SingleIndex& operator=(size_t pos);
+	SingleIndex& operator++();
+	SingleIndex& operator--();
+
+	int pp(std::intptr_t idxPtrNum);
+	int mm(std::intptr_t idxPtrNum);
+	
+	U meta();
+	SingleIndex& at(const U& metaPos);
+	
+	size_t dim(); // = 1
+	bool last();
+	bool first();
+
+	std::shared_ptr<RangeType> range();
+	
+	template <size_t N>
+	void getPtr();
+
+	std::shared_ptr<VIWB> getVPtr(size_t n);
+	size_t getStepSize(size_t n);
+
+	std::vector<IndexInfo> infoVec() const;
+	
+	std::string id() const;
+	void print(size_t offset);
+
+	template <class Expr>
+	auto ifor(Expr ex) const
+	    -> For<SingleIndex<U,TYPE>,Expr>;
+
+	template <class Expr>
+	auto iforh(Expr ex) const
+	    -> For<SingleIndex<U,TYPE>,Expr,ForType::HIDDEN>;
+
+    };
+
+    template <typename U, SpaceType TYPE>
+    class SingleRangeFactory : public RangeFactoryBase
+    {
+    public:
+	
+	typedef SingleRange<U,TYPE> oType;
+
+	SingleRangeFactory() = delete;
+	SingleRangeFactory(const std::vector<U>& space);
+	std::shared_ptr<RangeBase> create();
+	
+    };
+    
+    template <typename U, SpaceType TYPE>
+    class SingleRange : public RangeInterface<SingleIndex<U,TYPE> >
+    {
+    public:
+	typedef RangeBase RB;
+	typedef SingleIndex<U,TYPE> IndexType;
+	//typedef typename RangeInterface<SingleIndex<U,TYPE> >::IndexType IndexType;
+	
+	virtual size_t size() const override;
+	virtual size_t dim() const override;
+	
+	const U& get(size_t pos) const;
+	size_t getMeta(const U& metaPos) const;
+	
+	virtual IndexType begin() const override;
+	virtual IndexType end() const override;
+	virtual std::shared_ptr<VIWB> index() const override;
+	
+	friend SingleRangeFactory<U,TYPE>;
+
+	static constexpr bool defaultable = false;
+	static constexpr size_t ISSTATIC = 0;
+	static constexpr size_t SIZE = -1;
+
+	
+    protected:
+
+	SingleRange() = delete;
+	SingleRange(const SingleRange& in) = delete;
+	
+	SingleRange(const std::vector<U>& space);
+
+	std::vector<U> mSpace;
+    };
+
+}
+
+/* ========================= *
+ * ---   TEMPLATE CODE   --- *
+ * ========================= */
+
+namespace MultiArrayTools
+{
+    /******************
+     *  SingleIndex   *	     
+     ******************/
+
+    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>
+    IndexType SingleIndex<U,TYPE>::type() const
+    {
+	return IndexType::SINGLE;
+    }
+
+    template <typename U, SpaceType TYPE>
+    SingleIndex<U,TYPE>& SingleIndex<U,TYPE>::operator=(size_t pos)
+    {
+	IB::mPos = pos;
+	return *this;
+    }
+
+    template <typename U, SpaceType TYPE>
+    SingleIndex<U,TYPE>& SingleIndex<U,TYPE>::operator++()
+    {
+	++IB::mPos;
+	return *this;
+    }
+
+    template <typename U, SpaceType TYPE>
+    SingleIndex<U,TYPE>& SingleIndex<U,TYPE>::operator--()
+    {
+	--IB::mPos;
+	return *this;
+    }
+
+    template <typename U, SpaceType TYPE>
+    int SingleIndex<U,TYPE>::pp(std::intptr_t idxPtrNum)
+    {
+	++(*this);
+	return 1;
+    }
+
+    template <typename U, SpaceType TYPE>
+    int SingleIndex<U,TYPE>::mm(std::intptr_t idxPtrNum)
+    {
+	--(*this);
+	return 1;
+    }
+
+    template <typename U, SpaceType TYPE>
+    U SingleIndex<U,TYPE>::meta()
+    {
+	return std::dynamic_pointer_cast<SingleRange<U,TYPE> const>( IB::mRangePtr )->get( IB::pos() );
+    }
+
+    template <typename U, SpaceType TYPE>
+    SingleIndex<U,TYPE>& SingleIndex<U,TYPE>::at(const U& metaPos)
+    {
+	(*this) = std::dynamic_pointer_cast<SingleRange<U,TYPE> const>( IB::mRangePtr )->getMeta( metaPos );
+	return *this;
+    }
+
+    template <typename U, SpaceType TYPE>
+    size_t SingleIndex<U,TYPE>::dim() // = 1
+    {
+	return 1;
+    }
+
+    template <typename U, SpaceType TYPE>
+    bool SingleIndex<U,TYPE>::last()
+    {
+	return IB::mPos == IB::mMax - 1;
+    }
+
+    template <typename U, SpaceType TYPE>
+    bool SingleIndex<U,TYPE>::first()
+    {
+	return IB::mPos == 0;
+    }
+
+    template <typename U, SpaceType TYPE>
+    std::shared_ptr<typename SingleIndex<U,TYPE>::RangeType> SingleIndex<U,TYPE>::range()
+    {
+	return std::dynamic_pointer_cast<RangeType>( IB::mRangePtr );
+    }
+
+    template <typename U, SpaceType TYPE>
+    template <size_t N>
+    void SingleIndex<U,TYPE>::getPtr() {}
+
+    template <typename U, SpaceType TYPE>
+    std::shared_ptr<VIWB> SingleIndex<U,TYPE>::getVPtr(size_t n)
+    {
+	return std::shared_ptr<VIWB>();
+    }
+
+    template <typename U, SpaceType TYPE>
+    std::vector<IndexInfo> SingleIndex<U,TYPE>::infoVec() const
+    {
+	return std::move( std::vector<IndexInfo>() );
+    }
+    
+    template <typename U, SpaceType TYPE>
+    size_t SingleIndex<U,TYPE>::getStepSize(size_t n)
+    {
+	return 1;
+    }
+
+    template <typename U, SpaceType TYPE>
+    std::string SingleIndex<U,TYPE>::id() const
+    {
+	return std::string("sin") + std::to_string(IB::mId);
+    }
+
+    template <typename U, SpaceType TYPE>
+    void SingleIndex<U,TYPE>::print(size_t offset)
+    {
+	if(offset == 0){
+	    std::cout << " === " << std::endl;
+	}
+	for(size_t j = 0; j != offset; ++j) { std::cout << "\t"; }
+	std::cout << id() << "[" << reinterpret_cast<std::intptr_t>(this)
+		  << "](" << IB::mRangePtr << "): " << meta() << std::endl;
+    }
+
+    template <typename U, SpaceType TYPE>
+    template <class Expr>
+    auto SingleIndex<U,TYPE>::ifor(Expr ex) const
+	-> For<SingleIndex<U,TYPE>,Expr>
+    {
+	//static const size_t LAYER = typename Expr::LAYER; 
+	return For<SingleIndex<U,TYPE>,Expr>(this, ex);
+    }
+
+    template <typename U, SpaceType TYPE>
+    template <class Expr>
+    auto SingleIndex<U,TYPE>::iforh(Expr ex) const
+	-> For<SingleIndex<U,TYPE>,Expr,ForType::HIDDEN>
+    {
+	//static const size_t LAYER = typename Expr::LAYER; 
+	return For<SingleIndex<U,TYPE>,Expr,ForType::HIDDEN>(this, ex);
+    }
+
+    
+    /********************
+     *   SingleRange    *
+     ********************/
+
+    template <typename U, SpaceType TYPE>
+    SingleRangeFactory<U,TYPE>::SingleRangeFactory(const std::vector<U>& space)
+    {
+	mProd = std::shared_ptr<oType>( new SingleRange<U,TYPE>( space ) );
+    }
+
+    template <typename U, SpaceType TYPE>
+    std::shared_ptr<RangeBase> SingleRangeFactory<U,TYPE>::create()
+    {
+	setSelf();
+	return mProd;
+    }
+    
+    /********************
+     *   SingleRange    *
+     ********************/
+    
+    template <typename U, SpaceType TYPE>
+    SingleRange<U,TYPE>::SingleRange(const std::vector<U>& space) : RangeInterface<SingleIndex<U,TYPE> >(),
+	mSpace(space) {}
+    
+    template <typename U, SpaceType TYPE>
+    const U& SingleRange<U,TYPE>::get(size_t pos) const
+    {
+	return mSpace[pos];
+    }
+
+    template <typename U, SpaceType TYPE>
+    size_t SingleRange<U,TYPE>::getMeta(const U& metaPos) const
+    {
+	size_t cnt = 0;
+	for(auto& x: mSpace){
+	    if(x == metaPos){
+		return cnt;
+	    }
+	    ++cnt;
+	}
+	return cnt;
+    }
+
+    template <typename U, SpaceType TYPE>
+    size_t SingleRange<U,TYPE>::size() const
+    {
+	return mSpace.size();
+    }
+
+    template <typename U, SpaceType TYPE>
+    size_t SingleRange<U,TYPE>::dim() const
+    {
+	return 1;
+    }
+
+    template <typename U, SpaceType TYPE>
+    typename SingleRange<U,TYPE>::IndexType SingleRange<U,TYPE>::begin() const
+    {
+	SingleIndex<U,TYPE> i( std::dynamic_pointer_cast<SingleRange<U,TYPE> >
+			       ( 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
+    {
+	SingleIndex<U,TYPE> i( std::dynamic_pointer_cast<SingleRange<U,TYPE> >
+			       ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
+	i = size();
+	return i;
+    }
+
+    // put this in the interface class !!!
+    template <typename U, SpaceType TYPE>
+    std::shared_ptr<VIWB> SingleRange<U,TYPE>::index() const
+    {
+	typedef IndexWrapper<IndexType> IW;
+	return std::make_shared<IW>
+	    ( std::make_shared<IndexType>
+	      ( std::dynamic_pointer_cast<SingleRange<U,TYPE> >
+		( std::shared_ptr<RangeBase>( RB::mThis ) ) ) );
+    }
+
+}
+
+#include "range_types/header.h"
+
+#endif

src/include/ranges/vindex_base.h

@@ -0,0 +1,35 @@
+
+#ifndef __vindex_base_h__
+#define __vindex_base_h__
+
+//#include "ranges/range_base.h"
+#include "ranges/index_type.h"
+#include "base_def.h"
+
+namespace MultiArrayTools
+{
+
+#ifndef __range_base_h__
+    class RangeBase;
+#endif
+    
+    class VirtualIndexWrapperBase
+    {
+    public:
+	DEFAULT_MEMBERS(VirtualIndexWrapperBase);
+	
+	virtual IndexType type() const = 0;
+	virtual size_t dim() const = 0;
+	virtual size_t pos() const = 0;
+	virtual size_t max() const = 0;
+	virtual std::shared_ptr<RangeBase> rangePtr() const = 0;
+	virtual std::shared_ptr<VirtualIndexWrapperBase> getPtr(size_t n) const = 0;
+	virtual std::intptr_t getPtrNum() const = 0;
+	virtual size_t getStepSize(size_t n) const = 0;
+    };
+
+    typedef VirtualIndexWrapperBase VIWB;
+
+} // end namespace MultiArrayTools
+
+#endif

src/include/ranges/vindex_wrapper.h

@@ -0,0 +1,43 @@
+
+#ifndef __vindex_wrapper_h__
+#define __vindex_wrapper_h__
+
+#include "ranges/vindex_base.h"
+
+namespace MultiArrayTools
+{
+    
+    template <class I>
+    std::shared_ptr<IndexWrapper<I> > make_viwb(std::shared_ptr<I> idxPtr)
+    {
+	return std::make_shared<IndexWrapper<I> >(idxPtr);
+    }
+    
+    template <class I>
+    class IndexWrapper : public VirtualIndexWrapperBase
+    {
+    public:
+	
+	DEFAULT_MEMBERS(IndexWrapper);
+
+	IndexWrapper(std::shared_ptr<I> idxPtr) : mIdxPtr(idxPtr) {}
+
+	virtual IndexType type() const override { return mIdxPtr->type(); }
+	virtual size_t dim() const override { return mIdxPtr->dim(); }
+	virtual size_t pos() const override { return mIdxPtr->pos(); }
+	virtual size_t max() const override { return mIdxPtr->max(); }
+	virtual std::shared_ptr<RangeBase> rangePtr() const override { return mIdxPtr->vrange(); }
+	virtual std::shared_ptr<VirtualIndexWrapperBase> getPtr(size_t n) const override
+	{ return mIdxPtr->getVPtr(n); }
+	virtual std::intptr_t getPtrNum() const override { return reinterpret_cast<std::intptr_t>( mIdxPtr.get() ); }
+	virtual size_t getStepSize(size_t n) const override { return mIdxPtr->getStepSize(n); }
+
+	std::shared_ptr<I> get() const { return mIdxPtr; } // unwrap
+	
+    private:
+	std::shared_ptr<I> mIdxPtr;
+    };
+
+}
+
+#endif

src/include/utils.h

@@ -0,0 +1,30 @@
+
+#ifndef __utils_h__
+#define __utils_h__
+
+#include <ostream>
+#include "pack_num.h"
+
+namespace MultiArrayHelper
+{
+
+    using namespace MultiArrayTools;
+
+    template <typename... T>
+    std::ostream& operator<<(std::ostream& out, const std::tuple<T...>& tp);
+
+    /*==================*
+     *   TEMPLATE CODE  *
+     *==================*/
+    
+    template <typename... T>
+    std::ostream& operator<<(std::ostream& out, const std::tuple<T...>& tp)
+    {
+	PackNum<sizeof...(T)-1>::printTuple(out, tp);
+	return out;
+    }
+
+
+} // end namespace MultiArrayHelper
+
+#endif

src/include/xfor/exttype.h

@@ -0,0 +1,220 @@
+
+#ifndef __exttype_h__
+#define __exttype_h__
+
+#include <array>
+
+namespace MultiArrayHelper
+{
+
+    template <class X>
+    class MExt
+    {
+    private:
+
+	size_t mExt = 0u;
+	X mNext;
+	
+    public:
+	
+	static constexpr size_t NUM = X::NUM + 1;
+
+	MExt() = default;
+	MExt(const MExt& in) = default;
+	MExt& operator=(const MExt& in) = default;
+	MExt(MExt&& in) = default;
+	MExt& operator=(MExt&& in) = default;
+	
+	inline MExt(size_t ext, X next);
+
+	template <class Z>
+	inline MExt(size_t y, const Z& z);
+	
+	template <class Y, class Z>
+	inline MExt(const Y& y, const Z& z);
+
+	template <size_t N>
+	inline MExt(const std::array<size_t,N>& arr);
+
+	inline size_t val() const;
+	inline const X& next() const;
+	
+	inline MExt operator+(const MExt& in) const;
+	inline MExt operator*(size_t in) const;
+
+	template <class Y>
+	auto extend(const Y& y) const -> MExt<decltype(mNext.extend(y))>
+	{ return MExt<decltype(mNext.extend(y))>(mExt, mNext.extend(y)); }
+	
+    };    
+
+    template <>
+    class MExt<void>
+    {
+    private:
+
+	size_t mExt = 0u;
+	
+    public:
+	
+	static constexpr size_t NUM = 0;
+
+	MExt() = default;
+	MExt(const MExt& in) = default;
+	MExt& operator=(const MExt& in) = default;
+	MExt(MExt&& in) = default;
+	MExt& operator=(MExt&& in) = default;
+	
+	inline MExt(size_t ext);
+
+	template <class Z>
+	inline MExt(size_t y, const Z& z);
+	
+	template <class Y, class Z>
+	inline MExt(const Y& y, const Z& z);
+
+	template <size_t N>
+	inline MExt(const std::array<size_t,N>& arr);
+
+	inline size_t val() const;
+	inline size_t next() const { return 0; }
+	
+	inline MExt operator+(const MExt& in) const;
+	inline MExt operator*(size_t in) const;
+
+	template <class Y>
+	auto extend(const Y& y) const -> MExt<Y>
+	{ return MExt<Y>(mExt,y); }
+
+    };    
+
+
+    template <size_t I>
+    struct Getter
+    {
+	template <class ExtType>
+	static inline size_t get(const ExtType& et)
+	{
+	    return Getter<I-1>::get(et.next());
+	}
+
+	template <class ExtType>
+	static inline auto getX(const ExtType& et)
+	    -> decltype(Getter<I-1>::getX(et.next()))
+	{
+	    return Getter<I-1>::getX(et.next());
+	}	
+    };
+
+    template <>
+    struct Getter<0>
+    {
+	template <class ExtType>
+	static inline size_t get(const ExtType& et)
+	{
+	    return et.get();
+	}
+
+	template <class ExtType>
+	static inline auto getX(const ExtType& et)
+	    -> ExtType
+	{
+	    return et;
+	}	
+    };
+
+} // end namespace MultiArrayHelper
+
+/* ========================= *
+ * ---   TEMPLATE CODE   --- *
+ * ========================= */
+
+namespace MultiArrayHelper
+{
+
+    template <class X>
+    inline MExt<X>::MExt(size_t ext, X next) : mExt(ext), mNext(next) {}
+
+    template <class X>
+    template <size_t N>
+    inline MExt<X>::MExt(const std::array<size_t,N>& arr) :
+	mExt(std::get<NUM>(arr)), mNext(arr) {}
+
+    template <class X>
+    template <class Z>
+    inline MExt<X>::MExt(size_t y, const Z& z) :
+	mExt(z.val()), mNext(z.val(), z.next()) {}
+    
+    template <class X>
+    template <class Y, class Z>
+    inline MExt<X>::MExt(const Y& y, const Z& z) :
+	mExt(y.val()), mNext(y.next(), z) {}
+    
+    template <class X>
+    inline size_t MExt<X>::val() const
+    {
+	return mExt;
+    }
+
+    template <class X>
+    inline const X& MExt<X>::next() const
+    {
+	return mNext;
+    }
+    
+    template <class X>
+    inline MExt<X> MExt<X>::operator+(const MExt<X>& in) const
+    {
+	return MExt<X>(mExt + in.val(), mNext + in.next());
+    }
+
+    template <class X>
+    inline MExt<X> MExt<X>::operator*(size_t in) const
+    {
+	return MExt<X>(mExt * in, mNext * in);
+    }
+
+    
+    //template <>
+    inline MExt<void>::MExt(size_t ext) : mExt(ext) {}
+
+
+    //template <>
+    template <class Z>
+    inline MExt<void>::MExt(size_t y, const Z& z) :
+	mExt(z.val()) {}
+    
+    //template <>
+    template <class Y, class Z>
+    inline MExt<void>::MExt(const Y& y, const Z& z) :
+	mExt(y.val()) {}
+
+
+    //template <>
+    template <size_t N>
+    inline MExt<void>::MExt(const std::array<size_t,N>& arr) :
+	mExt(std::get<NUM>(arr)) {}
+
+    //template <>
+    inline size_t MExt<void>::val() const
+    {
+	return mExt;
+    }
+    
+    //template <>
+    inline MExt<void> MExt<void>::operator+(const MExt<void>& in) const
+    {
+	return MExt<void>(mExt + in.val());
+    }
+
+    //template <>
+    inline MExt<void> MExt<void>::operator*(size_t in) const
+    {
+	return MExt<void>(mExt * in);
+    }
+
+
+} // end namespace MultiArrayHelper
+
+
+#endif

src/include/xfor/for_type.h

@@ -0,0 +1,18 @@
+
+#ifndef __for_type_h__
+#define __for_type_h__
+
+namespace MultiArrayHelper
+{
+
+    
+    enum class ForType {
+	DEFAULT = 0,
+	HIDDEN = 1	
+    };
+
+    
+} // end namespace MultiArrayHelper
+
+
+#endif

src/include/xfor/for_utils.h

@@ -0,0 +1,39 @@
+
+#ifndef __for_utils_h__
+#define __for_utils_h__
+
+#include "ranges/rheader.h"
+#include <array>
+#include <tuple>
+
+namespace MultiArrayHelper
+{
+    namespace {
+        template <class Op>
+        using to_size_t = size_t;
+    }
+    
+    template <size_t N>
+    struct XFPackNum
+    {
+	template <class ETuple, typename... Args>
+	static inline ETuple mkPos(size_t pos, const ETuple& et, const ETuple& lt, const Args&... args)
+	{
+	    return std::move( XFPackNum<N-1>::mkPos(pos, et, lt, std::get<N>(lt) + pos * std::get<N>(et), args...) );
+	}
+    };
+
+    template <>
+    struct XFPackNum<0>
+    {
+	template <class ETuple, typename... Args>
+	static inline ETuple mkPos(size_t pos, const ETuple& et, const ETuple& lt, const Args&... args)
+	{
+	    return ETuple(std::get<0>(lt) + pos * std::get<0>(et), args...);
+	}
+
+    };
+    
+} // end namespace MultiArrayHelper
+
+#endif

src/include/xfor/xfor.h

@@ -0,0 +1,172 @@
+
+#ifndef __xfor_h__
+#define __xfor_h__
+
+#include <cstdlib>
+#include <memory>
+#include <tuple>
+#include "xfor/for_utils.h"
+#include "xfor/for_type.h"
+#include "xfor/exttype.h"
+
+namespace MultiArrayHelper
+{
+
+    // 'HIDDEN FOR' CLASS for nested for loops in contractions a.s.o.
+    // (NO COUNTING OF MASTER POSITION !!!!!)
+
+    template <ForType FT = ForType::DEFAULT>
+    struct PosForward
+    {
+	static inline size_t value(size_t last, size_t max, size_t pos)
+	{
+	    return last * max + pos;
+	}
+    };
+
+    template <>
+    struct PosForward<ForType::HIDDEN>
+    {
+	static inline size_t value(size_t last, size_t max, size_t pos)
+	{
+	    return last;
+	}
+    };
+
+    template <size_t ISSTATIC>
+    struct ForBound
+    {
+	template <size_t BOUND>
+	static inline size_t bound(size_t bound)
+	{
+	    return bound;
+	}
+    };
+
+    template <>
+    struct ForBound<1>
+    {
+	template <size_t BOUND>
+	static constexpr size_t bound(size_t bound)
+	{
+	    return BOUND;
+	}
+    };
+    
+    template <class IndexClass, class Expr, ForType FT = ForType::DEFAULT>
+    class For
+    {
+    private:
+	For() = default;
+	//For(const For& in) = default;
+	//For& operator=(const For& in) = default;
+		
+	const IndexClass* mIndPtr;
+	size_t mSPos;
+	size_t mMax;
+	Expr mExpr;
+
+	typedef decltype(mExpr.rootSteps()) ExtType;
+	ExtType mExt;
+
+    public:
+
+	static constexpr size_t LAYER = Expr::LAYER + 1;
+	static constexpr size_t SIZE = Expr::SIZE;
+
+	For(const For& in) = default;
+	For& operator=(const For& in) = default;
+	For(For&& in) = default;
+	For& operator=(For&& in) = default;
+	
+	For(const std::shared_ptr<IndexClass>& indPtr,
+	    Expr expr);
+
+	For(const IndexClass* indPtr,
+	    Expr expr);
+
+	
+	inline void operator()(size_t mlast, ExtType last) const;
+	inline void operator()(size_t mlast = 0) const;
+
+	auto rootSteps(std::intptr_t iPtrNum = 0) const -> ExtType;
+	    
+    };
+
+    template <size_t N>
+    size_t exceptMax(size_t max) { return max; }
+
+    template <>
+    size_t exceptMax<1>(size_t max) { return 1; }
+
+} // namespace MultiArrayHelper
+
+/* ========================= *
+ * ---   TEMPLATE CODE   --- *
+ * ========================= */
+
+#include <iostream>
+
+namespace MultiArrayHelper
+{
+
+    template <class IndexClass, class Expr, ForType FT>
+    For<IndexClass,Expr,FT>::For(const std::shared_ptr<IndexClass>& indPtr,
+				 Expr expr) :
+	mIndPtr(indPtr.get()), mSPos(mIndPtr->pos()), mMax(mIndPtr->max()), mExpr(expr),
+	mExt(expr.rootSteps( reinterpret_cast<std::intptr_t>( mIndPtr.get() )))
+    {
+	assert(mIndPtr != nullptr);
+	//VCHECK(mIndPtr->id());
+	//VCHECK(mIndPtr->max());
+    }
+
+    template <class IndexClass, class Expr, ForType FT>
+    For<IndexClass,Expr,FT>::For(const IndexClass* indPtr,
+				 Expr expr) :
+	mIndPtr(indPtr), mSPos(mIndPtr->pos()), mMax(mIndPtr->max()),
+	mExpr(std::forward<Expr>( expr )),
+	mExt(expr.rootSteps( reinterpret_cast<std::intptr_t>( mIndPtr ) ))
+    {
+	assert(mIndPtr != nullptr);
+	//VCHECK(mIndPtr->id());
+	//VCHECK(mIndPtr->max());
+    }
+    
+    template <class IndexClass, class Expr, ForType FT>
+    inline void For<IndexClass,Expr,FT>::operator()(size_t mlast,
+						    ExtType last) const
+    {
+	typedef typename IndexClass::RangeType RangeType;
+	for(size_t pos = 0u; pos != ForBound<RangeType::ISSTATIC>::template bound<RangeType::SIZE>(mMax); ++pos){
+	//for(size_t pos = mSPos; pos != mMax; ++pos){
+	    const size_t mnpos = PosForward<FT>::value(mlast, mMax, pos);
+	    const ExtType npos = last + mExt*pos;
+	    mExpr(mnpos, npos);
+	}
+    }
+
+    template <class IndexClass, class Expr, ForType FT>
+    inline void For<IndexClass,Expr,FT>::operator()(size_t mlast) const
+    {
+	typedef typename IndexClass::RangeType RangeType;
+	const ExtType last;
+	for(size_t pos = 0u; pos != ForBound<RangeType::ISSTATIC>::template bound<RangeType::SIZE>(mMax); ++pos){
+	//for(size_t pos = mSPos; pos != mMax; ++pos){
+	    const size_t mnpos = PosForward<FT>::value(mlast, mMax, pos);
+	    const ExtType npos = last + mExt*pos;
+	    mExpr(mnpos, npos);
+	}
+    }
+    
+    template <class IndexClass, class Expr, ForType FT>
+    auto For<IndexClass,Expr,FT>::rootSteps(std::intptr_t iPtrNum) const
+	-> ExtType
+    {
+	return mExpr.rootSteps(iPtrNum);
+    }
+
+    
+} // namespace MultiArrayHelper
+
+#endif

src/lib/CMakeLists.txt

@@ -0,0 +1,11 @@
+
+set(libmultiarray_a_SOURCES
+  ${CMAKE_SOURCE_DIR}/src/lib/ranges/range_base.cc
+  ${CMAKE_SOURCE_DIR}/src/lib/ranges/index_info.cc
+  #${CMAKE_SOURCE_DIR}/src/lib/operation_utils.cc
+  )
+
+add_library(multiarray STATIC ${libmultiarray_a_SOURCES})
+install(TARGETS multiarray
+  ARCHIVE DESTINATION ${INSTALL_PATH}/lib
+  LIBRARY DESTINATION ${INSTALL_PATH}/lib)

src/lib/ranges/index_info.cc

@@ -1,6 +1,6 @@
 
-#include "index_info.h"
+#include "ranges/index_info.h"
-#include "range_base.h"
+#include "ranges/range_base.h"
 
 namespace MultiArrayTools
 {

src/lib/ranges/range_base.cc

@@ -1,5 +1,5 @@
 
-#include "range_base.h"
+#include "ranges/range_base.h"
 
 namespace MultiArrayTools
 {

src/operation_utils.cc

@@ -1,176 +0,0 @@
-
-#include "operation_utils.h"
-#include "ranges/range_base.h"
-
-namespace MultiArrayTools
-{
-
-    namespace
-    {
-	using namespace MultiArrayHelper;
-    }
-    
-    void seekIndexInst(std::shared_ptr<VIWB> i, std::vector<std::shared_ptr<VIWB> >& ivec)
-    {
-	std::cout << __func__ << ":" << std::endl;
-	for(size_t inum = 0; inum != i->rangePtr()->dim(); ++inum){
-	    std::cout << i->getPtrNum() << std::endl;
-	    auto ii = i->getPtr(inum);
-	    if(ii->type() == IndexType::MULTI or
-	       ii->type() == IndexType::CONT){
-		seekIndexInst(ii, ivec);
-	    }
-	    ivec.push_back(ii);
-	}
-	std::cout << std::endl;
-    }
-
-    void seekIndexInst(const IndexInfo* i, std::vector<const IndexInfo*>& ivec)
-    {
-	std::cout << __func__ << ":" << std::endl;
-	for(size_t inum = 0; inum != i->dim(); ++inum){
-	    std::cout << i->getPtrNum() << std::endl;
-	    auto ii = i->getPtr(inum);
-	    if(ii->type() == IndexType::MULTI or
-	       ii->type() == IndexType::CONT){
-		seekIndexInst(ii, ivec);
-	    }
-	    ivec.push_back(ii);
-	}
-	std::cout << std::endl;
-    }
-    
-    BTSS getBlockType(std::shared_ptr<VIWB> i,
-		      std::shared_ptr<VIWB> j,
-		      bool first, size_t higherStepSize)
-    {
-	// returning BlockType and step size is redundant (change in the future)
-	// stepSize == 0 => VALUE
-	// stepSize == 1 => BLOCK
-	// stepSize > 1 => SPLIT :)
-	BTSS out(BlockType::VALUE, 0);
-	size_t lastNum = i->rangePtr()->dim();
-
-	for(size_t inum = 0; inum != lastNum; ++inum){
-	    auto ii = i->getPtr(inum);
-	    if(ii->getPtrNum() == j->getPtrNum()){
-		
-		if(inum == lastNum - 1 and first){
-		    out = BTSS(BlockType::BLOCK, 1);
-		}
-		else {
-		    first = false;
-		    out = BTSS(BlockType::SPLIT, i->getStepSize(inum) * higherStepSize + out.second);
-		}
-		continue;
-	    }
-	    
-	    if(ii->type() == IndexType::MULTI or
-	       ii->type() == IndexType::CONT){
-
-		BTSS tmp = getBlockType(ii, j, inum == lastNum - 1, i->getStepSize(inum) * higherStepSize);
-		if(tmp.first != BlockType::VALUE){
-		    out = tmp;
-		}
-	    }
-	}
-	return out;
-    }
-
-    BTSS getBlockType(const IndexInfo* i,
-		      const IndexInfo* j,
-		      bool first, size_t higherStepSize)
-    {
-	// returning BlockType and step size is redundant (change in the future)
-	// stepSize == 0 => VALUE
-	// stepSize == 1 => BLOCK
-	// stepSize > 1 => SPLIT :)
-	BTSS out(BlockType::VALUE, 0);
-	size_t lastNum = i->dim();
-
-	for(size_t inum = 0; inum != lastNum; ++inum){
-	    auto ii = i->getPtr(inum);
-	    if(ii->getPtrNum() == j->getPtrNum()){
-		
-		if(inum == lastNum - 1 and first){
-		    out = BTSS(BlockType::BLOCK, 1);
-		}
-		else {
-		    first = false;
-		    out = BTSS(BlockType::SPLIT, i->getStepSize(inum) * higherStepSize + out.second);
-		}
-		continue;
-	    }
-	    
-	    if(ii->type() == IndexType::MULTI or
-	       ii->type() == IndexType::CONT){
-
-		BTSS tmp = getBlockType(ii, j, inum == lastNum - 1, i->getStepSize(inum) * higherStepSize);
-		if(tmp.first != BlockType::VALUE){
-		    out = tmp;
-		}
-	    }
-	}
-	return out;
-    }
-    
-    size_t getBTNum(const std::vector<BTSS>& mp, BlockType bt)
-    {
-	size_t out = 0;
-	for(auto& xx: mp){
-	    if(xx.first == bt){
-		++out;
-	    }
-	}
-	return out;
-    }
-
-    void minimizeAppearanceOfType(std::map<std::shared_ptr<VIWB>, std::vector<BTSS> >& mp,
-				  BlockType bt)
-    {
-	size_t minNum = getBTNum( mp.begin()->second, bt );
-	for(auto& mm: mp){
-	    size_t tmp = getBTNum( mm.second, bt );
-	    if(tmp < minNum){
-		minNum = tmp;
-	    }
-	}
-
-	for(auto mit = mp.begin(); mit != mp.end(); ){
-	    size_t tmp = getBTNum( mit->second, bt );
-	    if(tmp > minNum){
-		mit = mp.erase(mit);
-	    }
-	    else {
-		++mit;
-	    }
-	}
-
-    }
-
-    void minimizeAppearanceOfType(std::map<const IndexInfo*, std::vector<BTSS> >& mp,
-				  BlockType bt)
-    {
-	size_t minNum = getBTNum( mp.begin()->second, bt );
-	for(auto& mm: mp){
-	    size_t tmp = getBTNum( mm.second, bt );
-	    if(tmp < minNum){
-		minNum = tmp;
-	    }
-	}
-
-	for(auto mit = mp.begin(); mit != mp.end(); ){
-	    size_t tmp = getBTNum( mit->second, bt );
-	    if(tmp > minNum){
-		mit = mp.erase(mit);
-	    }
-	    else {
-		++mit;
-	    }
-	}
-
-    }
-
-
-
-} // end namespace MultiArrayTools

src/operation_utils.h

@@ -1,52 +0,0 @@
-
-#ifndef __operation_utils_h__
-#define __operation_utils_h__
-
-#include <utility>
-#include <vector>
-#include <memory>
-#include <map>
-
-#include "block/block.h"
-#include "ranges/vindex_base.h"
-
-#include "ranges/index_info.h"
-
-namespace MultiArrayTools
-{
-
-    namespace
-    {
-	using namespace MultiArrayHelper;
-    }
-    
-    // <block type, step size within actual instance>
-    typedef std::pair<BlockType,size_t> BTSS;
-
-    
-    void seekIndexInst(std::shared_ptr<VIWB> i, std::vector<std::shared_ptr<VIWB> >& ivec);
-
-    void seekIndexInst(const IndexInfo* i, std::vector<const IndexInfo*>& ivec);
-    
-    BTSS getBlockType(std::shared_ptr<VIWB> i,
-		      std::shared_ptr<VIWB> j,
-		      bool first, size_t higherStepSize = 1);
-
-    BTSS getBlockType(const IndexInfo* i,
-		      const IndexInfo* j,
-		      bool first, size_t higherStepSize = 1);
-
-
-    size_t getBTNum(const std::vector<BTSS>& mp, BlockType bt);
-
-    void minimizeAppearanceOfType(std::map<std::shared_ptr<VIWB>, std::vector<BTSS> >& mp,
-				  BlockType bt);
-
-    void minimizeAppearanceOfType(std::map<const IndexInfo*, std::vector<BTSS> >& mp,
-				  BlockType bt);
-
-    
-    
-} // end namespace MultiArrayTools
-
-#endif

src/ranges/tests/index_unit_test.cc~

@@ -1,216 +0,0 @@
-// -*- C++ -*-
-
-#include <cstdlib>
-#include "gtest/gtest.h"
-#include <iostream>
-
-#include "rheader.h"
-
-//#include "multi_array_header.h"
-
-namespace MAT = MultiArrayTools;
-
-namespace {
-
-    using namespace MAT;
-
-    template <class Factory, typename T>
-    void swapFactory(std::shared_ptr<RangeFactoryBase>& fptr, std::initializer_list<T> ilist)
-    {
-	std::vector<T> tmp = ilist;
-	auto nptr = std::make_shared<Factory>( tmp );
-	fptr = nptr;
-    }
-
-    template <class Factory, class... Rs>
-    void swapMFactory(std::shared_ptr<RangeFactoryBase>& fptr, const Rs&... rs)
-    {
-	auto nptr = std::make_shared<Factory>( rs... );
-	fptr = nptr;
-    }
-
-    template <typename... Ts>
-    auto mkt(Ts&&... ts) -> decltype(std::make_tuple(ts...))
-    {
-	return std::make_tuple(ts...);
-    }
-    
-    class IndexTest : public ::testing::Test
-    {
-    protected:
-
-	typedef SingleRangeFactory<char,SpaceType::ANY> SRF;
-	typedef SRF::oType SRange;
-
-	typedef MultiRangeFactory<SRange,SRange,SRange> M3RF;
-	typedef M3RF::oType M3Range;
-
-	typedef MultiRangeFactory<SRange,M3Range,SRange> MasterRF;
-	typedef MasterRF::oType MasterRange;
-
-	typedef ContainerRangeFactory<M3Range,SRange> CRF;
-	typedef CRF::oType CRange;
-
-	IndexTest()
-	{
-	    swapFactory<SRF>(rfbptr, { 'e', 'b', 'n' } );
-	    sr1ptr = std::dynamic_pointer_cast<SRange>( rfbptr->create() );
-
-	    swapFactory<SRF>(rfbptr, { 'x', 'y', 'l', 'f' } );
-	    sr2ptr = std::dynamic_pointer_cast<SRange>( rfbptr->create() );
-
-	    swapFactory<SRF>(rfbptr, { 'a', 'b' } );
-	    std::shared_ptr<SRange> temp1 = std::dynamic_pointer_cast<SRange>( rfbptr->create() );
-	    swapFactory<SRF>(rfbptr, { '1' } );
-	    std::shared_ptr<SRange> temp2 = std::dynamic_pointer_cast<SRange>( rfbptr->create() );
-	    swapFactory<SRF>(rfbptr, { '0', '7' } );
-	    std::shared_ptr<SRange> temp3 = std::dynamic_pointer_cast<SRange>( rfbptr->create() );
-
-	    swapMFactory<M3RF>(rfbptr, temp1, temp2, temp3 );
-	    m3rptr = std::dynamic_pointer_cast<M3Range>( rfbptr->create() );
-
-	    swapMFactory<MasterRF>(rfbptr, sr1ptr, m3rptr, sr2ptr);
-	    mstrptr = std::dynamic_pointer_cast<MasterRange>( rfbptr->create() );
-
-	    swapMFactory<CRF>(rfbptr, m3rptr, sr2ptr);
-	    cr1ptr = std::dynamic_pointer_cast<CRange>( rfbptr->create() );
-
-	    swapMFactory<CRF>(rfbptr, m3rptr, sr1ptr);
-	    cr2ptr = std::dynamic_pointer_cast<CRange>( rfbptr->create() );
-	}
-
-	std::shared_ptr<RangeFactoryBase> rfbptr;
-	std::shared_ptr<SRange> sr1ptr;
-	std::shared_ptr<SRange> sr2ptr;
-	std::shared_ptr<M3Range> m3rptr;
-	std::shared_ptr<MasterRange> mstrptr;
-	std::shared_ptr<CRange> cr1ptr;
-	std::shared_ptr<CRange> cr2ptr;
-    };
-
-    TEST_F(IndexTest, SingleIndex_SimpleCall)
-    {
-	auto si = sr1ptr->begin();
-	EXPECT_EQ(si.max(), 3u);
-	EXPECT_EQ(si.pos(), 0u);
-	EXPECT_EQ(si.first(), true);
-	EXPECT_EQ(si.last(), false);
-	EXPECT_EQ(si.meta(), 'e');
-	si.at('n');
-	EXPECT_EQ(si.pos(), si.max()-1);
-	EXPECT_EQ(si.first(), false);
-	EXPECT_EQ(si.last(), true);
-	si = 1;
-	EXPECT_EQ(si.meta(), 'b');
-	++si;
-	EXPECT_EQ(si.meta(), 'n');
-	auto si2 = sr1ptr->end();
-	--si2;
-	EXPECT_EQ(si == si2, true);
-	EXPECT_EQ(si != si2, false);
-	auto si3 = sr2ptr->end();
-	--si3;
-	EXPECT_EQ(si == si3, false);
-	EXPECT_EQ(si != si3, true);
-    }
-
-    TEST_F(IndexTest, MultiIndex_SimpleCall)
-    {
-	auto mi = m3rptr->begin();
-	EXPECT_EQ(mi.max(), 4u);
-	EXPECT_EQ(mi.pos(), 0u);
-	EXPECT_EQ(mi.first(), true);
-	EXPECT_EQ(mi.last(), false);
-	EXPECT_EQ(mi.meta() == mkt('a','1','0'), true);
-	mi.at( mkt('b','1','7') );
-	EXPECT_EQ(mi.pos(), mi.max()-1);
-	EXPECT_EQ(mi.first(), false);
-	EXPECT_EQ(mi.last(), true);
-	mi = 1;
-	EXPECT_EQ(mi.meta() == mkt('a','1','7'), true);
-	++mi;
-	EXPECT_EQ(mi.meta() == mkt('b','1','0'), true);
-	++mi;
-	auto mi2 = m3rptr->end();
-	--mi2;
-	EXPECT_EQ(mi == mi2, true);
-	EXPECT_EQ(mi != mi2, false);
-	--mi2;
-	EXPECT_EQ(mi == mi2, false);
-	EXPECT_EQ(mi != mi2, true);
-	auto si = sr1ptr->begin();
-	EXPECT_EQ(si == mi, false);
-	EXPECT_EQ(si != mi, true);
-    }
-
-    TEST_F(IndexTest, MasterRange_Check)
-    {
-	EXPECT_EQ(mstrptr->size(), 48u);
-	EXPECT_EQ(mstrptr->template get<0>().size(), 3u);
-	EXPECT_EQ(mstrptr->template get<1>().size(), 4u);
-	EXPECT_EQ(mstrptr->template get<2>().size(), 4u);
-	EXPECT_EQ(mstrptr->dim(), 3);
-	EXPECT_EQ(mstrptr->template get<0>().dim(), 1u);
-	EXPECT_EQ(mstrptr->template get<1>().dim(), 3u);
-
-	auto mi = mstrptr->begin();
-	EXPECT_EQ(mi.meta() == mkt( 'e' , mkt('a', '1', '0') , 'x' ), true);
-	mi = mi.max()-1;
-	EXPECT_EQ(mi.meta() == mkt( 'n' , mkt('b', '1', '7') , 'f' ), true);
-	mi.template down<1>();
-	EXPECT_EQ(mi.meta() == mkt( 'n' , mkt('b', '1', '0') , 'f' ), true);
-	mi.template down<0>();
-	EXPECT_EQ(mi.meta() == mkt( 'b' , mkt('b', '1', '0') , 'f' ), true);
-	mi.template down<2>();
-	EXPECT_EQ(mi.meta() == mkt( 'b' , mkt('b', '1', '0') , 'l' ), true);
-	mi.template up<1>();
-	EXPECT_EQ(mi.meta() == mkt( 'b' , mkt('b', '1', '7') , 'l' ), true);
-
-	auto& subI = mi.template get<0>();
-	EXPECT_EQ(subI.meta(), 'b');
-	mi.template up<0>();
-	EXPECT_EQ(subI.meta(), 'n');
-	auto& subI2 = mi.template get<1>();
-	EXPECT_EQ(subI2.meta() == mkt('b', '1', '7'), true);
-	mi.template down<1>();
-	EXPECT_EQ(subI2.meta() == mkt('b', '1', '0'), true);
-    }
-
-    TEST_F(IndexTest, ContainerRange_Check)
-    {
-	EXPECT_EQ(cr1ptr->size(), 16u);
-	EXPECT_EQ(cr2ptr->size(), 12u);
-
-	auto mi = mstrptr->begin();
-	auto ci1 = cr1ptr->begin();
-	auto ci2 = cr2ptr->begin();
-
-	EXPECT_EQ(ci1.max(), 16u);
-	EXPECT_EQ(ci2.max(), 12u);
-	
-	ci1(mi.template getPtr<1>(), mi.template getPtr<2>());
-	ci2(mi.template getPtr<1>(), mi.template getPtr<0>());
-
-	EXPECT_EQ(ci1.pos(), 0u);
-	EXPECT_EQ(ci2.pos(), 0u);
-	++mi;
-	EXPECT_EQ(ci1().pos(), 1u);
-	EXPECT_EQ(ci2().pos(), 0u);
-	mi.template up<1>();
-	EXPECT_EQ(ci1().pos(), 5u);
-	EXPECT_EQ(ci2().pos(), 3u);
-	mi.template up<0>();
-	EXPECT_EQ(ci1().pos(), 5u);
-	EXPECT_EQ(ci2().pos(), 4u);
-	mi = mi.max()-1;
-	EXPECT_EQ(ci1().pos(), ci1.max()-1);
-	EXPECT_EQ(ci2().pos(), ci2.max()-1);
-    }
-    
-} // end namespace 
-
-int main(int argc, char** argv)
-{
-    ::testing::InitGoogleTest(&argc, argv);
-    return RUN_ALL_TESTS();
-}

src/tests/CMakeLists.txt

@@ -0,0 +1,21 @@
+
+#add_executable(iutest /ranges/index_unit_test.cc ${INDEX_CC_FILES})
+add_executable(iutest ranges/index_unit_test.cc)
+add_dependencies(iutest multiarray)
+target_link_libraries(iutest ${GTEST_BOTH_LIBRARIES} ${CMAKE_THREAD_LIBS_INIT} multiarray)
+add_test(NAME iutest COMMAND iutest)
+
+#add_executable(mautest src/tests/ma_unit_test.cc ${MA_CC_FILES})
+add_executable(mautest ma_unit_test.cc)
+add_dependencies(mautest multiarray)
+target_link_libraries(mautest ${GTEST_BOTH_LIBRARIES} ${CMAKE_THREAD_LIBS_INIT} multiarray)
+add_test(NAME mautest COMMAND mautest)
+
+add_executable(oputest op_unit_test.cc)
+add_dependencies(oputest multiarray)
+target_link_libraries(oputest ${GTEST_BOTH_LIBRARIES} ${CMAKE_THREAD_LIBS_INIT} multiarray)
+add_test(NAME oputest COMMAND oputest)
+
+add_executable(opptest op_perf_test.cc)
+add_dependencies(opptest multiarray)
+target_link_libraries(opptest multiarray)