... remaining fixes: hybrid dynamic multi-threaded operations work

src/include/dynamic_operation.cc.h

@@ -41,10 +41,89 @@ namespace MultiArrayTools
         return std::make_shared<DynamicOperation<T,Operation>>(*this);
     }
 
+    template <typename T, class Operation, class... Ranges>
+    DynamicOuterOp<T,Operation,Ranges...>::DynamicOuterOp(const DynamicOuterOp& in) :
+	mThreadId(omp_get_thread_num()), mOp(in.mOp),
+	mIndices(in.mIndices),
+	mMa((mThreadId != in.mThreadId) ? std::make_shared<MultiArray<T,Ranges...>>(*in.mMa) : in.mMa),
+	mProto((mThreadId != in.mThreadId) ? OperationRoot<T,Ranges...>(*mMa,mIndices) : in.mProto),
+	mL((mThreadId != in.mThreadId) ?
+	   mkILoop(std::make_tuple(*mProto.mOp,mOp), mIndices,
+		   std::make_tuple(mMa),
+		   std::make_tuple(mProto.mOp->assign( mOp, mkMIndex(mIndices) )),
+		   std::array<size_t,1>({1}), std::array<size_t,1>({0})) :
+	   in.mL)
+	   
+    {*mMa = 0;}
+
+    template <typename T, class Operation, class... Ranges>
+    DynamicOuterOp<T,Operation,Ranges...>::DynamicOuterOp(DynamicOuterOp&& in) :
+	mThreadId(omp_get_thread_num()), mOp(in.mOp),
+	mIndices(in.mIndices),
+	mMa((mThreadId != in.mThreadId) ? std::make_shared<MultiArray<T,Ranges...>>(*in.mMa) : in.mMa),
+	mProto((mThreadId != in.mThreadId) ? OperationRoot<T,Ranges...>(*mMa,mIndices) : in.mProto),
+	mL((mThreadId != in.mThreadId) ?
+	   mkILoop(std::make_tuple(*mProto.mOp,mOp), mIndices,
+		   std::make_tuple(mMa),
+		   std::make_tuple(mProto.mOp->assign( mOp, mkMIndex(mIndices) )),
+		   std::array<size_t,1>({1}), std::array<size_t,1>({0})) :
+	   in.mL)
+	   
+    {*mMa = 0;}
+
+    template <typename T, class Operation, class... Ranges>
+    DynamicOuterOp<T,Operation,Ranges...>&
+    DynamicOuterOp<T,Operation,Ranges...>::operator=(const DynamicOuterOp& in)
+    {
+	mThreadId = omp_get_thread_num();
+	mOp = in.mOp;
+	mIndices = in.mIndices;
+	if(mThreadId != in.mThreadId){
+	    mMa = std::make_shared<MultiArray<T,Ranges...>>(in.mMa);
+	    mProto = OperationRoot<T,Ranges...>(*mMa,mIndices);
+	    mL = mkILoop(std::make_tuple(*mProto.mOp,mOp), mIndices,
+			 std::make_tuple(mMa),
+			 std::make_tuple(mProto.mOp->assign( mOp, mkMIndex(mIndices) )),
+			 std::array<size_t,1>({1}), std::array<size_t,1>({0}));
+	}
+	else {
+	    mMa = in.mMa;
+	    mProto = in.mProto;
+	    mL = in.mL;
+	}
+	*mMa = 0;
+	return *this;
+    }
+
+    template <typename T, class Operation, class... Ranges>
+    DynamicOuterOp<T,Operation,Ranges...>&
+    DynamicOuterOp<T,Operation,Ranges...>::operator=(DynamicOuterOp&& in)
+    {
+	mThreadId = omp_get_thread_num();
+	mOp = in.mOp;
+	mIndices = in.mIndices;
+	if(mThreadId != in.mThreadId){
+	    mMa = std::make_shared<MultiArray<T,Ranges...>>(in.mMa);
+	    mProto = OperationRoot<T,Ranges...>(*mMa,mIndices);
+	    mL = mkILoop(std::make_tuple(*mProto.mOp,mOp), mIndices,
+			 std::make_tuple(mMa),
+			 std::make_tuple(mProto.mOp->assign( mOp, mkMIndex(mIndices) )),
+			 std::array<size_t,1>({1}), std::array<size_t,1>({0}));
+	}
+	else {
+	    mMa = in.mMa;
+	    mProto = in.mProto;
+	    mL = in.mL;
+	}
+	*mMa = 0;
+	return *this;
+    }
+
     template <typename T, class Operation, class... Ranges>
     DynamicOuterOp<T,Operation,Ranges...>::DynamicOuterOp(const Operation& op,
                                                           const std::shared_ptr<typename Ranges::IndexType>&... inds)
