some implementations reg FOR-utilities (step-array creation still missing)

src/multi_array_operation.h

@@ -50,7 +50,7 @@ namespace MultiArrayTools
 	
 	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,std::plus<T>,OperationClass,Second>;
@@ -84,7 +84,9 @@ namespace MultiArrayTools
 	class AssignmentExpr
 	{
 	public:
-	    AssignmentExpr(OperationMaster* mPtr, OpClass* secPtr);
+	    static size_t layer() { return 0; }
+	    
+	    AssignmentExpr(OperationMaster* mPtr, const OpClass* secPtr);
 
 	    AssignmentExpr(AssignmentExpr&& in) = default;
 	    AssignmentExpr& operator=(AssignmentExpr&& in) = default;
@@ -95,7 +97,7 @@ namespace MultiArrayTools
 	    AssignmentExpr() = default;
 	    
 	    OperationMaster* mMPtr;
-	    OpClass* mSecPtr;
+	    const OpClass* mSecPtr;
 	};
 	
 	typedef T value_type;
@@ -103,6 +105,8 @@ namespace MultiArrayTools
 	typedef ContainerRange<Ranges...> CRange;
 	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);
@@ -117,7 +121,7 @@ namespace MultiArrayTools
 	std::vector<BTSS> block(const IndexInfo* blockIndex, bool init = false) const;
 	const OperationMaster& block() const;
 	
-    protected:
+    private:
 
 	std::shared_ptr<IndexType> mkIndex(std::shared_ptr<typename CRange::IndexType>& index);
 	void performAssignment(std::intptr_t blockIndexNum);
@@ -141,6 +145,8 @@ namespace MultiArrayTools
 	typedef ContainerRange<Ranges...> CRange;
 	typedef typename CRange::IndexType IndexType;
 	typedef Block<T> bType;
+
+	static size_t rootNum() { return 1; }
 	
 	ConstOperationRoot(const MultiArrayBase<T,Ranges...>& ma,
 			   const std::shared_ptr<typename Ranges::IndexType>&... indices);
@@ -149,8 +155,10 @@ namespace MultiArrayTools
 
 	std::vector<BTSS> block(const IndexInfo* blockIndex, bool init = false) const;
 	const ConstOperationRoot& block() const;
+
+	std::tuple<size_t> rootSteps(const IndexInfo* ii = nullptr) const; // nullptr for simple usage with decltype
 	
-    protected:
+    private:
 
 	std::shared_ptr<IndexType>
 	mkIndex(const MultiArrayBase<T,Ranges...>& ma,
@@ -174,6 +182,8 @@ namespace MultiArrayTools
 	typedef typename CRange::IndexType IndexType;
 	typedef MBlock<T> bType;
 
+	static size_t rootNum() { return 1; }
+	
 	OperationRoot(MutableMultiArrayBase<T,Ranges...>& ma,
 		      const std::shared_ptr<typename Ranges::IndexType>&... indices);
 
@@ -186,8 +196,10 @@ namespace MultiArrayTools
 	OperationRoot& set(const IndexInfo* blockIndex);
 	std::vector<BTSS> block(const IndexInfo* blockIndex, bool init = false) const;
 	const OperationRoot& block() const;
+
+	std::tuple<size_t> rootSteps(const IndexInfo* ii = nullptr) const; // nullptr for simple usage with decltype
 	
-    protected:
+    private:
 
 	std::shared_ptr<IndexType>
 	mkIndex(const MultiArrayBase<T,Ranges...>& ma,
@@ -199,6 +211,18 @@ namespace MultiArrayTools
 	mutable bType mBlock;
 	const IndexInfo* mBlockII; // predefine to save time
     };
+
+    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 <typename T, class OpFunction, class... Ops>
     class Operation : public OperationTemplate<T,Operation<T,OpFunction,Ops...> >
@@ -210,17 +234,24 @@ namespace MultiArrayTools
 	typedef OperationTemplate<T,Operation<T,OpFunction,Ops...> > OT;
 	typedef OpFunction F;
 	typedef BlockResult<T> bType;
-	
+
+	static size_t rootNum() { return sumRootNum<Ops...>(); }
+
+    private:
+	std::tuple<Ops const&...> mOps;
+	mutable bType mRes;
+
+    public:
 	Operation(const Ops&... ops);
 	
 	const BlockResult<T>& get() const;
 
 	std::vector<BTSS> block(const IndexInfo* blockIndex, bool init = false) const;
 	const Operation& block() const;
+
+	auto rootSteps(const IndexInfo* ii = nullptr) const // nullptr for simple usage with decltype
+	    -> decltype(PackNum<sizeof...(Ops)-1>::mkStepTuple(ii, mOps));
 	
-    protected:
-	std::tuple<Ops const&...> mOps;
-	mutable bType mRes;
     };
     
     template <typename T, class Op, class IndexType>
