From e2405738c188ad119e317a61ba5900f51face073 Mon Sep 17 00:00:00 2001
From: Christian Zimmermann <christian.zimmermann@ur.de>
Date: Sun, 23 Oct 2022 18:29:07 +0200
Subject: [PATCH] derive COpInterface : XprInterface; allow return values for
 Xpr (VXpr+DXpr -> templates)

---
 src/include/operation/op_types.cc.h | 164 ++++++++++------------------
 src/include/operation/op_types.h    |  53 ++++-----
 src/include/operation/op_xpr.cc.h   |  49 ---------
 src/include/operation/op_xpr.h      |  37 -------
 src/include/ranges/crange.cc.h      |   2 +-
 src/include/ranges/dindex.h         |   2 +-
 src/include/ranges/xindex.cc.h      |   4 +-
 src/include/ranges/xindex.h         |   4 +-
 src/include/ranges/yrange.h         |   2 +-
 src/include/xpr/for.cc.h            | 139 +++++++++++++++--------
 src/include/xpr/for.h               |  57 ++++++----
 src/include/xpr/for_type.h          |  23 ----
 src/include/xpr/xpr_base.cc.h       |  42 +++----
 src/include/xpr/xpr_base.h          |  32 +++---
 src/lib/ranges/dindex.cc            |   4 +-
 src/lib/ranges/yrange.cc            |   4 +-
 16 files changed, 262 insertions(+), 356 deletions(-)
 delete mode 100644 src/include/operation/op_xpr.cc.h
 delete mode 100644 src/include/operation/op_xpr.h
 delete mode 100644 src/include/xpr/for_type.h

diff --git a/src/include/operation/op_types.cc.h b/src/include/operation/op_types.cc.h
index 4df3517..6b327ff 100644
--- a/src/include/operation/op_types.cc.h
+++ b/src/include/operation/op_types.cc.h
@@ -4,10 +4,8 @@
 
 #include "op_types.h"
 #include "xpr/pos_type.h"
-#include "ranges/range_helper.h"
 #include "xpr/op_xpr.h"
-#include "statics/static_for.h"
-
+#include "op_utility.h"
 
 namespace CNORXZ
 {
@@ -17,20 +15,29 @@ namespace CNORXZ
      **********************/
 
     template <class OpT>
-    template <class IndexT>
-    auto COpInterface<OpT>::c(const std::shared_ptr<IndexT>& ind) const
+    template <class F, class IndexT>
+    decltype(auto) COpInterface<OpT>::c(F&& f, const Sptr<IndexT>& ind) const
     {
-	return Contraction<T,OpT,IndexT>(THIS(), ind);
+	return mkContraction(std::forward<F>(f), THIS(), ind);
+    }
+
+    template <class OpT>
+    template <class IndexT>
+    decltype(auto) COpInterface<OpT>::c(const Sptr<IndexT>& ind) const
+    {
+	return mkContraction([](auto& a, const auto& b) { a += b; },
+			     THIS(), ind);
     }
 
     template <class OpT>
     template <class F, class... Args>
     constexpr auto a(F&& f, Args&&... args) const;
-    auto COpInterface<OpT>::o(F&& f, Args&&... args) const
+    decltype(auto) COpInterface<OpT>::o(F&& f, Args&&... args) const
     {
-        return Operation<R,F,OpT,Args...>(f, THIS(), args...);
+        return Operation<F,OpT,Args...>(std::forward<F>(f), THIS(), args...);
     }
 
+    
     /*********************
      *   OpInterface     *
      *********************/
@@ -39,7 +46,7 @@ namespace CNORXZ
     template <class IndexT, class F, class... Args>
     constexpr decltype(auto) OpInterface<OpT>::ax(const Sptr<IndexT>& ind, F&& f, Args&&... args)
     {
-	return ind->ifor( SPos<1>(), OpXpr<F,OpT,Args...>(f, THIS(), args...) );
+	return ind->ifor( Operation<F,OpT,Args...>(f, THIS(), args...) );
     }
 
     template <class OpT>
@@ -49,7 +56,7 @@ namespace CNORXZ
 	return ax(ind, f, args...)();
     }
 
-
+    
     /***************
      *   COpRoot   *
      ***************/
@@ -68,7 +75,7 @@ namespace CNORXZ
 
     template <typename T, class IndexT>
     template <class PosT>
-    constexpr decltype(auto) COpRoot<T,IndexT>::get(const PosT& pos) const
+    constexpr decltype(auto) COpRoot<T,IndexT>::operator()(const PosT& pos) const
     {
 	if constexpr(is_epos_type<PosT>::value){
 	    return vreg(mData,pos); // distinguish between consecutive/non-consecutive
@@ -78,6 +85,12 @@ namespace CNORXZ
 	}
     }
 
+    template <typename T, class IndexT>
+    constexpr decltype(auto) COpRoot<T,IndexT>::operator()() const
+    {
+	return mData[0];
+    }
+    
     template <typename T, class IndexT>
     template <SizeT I>
     constexpr decltype(auto) COpRoot<T,IndexT>::rootSteps(const IndexId<I>& id) const
@@ -85,13 +98,6 @@ namespace CNORXZ
 	return mIndex->stepSize(id);
     }
 
-    template <typename T, class IndexT>
-    template <class Expr>
-    constexpr decltype(auto) COpRoot<T,IndexT>::loop(Xpr&& xpr) const
-    {
-	return xpr;
-    }
-
     template <typename T, class IndexT>
     const T* COpRoot<T,IndexT>::data() const
     {
@@ -140,7 +146,7 @@ namespace CNORXZ
 
     template <typename T, class IndexT>
     template <class PosT>
-    constexpr decltype(auto) OpRoot<T,IndexT>::get(const PosT& pos) const
+    constexpr decltype(auto) OpRoot<T,IndexT>::operator()(const PosT& pos) const
     {
 	if constexpr(is_epos_type<PosT>::value){
 	    return vreg(mData,pos); // distinguish between consecutive/non-consecutive
@@ -150,6 +156,12 @@ namespace CNORXZ
 	}
     }
 
+    template <typename T, class IndexT>
+    constexpr decltype(auto) OpRoot<T,IndexT>::operator()() const
+    {
+	return mData[0];
+    }
+    
     template <typename T, class IndexT>
     template <SizeT I>
     constexpr decltype(auto) OpRoot<T,IndexT>::rootSteps(const IndexId<I>& id) const
@@ -157,100 +169,41 @@ namespace CNORXZ
 	return mIndex->stepSize(id);
     }
 
