clean up + switch loops -> auto vectorizing seems to work; TODO: re-enable step sizes...

src/arith.h

@@ -8,40 +8,39 @@ namespace MultiArrayHelper
     template <typename T>
     struct plus
     {
-	static inline T& acc(T& target, const T& arg)
+	static inline T apply(T a1, T a2)
 	{
-	    return target += arg;
+	    return a1 + a2;
 	}
     };
 
     template <typename T>
     struct minus
     {
-	static inline T& acc(T& target, const T& arg)
+	static inline T apply(T a1, T a2)
 	{
-	    return target -= arg;
+	    return a1 - a2;
 	}
     };
 
     template <typename T>
     struct multiplies
     {
-	static inline T& acc(T& target, const T& arg)
+	static inline T apply(T a1, T a2)
 	{
-	    return target *= arg;
+	    return a1 * a2;
 	}
     };
 
     template <typename T>
     struct divides
     {
-	static inline T& acc(T& target, const T& arg)
+	static inline T apply(T a1, T a2)
 	{
-	    return target /= arg;
+	    return a1 / a2;
 	}
     };
 
-
 } // end namespace MultiArrayHelper
     
 #endif

src/multi_array_operation.h

@@ -28,22 +28,6 @@ namespace MultiArrayTools
 	using namespace MultiArrayHelper;
     }
     
-    template <typename T>
-    Block<T> makeBlock(const T* vec, size_t stepSize, size_t blockSize);
-    
-    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
     {
@@ -86,14 +70,11 @@ namespace MultiArrayTools
 	{
 	private:
 	    AssignmentExpr() = default;
-	    //AssignmentExpr(const AssignmentExpr& in) = default;
-	    //AssignmentExpr& operator=(const AssignmentExpr& in) = default;
 	    	    
 	    OperationMaster& mM;
 	    const OpClass& mSec;
 
 	public:
-	    static size_t layer() { return 0; }
 
 	    static constexpr size_t LAYER = 0;
 	    static constexpr size_t SIZE = OpClass::SIZE;
@@ -103,7 +84,6 @@ namespace MultiArrayTools
 
 	    AssignmentExpr(const AssignmentExpr& in) = default;
 	    AssignmentExpr(AssignmentExpr&& in) = default;
-	    //AssignmentExpr& operator=(const AssignmentExpr&& in) = default;
 	    
 	    inline void operator()(size_t start = 0) const; 
 	    inline void operator()(size_t start, ExtType last) const;
@@ -118,8 +98,6 @@ namespace MultiArrayTools
 	typedef typename MultiRange<Ranges...>::IndexType IndexType;
 	typedef MBlock<T> bType;	
 
-	static size_t rootNum() { return 1; }
-	
 	OperationMaster(MutableMultiArrayBase<T,Ranges...>& ma, const OpClass& second,
 			std::shared_ptr<typename CRange::IndexType>& index);
 
@@ -127,16 +105,9 @@ namespace MultiArrayTools
 			std::shared_ptr<typename CRange::IndexType>& index,
 			const IndexInfo* blockIndex);
 
-	//MBlock<T>& get();
-	//const Block<T>& get() const;
-
 	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;
-	//inline const T& get(size_t pos) const;
-	
-	//std::vector<BTSS> block(const IndexInfo* blockIndex, bool init = false) const;
-	//const OperationMaster& block() const;
 
     private:
 
@@ -162,22 +133,15 @@ namespace MultiArrayTools
 	typedef OperationTemplate<T,ConstOperationRoot<T,Ranges...> > OT;
 	typedef ContainerRange<Ranges...> CRange;
 	typedef typename CRange::IndexType IndexType;
-	typedef Block<T> bType;
 
-	static size_t rootNum() { return 1; }
+	static constexpr size_t SIZE = 1;
-	static const size_t SIZE = 1;
 	
 	ConstOperationRoot(const MultiArrayBase<T,Ranges...>& ma,
 			   const std::shared_ptr<typename Ranges::IndexType>&... indices);
 	
-	//const Block<T>& get() const;
-	
 	template <class ET>
 	inline T get(ET pos) const;
 
-	//std::vector<BTSS> block(const IndexInfo* blockIndex, bool init = false) const;
-	//const ConstOperationRoot& block() const;
-
 	MExt<void> rootSteps(std::intptr_t iPtrNum = 0) const; // nullptr for simple usage with decltype
 
 	template <class Expr>
@@ -193,7 +157,6 @@ namespace MultiArrayTools
 	const T* mDataPtr;
 	std::shared_ptr<IndexType> mIndex;
 	IndexInfo mIInfo;
- 	//mutable bType mBlock;
     };
     
     template <typename T, class... Ranges>