@@ -231,19 +262,26 @@ 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(); }
+
+    private:
+
+	const Op& mOp;
+	std::shared_ptr<IndexType> mInd;
+	mutable bType mRes;
+
+    public:
 	Contraction(const Op& op, std::shared_ptr<IndexType> ind);
 	
 	const BlockResult<T>& get() const;
 
 	std::vector<BTSS> block(const IndexInfo* blockIndex, bool init = false) const;
 	const Contraction& block() const;
-	
-    protected:
 
-	const Op& mOp;
+	auto rootSteps(const IndexInfo* ii = nullptr) const  // nullptr for simple usage with decltype
-	std::shared_ptr<IndexType> mInd;
+	    -> decltype(mOp.rootSteps(ii));
-	mutable bType mRes;
+	
     };
     
 }
@@ -399,13 +437,21 @@ namespace MultiArrayTools
     template <typename T, class OpClass, class... Ranges>
     void OperationMaster<T,OpClass,Ranges...>::performAssignment(std::intptr_t blockIndexNum)
     {
-	//static auto loop = mkLoop(mIndex, *this, mSecond);
+#ifdef XX_USE_NEW_LOOP_ROUTINE_XX
-	//loop();
+	// === N E W ===
-	
+	AssignmentExpr ae(this, &mSecond);
+	typedef decltype(mSecond.rootSteps()) RootStepType;
+	std::array<RootStepType,/*sizeof single indices!!*/>
+	    ee(PackNum</*sizeof single indices!!*/-1>::mkExt(mIndex, mSecond));
+	static auto loop = mIndex->ifor(ee, ae);
+	loop();
+#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>
@@ -437,7 +483,7 @@ namespace MultiArrayTools
 	mBlock.set( mIndex->pos() );
 	return *this;
     }
-    
+
     /****************************
      *   ConstOperationRoot     *
      ****************************/
@@ -484,6 +530,12 @@ namespace MultiArrayTools
 	mBlock.set( (*mIndex)().pos() );
 	return *this;
     }
+
+    template <typename T, class... Ranges>
+    std::tuple<size_t> ConstOperationRoot<T,Ranges...>::rootSteps(const IndexInfo* ii) const
+    {
+	return std::tuple<size_t>(0ul); // !!!!!!
+    }
     
     /***********************
      *   OperationRoot     *
@@ -559,6 +611,12 @@ namespace MultiArrayTools
 	mBlock.set( (*mIndex)().pos() );
 	return *this;
     }
+
+    template <typename T, class... Ranges>
+    std::tuple<size_t> OperationRoot<T,Ranges...>::rootSteps(const IndexInfo* ii) const
+    {
+	return std::tuple<size_t>(0ul); // !!!!!!
+    }
     
     /*******************
      *   Operation     *
@@ -595,6 +653,14 @@ namespace MultiArrayTools
 	return *this;
     }
 
+    template <typename T, class OpFunction, class... Ops>
+    auto Operation<T,OpFunction,Ops...>::rootSteps(const IndexInfo* ii) const
+	-> decltype(PackNum<sizeof...(Ops)-1>::mkStepTuple(ii, mOps))
+    {
+	return PackNum<sizeof...(Ops)-1>::mkStepTuple(ii, mOps);
+    }
+
+    
     /*********************
      *   Contraction     *
      *********************/
@@ -629,7 +695,14 @@ namespace MultiArrayTools
     {
 	return *this;
     }
-    
+
+    template <typename T, class Op, class IndexType>
+    auto Contraction<T,Op,IndexType>::rootSteps(const IndexInfo* ii) const
+	-> decltype(mOp.rootSteps(ii))
+    {
+	return mOp.rootSteps(ii);
+    }
+
 }
 
 #endif