-    template <typename T, class IndexT>
-    template <class Expr>
-    constexpr decltype(auto) OpRoot<T,IndexT>::loop(Xpr&& xpr) const;
-    {
-	return xpr;
-    }
-
 
     /*******************
      *   Operation     *
      *******************/
     
-    template <typename T, class OpFunction, class... Ops>
-    Operation<T,OpFunction,Ops...>::Operation(const Ops&... ops) :
-	mOps(ops...)
-    {
-	static_assert( FISSTATIC, "need function instance for non-static function" );
-    }
-
-    template <typename T, class OpFunction, class... Ops>
-    Operation<T,OpFunction,Ops...>::Operation(std::shared_ptr<OpFunction> ff,
-					      const Ops&... ops) :
+    template <class F, class... Ops>
+    Operation<F,Ops...>::Operation(F&& f, Ops&&... ops) :
 	mOps(ops...),
 	mF(ff)
+    {}
+
+    template <class F, class... Ops>
+    template <class PosT>
+    inline decltype(auto) Operation<F,Ops...>::operator()(const PosT& pos) const
     {
-	static_assert( not FISSTATIC, "using instance of static function" );
+	return pos_unpack_args(mF, pos, mOps);
     }
 
-    template <size_t I, class OpFunction, class ETuple, class OpTuple, typename... Args>
-    inline auto
-    mkOpExpr(std::shared_ptr<OpFunction> f, const ETuple& pos, const OpTuple& ops, Args... args)
+    template <class F, class... Ops>
+    template <class PosT>
+    inline decltype(auto) Operation<F,Ops...>::operator()() const
     {
-	if constexpr(I == std::tuple_size<OpTuple>{}){
-	    if constexpr(OpFunction::FISSTATIC){
-		return OpFunction::apply(args...);
-	    }
-	    else {
-		(*f)(args...);
-	    }
-	}
-	else {
-	    typedef typename std::remove_reference<decltype(std::get<I>(ops))>::type NextOpType;
-	    return mkOpExpr<I+1>
-		( f, getX<NextOpType::SIZE>(pos), ops, args..., std::get<I>(ops).get(pos));
-	}
+	return exec(std::make_index_sequence<sizeof...(Ops)>{});
     }
 
-    template <typename T, class OpFunction, class... Ops>
-    template <class ET>
-    inline auto Operation<T,OpFunction,Ops...>::get(ET pos) const
+    template <class F, class... Ops>
+    constexpr decltype(auto) Operation<F,Ops...>::rootSteps(const IndexId<I>& id) const
     {
-	return mkOpExpr<0>(mF, pos, mOps.mOps);
+	return mOps.rootSteps(id);
     }
 
-    template <typename T, class OpFunction, class... Ops>
-    template <typename V, class ET>
-    inline auto Operation<T,OpFunction,Ops...>::vget(ET pos) const
+    template <class F, class... Ops>
+    constexpr decltype(auto) Operation<F,Ops...>::exec(std::index_sequence<Is...> is) const
     {
-	return mkVOpExpr<0,V>(mF, pos, mOps.mOps);
-    }
-
-    template <size_t I, class OpTuple, class ETuple>
-    static inline void setOpPos(OpTuple& ot, const ETuple& et)
-    {
-	if constexpr(I != std::tuple_size<OpTuple>{}){
-	    typedef typename std::remove_reference<decltype(std::get<I>(ot))>::type NextOpType;
-	    std::get<I>( ot ).set( et );
-	    setOpPos<I+1>(ot, getX<NextOpType::SIZE>(et));
-	}
-    }
-    
-    template <typename T, class OpFunction, class... Ops>
-    template <class ET>
-    inline Operation<T,OpFunction,Ops...>& Operation<T,OpFunction,Ops...>::set(ET pos)
-    {
-	setOpPos<0>(mOps.mOps,pos);
-	return *this;
-    }
-
-    template <typename T, class OpFunction, class... Ops>
-    auto Operation<T,OpFunction,Ops...>::rootSteps(std::intptr_t iPtrNum) const
-	-> ExtType
-    {
-	return mOps.rootSteps(iPtrNum);
-    }
-
-    template <typename T, class OpFunction, class... Ops>
-    template <class Expr>
-    auto Operation<T,OpFunction,Ops...>::loop(Expr exp) const
-    {
-	return sfor_m<sizeof...(Ops),0>
-	    ( [&](auto i){ return std::get<i>(mOps.mOps); },
-	      [&](auto f, auto next) { return f.loop(next); },
-	      exp );
+	return mF( std::get<Is>(mOps)() ... );
     }
 
     
@@ -265,11 +218,17 @@ namespace CNORXZ
 
     template <class CXpr>
     template <class PosT>
-    constexpr decltype(auto) Contraction<CXpr>::get(const PosT& pos) const
+    constexpr decltype(auto) Contraction<CXpr>::operator()(const PosT& pos) const
     {
 	return mCXpr(pos);
     }
 
+    template <class CXpr>
+    constexpr decltype(auto) Contraction<CXpr>::operator()() const
+    {
+	return mCXpr();
+    }
+
     template <class CXpr>
     template <SizeT I>
     constexpr decltype(auto) Contraction<CXpr>::rootSteps(const IndexId<I>& id) const;
