im com while fixing compile errors... not finished...

src/base_def.h

@@ -117,43 +117,7 @@ namespace MultiArrayTools
     // multi_array.h
     template <typename T, class... SRanges>
     class MultiArray;
-
+   
-    // block.h
-    template <typename T>
-    class BlockBase;
-
-    // block.h
-    template <typename T>
-    class MutableBlockBase;
-
-    // block.h
-    template <typename T>
-    class Block;
-
-    // block.h
-    template <typename T>
-    class MBlock;
-
-    // block.h
-    template <typename T>
-    class BlockValue;
-
-    // block.h
-    template <typename T>
-    class MBlockValue;
-
-    // block.h
-    template <typename T>
-    class SplitBlock;
-
-    // block.h
-    template <typename T>
-    class MSplitBlock;
-
-    // block.h
-    template <typename T>
-    class BlockResult;
-    
     // multi_array_operation.h
     template <typename T>
     class OperationBase;
@@ -163,7 +127,7 @@ namespace MultiArrayTools
     class MutableOperationBase;
 
     // multi_array_operation.h
-    template <class OperationClass>
+    template <typename T, class OperationClass>
     class OperationTemplate;
     
     // multi_array_operation.h
@@ -222,4 +186,47 @@ namespace MultiArrayTools
     */
 }
 
+namespace MultiArrayHelper
+{
+    // block.h
+    enum class BlockType;
+    
+    // block.h
+    template <typename T>
+    class BlockBase;
+
+    // block.h
+    template <typename T>
+    class MutableBlockBase;
+
+    // block.h
+    template <typename T>
+    class Block;
+
+    // block.h
+    template <typename T>
+    class MBlock;
+
+    // block.h
+    template <typename T>
+    class BlockValue;
+
+    // block.h
+    template <typename T>
+    class MBlockValue;
+
+    // block.h
+    template <typename T>
+    class SplitBlock;
+
+    // block.h
+    template <typename T>
+    class MSplitBlock;
+
+    // block.h
+    template <typename T>
+    class BlockResult;
+
+}
+
 #endif

src/block.cc

@@ -20,7 +20,7 @@ namespace MultiArrayHelper
 
     template <typename T>
     template <class OpFunction>
-    BlockResult<T> BlockBase<T>::operate(const BlockBase& in)
+    BlockResult<T> BlockBase<T>::operate(const BlockBase<T>& in)
     {
 	assert(mSize == in.size());
 	OpFunction f;
@@ -30,27 +30,28 @@ namespace MultiArrayHelper
 	}
 	return res;
     }
-
+    
+    
     template <typename T>
-    BlockResult<T> BlockBase<T>::operator+(const BlockBase& in)
+    BlockResult<T> BlockBase<T>::operator+(const BlockBase<T>& in)
     {
 	return operate<std::plus<T> >(in);
     }
 
     template <typename T>
-    BlockResult<T> BlockBase<T>::operator-(const BlockBase& in)
+    BlockResult<T> BlockBase<T>::operator-(const BlockBase<T>& in)
     {
 	return operate<std::minus<T> >(in);
     }
     
     template <typename T>
-    BlockResult<T> BlockBase<T>::operator*(const BlockBase& in)
+    BlockResult<T> BlockBase<T>::operator*(const BlockBase<T>& in)
     {
 	return operate<std::multiplies<T> >(in);
     }
     
     template <typename T>
-    BlockResult<T> BlockBase<T>::operator/(const BlockBase& in)
+    BlockResult<T> BlockBase<T>::operator/(const BlockBase<T>& in)
     {
 	return operate<std::divides<T> >(in);
     }
@@ -61,6 +62,15 @@ namespace MultiArrayHelper
     
     template <typename T>
     MutableBlockBase<T>::MutableBlockBase(size_t size) : BlockBase<T>(size) {}
+
+    template <typename T>
+    MutableBlockBase<T>& MutableBlockBase<T>::operator=(const BlockBase<T>& in)
+    {
+	for(size_t i = 0; i != BlockBase<T>::mSize; ++i){
+	    (*this)[i] = in[i];
+	}
+	return *this;
+    }
     
     /*************
      *   Block   *
@@ -96,11 +106,11 @@ namespace MultiArrayHelper
      **************/
 
     template <typename T>
-    MBlock<T>::MBlock(const std::vector<T>& data,
+    MBlock<T>::MBlock(std::vector<T>& data,
-		    size_t begPos, size_t size) :
+		      size_t begPos, size_t size) :
 	MutableBlockBase<T>(size),
 	mBegPtr(data.data() + begPos) {}
-
+   
     template <typename T>
     BlockType MBlock<T>::type() const
     {
@@ -159,7 +169,7 @@ namespace MultiArrayHelper
      *******************/
 
     template <typename T>
-    MBlockValue<T>::MBlockValue(const T& val, size_t size) :
+    MBlockValue<T>::MBlockValue(T& val, size_t size) :
 	BlockBase<T>(size),
 	mVal(val) {}
 
@@ -223,8 +233,8 @@ namespace MultiArrayHelper
      *******************/
 
     template <typename T>
-    MSplitBlock<T>::MSplitBlock(const std::vector<T>& data, size_t begPos,
+    MSplitBlock<T>::MSplitBlock(std::vector<T>& data, size_t begPos,
-			      size_t stepSize, size_t size) :
+				size_t stepSize, size_t size) :
 	BlockBase<T>(size),
 	mStepSize(stepSize),
 	mBegPtr(data.data() + begPos) {}

src/block.h

@@ -1,7 +1,7 @@
 // -*- C++ -*-
 
-#ifndef __block_h__
+#ifndef __ma_block_h__
-#define __block_h__
+#define __ma_block_h__
 
 #include <cstdlib>
 #include <vector>
@@ -16,7 +16,7 @@ namespace MultiArrayHelper
 	BLOCK = 1,
 	VALUE = 2,
 	SPLIT = 3,
-	RESULT = 4,
+	RESULT = 4
     };
     
     // manage vectorization in the future !!