@@ -206,10 +169,8 @@ namespace MultiArrayTools
 	typedef OperationTemplate<T,OperationRoot<T,Ranges...> > OT;
 	typedef ContainerRange<Ranges...> CRange;
 	typedef typename CRange::IndexType IndexType;
-	typedef MBlock<T> bType;
 
-	static size_t rootNum() { return 1; }
+	static constexpr size_t SIZE = 1;
-	static const size_t SIZE = 1;
 	
 	OperationRoot(MutableMultiArrayBase<T,Ranges...>& ma,
 		      const std::shared_ptr<typename Ranges::IndexType>&... indices);
@@ -217,19 +178,9 @@ namespace MultiArrayTools
 	template <class OpClass>
 	OperationMaster<T,OpClass,Ranges...> operator=(const OpClass& in);
 
-	//const MBlock<T>& get() const;
-	//MBlock<T>& get();
-
-	//template <class ET>
-	//inline const T& get(const ET& pos) const;
-
 	template <class ET>
 	inline T get(ET pos) const;
 
-	//OperationRoot& set(const IndexInfo* blockIndex);
-	//std::vector<BTSS> block(const IndexInfo* blockIndex, bool init = false) const;
-	//const OperationRoot& block() const;
-
 	MExt<void> rootSteps(std::intptr_t iPtrNum = 0) const; // nullptr for simple usage with decltype
 
 	template <class Expr>
@@ -245,8 +196,6 @@ namespace MultiArrayTools
 	T* mDataPtr;
 	std::shared_ptr<IndexType> mIndex;
 	IndexInfo mIInfo;
-	//mutable bType mBlock;
-	//const IndexInfo* mBlockII; // predefine to save time
     };
 
     template <class Op>
@@ -298,29 +247,20 @@ namespace MultiArrayTools
 	typedef OperationBase<T> OB;
 	typedef OperationTemplate<T,Operation<T,OpFunction,Ops...> > OT;
 	typedef OpFunction F;
-	typedef BlockResult<T> bType;
 	
-	static size_t rootNum() { return sumRootNum<Ops...>(); }
 	static constexpr size_t SIZE = RootSum<Ops...>::SIZE;
 
     private:
 	std::tuple<Ops...> mOps;
-	mutable bType mRes;
 
     public:
 	typedef decltype(PackNum<sizeof...(Ops)-1>::mkSteps(0, mOps)) ETuple;
-	//typedef decltype(PackNum<sizeof...(Ops)-1>::template mkLoopType<Ops...>) LType;
 		
 	Operation(const Ops&... ops);
 	
-	//const BlockResult<T>& get() const;
-	
 	template <class ET>
 	inline T get(ET pos) const;
 
-	//std::vector<BTSS> block(const IndexInfo* blockIndex, bool init = false) const;
-	//const Operation& block() const;
-
 	auto rootSteps(std::intptr_t iPtrNum = 0) const // nullptr for simple usage with decltype
 	    -> decltype(PackNum<sizeof...(Ops)-1>::mkSteps(iPtrNum, mOps));
 
@@ -337,30 +277,22 @@ namespace MultiArrayTools
 
 	typedef T value_type;
 	typedef OperationTemplate<T,Contraction<T,Op,IndexType> > OT;
-	typedef BlockResult<T> bType;
 
-	static size_t rootNum() { return typename Op::rootNum(); }
+	static constexpr size_t SIZE = Op::SIZE;
-	static const size_t SIZE = Op::SIZE;
 	
     private:
 
 	const Op& mOp;
 	std::shared_ptr<IndexType> mInd;
-	mutable bType mRes;
 
     public:
 	typedef decltype(mOp.rootSteps(0)) ETuple;
 	
 	Contraction(const Op& op, std::shared_ptr<IndexType> ind);
 
-	//const BlockResult<T>& get() const;
-
 	template <class ET>
 	inline T get(ET pos) const;
 
-	//std::vector<BTSS> block(const IndexInfo* blockIndex, bool init = false) const;
-	//const Contraction& block() const;
-
 	auto rootSteps(std::intptr_t iPtrNum = 0) const  // nullptr for simple usage with decltype
 	    -> decltype(mOp.rootSteps(iPtrNum));
 
@@ -381,53 +313,6 @@ namespace MultiArrayTools
 	using namespace MultiArrayHelper;
     }
     
-    template <typename T>
-    Block<T> makeBlock(const T* vec, size_t stepSize, size_t blockSize)
-    {
-	return Block<T>(vec, 0, blockSize, stepSize);
-    }
-
-    template <typename T>
-    MBlock<T> makeBlock(T* vec, size_t stepSize, size_t blockSize)
-    {
-	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;
-	
-	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;
-    }
-
-    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     *
      ***************************/