src/pack_num.h

@@ -44,11 +44,24 @@ namespace MultiArrayHelper
 	}
 
 	template <typename... T>
-	static void printTuple(std::ostream& out, const std::tuple<T...>& tp){
+	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 mkStepTuple(const IndexInfo* ii, std::tuple<Ops const&...> otp)
+	    -> decltype(std::tuple_cat( PackNum<N-1>::mkStepTuple(ii, otp), std::get<N>(otp).rootSteps(ii) ))
+	{
+	    return std::tuple_cat( PackNum<N-1>::mkStepTuple(ii, otp), std::get<N>(otp).rootSteps(ii) );
+	}
+	/*
+	template <class IndexClass, class OpClass>
+	static auto mkExt(const std::shared_ptr<IndexClass>& idxPtr, const OpClass& second)
+	{
+	// !!!!!
+	}*/
     };
     
     template<>
@@ -83,10 +96,18 @@ namespace MultiArrayHelper
 	}
 
 	template <typename... T>
-	static void printTuple(std::ostream& out, const std::tuple<T...>& tp){
+	static void printTuple(std::ostream& out, const std::tuple<T...>& tp)
+	{	    
 	    out << std::get<sizeof...(T)-1>(tp);
 	}
-	
+
+	template <class... Ops>
+	static auto mkStepTuple(const IndexInfo* ii, std::tuple<Ops const&...> otp)
+	    -> decltype(std::get<0>(otp).rootSteps(ii))
+	{
+	    return std::get<0>(otp).rootSteps(ii);
+	}
+
     };
 
 

src/ranges/container_range.h

@@ -87,10 +87,10 @@ namespace MultiArrayTools
 	
 	std::string id();
 	void print(size_t offset);
-
+	
-	template <class... Exprs>
+	template <class Ext, class Exprs>
-	auto ifor(Exprs&&... exs) const
+	auto ifor(const Ext& ext, Exprs&& exs) const
-	    -> decltype(RPackNum<sizeof...(Indices)-1>::mkFor(mIPack, exs...));
+	    -> decltype(RPackNum<sizeof...(Indices)-1>::mkFor(mIPack, ext, exs));
 	
     };
 
@@ -371,11 +371,11 @@ namespace MultiArrayTools
     }
 
     template <class... Indices>
-    template <class... Exprs>
+    template <class Ext, class Exprs>
-    auto ContainerIndex<Indices...>::ifor(Exprs&&... exs) const
+    auto ContainerIndex<Indices...>::ifor(const Ext& ext, Exprs&& exs) const
-	-> decltype(RPackNum<sizeof...(Indices)-1>::mkFor(mIPack, exs...))
+	-> decltype(RPackNum<sizeof...(Indices)-1>::mkFor(mIPack, ext, exs))
     {
-	return RPackNum<sizeof...(Indices)-1>::mkFor(mIPack, exs...);
+	return RPackNum<sizeof...(Indices)-1>::mkFor(mIPack, ext, exs);
     }
 
     

src/ranges/multi_range.h

@@ -104,9 +104,9 @@ namespace MultiArrayTools
 	std::string id();
 	void print(size_t offset);
 
