remove master class + add separate parallel op root + pfor expression in parallel scope

src/include/expressions.h

@@ -30,10 +30,10 @@ namespace MultiArrayTools
     using OX = Operation<double,OpF<double>,oo<EC,MAs>...>;
 
     template <class EC, class Second>
-    using AEXT = typename OperationMaster<double,SelfIdentity<double>,Second,DynamicRange<EC>>::AssignmentExpr;
+    using AEXT = AssignmentExpr<double,Second>;
 
     template <class EC, class Second>
-    using AEXT_P = typename OperationMaster<double,plus<double>,Second,DynamicRange<EC>>::AssignmentExpr;
+    using AEXT_P = AssignmentExpr<double,Second>;
 
     template <class EC, template <class> class OpF, class... MAs>
     using AEX = AEXT<EC,OX<EC,OpF,MAs...>>;

src/include/mbase_def.h

@@ -35,14 +35,14 @@ namespace MultiArrayTools
     template <typename T, class OperationClass>
     class OperationTemplate;
     
-    // multi_array_operation.h
-    template <typename T, class AOp, class OpClass, class... Ranges>
-    class OperationMaster;
-    
     // multi_array_operation.h
     template <typename T, class... Ranges>
     class OperationRoot;
 
+    // multi_array_operation.h
+    template <typename T, class... Ranges>
+    class ParallelOperationRoot;
+
     // multi_array_operation.h
     template <typename T>
     class OperationValue;

src/include/multi_array_operation.cc.h

@@ -117,108 +117,42 @@ namespace MultiArrayTools
         return Operation<R,function<R,T,typename Args::value_type...>,OperationClass, Args...>(ll, THIS(), args...);
     }
 
-    /*****************************************
+    /************************
-     *   OperationMaster::AssignmentExpr     *
+     *   AssignmentExpr     *
-     *****************************************/
+     ************************/
-    /*
+
-    template <typename T, class AOp, class OpClass, class... Ranges>
+    template <typename T, class OpClass>
-    OperationMaster<T,AOp,OpClass,Ranges...>::AssignmentExpr::
+    AssignmentExpr<T,OpClass>::AssignmentExpr(T* dataPtr, const OpClass& sec) :
-    AssignmentExpr(OperationMaster& m, const OpClass& sec) :
-	mM(m), mSec(sec) {}
-    */
-    template <typename T, class AOp, class OpClass, class... Ranges>
-    OperationMaster<T,AOp,OpClass,Ranges...>::AssignmentExpr::
-    AssignmentExpr(T* dataPtr, const OpClass& sec) :
         mSec(sec), mDataPtr(dataPtr) {}
     
-    template <typename T, class AOp, class OpClass, class... Ranges>
+    template <typename T, class OpClass>
-    inline void OperationMaster<T,AOp,OpClass,Ranges...>::AssignmentExpr::
+    inline void AssignmentExpr<T,OpClass>::operator()(size_t start, ExtType last)
-    operator()(size_t start, ExtType last)
     {
-	//VCHECK(mSec.template get<ExtType>(last));
-	//mM.set(start, mSec.template get<ExtType>(last) );
         mDataPtr[start] = mSec.template get<ExtType>(last);
-        //AOp::sapply(mDataPtr[start], mSec.template get<ExtType>(last));
     }
     
-    template <typename T, class AOp, class OpClass, class... Ranges>
+    template <typename T, class OpClass>
-    typename OperationMaster<T,AOp,OpClass,Ranges...>::AssignmentExpr::ExtType
+    typename AssignmentExpr<T,OpClass>::ExtType AssignmentExpr<T,OpClass>::rootSteps(std::intptr_t iPtrNum) const
-    OperationMaster<T,AOp,OpClass,Ranges...>::AssignmentExpr::
-    rootSteps(std::intptr_t iPtrNum) const
     {
 	return mSec.rootSteps(iPtrNum);
     }
 
