replace VIWB by IndexInfo in blocking routines -> remaining performance lack in get()

src/multi_array_operation.h

@@ -47,10 +47,16 @@ namespace MultiArrayTools
     template <typename T>
     MBlock<T> makeBlock(T* vec, size_t stepSize, size_t blockSize);
 
+    // dont use this for now !!
     template <class OpClass>
     std::shared_ptr<VIWB> seekBlockIndex(std::shared_ptr<VIWB> ownIdx,
 					 const OpClass& second);
 
+    template <class OpClass>
+    const IndexInfo* seekBlockIndex(const IndexInfo* ownII,
+				    const OpClass& second);
+
+    
     template <typename T, class OperationClass>
     class OperationTemplate
     {
@@ -99,12 +105,13 @@ namespace MultiArrayTools
 
 	OperationMaster(MutableMultiArrayBase<T,Ranges...>& ma, const OpClass& second,
 			std::shared_ptr<typename CRange::IndexType>& index,
-			std::shared_ptr<VIWB> blockIndex);
+			const IndexInfo* blockIndex);
-
+	
 	MBlock<T>& get();
 	const Block<T>& get() const;
 
 	std::vector<BTSS> block(const std::shared_ptr<VIWB> blockIndex, bool init = false) const;
+	std::vector<BTSS> block(const IndexInfo* blockIndex, bool init = false) const;
 	const OperationMaster& block() const;
 	
     protected:
@@ -137,6 +144,7 @@ namespace MultiArrayTools
 	const Block<T>& get() const;
 
 	std::vector<BTSS> block(const std::shared_ptr<VIWB> blockIndex, bool init = false) const;
+	std::vector<BTSS> block(const IndexInfo* blockIndex, bool init = false) const;
 	const ConstOperationRoot& block() const;
 	
     protected:
@@ -171,10 +179,12 @@ namespace MultiArrayTools
 	
 	const MBlock<T>& get() const;
 	MBlock<T>& get();
-	
+
+	OperationRoot& set(const IndexInfo* blockIndex);
 	OperationRoot& set(std::shared_ptr<VIWB> blockIndex);
 	
 	std::vector<BTSS> block(const std::shared_ptr<VIWB> blockIndex, bool init = false) const;
+	std::vector<BTSS> block(const IndexInfo* blockIndex, bool init = false) const;
 	const OperationRoot& block() const;
 	
     protected:
@@ -188,6 +198,7 @@ namespace MultiArrayTools
 	IndexInfo mIInfo;
 	mutable MBlock<T> mBlock;
 	std::shared_ptr<VIWB> mBlockIndex; // predefine to save time
+	const IndexInfo* mBlockII; // predefine to save time
     };
     
     template <typename T, class OpFunction, class... Ops>
@@ -206,6 +217,7 @@ namespace MultiArrayTools
 	const BlockResult<T>& get() const;
 
 	std::vector<BTSS> block(const std::shared_ptr<VIWB> blockIndex, bool init = false) const;
+	std::vector<BTSS> block(const IndexInfo* blockIndex, bool init = false) const;
 	const Operation& block() const;
 	
     protected:
@@ -226,6 +238,7 @@ namespace MultiArrayTools
 	const BlockResult<T>& get() const;
 
 	std::vector<BTSS> block(const std::shared_ptr<VIWB> blockIndex, bool init = false) const;
+	std::vector<BTSS> block(const IndexInfo* blockIndex, bool init = false) const;
 	const Contraction& block() const;
 	
     protected:
@@ -260,10 +273,12 @@ namespace MultiArrayTools
 	return MBlock<T>(vec, 0, blockSize, stepSize);
     }
 
+    // dont use this for now !!
     template <class OpClass>
     std::shared_ptr<VIWB> seekBlockIndex(std::shared_ptr<VIWB> ownIdx,
 					 const OpClass& second)
     {
+	assert(0); // dont use this for now !!
 	std::vector<std::shared_ptr<VIWB> > ivec;
 	seekIndexInst(ownIdx, ivec);
 	std::map<std::shared_ptr<VIWB>, std::vector<BTSS> > mp;
@@ -279,6 +294,20 @@ namespace MultiArrayTools
 	return mp.begin()->first;
     }
 
+    template <class OpClass>
+    const IndexInfo* seekBlockIndex(const IndexInfo* ownII,
+				    const OpClass& second)
+	
+    {
+	const IndexInfo* ii = ownII;
+	while(ii->type() == IndexType::CONT or
+	      ii->type() == IndexType::MULTI){
+	    ii = ii->getPtr(ii->dim()-1);
+	}
+	return ii;
+    }
+
+    
     /***************************
      *   OperationTemplate     *
      ***************************/