-	template <class... Exprs>
+	template <class Ext, class Exprs>
-	auto ifor(Exprs&&... exs) const
+	auto ifor(const Ext& ext, Exprs&& exs) const
-	    -> decltype(RPackNum<sizeof...(Indices)-1>::mkFor(mIPack, exs...));
+	    -> decltype(RPackNum<sizeof...(Indices)-1>::mkFor(mIPack, ext, exs));
     };
 
     /*************************
@@ -415,11 +415,11 @@ namespace MultiArrayTools
     }
 
     template <class... Indices>
-    template <class... Exprs>
+    template <class Ext, class Exprs>
-    auto MultiIndex<Indices...>::ifor(Exprs&&... exs) const
+    auto MultiIndex<Indices...>::ifor(const Ext& ext, Exprs&& exs) const
-	-> decltype(RPackNum<sizeof...(Indices)-1>::mkFor(mIPack, exs...))
+	-> decltype(RPackNum<sizeof...(Indices)-1>::mkFor(mIPack, ext, exs))
     {
-	return RPackNum<sizeof...(Indices)-1>::mkFor(mIPack, exs...);
+	return RPackNum<sizeof...(Indices)-1>::mkFor(mIPack, ext, exs);
     }
     
     /*************************

src/ranges/rpack_num.h

@@ -219,14 +219,13 @@ namespace MultiArrayHelper
 	    RPackNum<N-1>::buildInfoVec(out, ip, bs);
 	}
 
-	template <class IndexPack, class... Exprs>
+	template <class IndexPack, class Ext, class Exprs>
-	static auto mkFor(const IndexPack& ipack, Exprs&&... exs)
+	static auto mkFor(const IndexPack& ipack, const Ext& ext, Exprs&& exs)
 	    -> decltype(std::get<std::tuple_size<IndexPack>::value-N>(ipack)
-			->ifor(RPackNum<N-1>::template mkFor<IndexPack,Exprs...>
+			->ifor(ext, RPackNum<N-1>::mkFor(ipack, ext, exs) ) )
-			       (ipack, exs...) ) )
 	{
 	    return std::get<std::tuple_size<IndexPack>::value-N>(ipack)
-		->ifor( RPackNum<N-1>::template mkFor<IndexPack,Exprs...>(ipack, exs...) );
+		->ifor( ext, RPackNum<N-1>::mkFor(ipack, ext, exs) );
 	}
 	
     };
@@ -392,12 +391,12 @@ namespace MultiArrayHelper
 	    out.emplace_back(*std::get<POS>(ip), std::get<POS>(bs));
 	}
 
-	template <class IndexPack, class... Exprs>
+	template <class IndexPack, class Ext, class Exprs>
-	static auto mkFor(const IndexPack& ipack, Exprs&&... exs)
+	static auto mkFor(const IndexPack& ipack, const Ext& ext, Exprs&& exs)
-	    -> decltype(std::get<std::tuple_size<IndexPack>::value-1>(ipack)
+	    -> decltype(std::get<std::tuple_size<IndexPack>::value-1>(ipack, ext)
-			->ifor(exs...) )
+			->ifor(ext, exs) )
 	{
-	    return std::get<std::tuple_size<IndexPack>::value-1>(ipack)->ifor(exs...);
+	    return std::get<std::tuple_size<IndexPack>::value-1>(ipack)->ifor(ext, exs);
 	}
     };
 	

src/ranges/single_range.h

@@ -64,8 +64,8 @@ namespace MultiArrayTools
 	std::string id();
 	void print(size_t offset);
 
-	template <class Expr, class... Ops>
+	template <class Expr, size_t DIM, class... Ops>
-	auto ifor(std::tuple<to_size_t<Ops>...>&& ee, Expr&& ex) const
+	auto ifor(const std::array<std::tuple<to_size_t<Ops>...>,DIM>& ee, Expr&& ex) const
 	    -> For<SingleIndex<U,TYPE>,Expr,Ops...>;
     };
 
@@ -248,12 +248,12 @@ namespace MultiArrayTools
     }
 
     template <typename U, SpaceType TYPE>
-    template <class Expr, class... Ops>
+    template <class Expr, size_t DIM, class... Ops>
-    auto SingleIndex<U,TYPE>::ifor(std::tuple<to_size_t<Ops>...>&& ee, Expr&& ex) const
+    auto SingleIndex<U,TYPE>::ifor(const std::array<std::tuple<to_size_t<Ops>...>,DIM>& ee, Expr&& ex) const
 	-> For<SingleIndex<U,TYPE>,Expr,Ops...>
     {
-	return For<SingleIndex<U,TYPE>,Expr>
+	static const size_t layer = typename Expr::layer(); 
-	    ( std::make_shared<SingleIndex<U,TYPE> >(*this), ex, ee );
+	return For<SingleIndex<U,TYPE>,Expr,Ops...>( this, ex, std::get<layer>( ee ) );
     }
 
     

src/xfor/xfor.h

@@ -15,6 +15,8 @@ namespace MultiArrayHelper
     {
     public:
 
+	static size_t layer() { return typename Expr::layer() + 1; }
+	
 	For(For&& in) = default;
 	For& operator=(For&& in) = default;
 	
@@ -24,18 +26,32 @@ namespace MultiArrayHelper
 
 	For(const std::shared_ptr<IndexClass>& indPtr,
 	    Expr&& expr, std::tuple<to_size_t<Ops>...>&& ext);
+
+	template <typename... Args>
+	For(IndexClass* indPtr,
+	    std::tuple<to_size_t<Ops>...>&& ext, const Args&... args);
+
+	For(IndexClass* indPtr,
+	    Expr&& expr, std::tuple<to_size_t<Ops>...>&& ext);
+
 	
 	inline void operator()(size_t mlast, const std::tuple<to_size_t<Ops>...>& last) const;
+	inline void operator()(size_t mlast = 0) const;
 	
     private:
 	For() = default;
 	
-	mutable std::shared_ptr<IndexClass> mIndPtr;
+	IndexClass* mIndPtr;
 	const Expr mExpr;
 	const std::tuple<to_size_t<Ops>...> mExt;
     };
 
-    
+    template <size_t N>
+    size_t exceptMax(size_t max) { return max; }
+
+    template <>
+    size_t exceptMax<1>(size_t max) { return 1; }
+
     
 } // namespace MultiArrayHelper
 
@@ -51,11 +67,23 @@ namespace MultiArrayHelper
     For<IndexClass,Expr,Ops...>::For(const std::shared_ptr<IndexClass>& indPtr,
 				     std::tuple<to_size_t<Ops>...>&& ext,
 				     const Args&... args) :
-	mIndPtr(indPtr), mExpr(args...), mExt(ext) {}
+	mIndPtr(indPtr.get()), mExpr(args...), mExt(ext) {}
 
     template <class IndexClass, class Expr, class... Ops>
     For<IndexClass,Expr,Ops...>::For(const std::shared_ptr<IndexClass>& indPtr,
 				     Expr&& expr, std::tuple<to_size_t<Ops>...>&& ext) :
+	mIndPtr(indPtr.get()), mExpr(expr), mExt(ext) {}
+
+    template <class IndexClass, class Expr, class... Ops>
+    template <typename... Args>
+    For<IndexClass,Expr,Ops...>::For(IndexClass* indPtr,
+				     std::tuple<to_size_t<Ops>...>&& ext,
+				     const Args&... args) :
+	mIndPtr(indPtr), mExpr(args...), mExt(ext) {}
+
+    template <class IndexClass, class Expr, class... Ops>
+    For<IndexClass,Expr,Ops...>::For(IndexClass* indPtr,
+				     Expr&& expr, std::tuple<to_size_t<Ops>...>&& ext) :
 	mIndPtr(indPtr), mExpr(expr), mExt(ext) {}
     
     template <class IndexClass, class Expr, class... Ops>
@@ -64,7 +92,21 @@ namespace MultiArrayHelper
     {
 	static const size_t opNum = sizeof...(Ops);
 	auto& ind = *mIndPtr;
-	const size_t max = ind.max();
+	const size_t max = exceptMax<For<IndexClass,Expr,Ops...>::layer()>( ind.max() ); // blocking
+	for(ind = 0; ind.pos() != max; ++ind){
+	    const size_t mnpos = mlast * max + ind.pos();
+	    const std::tuple<to_size_t<Ops>...> npos = std::move( XFPackNum<opNum-1>::mkPos(ind, mExt, last) );
+	    mExpr(mnpos, npos);
+	}
+    }
+
+    template <class IndexClass, class Expr, class... Ops>
+    inline void For<IndexClass,Expr,Ops...>::operator()(size_t mlast) const
+    {
+	static const size_t opNum = sizeof...(Ops);
+	std::tuple<to_size_t<Ops>...> last(to_size_t<Ops>(0)...);
+	auto& ind = *mIndPtr;
+	const size_t max = exceptMax<For<IndexClass,Expr,Ops...>::layer()>( ind.max() ); // blocking
 	for(ind = 0; ind.pos() != max; ++ind){
 	    const size_t mnpos = mlast * max + ind.pos();
 	    const std::tuple<to_size_t<Ops>...> npos = std::move( XFPackNum<opNum-1>::mkPos(ind, mExt, last) );