-    template <typename T, class AOp, class OpClass, class... Ranges>
+    template <typename T, class OpClass>
-    OperationMaster<T,AOp,OpClass,Ranges...>::AddExpr::
+    AddExpr<T,OpClass>::AddExpr(T* dataPtr, const OpClass& sec) :
-    AddExpr(T* dataPtr, const OpClass& sec) :
         mSec(sec), mDataPtr(dataPtr) {}
     
-    template <typename T, class AOp, class OpClass, class... Ranges>
+    template <typename T, class OpClass>
-    inline void OperationMaster<T,AOp,OpClass,Ranges...>::AddExpr::
+    inline void AddExpr<T,OpClass>::operator()(size_t start, ExtType last)
-    operator()(size_t start, ExtType last)
     {
-	//VCHECK(mSec.template get<ExtType>(last));
-	//mM.set(start, mSec.template get<ExtType>(last) );
         mDataPtr[start] += mSec.template get<ExtType>(last);
-        //AOp::sapply(mDataPtr[start], mSec.template get<ExtType>(last));
     }
     
-    template <typename T, class AOp, class OpClass, class... Ranges>
+    template <typename T, class OpClass>
-    typename OperationMaster<T,AOp,OpClass,Ranges...>::AddExpr::ExtType
+    typename AddExpr<T,OpClass>::ExtType AddExpr<T,OpClass>::rootSteps(std::intptr_t iPtrNum) const
-    OperationMaster<T,AOp,OpClass,Ranges...>::AddExpr::
-    rootSteps(std::intptr_t iPtrNum) const
     {
 	return mSec.rootSteps(iPtrNum);
     }
 