@@ -277,18 +236,11 @@ namespace CNORXZ
 	return mCXpr.stepSize(id);
     }
 
-    template <class CXpr>
-    template <class Expr>
-    constexpr decltype(auto) Contraction<CXpr>::loop(Xpr&& xpr) const
-    {
-	return xpr;
-    }
-
     template <class F, class Op, class IndexT>
     constexpr decltype(auto) mkContracion(F&& f, Op&& op, const Sptr<IndexT>& i)
     {
-	typedef decltype(i->ifor( mkOpXpr( f, op ) )) CXprT;
-	return Contraction<CXprT>( i->ifor( mkOpXpr( f, op ) ) );
+	typedef decltype(i->ifor( op, f )) CXprT; // TODO: implement ifor with func arg!!!
+	return Contraction<CXprT>( i->ifor( op, f ) );
     }
 }
 
diff --git a/src/include/operation/op_types.h b/src/include/operation/op_types.h
index 9b25aab..e5ef07c 100644
--- a/src/include/operation/op_types.h
+++ b/src/include/operation/op_types.h
@@ -10,30 +10,22 @@ namespace CNORXZ
 
 
     template <class OpT>
-    class COpInterface
+    class COpInterface : public XprInterface<OpT>
     {
     public:
 	
 	OpT& THIS() { return static_cast<OpT&>(*this); }
 	const OpT& THIS() const { return static_cast<const OpT&>(*this); }
 
+	template <class F, class IndexT>
+	constexpr decltype(auto) c(F&& f, const Sptr<IndexT>& ind) const;
+
 	template <class IndexT>
-	constexpr auto c(const Sptr<IndexT>& ind) const;
+	constexpr decltype(auto) c(const Sptr<IndexT>& ind) const;
 
 	template <class F, class... Args>
-	constexpr auto o(F&& f, Args&&... args) const;
+	constexpr decltype(auto) o(F&& f, Args&&... args) const;
 
-	template <class PosT>
-	constexpr decltype(auto) get(const PosT& pos) const
-	{ return THIS().get(pos); }
-
-	template <SizeT I>
-	constexpr decltype(auto) rootSteps(const IndexId<I>& id) const
-	{ return THIS().rootSteps(id); }
-	
-	template <class Xpr>
-	constexpr decltype(auto) loop(Xpr&& exp) const
-	{ return THIS().loop(exp); }
     };
 
     template <class OpT>
@@ -72,14 +64,13 @@ namespace CNORXZ
 	constexpr COpRoot(const T* data, const Sptr<IndexT>& ind);
 
 	template <class PosT>
-	constexpr decltype(auto) get(const PosT& pos) const;
+	constexpr decltype(auto) operator()(const PosT& pos) const;
+
+	constexpr decltype(auto) operator()() const;
 
 	template <SizeT I>
 	constexpr decltype(auto) rootSteps(const IndexId<I>& id) const;
 	
-	template <class Xpr>
-	constexpr decltype(auto) loop(Xpr&& xpr) const;
-
 	const T* data() const;
 	
     private:
@@ -114,14 +105,13 @@ namespace CNORXZ
         constexpr OpRoot& operator=(const OpRoot& in);
         
 	template <class PosT>
-	constexpr decltype(auto) get(const PosT& pos) const;
+	constexpr decltype(auto) operator()(const PosT& pos) const;
+
+	constexpr decltype(auto) operator()() const;
 
 	template <SizeT I>
 	constexpr decltype(auto) rootSteps(const IndexId<I>& id) const;
 
-	template <class Xpr>
-	constexpr decltype(auto) loop(Xpr&& exp) const;
-        
 	T* data() const;
 
     private:
@@ -140,15 +130,17 @@ namespace CNORXZ
 	constexpr Operation(F&& f, Ops&&... ops);
 
 	template <class PosT>
-	constexpr decltype(auto) get(const PosT& pos) const;
+	constexpr decltype(auto) operator()(const PosT& pos) const;
+
+	constexpr decltype(auto) operator()() const;
 
 	template <SizeT I>
-	constexpr decltype(auto) rootSteps(const IndexId<I>& id) const
-
-	template <class Xpr>
-	constexpr decltype(auto) loop(Xpr&& xpr) const;
+	constexpr decltype(auto) rootSteps(const IndexId<I>& id) const;
 
     private:
+	template <SizeT... Is>
+	constexpr decltype(auto) exec(std::index_sequence<Is...> is) const;
+     
 	Tuple<Ops...> mOps;
 	F mF;
 	
@@ -163,14 +155,13 @@ namespace CNORXZ
 	constexpr Contraction(CXpr&& cxpr);
 
 	template <class PosT>
-	constexpr decltype(auto) get(const PosT& pos) const;
+	constexpr decltype(auto) operator()(const PosT& pos) const;
+
+	constexpr decltype(auto) operator()() const;
 
 	template <SizeT I>
 	constexpr decltype(auto) rootSteps(const IndexId<I>& id) const;
 
-	template <class Xpr>
-	constexpr decltype(auto) loop(Xpr&& xpr) const;
-
     private:
 	CXpr mCXpr;
 	Sptr<IndexType> mInd;