@@ -28,6 +28,9 @@ namespace MultiArrayHelper
 	BlockBase() = default;
 	BlockBase(size_t size);
 
+	BlockBase(BlockBase&& res) = default;
+	BlockBase& operator=(BlockBase&& res) = default;
+	
 	virtual BlockType type() const = 0;
 	
 	virtual size_t size() const;
@@ -37,7 +40,7 @@ namespace MultiArrayHelper
 	
 	template <class OpFunction>
 	BlockResult<T> operate(const BlockBase& in);
-
+	
 	BlockResult<T> operator+(const BlockBase& in);
 	BlockResult<T> operator-(const BlockBase& in);
 	BlockResult<T> operator*(const BlockBase& in);
@@ -54,8 +57,13 @@ namespace MultiArrayHelper
 
 	MutableBlockBase() = default;
 	MutableBlockBase(size_t size);
+
+	MutableBlockBase& operator=(const BlockBase<T>& in);
 	
-	virtual T& operator[](size_t pos) override;
+	MutableBlockBase(MutableBlockBase&& res) = default;
+	MutableBlockBase& operator=(MutableBlockBase&& res) = default;
+	
+	virtual T& operator[](size_t pos) = 0;
 	
     };
     
@@ -84,7 +92,7 @@ namespace MultiArrayHelper
 	virtual BlockType type() const override;
 	virtual const T& operator[](size_t pos) const override;
 	virtual T& operator[](size_t pos) override;
-	virtual Block& set(const T* nbeg) override;
+	virtual MBlock& set(const T* nbeg) override;
 	
     protected:
 	T* mBegPtr;
@@ -115,7 +123,7 @@ namespace MultiArrayHelper
 	virtual BlockType type() const override;
 	virtual const T& operator[](size_t pos) const override;
 	virtual T& operator[](size_t pos) override;
-	virtual BlockValue& set(const T* nbeg) override;
+	virtual MBlockValue& set(const T* nbeg) override;
 	
     protected:
 	T& mVal;
@@ -151,7 +159,7 @@ namespace MultiArrayHelper
 	virtual BlockType type() const override;
 	virtual const T& operator[](size_t pos) const override;
 	virtual T& operator[](size_t pos) override;
-	virtual SplitBlock& set(const T* nbeg) override;
+	virtual MSplitBlock& set(const T* nbeg) override;
 	
     protected:
 	size_t mStepSize;
@@ -166,6 +174,9 @@ namespace MultiArrayHelper
 	BlockResult() = default;
 	BlockResult(size_t size);
 
+	BlockResult(BlockResult&& res) = default;
+	BlockResult& operator=(BlockResult&& res) = default;
+	
 	virtual BlockType type() const override;
 	virtual const T& operator[](size_t pos) const override;
 	virtual T& operator[](size_t i) override;

src/container_range.cc

@@ -79,17 +79,17 @@ namespace MultiArrayTools
     }
 
     template <class... Indices>
-    size_t ContainerIndex<Indices...>::pp(std::shared_ptr<const IndexBase>& idxPtr)
+    size_t ContainerIndex<Indices...>::pp(std::shared_ptr<IndexBase>& idxPtr)
     {
-	size_t tmp = pp(mIPack, mBlockSizes, idxPtr);
+	size_t tmp = PackNum<sizeof...(Indices)-1>::pp(mIPack, mBlockSizes, idxPtr);
 	IB::mPos += tmp;
 	return tmp;
     }
 
     template <class... Indices>
-    size_t ContainerIndex<Indices...>::mm(std::shared_ptr<const IndexBase>& idxPtr)
+    size_t ContainerIndex<Indices...>::mm(std::shared_ptr<IndexBase>& idxPtr)
     {
-	size_t tmp = mm(mIPack, mBlockSizes, idxPtr);
+	size_t tmp = PackNum<sizeof...(Indices)-1>::mm(mIPack, mBlockSizes, idxPtr);
 	IB::mPos -= tmp;
 	return tmp;
     }
@@ -140,7 +140,7 @@ namespace MultiArrayTools
     }
 
     template <class... Indices>