-    
-    /*************************
-     *   OperationMaster     *
-     *************************/
-
-    template <typename T, class AOp, class OpClass, class... Ranges>
-    OperationMaster<T,AOp,OpClass,Ranges...>::
-    OperationMaster(MutableMultiArrayBase<T,Ranges...>& ma, const OpClass& second,
-		    IndexType& index, bool doParallel) :
-	mSecond(second), mDataPtr(ma.data()),
-	mIndex(index), mDoParallel(doParallel)
-    {
-	performAssignment(0);
-    }
-
-    template <typename T, class AOp, class OpClass, class... Ranges>
-    OperationMaster<T,AOp,OpClass,Ranges...>::
-    OperationMaster(T* data, const OpClass& second,
-		    IndexType& index, bool doParallel) :
-	mSecond(second), mDataPtr(data),
-	mIndex(index), mDoParallel(doParallel)
-    {
-	performAssignment(0);
-    }
-
-    template <typename T, class AOp, class OpClass, class... Ranges>
-    void OperationMaster<T,AOp,OpClass,Ranges...>::performAssignment(std::intptr_t blockIndexNum)
-    {
-        /*
-	AssignmentExpr ae(*this, mSecond); // Expression to be executed within loop
-        if(mDoParallel){
-            auto ploop = mIndex.pifor( 1, mSecond.loop(ae) );
-            ploop(); // execute overall loop(s) and so internal hidden loops and so the inherited expressions
-        }
-        else {
-            auto loop = mIndex.ifor( 1, mSecond.loop(ae) );
-            loop(); // execute overall loop(s) and so internal hidden loops and so the inherited expressions
-        }
-        */
-    }
-
-    template <typename T, class AOp, class OpClass, class... Ranges>
-    inline T OperationMaster<T,AOp,OpClass,Ranges...>::get(size_t pos) const
-    {
-	return mDataPtr[pos];
-    }
-
    
     /****************************
      *   ConstOperationRoot     *
@@ -232,10 +166,8 @@ namespace MultiArrayTools
         mOrigDataPtr(ma.data()),
 	mIndex( ma.begin() )
     {
-	//VCHECK(ma.data());
 	mIndex(indices...);
         mDataPtr = mOrigDataPtr + mIndex.pos();
-	//mOff = mIndex.pos();
     }
 
     template <typename T, class... Ranges>
@@ -397,7 +329,6 @@ namespace MultiArrayTools
     {
 	mIndex(indices...);
         mDataPtr = mOrigDataPtr + mIndex.pos();
-	//mOff = mIndex.pos();
     }
 
     template <typename T, class... Ranges>
@@ -408,67 +339,57 @@ namespace MultiArrayTools
 	mIndex( ind )
     {
         mDataPtr = mOrigDataPtr + mIndex.pos();
-	//mOff = mIndex.pos();
+    }
+
+    template <typename T, class... Ranges>
+    template <class OpClass>
+    auto OperationRoot<T,Ranges...>::assign(const OpClass& in)
+        -> decltype(mIndex.ifor(1,in.loop(AssignmentExpr<T,OpClass>(mOrigDataPtr,in))))
+    {
+        return mIndex.ifor(1,in.loop(AssignmentExpr<T,OpClass>(mOrigDataPtr,in)));
+    }
+
+    template <typename T, class... Ranges>
+    template <class OpClass>
+    auto OperationRoot<T,Ranges...>::plus(const OpClass& in)
+        -> decltype(mIndex.ifor(1,in.loop(AddExpr<T,OpClass>(mOrigDataPtr,in))))
+    {
+        return mIndex.ifor(1,in.loop(AddExpr<T,OpClass>(mOrigDataPtr,in)));
     }
 
     template <typename T, class... Ranges>
     template <class OpClass>
     OperationRoot<T,Ranges...>& OperationRoot<T,Ranges...>::operator=(const OpClass& in)
-    //OperationMaster<T,SelfIdentity<T>,OpClass,Ranges...> OperationRoot<T,Ranges...>::operator=(const OpClass& in)
     {
-	typename OperationMaster<T,SelfIdentity<T>,OpClass,Ranges...>::AssignmentExpr ae(mOrigDataPtr, in);
+        assign(in)();
-        // Expression to be executed within loop
-        if(mDoParallel){
-            auto ploop = mIndex.pifor( 1, in.loop(ae) );
-            ploop(); // execute overall loop(s) and so internal hidden loops and so the inherited expressions
-        }
-        else {
-            auto loop = mIndex.ifor( 1, in.loop(ae) );
-            loop(); // execute overall loop(s) and so internal hidden loops and so the inherited expressions
-        }
         return *this;
-	//return OperationMaster<T,SelfIdentity<T>,OpClass,Ranges...>(mDataPtr, in, mIndex, mDoParallel);
     }
 
     template <typename T, class... Ranges>
     template <class OpClass>
     OperationRoot<T,Ranges...>& OperationRoot<T,Ranges...>::operator+=(const OpClass& in)
-    //OperationMaster<T,plus<T>,OpClass,Ranges...> OperationRoot<T,Ranges...>::operator+=(const OpClass& in)
     {
-        typename OperationMaster<T,plus<T>,OpClass,Ranges...>::AddExpr ae(mOrigDataPtr, in);
+        plus(in)();
-        // Expression to be executed within loop
-        if(mDoParallel){
-            auto ploop = mIndex.pifor( 1, in.loop(ae) );
-            ploop(); // execute overall loop(s) and so internal hidden loops and so the inherited expressions
-        }
-        else {
-            auto loop = mIndex.ifor( 1, in.loop(ae) );
-            loop(); // execute overall loop(s) and so internal hidden loops and so the inherited expressions
-        }
         return *this;
-	//return OperationMaster<T,plus<T>,OpClass,Ranges...>(mDataPtr, in, mIndex, mDoParallel);
     }
 
     template <typename T, class... Ranges>
     OperationRoot<T,Ranges...>& OperationRoot<T,Ranges...>::operator=(const OperationRoot<T,Ranges...>& in)
-    //OperationMaster<T,SelfIdentity<T>,OperationRoot<T,Ranges...>,Ranges...>
-    //OperationRoot<T,Ranges...>::operator=(const OperationRoot<T,Ranges...>& in)
     {
 	return operator=<OperationRoot<T,Ranges...> >(in);
     }
 
     template <typename T, class... Ranges>
-    OperationRoot<T,Ranges...>& OperationRoot<T,Ranges...>::par()
+    ParallelOperationRoot<T,Ranges...> OperationRoot<T,Ranges...>::par()
     {
-        mDoParallel = true;
+        return ParallelOperationRoot<T,Ranges...>(mOrigDataPtr, mIndex);
-        return *this;
     }
     
     template <typename T, class... Ranges>
     template <class ET>
     inline T& OperationRoot<T,Ranges...>::get(ET pos) const
     {
-	return mDataPtr[pos.val()/*+mOff*/];
+	return mDataPtr[pos.val()];
     }
 
     template <typename T, class... Ranges>