diff --git a/src/include/operation/op_xpr.cc.h b/src/include/operation/op_xpr.cc.h
deleted file mode 100644
index 8488766..0000000
--- a/src/include/operation/op_xpr.cc.h
+++ /dev/null
@@ -1,49 +0,0 @@
-
-#ifndef __cxz_op_expr_cc_h__
-#define __cxz_op_expr_cc_h__
-
-#include "op_expr.h"
-#include "op_utility.h"
-
-namespace CNORXZ
-{
-    template <class F, class... Ops>
-    constexpr OpXpr<F,Ops...>::OpXpr(F&& f, Ops&&... ops) :
-	mF(std::forward<F>(f)),
-	mOps(ops...)
-    {
-	static_assert((is_operation<Ops>::value and ...), "got non-op type");
-    }
-	
-    template <class F, class... Ops>
-    template <class PosT1, class PosT2>
-    inline SizeT OpXpr<F,Ops...>::operator()(const PosT& last) const
-    {
-	pos_unpack_args(mF,last,mOps); // utility function (to be implemented)
-	// depending on whether Ops[N] is static or not call statically or dynamically .next()
-	return 0;
-    }
-
-    template <class F, class... Ops>
-    inline SizeT OpXpr<F,Ops...>::operator()() const
-    {
-	return 0;
-    }
-
-    template <class F, class... Ops>
-    template <SizeT I>
-    constexpr decltype(auto) OpXpr<F,Ops...>::rootSteps(const IndexId<I>& id) const
-    {
-	return rootStepsI(id, std::make_index_sequence<sizeof...(Ops)>{});
-    }
-
-    template <class F, class... Ops>
-    template <SizeT I, SizeT... Js>
-    constexpr decltype(auto) OpXpr<F,Ops...>::rootStepsI(const IndexId<I>& id, std::index_sequence<Js...> is) const
-    {
-	return (std::get<Js>(mOps).rootSteps(id) << ...);
-	// TODO: implement a << b which is a.extend(b)!!  
-    }
-}
-
-#endif
diff --git a/src/include/operation/op_xpr.h b/src/include/operation/op_xpr.h
deleted file mode 100644
index 727f283..0000000
--- a/src/include/operation/op_xpr.h
+++ /dev/null
@@ -1,37 +0,0 @@
-
-#ifndef __cxz_op_expr_h__
-#define __cxz_op_expr_h__
-
-#include "base/base.h"
-#include "operation/operation_base.h"
-#include "xpr/xpr_base.h"
-
-namespace CNORXZ
-{
-
-    template <class F, class... Ops>
-    class OpXpr
-    {
-    private:
-	F mF;
-	Tuple<Ops...> mOps;
-
-	template <SizeT I, SizeT... Js>
-	constexpr decltype(auto) rootStepsI(const IndexId<I>& id, std::index_sequence<Js...> is) const;
- 	
-    public:
-	DEFAULT_MEMBERS(OpXpr);
-
-	constexpr OpXpr(F&& f, Ops&&... ops);
-
-	template <class PosT>
-	inline SizeT operator()(const PosT& last) const;
-
-	inline SizeT operator()() const;
-
-	template <SizeT I>
-	constexpr decltype(auto) rootSteps(const IndexId<I>& id) const;
-    };
-}
-
-#endif
diff --git a/src/include/ranges/crange.cc.h b/src/include/ranges/crange.cc.h
index e07ad5f..6aa0e12 100644
--- a/src/include/ranges/crange.cc.h
+++ b/src/include/ranges/crange.cc.h
@@ -9,7 +9,7 @@ namespace CNORXZ
     template <class Xpr>
     decltype(auto) CIndex::ifor(const Xpr& xpr) const
     {
-	return For<0,Xpr>(this->max(), this->id(), xpr);
+	return For<0,Xpr>(this->max(), this->id(), xpr, NoF());
     }
 }
 
diff --git a/src/include/ranges/dindex.h b/src/include/ranges/dindex.h
index e4231d2..9d4817a 100644
--- a/src/include/ranges/dindex.h
+++ b/src/include/ranges/dindex.h
@@ -50,7 +50,7 @@ namespace CNORXZ
 	DType meta() const;
 	DIndex& at(const DType& meta);
 
-	DXpr ifor(const DXpr& xpr) const;
+	DXpr<SizeT> ifor(const DXpr<SizeT>& xpr) const;
 
     private:
 	XIndexPtr mI;
diff --git a/src/include/ranges/xindex.cc.h b/src/include/ranges/xindex.cc.h
index b445046..655d6c0 100644
--- a/src/include/ranges/xindex.cc.h
+++ b/src/include/ranges/xindex.cc.h
@@ -155,9 +155,9 @@ namespace CNORXZ
     }
 
     template <class Index, typename Meta>
-    DXpr XIndex<Index,Meta>::ifor(const DXpr& xpr) const
+    DXpr<SizeT> XIndex<Index,Meta>::ifor(const DXpr<SizeT>& xpr) const
     {
-	return DXpr(mI->ifor(xpr));
+	return DXpr<SizeT>(mI->ifor(xpr));
     }
 
 }
diff --git a/src/include/ranges/xindex.h b/src/include/ranges/xindex.h
index 38b359b..1fe2a89 100644
--- a/src/include/ranges/xindex.h
+++ b/src/include/ranges/xindex.h
@@ -42,7 +42,7 @@ namespace CNORXZ
 	virtual DType meta() const = 0;
 	virtual XIndexBase& at(const DType& meta) = 0;
 
-	virtual DXpr ifor(const DXpr& xpr) const = 0;
+	virtual DXpr<SizeT> ifor(const DXpr<SizeT>& xpr) const = 0;
     };
 
     //Sptr<XIndexBase>& operator++(Sptr<XIndexBase>& i);
@@ -91,7 +91,7 @@ namespace CNORXZ
 	virtual DType meta() const override final;
 	virtual XIndexBase& at(const DType& meta) override final;
 
-	virtual DXpr ifor(const DXpr& xpr) const override final;
+	virtual DXpr<SizeT> ifor(const DXpr<SizeT>& xpr) const override final;
 
     private:
 	IndexPtr<Index,Meta> mI;
diff --git a/src/include/ranges/yrange.h b/src/include/ranges/yrange.h
index 58ae0a8..95f9b1d 100644
--- a/src/include/ranges/yrange.h
+++ b/src/include/ranges/yrange.h
@@ -51,7 +51,7 @@ namespace CNORXZ
 	DType meta() const;
 	YIndex& at(const DType& meta);
 