-    std::shared_ptr<const IndexBase> ContainerIndex<Indices...>::getPtr(size_t n) const
+    std::shared_ptr<IndexBase> ContainerIndex<Indices...>::getPtr(size_t n) const
     {
 	if(n >= sizeof...(Indices)){
 	    assert(0);
@@ -149,7 +149,18 @@ namespace MultiArrayTools
 	ContainerIndex<Indices...> const* t = this;
 	return PackNum<sizeof...(Indices)-1>::getIndexPtr(*t, n);
     }
+
     
+    template <class... Indices>
+    size_t ContainerIndex<Indices...>::getStepSize(size_t n) const
+    {
+	if(n >= sizeof...(Indices)){
+	    assert(0);
+	    // throw !!
+	}
+	return mBlockSizes[n+1];
+    }
+	
     template <class... Indices>
     bool ContainerIndex<Indices...>::first() const
     {

src/container_range.h

@@ -28,7 +28,8 @@ namespace MultiArrayTools
     protected:
 	bool mExternControl = false;
 	IndexPack mIPack;
-    
+	std::array<size_t,sizeof...(Indices)+1> mBlockSizes; 
+	
     public:
 	ContainerIndex() = delete;
 
@@ -44,8 +45,8 @@ namespace MultiArrayTools
 	virtual ContainerIndex& operator--() override;
 	virtual ContainerIndex& operator=(size_t pos) override;
 
-	virtual size_t pp(std::shared_ptr<const IndexBase>& idxPtr) override;
+	virtual size_t pp(std::shared_ptr<IndexBase>& idxPtr) override;
-	virtual size_t mm(std::shared_ptr<const IndexBase>& idxPtr) override;
+	virtual size_t mm(std::shared_ptr<IndexBase>& idxPtr) override;
 	
 	virtual MetaType meta() const override;
 	virtual ContainerIndex& at(const MetaType& metaPos) override;
@@ -63,7 +64,8 @@ namespace MultiArrayTools
 	template <size_t N>
 	auto getPtr() const -> decltype( std::get<N>( mIPack ) )&;
 
-	virtual std::shared_ptr<const IndexBase> getPtr(size_t n) const override;
+	virtual std::shared_ptr<IndexBase> getPtr(size_t n) const override;
+	virtual size_t getStepSize(size_t n) const override;
 	
 	ContainerIndex& operator()(const std::shared_ptr<Indices>&... inds); // control via external indices
 

src/index_base.cc

@@ -39,12 +39,17 @@ namespace MultiArrayTools
 	return mLocked;
     }
     
-    IndexBase& IndexBase::lock(std::shared_ptr<const IndexBase>& idx)
+    IndexBase& IndexBase::lock(std::shared_ptr<IndexBase>& idx)
     {
 	mLocked = (idx.get() == this);
 	return *this;
     }
     */
+
+    std::shared_ptr<RangeBase> IndexBase::rangePtr() const
+    {
+	return mRangePtr;
+    }
     
     IndexBase::operator size_t() const
     {

src/index_base.h

@@ -44,8 +44,8 @@ namespace MultiArrayTools
 	virtual IndexBase& operator++() = 0;
 	virtual IndexBase& operator--() = 0;
 
-	virtual size_t pp(std::shared_ptr<const IndexBase>& idxPtr) = 0;
+	virtual size_t pp(std::shared_ptr<IndexBase>& idxPtr) = 0;
-	virtual size_t mm(std::shared_ptr<const IndexBase>& idxPtr) = 0;
+	virtual size_t mm(std::shared_ptr<IndexBase>& idxPtr) = 0;
 	
 	bool operator==(const IndexBase& in) const;
 	bool operator!=(const IndexBase& in) const;
@@ -58,9 +58,12 @@ namespace MultiArrayTools
 	virtual bool first() const = 0;
 	
 	//virtual bool locked() const;
-	//virtual IndexBase& lock(std::shared_ptr<const IndexBase>& idx);
+	//virtual IndexBase& lock(std::shared_ptr<IndexBase>& idx);
-	
+
-	virtual std::shared_ptr<const IndexBase> getPtr(size_t n) const = 0;
+	virtual std::shared_ptr<RangeBase> rangePtr() const;
+	virtual std::shared_ptr<IndexBase> getPtr(size_t n) const = 0;
+
+	virtual size_t getStepSize(size_t n) const = 0;
 	
 	virtual operator size_t() const;
 

src/multi_array_header.h

@@ -10,6 +10,7 @@
 #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 "slice.h"

src/multi_array_operation.cc

@@ -9,9 +9,9 @@ namespace MultiArrayTools
 	using namespace MultiArrayHelper;
     }
 
-    void seekIndexInst(std::shared_ptr<const IndexBase> i, std::vector<std::shared_ptr<const IndexBase> >& ivec)
+    void seekIndexInst(std::shared_ptr<IndexBase> i, std::vector<std::shared_ptr<IndexBase> >& ivec)
     {
-	for(size_t inum = 0; inum != i->range()->dim(); ++inum){
+	for(size_t inum = 0; inum != i->rangePtr()->dim(); ++inum){
 	    auto ii = i->getPtr(inum);
 	    if(ii->type() == IndexType::MULTI){
 		seekIndexInst(ii, ivec);
@@ -20,16 +20,15 @@ namespace MultiArrayTools
 	}
     }
 
-    // !!!
+    BTSS getBlockType(std::shared_ptr<IndexBase> i,
-    BTSS getBlockType(std::shared_ptr<const IndexBase> i,
+		      std::shared_ptr<IndexBase> j, bool first)
-					     std::shared_ptr<const IndexBase> j, bool first)
     {
 	// 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);
-	for(size_t inum = 0; inum != i->range()->dim(); ++inum){
+	for(size_t inum = 0; inum != i->rangePtr()->dim(); ++inum){
 	    auto ii = i->getPtr(inum);
 	    if(ii == j){
 		if(inum == 0 and first){
@@ -79,6 +78,62 @@ namespace MultiArrayTools
 	    return std::make_shared<MSplitBlock<T> >(vec, 0, stepSize, blockSize);
 	}
     }
+
+    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<IndexBase>, 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;
+	    }
+	}
+
+    }
+    
+    template <typename T>
+    std::shared_ptr<IndexBase> seekBlockIndex(std::shared_ptr<IndexBase> ownIdx,
+						    const OperationBase<T>& second)
+    {
+	std::vector<std::shared_ptr<IndexBase> > ivec;
+	seekIndexInst(ownIdx, ivec);
+	std::map<std::shared_ptr<IndexBase>, std::vector<BTSS> > mp;
+	
+	for(auto& xx: ivec){
+	    mp[xx] = second.block(xx);
+	}
+
+	// seek minimal number of VALUEs => guarantees absence of conflicting blocks
+	minimizeAppearanceOfType(mp, BlockType::VALUE);
+
+	// seek mininmal number of SPLITs => maximal vectorization possible
+	minimizeAppearanceOfType(mp, BlockType::SPLIT);
+	
+	return mp.begin()->first;
+    }
+
     
     /*********************************
      *   MultiArrayOperationBase     *
@@ -91,39 +146,39 @@ namespace MultiArrayTools
      *   OperationTemplate     *
      ***************************/
 
-    template <class OperationClass>
+    template <typename T, class OperationClass>
-    OperationTemplate<OperationClass>::OperationTemplate(OperationClass* oc) : mOc(oc) {}
+    OperationTemplate<T,OperationClass>::OperationTemplate(OperationClass* oc) : mOc(oc) {}
     
-    template <class OperationClass>
+    template <typename T, class OperationClass>
     template <class Second>
-    auto OperationTemplate<OperationClass>::operator+(const Second& in) const
+    auto OperationTemplate<T,OperationClass>::operator+(const Second& in) const
-	-> Operation<value_type,std::plus<value_type>,OperationClass,Second>
+	-> Operation<T,std::plus<T>,OperationClass,Second>
     {
-    	return Operation<value_type,std::plus<value_type>,OperationClass,Second>(*mOc, in);
+    	return Operation<T,std::plus<T>,OperationClass,Second>(*mOc, in);
     }
 
-    template <class OperationClass>
+    template <typename T, class OperationClass>
     template <class Second>
-    auto OperationTemplate<OperationClass>::operator-(const Second& in) const
+    auto OperationTemplate<T,OperationClass>::operator-(const Second& in) const
-	-> Operation<value_type,std::minus<value_type>,OperationClass,Second>
+	-> Operation<T,std::minus<T>,OperationClass,Second>
     {
-    	return Operation<value_type,std::minus<value_type>,OperationClass,Second>(*mOc, in);
+    	return Operation<T,std::minus<T>,OperationClass,Second>(*mOc, in);
     }
     
-    template <class OperationClass>
+    template <typename T, class OperationClass>
     template <class Second>
-    auto OperationTemplate<OperationClass>::operator*(const Second& in) const
+    auto OperationTemplate<T,OperationClass>::operator*(const Second& in) const
-	-> Operation<value_type,std::multiplies<value_type>,OperationClass,Second>
+	-> Operation<T,std::multiplies<T>,OperationClass,Second>
     {
-    	return Operation<value_type,std::multiplies<value_type>,OperationClass,Second>(*mOc, in);
+    	return Operation<T,std::multiplies<T>,OperationClass,Second>(*mOc, in);
     }
 
-    template <class OperationClass>
+    template <typename T, class OperationClass>
     template <class Second>
-    auto OperationTemplate<OperationClass>::operator/(const Second& in) const
+    auto OperationTemplate<T,OperationClass>::operator/(const Second& in) const
-	-> Operation<double,std::divides<value_type>,OperationClass,Second>
+	-> Operation<T,std::divides<T>,OperationClass,Second>
     {
-    	return Operation<value_type,std::divides<value_type>,OperationClass,Second>(*mOc, in);
+    	return Operation<T,std::divides<T>,OperationClass,Second>(*mOc, in);
     }
 
     /*************************
@@ -142,16 +197,16 @@ namespace MultiArrayTools
 	mIndex = std::make_shared<IndexType>( mr->begin() );
 	(*mIndex) = *index;
 
-	auto blockIndex = seekBlockIndex(mIndex, second);
+	auto blockIndex = seekBlockIndex( mIndex, second);
 	block(blockIndex);
 	second.block(blockIndex);
-	for(*mIndex = 0; mIndex->pos() != mIndex->max(); mIndex->pp(blockIndex) )){
+	for(*mIndex = 0; mIndex->pos() != mIndex->max(); mIndex->pp(blockIndex) ){
 	    get() = mSecond.get();
 	}
     }
     
     template <typename T, class... Ranges>
-    BlockBase<T>& OperationMaster<T,Ranges...>::get()
+    MutableBlockBase<T>& OperationMaster<T,Ranges...>::get()
     {
 	block();
 	return *mBlockPtr;
@@ -165,7 +220,7 @@ namespace MultiArrayTools
     }
 
     template <typename T, class... Ranges>
-    std::vector<BTSS> OperationMaster<T,Ranges...>::block(const std::shared_ptr<IndexBase>& blockIndex) const
+    std::vector<BTSS> OperationMaster<T,Ranges...>::block(const std::shared_ptr<IndexBase> blockIndex) const
     {
 	std::vector<BTSS> btv(1, getBlockType(mIndex, blockIndex, true) );
 	mBlockPtr = makeBlock(mArrayRef.data(), btv[0].second, blockIndex->max());
@@ -187,7 +242,7 @@ namespace MultiArrayTools
     ConstOperationRoot<T,Ranges...>::
     ConstOperationRoot(const MultiArrayBase<T,Ranges...>& ma,
 		       const std::shared_ptr<typename Ranges::IndexType>&... indices) :
-	OperationBase<T>(), OperationTemplate<ConstOperationRoot<T,Ranges...> >(this),
+	OperationBase<T>(), OperationTemplate<T,ConstOperationRoot<T,Ranges...> >(this),
 	mArrayRef(ma), mIndex( std::make_shared<IndexType>( mArrayRef.range() ) )
     {
 	mIndex(indices...);
@@ -201,7 +256,7 @@ namespace MultiArrayTools
     }
 
     template <typename T, class... Ranges>
-    std::vector<BTSS> ConstOperationRoot<T,Ranges...>::block(const std::shared_ptr<IndexBase>& blockIndex) const
+    std::vector<BTSS> ConstOperationRoot<T,Ranges...>::block(const std::shared_ptr<IndexBase> blockIndex) const
     {
 	std::vector<BTSS> btv(1, getBlockType(mIndex, blockIndex, true) );
 	mBlockPtr = makeBlock(mArrayRef.data(), btv[0].second, blockIndex->max());
@@ -223,7 +278,7 @@ namespace MultiArrayTools
     OperationRoot<T,Ranges...>::
     OperationRoot(MutableMultiArrayBase<T,Ranges...>& ma,
 		  const std::shared_ptr<typename Ranges::IndexType>&... indices) :
-	MutableOperationBase<T>(), OperationTemplate<OperationRoot<T,Ranges...> >(this),
+	MutableOperationBase<T>(), OperationTemplate<T,OperationRoot<T,Ranges...> >(this),
 	mArrayRef(ma), mIndex( std::make_shared<IndexType>( mArrayRef.range() ) )
     {
 	(*mIndex)(indices...);
@@ -243,14 +298,14 @@ namespace MultiArrayTools
     }
 
     template <typename T, class... Ranges>
-    BlockBase<T>& OperationRoot<T,Ranges...>::get()
+    MutableBlockBase<T>& OperationRoot<T,Ranges...>::get()
     {
 	block();
 	return *mBlockPtr;
     }
 
     template <typename T, class... Ranges>
-    std::vector<BTSS> OperationRoot<T,Ranges...>::block(const std::shared_ptr<IndexBase>& blockIndex) const
+    std::vector<BTSS> OperationRoot<T,Ranges...>::block(const std::shared_ptr<IndexBase> blockIndex) const
     {
 	std::vector<BTSS> btv(1, getBlockType(mIndex, blockIndex, true) );
 	mBlockPtr = makeBlock(mArrayRef.data(), btv[0].second, blockIndex->max());
@@ -270,7 +325,7 @@ namespace MultiArrayTools
     
     template <typename T, class OpFunction, class... Ops>
     Operation<T,OpFunction,Ops...>::Operation(const Ops&... ops) :
-	OperationBase<T>(), OperationTemplate<Operation<T,OpFunction,Ops...> >(this),
+	OperationBase<T>(), OperationTemplate<T,Operation<T,OpFunction,Ops...> >(this),
 	mOps(ops...) {}
 
     template <typename T, class OpFunction, class... Ops>
@@ -280,18 +335,18 @@ namespace MultiArrayTools
 	return mRes;
     }
 
-    template <typename T, class... Ranges>
+    template <typename T, class OpFunction, class... Ops>
-    std::vector<BTSS> Operation<T,Ranges...>::block(const std::shared_ptr<IndexBase>& blockIndex) const
+    std::vector<BTSS> Operation<T,OpFunction,Ops...>::block(const std::shared_ptr<IndexBase> blockIndex) const
     {
 	std::vector<BTSS> btv;
-	PackNum<sizeof...(Ranges)-1>::makeBlockTypeVec(btv, mOps, blockIndex);
+	PackNum<sizeof...(Ops)-1>::makeBlockTypeVec(btv, mOps, blockIndex);
 	return btv;
     }
 
-    template <typename T, class... Ranges>
+    template <typename T, class OpFunction, class... Ops>
-    Operation<T,Ranges...>& Operation<T,Ranges...>::block() const
+    Operation<T,OpFunction,Ops...>& Operation<T,OpFunction,Ops...>::block() const
     {
-	mBlockPtr->set( &mArrayRef[ (*mIndex)() ] );
+	//mBlockPtr->set( &mArrayRef[ (*mIndex)() ] );
 	return *this;
     }
 

src/multi_array_operation.h

@@ -32,14 +32,14 @@ namespace MultiArrayTools
      * 
      */
 
-    void seekIndexInst(std::shared_ptr<const IndexBase> i, std::vector<std::shared_ptr<const IndexBase> >& ivec);
+    void seekIndexInst(std::shared_ptr<IndexBase> i, std::vector<std::shared_ptr<IndexBase> >& ivec);
 
     
     // <block type, step size within actual instance>
     typedef std::pair<BlockType,size_t> BTSS;
     
-    BTSS getBlockType(std::shared_ptr<const IndexBase> i,
+    BTSS getBlockType(std::shared_ptr<IndexBase> i,
-		      std::shared_ptr<const IndexBase> j, bool first);
+		      std::shared_ptr<IndexBase> j, bool first);
 
     template <typename T>
     std::shared_ptr<BlockBase<T> > makeBlock(const std::vector<T>& vec, size_t stepSize, size_t blockSize);
@@ -47,73 +47,26 @@ namespace MultiArrayTools
     template <typename T>
     std::shared_ptr<MutableBlockBase<T> > makeBlock(std::vector<T>& vec, size_t stepSize, size_t blockSize);
 
-    size_t getBTNum(const std::vector<BTSS>& mp, BlockType bt)
+    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<const IndexBase>, std::vector<BTSS> > mp,
+    void minimizeAppearanceOfType(std::map<std::shared_ptr<IndexBase>, std::vector<BTSS> > mp,
-				  BlockType bt)
+				  BlockType bt);
-    {
-	size_t minNum = getBTNum( *mp.begin(), 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;
-	    }
-	}
-
-    }
-    
     template <typename T>
-    std::shared_ptr<const IndexBase> seekBlockIndex(std::shared_ptr<const IndexBase>& ownIdx,
+    std::shared_ptr<IndexBase> seekBlockIndex(std::shared_ptr<IndexBase> ownIdx,
-						    const OperationBase<T>& second)
+					      const OperationBase<T>& second);
-    {
-	std::vector<std::shared_ptr<const IndexBase> > ivec;
-	seekIndexInst(ownIdx, ivec);
-	std::map<std::shared_ptr<const IndexBase>, std::vector<BTSS> > mp;
-	
-	for(auto& xx: ivec){
-	    mp[xx] = second.block(xx);
-	}
-
-	// seek minimal number of VALUEs => guarantees absence of conflicting blocks
-	minimizeAppearanceOfType(mp, BlockType::VALUE);
-
-	// seek mininmal number of SPLITs => maximal vectorization possible
-	minimizeAppearanceOfType(mp, BlockType::SPLIT);
-	
-	return *mp.begin();
-    }
     
     template <typename T>
     class OperationBase
     {
     public:
-
 	typedef T value_type;
-	
+
 	OperationBase() = default;
 	virtual ~OperationBase() = default;
 
 	 // init block, return resulting type (BLOCK, VALUE, SPLIT)
-	virtual std::vector<BlockType> block(const std::shared_ptr<IndexBase>& blockIndex) const = 0;
+	virtual std::vector<BTSS> block(const std::shared_ptr<IndexBase> blockIndex) const = 0;
 	virtual OperationBase& block() const = 0; // update block
 	
 	//virtual size_t argNum() const = 0;
@@ -127,33 +80,33 @@ namespace MultiArrayTools
 	typedef T value_type;
 	
 	MutableOperationBase() = default;
-	virtual BlockBase<T>& get() = 0;
+
+	virtual MutableBlockBase<T>& get() = 0;
+
     };
 
-    template <class OperationClass>
+    template <typename T, class OperationClass>
     class OperationTemplate
     {
     public:
-
-	typedef typename OperationClass::value_type value_type;
 	
 	OperationTemplate(OperationClass* oc);
 	
 	template <class Second>
 	auto operator+(const Second& in) const
-	    -> Operation<value_type,std::plus<value_type>,OperationClass,Second>;
+	    -> Operation<T,std::plus<T>,OperationClass,Second>;
 	
 	template <class Second>
 	auto operator-(const Second& in) const
-	    -> Operation<value_type,std::minus<value_type>,OperationClass,Second>;
+	    -> Operation<T,std::minus<T>,OperationClass,Second>;
 	
 	template <class Second>
 	auto operator*(const Second& in) const
-	    -> Operation<value_type,std::multiplies<value_type>,OperationClass,Second>;
+	    -> Operation<T,std::multiplies<T>,OperationClass,Second>;
 	
 	template <class Second>
 	auto operator/(const Second& in) const
-	    -> Operation<value_type,std::divides<value_type>,OperationClass,Second>;
+	    -> Operation<T,std::divides<T>,OperationClass,Second>;
 	
     private:
 	OperationClass* mOc;
@@ -175,7 +128,7 @@ namespace MultiArrayTools
 	virtual MutableBlockBase<T>& get() override;
 	virtual const BlockBase<T>& get() const override;
 
-	virtual std::vector<BlockType> block(const std::shared_ptr<IndexBase>& blockIndex) const override;
+	virtual std::vector<BTSS> block(const std::shared_ptr<IndexBase> blockIndex) const override;
 	virtual OperationMaster& block() const override;
 	
     protected:
@@ -184,19 +137,19 @@ namespace MultiArrayTools
 	OperationBase<T> const& mSecond;
 	MutableMultiArrayBase<T,Ranges...>& mArrayRef;
 	std::shared_ptr<IndexType> mIndex;
-	mutable std::shared_ptr<MutableBlockBase> mBlockPtr;
+	mutable std::shared_ptr<MutableBlockBase<T> > mBlockPtr;
     };
 
     
     template <typename T, class... Ranges>
     class ConstOperationRoot : public OperationBase<T>,
-			       public OperationTemplate<ConstOperationRoot<T,Ranges...> >
+			       public OperationTemplate<T,ConstOperationRoot<T,Ranges...> >
     {
     public:
 
 	typedef T value_type;
 	typedef OperationBase<T> OB;
-	typedef OperationTemplate<ConstOperationRoot<T,Ranges...> > OT;
+	typedef OperationTemplate<T,ConstOperationRoot<T,Ranges...> > OT;
 	typedef ContainerRange<Ranges...> CRange;
 	typedef typename CRange::IndexType IndexType;
 	
@@ -205,25 +158,25 @@ namespace MultiArrayTools
 	
 	virtual const BlockBase<T>& get() const override;
 
-	virtual std::vector<BlockType> block(const std::shared_ptr<IndexBase>& blockIndex) const override;
+	virtual std::vector<BTSS> block(const std::shared_ptr<IndexBase> blockIndex) const override;
 	virtual ConstOperationRoot& block() const override;
 	
     protected:
 	
 	MultiArrayBase<T,Ranges...> const& mArrayRef;
 	std::shared_ptr<IndexType> mIndex;
- 	mutable std::shared_ptr<BlockBase> mBlockPtr;
+ 	mutable std::shared_ptr<BlockBase<T> > mBlockPtr;
     };
 
     template <typename T, class... Ranges>
     class OperationRoot : public MutableOperationBase<T>,
-			  public OperationTemplate<OperationRoot<T,Ranges...> >
+			  public OperationTemplate<T,OperationRoot<T,Ranges...> >
     {
     public:
 
 	typedef T value_type;
 	typedef OperationBase<T> OB;
-	typedef OperationTemplate<OperationRoot<T,Ranges...> > OT;
+	typedef OperationTemplate<T,OperationRoot<T,Ranges...> > OT;
 	typedef ContainerRange<Ranges...> CRange;
 	typedef typename CRange::IndexType IndexType;
 	
@@ -235,32 +188,32 @@ namespace MultiArrayTools
 	virtual const BlockBase<T>& get() const override;
 	virtual MutableBlockBase<T>& get() override;
 
-	virtual std::vector<BlockType> block(const std::shared_ptr<IndexBase>& blockIndex) const override;
+	virtual std::vector<BTSS> block(const std::shared_ptr<IndexBase> blockIndex) const override;
 	virtual OperationRoot& block() const override;
 	
     protected:
 	
 	MutableMultiArrayBase<T,Ranges...>& mArrayRef;
 	std::shared_ptr<IndexType> mIndex;
-	mutable std::shared_ptr<MutableBlockBase> mBlockPtr;
+	mutable std::shared_ptr<MutableBlockBase<T> > mBlockPtr;
     };
 
     template <typename T, class OpFunction, class... Ops>
     class Operation : public OperationBase<T>,
-		      public OperationTemplate<Operation<T,OpFunction,Ops...> >
+		      public OperationTemplate<T,Operation<T,OpFunction,Ops...> >
     {
     public:
 
 	typedef T value_type;
 	typedef OperationBase<T> OB;
-	typedef OperationTemplate<Operation<T,OpFunction,Ops...> > OT;
+	typedef OperationTemplate<T,Operation<T,OpFunction,Ops...> > OT;
 	typedef OpFunction F;
 	
 	Operation(const Ops&... ops);
 	
 	virtual const BlockBase<T>& get() const override;
 
-	virtual std::vector<BlockType> block(const std::shared_ptr<IndexBase>& blockIndex) const override;
+	virtual std::vector<BTSS> block(const std::shared_ptr<IndexBase> blockIndex) const override;
 	virtual Operation& block() const override;
 	
     protected:

src/multi_range.cc

@@ -81,17 +81,17 @@ namespace MultiArrayTools
     }
 
     template <class... Indices>
-    size_t MultiIndex<Indices...>::pp(std::shared_ptr<const IndexBase>& idxPtr)
+    size_t MultiIndex<Indices...>::pp(std::shared_ptr<IndexBase>& idxPtr)
     {
-	size_t tmp = pp(mIPack, mBlockSizes, idxPtr);
+	size_t tmp = PackNum<sizeof...(Indices)-1>::pp(mIPack, mBlockSizes, idxPtr);
 	IB::mPos += tmp;
 	return tmp;
     }
 
     template <class... Indices>
-    size_t MultiIndex<Indices...>::mm(std::shared_ptr<const IndexBase>& idxPtr)
+    size_t MultiIndex<Indices...>::mm(std::shared_ptr<IndexBase>& idxPtr)
     {
-	size_t tmp = mm(mIPack, mBlockSizes, idxPtr);
+	size_t tmp = PackNum<sizeof...(Indices)-1>::mm(mIPack, mBlockSizes, idxPtr);
 	IB::mPos -= tmp;
 	return tmp;
     }
@@ -148,7 +148,7 @@ namespace MultiArrayTools
     }
 
     template <class... Indices>