@@ -484,7 +405,6 @@ namespace MultiArrayTools
     MExt<void> OperationRoot<T,Ranges...>::rootSteps(std::intptr_t iPtrNum) const
     {
 	return MExt<void>(getStepSize( mIndex, iPtrNum ));
-	//return MExt<void>(getStepSize( mIndex.info(), iPtrNum ));
     }
 
     template <typename T, class... Ranges>
@@ -511,6 +431,108 @@ namespace MultiArrayTools
 	return out;
     }
 
+    /*******************************
+     *   ParallelOperationRoot     *
+     *******************************/
+    
+    template <typename T, class... Ranges>
+    ParallelOperationRoot<T,Ranges...>::
+    ParallelOperationRoot(MutableMultiArrayBase<T,Ranges...>& ma,
+		  const std::shared_ptr<typename Ranges::IndexType>&... indices) :
+	mDataPtr(ma.data()),
+        mOrigDataPtr(ma.data()),
+	mIndex( ma.begin() )
+    {
+	mIndex(indices...);
+        mDataPtr = mOrigDataPtr + mIndex.pos();
+    }
+
+    template <typename T, class... Ranges>
+    ParallelOperationRoot<T,Ranges...>::
+    ParallelOperationRoot(T* data, const IndexType& ind) :
+	mDataPtr(data),
+        mOrigDataPtr(data),
+	mIndex( ind )
+    {
+        mDataPtr = mOrigDataPtr + mIndex.pos();
+    }
+
+    template <typename T, class... Ranges>
+    template <class OpClass>
+    auto ParallelOperationRoot<T,Ranges...>::assign(const OpClass& in)
+        -> decltype(mIndex.pifor(1,in.loop(AssignmentExpr<T,OpClass>(mOrigDataPtr,in))))
+    {
+        return mIndex.pifor(1,in.loop(AssignmentExpr<T,OpClass>(mOrigDataPtr,in)));
+    }
+
+    template <typename T, class... Ranges>
+    template <class OpClass>
+    auto ParallelOperationRoot<T,Ranges...>::plus(const OpClass& in)
+        -> decltype(mIndex.pifor(1,in.loop(AddExpr<T,OpClass>(mOrigDataPtr,in))))
+    {
+        return mIndex.pifor(1,in.loop(AddExpr<T,OpClass>(mOrigDataPtr,in)));
+    }
+
+    template <typename T, class... Ranges>
+    template <class OpClass>
+    ParallelOperationRoot<T,Ranges...>& ParallelOperationRoot<T,Ranges...>::operator=(const OpClass& in)
+    {
+        assign(in)();
+        return *this;
+    }
+
+    template <typename T, class... Ranges>
+    template <class OpClass>
+    ParallelOperationRoot<T,Ranges...>& ParallelOperationRoot<T,Ranges...>::operator+=(const OpClass& in)
+    {
+        plus(in)();
+        return *this;
+    }
+
+    template <typename T, class... Ranges>
+    ParallelOperationRoot<T,Ranges...>&
+    ParallelOperationRoot<T,Ranges...>::operator=(const ParallelOperationRoot<T,Ranges...>& in)
+    {
+	return operator=<ParallelOperationRoot<T,Ranges...> >(in);
+    }
+    
+    template <typename T, class... Ranges>
+    template <class ET>
+    inline T& ParallelOperationRoot<T,Ranges...>::get(ET pos) const
+    {
+	return mDataPtr[pos.val()/*+mOff*/];
+    }
+
+    template <typename T, class... Ranges>
+    template <class ET>
+    inline ParallelOperationRoot<T,Ranges...>& ParallelOperationRoot<T,Ranges...>::set(ET pos)
+    {
+	mIndex = pos.val();
+	mDataPtr = mOrigDataPtr + mIndex.pos();
+	return *this;
+    }
+
+    template <typename T, class... Ranges>
+    MExt<void> ParallelOperationRoot<T,Ranges...>::rootSteps(std::intptr_t iPtrNum) const
+    {
+	return MExt<void>(getStepSize( mIndex, iPtrNum ));
+	//return MExt<void>(getStepSize( mIndex.info(), iPtrNum ));
+    }
+
+    template <typename T, class... Ranges>
+    template <class Expr>
+    Expr ParallelOperationRoot<T,Ranges...>::loop(Expr exp) const
+    {
+	return exp;
+    }
+
+    template <typename T, class... Ranges>
+    T* ParallelOperationRoot<T,Ranges...>::data() const
+    {
+        auto i = mIndex;
+	return mOrigDataPtr + i().pos();
+    }
+
 
     /************************
      *   OperationValue     *

src/include/multi_array_operation.h

@@ -101,95 +101,54 @@ namespace MultiArrayTools
         }
     };
 
-    template <typename T, class AOp, class OpClass, class... Ranges>
+    template <typename T, class OpClass>
-    class OperationMaster
+    class AssignmentExpr
     {
+    private:
+        AssignmentExpr() = default;
+	    	    
+        OpClass mSec;
+        T* mDataPtr;
+            
     public:
 
-	class AssignmentExpr
+        static constexpr size_t LAYER = 0;
-	{
+        static constexpr size_t SIZE = OpClass::SIZE;
-	private:
+        typedef decltype(mSec.rootSteps()) ExtType;
-	    AssignmentExpr() = default;
-	    	    
-	    //OperationMaster mM;
-	    OpClass mSec;
-            T* mDataPtr;
-            
-	public:
-
-	    static constexpr size_t LAYER = 0;
-	    static constexpr size_t SIZE = OpClass::SIZE;
-	    typedef decltype(mSec.rootSteps()) ExtType;
 	    
-	    //AssignmentExpr(OperationMaster& m, const OpClass& sec);
+        AssignmentExpr(T* dataPtr, const OpClass& sec);
-            AssignmentExpr(T* dataPtr, const OpClass& sec);
+        AssignmentExpr(const AssignmentExpr& in) = default;
-
+        AssignmentExpr(AssignmentExpr&& in) = default;
-	    AssignmentExpr(const AssignmentExpr& in) = default;
-	    AssignmentExpr(AssignmentExpr&& in) = default;
 	    
-	    inline void operator()(size_t start = 0); 
+        inline void operator()(size_t start = 0); 
-	    inline void operator()(size_t start, ExtType last);
+        inline void operator()(size_t start, ExtType last);
-
+        auto rootSteps(std::intptr_t iPtrNum = 0) const -> ExtType;
-	    auto rootSteps(std::intptr_t iPtrNum = 0) const -> ExtType;
+    };
-	    
+
-	};
+    template <typename T, class OpClass>
-
+    class AddExpr
-        class AddExpr
+    {
-	{
+    private:
-	private:
+        AddExpr() = default;
-	    AddExpr() = default;
+	    	    
-	    	    
+        OpClass mSec;
-	    //OperationMaster mM;
+        T* mDataPtr;
-	    OpClass mSec;
+            
-            T* mDataPtr;
+    public:
-            
+
-	public:
+        static constexpr size_t LAYER = 0;
-
+        static constexpr size_t SIZE = OpClass::SIZE;
-	    static constexpr size_t LAYER = 0;
+        typedef decltype(mSec.rootSteps()) ExtType;
-	    static constexpr size_t SIZE = OpClass::SIZE;
+	    
-	    typedef decltype(mSec.rootSteps()) ExtType;
+        AddExpr(T* dataPtr, const OpClass& sec);
-	    
+        AddExpr(const AddExpr& in) = default;
-	    //AssignmentExpr(OperationMaster& m, const OpClass& sec);
+        AddExpr(AddExpr&& in) = default;
-            AddExpr(T* dataPtr, const OpClass& sec);
+	    
-
+        inline void operator()(size_t start = 0); 
-	    AddExpr(const AddExpr& in) = default;
+        inline void operator()(size_t start, ExtType last);
-	    AddExpr(AddExpr&& in) = default;
+        auto rootSteps(std::intptr_t iPtrNum = 0) const -> ExtType;
-	    
-	    inline void operator()(size_t start = 0); 
-	    inline void operator()(size_t start, ExtType last);
-
-	    auto rootSteps(std::intptr_t iPtrNum = 0) const -> ExtType;
-	    
-	};
-
-	typedef T value_type;
-	//typedef OperationBase<T> OB;
-	typedef ContainerRange<T,Ranges...> CRange;
-	typedef ContainerIndex<T,typename Ranges::IndexType...> IndexType;
-	//typedef typename MultiRange<Ranges...>::IndexType IndexType;
-	
-	OperationMaster(MutableMultiArrayBase<T,Ranges...>& ma, const OpClass& second,
-			IndexType& index, bool doParallel = false);
-
-	OperationMaster(T* data, const OpClass& second,
-			IndexType& index, bool doParallel = false);
-
-	inline void set(size_t pos, T val) { AOp::sapply(mDataPtr[pos],val); }
-        
-	//inline void add(size_t pos, T val) { mDataPtr[pos] += val; }
-	inline T get(size_t pos) const;
-
-    private:
-
-	void performAssignment(std::intptr_t blockIndexNum);
-	OpClass const& mSecond;
-	//MutableMultiArrayBase<T,Ranges...>& mArrayRef;
-	T* mDataPtr;
-	IndexType mIndex;
-        bool mDoParallel;
     };
 
-    
     template <typename T, class... Ranges>
     class ConstOperationRoot : public OperationTemplate<T,ConstOperationRoot<T,Ranges...> >
     {
@@ -226,11 +185,9 @@ namespace MultiArrayTools
 	
     private:
 
-	//MultiArrayBase<T,Ranges...> const& mArrayRef;
 	const T* mDataPtr;
         const T* mOrigDataPtr;
 	IndexType mIndex;
-	//size_t mOff = 0;
 	std::shared_ptr<MultiArrayBase<T,Ranges...> > mMaPtr; // never remove this ptr, otherwise we lose temporary container instances!
     };
 
@@ -300,8 +257,6 @@ namespace MultiArrayTools
 		
     private:
 
-	//MultiArrayBase<T,Ranges...> const& mArrayRef;
-	//const T* mDataPtr;
         mutable IndexType mWorkIndex;
         std::shared_ptr<IndexType> mIndex;
     };
@@ -318,21 +273,27 @@ namespace MultiArrayTools
 
 	static constexpr size_t SIZE = 1;
         static constexpr bool CONT = true;
-	
+
+    private:
+
+	T* mDataPtr;
+        T* mOrigDataPtr;
+	IndexType mIndex;
+
+    public:
 	OperationRoot(MutableMultiArrayBase<T,Ranges...>& ma,
 		      const std::shared_ptr<typename Ranges::IndexType>&... indices);
 
 	OperationRoot(T* data, const IndexType& ind);
-        /*
-	template <class OpClass>
-	OperationMaster<T,SelfIdentity<T>,OpClass,Ranges...> operator=(const OpClass& in);
 
         template <class OpClass>
-	OperationMaster<T,plus<T>,OpClass,Ranges...> operator+=(const OpClass& in);
+        auto assign(const OpClass& in)
-
+            -> decltype(mIndex.ifor(1,in.loop(AssignmentExpr<T,OpClass>(mOrigDataPtr,in))));
-	OperationMaster<T,SelfIdentity<T>,OperationRoot,Ranges...> operator=(const OperationRoot& in);
-        */
 