-	DXpr ifor(const DXpr& xpr) const;
+	DXpr<SizeT> ifor(const DXpr<SizeT>& xpr) const;
 
     private:
 
diff --git a/src/include/xpr/for.cc.h b/src/include/xpr/for.cc.h
index 58ec16e..14b3948 100644
--- a/src/include/xpr/for.cc.h
+++ b/src/include/xpr/for.cc.h
@@ -8,42 +8,63 @@
 
 namespace CNORXZ
 {
+    /*************
+     *   ZeroF   *
+     *************/
+
+    template <typename... T>
+    constexpr decltype(auto) ZeroF::operator()(const T&... as) const
+    {
+	return 0;
+    }
+
+    template <typename... T>
+    constexpr decltype(auto) NoF::operator()(const T&... as) const
+    {
+	return;
+    }
+
     /***********
      *   For   *
      ***********/
     
-    template <SizeT L, class Xpr>
-    constexpr For<L,Xpr>::For(SizeT size, const IndexId<L>& id, const Xpr& xpr) :
+    template <SizeT L, class Xpr, class F>
+    constexpr For<L,Xpr,F>::For(SizeT size, const IndexId<L>& id, const Xpr& xpr, F&& f) :
 	mSize(size),
 	mId(id),
 	mXpr(xpr),
-	mExt(mXpr.rootSteps(mId))
+	mExt(mXpr.rootSteps(mId)),
+	mF(f)
     {}
 
-    template <SizeT L, class Xpr>
+    template <SizeT L, class Xpr, class F>
     template <class PosT>
-    inline SizeT For<L,Xpr>::operator()(const PosT& last) const
+    inline decltype(auto) For<L,Xpr,F>::operator()(const PosT& last) const
     {
+	typedef typename std::remove_reference<decltype(mXpr(last + mExt * UPos(0)))>::type OutT;
+	auto o = OutT();
 	for(SizeT i = 0; i != mSize; ++i){
 	    const auto pos = last + mExt * UPos(i);
-	    mXpr(pos);
+	    mF(o, mXpr(pos));
 	}
-	return 0;
+	return o;
     }
 	
-    template <SizeT L, class Xpr>
-    inline SizeT For<L,Xpr>::operator()() const
+    template <SizeT L, class Xpr, class F>
+    inline decltype(auto) For<L,Xpr,F>::operator()() const
     {
+	typedef typename std::remove_reference<decltype(mXpr(mExt * UPos(0)))>::type OutT;
+	auto o = OutT();
 	for(SizeT i = 0; i != mSize; ++i){
 	    const auto pos = mExt * UPos(i);
-	    mXpr(pos);
+	    mF(o, mXpr(pos));
 	}
-	return 0;
+	return o;
     }
 
-    template <SizeT L, class Xpr>
+    template <SizeT L, class Xpr, class F>
     template <SizeT I>
-    inline decltype(auto) For<L,Xpr>::rootSteps(const IndexId<I>& id) const
+    inline decltype(auto) For<L,Xpr,F>::rootSteps(const IndexId<I>& id) const
     {
 	return mXpr.rootSteps(id);
     }
@@ -53,60 +74,61 @@ namespace CNORXZ
      *   SFor   *
      ************/
 
-    template <SizeT N, SizeT L, class Xpr>
-    constexpr SFor<N,L,Xpr>::SFor(const IndexId<L>& id, const Xpr& xpr) :
+    template <SizeT N, SizeT L, class Xpr, class F>
+    constexpr SFor<N,L,Xpr,F>::SFor(const IndexId<L>& id, const Xpr& xpr, F&& f) :
 	mId(id),
 	mXpr(xpr),
-	mExt(mXpr.RootSteps(mId))
+	mExt(mXpr.RootSteps(mId)),
+	mF(f)
     {}
 
-    template <SizeT N, SizeT L, class Xpr>
+    template <SizeT N, SizeT L, class Xpr, class F>
     template <class PosT>
-    constexpr SizeT SFor<N,L,Xpr>::operator()(const PosT& last) const
+    constexpr decltype(auto) SFor<N,L,Xpr,F>::operator()(const PosT& last) const
     {
 	return exec<0>(last);
     }
 	
-    template <SizeT N, SizeT L, class Xpr>
-    constexpr SizeT SFor<N,L,Xpr>::operator()() const
+    template <SizeT N, SizeT L, class Xpr, class F>
+    constexpr decltype(auto) SFor<N,L,Xpr,F>::operator()() const
     {
 	return exec<0>();
     }
 
-    template <SizeT N, SizeT L, class Xpr>
+    template <SizeT N, SizeT L, class Xpr, class F>
     template <SizeT I>
-    constexpr decltype(auto) SFor<N,L,Xpr>::rootSteps(const IndexId<I>& id) const
+    constexpr decltype(auto) SFor<N,L,Xpr,F>::rootSteps(const IndexId<I>& id) const
     {
 	return mXpr.rootSteps(id);
     }
 
-    template <SizeT N, SizeT L, class Xpr>
+    template <SizeT N, SizeT L, class Xpr, class F>
     template <SizeT I, class PosT>
-    constexpr SizeT SFor<N,L,Xpr>::exec(const PosT& last) const
+    constexpr decltype(auto) SFor<N,L,Xpr,F>::exec(const PosT& last) const
     {
 	constexpr SPos<I> i;
 	const auto pos = last + mExt * i;
-	mXpr(pos);
+	auto o = mXpr(pos);
 	if constexpr(I < N-1){
-	    return exec<I+1>(last);
+	    return mF(o,exec<I+1>(last));
 	}
 	else {
-	    return 0;
+	    return o;
 	}
     }
 
-    template <SizeT N, SizeT L, class Xpr>
+    template <SizeT N, SizeT L, class Xpr, class F>
     template <SizeT I>
-    constexpr SizeT SFor<N,L,Xpr>::exec() const
+    constexpr decltype(auto) SFor<N,L,Xpr,F>::exec() const
     {
 	constexpr SPos<I> i;
 	const auto pos = mExt * i;
-	mXpr(pos);
+	auto o = mXpr(pos);
 	if constexpr(I < N-1){
-	    return exec<I+1>();
+	    return mF(o,exec<I+1>());
 	}
 	else {
-	    return 0;
+	    return o;
 	}
     }
 
@@ -114,12 +136,13 @@ namespace CNORXZ
      *   TFor   *
      ************/
 
-    template <SizeT L, class Xpr>
-    constexpr TFor<L,Xpr>::TFor(SizeT size, const IndexId<L>& id, const Xpr& xpr) :
+    template <SizeT L, class Xpr, class F>
+    constexpr TFor<L,Xpr,F>::TFor(SizeT size, const IndexId<L>& id, const Xpr& xpr, F&& f) :
 	mSize(size),
 	mId(id),
 	mXpr(xpr),
-	mExt(mXpr.rootSteps(mId))
+	mExt(mXpr.rootSteps(mId)),
+	mF(f)
     {
 	// check for write access!!!
 	/*
@@ -130,12 +153,18 @@ namespace CNORXZ
 	*/
     }
 
-    template <SizeT L, class Xpr>
+    template <SizeT L, class Xpr, class F>
     template <class PosT>
-    inline SizeT TFor<L,Xpr>::operator()(const PosT& last) const
+    inline decltype(auto) TFor<L,Xpr,F>::operator()(const PosT& last) const
     {
+	typedef typename std::remove_reference<decltype(mXpr(last + mExt * UPos(0)))>::type OutT;
 	int i = 0;
 	const int size = static_cast<int>(mSize);
+	Vector<OutT> ov;
+	if constexpr(not std::is_same<F,ZeroF>::value){
+	    // replace if statement by "does s.th. non-trivial"
+	    ov.resize(mSize);
+	}
 #pragma omp parallel private(i)
 	{
 	    auto xpr = mXpr;
@@ -145,14 +174,27 @@ namespace CNORXZ
 		xpr(pos);
 	    }
 	}
-	return 0;
+	OutT o;
+	if constexpr(not std::is_same<F,ZeroF>::value){
+	    // replace if statement by "does s.th. non-trivial"
+	    for(SizeT i = 0; i != mSize; ++i){
+		mF(o, ov[i]);
+	    }
+	}
+	return o;
     }
 	