-    std::shared_ptr<const IndexBase> MultiIndex<Indices...>::getPtr(size_t n) const
+    std::shared_ptr<IndexBase> MultiIndex<Indices...>::getPtr(size_t n) const
     {
 	if(n >= sizeof...(Indices)){
 	    assert(0);
@@ -158,7 +158,8 @@ namespace MultiArrayTools
 	return PackNum<sizeof...(Indices)-1>::getIndexPtr(*t, n);
     }
 
-    size_t getStepSize(size_t n) const
+    template <class... Indices>
+    size_t MultiIndex<Indices...>::getStepSize(size_t n) const
     {
 	if(n >= sizeof...(Indices)){
 	    assert(0);
@@ -202,7 +203,7 @@ namespace MultiArrayTools
     }
     /*
     template <class... Indices>
-    MultiIndex<Indices...>& MultiIndex<Indices...>::lock(std::shared_ptr<const IndexBase>& idx)
+    MultiIndex<Indices...>& MultiIndex<Indices...>::lock(std::shared_ptr<IndexBase>& idx)
     {
 	IB::mLocked = (idx.get() == this);
 	PackNum<sizeof...(Indices)-1>::lock(mIPack, idx);

src/multi_range.h

@@ -46,8 +46,8 @@ namespace MultiArrayTools
 	virtual MultiIndex& operator--() override;
 	virtual MultiIndex& operator=(size_t pos) override;
 
-	virtual size_t pp(std::shared_ptr<const IndexBase>& idxPtr) override;
+	virtual size_t pp(std::shared_ptr<IndexBase>& idxPtr) override;
-	virtual size_t mm(std::shared_ptr<const IndexBase>& idxPtr) override;
+	virtual size_t mm(std::shared_ptr<IndexBase>& idxPtr) override;
 	
 	template <size_t DIR>
 	MultiIndex& up();
@@ -62,7 +62,7 @@ namespace MultiArrayTools
 	auto getPtr() const -> decltype( std::get<N>( mIPack ) )&;
 	
 	const IndexBase& get(size_t n) const;
-	virtual std::shared_ptr<const IndexBase> getPtr(size_t n) const override;
+	virtual std::shared_ptr<IndexBase> getPtr(size_t n) const override;
 	virtual size_t getStepSize(size_t n) const override;
 	
         virtual MetaType meta() const override;
@@ -75,7 +75,7 @@ namespace MultiArrayTools
 
 	std::shared_ptr<RangeType> range() const;
 
-	//virtual MultiIndex& lock(std::shared_ptr<const IndexBase>& idx) override;
+	//virtual MultiIndex& lock(std::shared_ptr<IndexBase>& idx) override;
 	
 	// raplace instances (in contrast to its analogon in ContainerIndex
 	// MultiIndices CANNOT be influences be its subindices, so there is

src/pack_num.h

@@ -52,7 +52,7 @@ namespace MultiArrayHelper
 
 	template <class... Indices>
 	static void lock(std::tuple<std::shared_ptr<Indices>...>& ip,
-			 std::shared_ptr<const IndexBase>& toLock)
+			 std::shared_ptr<IndexBase> toLock)
 	{
 	    std::get<N>(ip)->lock(toLock);
 	    PackNum<N-1>::lock(ip, toLock);
@@ -87,7 +87,7 @@ namespace MultiArrayHelper
 	template <class... Indices>
 	static inline size_t pp(std::tuple<std::shared_ptr<Indices>...>& ip,
 				std::array<size_t,sizeof...(Indices)>& bs,
-				std::shared_ptr<const IndexBase>& idxPtr)
+				std::shared_ptr<IndexBase> idxPtr)
 	{
 	    auto siPtr = std::get<N>(ip);
 	    if(siPtr.get() == idxPtr.get()){
@@ -109,7 +109,7 @@ namespace MultiArrayHelper
 	{
 	    auto& si = *std::get<N>(ip);
 	    if(si.first()){
-		si = si.max();
+		si = si.max() - 1;
 		PackNum<N-1>::mm(ip);
 	    }
 	    else {
@@ -120,7 +120,7 @@ namespace MultiArrayHelper
 	template <class... Indices>
 	static inline size_t mm(std::tuple<std::shared_ptr<Indices>...>& ip,
 				std::array<size_t,sizeof...(Indices)>& bs,
-				std::shared_ptr<const IndexBase>& idx)
+				std::shared_ptr<IndexBase> idxPtr)
 	{
 	    auto siPtr = std::get<N>(ip);
 	    if(siPtr.get() == idxPtr.get()){
@@ -128,10 +128,11 @@ namespace MultiArrayHelper
 	    }
 	    else {
 		if(siPtr->first()){
+		    (*siPtr) = siPtr->max() - 1;
 		    return PackNum<N-1>::mm(ip, bs) - siPtr->max() + 1;
 		}
 		else {
-		    return siPtr->pp(idx);
+		    return siPtr->mm(idxPtr);
 		}
 	    }
 	}
@@ -229,7 +230,7 @@ namespace MultiArrayHelper
 	template <class... Ops>
 	static void makeBlockTypeVec(std::vector<std::pair<BlockType,size_t> >& btv,
 				     const std::tuple<Ops...>& ops,
-				     std::shared_ptr<const IndexBase>& idxPtr)
+				     std::shared_ptr<IndexBase> idxPtr)
 	{
 	    auto& subvec = std::get<N>(ops).block(idxPtr);
 	    btv.insert(btv.end(), subvec.begin(), subvec.end() );
@@ -261,7 +262,7 @@ namespace MultiArrayHelper
 
 	template <class... Indices>
 	static void lock(std::tuple<std::shared_ptr<Indices>...>& ip,
-			 std::shared_ptr<const IndexBase>& toLock)
+			 std::shared_ptr<IndexBase> toLock)
 	{
 	    std::get<0>(ip)->lock(toLock);
 	}
@@ -282,21 +283,20 @@ namespace MultiArrayHelper
 	
 	template <class... Indices>
 	static inline size_t pp(std::tuple<std::shared_ptr<Indices>...>& ip,
-			      std::array<size_t,sizeof...(Indices)+1>& bs)
+				std::array<size_t,sizeof...(Indices)>& bs,
+				std::shared_ptr<IndexBase> idxPtr)
 	{
-	    auto& si = *std::get<0>(ip);
+	    auto siPtr = std::get<0>(ip);
-	    if(si.locked()){
+	    if(siPtr.get() == idxPtr.get()){
-		si = si.max()+1;
+		return std::get<0>(bs);
-		return std::get<0>(bs) - std::get<1>(bs) + 1;
 	    }
-	    else {		
+	    else {
-		++si;
+		return siPtr->pp(idxPtr);
-		return 1;
 	    }
 	}
 
 	template <class... Indices>
-	static inline size_t mm(std::tuple<std::shared_ptr<Indices>...>& ip)
+	static inline void mm(std::tuple<std::shared_ptr<Indices>...>& ip)
 	{
 	    auto& si = *std::get<0>(ip);
 	    --si;
@@ -304,17 +304,15 @@ namespace MultiArrayHelper
 	
 	template <class... Indices>
 	static inline size_t mm(std::tuple<std::shared_ptr<Indices>...>& ip,
-			      std::array<size_t,sizeof...(Indices)>& bs)
+				std::array<size_t,sizeof...(Indices)>& bs,
+				std::shared_ptr<IndexBase> idxPtr)
 	{
-	    auto& si = *std::get<0>(ip);
+	    auto siPtr = std::get<0>(ip);
-	    if(si.locked()){
+	    if(siPtr.get() == idxPtr.get()){
-		si = 0;
+		return std::get<0>(bs);
-		--si;
-		return std::get<0>(bs) - std::get<1>(bs) + 1;
 	    }
-	    else {		
+	    else {
-		--si;
+		return siPtr->mm(idxPtr);
-		return 1;
 	    }
 	}
 	
@@ -405,7 +403,7 @@ namespace MultiArrayHelper
 	template <class... Ops>
 	static void makeBlockTypeVec(std::vector<std::pair<BlockType,size_t> >& btv,
 				     const std::tuple<Ops...>& ops,
-				     std::shared_ptr<const IndexBase>& idxPtr)
+				     std::shared_ptr<IndexBase> idxPtr)
 	{
 	    auto& subvec = std::get<0>(ops).block(idxPtr);
 	    btv.insert(btv.end(), subvec.begin(), subvec.end() );

src/single_range.cc

@@ -39,14 +39,14 @@ namespace MultiArrayTools
     }
 
     template <typename U, RangeType TYPE>
-    size_t SingleIndex<U,TYPE>::pp(std::shared_ptr<const IndexBase>& idxPtr)
+    size_t SingleIndex<U,TYPE>::pp(std::shared_ptr<IndexBase>& idxPtr)
     {
 	++(*this);
 	return 1;
     }
 
     template <typename U, RangeType TYPE>
-    size_t SingleIndex<U,TYPE>::mm(std::shared_ptr<const IndexBase>& idxPtr)
+    size_t SingleIndex<U,TYPE>::mm(std::shared_ptr<IndexBase>& idxPtr)
     {
 	--(*this);
 	return 1;
@@ -84,9 +84,15 @@ namespace MultiArrayTools
     }
 
     template <typename U, RangeType TYPE>
-    std::shared_ptr<const IndexBase> SingleIndex<U,TYPE>::getPtr(size_t n) const
+    std::shared_ptr<IndexBase> SingleIndex<U,TYPE>::getPtr(size_t n) const
     {
-	return std::shared_ptr<const IndexBase>();
+	return std::shared_ptr<IndexBase>();
+    }
+
+    template <typename U, RangeType TYPE>    
+    size_t SingleIndex<U,TYPE>::getStepSize(size_t n) const
+    {
+	return 1;
     }
     
     /********************

src/single_range.h

@@ -33,8 +33,8 @@ namespace MultiArrayTools
 	virtual SingleIndex& operator++() override;
 	virtual SingleIndex& operator--() override;
 
-	virtual size_t pp(std::shared_ptr<const IndexBase>& idxPtr) override;
+	virtual size_t pp(std::shared_ptr<IndexBase>& idxPtr) override;
-	virtual size_t mm(std::shared_ptr<const IndexBase>& idxPtr) override;
+	virtual size_t mm(std::shared_ptr<IndexBase>& idxPtr) override;
 	
 	virtual U meta() const override;
 	virtual SingleIndex& at(const U& metaPos) override;
@@ -43,7 +43,8 @@ namespace MultiArrayTools
 	virtual bool last() const override;
 	virtual bool first() const override;
 
-	virtual std::shared_ptr<const IndexBase> getPtr(size_t n) const override;
+	virtual std::shared_ptr<IndexBase> getPtr(size_t n) const override;
+	virtual size_t getStepSize(size_t n) const override;
 	
 	virtual std::string id() const override { return std::string("sin") + std::to_string(IB::mId); }
     };