@@ -334,7 +363,7 @@ namespace MultiArrayTools
 		    std::shared_ptr<typename CRange::IndexType>& index) :
 	mSecond(second), mArrayRef(ma), mIndex(mkIndex(index)), mIInfo(*mIndex)
     {
-	auto blockIndex = seekBlockIndex( make_viwb( mIndex ), second);
+	auto blockIndex = seekBlockIndex( &mIInfo, second);
 	std::intptr_t blockIndexNum = blockIndex->getPtrNum();
 	
 	block(blockIndex, true);
@@ -347,7 +376,7 @@ namespace MultiArrayTools
     OperationMaster<T,OpClass,Ranges...>::
     OperationMaster(MutableMultiArrayBase<T,Ranges...>& ma, const OpClass& second,
 		    std::shared_ptr<typename CRange::IndexType>& index,
-		    std::shared_ptr<VIWB> blockIndex) :
+		    const IndexInfo* blockIndex) :
 	mSecond(second), mArrayRef(ma), mIndex(mkIndex(index)), mIInfo(*mIndex)
     {
 	std::intptr_t blockIndexNum = blockIndex->getPtrNum();
@@ -356,7 +385,6 @@ namespace MultiArrayTools
 	performAssignment(blockIndexNum);
     }
 
-
     template <typename T, class OpClass, class... Ranges>
     std::shared_ptr<typename OperationMaster<T,OpClass,Ranges...>::IndexType>
     OperationMaster<T,OpClass,Ranges...>::
@@ -411,6 +439,17 @@ namespace MultiArrayTools
 	return btv;
     }
 