-        : mOp(op),
+        : mThreadId(omp_get_thread_num()), mOp(op),
+	  mIndices(inds...),
 	  mMa(std::make_shared<MultiArray<T,Ranges...>>(mkArray<T>(inds->range()...))),
 	  mProto(OperationRoot<T,Ranges...>(*mMa,inds...)),
           mL(std::make_tuple(*mProto.mOp,mOp), std::make_tuple(inds...),

src/include/dynamic_operation.h

@@ -70,8 +70,10 @@ namespace MultiArrayTools
     class DynamicOuterOp : public DynamicOperationBase<OpH<OperationRoot<T,Ranges...>>>
     {
     private:
+	size_t mThreadId;
 	Operation mOp;
         //OperationRoot<T,Ranges...> mProto;
+	std::tuple<std::shared_ptr<typename Ranges::IndexType>...> mIndices;
 	std::shared_ptr<MultiArray<T,Ranges...>> mMa;
         OpH<OperationRoot<T,Ranges...>> mProto;
 
@@ -90,11 +92,11 @@ namespace MultiArrayTools
         typedef decltype(mL.rootSteps()) ET;
 	//typedef decltype(std::declval<Operation>().rootSteps()) ET;
 	
-	DynamicOuterOp() = default;
+	DynamicOuterOp() : mThreadId(omp_get_thread_num()) {}
-	DynamicOuterOp(const DynamicOuterOp& in) = default;
+	DynamicOuterOp(const DynamicOuterOp& in);
-	DynamicOuterOp(DynamicOuterOp&& in) = default;
+	DynamicOuterOp(DynamicOuterOp&& in);
-	DynamicOuterOp& operator=(const DynamicOuterOp& in) = default;
+	DynamicOuterOp& operator=(const DynamicOuterOp& in);
-	DynamicOuterOp& operator=(DynamicOuterOp&& in) = default;
+	DynamicOuterOp& operator=(DynamicOuterOp&& in);
 
 	DynamicOuterOp(const Operation& op, const std::shared_ptr<typename Ranges::IndexType>&... inds);
 	

src/include/helper_tools.cc.h

@@ -36,6 +36,41 @@ namespace MultiArrayTools
 	return std::dynamic_pointer_cast<MultiRange<RangeTypes...> >( mrf.create() );
     }
 
+    namespace
+    {
+	template <size_t N>
+	struct IndexToRangeTuple
+	{
+	    template <class... IndexTypes>
+	    static inline void set(std::tuple<std::shared_ptr<typename IndexTypes::RangeType>...>& out,
+		     const std::tuple<std::shared_ptr<IndexTypes>...>& indices)
+	    {
+		std::get<N>(out) = std::get<N>(indices)->range();
+		IndexToRangeTuple<N-1>::set(out,indices);
+	    }
+	};
+
+	template <>
+	struct IndexToRangeTuple<0>
+	{
+	    template <class... IndexTypes>
+	    static inline void set(std::tuple<std::shared_ptr<typename IndexTypes::RangeType>...>& out,
+		     const std::tuple<std::shared_ptr<IndexTypes>...>& indices)
+	    {
+		std::get<0>(out) = std::get<0>(indices)->range();
+	    }
+	};
+    }
+    
+    template <class... IndexTypes>
+    auto indexToRangeTuple(const std::tuple<std::shared_ptr<IndexTypes>...>& indices)
+	-> std::tuple<std::shared_ptr<typename IndexTypes::RangeType>...>
+    {
+	std::tuple<std::shared_ptr<typename IndexTypes::RangeType>...> out;
+	IndexToRangeTuple<sizeof...(IndexTypes)-1>::set(out, indices);
+	return out;
+    }
+
     template <class... IndexTypes>
     auto mkMIndex(std::shared_ptr<IndexTypes>... indices)
 	-> decltype( getIndex( mkMulti( indices->range()... ) ) )
@@ -45,6 +80,16 @@ namespace MultiArrayTools
 	return mi;
     }
 
+    template <class... IndexTypes>
+    auto mkMIndex(const std::tuple<std::shared_ptr<IndexTypes>...>& indices)
+	-> decltype( getIndex( mkMulti( indexToRangeTuple(indices) ) ) )
+    {
+	auto mi = getIndex( mkMulti( indexToRangeTuple(indices) ) );	    
+	(*mi)( indices );
+	return mi;
+    }
+
+
     template <class Index>
     auto mkIndexW(const std::shared_ptr<Index>& ind)
 	-> std::shared_ptr<IndexW>
@@ -102,7 +147,7 @@ namespace MultiArrayTools
 
     template <class... RangeTypes>
     auto mkMulti(std::tuple<std::shared_ptr<RangeTypes>...> rangesTuple)