@@ -481,14 +366,6 @@ namespace MultiArrayTools
     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)
-    {
-	
-    }*/
-
     template <typename T, class OpClass, class... Ranges>
     inline void OperationMaster<T,OpClass,Ranges...>::AssignmentExpr::
     operator()(size_t start, ExtType last) const
@@ -516,12 +393,6 @@ namespace MultiArrayTools
 	mSecond(second), mArrayRef(ma), mDataPtr(mArrayRef.data()),
 	mIndex(mkIndex(index)), mIInfo(*mIndex)
     {
-	//auto blockIndex = seekBlockIndex( &mIInfo, second);
-	//std::intptr_t blockIndexNum = blockIndex->getPtrNum();
-	
-	//block(blockIndex, true);
-	//second.block(blockIndex, true);
-
 	performAssignment(0);
     }
 
@@ -533,9 +404,6 @@ namespace MultiArrayTools
 	mSecond(second), mArrayRef(ma), mDataPtr(mArrayRef.data()),
 	mIndex(mkIndex(index)), mIInfo(*mIndex)
     {
-	//std::intptr_t blockIndexNum = blockIndex->getPtrNum();
-	//second.block(blockIndex, true);
-
 	performAssignment(0);
     }
 
@@ -555,71 +423,18 @@ namespace MultiArrayTools
     template <typename T, class OpClass, class... Ranges>
     void OperationMaster<T,OpClass,Ranges...>::performAssignment(std::intptr_t blockIndexNum)
     {
-#define XX_USE_NEW_LOOP_ROUTINE_XX
-#ifdef XX_USE_NEW_LOOP_ROUTINE_XX
-	// === N E W ===
 	AssignmentExpr ae(*this, mSecond); // Expression to be executed within loop
-	//const auto hiddenLoop = mSecond.loop(AssignmentExpr(*this, mSecond)); // hidden loop within 'mSecond' e.g. contractions
+	const auto loop = mSecond.template loop<decltype(mIndex->ifor(ae))>( mIndex->ifor(ae) );
-	auto loop = mIndex->ifor
+	// hidden Loops outside ! -> auto vectorizable
-	    ( mSecond.template loop<AssignmentExpr>
-	      ( std::move(ae) ) ); // init overall loop(s) 
 	loop(); // execute overall loop(s) and so internal hidden loops and so the inherited expressions
-#else
-	// === O L D ===
-	for(*mIndex = 0; mIndex->pos() != mIndex->max(); mIndex->pp(blockIndexNum) ){
-	    block();
-	    get() = mSecond.get();
-	}
-#endif
-    }
-    /*
-    template <typename T, class OpClass, class... Ranges>
-    MBlock<T>& OperationMaster<T,OpClass,Ranges...>::get()
-    {
-	return mBlock;
     }
 
-    template <typename T, class OpClass, class... Ranges>
-    const Block<T>& OperationMaster<T,OpClass,Ranges...>::get() const
-    {
-	return mBlock;
-    }
-    */
     template <typename T, class OpClass, class... Ranges>
     inline T OperationMaster<T,OpClass,Ranges...>::get(size_t pos) const
     {
-	//assert(pos < mIndex->max());
-	//if(pos >= mIndex->max()) { VCHECK(pos); VCHECK(mIndex->max()); assert(0); }
-	//VCHECK(pos);
 	return mDataPtr[pos];
     }
 
-    /*
-    template <typename T, class OpClass, class... Ranges>
-    inline const T& OperationMaster<T,OpClass,Ranges...>::get(size_t pos) const
-    {
-	//VCHECK(pos);
-	return mDataPtr[pos];
-	}*/
-    /*
-    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
-    {
-	mBlock.set( mIndex->pos() );
-	return *this;
-    }
-    */
    
     /****************************
      *   ConstOperationRoot     *
@@ -629,7 +444,6 @@ namespace MultiArrayTools
     ConstOperationRoot<T,Ranges...>::
     ConstOperationRoot(const MultiArrayBase<T,Ranges...>& ma,
 		       const std::shared_ptr<typename Ranges::IndexType>&... indices) :
-	//OperationTemplate<T,ConstOperationRoot<T,Ranges...> >(this),
 	mArrayRef(ma), mDataPtr(mArrayRef.data()),
 	mIndex( mkIndex(ma,indices...) ), mIInfo(*mIndex)
     {}
@@ -644,38 +458,14 @@ namespace MultiArrayTools
 	(*mIndex)(indices...);
 	return i;
     }
-    /*
+
-    template <typename T, class... Ranges>
-    const Block<T>& ConstOperationRoot<T,Ranges...>::get() const
-    {
-	block();
-	return mBlock;
-    }
-    */
     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>