+        template <class OpClass>
+        auto plus(const OpClass& in)
+            -> decltype(mIndex.ifor(1,in.loop(AddExpr<T,OpClass>(mOrigDataPtr,in))));
+        
         template <class OpClass>
         OperationRoot& operator=(const OpClass& in);
 
@@ -341,7 +302,7 @@ namespace MultiArrayTools
 
         OperationRoot& operator=(const OperationRoot& in);
         
-        OperationRoot& par();
+        ParallelOperationRoot<T,Ranges...> par();
         
 	template <class ET>
 	inline T& get(ET pos) const;
@@ -353,21 +314,69 @@ namespace MultiArrayTools
 
 	template <class Expr>
 	Expr loop(Expr exp) const;
-
+        
 	T* data() const;
 	
 	template <class... Indices>
 	auto sl(const std::shared_ptr<Indices>&... inds)
 	    -> Slice<T,typename Indices::RangeType...>;
 
+    };
+
+    template <typename T, class... Ranges>
+    class ParallelOperationRoot : public OperationTemplate<T,ParallelOperationRoot<T,Ranges...> >
+    {
+    public:
+
+	typedef T value_type;
+	typedef OperationBase<T,ParallelOperationRoot<T,Ranges...> > OT;
+	typedef ContainerRange<T,Ranges...> CRange;
+	typedef ContainerIndex<T,typename Ranges::IndexType...> IndexType;
+
+	static constexpr size_t SIZE = 1;
+        static constexpr bool CONT = true;
+
     private:
 
-	//MutableMultiArrayBase<T,Ranges...>& mArrayRef;
 	T* mDataPtr;
         T* mOrigDataPtr;
 	IndexType mIndex;
- 	//size_t mOff = 0;
+
-        bool mDoParallel = false;
+    public:
+	ParallelOperationRoot(MutableMultiArrayBase<T,Ranges...>& ma,
+                              const std::shared_ptr<typename Ranges::IndexType>&... indices);
+
+	ParallelOperationRoot(T* data, const IndexType& ind);
+
+        template <class OpClass>
+        auto assign(const OpClass& in)
+            -> decltype(mIndex.pifor(1,in.loop(AssignmentExpr<T,OpClass>(mOrigDataPtr,in))));
+
+        template <class OpClass>
+        auto plus(const OpClass& in)
+            -> decltype(mIndex.pifor(1,in.loop(AddExpr<T,OpClass>(mOrigDataPtr,in))));
+
+        template <class OpClass>
+        ParallelOperationRoot& operator=(const OpClass& in);
+
+        template <class OpClass>
+        ParallelOperationRoot& operator+=(const OpClass& in);
+
+        ParallelOperationRoot& operator=(const ParallelOperationRoot& in);
+        
+	template <class ET>
+	inline T& get(ET pos) const;
+
+	template <class ET>
+	inline ParallelOperationRoot& set(ET pos);
+
+	MExt<void> rootSteps(std::intptr_t iPtrNum = 0) const; // nullptr for simple usage with decltype
+
+	template <class Expr>
+	Expr loop(Expr exp) const;
+        
+	T* data() const;
+
     };
 
     template <typename T>