-    template <SizeT L, class Xpr>
-    inline SizeT TFor<L,Xpr>::operator()() const
+    template <SizeT L, class Xpr, class F>
+    inline decltype(auto) TFor<L,Xpr,F>::operator()() const
     {
+	typedef typename std::remove_reference<decltype(mXpr(mExt * UPos(0)))>::type OutT;
 	int i = 0;
 	const int size = static_cast<int>(mSize);
+	Vector<OutT> ov;
+	if constexpr(not std::is_same<F,ZeroF>::value){
+	    // replace if statement by "does s.th. non-trivial"
+	    ov.resize(mSize);
+	}
 #pragma omp parallel private(i)
 	{
 	    auto xpr = mXpr;
@@ -162,12 +204,19 @@ namespace CNORXZ
 		xpr(pos);
 	    }
 	}
-	return 0;
+	OutT o;
+	if constexpr(not std::is_same<F,ZeroF>::value){
+	    // replace if statement by "does s.th. non-trivial"
+	    for(SizeT i = 0; i != mSize; ++i){
+		mF(o, ov[i]);
+	    }
+	}
+	return o;
     }
 
-    template <SizeT L, class Xpr>
+    template <SizeT L, class Xpr, class F>
     template <SizeT I>
-    inline decltype(auto) TFor<L,Xpr>::rootSteps(const IndexId<I>& id) const
+    inline decltype(auto) TFor<L,Xpr,F>::rootSteps(const IndexId<I>& id) const
     {
 	return mXpr.rootSteps(id);
     }
@@ -186,14 +235,14 @@ namespace CNORXZ
 
     template <SizeT N, SizeT L, class Xpr>
     template <class PosT>
-    constexpr SizeT EFor<N,L,Xpr>::operator()(const PosT& last) const
+    constexpr decltype(auto) EFor<N,L,Xpr>::operator()(const PosT& last) const
     {
 	auto pos = mkEPos<N>(last, mExt);
 	return mXpr(pos);
     }
 	
     template <SizeT N, SizeT L, class Xpr>