-	-> MultiRange<RangeTypes...>
+	-> std::shared_ptr<MultiRange<RangeTypes...>>
     {
 	MultiRangeFactory<RangeTypes...> mrf( rangesTuple );
 	return std::dynamic_pointer_cast<MultiRange<RangeTypes...> >( mrf.create() );

src/include/helper_tools.h

@@ -56,18 +56,27 @@ namespace MultiArrayTools
     template <class Op, class MA, class... IndexTypes>
     auto mkMapI(const std::tuple<Op,MA>& f, std::shared_ptr<IndexTypes>... indices)
 	-> decltype( mkGenMapI<SpaceType::ANY>(f, indices... ) );
+
+    template <class... IndexTypes>
+    auto indexToRangeTuple(const std::tuple<std::shared_ptr<IndexTypes>...>& indices)
+	-> std::tuple<std::shared_ptr<typename IndexTypes::RangeType>...>;
+
+    template <class... RangeTypes>
+    auto mkMulti(std::tuple<std::shared_ptr<RangeTypes>...> rangesTuple)
+	-> std::shared_ptr<MultiRange<RangeTypes...>>;
     
     template <class... IndexTypes>
     auto mkMIndex(std::shared_ptr<IndexTypes>... indices)
 	-> decltype( getIndex( mkMulti( indices.range()... ) ) );
 
+    template <class... IndexTypes>
+    auto mkMIndex(const std::tuple<std::shared_ptr<IndexTypes>...>& indices)
+	-> decltype( getIndex( mkMulti( indexToRangeTuple(indices) ) ) );
+
     template <class Index>
     auto mkIndexW(const std::shared_ptr<Index>& ind)
 	-> std::shared_ptr<IndexW>;
 
-    template <class... RangeTypes>
-    auto mkMulti(std::tuple<std::shared_ptr<RangeTypes>...> rangesTuple)
-	-> MultiRange<RangeTypes...>;
 
     template <class RangeFactory>
     auto createExplicit(RangeFactory& rf)

src/include/multi_array_operation.cc.h

@@ -481,6 +481,18 @@ namespace MultiArrayTools
         mDataPtr = mOrigDataPtr + mIndex.pos();
     }
 
+    template <typename T, class... Ranges>
+    OperationRoot<T,Ranges...>::
+    OperationRoot(MutableMultiArrayBase<T,Ranges...>& ma,
+		  const std::tuple<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>
     OperationRoot<T,Ranges...>::
     OperationRoot(T* data, const IndexType& ind) :

src/include/multi_array_operation.h

@@ -368,6 +368,9 @@ namespace MultiArrayTools
 	OperationRoot(MutableMultiArrayBase<T,Ranges...>& ma,
 		      const std::shared_ptr<typename Ranges::IndexType>&... indices);
 
+	OperationRoot(MutableMultiArrayBase<T,Ranges...>& ma,
+		      const std::tuple<std::shared_ptr<typename Ranges::IndexType>...>& indices);
+
 	OperationRoot(T* data, const IndexType& ind);
 
         template <class OpClass>

src/include/ranges/multi_range.h

@@ -85,6 +85,7 @@ namespace MultiArrayTools
 	// NO foreign/external controll)
 	// Do NOT share index instances between two or more MultiIndex instances
 	MultiIndex& operator()(std::shared_ptr<Indices>&... indices);
+	MultiIndex& operator()(const std::tuple<std::shared_ptr<Indices>...>& indices);
 
 	// ==== >>>>> STATIC POLYMORPHISM <<<<< ====
 	
@@ -324,6 +325,14 @@ namespace MultiArrayTools
 	return *this;
     }
 
+    template <class... Indices>
+    MultiIndex<Indices...>& MultiIndex<Indices...>::operator()(const std::tuple<std::shared_ptr<Indices>...>& indices)
+    {
+	RPackNum<sizeof...(Indices)-1>::swapIndices(mIPack, indices);
+	RPackNum<sizeof...(Indices)-1>::setIndexPack(mIPack, IB::mPos);
+	return *this;
+    }
+
     template <class... Indices>
     IndexType MultiIndex<Indices...>::type() const
     {

src/tests/op4_unit_test.cc

@@ -137,14 +137,12 @@ namespace
         resx2(i1,di4) = mkDynOp(ma1(i1,di1) * ma2(i1,di2));
 	resx3(i1,di4) = mkDynOp(mkDynOp(ma1(i1,di1)) * mkDynOp(ma2(i1,di2)));
 
-	//auto op1 = mkDynOutOp((ma1(i1,di1) * ma2(i1,di2)), ci4_1, ci4_2);
+	auto op1x = mkDynOutOp((ma1(i1,di1) * ma2(i1,di2)), ci4_1, ci4_2);
-	//auto opr = resx4(i1,di4);
-	auto op1x = (ma1(i1,di1) * ma2(i1,di2));
 	auto opr = resx4(i1,di4);
 	
 	auto loop = mkPILoop
 	    ( [&op1x,&opr,&xx,this](){
-		auto op1 = mkDynOutOp(op1x, ci4_1, ci4_2);
+		auto op1 = op1x;
 		return mkGetExpr(op1,mkILoop(std::make_tuple(opr,op1,*op1.data()->mOp), std::make_tuple(ci4_1, ci4_2),
 					     std::make_tuple(xx),
 					     std::make_tuple(opr.assign( *op1.data()->mOp, mkMIndex(ci4_1, ci4_2) )),