src/include/xfor/exttype.h

@@ -7,6 +7,41 @@
 namespace MultiArrayHelper
 {
 
+    template <size_t I>
+    struct Getter
+    {
+	template <class ExtType>
+	static inline size_t get(const ExtType& et)
+	{
+	    return Getter<I-1>::get(et.next());
+	}
+
+	template <class ExtType>
+	static inline auto getX(const ExtType& et)
+	    -> decltype(Getter<I-1>::getX(et.next()))
+	{
+	    return Getter<I-1>::getX(et.next());
+	}	
+    };
+
+    template <>
+    struct Getter<0>
+    {
+	template <class ExtType>
+	static inline size_t get(const ExtType& et)
+	{
+	    return et.get();
+	}
+
+	template <class ExtType>
+	static inline auto getX(const ExtType& et)
+	    -> ExtType
+	{
+	    return et;
+	}	
+    };
+
+    
     template <class X>
     class MExt
     {
@@ -38,7 +73,12 @@ namespace MultiArrayHelper
 
 	inline const size_t& val() const;
 	inline const X& next() const;
-	
+
+        template <size_t N>
+        inline auto nn() const
+            -> decltype(Getter<N>::getX(*this))
+        { return Getter<N>::getX(*this); }
+        
 	inline MExt operator+(const MExt& in) const;
 	inline MExt operator*(size_t in) const;
 
@@ -78,7 +118,12 @@ namespace MultiArrayHelper
 
 	inline const size_t& val() const;
 	inline size_t next() const { return 0; }
-	
+
+        template <size_t N>
+        inline auto nn() const
+            -> decltype(Getter<N>::getX(*this))
+        { return Getter<N>::getX(*this); }
+        
 	inline MExt operator+(const MExt& in) const;
 	inline MExt operator*(size_t in) const;
 
@@ -89,39 +134,6 @@ namespace MultiArrayHelper
     };    
 
 
-    template <size_t I>
-    struct Getter
-    {
-	template <class ExtType>
-	static inline size_t get(const ExtType& et)
-	{
-	    return Getter<I-1>::get(et.next());
-	}
-
-	template <class ExtType>
-	static inline auto getX(const ExtType& et)
-	    -> decltype(Getter<I-1>::getX(et.next()))
-	{
-	    return Getter<I-1>::getX(et.next());
-	}	
-    };
-
-    template <>
-    struct Getter<0>
-    {
-	template <class ExtType>
-	static inline size_t get(const ExtType& et)
-	{
-	    return et.get();
-	}
-
-	template <class ExtType>
-	static inline auto getX(const ExtType& et)
-	    -> ExtType
-	{
-	    return et;
-	}	
-    };
 
 } // end namespace MultiArrayHelper
 