-    constexpr SizeT EFor<N,L,Xpr>::operator()() const
+    constexpr decltype(auto) EFor<N,L,Xpr>::operator()() const
     {
 	auto pos = mkEPos<N>(SPos<0>(), mExt);
 	return mXpr(pos);
diff --git a/src/include/xpr/for.h b/src/include/xpr/for.h
index 8b7771e..0759ad8 100644
--- a/src/include/xpr/for.h
+++ b/src/include/xpr/for.h
@@ -7,19 +7,32 @@
 
 namespace CNORXZ
 {
+    class ZeroF
+    {
+    public:
+	template <typename... T>
+	constexpr decltype(auto) operator()(const T&... as) const;
+    };
 
-    template <SizeT L, class Xpr>
-    class For : public XprInterface<For<L,Xpr>>
+    class NoF
+    {
+    public:
+	template <typename... T>
+	constexpr decltype(auto) operator()(const T&... as) const;
+    };
+
+    template <SizeT L, class Xpr, class F = NoF>
+    class For : public XprInterface<For<L,Xpr,F>>
     {
     public:
 	DEFAULT_MEMBERS(For);
 
-	constexpr For(SizeT size, const IndexId<L>& id, const Xpr& xpr);
+	constexpr For(SizeT size, const IndexId<L>& id, const Xpr& xpr, F&& f);
 
 	template <class PosT>
-	inline SizeT operator()(const PosT& last) const;
+	inline decltype(auto) operator()(const PosT& last) const;
 	
-	inline SizeT operator()() const;
+	inline decltype(auto) operator()() const;
 
 	template <SizeT I>
 	inline decltype(auto) rootSteps(const IndexId<I>& id) const;
@@ -30,22 +43,23 @@ namespace CNORXZ
 	Xpr mXpr;
 	typedef decltype(mXpr.rootSteps(mId)) XPosT;
 	XPosT mExt;
-	
+	F mF;
     };
 
+    
     // unrolled loop:
-    template <SizeT N, SizeT L, class Xpr>
-    class SFor : public XprInterface<SFor<N,L,Xpr>>
+    template <SizeT N, SizeT L, class Xpr, class F = NoF>
+    class SFor : public XprInterface<SFor<N,L,Xpr,F>>
     {
     public:
 	DEFAULT_MEMBERS(SFor);
 
-	constexpr SFor(const IndexId<L>& id, const Xpr& xpr);
+	constexpr SFor(const IndexId<L>& id, const Xpr& xpr, F&& f);
 
 	template <class PosT>
-	constexpr SizeT operator()(const PosT& last) const;
+	constexpr decltype(auto) operator()(const PosT& last) const;
 	
-	constexpr SizeT operator()() const;
+	constexpr decltype(auto) operator()() const;
 
 	template <SizeT I>
 	constexpr decltype(auto) rootSteps(const IndexId<I>& id) const;
@@ -53,32 +67,33 @@ namespace CNORXZ
     private:
 
 	template <SizeT I, class PosT>
-	constexpr SizeT exec(const PosT& last) const;
+	constexpr decltype(auto) exec(const PosT& last) const;
 
 	template <SizeT I>
-	constexpr SizeT exec() const;
+	constexpr decltype(auto) exec() const;
 
 	IndexId<L> mId;
 	Xpr mXpr;
 	typedef decltype(mXpr.RootSteps(mId)) XPosT;
 	XPosT mExt;
+	F mF;
 	
     };
 
 
     // multi-threading
-    template <SizeT L, class Xpr>
-    class TFor : public XprInterface<TFor<L,Xpr>>
+    template <SizeT L, class Xpr, class F = NoF>
+    class TFor : public XprInterface<TFor<L,Xpr,F>>
     {
     public:
 	DEFAULT_MEMBERS(TFor);
 
-	constexpr TFor(SizeT size, const IndexId<L>& id, const Xpr& xpr);
+	constexpr TFor(SizeT size, const IndexId<L>& id, const Xpr& xpr, F&& f);
 
 	template <class PosT>
-	inline SizeT operator()(const PosT& last) const;
+	inline decltype(auto) operator()(const PosT& last) const;
 	
-	inline SizeT operator()() const;
+	inline decltype(auto) operator()() const;
 
 	template <SizeT I>
 	inline decltype(auto) rootSteps(const IndexId<I>& id) const;
@@ -89,7 +104,7 @@ namespace CNORXZ
 	Xpr mXpr;
 	typedef decltype(mXpr.rootSteps(mId)) XPosT;
 	XPosT mExt;
-	
+	F mF;
     };
 
     // Extension For (Vectorization)
@@ -102,9 +117,9 @@ namespace CNORXZ
 	constexpr EFor(const IndexId<L>& id, const Xpr& xpr);
 
 	template <class PosT>
-	constexpr SizeT operator()(const PosT& last) const;
+	constexpr decltype(auto) operator()(const PosT& last) const;
 
-	constexpr SizeT operator()() const;
+	constexpr decltype(auto) operator()() const;
 
 	template <SizeT I>
 	constexpr decltype(auto) rootSteps(const IndexId<I>& id) const;
diff --git a/src/include/xpr/for_type.h b/src/include/xpr/for_type.h
deleted file mode 100644
index 6833c99..0000000
--- a/src/include/xpr/for_type.h
+++ /dev/null
@@ -1,23 +0,0 @@
-
-#ifndef __cxz_for_type_h__
-#define __cxz_for_type_h__
-
-namespace CNORXZInternal
-{
-    
-    enum class ForType {
-	DEFAULT = 0,
-	HIDDEN = 1	
-    };
-
-    template <class IndexClass, class Expr, ForType FT = ForType::DEFAULT, size_t DIV = 1>
-    class For;
-
-    template <class IndexClass, class Expr, size_t DIV = 1>
-    class PFor;
-
-    
-} // end namespace CNORXZInternal
-
-
-#endif
diff --git a/src/include/xpr/xpr_base.cc.h b/src/include/xpr/xpr_base.cc.h
index 99da39a..dcc24c5 100644
--- a/src/include/xpr/xpr_base.cc.h
+++ b/src/include/xpr/xpr_base.cc.h
@@ -11,31 +11,31 @@ namespace CNORXZ
      *   VXpr   *
      ************/
     
-    template <class Xpr>
-    VXpr<Xpr>::VXpr(const XprInterface<Xpr>& a) :
+    template <typename T, class Xpr>
+    VXpr<T,Xpr>::VXpr(const XprInterface<Xpr>& a) :
 	Xpr(a.THIS())
     {}
 
-    template <class Xpr>
-    Uptr<VXprBase> VXpr<Xpr>::copy() const
+    template <typename T, class Xpr>
+    Uptr<VXprBase<T>> VXpr<T,Xpr>::copy() const
     {
-	return std::make_unique<VXpr<Xpr>>(*this);
+	return std::make_unique<VXpr<T,Xpr>>(*this);
     }
 
-    template <class Xpr>
-    SizeT VXpr<Xpr>::vexec(const DPos& last) const
+    template <typename T, class Xpr>
+    T VXpr<T,Xpr>::vexec(const DPos& last) const
     {
 	return (*this)(last);
     }
 
-    template <class Xpr>
-    SizeT VXpr<Xpr>::vexec() const
+    template <typename T, class Xpr>
+    T VXpr<T,Xpr>::vexec() const
     {
 	return (*this)();
     }
     
-    template <class Xpr>
-    DPos VXpr<Xpr>::vrootSteps(const IndexId<0>& id) const
+    template <typename T, class Xpr>
+    DPos VXpr<T,Xpr>::vrootSteps(const IndexId<0>& id) const
     {
 	return DPos(this->rootSteps(id));
     }
@@ -44,25 +44,29 @@ namespace CNORXZ
      *   DXpr   *
      ************/
 
+    template <typename T>
     template <class Xpr>
-    DXpr::DXpr(const Xpr& a) :
-	ObjHandle<VXprBase>(std::make_unique<VXpr<Xpr>>(a))
+    DXpr<T>::DXpr(const Xpr& a) :
+	ObjHandle<VXprBase<T>>(std::make_unique<VXpr<T,Xpr>>(a))
     {}
 
-    inline SizeT DXpr::operator()(const DPos& last) const
+    template <typename T>
+    inline T DXpr<T>::operator()(const DPos& last) const
     {
-	return mC->vexec(last);
+	return VB::mC->vexec(last);
     }
     
-    inline SizeT DXpr::operator()() const
+    template <typename T>
+    inline T DXpr<T>::operator()() const
     {
-	return mC->vexec();
+	return VB::mC->vexec();
     }
 
+    template <typename T>
     template <SizeT I>
-    inline DPos DXpr::rootSteps(const IndexId<I>& id) const
+    inline DPos DXpr<T>::rootSteps(const IndexId<I>& id) const
     {
-	return mC->vrootSteps(IndexId<0>(id.id()));
+	return VB::mC->vrootSteps(IndexId<0>(id.id()));
     }
 }
 
diff --git a/src/include/xpr/xpr_base.h b/src/include/xpr/xpr_base.h
index aad448d..53bf760 100644
--- a/src/include/xpr/xpr_base.h
+++ b/src/include/xpr/xpr_base.h
@@ -19,53 +19,57 @@ namespace CNORXZ
 	inline const Xpr& THIS() const { return static_cast<const Xpr&>(*this); }
 
 	template <class PosT>
-	inline SizeT operator()(const PosT& last) const { return THIS()(last); } 
+	constexpr decltype(auto) operator()(const PosT& last) const { return THIS()(last); } 
 
-	inline SizeT operator()() const { return THIS()(); }
+	constexpr decltype(auto) operator()() const { return THIS()(); }
 
 	template <SizeT I>
 	constexpr decltype(auto) rootSteps(const IndexId<I>& id) const { return THIS().rootSteps(id); }
     };
 