-    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
-    {
-	mBlock.set( (*mIndex)().pos() );
-	return *this;
-    }
-    */
     template <typename T, class... Ranges>
     MExt<void> ConstOperationRoot<T,Ranges...>::rootSteps(std::intptr_t iPtrNum) const
     {
@@ -698,10 +488,8 @@ namespace MultiArrayTools
     OperationRoot<T,Ranges...>::
     OperationRoot(MutableMultiArrayBase<T,Ranges...>& ma,
 		  const std::shared_ptr<typename Ranges::IndexType>&... indices) :
-	//OperationTemplate<T,OperationRoot<T,Ranges...> >(this),
 	mArrayRef(ma), mDataPtr(mArrayRef.data()),
 	mIndex( mkIndex( ma, indices... ) ), mIInfo(*mIndex)
-	//mBlockII(nullptr)
     {}
 
     template <typename T, class... Ranges>
@@ -719,67 +507,16 @@ namespace MultiArrayTools
     template <class OpClass>
     OperationMaster<T,OpClass,Ranges...> OperationRoot<T,Ranges...>::operator=(const OpClass& in)
     {
-	//if(mBlockII != nullptr){
-	//    return OperationMaster<T,OpClass,Ranges...>(mArrayRef, in, mIndex, mBlockII);
-	//}
-	//else {
 	return OperationMaster<T,OpClass,Ranges...>(mArrayRef, in, mIndex);
-	    //}
-    }
-    /*
-    template <typename T, class... Ranges>
-    const MBlock<T>& OperationRoot<T,Ranges...>::get() const
-    {
-	block();
-	return mBlock;
     }
 
-    template <typename T, class... Ranges>
-    MBlock<T>& OperationRoot<T,Ranges...>::get()
-    {
-	block();
-	return mBlock;
-    }
-    
-    template <typename T, class... Ranges>
-    template <class ET>
-    inline const T& OperationRoot<T,Ranges...>::get(const ET& pos) const
-    {
-	return mDataPtr[pos.val()];
-    }
-    */
     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>
-    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 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
-    {
-	mBlock.set( (*mIndex)().pos() );
-	return *this;
-    }
-    */
     template <typename T, class... Ranges>
     MExt<void> OperationRoot<T,Ranges...>::rootSteps(std::intptr_t iPtrNum) const
     {
@@ -799,16 +536,8 @@ namespace MultiArrayTools
     
     template <typename T, class OpFunction, class... Ops>
     Operation<T,OpFunction,Ops...>::Operation(const Ops&... ops) :
-	//OperationTemplate<T,Operation<T,OpFunction,Ops...> >(this),
 	mOps(ops...) {}
-    /*
+
-    template <typename T, class OpFunction, class... Ops>
-    const BlockResult<T>& Operation<T,OpFunction,Ops...>::get() const
-    {
-	PackNum<sizeof...(Ops)-1>::template unpackArgs<T,OpFunction>(mRes, mOps);
-	return mRes;
-    }
-    */
     template <typename T, class OpFunction, class... Ops>
     template <class ET>
     inline T Operation<T,OpFunction,Ops...>::get(ET pos) const
@@ -817,26 +546,7 @@ namespace MultiArrayTools
 	return PackNum<sizeof...(Ops)-1>::
 	    template mkOpExpr<SIZE,T,ET,OpTuple,OpFunction>(pos, mOps);
     }
-    /*
-    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
-    {
-	//mBlock.set( mIndex->pos() );
-	return *this;
-    }
-    */
     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))
@@ -859,20 +569,9 @@ namespace MultiArrayTools
 
     template <typename T, class Op, class IndexType>
     Contraction<T,Op,IndexType>::Contraction(const Op& op, std::shared_ptr<IndexType> ind) :
-	//OperationTemplate<T,Contraction<T,Op,IndexType> >(this),
 	mOp(op),
 	mInd(ind) {}
-    /*
+
-    template <typename T, class Op, class IndexType>
-    const BlockResult<T>& Contraction<T,Op,IndexType>::get() const
-    {
-	BlockBinaryOpSelf<T,std::plus<T>,decltype(mOp.get())> f(mRes);
-	for(*mInd = 0; mInd->pos() != mInd->max(); ++(*mInd)){
-	    f(mOp.get());
-	}
-	return mRes;
-    }
-    */
     // forward loop !!!!
     template <typename T, class Op, class IndexType>
     template <class ET>
@@ -880,22 +579,7 @@ namespace MultiArrayTools
     {
 	return mOp.template get<ET>(pos);
     }
-    /*
-    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
-    {
-	return *this;
-    }
-    */
     template <typename T, class Op, class IndexType>
     auto Contraction<T,Op,IndexType>::rootSteps(std::intptr_t iPtrNum) const
 	-> decltype(mOp.rootSteps(iPtrNum))