src/include/xfor/xfor.h

@@ -402,8 +402,6 @@ namespace MultiArrayHelper
     {
 	typedef typename IndexClass::RangeType RangeType;
 	for(size_t pos = 0u; pos != ForBound<RangeType::ISSTATIC>::template bound<RangeType::SIZE>(mMax); ++pos){
-	//for(size_t pos = mSPos; pos != mMax; ++pos){
-	    //const size_t mnpos = PosForward<FT>::value(mlast, mMax, pos);
 	    const size_t mnpos = PosForward<FT>::valuex(mlast, mStep, pos);
 	    const ExtType npos = last + mExt*pos;
 	    mExpr(mnpos, npos);
@@ -416,8 +414,6 @@ namespace MultiArrayHelper
 	typedef typename IndexClass::RangeType RangeType;
 	const ExtType last;
 	for(size_t pos = 0u; pos != ForBound<RangeType::ISSTATIC>::template bound<RangeType::SIZE>(mMax); ++pos){
-	//for(size_t pos = mSPos; pos != mMax; ++pos){
-	    //const size_t mnpos = PosForward<FT>::value(mlast, mMax, pos);
 	    const size_t mnpos = PosForward<FT>::valuex(mlast, mStep, pos);
 	    const ExtType npos = last + mExt*pos;
 	    mExpr(mnpos, npos);
@@ -498,9 +494,9 @@ namespace MultiArrayHelper
         int pos = 0;
         size_t mnpos = 0;
         ExtType npos;
-        auto expr = mExpr;
+#pragma omp parallel shared(mExpr) private(pos,mnpos,npos)
-#pragma omp parallel shared(expr) private(pos,mnpos,npos)
         {
+            auto expr = mExpr;
 #pragma omp for nowait
             for(pos = 0; pos < static_cast<int>(ForBound<RangeType::ISSTATIC>::template bound<RangeType::SIZE>(mMax)); pos++){
                 mnpos = PosForward<ForType::DEFAULT>::valuex(mlast, mStep, pos);
@@ -519,9 +515,9 @@ namespace MultiArrayHelper
         int pos = 0;
         size_t mnpos = 0;
         ExtType npos;
-        auto expr = mExpr;
+#pragma omp parallel shared(mExpr) private(pos,mnpos,npos)
-#pragma omp parallel shared(expr) private(pos,mnpos,npos)
         {
+            auto expr = mExpr;
 #pragma omp for nowait
             for(pos = 0; pos < static_cast<int>(ForBound<RangeType::ISSTATIC>::template bound<RangeType::SIZE>(mMax)); pos++){
                 mnpos = PosForward<ForType::DEFAULT>::valuex(mlast, mStep, pos);