+    template <typename T>
     class VXprBase
     {
     public:
 	DEFAULT_MEMBERS(VXprBase);
 
-	virtual Uptr<VXprBase> copy() const = 0;
+	virtual Uptr<VXprBase<T>> copy() const = 0;
 
-	virtual SizeT vexec(const DPos& last) const = 0;
-	virtual SizeT vexec() const = 0;
+	virtual T vexec(const DPos& last) const = 0;
+	virtual T vexec() const = 0;
 
 	virtual DPos vrootSteps(const IndexId<0>& id) const = 0;
     };
 
-    template <class Xpr>
-    class VXpr : public VXprBase, public Xpr
+    template <typename T, class Xpr>
+    class VXpr : public VXprBase<T>, public Xpr
     {
     public:
 	DEFAULT_MEMBERS(VXpr);
 	VXpr(const XprInterface<Xpr>& a);
 
-	virtual Uptr<VXprBase> copy() const override final;
+	virtual Uptr<VXprBase<T>> copy() const override final;
 
-	virtual SizeT vexec(const DPos& last) const override final;
-	virtual SizeT vexec() const override final;
+	virtual T vexec(const DPos& last) const override final;
+	virtual T vexec() const override final;
 
 	virtual DPos vrootSteps(const IndexId<0>& id) const override final;
     };
 
-    class DXpr : public ObjHandle<VXprBase>,
-		 public XprInterface<DXpr>
+    template <typename T>
+    class DXpr : public ObjHandle<VXprBase<T>>,
+		 public XprInterface<DXpr<T>>
     {
     public:
+	typedef ObjHandle<VXprBase<T>> VB;
+	
 	DEFAULT_MEMBERS(DXpr);
 
 	template <class Xpr>
 	explicit DXpr(const Xpr& a);
 
-	inline SizeT operator()(const DPos& last) const;
-	inline SizeT operator()() const;
+	inline T operator()(const DPos& last) const;
+	inline T operator()() const;
 
 	template <SizeT I>
 	inline DPos rootSteps(const IndexId<I>& id) const;
diff --git a/src/lib/ranges/dindex.cc b/src/lib/ranges/dindex.cc
index e2ae9fb..c31d910 100644
--- a/src/lib/ranges/dindex.cc
+++ b/src/lib/ranges/dindex.cc
@@ -138,9 +138,9 @@ namespace CNORXZ
 	return *this;
     }
 
-    DXpr DIndex::ifor(const DXpr& xpr) const
+    DXpr<SizeT> DIndex::ifor(const DXpr<SizeT>& xpr) const
     {
-	return DXpr(mI->ifor(xpr));
+	return DXpr<SizeT>(mI->ifor(xpr));
     }
     
 }
diff --git a/src/lib/ranges/yrange.cc b/src/lib/ranges/yrange.cc
index 3789c99..f1f09f8 100644
--- a/src/lib/ranges/yrange.cc
+++ b/src/lib/ranges/yrange.cc
@@ -171,10 +171,10 @@ namespace CNORXZ
 	return *this;
     }
 
-    DXpr YIndex::ifor(const DXpr& xpr) const
+    DXpr<SizeT> YIndex::ifor(const DXpr<SizeT>& xpr) const
     {
 	assert(0);
-	return DXpr();
+	return DXpr<SizeT>();
     }
     
     /**********************