+    template <typename T, class OpClass, class... Ranges>
+    std::vector<BTSS> OperationMaster<T,OpClass,Ranges...>::block(const IndexInfo* blockIndex, bool init) const
+    {
+	std::vector<BTSS> btv(1, getBlockType( &mIInfo, blockIndex, true) );
+	if(init){
+	    mBlock = makeBlock(mArrayRef.data(), btv[0].second, blockIndex->max());
+	}
+	return btv;
+    }
+
+    
     template <typename T, class OpClass, class... Ranges>
     const OperationMaster<T,OpClass,Ranges...>& OperationMaster<T,OpClass,Ranges...>::block() const
     {
@@ -451,6 +490,7 @@ namespace MultiArrayTools
     template <typename T, class... Ranges>
     std::vector<BTSS> ConstOperationRoot<T,Ranges...>::block(const std::shared_ptr<VIWB> blockIndex, bool init) const
     {
+	assert(0);
 	std::vector<BTSS> btv(1, getBlockType( make_viwb( mIndex ), blockIndex, true) );
 	if(init){
 	    mBlock = makeBlock(mArrayRef.data(), btv[0].second, blockIndex->max());
@@ -458,6 +498,16 @@ namespace MultiArrayTools
 	return btv;
     }
 
+    template <typename T, class... Ranges>
+    std::vector<BTSS> ConstOperationRoot<T,Ranges...>::block(const IndexInfo* blockIndex, bool init) const
+    {
+	std::vector<BTSS> btv(1, getBlockType( &mIInfo, blockIndex, true) );
+	if(init){
+	    mBlock = makeBlock(mArrayRef.data(), btv[0].second, blockIndex->max());
+	}
+	return btv;
+    }
+    
     template <typename T, class... Ranges>
     const ConstOperationRoot<T,Ranges...>& ConstOperationRoot<T,Ranges...>::block() const
     {
@@ -475,7 +525,7 @@ namespace MultiArrayTools
 		  const std::shared_ptr<typename Ranges::IndexType>&... indices) :
 	//OperationTemplate<T,OperationRoot<T,Ranges...> >(this),
 	mArrayRef(ma), mIndex( mkIndex( ma, indices... ) ), mIInfo(*mIndex),
-	mBlockIndex(nullptr)
+	mBlockIndex(nullptr), mBlockII(nullptr)
     {}
 
     template <typename T, class... Ranges>
@@ -493,8 +543,8 @@ namespace MultiArrayTools
     template <class OpClass>
     OperationMaster<T,OpClass,Ranges...> OperationRoot<T,Ranges...>::operator=(const OpClass& in)
     {
-	if(mBlockIndex){
+	if(mBlockII != nullptr){
-	    return OperationMaster<T,OpClass,Ranges...>(mArrayRef, in, mIndex, mBlockIndex);
+	    return OperationMaster<T,OpClass,Ranges...>(mArrayRef, in, mIndex, mBlockII);
 	}
 	else {
 	    return OperationMaster<T,OpClass,Ranges...>(mArrayRef, in, mIndex);
@@ -522,10 +572,20 @@ namespace MultiArrayTools
 	mBlockIndex = blockIndex;
 	return *this;
     }
+
+    template <typename T, class... Ranges>
+    OperationRoot<T,Ranges...>&
+    OperationRoot<T,Ranges...>::set(const IndexInfo* blockIndex)
+    {
+	mBlockII = blockIndex;
+	return *this;
+    }
+
     
     template <typename T, class... Ranges>
     std::vector<BTSS> OperationRoot<T,Ranges...>::block(const std::shared_ptr<VIWB> blockIndex, bool init) const
     {
+	assert(0);
 	std::vector<BTSS> btv(1, getBlockType( make_viwb( mIndex ), blockIndex, true) );
 	if(init){
 	    mBlock = makeBlock(mArrayRef.data(), btv[0].second, blockIndex->max());
@@ -533,6 +593,16 @@ namespace MultiArrayTools
 	return btv;
     }
 
+    template <typename T, class... Ranges>
+    std::vector<BTSS> OperationRoot<T,Ranges...>::block(const IndexInfo* blockIndex, bool init) const
+    {
+	std::vector<BTSS> btv(1, getBlockType( &mIInfo, blockIndex, true) );
+	if(init){
+	    mBlock = makeBlock(mArrayRef.data(), btv[0].second, blockIndex->max());
+	}
+	return btv;
+    }
+    
     template <typename T, class... Ranges>
     const OperationRoot<T,Ranges...>& OperationRoot<T,Ranges...>::block() const
     {
@@ -559,6 +629,7 @@ namespace MultiArrayTools
     template <typename T, class OpFunction, class... Ops>
     std::vector<BTSS> Operation<T,OpFunction,Ops...>::block(const std::shared_ptr<VIWB> blockIndex, bool init) const
     {
+	assert(0);
 	std::vector<BTSS> btv;
 	PackNum<sizeof...(Ops)-1>::makeBlockTypeVec(btv, mOps, blockIndex, init);
 	if(init){
@@ -567,6 +638,18 @@ namespace MultiArrayTools
 	return btv;
     }
 
+    template <typename T, class OpFunction, class... Ops>
+    std::vector<BTSS> Operation<T,OpFunction,Ops...>::block(const IndexInfo* blockIndex, bool init) const
+    {
+	std::vector<BTSS> btv;
+	PackNum<sizeof...(Ops)-1>::makeBlockTypeVec(btv, mOps, blockIndex, init);
+	if(init){
+	    mRes.init(blockIndex->max());
+	}
+	return btv;
+    }
+
+    
     template <typename T, class OpFunction, class... Ops>
     const Operation<T,OpFunction,Ops...>& Operation<T,OpFunction,Ops...>::block() const
     {
@@ -603,6 +686,15 @@ namespace MultiArrayTools
 	return mOp.block(blockIndex, init);
     }
 
+    template <typename T, class Op, class IndexType>
+    std::vector<BTSS> Contraction<T,Op,IndexType>::block(const IndexInfo* blockIndex, bool init) const
+    {
+	if(init){
+	    mRes.init(blockIndex->max());
+	}
+	return mOp.block(blockIndex, init);
+    }
+    
     template <typename T, class Op, class IndexType>
     const Contraction<T,Op,IndexType>& Contraction<T,Op,IndexType>::block() const
     {

src/operation_utils.cc

@@ -12,7 +12,9 @@ namespace MultiArrayTools
     
     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){
@@ -20,8 +22,24 @@ namespace MultiArrayTools
 	    }
 	    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)
@@ -59,6 +77,43 @@ namespace MultiArrayTools
 	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;
@@ -92,7 +147,30 @@ namespace MultiArrayTools
 	}
 
     }
-    
+
+    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

@@ -10,6 +10,8 @@
 #include "block/block.h"
 #include "ranges/vindex_base.h"
 
+#include "ranges/index_info.h"
+
 namespace MultiArrayTools
 {
 
@@ -24,18 +26,26 @@ namespace MultiArrayTools
     
     void seekIndexInst(std::shared_ptr<VIWB> i, std::vector<std::shared_ptr<VIWB> >& ivec);
 
-    //void seekIndexInst(const IndexInfo& i, std::vector<IndexInfo>& 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
 

src/pack_num.h

@@ -27,6 +27,16 @@ namespace MultiArrayHelper
 	    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)
 	{
@@ -62,7 +72,16 @@ namespace MultiArrayHelper
 	    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){
 	    out << std::get<sizeof...(T)-1>(tp);

src/ranges/container_range.h

@@ -53,7 +53,7 @@ namespace MultiArrayTools
 
 	// ==== >>>>> STATIC POLYMORPHISM <<<<< ====
 	
-	IndexType type();
+	IndexType type() const;
 	
 	ContainerIndex& operator++();
 	ContainerIndex& operator--();
@@ -216,7 +216,7 @@ namespace MultiArrayTools
     }
 
     template <class... Indices>
-    IndexType ContainerIndex<Indices...>::type() { return IndexType::CONT; }
+    IndexType ContainerIndex<Indices...>::type() const { return IndexType::CONT; }
 
     template <class... Indices>
     ContainerIndex<Indices...>& ContainerIndex<Indices...>::operator++()

src/ranges/index_info.h

@@ -64,7 +64,8 @@ namespace MultiArrayTools
 	mPtrNum( reinterpret_cast<std::intptr_t>( &ind ) ),
 	mDim(ind.vrange()->dim()),
 	mMax(ind.max()),
-	mStepSize(stepSize)
+	mStepSize(stepSize),
+	mType(ind.type())
     {}
 
     template <class IndexClass>

src/ranges/multi_range.h

@@ -70,7 +70,7 @@ namespace MultiArrayTools
 
 	// ==== >>>>> STATIC POLYMORPHISM <<<<< ====
 	
-	IndexType type();
+	IndexType type() const;
 
 	MultiIndex& operator=(size_t pos);
 
@@ -264,7 +264,7 @@ namespace MultiArrayTools
     }
 
     template <class... Indices>
-    IndexType MultiIndex<Indices...>::type()
+    IndexType MultiIndex<Indices...>::type() const
     {
 	return IndexType::MULTI;
     }

src/ranges/rpack_num.h

@@ -166,10 +166,10 @@ namespace MultiArrayHelper
 	}
 	
 	template <class... Indices>
-	static size_t makePos(const std::tuple<std::shared_ptr<Indices>...>& iPtrTup)
+	static inline size_t makePos(const std::tuple<std::shared_ptr<Indices>...>& iPtrTup)
 	{
-	    const auto& idx = *std::get<N>(iPtrTup);
+	    //const auto& idx = *std::get<N>(iPtrTup);
-	    return idx.pos() + RPackNum<N-1>::makePos(iPtrTup) * idx.max();
+	    return std::get<N>(iPtrTup)->pos() + RPackNum<N-1>::makePos(iPtrTup) * std::get<N>(iPtrTup)->max();
 	}
 
 	template <class Pack, class IndexType, class... Indices>
@@ -337,7 +337,7 @@ namespace MultiArrayHelper
 	}
 	
 	template <class... Indices>
-	static size_t makePos(const std::tuple<std::shared_ptr<Indices>...>& iPtrTup)
+	static inline size_t makePos(const std::tuple<std::shared_ptr<Indices>...>& iPtrTup)
 	{
 	    return std::get<0>(iPtrTup)->pos();
 	}

src/ranges/single_range.h

@@ -29,7 +29,7 @@ namespace MultiArrayTools
 
 	// ==== >>>>> STATIC POLYMORPHISM <<<<< ====
 	
-	IndexType type();
+	IndexType type() const;
 	
 	SingleIndex& operator=(size_t pos);
 	SingleIndex& operator++();
@@ -121,7 +121,7 @@ namespace MultiArrayTools
 	IndexInterface<SingleIndex<U,TYPE>,U>(range, 0) {}
 
     template <typename U, SpaceType TYPE>
-    IndexType SingleIndex<U,TYPE>::type()
+    IndexType SingleIndex<U,TYPE>::type() const
     {
 	return IndexType::SINGLE;
     }

src/tests/op_perf_test.cc

@@ -167,7 +167,8 @@ namespace {
 	auto mix = MAT::mkMIndex( alpha, beta, gamma );
 
 	std::clock_t begin = std::clock();
-	res1(delta, deltap).set(vdeltap) = ma(delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix);
+	//res1(delta, deltap).set(vdeltap) = ma(delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix);
+	res1(delta, deltap) = ma(delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix);
 	std::clock_t end = std::clock();
 	std::cout << "MultiArray time: " << static_cast<double>( end - begin ) / CLOCKS_PER_SEC
 		  << std::endl;

src/tests/op_unit_test.cc

@@ -220,7 +220,8 @@ namespace {
 	auto mix = MAT::mkMIndex( alpha, beta, gamma );
 
 	std::clock_t begin = std::clock();
-	res1(delta, deltap).set(vdeltap) = ma(delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix);
+	//res1(delta, deltap).set(vdeltap) = ma(delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix);
+	res1(delta, deltap) = ma(delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix);
 	std::clock_t end = std::clock();
 	std::cout << "MultiArray time: " << static_cast<double>( end - begin ) / CLOCKS_PER_SEC
 		  << std::endl;