.woodpecker.yml: try fix

.woodpecker.yml

@@ -3,13 +3,20 @@
 steps:
     build:
         image: chizeta/centos:8
+        pull: true
         commands:
             - source /opt/rh/gcc-toolset-9/enable
+            - source /etc/profile.d/modules.sh
+            - module use --append /usr/share/modulefiles
+            - module load mpi
+            - chown -R test:test .
+            - su test
             - mkdir build-gcc
             - cd build-gcc
             - cmake3 -DSCALAR_BUILD=on ..
             - make -j2
             - make test
+            - cat /woodpecker/src/git.f3l.de/chizeta/cnorxz/build-gcc/Testing/Temporary/LastTest.log
             - cd ..
             - mkdir build-clang
             - cd build-clang

CMakeLists.txt

@@ -123,6 +123,14 @@ else()
   endif()
 endif()
 
+# CHECK LIBRARIES : mpi
+
+if(DEFINED ENABLE_mpi)
+  set(ENABLE_mpi ${ENABLE_mpi} CACHE BOOL "enable mpi")
+else()
+  set(ENABLE_mpi TRUE CACHE BOOL "enable mpi")
+endif()
+
 # DEFINES
 
 add_definitions(-DVERSION="${VERSION}")

src/CMakeLists.txt

@@ -12,4 +12,8 @@ if(ENABLE_cereal)
   add_subdirectory(opt/cereal)
 endif()
 
+if(ENABLE_mpi)
+  add_subdirectory(opt/mpi)
+endif()
+
 install(DIRECTORY include/ DESTINATION ${INSTALL_PATH}/include/cnorxz)

src/bin/config.cc

@@ -25,6 +25,7 @@ auto configMap()
     m["--version"] = CNORXZ::Config::version;
     m["--commit"] = CNORXZ::Config::commit;
     m["--flags"] = CNORXZ::Config::flags;
+    m["--definitions"] = CNORXZ::Config::definitions;
     return m;
 }
 

src/include/array/array_base.cc.h

@@ -30,8 +30,8 @@ namespace CNORXZ
     template <typename I, typename M>
     const T& CArrayBase<T>::operator[](const IndexInterface<I,M>& i) const
     {
-	if(formatIsTrivial()){
-	    return data()[i.lex()];
+	if(formatIsTrivial() or not i.formatIsTrivial()){
+	    return data()[i.pos()];
 	}
 	else {
 	    auto ai = itLex(i);
@@ -130,19 +130,12 @@ namespace CNORXZ
 	else {
 	    if(i->formatIsTrivial()){
 		// try to apply container format.
-		// if the reformat changes the index type in any manner
-		// the format is not applicable:
-		if constexpr(std::is_same<decltype(i->reformat( Vector<SizeT>(), Vector<SizeT>() )),Index>::value){
-		    auto beg = begin();
-		    auto aformat = beg.deepFormat();
-		    auto amax = beg.deepMax();
-		    auto fi = i->reformat( aformat, amax );
-		    return coproot(*this, moveToPtr( fi ) );
-		}
-		else {
-		    this->checkFormatCompatibility(*i);
-		    return coproot(*this, i);
-		}
+		auto beg = begin();
+		auto aformat = beg.deepFormat();
+		auto amax = beg.deepMax();
+		//auto j = std::make_shared<Index>(*i); // use copy, otherwise would change i
+		i->reformat( aformat, amax );
+		return coproot(*this, i );
 	    }
 	    else {
 		// check if format is compatible
@@ -192,7 +185,7 @@ namespace CNORXZ
 	auto j = begin();
 	CXZ_ASSERT(acc.lmax().val() == j.lmax().val(),
 		   "got index of iteration space size = " << acc.lmax().val()
-		   << ", expected size = " << acc.lmax().val());
+		   << ", expected size = " << j.lmax().val());
 	Vector<SizeT> f1 = toVec(acc.deepFormat());
 	Vector<SizeT> f2 = j.deepFormat();
 	std::sort(f1.begin(),f1.end());
@@ -234,8 +227,8 @@ namespace CNORXZ
     template <typename I, typename M>
     T& ArrayBase<T>::operator[](const IndexInterface<I,M>& i)
     {
-	if(this->formatIsTrivial()){
-	    return data()[i.lex()];
+	if(this->formatIsTrivial() or not i.formatIsTrivial()){
+	    return data()[i.pos()];
 	}
 	else {
 	    auto ai = itLex(i);
@@ -322,19 +315,12 @@ namespace CNORXZ
 	else {
 	    if(i->formatIsTrivial()){
 		// try to apply container format.
-		// if the reformat changes the index type in any manner
-		// the format is not applicable:
-		if constexpr(std::is_same<decltype(i->reformat( Vector<SizeT>(), Vector<SizeT>() )),Index>::value){
-		    auto beg = begin();
-		    auto aformat = beg.deepFormat();
-		    auto amax = beg.deepMax();
-		    auto fi = i->reformat( aformat, amax );
-		    return oproot(*this, moveToPtr( fi ) );
-		}
-		else {
-		    this->checkFormatCompatibility(*i);
-		    return oproot(*this, i);
-		}
+		auto beg = begin();
+		auto aformat = beg.deepFormat();
+		auto amax = beg.deepMax();
+		//auto j = std::make_shared<Index>(*i); // use copy, otherwise would change i
+		i->reformat( aformat, amax );
+		return oproot(*this, i );
 	    }
 	    else {
 		// check if format is compatible

src/include/array/array_base.h

@@ -45,7 +45,12 @@ namespace CNORXZ
 	/** default destructor */
 	virtual ~CArrayBase() = default;
 
-	/** const data element access
+	/** Return unique pointer copy of itself. */
+	virtual Uptr<CArrayBase> copy() const = 0;
+	
+	/** const data element access.
+	    If the array's format is trivial or the index is already non-trivially formatted,
+	    this is equivalent to data()[i.pos()]. Otherwise, the array's format is applied.
 	    @tparam I index type
 	    @tparam M meta data type
 	    @param i index
@@ -100,7 +105,8 @@ namespace CNORXZ
 	Sptr<CArrayBase<T>> sl(const IndexInterface<I,M>& begin,
 			       const IndexInterface<I,M>& end) const;
 
-	/** create operation on this container
+	/** create operation on this container.
+	    Caution: might modify the index format.
 	    @tparam Index type of operation index
 	    @param i operation index
 	 */
@@ -143,6 +149,19 @@ namespace CNORXZ
 	/** check if container views the data, i.e. it does not own it */
 	virtual bool isView() const = 0;
 
+	/** Perform compatibility checks
+	    @tparam Acc index type or index pack type
+	    @param acc index or index pack.
+	 */
+	template <class Acc>
+	void checkFormatCompatibility(const Acc& acc) const;
+
+	/** check if format is trivial
+	    @return true if container is data owning array, else return
+	    result of the corresponding container index
+	 */
+	virtual bool formatIsTrivial() const = 0;
+
     protected:
 	RangePtr mRange; /**< the container range */
 
@@ -165,18 +184,6 @@ namespace CNORXZ
 	template <class Acc>
 	const_iterator itLexSave(const Acc& acc) const;
 
-	/** Perform compatibility checks
-	    @tparam Acc index type or index pack type
-	    @param acc index or index pack.
-	 */
-	template <class Acc>
-	void checkFormatCompatibility(const Acc& acc) const;
-
-	/** check if format is trivial
-	    @return true if container is data owning array, else return
-	    result of the corresponding container index
-	 */
-	virtual bool formatIsTrivial() const = 0;
     };
 
     /** ****

src/include/array/marray.cc.h

@@ -61,6 +61,12 @@ namespace CNORXZ
 	return *this;
     }
 
+    template <typename T>
+    Uptr<CArrayBase<T>> MArray<T>::copy() const
+    {
+	return std::make_unique<MArray<T>>(*this);
+    }
+
     template <typename T>
     const T* MArray<T>::data() const
     {

src/include/array/marray.h

@@ -38,16 +38,15 @@ namespace CNORXZ
 	MArray& init(const RangePtr& range);
 	MArray& extend(const RangePtr& range);
 
+	virtual Uptr<CArrayBase<T>> copy() const override;
 	virtual const T* data() const override;
 	virtual T* data() override;
 	virtual const_iterator cbegin() const override;
 	virtual const_iterator cend() const override;
 	virtual bool isView() const override;
+	virtual bool formatIsTrivial() const override final;
 
 	SERIALIZATION_FUNCTIONS;
-
-    protected:
-	virtual bool formatIsTrivial() const override final;
 	
     private:
 	Vector<T> mCont;

src/include/array/slice.cc.h

@@ -29,6 +29,12 @@ namespace CNORXZ
 	mOff(off)
     {}
 
+    template <typename T>
+    Uptr<CArrayBase<T>> CSlice<T>::copy() const
+    {
+	return std::make_unique<CSlice<T>>(*this);
+    }
+
     template <typename T>
     const T* CSlice<T>::data() const
     {
@@ -76,6 +82,12 @@ namespace CNORXZ
 	mOff(off)
     {}
 
+    template <typename T>
+    Uptr<CArrayBase<T>> Slice<T>::copy() const
+    {
+	return std::make_unique<Slice<T>>(*this);
+    }
+
     template <typename T>
     T* Slice<T>::data()
     {

src/include/array/slice.h

@@ -41,16 +41,16 @@ namespace CNORXZ
 	 */
 	CSlice(const RangePtr& range, const CArrayBase<T>* parent,
 	       const YFormat& blockSizes, SizeT off);
-	
+
+	virtual Uptr<CArrayBase<T>> copy() const override;
 	virtual const T* data() const override;
 	virtual const_iterator cbegin() const override;
 	virtual const_iterator cend() const override;
 	virtual bool isView() const override final;
-
-    protected:
-
 	virtual bool formatIsTrivial() const override final;
 	
+    protected:
+
 	const CArrayBase<T>* mCParent = nullptr; /**< pointer to the original container */
 	YFormat mBlockSizes; /**< the format */
 	SizeT mOff = 0; /** pointer offset w.r.t. the original pointer */
@@ -82,13 +82,12 @@ namespace CNORXZ
 	Slice(const RangePtr& range, ArrayBase<T>* parent,
 	      const YFormat& blockSizes, SizeT off);
 
+	virtual Uptr<CArrayBase<T>> copy() const override;
 	virtual const T* data() const override;
 	virtual T* data() override;
 	virtual const_iterator cbegin() const override;
 	virtual const_iterator cend() const override;
 	virtual bool isView() const override final;
-
-    protected:
 	virtual bool formatIsTrivial() const override final;
 
     private:

src/include/base/config.h

@@ -29,6 +29,9 @@ namespace CNORXZ
 
 	/** return compile flags */
 	String flags();
+
+	/** return compile definitions */
+	String definitions();
     }
 }
 

src/include/base/obj_handle.cc.h

@@ -41,7 +41,7 @@ namespace CNORXZ
     template <typename T>
     ObjHandle<T>& ObjHandle<T>::operator=(ObjHandle&& a)
     {
-	mC = a.mC;
+	mC = std::move(a.mC);
 	return *this;
     }
     
@@ -69,6 +69,18 @@ namespace CNORXZ
 	return &*mC;
     }
 
+    template <typename T>
+    T* ObjHandle<T>::get()
+    {
+	return &*mC;
+    }
+    
+    template <typename T>
+    const T* ObjHandle<T>::get() const
+    {
+	return &*mC;
+    }
+
 }
 
 #endif

src/include/base/obj_handle.h

@@ -77,6 +77,12 @@ namespace CNORXZ
 
 	/** get pointer to data (const) */
 	const T* operator->() const;
+
+	/** get object. */
+	T* get();
+
+	/** get object. */
+	const T* get() const;
     };
     
 }

src/include/base/types.h

@@ -171,12 +171,20 @@ namespace CNORXZ
     
     // definition: ranges/urange.h
     template <typename Meta>
-    class URange; // generic simple range (uni-dimensional)
+    class URange; // generic simple range (uni-dimensional, ordered)
 
     // definition: ranges/urange.h
     template <typename Meta>
     class UIndex;
 
+    // definition: ranges/wrange.h
+    template <typename Meta>
+    class WRange; // generic simple range (uni-dimensional)
+
+    // definition: ranges/wrange.h
+    template <typename Meta>
+    class WIndex;
+
     // definition: ranges/prange.h
     template <class Range>
     class PRange;

src/include/base/utils.h

@@ -36,78 +36,113 @@ namespace CNORXZ
     {
 	return Vector<T> { a };
     }
-    
+
+    template <typename T, typename U>
+    struct Concat
+    {
+	static constexpr Arr<T,2> cat(const T& a1, const U& a2)
+	{
+	    static_assert( std::is_same<T,U>::value, "types have to be vector, array or equal" );
+	    return Arr<T,2> { a1, a2 };
+	}
+    };
+
     template <typename T, SizeT N1, SizeT N2>
-    constexpr Arr<T,N1+N2> cat2(const Arr<T,N1>& a1, const Arr<T,N2>& a2)
+    struct Concat<Arr<T,N1>,Arr<T,N2>>
     {
-	return iter<0,N1+N2>
-	    ( [&](auto i) { if constexpr(i < N1) { return std::get<i>(a1); } else { return std::get<i-N1>(a2); } },
-	      [](const auto&... e) { return Arr<T,N1+N2> { e... }; } );
-    }
+	static constexpr Arr<T,N1+N2> cat(const Arr<T,N1>& a1, const Arr<T,N2>& a2)
+	{
+	    return iter<0,N1+N2>
+		( [&](auto i) { if constexpr(i < N1) { return std::get<i>(a1); } else { return std::get<i-N1>(a2); } },
+		  [](const auto&... e) { return Arr<T,N1+N2> { e... }; } );
+	}
+    };
 
     template <typename T, SizeT N1>
-    constexpr Arr<T,N1+1> cat2(const Arr<T,N1>& a1, const T& a2)
+    struct Concat<Arr<T,N1>,T>
     {
-	return iter<0,N1>
-	    ( [&](auto i) { return std::get<i>(a1); },
-	      [&](const auto&... e) { return Arr<T,N1+1> { e..., a2 }; } );
-    }
+	static constexpr Arr<T,N1+1> cat(const Arr<T,N1>& a1, const T& a2)
+	{
+	    return iter<0,N1>
+		( [&](auto i) { return std::get<i>(a1); },
+		  [&](const auto&... e) { return Arr<T,N1+1> { e..., a2 }; } );
+	}
+    };
 
     template <typename T, SizeT N1>
-    constexpr Arr<T,N1+1> cat2(const T& a1, const Arr<T,N1>& a2)
+    struct Concat<T,Arr<T,N1>>
     {
-	return iter<0,N1>
-	    ( [&](auto i) { return std::get<i>(a2); },
-	      [&](const auto&... e) { return Arr<T,N1+1> { a1, e... }; } );
-    }
+	static constexpr Arr<T,N1+1> cat(const T& a1, const Arr<T,N1>& a2)
+	{
+	    return iter<0,N1>
+		( [&](auto i) { return std::get<i>(a2); },
+		  [&](const auto&... e) { return Arr<T,N1+1> { a1, e... }; } );
+	}
+    };
 
-    template <typename T>
-    constexpr Arr<T,2> cat2(const T& a1, const T& a2)
-    {
-	return Arr<T,2> { a1, a2 };
-    }
 
     template <typename T, SizeT N2>
-    Vector<T> cat2(const Vector<T>& a1, const Arr<T,N2>& a2)
+    struct Concat<Vector<T>,Arr<T,N2>>
     {
-	Vector<T> o(a1.size()+N2);
-	std::copy(a1.begin(), a1.end(), o.begin());
-	std::copy(a2.begin(), a2.end(), o.begin()+a1.size());
-	return o;
-    }
+	static Vector<T> cat(const Vector<T>& a1, const Arr<T,N2>& a2)
+	{
+	    Vector<T> o(a1.size()+N2);
+	    std::copy(a1.begin(), a1.end(), o.begin());
+	    std::copy(a2.begin(), a2.end(), o.begin()+a1.size());
+	    return o;
+	}
+    };
 
     template <typename T, SizeT N1>
-    Vector<T> cat2(const Arr<T,N1>& a1, const Vector<T>& a2)
+    struct Concat<Arr<T,N1>,Vector<T>>
     {
-	Vector<T> o(N1+a2.size());
-	std::copy(a1.begin(), a1.end(), o.begin());
-	std::copy(a2.begin(), a2.end(), o.begin()+N1);
-	return o;
-    }
+	static Vector<T> cat(const Arr<T,N1>& a1, const Vector<T>& a2)
+	{
+	    Vector<T> o(N1+a2.size());
+	    std::copy(a1.begin(), a1.end(), o.begin());
+	    std::copy(a2.begin(), a2.end(), o.begin()+N1);
+	    return o;
+	}
+    };
 
     template <typename T>
-    Vector<T> cat2(const Vector<T>& a1, const Vector<T>& a2)
+    struct Concat<Vector<T>,Vector<T>>
     {
-	Vector<T> o(a1.size()+a2.size());
-	std::copy(a1.begin(), a1.end(), o.begin());
-	std::copy(a2.begin(), a2.end(), o.begin()+a1.size());
-	return o;
-    }
+	static Vector<T> cat(const Vector<T>& a1, const Vector<T>& a2)
+	{
+	    Vector<T> o(a1.size()+a2.size());
+	    std::copy(a1.begin(), a1.end(), o.begin());
+	    std::copy(a2.begin(), a2.end(), o.begin()+a1.size());
+	    return o;
+	}
+    };
 
     template <typename T>
-    Vector<T> cat2(const Vector<T>& a1, const T& a2)
+    struct Concat<Vector<T>,T>
     {
-	Vector<T> o(a1);
-	o.push_back(a2);
-	return o;
-    }
+	static Vector<T> cat(const Vector<T>& a1, const T& a2)
+	{
+	    Vector<T> o(a1);
+	    o.push_back(a2);
+	    return o;
+	}
+    };
 
-    template <typename T, SizeT N1>
-    Vector<T> cat2(const T& a1, const Vector<T>& a2)
+    template <typename T>
+    struct Concat<T,Vector<T>>
     {
-	Vector<T> o { a1 };
-	o.insert(o.end(), a2.begin(), a2.end());
-	return o;
+	static Vector<T> cat(const T& a1, const Vector<T>& a2)
+	{
+	    Vector<T> o { a1 };
+	    o.insert(o.end(), a2.begin(), a2.end());
+	    return o;
+	}
+    };
+
+    template <typename T1, typename T2>
+    decltype(auto) cat2(const T1& a1, const T2& a2)
+    {
+	return Concat<T1,T2>::cat(a1, a2);
     }
 
     template <typename T1, typename T2, typename... Ts>

src/include/operation/op_types.cc.h

@@ -379,6 +379,12 @@ namespace CNORXZ
 	return OpRoot<T,IndexT>(a, ind);
     }
 
+    template <typename T, class IndexT>
+    constexpr decltype(auto) oproot(T* a, const Sptr<IndexT>& ind)
+    {
+	return OpRoot<T,IndexT>(a, ind);
+    }
+
     /*=================+
      |   Operation     |
      +=================*/
@@ -509,7 +515,35 @@ namespace CNORXZ
 	return i->xpr(i);
     }
 
-    
+    /*==============+
+     |    PosOp     |
+     +==============*/
+
+    template <class Index>
+    constexpr PosOp<Index>::PosOp(const Sptr<Index>& i) :
+	mI(i)
+    {}
+
+    template <class Index>
+    template <class PosT>
+    constexpr decltype(auto) PosOp<Index>::operator()(const PosT& pos) const
+    {
+	return static_cast<SizeT>(pos);
+    }
+
+    template <class Index>
+    constexpr decltype(auto) PosOp<Index>::operator()() const
+    {
+	return static_cast<SizeT>(0);
+    }
+
+    template <class Index>
+    template <SizeT I>
+    constexpr decltype(auto) PosOp<Index>::rootSteps(const IndexId<I>& id) const
+    {
+	return mI->stepSize(id);
+    }
+
 }
 
 #endif

src/include/operation/op_types.h

@@ -217,6 +217,9 @@ namespace CNORXZ
     template <typename T, class IndexT>
     constexpr decltype(auto) oproot(ArrayBase<T>& a, const Sptr<IndexT>& ind);
 
+    template <typename T, class IndexT>
+    constexpr decltype(auto) oproot(T* a, const Sptr<IndexT>& ind);
+
     template <class F, class... Ops>
     class Operation : public OpInterface<Operation<F,Ops...>>
     {
@@ -290,9 +293,37 @@ namespace CNORXZ
     template <class F, class Op, class IndexT>
     constexpr decltype(auto) contraction(F&& f, Op&& op, const Sptr<IndexT>& i);
 
+    // redundant? see xpr()
     template <class IndexT>
     constexpr decltype(auto) indexOp(const Sptr<IndexT>& i);
 
+    template <class Index>
+    class PosOp : public COpInterface<PosOp<Index>>
+    {
+    public:
+	typedef COpInterface<PosOp<Index>> OI;
+
+	constexpr PosOp() = default;
+	constexpr PosOp(const Sptr<Index>& i);
+
+	template <class PosT>
+	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;
+
+    private:
+	Sptr<Index> mI;
+    };
+
+    template <class Index>
+    constexpr decltype(auto) posop(const Sptr<Index>& i)
+    {
+	return PosOp<Index>(i);
+    }
+
 }
 
 #endif

src/include/operation/op_utility.cc.h

@@ -52,7 +52,7 @@ namespace CNORXZ
     inline auto pos_unpack_args(const F& f, const PosT& pos, const Tuple<Ops...>& args)
     {
 	static_assert(is_pos_type<PosT>::value, "got non-pos-type");
-	static_assert((is_operation<Ops>::value and ...), "got non-operation type");
+	static_assert((is_xpr<Ops>::value and ...), "got non-operation type");
 	typedef std::make_index_sequence<sizeof...(Ops)> Idxs;
 	typedef std::index_sequence<op_size<Ops>::value...> OpSizes;
 	return pos_unpack_args_i(f, pos, args, OpSizes{}, Idxs{});

src/include/ranges/crange.h

@@ -185,7 +185,8 @@ namespace CNORXZ
 	 */
 	SizeT get(SizeT pos) const;
 
-	/** return position for given meta data
+	/** return position for given meta data.
+	    Returns size() if metaPos is out of scope.
 	    @param metaPos meta data, size type
 	 */
 	SizeT getMeta(SizeT metaPos) const;

src/include/ranges/mrange.cc.h

@@ -60,8 +60,12 @@ namespace CNORXZ
     template <class FormatT, class... Indices>
     inline void GMIndex<FormatT,Indices...>::mkPos()
     {
+	bool outOfScope = false;
 	mLex = iter<0,NI>
-	    ([&](auto i) { return mIPack[i]->lex() * mLexFormat[i].val(); },
+	    ([&](auto i) {
+		outOfScope |= mIPack[i]->lex() >= mIPack[i]->lmax().val();
+		return mIPack[i]->lex() * mLexFormat[i].val();
+	    },
 	     [](const auto&... e) { return (e + ...); });
 	if constexpr(not std::is_same<FormatT,None>::value){
 	    IB::mPos = iter<0,NI>
@@ -71,6 +75,10 @@ namespace CNORXZ
 	else {
 	    IB::mPos = mLex;
 	}
+	if(outOfScope){
+	    IB::mPos = pmax().val();
+	    mLex = lmax().val();
+	}
     }
 
     template <class FormatT, class... Indices>
@@ -79,7 +87,7 @@ namespace CNORXZ
     {
 	return gmformat
 	    ( iter<Is,NI>
-	      ( [&](auto i) { return ipack[i]->pmax(); },
+	      ( [&](auto i) { return ipack[i]->lmax(); },
 		[](const auto&... as) { return (as * ...); } )...,
 	      SPos<1>() );
     }
@@ -404,9 +412,14 @@ namespace CNORXZ
     template <SizeT I>
     decltype(auto) GMIndex<FormatT,Indices...>::stepSize(const IndexId<I>& id) const
     {
+	// TODO: also for YRange!!!
+	const auto own = [&]() {
+	    if constexpr(I != 0){ return SPos<0> {}; }
+	    else { return UPos(id == this->id() ? 1 : 0); }
+	};
 	return iter<0,NI>
 	    ( [&](auto i) { return mIPack[i]->stepSize(id) * format()[i]; },
-	      [](const auto&... ss) { return ( ss + ... ); });
+	      [](const auto&... ss) { return ( ss + ... ); }) + own();
     }
 
     template <class FormatT, class... Indices>
@@ -415,11 +428,16 @@ namespace CNORXZ
 	const String blim = "(";
 	const String elim = ")";
 	const String dlim = ",";
-	return iter<1,NI>
-	    ( [&](auto i) { return mIPack[i]->stringMeta(); },
-	      [&](const auto&... xs) {
-		  return blim + mIPack[CSizeT<0>{}]->stringMeta() + ( (dlim + xs) + ... ) + elim;
-	      } );
+	if constexpr(NI == 1){
+	    return blim + mIPack[CSizeT<0>{}]->stringMeta() + elim;
+	}
+	else {
+	    return iter<1,NI>
+		( [&](auto i) { return mIPack[i]->stringMeta(); },
+		  [&](const auto&... xs) {
+		      return blim + mIPack[CSizeT<0>{}]->stringMeta() + ( (dlim + xs) + ... ) + elim;
+		  } );
+	}
     }
     
     template <class FormatT, class... Indices>
@@ -433,9 +451,7 @@ namespace CNORXZ
     GMIndex<FormatT,Indices...>& GMIndex<FormatT,Indices...>::at(const MetaType& metaPos)
     {
 	iter<0,NI>( [&](auto i) { mIPack[i]->at( std::get<i>(metaPos) ); }, NoF {} );
-	IB::mPos = iter<0,NI>
-	    ( [&](auto i) { return mIPack[i]->pos()*format()[i].val(); },
-	      [](const auto&... xs) { return (xs + ...); });
+	mkPos();
 	return *this;
     }
 
@@ -625,6 +641,12 @@ namespace CNORXZ
 	return MIndex<Indices...>(pack);
     }
 
+    template <class... Indices>
+    constexpr decltype(auto) mindexPtr(const Sptr<Indices>&... is)
+    {
+	return std::make_shared<MIndex<Indices...>>(is...);
+    }
+
     template <class... Indices>
     constexpr decltype(auto) mindexPtr(const SPack<Indices...>& pack)
     {
@@ -637,6 +659,12 @@ namespace CNORXZ
 	return std::make_shared<GMIndex<FormatT,Indices...>>(bs, pack);
     }    
     
+    template <class... Indices>
+    constexpr decltype(auto) gmindexPtr(const SPack<Indices...>& pack)
+    {
+	return std::make_shared<GMIndex<MFormat<sizeof...(Indices)>,Indices...>>(pack);
+    }
+
     template <class I1, class FormatT, class... Indices>
     decltype(auto) operator*(const Sptr<GMIndex<FormatT,Indices...>>& a, const Sptr<I1>& b)
     {
@@ -801,6 +829,12 @@ namespace CNORXZ
     {
 	return MRangeFactory<Ranges...>(std::make_tuple(rs...)).create();
     }
+
+    template <class... Ranges>
+    Sptr<MRange<Ranges...>> xplMrange(const Sptr<Ranges>&... rs)
+    {
+	return std::dynamic_pointer_cast<MRange<Ranges...>>( mrange( rs... ) );
+    }
     
     template <class... Ranges>
     Sptr<MRange<Ranges...>> RangeCast<MRange<Ranges...>>::func(const RangePtr& r)

src/include/ranges/mrange.h

@@ -163,7 +163,7 @@ namespace CNORXZ
 
 	// replace sub-index instances; only use if you know what you are doing!
 	/** Replace sub-index instances and update index position correspondingly.
-	    @param new index instances.
+	    @param mi New index instances.
 	 */
 	GMIndex& operator()(const Sptr<MIndex<Indices...>>& mi);
 
@@ -286,6 +286,12 @@ namespace CNORXZ
     template <class... Indices>
     constexpr decltype(auto) mindex(const SPack<Indices...>& pack);
 
+    /** Create pointer to MIndex from index pack.
+	@param is Input index pointers.
+     */
+    template <class... Indices>
+    constexpr decltype(auto) mindexPtr(const Sptr<Indices>&... is);
+
     /** Create pointer to MIndex from index pack.
 	@param pack Pack of input indices.
      */
@@ -299,6 +305,12 @@ namespace CNORXZ
     template <class FormatT, class... Indices>
     constexpr decltype(auto) gmindexPtr(const FormatT& bs, const SPack<Indices...>& pack);
 
+    /** Create pointer to GMIndex from index pack (default format).
+	@param pack Pack of input indices.
+     */
+    template <class... Indices>
+    constexpr decltype(auto) gmindexPtr(const SPack<Indices...>& pack);
+
     /** Specialization for index multiplication with GMIndex on the l.h.s.
 	@param a First index of type GMIndex.
 	@param b Second index of arbitrary type.
@@ -361,8 +373,8 @@ namespace CNORXZ
 	/** Get sub-ranges. */
 	decltype(auto) space() const;
 
-	/** Get sub-range.
-	    @param pos Position of the sub-range.
+	/** Get meta data for given lexicographic position.
+	    @param pos Lexicographic position.
 	 */
 	const MetaType get(SizeT pos) const;
 
@@ -399,6 +411,12 @@ namespace CNORXZ
     template <class... Ranges>
     RangePtr mrange(const Sptr<Ranges>&... rs);
 
+    /** Create MRange pointer.
+	@param rs Pointer to sub-ranges. 
+     */
+    template <class... Ranges>
+    Sptr<MRange<Ranges...>> xplMrange(const Sptr<Ranges>&... rs);
+
     /** ****
 	Specialization of RangeCast for MRange.
 	@see RangeCast.

src/include/ranges/prange.cc.h

@@ -367,9 +367,8 @@ namespace CNORXZ
     {
 	Vector<MetaType> meta(this->size());
 	auto i = mRange->begin();
-	SizeT j = 0;
-	for(const auto& p: mParts){
-	    meta[j++] = *(i = p);
+	for(SizeT j = 0; j != meta.size(); ++j){
+	    meta[j] = *(i = mParts[j]);
 	}
 	return URangeFactory<MetaType>( meta ).create();
     }

src/include/ranges/ranges.cc.h

@@ -14,6 +14,7 @@
 #include "mrange.cc.h"
 #include "xindex.cc.h"
 #include "urange.cc.h"
+#include "wrange.cc.h"
 #include "srange.cc.h"
 #include "crange.cc.h"
 #include "prange.cc.h"

src/include/ranges/ranges.h

@@ -16,6 +16,7 @@
 #include "xindex.h"
 #include "yrange.h"
 #include "urange.h"
+#include "wrange.h"
 #include "srange.h"
 #include "crange.h"
 #include "prange.h"

src/include/ranges/urange.cc.h

@@ -299,8 +299,12 @@ namespace CNORXZ
 	auto b = mSpace.begin();
 	auto e = mSpace.end();
 	auto i = std::lower_bound(b, e, meta, std::less<MetaT>());
-	CXZ_ASSERT(i != e, "element with meta data = " << toString(meta) << " not in range"); // check this first, otherwise the next test may potentially result in a seg fault!
-	CXZ_ASSERT(*i == meta, "element with meta data = " << toString(meta) << " not in range");
+	if(i == e){
+	    return size();
+	}
+	if(*i != meta){
+	    return size();
+	}
 	return i - b;
     }
 
@@ -359,7 +363,21 @@ namespace CNORXZ
 	template <typename T>
 	static inline Sptr<URange<MetaT>> transform(const RangePtr& r)
 	{
-	    if(r->type() == typeid(URange<T>)){
+	    if(r->type() == typeid(CRange)){
+		auto rr = std::dynamic_pointer_cast<CRange>(r);
+		Vector<MetaT> v(rr->size());
+		std::transform(rr->begin(), rr->end(), v.begin(),
+			       [](const T& x) { return static_cast<MetaT>(x); } );
+		return std::dynamic_pointer_cast<URange<MetaT>>
+		    ( URangeFactory<MetaT>(std::move(v)).create() );
+	    }
+	    else if(r->type() == typeid(PRange<CRange>)){
+		return transform<T>( std::dynamic_pointer_cast<PRange<CRange>>(r)->derive() );
+	    }
+	    else if(r->type() == typeid(PRange<URange<T>>)){
+		return transform<T>( std::dynamic_pointer_cast<PRange<URange<T>>>(r)->derive() );
+	    }
+	    else if(r->type() == typeid(URange<T>)){
 		auto rr = std::dynamic_pointer_cast<URange<T>>(r);
 		Vector<MetaT> v(rr->size());
 		std::transform(rr->begin(), rr->end(), v.begin(),
@@ -377,6 +395,7 @@ namespace CNORXZ
 	    static_assert(std::is_fundamental<MetaT>::value, "got non-fundamental type");
 	    CXZ_ASSERT(r->dim() == 1, "range cast into URange<Int>: source range must have dim = 1, got " << r->dim());
 	    Sptr<URange<MetaT>> o = nullptr;
+	    // TODO: cast from CRange!!!
 	    o = transform<SizeT>(r); if(o) return o;
 	    o = transform<Int>(r); if(o) return o;
 	    o = transform<LInt>(r); if(o) return o;
@@ -387,9 +406,6 @@ namespace CNORXZ
 	}
     };
 
-    template <typename T> struct is_vector : std::false_type {};
-    template <typename T, typename A> struct is_vector<std::vector<T,A>> : std::true_type {};
-
     template <typename U>
     struct URangeCast<Vector<U>>
     {

src/include/ranges/urange.h

@@ -162,7 +162,10 @@ namespace CNORXZ
 
     /** ****
 	Uni-(1-)dimensional range with non-trivial meta data space
-	i.e. the parameter space can be arbitrary.
+	i.e. the parameter space can be arbitrary. However, the parameter space
+	is assumed to be ordered, i.e. the i-1-th element is assumed to be smaller
+	(according to std::less) than the i-th element.
+	If you require arbitrary ordering, use WRange / WIndex instead!
 	@tparam MetaT Meta data type.
      */
     template <typename MetaT>
@@ -194,6 +197,7 @@ namespace CNORXZ
 	const MetaType* get() const;
 
 	/** Get range position for given meta data.
+	    Returns size() if metaPos is not part of the range.
 	    @param metaPos Meta data.
 	    @return Position of the given meta data if it is contained by the range.
 	 */
@@ -240,6 +244,11 @@ namespace CNORXZ
     {
 	static constexpr bool value = true;
     };
+
+    // TODO: move to utils file!!!
+    template <typename T> struct is_vector : std::false_type {};
+    template <typename T, typename A> struct is_vector<std::vector<T,A>> : std::true_type {};
+
 }
 
 #endif

src/include/ranges/wrange.cc.h

@@ -0,0 +1,515 @@
+// -*- C++ -*-
+/**
+   
+   @file include/ranges/wrange.cc.h
+   @brief WRange, WRangeFactory and WIndex implementations.
+
+   Copyright (c) 2024 Christian Zimmermann. All rights reserved.
+   Mail: chizeta@f3l.de
+   
+**/
+
+#ifndef __cxz_wrange_cc_h__
+#define __cxz_wrange_cc_h__
+
+#include <functional>
+#include <algorithm>
+
+#include "wrange.h"
+#include "prange.h"
+#include "index_mul.h"
+#include "xpr/for.h"
+#include "operation/op_types.h"
+#include "operation/op_types.cc.h"
+
+namespace CNORXZ
+{
+    /*===============+
+     |    WIndex     |	     
+     +===============*/
+
+    template <typename MetaT>
+    WIndex<MetaT>::WIndex(const RangePtr& range, SizeT pos) :
+	IndexInterface<WIndex<MetaT>,MetaT>(pos),
+	mRangePtr(rangeCast<RangeType>(range)),
+	mMetaPtr(mRangePtr->get()),
+	mSpaceToPos(mRangePtr->spaceToPos()),
+	mPosToSpace(mRangePtr->posToSpace())
+    {}
+
+    template <typename MetaT>
+    WIndex<MetaT>& WIndex<MetaT>::operator=(size_t lexpos)
+    {
+	IB::mPos = lexpos;
+	return *this;
+    }
+
+    template <typename MetaT>
+    WIndex<MetaT>& WIndex<MetaT>::operator++()
+    {
+	++IB::mPos;
+	return *this;
+    }
+
+    template <typename MetaT>
+    WIndex<MetaT>& WIndex<MetaT>::operator--()
+    {
+	--IB::mPos;
+	return *this;
+    }
+
+    template <typename MetaT>
+    WIndex<MetaT> WIndex<MetaT>::operator+(Int n) const
+    {
+	return WIndex(mRangePtr, IB::mPos + n);
+    }
+    
+    template <typename MetaT>
+    WIndex<MetaT> WIndex<MetaT>::operator-(Int n) const
+    {
+	return WIndex(mRangePtr, IB::mPos - n);
+    }
+
+    template <typename MetaT>
+    SizeT WIndex<MetaT>::operator-(const WIndex& i) const
+    {
+	return lex() - i.lex();
+    }
+    
+    template <typename MetaT>
+    WIndex<MetaT>& WIndex<MetaT>::operator+=(Int n)
+    {
+	IB::mPos += n;
+	return *this;
+    }
+    
+    template <typename MetaT>
+    WIndex<MetaT>& WIndex<MetaT>::operator-=(Int n)
+    {
+	IB::mPos -= n;
+	return *this;
+    }
+
+    template <typename MetaT>
+    SizeT WIndex<MetaT>::lex() const
+    {
+	return IB::mPos;
+    }
+    
+    template <typename MetaT>
+    UPos WIndex<MetaT>::pmax() const
+    {
+	return UPos(mRangePtr->size());
+    }
+
+    template <typename MetaT>
+    UPos WIndex<MetaT>::lmax() const
+    {
+	return UPos(mRangePtr->size());
+    }
+
+    template <typename MetaT>
+    IndexId<0> WIndex<MetaT>::id() const
+    {
+	return IndexId<0>(this->ptrId());
+    }
+
+    template <typename MetaT>
+    const MetaT& WIndex<MetaT>::operator*() const
+    {
+	return meta();
+    }
+    
+    template <typename MetaT>
+    String WIndex<MetaT>::stringMeta() const
+    {
+	return toString(this->meta());
+    }
+    
+    template <typename MetaT>
+    const MetaT& WIndex<MetaT>::meta() const
+    {
+	return mMetaPtr[mPosToSpace[IB::mPos]];
+    }
+    
+    template <typename MetaT>
+    WIndex<MetaT>& WIndex<MetaT>::at(const MetaT& metaPos)
+    {
+	(*this) = mRangePtr->getMeta(metaPos);
+	return *this;
+    }
+
+    template <typename MetaT>
+    decltype(auto) WIndex<MetaT>::xpr(const Sptr<WIndex<MetaT>>& _this) const
+    {
+	CXZ_ERROR("not implemented"); // TODO: Operation root with reordered meta data access!!!
+	return coproot(mMetaPtr,_this);
+    }
+
+    template <typename MetaT>
+    RangePtr WIndex<MetaT>::prange(const WIndex<MetaT>& last) const
+    {
+	CXZ_ASSERT(last >= *this, "got last index position (" << last.lex()
+		   << ") smaller than begin index position (" << lex() << ")");
+	const SizeT beginPos = lex();
+	Vector<SizeT> parts(last.lex()-beginPos+1);
+	for(auto i = *this; i != last+1; ++i){
+	    parts[i.lex()-beginPos] = i.lex();
+	}
+	return CNORXZ::prange(mRangePtr, parts);
+    }
+
+    template <typename MetaT>
+    SizeT WIndex<MetaT>::deepFormat() const
+    {
+	return 1;
+    }
+
+    template <typename MetaT>
+    SizeT WIndex<MetaT>::deepMax() const
+    {
+	return lmax().val();
+    }
+
+    template <typename MetaT>
+    WIndex<MetaT>& WIndex<MetaT>::reformat(const Vector<SizeT>& f, const Vector<SizeT>& s)
+    {
+	CXZ_ASSERT(f[0]*s[0] == lmax().val(), "got wrong extension: " << f[0]*s[0]
+		   << " vs " << lmax().val());
+	CXZ_ASSERT(CNORXZ::formatIsTrivial(f,s), "format is not trivial: f = " << toString(f)
+		   << ", s = " << toString(s));
+	return *this;
+    }
+
+    template <typename MetaT>
+    size_t WIndex<MetaT>::dim() const // = 1
+    {
+	return 1;
+    }
+
+    template <typename MetaT>
+    Sptr<WRange<MetaT>> WIndex<MetaT>::range() const
+    {
+	return mRangePtr;
+    }
+
+    template <typename MetaT>
+    template <SizeT I>
+    decltype(auto) WIndex<MetaT>::stepSize(const IndexId<I>& id) const
+    {
+	if constexpr(I != 0){
+	    return SPos<0>();
+	}
+	else {
+	    return UPos(id == this->id() ? 1 : 0);
+	}
+    }
+
+    template <typename MetaT>
+    template <class Xpr, class F>
+    decltype(auto) WIndex<MetaT>::ifor(const Xpr& xpr, F&& f) const
+    {
+	return For<0,Xpr,F>(this->pmax().val(), this->id(), xpr, std::forward<F>(f));
+    }
+
+    template <typename MetaT>
+    bool WIndex<MetaT>::formatIsTrivial() const
+    {
+	return true;
+    }
+    
+    template <typename MetaT, class I1>
+    decltype(auto) operator*(const Sptr<WIndex<MetaT>>& a, const Sptr<I1>& b)
+    {
+	return iptrMul(a, b);
+    }
+
+    /*====================+
+     |   WRangeFactory    |
+     +====================*/
+
+    template <typename MetaT>
+    WRangeFactory<MetaT>::WRangeFactory(const Vector<MetaT>& space) :
+	mSpace(space) {}
+
+    template <typename MetaT>
+    WRangeFactory<MetaT>::WRangeFactory(Vector<MetaT>&& space) :
+	mSpace(space) {}
+
+    template <typename MetaT>
+    WRangeFactory<MetaT>::WRangeFactory(const Vector<MetaT>& space, const RangePtr& ref) :
+	mSpace(space), mRef(ref) {}
+
+    template <typename MetaT>
+    WRangeFactory<MetaT>::WRangeFactory(Vector<MetaT>&& space, const RangePtr& ref) :
+	mSpace(space), mRef(ref) {}
+
+    template <typename MetaT>
+    void WRangeFactory<MetaT>::make()
+    {
+	const auto& info = typeid(WRange<MetaT>);
+	if(mRef != nullptr) {
+	    mProd = this->fromCreated(info, {mRef->id()});
+	}
+	if(mProd == nullptr){
+	    RangePtr key = mProd = std::shared_ptr<WRange<MetaT>>
+		( new WRange<MetaT>( std::move(mSpace) ) );
+	    if(mRef != nullptr) { key = mRef; }
+	    this->addToCreated(info, { key->id() }, mProd);
+	}
+    }
+    
+    /*=============+
+     |   WRange    |
+     +=============*/
+
+    template <typename MetaT>
+    void WRange<MetaT>::setupSpace(const Vector<MetaT>& space)
+    {
+	mSpace = space;
+	const SizeT s = mSpace.size();
+	mSpaceToPos.resize(s);
+	mPosToSpace.resize(s);
+	std::sort(mSpace.begin(), mSpace.end(), std::less<MetaT>());
+	auto itdupl = std::adjacent_find(mSpace.begin(), mSpace.end());
+	CXZ_ASSERT(itdupl == mSpace.end(), "found duplicate: " << *itdupl);
+	for(SizeT lx = 0; lx != s; ++lx){
+	    for(SizeT sp = 0; sp != s; ++sp){
+		if(space[lx] == mSpace[sp]){
+		    mPosToSpace[lx] = sp;
+		    mSpaceToPos[sp] = lx;
+		    break;
+		}
+	    }
+	}
+    }
+    
+    template <typename MetaT>
+    WRange<MetaT>::WRange(const Vector<MetaT>& space) :
+	RangeInterface<WRange<MetaT>>()
+    {
+	setupSpace(space);
+    }
+
+    template <typename MetaT>
+    WRange<MetaT>::WRange(Vector<MetaT>&& space) :
+	RangeInterface<WRange<MetaT>>()
+    {
+	setupSpace(space);
+    }
+
+        
+    template <typename MetaT>
+    const MetaT& WRange<MetaT>::get(SizeT pos) const
+    {
+	return mSpace[mPosToSpace[pos]];
+    }
+
+    template <typename MetaT>
+    const MetaT* WRange<MetaT>::get() const
+    {
+	return mSpace.data();
+    }
+
+    template <typename MetaT>
+    const SizeT* WRange<MetaT>::spaceToPos() const
+    {
+	return mSpaceToPos.data();
+    }
+
+    template <typename MetaT>
+    const SizeT* WRange<MetaT>::posToSpace() const
+    {
+	return mPosToSpace.data();
+    }
+
+    template <typename MetaT>
+    SizeT WRange<MetaT>::getMeta(const MetaT& meta) const
+    {
+	auto b = mSpace.begin();
+	auto e = mSpace.end();
+	auto i = std::lower_bound(b, e, meta, std::less<MetaT>());
+	if(i == e){
+	    return size();
+	}
+	if(*i != meta){
+	    return size();
+	}
+	return mSpaceToPos[i - b];
+    }
+
+    template <typename MetaT>
+    SizeT WRange<MetaT>::size() const
+    {
+	return mSpace.size();
+    }
+
+    template <typename MetaT>
+    SizeT WRange<MetaT>::dim() const
+    {
+	return 1;
+    }
+
+    template <typename MetaT>
+    String WRange<MetaT>::stringMeta(SizeT pos) const
+    {
+	return toString(this->get(pos));
+    }
+    
+    template <typename MetaT>
+    const TypeInfo& WRange<MetaT>::type() const
+    {
+	return typeid(WRange<MetaT>);
+    }
+    
+    template <typename MetaT>
+    const TypeInfo& WRange<MetaT>::metaType() const
+    {
+	return typeid(MetaT);
+    }
+
+    template <typename MetaT>
+    RangePtr WRange<MetaT>::extend(const RangePtr& r) const
+    {
+	auto rx = rangeCast<WRange<MetaT>>(r);
+	Vector<MetaT> nspace(mSpace.size()+rx->size());
+	SizeT i = 0;
+	for(; i != mSpace.size(); ++i){
+	    nspace[i] = mSpace[mPosToSpace[i]];
+	}
+	for(SizeT j = 0; j != rx->mSpace.size(); ++j, ++i){
+	    nspace[i] = rx->mSpace[rx->mPosToSpace[j]];
+	}
+	return WRangeFactory<MetaT>( nspace ).create();
+    }
+
+    template <typename MetaT>
+    Vector<Uuid> WRange<MetaT>::key() const
+    {
+	return Vector<Uuid> { this->id() };
+    }
+
+    /*=================+
+     |   Range Casts   |
+     +=================*/
+
+    template <typename MetaT>
+    struct WRangeCast
+    {
+	template <typename T>
+	static inline Sptr<WRange<MetaT>> transform(const RangePtr& r)
+	{
+	    if(r->type() == typeid(PRange<WRange<T>>)){
+		return transform<T>( std::dynamic_pointer_cast<PRange<WRange<T>>>(r)->derive() );
+	    }
+	    else if(r->type() == typeid(WRange<T>)){
+		auto rr = std::dynamic_pointer_cast<WRange<T>>(r);
+		Vector<MetaT> v(rr->size());
+		std::transform(rr->begin(), rr->end(), v.begin(),
+			       [](const T& x) { return static_cast<MetaT>(x); } );
+		return std::dynamic_pointer_cast<WRange<MetaT>>
+		    ( WRangeFactory<MetaT>(std::move(v)).create() );
+	    }
+	    else {
+		return nullptr;
+	    }
+	}
+	
+	static inline Sptr<WRange<MetaT>> cast(const RangePtr& r)
+	{
+	    static_assert(std::is_fundamental<MetaT>::value, "got non-fundamental type");
+	    CXZ_ASSERT(r->dim() == 1, "range cast into WRange<Int>: source range must have dim = 1, got " << r->dim());
+	    Sptr<WRange<MetaT>> o = nullptr;
+	    // TODO: cast from CRange!!!
+	    o = transform<SizeT>(r); if(o) return o;
+	    o = transform<Int>(r); if(o) return o;
+	    o = transform<LInt>(r); if(o) return o;
+	    o = transform<Double>(r); if(o) return o;
+	    // else general transform using DType (better than nothing), to be implemented!!!
+	    CXZ_ERROR("no range cast available for input range '" << r->type().name() << "'");
+	    return nullptr;
+	}
+    };
+
+    template <typename U>
+    struct WRangeCast<Vector<U>>
+    {
+	template <typename T>
+	static inline Sptr<WRange<Vector<U>>> transform(const RangePtr& r)
+	{
+	    if(r->type() == typeid(WRange<T>)){
+		auto rr = std::dynamic_pointer_cast<WRange<T>>(r);
+		Vector<Vector<U>> v(rr->size());
+		std::transform(rr->begin(), rr->end(), v.begin(),
+			       [](const T& x) { return Vector<U> { static_cast<U>(x) }; } );
+		return std::dynamic_pointer_cast<WRange<Vector<U>>>
+		    ( WRangeFactory<Vector<U>>(std::move(v)).create() );
+	    }
+	    else {
+		return nullptr;
+	    }
+	}
+
+	template <typename T, SizeT N>
+	static inline Sptr<WRange<Vector<U>>> atransform(const RangePtr& r)
+	{
+	    if(r->type() == typeid(WRange<Arr<T,N>>)){
+		auto rr = std::dynamic_pointer_cast<WRange<Arr<T,N>>>(r);
+		Vector<Vector<U>> v(rr->size());
+		std::transform(rr->begin(), rr->end(), v.begin(),
+			       [](const Arr<T,N>& x) {
+				   return iter<0,N>( [&](auto i) { return static_cast<U>(x[i]); },
+						     [](const auto&... e) { return Vector<U>{ e... }; });
+			       } );
+		return std::dynamic_pointer_cast<WRange<Vector<U>>>
+		    ( WRangeFactory<Vector<U>>(std::move(v)).create() );
+	    }
+	    else {
+		return nullptr;
+	    }
+	}
+
+	static inline Sptr<WRange<Vector<U>>> cast(const RangePtr& r)
+	{
+	    Sptr<WRange<Vector<U>>> o = nullptr;
+	    if constexpr(std::is_fundamental<U>::value){
+		o = transform<SizeT>(r); if(o) return o;
+		o = transform<Int>(r); if(o) return o;
+		o = transform<LInt>(r); if(o) return o;
+		o = transform<Double>(r); if(o) return o;
+		o = atransform<SizeT,2>(r); if(o) return o;
+		o = atransform<Int,2>(r); if(o) return o;
+		o = atransform<LInt,2>(r); if(o) return o;
+		o = atransform<Double,2>(r); if(o) return o;
+		o = atransform<SizeT,3>(r); if(o) return o;
+		o = atransform<Int,3>(r); if(o) return o;
+		o = atransform<LInt,3>(r); if(o) return o;
+		o = atransform<Double,3>(r); if(o) return o;
+	    }
+	    // else general transform using DType (better than nothing), to be implemented!!!
+	    CXZ_ERROR("no range cast available for input range '" << r->type().name() << "'");
+	}
+    };
+
+    
+    template <typename MetaT>
+    Sptr<WRange<MetaT>> RangeCast<WRange<MetaT>>::func(const RangePtr& r)
+    {
+	if constexpr(std::is_fundamental<MetaT>::value or is_vector<MetaT>::value){
+	    return WRangeCast<MetaT>::cast(r);
+	}
+	else {
+	    CXZ_ERROR("no range cast available for input range '" << r->type().name() << "'");
+	    return nullptr;
+	}
+    }
+
+    template <typename MetaT>
+    RangePtr wrange(const Vector<MetaT>& space)
+    {
+	return WRangeFactory<MetaT>(space).create();
+    }
+    
+}
+
+#endif

src/include/ranges/wrange.h

@@ -0,0 +1,266 @@
+// -*- C++ -*-
+/**
+   
+   @file include/ranges/wrange.h
+   @brief WRange, WRangeFactory and WIndex declaration.
+
+   Copyright (c) 2024 Christian Zimmermann. All rights reserved.
+   Mail: chizeta@f3l.de
+   
+**/
+
+#ifndef __cxz_wrange_h__
+#define __cxz_wrange_h__
+
+#include "base/base.h"
+#include "ranges/index_base.h"
+#include "ranges/range_base.h"
+#include "xpr/xpr.h"
+
+#include <iterator>
+
+namespace CNORXZ
+{
+    
+    /** ****
+	Specific index for WRange.
+	@tparam MetaT Meta data type.
+     */
+    template <typename MetaT>
+    class WIndex : public IndexInterface<WIndex<MetaT>,MetaT>
+    {
+    public:
+
+	typedef IndexInterface<WIndex<MetaT>,MetaT> IB;
+	typedef WRange<MetaT> RangeType;
+	typedef MetaT MetaType;
+
+	INDEX_RANDOM_ACCESS_ITERATOR_DEFS(MetaType);
+	DEFAULT_MEMBERS(WIndex); /**< default constructors and assignments */
+	
+	/** Construct index from range and position.
+	    @param range Range to iterate over.
+	    @param pos lexicographic position.
+	 */
+	WIndex(const RangePtr& range, SizeT pos = 0);
+
+	/** @copydoc IndexInterface::operator=(SizeT) */
+	WIndex& operator=(SizeT lexpos);
+
+	/** @copydoc IndexInterface::operator++() */
+	WIndex& operator++();
+
+	/** @copydoc IndexInterface::operator--() */
+	WIndex& operator--();
+
+	/** @copydoc IndexInterface::operator+() */
+	WIndex operator+(Int n) const;
+
+	/** @copydoc IndexInterface::operator-() */
+	WIndex operator-(Int n) const;
+
+	/** @copydoc IndexInterface::operator-(WIndex) */
+	SizeT operator-(const WIndex& i) const;
+
+	/** @copydoc IndexInterface::operator+=() */
+	WIndex& operator+=(Int n);
+
+	/** @copydoc IndexInterface::operator-=() */
+	WIndex& operator-=(Int n);
+
+	/** @copydoc IndexInterface::lex() */
+	SizeT lex() const;
+
+	/** @copydoc IndexInterface::pmax() */
+	UPos pmax() const;
+
+	/** @copydoc IndexInterface::lmax() */
+	UPos lmax() const;
+
+	/** @copydoc IndexInterface::id() */
+	IndexId<0> id() const;
+	
+	/** @copydoc IndexInterface::operator*() */
+	const MetaT& operator*() const;
+	
+	/** @copydoc IndexInterface::dim() */
+	SizeT dim() const; // = 1
+
+	/** @copydoc IndexInterface::range() */
+	Sptr<RangeType> range() const;
+
+	/** @copydoc IndexInterface::stepSize() */
+	template <SizeT I>
+	decltype(auto) stepSize(const IndexId<I>& id) const;
+
+	/** @copydoc IndexInterface::stringMeta() */
+	String stringMeta() const;
+
+	/** @copydoc IndexInterface::meta() */
+	const MetaT& meta() const;
+
+	/** @copydoc IndexInterface::at() */
+	WIndex& at(const MetaT& metaPos);
+
+	/** @copydoc IndexInterface::prange() */
+	RangePtr prange(const WIndex<MetaType>& last) const;
+
+	/** @copydoc IndexInterface::deepFormat() */
+	SizeT deepFormat() const;
+
+	/** @copydoc IndexInterface::deepMax() */
+	SizeT deepMax() const;
+
+	/** @copydoc IndexInterface::reformat() */
+	WIndex& reformat(const Vector<SizeT>& f, const Vector<SizeT>& s);
+
+	/** @copydoc IndexInterface::ifor() */
+	template <class Xpr, class F>
+	decltype(auto) ifor(const Xpr& xpr, F&& f) const;
+
+	/** @copydoc IndexInterface::formatIsTrivial() */
+	bool formatIsTrivial() const;
+
+	/** @copydoc IndexInterface::xpr() */
+	decltype(auto) xpr(const Sptr<WIndex<MetaType>>& _this) const;
+
+    private:
+	Sptr<RangeType> mRangePtr;
+	const MetaT* mMetaPtr;
+	const SizeT* mSpaceToPos;
+	const SizeT* mPosToSpace;
+    };
+
+    template <typename MetaT>
+    void swap(WIndex<MetaT>& a, WIndex<MetaT>& b) { a.swap(b); } 
+    
+    /** Make index pack of a WIndex and another index.
+	@param a pointer to WIndex.
+	@param b pointer to another index.
+     */
+    template <typename MetaT, class I1>
+    decltype(auto) operator*(const Sptr<WIndex<MetaT>>& a, const Sptr<I1>& b);
+    
+    /** ****
+	Specific factory for WRange.
+	@tparam MetaT Meta data type.
+     */
+    template <typename MetaT>
+    class WRangeFactory : public RangeFactoryBase
+    {
+    public:
+	WRangeFactory(const Vector<MetaT>& space);
+        WRangeFactory(Vector<MetaT>&& space);
+	WRangeFactory(const Vector<MetaT>& space, const RangePtr& ref);
+        WRangeFactory(Vector<MetaT>&& space, const RangePtr& ref);
+	
+    private:
+	WRangeFactory() = default;
+	virtual void make() override final;
+
+	Vector<MetaT> mSpace;
+	RangePtr mRef;
+    };
+
+    /** ****
+	Uni-(1-)dimensional range with non-trivial meta data space
+	i.e. the parameter space can be arbitrary. In contrast to URange
+	the meta data can have arbitrary ordering.
+	Note: Internally, the meta data is still ordered for the sake of performance during look-up.
+	The externally visible non-ordering is realized through maps between the internal (meta space)
+	ordering, and the externally visible lexicographic ordering.
+	@tparam MetaT Meta data type.
+     */
+    template <typename MetaT>
+    class WRange : public RangeInterface<WRange<MetaT>>
+    {
+    public:
+	typedef RangeBase RB;
+	typedef WIndex<MetaT> IndexType;
+	typedef MetaT MetaType;
+	
+	friend WRangeFactory<MetaType>;
+
+	virtual SizeT size() const override final;
+	virtual SizeT dim() const override final;
+	virtual String stringMeta(SizeT pos) const override final;
+	virtual const TypeInfo& type() const override final;
+	virtual const TypeInfo& metaType() const override final;
+	virtual RangePtr extend(const RangePtr& r) const override final;
+
+	/** Get meta data at given range position.
+	    @param pos Integer indicating requested position.
+	    @return Meta data at given postion.
+	 */
+	const MetaType& get(SizeT pos) const;
+
+	/** Get range position for given meta data.
+	    Returns size() if metaPos is not part of the range.
+	    @param metaPos Meta data.
+	    @return Position of the given meta data if it is contained by the range.
+	 */
+	SizeT getMeta(const MetaType& metaPos) const;
+
+	/** Get meta data array.
+	    @return Pointer to first element (memory ordering) of the underlying meta data array.
+	 */
+	const MetaType* get() const;
+
+	/** Get map from meta space memory position to lexicographic position.
+	    @return Pointer to first element of the map.
+	 */ 
+	const SizeT* spaceToPos() const;
+
+	/** Get map from lexicographic position to meta space memory position.
+	    @return Pointer to first element of the map.
+	 */ 
+	const SizeT* posToSpace() const;
+
+    private:
+
+	void setupSpace(const Vector<MetaT>& space);
+	
+	WRange() = default;
+	WRange(const WRange& in) = delete;
+	WRange(const Vector<MetaType>& space);
+        WRange(Vector<MetaType>&& space);
+
+	Vector<MetaType> mSpace;
+	Vector<SizeT> mSpaceToPos;
+	Vector<SizeT> mPosToSpace;
+
+	virtual Vector<Uuid> key() const override final;
+	
+	SERIALIZATION_FUNCTIONS_NOPUB;
+    };
+
+    /** ***
+	Specialize RangeCast for casts to WRange
+	@see RangeCast
+     */    
+    template <typename MetaType>
+    struct RangeCast<WRange<MetaType>>
+    {
+	/** cast the range */
+	static Sptr<WRange<MetaType>> func(const RangePtr& r);
+    };
+
+    /** Create an WRange, calls WRangeFactory.
+	@param space Meta data space to create an WRange on.
+	@return Created range.
+     */
+    template <typename MetaT>
+    RangePtr wrange(const Vector<MetaT>& space);
+
+    /** ***
+	WIndex can be used as expression
+	@see index_expression_exists
+     */    
+    template <typename MetaT>
+    struct index_expression_exists<WIndex<MetaT>>
+    {
+	static constexpr bool value = true;
+    };
+}
+
+#endif

src/include/ranges/xindex.cc.h

@@ -142,7 +142,7 @@ namespace CNORXZ
     template <class Index, typename Meta>
     UPos XIndex<Index,Meta>::stepSize(const IndexId<0>& id) const
     {
-	return mI->stepSize(id);
+	return UPos(mI->stepSize(id).val());
     }
 
     template <class Index, typename Meta>

src/include/ranges/yrange.h

@@ -311,6 +311,7 @@ namespace CNORXZ
     public:
 	typedef RangeBase RB;
 	typedef YIndex IndexType;
+	typedef typename IndexType::MetaType MetaType;
 
 	friend YRangeFactory;
 
@@ -343,6 +344,12 @@ namespace CNORXZ
      */
     RangePtr yrange(const Vector<RangePtr>& rs);
 
+    /** Create YRange from sub-ranges.
+	@param rs Vector of pointers to the sub-ranges used by the YRange.
+	@return A shared pointer to the created YRange.
+     */
+    Sptr<YRange> xplYrange(const Vector<RangePtr>& rs);
+
     /** ****
 	Specialize RangeCast for casts to YRange.
 	@see RangeCast

src/include/xpr/for.cc.h

@@ -150,7 +150,9 @@ namespace CNORXZ
 	constexpr SPos<I> i;
 	const auto pos = last + mExt( i );
 	if constexpr(I < N-1){
-	    return mF(mXpr(pos),exec<I+1>(last));
+	    auto x = mXpr(pos);
+	    mF(x,exec<I+1>(last));
+	    return x;
 	}
 	else {
 	    return mXpr(pos);
@@ -164,7 +166,9 @@ namespace CNORXZ
 	constexpr SPos<I> i;
 	const auto pos = mExt( i );
 	if constexpr(I < N-1){
-	    return mF(mXpr(pos),exec<I+1>());
+	    auto x = mXpr(pos);
+	    mF(x,exec<I+1>());
+	    return x;
 	}
 	else {
 	    return mXpr(pos);

src/include/xpr/map_xpr.cc.h

@@ -0,0 +1,113 @@
+// -*- C++ -*-
+/**
+   
+   @file include/xpr/map_xpr.cc.h
+   @brief Map xpression template implementations.
+
+   Copyright (c) 2024 Christian Zimmermann. All rights reserved.
+   Mail: chizeta@f3l.de
+   
+**/
+
+#ifndef __cxz_map_xpr_cc_h__
+#define __cxz_map_xpr_cc_h__
+
+#include "map_xpr.h"
+
+namespace CNORXZ
+{
+    template <class TarIndex, class SrcIndex, class F>
+    void
+    MapSetup<TarIndex,SrcIndex,F>::setup(const Sptr<TarIndex>& ti, const Sptr<SrcIndex>& si,
+					 const F& f, const Sptr<Vector<SizeT>>& m)
+    {
+	auto six = *si;
+	auto sie = si->range()->end();
+	auto tix = *ti;
+	for(six = 0; six != sie; ++six){
+	    tix.at( f(*six) );
+	    (*m)[six.pos()] = tix.pos();
+	}
+    }
+
+    template <class TarIndex, class SrcIndex, class F>
+    Sptr<Vector<SizeT>>
+    MapSetup<TarIndex,SrcIndex,F>::setup(const Sptr<TarIndex>& ti, const Sptr<SrcIndex>& si,
+					 const F& f)
+    {
+	auto o = std::make_shared<Vector<SizeT>>(si->lmax().val());
+	setupMap(ti,si,f,o);
+	return o;
+    }
+
+    template <class TarIndex, class SrcIndex, class F>
+    void setupMap(const Sptr<TarIndex>& ti, const Sptr<SrcIndex>& si,
+		  const F& f, const Sptr<Vector<SizeT>>& m)
+    {
+	MapSetup<TarIndex,SrcIndex,F>::setup(ti,si,f,m);
+    }
+
+    template <class TarIndex, class SrcIndex, class F>
+    Sptr<Vector<SizeT>> setupMap(const Sptr<TarIndex>& ti, const Sptr<SrcIndex>& si,
+				 const F& f)
+    {
+	return MapSetup<TarIndex,SrcIndex,F>::setup(ti,si,f);
+    }
+
+    template <class TarIndex, class SrcIndex, class Xpr>
+    template <class F>
+    MapXpr<TarIndex,SrcIndex,Xpr>::MapXpr(const Sptr<TarIndex>& ti, const Sptr<SrcIndex>& si,
+					  const F& f, Xpr&& xpr) :
+	mTi(ti), mSi(si),
+	mMap(nullptr),
+	mXpr(std::forward<Xpr>(xpr)),
+	mExt(mkFPos( mXpr.rootSteps(mTi->id()), mMap->data() ))
+    {
+	mMap = setupMap(ti,si,f);
+    }
+
+    template <class TarIndex, class SrcIndex, class Xpr>
+    MapXpr<TarIndex,SrcIndex,Xpr>::MapXpr(const Sptr<TarIndex>& ti, const Sptr<SrcIndex>& si,
+					  const Sptr<Vector<SizeT>>& m, Xpr&& xpr) :
+	mTi(ti), mSi(si),
+	mMap(m), mXpr(std::forward<Xpr>(xpr)),
+	mExt(mkFPos( mXpr.rootSteps(mTi->id()), mMap->data() ))
+    {}
+	    
+
+    template <class TarIndex, class SrcIndex, class Xpr>
+    template <class PosT>
+    decltype(auto) MapXpr<TarIndex,SrcIndex,Xpr>::operator()(const PosT& last) const
+    {
+	return mXpr( last.next() + mExt( last ) );
+    }
+
+    template <class TarIndex, class SrcIndex, class Xpr>
+    decltype(auto) MapXpr<TarIndex,SrcIndex,Xpr>::operator()() const
+    {
+	return mXpr( mExt( UPos(0) ) );
+    }
+
+    template <class TarIndex, class SrcIndex, class Xpr>
+    template <SizeT I>
+    decltype(auto) MapXpr<TarIndex,SrcIndex,Xpr>::rootSteps(const IndexId<I>& id) const
+    {
+	return mSi->stepSize(id) << mXpr.rootSteps(id);
+    }
+
+    template <class TarIndex, class SrcIndex, class F, class Xpr>
+    decltype(auto) mapXpr(const Sptr<TarIndex>& ti, const Sptr<SrcIndex>& si, const F& f, Xpr&& xpr)
+    {
+	return MapXpr<TarIndex,SrcIndex,Xpr>(ti,si,f,std::forward<Xpr>(xpr));
+    }
+
+    template <class TarIndex, class SrcIndex, class Xpr>
+    decltype(auto) mapXpr(const Sptr<TarIndex>& ti, const Sptr<SrcIndex>& si,
+			  const Sptr<Vector<SizeT>>& m, Xpr&& xpr)
+    {
+	return MapXpr<TarIndex,SrcIndex,Xpr>(ti,si,m,std::forward<Xpr>(xpr));
+    }
+
+}
+
+#endif

src/include/xpr/map_xpr.h

@@ -0,0 +1,80 @@
+// -*- C++ -*-
+/**
+
+   @file include/xpr/map_xpr.h
+   @brief Map expression declarations.
+
+   Copyright (c) 2024 Christian Zimmermann. All rights reserved.
+   Mail: chizeta@f3l.de
+   
+ **/
+
+#ifndef __cxz_map_xpr_h__
+#define __cxz_map_xpr_h__
+
+#include "xpr_base.h"
+
+namespace CNORXZ
+{
+
+    template <class TarIndex, class SrcIndex, class F>
+    struct MapSetup
+    {
+	static void setup(const Sptr<TarIndex>& ti, const Sptr<SrcIndex>& si,
+			  const F& f, const Sptr<Vector<SizeT>>& m);
+	
+	static Sptr<Vector<SizeT>> setup(const Sptr<TarIndex>& ti, const Sptr<SrcIndex>& si,
+					 const F& f);
+    };
+
+    template <class TarIndex, class SrcIndex, class F>
+    void setupMap(const Sptr<TarIndex>& ti, const Sptr<SrcIndex>& si,
+		  const F& f, const Sptr<Vector<SizeT>>& m);
+
+    template <class TarIndex, class SrcIndex, class F>
+    Sptr<Vector<SizeT>> setupMap(const Sptr<TarIndex>& ti, const Sptr<SrcIndex>& si,
+				 const F& f);
+    
+    template <class TarIndex, class SrcIndex, class Xpr>
+    class MapXpr : public XprInterface<MapXpr<TarIndex,SrcIndex,Xpr>>
+    {
+    private:
+	Sptr<TarIndex> mTi;
+	Sptr<SrcIndex> mSi;
+	Sptr<Vector<SizeT>> mMap;
+	Xpr mXpr;
+	typedef decltype(mkFPos( mXpr.rootSteps(mTi->id()), mMap->data() )) Ext;
+	Ext mExt;
+    
+    public:
+	
+	MapXpr() = default;
+
+	// src local
+	template <class F>
+	MapXpr(const Sptr<TarIndex>& ti, const Sptr<SrcIndex>& si, const F& f, Xpr&& xpr);
+
+	MapXpr(const Sptr<TarIndex>& ti, const Sptr<SrcIndex>& si,
+	       const Sptr<Vector<SizeT>>& m, Xpr&& xpr);
+	    
+
+	template <class PosT>
+	decltype(auto) operator()(const PosT& last) const;
+
+	decltype(auto) operator()() const;
+
+	template <SizeT I>
+	decltype(auto) rootSteps(const IndexId<I>& id) const;
+    };
+
+    template <class TarIndex, class SrcIndex, class F, class Xpr>
+    decltype(auto) mapXpr(const Sptr<TarIndex>& ti, const Sptr<SrcIndex>& si, const F& f, Xpr&& xpr);
+
+    template <class TarIndex, class SrcIndex, class Xpr>
+    decltype(auto) mapXpr(const Sptr<TarIndex>& ti, const Sptr<SrcIndex>& si,
+			  const Sptr<Vector<SizeT>>& m, Xpr&& xpr);
+    
+   
+}
+
+#endif

src/include/xpr/pos_type.cc.h

@@ -187,7 +187,8 @@ namespace CNORXZ
      |   FPos   |
      +==========*/
 
-    inline FPos::FPos(SizeT ext, const SizeT* map) : mExt(ext), mMap(map) {}
+    inline FPos::FPos(SizeT ext, const SizeT* map, SizeT max, SizeT max2) :
+	mExt(ext), mMap(map), mMax(max), mMax2(max2) {}
     
     constexpr SizeT FPos::size() const
     {
@@ -234,6 +235,36 @@ namespace CNORXZ
 	return val();
     }
 
+    constexpr SizeT FPos::max() const
+    {
+	return mMax;
+    }
+
+    template <class PosT>
+    struct MkFPos
+    {
+	static constexpr decltype(auto) mk(const PosT& pos, const SizeT* map)
+	{
+	    return FPos(pos.val(), map);
+	}
+    };
+
+    template <class BPosT, class NPosT>
+    struct MkFPos<MPos<BPosT,NPosT>>
+    {
+	static constexpr decltype(auto) mk(const MPos<BPosT,NPosT>& pos, const SizeT* map)
+	{
+	    return mkMPos( MkFPos<BPosT>::mk( pos, map ), MkFPos<NPosT>::mk( pos.next(), map ) );
+	}
+    };
+
+    template <class PosT>
+    constexpr decltype(auto) mkFPos(const PosT& pos, const SizeT* map)
+    {
+	return MkFPos<PosT>::mk(pos, map);
+    }
+
+    
     /*===========+
      |   SFPos   |
      +===========*/
@@ -320,6 +351,12 @@ namespace CNORXZ
 	return val();
     }
 
+    template <SizeT N, SizeT... Ms>
+    constexpr SizeT SFPos<N,Ms...>::max() const
+    {
+	return sizeof...(Ms); // CHECK!!!
+    }
+    
     /*==========+
      |   MPos   |
      +==========*/
@@ -411,6 +448,12 @@ namespace CNORXZ
 	return extend(a);
     }
     
+    template <class BPosT, class NPosT>
+    constexpr decltype(auto) mkMPos(const BPosT& bpos, const NPosT& npos)
+    {
+	return MPos<BPosT,NPosT>(bpos, npos);
+    }
+
     /*==========+
      |   DPos   |
      +==========*/

src/include/xpr/pos_type.h

@@ -93,11 +93,13 @@ namespace CNORXZ
     private:
 	SizeT mExt = 0;
 	const SizeT* mMap = nullptr;
-
+	SizeT mMax = 0;
+	SizeT mMax2 = 0; // !!!!
+	
     public:
 	DEFAULT_MEMBERS(FPos);
 
-	inline FPos(SizeT ext, const SizeT* map);
+	inline FPos(SizeT ext, const SizeT* map, SizeT max = 0, SizeT max2 = 0);
 
 	constexpr SizeT size() const;
 	constexpr const SizeT& val() const;
@@ -118,8 +120,15 @@ namespace CNORXZ
 	constexpr decltype(auto) operator<<(const PosT& a) const;
 
 	explicit constexpr operator SizeT() const;
+
+	constexpr SizeT max() const;
+	constexpr SizeT max2() const { return mMax2; }
+	constexpr const SizeT* map() const { return mMap; }
     };
 
+    template <class PosT>
+    constexpr decltype(auto) mkFPos(const PosT& pos, const SizeT* map);
+    
     template <SizeT N, SizeT... Ms>
     class SFPos
     {
@@ -151,6 +160,8 @@ namespace CNORXZ
 	explicit constexpr operator FPos() const;
 
 	explicit constexpr operator SizeT() const;
+
+	constexpr SizeT max() const;
     };
 
     template <class BPosT, class NPosT>
@@ -189,6 +200,9 @@ namespace CNORXZ
 	constexpr decltype(auto) operator<<(const PosT& a) const;
     };
 
+    template <class BPosT, class NPosT>
+    constexpr decltype(auto) mkMPos(const BPosT& bpos, const NPosT& npos);
+
     // treat as scalar pos!
     class DPos : public ObjHandle<VPosBase>
     {

src/include/xpr/xpr.cc.h

@@ -15,3 +15,4 @@
 #include "for.cc.h"
 #include "index_id.cc.h"
 #include "func.cc.h"
+#include "map_xpr.cc.h"

src/include/xpr/xpr.h

@@ -15,5 +15,6 @@
 #include "for.h"
 #include "index_id.h"
 #include "func.h"
+#include "map_xpr.h"
 
 #include "xpr.cc.h"

src/lib/base/config.cc

@@ -31,5 +31,17 @@ namespace CNORXZ
 	{
 	    return String(CXX_FLAGS);
 	}
+
+	String definitions()
+	{
+	    String o;
+#ifdef HAVE_CEREAL
+	    o += "HAVE_CEREAL ";
+#endif
+#ifdef CXZ_HAVE_AVX
+	    o += "CXZ_HAVE_AVX ";
+#endif
+	    return o;
+	}
     }
 }

src/lib/ranges/crange.cc

@@ -115,7 +115,7 @@ namespace CNORXZ
     
     CIndex& CIndex::at(const SizeT& metaPos)
     {
-	IB::mPos = metaPos;
+	IB::mPos = metaPos < lmax().val() ? metaPos : lmax().val();
 	return *this;
     }
 
@@ -204,6 +204,9 @@ namespace CNORXZ
     
     SizeT CRange::getMeta(SizeT metaPos) const
     {
+	if(metaPos >= size()){
+	    return size();
+	}
 	return metaPos;
     }
     

src/lib/ranges/yrange.cc

@@ -64,10 +64,16 @@ namespace CNORXZ
     {
 	mLex = 0;
 	IB::mPos = 0;
+	bool outOfScope = false;
 	for(SizeT i = 0; i != dim(); ++i){
+	    outOfScope |= mIs[i]->lex() >= mIs[i]->lmax().val();
 	    mLex += mIs[i]->lex() * mLexFormat[i].val();
 	    IB::mPos += mIs[i]->pos() * mFormat[i].val();
 	}
+	if(outOfScope){
+	    IB::mPos = pmax().val();
+	    mLex = lmax().val();
+	}
     }
     
     inline void YIndex::up(SizeT i)
@@ -668,6 +674,11 @@ namespace CNORXZ
 	return YRangeFactory(rs).create();
     }
 
+    Sptr<YRange> xplYrange(const Vector<RangePtr>& rs)
+    {
+	return std::dynamic_pointer_cast<YRange>( yrange( rs ) );
+    }
+
     /*=================+
      |   Range Casts   |
      +=================*/

src/opt/hdf5/COMMENT

@@ -1 +0,0 @@
-coming soon

src/opt/hdf5/lib/h5_group.cc

@@ -93,6 +93,7 @@ namespace CNORXZ
 		const String next = name.substr(delimpos+1);
 		auto g = getGroup(thisname);
 		g->open();
+		CHECK;
 		return g->get(next);
 	    }
 	    auto i = this->getIndexTo(thisname);
@@ -101,7 +102,7 @@ namespace CNORXZ
 
 	Sptr<Group> Group::getGroup(const String& name) const
 	{
-	    auto group = this->get(name);
+	    ContentPtr group = this->get(name);
 	    CXZ_ASSERT(group->type() == ContentType::GROUP,
 		       "element '" << name << "' is not of type GROUP");
 	    return std::dynamic_pointer_cast<Group>( group );
@@ -236,6 +237,7 @@ namespace CNORXZ
 	    auto dvec = [](const String& n) { return Vector<DType>({DType(n)}); };
 	    auto i = mCont.begin();
 	    i.at(dvec(name));
+	    CXZ_ASSERT(i != mCont.end(), "no element '" << name << "' in group " << path());
 	    return i;
 	}
 
@@ -245,6 +247,7 @@ namespace CNORXZ
 	    auto dvec = [](const String& n) { return Vector<DType>({DType(n)}); };
 	    auto i = mCont.begin();
 	    i.at(dvec(name));
+	    CXZ_ASSERT(i != mCont.end(), "no element '" << name << "' in group " << path());
 	    return i;
 	}
 	

src/opt/mpi/CMakeLists.txt

@@ -0,0 +1,16 @@
+
+find_package(MPI REQUIRED)
+if(MPI_FOUND)
+  include_directories(${MPI_C_INCLUDE_DIRS})
+else()
+  message(FATAL_ERROR "HDF5 not found")
+endif()
+message(STATUS "mpi lib = ${MPI_C_LIBRARIES}")
+set(MPI_LIBS ${MPI_LIBRARIES})
+
+include_directories(${CMAKE_CURRENT_SOURCE_DIR}/include)
+
+add_subdirectory(lib)
+add_subdirectory(tests)
+
+install(DIRECTORY include/ DESTINATION ${INSTALL_PATH}/include/cnorxz/mpi)

src/opt/mpi/include/cnorxz_mpi.cc.h

@@ -0,0 +1,17 @@
+// -*- C++ -*-
+/**
+
+   @file opt/mpi/include/cnorxz_mpi.cc.h
+   @brief CNORXZ MPI template sources header
+
+   Copyright (c) 2024 Christian Zimmermann. All rights reserved.
+   Mail: chizeta@f3l.de
+
+ **/
+
+#include "mpi_wrappers.cc.h"
+#include "rop_types.cc.h"
+#include "rrange.cc.h"
+#include "raindex.cc.h"
+#include "rarray.cc.h"
+#include "rmap_xpr.cc.h"

src/opt/mpi/include/cnorxz_mpi.h

@@ -0,0 +1,21 @@
+// -*- C++ -*-
+/**
+
+   @file opt/mpi/include/cnorxz_mpi.h
+   @brief CNORXZ MPI main header
+
+   Copyright (c) 2024 Christian Zimmermann. All rights reserved.
+   Mail: chizeta@f3l.de
+
+ **/
+
+#include "mpi_base.h"
+#include "mpi_wrappers.h"
+#include "rrange.h"
+#include "raindex.h"
+#include "rarray.h"
+#include "rop_types.h"
+#include "typemap.h"
+#include "rmap_xpr.h"
+
+#include "cnorxz_mpi.cc.h"

src/opt/mpi/include/mpi_base.h

@@ -0,0 +1,48 @@
+// -*- C++ -*-
+
+/**
+   
+   @file
+   @brief
+
+   Copyright (c) 2024 Christian Zimmermann. All rights reserved.
+   Mail: chizeta@f3l.de
+   
+**/
+
+#ifndef __cxz_mpi_base_h__
+#define __cxz_mpi_base_h__
+
+#include "cnorxz.h"
+#include "mpi.h"
+
+namespace CNORXZ
+{
+    namespace mpi
+    {
+	// class declarations
+
+	template <class RangeI, class RangeK>
+	class RRange;
+
+	template <class IndexI, class IndexK>
+	class RIndex;
+
+	template <typename T>
+	class RAIndex;
+
+	template <typename T>
+	class RArray;
+	
+	// wrapper functions
+	
+	/** Get number of THIS rank. */
+	SizeT getRankNumber();
+
+	/** Get total number of ranks. */
+	SizeT getNumRanks();
+	
+    } // namespace mpi
+} // namespace CNORXZ
+
+#endif

src/opt/mpi/include/mpi_wrappers.cc.h

@@ -0,0 +1,107 @@
+// -*- C++ -*-
+/**
+
+   @file opt/mpi/include/mpi_wrappers.cc.h
+   @brief MPI wrapper template implementations.
+
+   Copyright (c) 2024 Christian Zimmermann. All rights reserved.
+   Mail: chizeta@f3l.de
+   
+**/
+
+#ifndef __cxz_mpi_wrappers_cc_h__
+#define __cxz_mpi_wrappers_cc_h__
+
+#include "mpi_wrappers.h"
+
+namespace CNORXZ
+{
+    namespace mpi
+    {
+	
+	template <typename T>
+	void BCast<T>::bcast(T& d, SizeT root)
+	{
+	    static_assert( Typemap<T>::exists, "no bcast implementation for given type" );
+	    const int ret = MPI_Bcast( reinterpret_cast<void*>(&d), 1,
+				       Typemap<T>::value(), root, MPI_COMM_WORLD );
+	    CXZ_ASSERT(ret == MPI_SUCCESS, "got bcast error = " << ret);
+	    return;
+	}
+
+	// preliminary solutions, use cereal serialization in the future!!!
+	template <typename T>
+	void BCast<Vector<T>>::bcast(Vector<T>& d, SizeT root)
+	{
+	    SizeT size = d.size();
+	    const int ret = MPI_Bcast( reinterpret_cast<void*>(&size), 1,
+				       MPI_UNSIGNED_LONG, root, MPI_COMM_WORLD );
+	    CXZ_ASSERT(ret == MPI_SUCCESS, "got bcast error = " << ret);
+	    if(size != d.size()){
+		d.resize(size);
+	    }
+	    if constexpr( BCast<T>::special ){
+		for(auto& x: d){
+		    bcast(x, root);
+		}
+	    }
+	    else {
+		const int ret2 = MPI_Bcast( reinterpret_cast<void*>(d.data()), size,
+					    Typemap<T>::value(), root, MPI_COMM_WORLD );
+		CXZ_ASSERT(ret2 == MPI_SUCCESS, "got bcast error = " << ret2);
+	    }
+	}
+
+	template <typename T, SizeT N>
+	void BCast<Arr<T,N>>::bcast(Arr<T,N>& d, SizeT root)
+	{
+	    if constexpr( BCast<T>::special ){
+		for(auto& x: d){
+		    bcast(x, root);
+		}
+	    }
+	    else {
+		const int ret = MPI_Bcast( reinterpret_cast<void*>(d.data()), N,
+					   Typemap<T>::value(), root, MPI_COMM_WORLD );
+		CXZ_ASSERT(ret == MPI_SUCCESS, "got bcast error = " << ret);
+	    }
+	}
+
+	inline void BCast<String>::bcast(String& d, SizeT root)
+	{
+	    SizeT size = d.size();
+	    const int ret = MPI_Bcast( reinterpret_cast<void*>(&size), 1,
+				       MPI_UNSIGNED_LONG, root, MPI_COMM_WORLD );
+	    CXZ_ASSERT(ret == MPI_SUCCESS, "got bcast error = " << ret);
+	    if(size != d.size()){
+		d.resize(size);
+	    }
+	    const int ret2 = MPI_Bcast( reinterpret_cast<void*>(d.data()), size,
+				       MPI_CHAR, root, MPI_COMM_WORLD );
+	    CXZ_ASSERT(ret2 == MPI_SUCCESS, "got bcast error = " << ret2);
+	}
+
+	template <typename... Ts>
+	void BCast<Tuple<Ts...>>::bcast(Tuple<Ts...>& d, SizeT root)
+	{
+	    if constexpr( ( BCast<Ts>::special or ... ) ){
+		iter<0,sizeof...(Ts)>( [&](auto i) { bcast( std::get<i>(d), root ); }, NoF {} );
+	    }
+	    else {
+		const int ret = MPI_Bcast( reinterpret_cast<void*>(&d), sizeof(d),
+					   MPI_BYTE, root, MPI_COMM_WORLD );
+		CXZ_ASSERT(ret == MPI_SUCCESS, "got bcast error = " << ret);
+	    }
+	}
+
+	template <typename T>
+	void bcast(T& d, SizeT root)
+	{
+	    BCast<T>::bcast(d, root);
+	    return;
+	}
+	
+    } // namespace mpi
+} // namespace CNORXZ
+
+#endif

src/opt/mpi/include/mpi_wrappers.h

@@ -0,0 +1,65 @@
+// -*- C++ -*-
+/**
+
+   @file opt/mpi/include/mpi_wrappers.h
+   @brief MPI wrapper template delarations.
+
+   Copyright (c) 2024 Christian Zimmermann. All rights reserved.
+   Mail: chizeta@f3l.de
+   
+**/
+
+#ifndef __cxz_mpi_wrappers_h__
+#define __cxz_mpi_wrappers_h__
+
+#include "mpi_base.h"
+#include "typemap.h"
+
+namespace CNORXZ
+{
+    namespace mpi
+    {
+	
+	template <typename T>
+	struct BCast
+	{
+	    static constexpr bool special = false;
+	    static void bcast(T& d, SizeT root);
+	};
+
+	template <typename T>
+	struct BCast<Vector<T>>
+	{
+	    static constexpr bool special = true;
+	    static void bcast(Vector<T>& d, SizeT root);
+	};
+
+	template <typename T, SizeT N>
+	struct BCast<Arr<T,N>>
+	{
+	    static constexpr bool special = true;
+	    static void bcast(Arr<T,N>& d, SizeT root);
+	};
+
+	template <>
+	struct BCast<String>
+	{
+	    static constexpr bool special = true;
+	    static inline void bcast(String& d, SizeT root);
+	};
+
+	template <typename... Ts>
+	struct BCast<Tuple<Ts...>>
+	{
+	    static constexpr bool special = true;
+	    static void bcast(Tuple<Ts...>& d, SizeT root);
+	};
+
+	template <typename T>
+	void bcast(T& d, SizeT root);
+
+
+    } // namespace mpi
+} // namespace CNORXZ
+
+#endif

src/opt/mpi/include/raindex.cc.h

@@ -0,0 +1,174 @@
+// -*- C++ -*-
+/**
+   
+   @file opt/mpi/include/raindex.cc.h
+   @brief RAIndex template implementations.
+
+   Copyright (c) 2024 Christian Zimmermann. All rights reserved.
+   Mail: chizeta@f3l.de
+   
+**/
+
+#ifndef __cxz_mpi_raindex_cc_h__
+#define __cxz_mpi_raindex_cc_h__
+
+#include <cstring>
+#include "raindex.h"
+
+namespace CNORXZ
+{
+    namespace mpi
+    {
+	/*===============+
+	 |    RAIndex    |
+	 +===============*/
+	 
+	template <typename T>
+	RAIndex<T>::RAIndex(const T* loc, const RangePtr& range, SizeT lexpos) :
+	    RIndex<YIndex,YIndex>(range, lexpos),
+	    mLoc(loc),
+	    mMin(0),
+	    mMax(0),
+	    mThisRank(getRankNumber())
+	{
+	    setBufferSize();
+	    setBuffer();
+	}
+	    
+
+	template <typename T>
+	RAIndex<T>::RAIndex(const T* loc, const RIndex<YIndex,YIndex>& i) :
+	    RIndex<YIndex,YIndex>(i),
+	    mLoc(loc),
+	    mMin(0),
+	    mMax(0),
+	    mThisRank(getRankNumber())
+	{
+	    setBufferSize();
+	    setBuffer();
+	}
+	    
+
+	template <typename T>
+	RAIndex<T>::RAIndex(const T* loc, const RAIndex<T>& i) :
+	    RAIndex(i)
+	{
+	    mLoc = loc;
+	    setBuffer();
+	}
+
+	template <typename T>
+	RAIndex<T> RAIndex<T>::operator+(Int n) const
+	{
+	    RAIndex<T> o = *this;
+	    return o += n;
+	}
+
+	template <typename T>
+	RAIndex<T> RAIndex<T>::operator-(Int n) const
+	{
+	    RAIndex<T> o = *this;
+	    return o -= n;
+	}
+
+	template <typename T>
+	const T& RAIndex<T>::operator*() const
+	{
+	    if(this->pos() < mMin or this->pos() >= mMax){
+		setBuffer();
+	    }
+	    if(rank() != mThisRank){
+		return mBuf[local()->pos() % mBufSize];
+	    }
+	    else {
+		return mLoc[local()->pos()];
+	    }
+	}
+
+	template <typename T>
+	const T* RAIndex<T>::operator->() const
+	{
+	    if(this->pos() < mMin or this->pos() >= mMax){
+		setBuffer();
+	    }
+	    if(rank() != mThisRank){
+		return mBuf.data() + local()->pos() % mBufSize;
+	    }
+	    else {
+		return mLoc + local()->pos();
+	    }
+	}
+	
+	template <typename T>
+	void RAIndex<T>::setBufferSize()
+	{
+	    // use the contiguous stride's size as buffer size:
+	    mBufSize = 1;
+	    const auto& df = deepFormat();
+	    for(SizeT i = 0; i != df.size(); ++i){
+		if(df[i] == 1){
+		    mBufSize = deepMax()[i];
+		}
+	    }
+	}
+
+	template <typename T>
+	void RAIndex<T>::setBuffer() const
+	{
+	    if(mBuf.size() != mBufSize){
+		mBuf.resize(mBufSize);
+	    }
+	    // A Bcast alternative with const pointer to source would be better...
+	    const T* d = mLoc + (local()->pos() / mBufSize) * mBufSize;
+	    std::memcpy(mBuf.data(), d, mBufSize*sizeof(T));
+	    MPI_Bcast(mBuf.data(), mBufSize*sizeof(T), MPI_BYTE, static_cast<int>(rank()),
+		      MPI_COMM_WORLD );
+	    mMin = (this->pos() / mBufSize) * mBufSize;
+	    mMax = (this->pos() / mBufSize + 1) * mBufSize;
+	}
+
+	/*===============+
+	 |    RBIndex    |
+	 +===============*/
+
+	template <typename T>
+	RBIndex<T>::RBIndex(T* loc, const RangePtr& range, SizeT lexpos) :
+	    RAIndex<T>(loc, range, lexpos),
+	    mLoc(loc)
+	{}
+
+	template <typename T>
+	RBIndex<T>::RBIndex(T* loc, const RAIndex<T>& i) :
+	    RAIndex<T>(loc, i),
+	    mLoc(loc)
+	{}
+
+	template <typename T>
+	RBIndex<T> RBIndex<T>::operator+(Int n) const
+	{
+	    RBIndex<T> o = *this;
+	    return o += n;
+	}
+
+	template <typename T>
+	RBIndex<T> RBIndex<T>::operator-(Int n) const
+	{
+	    RBIndex<T> o = *this;
+	    return o -= n;
+	}
+
+	template <typename T>
+	RBIndex<T>& RBIndex<T>::set(const T& val)
+	{
+	    if(AI::rank() == AI::mThisRank){
+		mLoc[AI::local()->pos()] = val;
+	    }
+	    AI::setBuffer();
+	    return *this;
+	}
+
+
+    } // namespace mpi
+} // namespace CNOXRZ
+
+#endif

src/opt/mpi/include/raindex.h

@@ -0,0 +1,86 @@
+// -*- C++ -*-
+/**
+   
+   @file opt/mpi/include/raindex.h
+   @brief RAIndex declarations.
+
+   Copyright (c) 2024 Christian Zimmermann. All rights reserved.
+   Mail: chizeta@f3l.de
+   
+**/
+
+#ifndef __cxz_mpi_raindex_h__
+#define __cxz_mpi_raindex_h__
+
+#include "cnorxz.h"
+#include "rrange.h"
+
+namespace CNORXZ
+{
+    namespace mpi
+    {
+
+	/** ****
+	    Index for multi-ranked array.
+	    This is meant to be a global index, its position must be equal on each rank,
+	    otherwise deadlocks will occur. This kind of index exists for
+	    completeness (global array element access), it should not be used in
+	    operations that require high performance.
+	    For local iterations (which is what one usually does) a YIndex is sufficient.
+	    @tparam T data type.
+	*/
+	template <typename T>
+	class RAIndex : public RIndex<YIndex,YIndex>
+	{
+	public:
+	    typedef typename RIndex<YIndex,YIndex>::IB IB;
+
+	    DEFAULT_MEMBERS(RAIndex);
+	    RAIndex(const T* loc, const RangePtr& range, SizeT lexpos = 0);
+	    RAIndex(const T* loc, const RIndex<YIndex,YIndex>& i);
+	    RAIndex(const T* loc, const RAIndex<T>& i);
+
+	    RAIndex operator+(Int n) const;
+	    RAIndex operator-(Int n) const;
+
+	    const T& operator*() const;
+	    const T* operator->() const;
+
+	    void setBufferSize();
+	    
+	private:
+
+	    void setBuffer() const;
+	    
+	    const T* mLoc = nullptr;
+	    mutable Vector<T> mBuf; // used if iterating over content on different rank
+	    mutable SizeT mMin; // current global minimum position in the buffer
+	    mutable SizeT mMax; // current global maximum position in the buffer
+	    SizeT mBufSize;
+	    SizeT mThisRank;
+	};
+
+	template <typename T>
+	class RBIndex : public RAIndex<T>
+	{
+	public:
+	    typedef RAIndex<T> AI;
+	    typedef typename AI::IB IB;
+
+	    DEFAULT_MEMBERS(RBIndex);
+	    RBIndex(T* loc, const RangePtr& range, SizeT lexpos = 0);
+	    RBIndex(T* loc, const RAIndex<T>& i);
+
+	    RBIndex operator+(Int n) const;
+	    RBIndex operator-(Int n) const;
+
+	    RBIndex& set(const T& val);
+	    
+	private:
+	    T* mLoc = nullptr;
+	};
+	
+    } // namespace mpi
+} // namespace CNOXRZ
+
+#endif

src/opt/mpi/include/rarray.cc.h

@@ -0,0 +1,524 @@
+// -*- C++ -*-
+/**
+   
+   @file opt/mpi/include/rarray.cc.h
+   @brief RArray template implementations.
+
+   Copyright (c) 2024 Christian Zimmermann. All rights reserved.
+   Mail: chizeta@f3l.de
+   
+**/
+
+#ifndef __cxz_mpi_rarray_cc_h__
+#define __cxz_mpi_rarray_cc_h__
+
+#include "rarray.h"
+#include "raindex.h"
+#include "rop_types.h"
+
+namespace CNORXZ
+{
+    namespace mpi
+    {
+
+	/*===============+
+	 |    RCArray    |
+	 +===============*/
+	
+	template <typename T>
+	RCArray<T>::RCArray(const RCArray& a) :
+	    mA(a.mA->copy()),
+	    mGeom(a.mGeom),
+	    mGlobal(a.mGlobal)
+	{}
+
+	template <typename T>
+	RCArray<T>& RCArray<T>::operator=(const RCArray& a)
+	{
+	    mA = ObjHandle<CArrayBase<T>>(a.mA->copy());
+	    mGeom = a.mGeom;
+	    mGlobal = a.mGlobal;
+	}
+
+	template <typename T>
+	RCArray<T>::RCArray(const CArrayBase<T>& a, const RangePtr& geom) :
+	    mA(a.copy()),
+	    mGeom(geom),
+	    mGlobal(RRangeFactory(rangeCast<YRange>(a.range()),rangeCast<YRange>(mGeom)).create())
+	{}
+
+	template <typename T>
+	template <class IndexI, class IndexK>
+	RCArray<T>::RCArray(const Sptr<RRange<IndexI,IndexK>>& range) :
+	    mA(std::make_unique<T>(range->local())),
+	    mGeom(range->geom()),
+	    mGlobal(range)
+	{}
+
+	template <typename T>
+	template <typename I, typename M>
+	T RCArray<T>::operator[](const IndexInterface<I,M>& i) const
+	{
+	    if constexpr(is_rank_index<I>::value){
+		return *(begin() + i.lex());
+	    }
+	    else {
+		return (*mA)[i];
+	    }
+	}
+
+	template <typename T>
+	template <typename I, typename M>
+	T RCArray<T>::at(const IndexInterface<I,M>& i) const
+	{
+	    if constexpr(is_rank_index<I>::value){
+		CXZ_ASSERT(i.rank() == getNumRanks(), "rank number out of scope");
+		// consistency checks!!!
+		return *(begin() + i.lex());
+	    }
+	    else {
+		return mA->at(i);
+	    }
+	}
+
+	template <typename T>
+	template <class... Indices>
+	T RCArray<T>::operator[](const SPack<Indices...>& pack) const
+	{
+	    typedef typename std::remove_reference<decltype(*pack[CSizeT<0>{}])>::type I0;
+	    if constexpr(is_rank_index<I0>::value){
+		return *(begin() + pack.lex());
+	    }
+	    else {
+		return (*mA)[pack];
+	    }
+	}
+
+	template <typename T>
+	template <class... Indices>
+	T RCArray<T>::at(const SPack<Indices...>& pack) const
+	{
+	    typedef typename std::remove_reference<decltype(*pack[CSizeT<0>{}])>::type I0;
+	    if constexpr(is_rank_index<I0>::value){
+		CXZ_ASSERT(pack[CSizeT<0>{}]->rank() == getNumRanks(),
+			   "rank number out of scope");
+		return *(begin() + pack.lex());
+	    }
+	    else {
+		return mA->at(pack);
+	    }
+	}
+
+	template <typename T>
+	T RCArray<T>::operator[](const DPack& pack) const
+	{
+	    // TODO: assert that none of the indices is rank index
+	    return (*mA)[pack];
+	}
+
+	template <typename T>
+	T RCArray<T>::at(const DPack& pack) const
+	{
+	    // TODO: assert that none of the indices is rank index
+	    return mA->at(pack);
+	}
+
+	template <typename T>
+	template <typename I, typename M>
+	Sptr<CArrayBase<T>> RCArray<T>::sl(const IndexInterface<I,M>& begin,
+					   const IndexInterface<I,M>& end) const
+	{
+	    if constexpr(is_rank_index<I>::value){
+		CXZ_ERROR("not implemented");
+		return nullptr;
+	    }
+	    else {
+		return mA->sl(begin, end);
+	    }
+	}
+
+	template <typename T>
+	template <class Index>
+	COpRoot<T,Index> RCArray<T>::operator()(const Sptr<Index>& i) const
+	{
+	    CXZ_ERROR("not implemented");
+	    return COpRoot<T,Index>();
+	}
+
+	template <typename T>
+	template <class... Indices>
+	inline decltype(auto) RCArray<T>::operator()(const SPack<Indices...>& pack) const
+	{
+	    typedef typename std::remove_reference<decltype(*pack[CSizeT<0>{}])>::type I0;
+	    if constexpr(is_rank_index<I0>::value){
+		// preliminary:
+		CXZ_ASSERT(mA->formatIsTrivial(),
+			   "array has non-trivial format, rank operations require trivial format");
+		auto ri = pack[CSizeT<0>{}];
+		auto li = iter<1,sizeof...(Indices)>
+		    ( [&](auto i) { return pack[CSizeT<i>{}]; },
+		      [](const auto&... x) { return mindexPtr( (x * ...) ); } );
+		return croproot(*this, ri, li);
+	    }
+	    else {
+		return (*mA)(pack);
+	    }
+	}
+
+	template <typename T>
+	inline decltype(auto) RCArray<T>::operator()(const DPack& pack) const
+	{
+	    // TODO: assert that none of the indices is rank index
+	    return (*mA)(pack);
+	}
+
+	template <typename T>
+	const T* RCArray<T>::data() const
+	{
+	    return mA->data();
+	}
+
+	template <typename T>
+	SizeT RCArray<T>::size() const
+	{
+	    return mA->size() * mGeom->size();
+	}
+
+	template <typename T>
+	RangePtr RCArray<T>::range() const
+	{
+	    return mGlobal;
+	}
+
+	template <typename T>
+	typename RCArray<T>::const_iterator RCArray<T>::begin() const
+	{
+	    return const_iterator(mA->data(), mGlobal);
+	}
+
+	template <typename T>
+	typename RCArray<T>::const_iterator RCArray<T>::end() const
+	{
+	    return const_iterator(mA->data(), mGlobal, mGlobal->size());
+	}
+
+	template <typename T>
+	typename RCArray<T>::const_iterator RCArray<T>::cbegin() const
+	{
+	    return const_iterator(mA->data(), mGlobal);
+	}
+
+	template <typename T>
+	typename RCArray<T>::const_iterator RCArray<T>::cend() const
+	{
+	    return const_iterator(mA->data(), mGlobal, mGlobal->size());
+	}
+
+	template <typename T>
+	bool RCArray<T>::isView() const
+	{
+	    return mA->isView();
+	}
+
+	template <typename T>
+	const CArrayBase<T>& RCArray<T>::local() const
+	{
+	    return *mA;
+	}
+
+	template <typename T>
+	RangePtr RCArray<T>::geom() const
+	{
+	    return mGeom;
+	}
+
+	template <typename T>
+	const Vector<const T*>& RCArray<T>::buffermap() const
+	{
+	    return mMap;
+	}
+
+	template <typename T>
+	template <class SrcI, class RSrcI, class I>
+	void RCArray<T>::load(const Sptr<RIndex<SrcI,RSrcI>>& ai,
+			      const Sptr<I>& i, const Vector<bool>& required) const
+	{
+	    mA->checkFormatCompatibility(mindex(ai->local()*i));
+	    const SizeT blocks = i->pmax().val();
+	    setupBuffer(ai, required, *mA, mBuf, mMap, blocks);
+	}
+	
+	/*==============+
+	 |    RArray    |
+	 +==============*/
+
+	template <typename T>
+	RArray<T>::RArray(const RArray& a) :
+	    RCArray<T>(a),
+	    mB(dynamic_cast<ArrayBase<T>*>(RCA::mA.get()))
+	{}
+
+	template <typename T>
+	RArray<T>& RArray<T>::operator=(const RArray& a)
+	{
+	    RCArray<T>::operator=(a);
+	    mB = dynamic_cast<ArrayBase<T>*>(RCA::mA.get());
+	    return *this;
+	}
+
+	template <typename T>
+	template <class IndexI, class IndexK>
+	RArray<T>::RArray(const Sptr<RRange<IndexI,IndexK>>& range) :
+	    RCArray<T>(range),
+	    mB(dynamic_cast<ArrayBase<T>*>(RCA::mA.get()))
+	{}
+	    
+	template <typename T>
+	RArray<T>::RArray(const ArrayBase<T>& a, const RangePtr& geom) :
+	    RCArray<T>(a, geom),
+	    mB(dynamic_cast<ArrayBase<T>*>(RCA::mA.get()))
+	{}
+
+	template <typename T>
+	template <typename I, typename M>
+	RArray<T>& RArray<T>::set(const IndexInterface<I,M>& i, const T& val)
+	{
+	    auto it = begin() + i.lex();
+	    it.set(val);
+	    return *this;
+	}
+
+	template <typename T>
+	template <class... Indices>
+	RArray<T>& RArray<T>::set(const SPack<Indices...>& pack, const T& val)
+	{
+	    auto it = begin() + pack.lex();
+	    it.set(val);
+	    return *this;
+	}
+
+	template <typename T>
+	RArray<T>& RArray<T>::set(const DPack& pack, const T& val)
+	{
+	    auto it = begin() + pack.lex();
+	    it.set(val);
+	    return *this;
+	}
+
+	template <typename T>
+	template <typename I, typename M>
+	Sptr<ArrayBase<T>> RArray<T>::sl(const IndexInterface<I,M>& begin,
+					 const IndexInterface<I,M>& end)
+	{
+	    if constexpr(is_rank_index<I>::value){
+		CXZ_ERROR("not implemented");
+		return nullptr;
+	    }
+	    else {
+		return mB->sl(begin, end);
+	    }
+	}
+
+	template <typename T>
+	template <class Index>
+	COpRoot<T,Index> RArray<T>::operator()(const Sptr<Index>& i) const
+	{
+	    CXZ_ERROR("not implemented");
+	    return COpRoot<T,Index>();
+	}
+
+	template <typename T>
+	template <class... Indices>
+	inline decltype(auto) RArray<T>::operator()(const SPack<Indices...>& pack) const
+	{
+	    typedef typename std::remove_reference<decltype(*pack[CSizeT<0>{}])>::type I0;
+	    if constexpr(is_rank_index<I0>::value){
+		// preliminary:
+		CXZ_ASSERT(mB->formatIsTrivial(),
+			   "array has non-trivial format, rank operations require trivial format");
+		auto ri = pack[CSizeT<0>{}];
+		auto li = iter<1,sizeof...(Indices)>
+		    ( [&](auto i) { return pack[CSizeT<i>{}]; },
+		      [](const auto&... x) { return mindexPtr( (x * ...) ); } );
+		return croproot(*this, ri, li);
+	    }
+	    else {
+		return (*mB)(pack);
+	    }
+	}
+
+	template <typename T>
+	inline decltype(auto) RArray<T>::operator()(const DPack& pack) const
+	{
+	    // TODO: assert that none of the indices is rank index
+	    return (*mB)(pack);
+	}
+
+	template <typename T>
+	template <class Index>
+	OpRoot<T,Index> RArray<T>::rop(const Sptr<Index>& i)
+	{
+	    return (*mB)(i);
+	}
+	
+	template <typename T>
+	template <class... Indices>
+	inline decltype(auto) RArray<T>::rop(const SPack<Indices...>& pack)
+	{
+	    typedef typename std::remove_reference<decltype(*pack[CSizeT<0>{}])>::type I0;
+	    if constexpr(is_rank_index<I0>::value){
+		// preliminary:
+		CXZ_ASSERT(mB->formatIsTrivial(),
+			   "array has non-trivial format, rank operations require trivial format");
+		/*
+		auto ri = pack[CSizeT<0>{}];
+		auto li = iter<1,sizeof...(Indices)>
+		    ( [&](auto i) { return pack[CSizeT<i>{}]; },
+		      [](const auto&... x) { return mindexPtr( (x * ...) ); } );
+		*/
+		return oproot(*mB, mindexPtr(pack));
+	    }
+	    else {
+		return (*mB)(pack);
+	    }
+	}
+	
+	template <typename T>
+	inline decltype(auto) RArray<T>::rop(const DPack& pack)
+	{
+	    return (*mB)(pack);
+	}
+
+	template <typename T>
+	T* RArray<T>::data()
+	{
+	    return mB->data();
+	}
+
+	template <typename T>
+	typename RArray<T>::iterator RArray<T>::begin()
+	{
+	    return iterator(mB->data(), RCA::mGlobal);
+	}
+
+	template <typename T>
+	typename RArray<T>::iterator RArray<T>::end()
+	{
+	    return iterator(mB->data(), RCA::mGlobal, RCA::mGlobal->size());
+	}
+	    
+	template <typename T>
+	ArrayBase<T>& RArray<T>::local()
+	{
+	    return *mB;
+	}
+	
+	/*============================+
+	 |    non-member functions    |
+	 +============================*/
+
+	template <class SrcI, class RSrcI, typename T>
+	void setupBuffer(const Sptr<RIndex<SrcI,RSrcI>>& rgj,
+			 const Vector<bool>& required, const CArrayBase<T>& data,
+			 Vector<T>& buf, Vector<const T*>& map, const SizeT blocks)
+	{
+	    const SizeT myrank = getRankNumber();
+	    const SizeT Nranks = getNumRanks();
+	    const SizeT locsz = rgj->local()->lmax().val();
+	    const SizeT mapsize = rgj->range()->size();
+	    map = Vector<const T*>(mapsize,nullptr);
+	    Vector<Vector<T>> sendbuf(Nranks);
+	    for(auto& sb: sendbuf){
+		sb.reserve(data.size());
+	    }
+	    Vector<Vector<SizeT>> request(Nranks);
+
+	    // Second loop: setup send buffer
+	    auto mi = mindexPtr(rgj->rankI(), rgj->local());
+	    mi->ifor( operation
+		      ( [&](SizeT p) {
+			  const SizeT r = p / locsz;
+			  if(myrank != r and required[p]){
+			      request[r].push_back(p % locsz);
+			  }
+		      } , posop(mi) ) ,
+		NoF {} )();
+
+	    // transfer:
+	    Vector<SizeT> reqsizes(Nranks);
+	    SizeT bufsize = 0;
+	    Vector<Vector<SizeT>> ext(Nranks);
+	    for(auto& e: ext){
+		e.resize(Nranks);
+	    }
+	    for(SizeT i = 0; i != Nranks; ++i){
+		reqsizes[i] = request[i].size();
+		bufsize += reqsizes[i]*blocks;
+		ext[myrank][i] = reqsizes[i];
+	    }
+	    buf.resize(bufsize);
+	    MPI_Status stat;
+
+	    // transfer requests:
+	    for(SizeT o = 1; o != Nranks; ++o){
+		const SizeT dstr = (myrank + o) % Nranks;
+		const SizeT srcr = (myrank - o + Nranks) % Nranks;
+		SizeT sendsize = 0;
+		MPI_Sendrecv(reqsizes.data()+dstr, 1, MPI_UNSIGNED_LONG, dstr, 0,
+			     &sendsize, 1, MPI_UNSIGNED_LONG, srcr, 0, MPI_COMM_WORLD, &stat);
+		ext[srcr][myrank] = sendsize;
+		Vector<SizeT> sendpos(sendsize);
+		MPI_Sendrecv(request[dstr].data(), reqsizes[dstr], MPI_UNSIGNED_LONG, dstr, 0,
+			     sendpos.data(), sendsize, MPI_UNSIGNED_LONG, srcr, 0, MPI_COMM_WORLD, &stat);
+		sendbuf[srcr].resize(sendsize*blocks);
+		for(SizeT i = 0; i != sendsize; ++i){
+		    std::memcpy( sendbuf[srcr].data()+i*blocks, data.data()+sendpos[i]*blocks, blocks*sizeof(T) );
+		}
+	    }
+
+	    const MPI_Datatype dt = Typemap<T>::value();
+	
+	    // transfer data:
+	    for(SizeT o = 1; o != Nranks; ++o){
+		const SizeT dstr = (myrank + o) % Nranks;
+		const SizeT srcr = (myrank - o + Nranks) % Nranks;
+		SizeT off = 0;
+		for(SizeT p = 0; p != srcr; ++p){
+		    off += ext[myrank][p];
+		}
+	
+		MPI_Sendrecv(sendbuf[dstr].data(), ext[dstr][myrank]*blocks, dt, dstr, 0,
+			     buf.data()+off*blocks, ext[myrank][srcr]*blocks, dt, srcr, 0,
+			     MPI_COMM_WORLD, &stat);
+	
+	    }
+
+	    // Third loop: Assign map to target buffer positions:
+	    const SizeT myrankoff = myrank*locsz;
+	    assert(mapsize == Nranks*locsz);
+	    Vector<SizeT> cnt(Nranks);
+	    mi->ifor( operation
+		      ( [&](SizeT p) {
+			  const SizeT r = p / locsz;
+			  const SizeT l = p % locsz;
+			  const SizeT mpidx = (p - myrankoff + mapsize) % mapsize;
+			  if(myrank != r and required[p]){
+			      SizeT off = 0;
+			      for(SizeT s = 0; s != r; ++s){
+				  off += ext[myrank][s];
+			      }
+			      map[mpidx] = buf.data() + off*blocks + cnt[r]*blocks;
+			      ++cnt[r];
+			  }
+			  if(myrank == r){
+			      assert(mpidx < locsz);
+			      map[mpidx] = data.data() + l*blocks;
+			  }
+		      } , posop(mi) ), NoF {} )();
+
+	}
+	
+    } // namespace mpi
+} // namespace CNORXZ
+
+#endif

src/opt/mpi/include/rarray.h

@@ -0,0 +1,259 @@
+// -*- C++ -*-
+/**
+   
+   @file opt/mpi/include/rarray.h
+   @brief RArray declarations.
+
+   Copyright (c) 2024 Christian Zimmermann. All rights reserved.
+   Mail: chizeta@f3l.de
+   
+**/
+
+#ifndef __cxz_mpi_rarray_h__
+#define __cxz_mpi_rarray_h__
+
+#include "cnorxz.h"
+#include "raindex.h"
+
+namespace CNORXZ
+{
+    namespace mpi
+    {
+	/** ****
+	    Wrapper for cnorxz array types for handling MPI ranks.
+	    @tparam T data type.
+	*/
+	template <typename T>
+	class RCArray
+	{
+	public:
+	    typedef RAIndex<T> const_iterator;
+
+	    DEFAULT_C(RCArray);
+	    DEFAULT_MOVE(RCArray);
+	    RCArray(const RCArray& a);
+	    RCArray& operator=(const RCArray& a);
+
+	    /** Construct from a rank range.
+		@param range The range.
+	     */
+	    template <class IndexI, class IndexK>
+	    RCArray(const Sptr<RRange<IndexI,IndexK>>& range);
+	    
+	    /** Construct from local array object.
+		@param a Local array.
+		@param geom Rank geometry.
+	     */
+	    RCArray(const CArrayBase<T>& a, const RangePtr& geom);
+
+	    /** @copydoc CArrayBase::operator[] */
+	    template <typename I, typename M>
+	    T operator[](const IndexInterface<I,M>& i) const;
+
+	    /** @copydoc CArrayBase::at()
+		Return copy instead of reference.
+	     */
+	    template <typename I, typename M>
+	    T at(const IndexInterface<I,M>& i) const;
+
+	    /** @copydoc CArrayBase::operator[] */
+	    template <class... Indices>
+	    T operator[](const SPack<Indices...>& pack) const;
+
+	    /** @copydoc CArrayBase::at() */
+	    template <class... Indices>
+	    T at(const SPack<Indices...>& pack) const;
+
+	    /** @copydoc CArrayBase::operator[] */
+	    T operator[](const DPack& pack) const;
+
+	    /** @copydoc CArrayBase::at() */
+	    T at(const DPack& pack) const;
+
+	    /** @copydoc CArrayBase::sl() */
+	    template <typename I, typename M>
+	    Sptr<CArrayBase<T>> sl(const IndexInterface<I,M>& begin,
+				   const IndexInterface<I,M>& end) const;
+
+	    /** @copydoc CArrayBase::operator() */
+	    template <class Index>
+	    COpRoot<T,Index> operator()(const Sptr<Index>& i) const;
+
+	    /** @copydoc CArrayBase::operator() */
+	    template <class... Indices>
+	    inline decltype(auto) operator()(const SPack<Indices...>& pack) const;
+
+	    /** @copydoc CArrayBase::operator() */
+	    inline decltype(auto) operator()(const DPack& pack) const;
+
+	    /** @copydoc CArrayBase::data() */
+	    const T* data() const;
+
+	    /** @copydoc CArrayBase::size() */
+	    SizeT size() const;
+
+	    /** @copydoc CArrayBase:range() */
+	    RangePtr range() const;
+
+	    /** @copydoc CArrayBase::begin() */
+	    const_iterator begin() const;
+
+	    /** @copydoc CArrayBase::end() */
+	    const_iterator end() const;
+
+	    /** @copydoc CArrayBase::cbegin() */
+	    const_iterator cbegin() const;
+
+	    /** @copydoc CArrayBase::cend() */
+	    const_iterator cend() const;
+
+	    /** @copydoc CArrayBase::isView() */
+	    bool isView() const;
+
+	    /** Get local array object. */
+	    const CArrayBase<T>& local() const;
+
+	    /** Get rank geometry. */
+	    RangePtr geom() const;
+
+	    /** Get the buffer map.
+		buffermap()[r (ranked)][i] returns the same as
+		loc[r * i] if 'loc' is a purely local array with the same content.
+		The buffer needs to be initialized before according to required range.
+	     */
+	    const Vector<const T*>& buffermap() const;
+
+	    /** Load all data from other ranks that is accessed by ai in a loop over lpi.
+		imap indicates the global (!) position of ai for a given (local) position of lpi.
+		@param ai Access rank index.
+		@param i Index for non-ranked dimensions.
+		@param required Vector of booleans indicating global positions to load.
+	     */
+	    template <class SrcI, class RSrcI, class I>
+	    void load(const Sptr<RIndex<SrcI,RSrcI>>& ai,
+		      const Sptr<I>& i, const Vector<bool>& required) const;
+
+
+	protected:
+	    ObjHandle<CArrayBase<T>> mA;
+	    RangePtr mGeom;
+	    RangePtr mGlobal;
+
+	private:
+	    mutable Vector<T> mBuf; // data from other ranks
+	    mutable Vector<const T*> mMap; // point to positions in mA or mBuf
+	};
+
+	/** ****
+	    Wrapper for cnorxz array types for handling MPI ranks.
+	    @tparam T data type.
+	*/
+	template <typename T>
+	class RArray : public RCArray<T>
+	{
+	public:
+	    typedef RCArray<T> RCA;
+	    typedef typename RCA::const_iterator const_iterator;
+	    typedef RBIndex<T> iterator;
+	
+	    using RCA::operator[];
+	    using RCA::operator();
+	    using RCA::at;
+	    using RCA::data;
+	    using RCA::begin;
+	    using RCA::end;
+	    using RCA::cbegin;
+	    using RCA::cend;
+	    using RCA::sl;
+
+	    DEFAULT_C(RArray);
+	    DEFAULT_MOVE(RArray);
+	    RArray(const RArray& a);
+	    RArray& operator=(const RArray& a);
+
+	    /** Construct from a rank range.
+		@param range The range.
+	     */
+	    template <class IndexI, class IndexK>
+	    RArray(const Sptr<RRange<IndexI,IndexK>>& range);
+	    
+	    /** Construct from local array object.
+		@param a Local array.
+		@param geom Rank geometry.
+	     */
+	    RArray(const ArrayBase<T>& a, const RangePtr& geom);
+
+	    /** Assign value at a position indicated by an index.
+		@param i Index indicating the position.
+		@param val Value.
+	     */
+	    template <typename I, typename M>
+	    RArray& set(const IndexInterface<I,M>& i, const T& val);
+
+	    /** Assign value at a position indicated by an index pack.
+		@param pack Index pack indicating the position.
+		@param val Value.
+	     */
+	    template <class... Indices>
+	    RArray& set(const SPack<Indices...>& pack, const T& val);
+
+	    /** Assign value at a position indicated by an index pack.
+		@param pack Index pack indicating the position.
+		@param val Value.
+	     */
+	    RArray& set(const DPack& pack, const T& val);
+
+	    /** @copydoc ArrayBase::sl() */
+	    template <typename I, typename M>
+	    Sptr<ArrayBase<T>> sl(const IndexInterface<I,M>& begin,
+				  const IndexInterface<I,M>& end);
+
+	    /** @copydoc ArrayBase::operator() */
+	    template <class Index>
+	    COpRoot<T,Index> operator()(const Sptr<Index>& i) const;
+
+	    /** @copydoc ArrayBase::operator() */
+	    template <class... Indices>
+	    inline decltype(auto) operator()(const SPack<Indices...>& pack) const;
+
+	    /** @copydoc ArrayBase::operator() */
+	    inline decltype(auto) operator()(const DPack& pack) const;
+
+	    /** @copydoc ArrayBase::operator() */
+	    template <class Index>
+	    OpRoot<T,Index> rop(const Sptr<Index>& i);
+
+	    /** @copydoc ArrayBase::operator() */
+	    template <class... Indices>
+	    inline decltype(auto) rop(const SPack<Indices...>& pack);
+
+	    /** @copydoc ArrayBase::operator() */
+	    inline decltype(auto) rop(const DPack& pack);
+
+	    /** @copydoc ArrayBase::data() */
+	    T* data();
+
+	    /** @copydoc ArrayBase::begin() */
+	    iterator begin();
+
+	    /** @copydoc ArrayBase::end() */
+	    iterator end();
+	    
+	    /** Get local array object. */
+	    ArrayBase<T>& local();
+
+	private:
+	    ArrayBase<T>* mB = nullptr;
+
+	};
+
+	template <class SrcI, class RSrcI, typename T>
+	void setupBuffer(const Sptr<RIndex<SrcI,RSrcI>>& rgj,
+			 const Vector<bool>& required, const CArrayBase<T>& data,
+			 Vector<T>& buf, Vector<const T*>& map, const SizeT blocks);
+
+	
+    } // namespace mpi
+} // namespace CNORXZ
+
+#endif

src/opt/mpi/include/rmap_xpr.cc.h

@@ -0,0 +1,75 @@
+// -*- C++ -*-
+/**
+
+   @file opt/include/rmap_xpr.cc.h
+   @brief MPI specific specializations for MapXpr
+
+   Copyright (c) 2024 Christian Zimmermann. All rights reserved.
+   Mail: chizeta@f3l.de
+   
+ **/
+
+#ifndef __cxz_mpi_rmap_xpr_cc_h__
+#define __cxz_mpi_rmap_xpr_cc_h__
+
+#include "rmap_xpr.h"
+#include "mpi_base.h"
+
+namespace CNORXZ
+{
+    template <class TarI, class RTarI, class SrcIndex, class F>
+    void
+    MapSetup<mpi::RIndex<TarI,RTarI>,SrcIndex,F>::setup(const Sptr<mpi::RIndex<TarI,RTarI>>& ti,
+							const Sptr<SrcIndex>& si,
+							const F& f, const Sptr<Vector<SizeT>>& m)
+    {
+	auto six = *si;
+	auto sie = si->range()->end();
+	auto tix = *ti;
+	const SizeT locsz = tix.local()->pmax().val();
+	const SizeT tarsize = locsz*mpi::getNumRanks();
+	const SizeT mapsize = m->size();
+	const SizeT myrank = mpi::getRankNumber();
+	if constexpr(mpi::is_rank_index<SrcIndex>::value){
+	    CXZ_ASSERT(mapsize == six.local()->pmax().val(), "map not well-formatted: size = "
+		       << mapsize << ", expected " << six.local()->pmax().val());
+	    for(six = 0; six != sie; ++six){
+		tix.at( f(*six) );
+		if(six.rank() == myrank){
+		    const SizeT idx = (tix.pos() - locsz*tix.rank() + tarsize) % tarsize;
+		    (*m)[six.local()->pos()] = idx;
+		}
+	    }
+	}
+	else {
+	    CXZ_ASSERT(mapsize == six.pmax().val(), "map not well-formatted: size = "
+		       << mapsize << ", expected " << six.pmax().val());
+	    for(six = 0; six != sie; ++six){
+		tix.at( f(*six) );
+		const SizeT idx = (tix.pos() - locsz*tix.rank() + tarsize) % tarsize;
+		(*m)[six.pos()] = idx;
+	    }
+	}
+    }
+
+    template <class TarI, class RTarI, class SrcIndex, class F>
+    Sptr<Vector<SizeT>>
+    MapSetup<mpi::RIndex<TarI,RTarI>,SrcIndex,F>::setup(const Sptr<mpi::RIndex<TarI,RTarI>>& ti,
+							const Sptr<SrcIndex>& si,
+							const F& f)
+    {
+	SizeT mapsize = 0;
+	if constexpr(mpi::is_rank_index<SrcIndex>::value){
+	    mapsize = si->local()->lmax().val();
+	}
+	else {
+	    mapsize = si->lmax().val();
+	}
+	auto o = std::make_shared<Vector<SizeT>>(mapsize);
+	setup(ti,si,f,o);
+	return o;
+    }
+    
+}
+
+#endif

src/opt/mpi/include/rmap_xpr.h

@@ -0,0 +1,33 @@
+// -*- C++ -*-
+/**
+
+   @file opt/include/rmap_xpr.h
+   @brief MPI specific specializations for MapXpr
+
+   Copyright (c) 2024 Christian Zimmermann. All rights reserved.
+   Mail: chizeta@f3l.de
+   
+ **/
+
+#ifndef __cxz_mpi_rmap_xpr_h__
+#define __cxz_mpi_rmap_xpr_h__
+
+#include "xpr/map_xpr.h"
+#include "rrange.h"
+
+namespace CNORXZ
+{
+    template <class TarI, class RTarI, class SrcIndex, class F>
+    struct MapSetup<mpi::RIndex<TarI,RTarI>,SrcIndex,F>
+    {
+	static void setup(const Sptr<mpi::RIndex<TarI,RTarI>>& ti,
+			  const Sptr<SrcIndex>& si,
+			  const F& f, const Sptr<Vector<SizeT>>& m);
+
+	static Sptr<Vector<SizeT>> setup(const Sptr<mpi::RIndex<TarI,RTarI>>& ti,
+					 const Sptr<SrcIndex>& si,
+					 const F& f);
+    };
+}
+
+#endif

src/opt/mpi/include/rop_types.cc.h

@@ -0,0 +1,146 @@
+
+#include "rop_types.h"
+
+namespace CNORXZ
+{
+    namespace mpi
+    {
+
+	/*================+
+	 |    CROpRoot    |
+	 +================*/
+    
+	template <typename T, class RIndexT, class IndexT>
+	constexpr CROpRoot<T,RIndexT,IndexT>::CROpRoot(const RCArray<T>& a, const Sptr<RIndexT>& ri,
+						       const Sptr<IndexT>& li) :
+	    mData(a.buffermap().data()),
+	    mRIndex(ri),
+	    mIndex(li)
+	{
+	    CXZ_ASSERT(a.buffermap().size() == ri->lmax().val(),
+		       "data map not properly initialized: map size = " << a.buffermap().size()
+		       << ", rank index range size = " << ri->lmax().val());
+	}
+
+	template <typename T, class RIndexT, class IndexT>
+	template <class PosT>
+	constexpr decltype(auto) CROpRoot<T,RIndexT,IndexT>::operator()(const PosT& pos) const
+	{
+	    return (mData[pos.val()])[pos.next().val()];
+	}
+
+	template <typename T, class RIndexT, class IndexT>
+	constexpr decltype(auto) CROpRoot<T,RIndexT,IndexT>::operator()() const
+	{
+	    return (mData[0])[0];
+	}
+
+	template <typename T, class RIndexT, class IndexT>
+	template <SizeT I>
+	constexpr decltype(auto) CROpRoot<T,RIndexT,IndexT>::rootSteps(const IndexId<I>& id) const
+	{
+	    return mRIndex->stepSize(id) << mIndex->stepSize(id);
+	}
+
+	template <typename T, class RIndexT, class IndexT>
+	constexpr decltype(auto) croproot(const RCArray<T>& a, const Sptr<RIndexT>& ri,
+					  const Sptr<IndexT>& li)
+	{
+	    return CROpRoot<T,RIndexT,IndexT>(a,ri,li);
+	}
+    
+	/*===============+
+	 |    ROpRoot    |
+	 +===============*/
+
+	template <typename T, class RIndexT, class IndexT>
+	constexpr ROpRoot<T,RIndexT,IndexT>::ROpRoot(RArray<T>& a, const Sptr<RIndexT>& ri,
+						     const Sptr<IndexT>& li) :
+	    mData(a.data()),
+	    mRIndex(ri),
+	    mIndex(li)
+	{
+	    CXZ_ERROR("nope");
+	}
+	
+	template <typename T, class RIndexT, class IndexT>
+	template <class Op>
+	constexpr ROpRoot<T,RIndexT,IndexT>& ROpRoot<T,RIndexT,IndexT>::operator=(const Op& in)
+	{
+	    OI::a(mIndex, [](auto& a, const auto& b) { a = b; }, in);
+	    return *this;
+	}
+	
+	template <typename T, class RIndexT, class IndexT>
+	template <class Op>
+	constexpr ROpRoot<T,RIndexT,IndexT>& ROpRoot<T,RIndexT,IndexT>::operator+=(const Op& in)
+	{
+	    OI::a(mIndex, [](auto& a, const auto& b) { a += b; }, in);
+	    return *this;
+	}
+	
+	template <typename T, class RIndexT, class IndexT>
+	constexpr ROpRoot<T,RIndexT,IndexT>& ROpRoot<T,RIndexT,IndexT>::operator=(const ROpRoot& in)
+	{
+	    OI::a(mIndex, [](auto& a, const auto& b) { a = b; }, in);
+	    return *this;
+	}
+
+	template <typename T, class RIndexT, class IndexT>
+	template <class PosT>
+	constexpr decltype(auto) ROpRoot<T,RIndexT,IndexT>::operator()(const PosT& pos) const
+	{
+	    return mData[pos.val()];
+	}
+
+	template <typename T, class RIndexT, class IndexT>
+	constexpr decltype(auto) ROpRoot<T,RIndexT,IndexT>::operator()() const
+	{
+	    return mData[0];
+	}
+
+	template <typename T, class RIndexT, class IndexT>
+	template <SizeT I>
+	constexpr decltype(auto) ROpRoot<T,RIndexT,IndexT>::rootSteps(const IndexId<I>& id) const
+	{
+	    return mIndex->stepSize(id);
+	}
+    
+	template <typename T, class RIndexT, class IndexT>
+	constexpr decltype(auto) roproot(RArray<T>& a, const Sptr<RIndexT>& ri,
+					 const Sptr<IndexT>& li)
+	{
+	    return ROpRoot<T,RIndexT,IndexT>(a,ri,li);
+	}
+
+	/*====================+
+	 |    RContraction    |
+	 +====================*/
+	/*
+	template <class CXpr>
+	constexpr RContraction<CXpr>::RContraction(CXpr&& cxpr) :
+	    mCXpr(std::forward<CXpr>(cxpr))
+	{}
+
+	template <class CXpr>
+	template <class PosT>
+	constexpr decltype(auto) RContraction<CXpr>::operator()(const PosT& pos) const
+	{
+
+	}
+
+	template <class CXpr>
+	constexpr decltype(auto) RContraction<CXpr>::operator()() const
+	{
+
+	}
+
+	template <class CXpr>
+	template <SizeT I>
+	constexpr decltype(auto) RContraction<CXpr>::rootSteps(const IndexId<I>& id) const
+	{
+
+	}
+	*/
+    } // namespace mpi    
+} // namespace CNORXZ

src/opt/mpi/include/rop_types.h

@@ -0,0 +1,120 @@
+// -*- C++ -*-
+/**
+   
+   @file opt/include/rop_types.h
+   @brief Ranked Operation types declarations.
+
+   Copyright (c) 2024 Christian Zimmermann. All rights reserved.
+   Mail: chizeta@f3l.de
+
+**/
+
+#ifndef __cxz_mpi_rop_types_h__
+#define __cxz_mpi_rop_types_h__
+
+#include "operation/operation.h"
+
+namespace CNORXZ
+{
+    namespace mpi
+    {
+	template <typename T, class RIndexT, class IndexT>
+	class CROpRoot : public COpInterface<CROpRoot<T,RIndexT,IndexT>>
+	{
+	public:
+	    typedef COpInterface<CROpRoot<T,RIndexT,IndexT>> OI;
+
+	    constexpr CROpRoot() = default;
+	    constexpr CROpRoot(const RCArray<T>& a, const Sptr<RIndexT>& ri,
+			       const Sptr<IndexT>& li);
+	    
+	    template <class PosT>
+	    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;
+
+	private:
+
+	    const T* const* mData;
+	    Sptr<RIndexT> mRIndex;
+	    Sptr<IndexT> mIndex;
+	};
+
+
+	template <typename T, class RIndexT, class IndexT>
+	constexpr decltype(auto) croproot(const RCArray<T>& a, const Sptr<RIndexT>& ri,
+					  const Sptr<IndexT>& li);
+
+	template <typename T, class RIndexT, class IndexT>
+	class ROpRoot : public OpInterface<ROpRoot<T,RIndexT,IndexT>>
+	{
+	public:
+	    typedef OpInterface<ROpRoot<T,RIndexT,IndexT>> OI;
+
+	    constexpr ROpRoot() = default;
+	    constexpr ROpRoot(RArray<T>& a, const Sptr<RIndexT>& ri,
+			      const Sptr<IndexT>& li);
+
+	    template <class Op>
+	    constexpr ROpRoot& operator=(const Op& in);
+
+	    template <class Op>
+	    constexpr ROpRoot& operator+=(const Op& in);
+
+	    constexpr ROpRoot& operator=(const ROpRoot& in);
+
+	    template <class PosT>
+	    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;
+
+	private:
+
+	    T* mData;
+	    Sptr<RIndexT> mRIndex;
+	    Sptr<IndexT> mIndex;
+	};
+
+	template <typename T, class RIndexT, class IndexT>
+	constexpr decltype(auto) roproot(RArray<T>& a, const Sptr<RIndexT>& ri,
+					 const Sptr<IndexT>& li);
+
+	/*
+	template <class CXpr>
+	class RContraction : public OpInterfacte<RContraction<CXpr>>
+	{
+	public:
+	    typedef OpInterfacte<RContraction<CXpr>> OI;
+
+	    constexpr RContraction() = default;
+	    constexpr RContraction(CXpr&& cxpr);
+
+	    template <class PosT>
+	    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;
+
+	private:
+	    CXpr mCXpr;
+	};
+	*/
+    } // namespace mpi
+
+    template <typename T, class RIndexT, class IndexT>
+    struct op_size<mpi::CROpRoot<T,RIndexT,IndexT>>
+    {
+	static constexpr SizeT value = 2;
+    };
+    
+} // namespace CNORXZ
+
+#endif

src/opt/mpi/include/rrange.cc.h

@@ -0,0 +1,559 @@
+// -*- C++ -*-
+/**
+   
+   @file opt/mpi/include/rrange.cc.h
+   @brief RRange and RIndex template implementation.
+   
+   Copyright (c) 2024 Christian Zimmermann. All rights reserved.
+   Mail: chizeta@f3l.de
+   
+**/
+
+#ifndef __cxz_mpi_rrange_cc_h__
+#define __cxz_mpi_rrange_cc_h__
+
+#include "rrange.h"
+#include "mpi_wrappers.h"
+#include "mpi_wrappers.cc.h"
+
+namespace CNORXZ
+{
+    namespace mpi
+    {
+
+	/*==============+
+	 |    RIndex    |
+	 +==============*/
+
+	template <class IndexI, class IndexK>
+	RIndex<IndexI,IndexK>::RIndex(const RIndex& in) :
+	    mRange(in.mRange),
+	    mI(std::make_shared<IndexI>(mRange->local())),
+	    mK(std::make_shared<IndexK>(mRange->geom()))
+	{
+	    *this = in.lex();
+	}
+
+	template <class IndexI, class IndexK>
+	RIndex<IndexI,IndexK>& RIndex<IndexI,IndexK>::operator=(const RIndex& in)
+	{
+	    mRange = in.mRange;
+	    mI = std::make_shared<IndexI>(mRange->local());
+	    mK = std::make_shared<IndexK>(mRange->geom());
+	    *this = in.lex();
+	    return *this;
+	}
+
+	template <class IndexI, class IndexK>
+	RIndex<IndexI,IndexK>::RIndex(const RangePtr& global, SizeT lexpos) :
+	    mRange(rangeCast<RangeType>(global)),
+	    mI(std::make_shared<IndexI>(mRange->local())),
+	    mK(std::make_shared<IndexK>(mRange->geom()))
+	{
+	    *this = lexpos;
+	}
+
+	template <class IndexI, class IndexK>
+	RIndex<IndexI,IndexK>::RIndex(const Sptr<IndexI>& i, const Sptr<IndexK>& k) :
+	    mRange(rangeCast<RangeType>( RRangeFactory(i->range(), k->range()).create() )),
+	    mI(i),
+	    mK(k)
+	{
+	    (*this)();
+	}
+
+	template <class IndexI, class IndexK>
+	RIndex<IndexI,IndexK>& RIndex<IndexI,IndexK>::operator=(SizeT pos)
+	{
+	    mLex = pos; // = lex
+	    if(mLex >= lmax().val()){
+		mLex = lmax().val();
+		return *this;
+	    }
+	    // pos is the lexicographic position of the global range.
+	    // Hence, have to consider the rank geometry.
+	    const auto& i = mI->pack();
+	    const auto& k = mK->pack();
+	    const auto& ilf = mI->lexFormat();
+	    const auto& klf = mK->lexFormat();
+	    SizeT r = 0;
+	    SizeT l = 0;
+	    if constexpr(has_static_sub<IndexI>::value){
+		constexpr SizeT NI = index_dim<IndexI>::value;
+		iter<0,NI>( [&](auto mu) {
+		    const SizeT jmu = (mLex / (ilf[mu].val()*klf[mu].val())) %
+			(i[mu]->lmax().val() * k[mu]->lmax().val());
+		    r += ( jmu / i[mu]->lmax().val() ) * klf[mu].val();
+		    l += ( jmu % i[mu]->lmax().val() ) * ilf[mu].val();
+		}, NoF{} );
+	    }
+	    else if constexpr( has_static_sub<IndexK>::value){
+		constexpr SizeT NI = index_dim<IndexK>::value;
+		iter<0,NI>( [&](auto mu) {
+		    const SizeT jmu = (mLex / (ilf[mu].val()*klf[mu].val())) %
+			(i[mu]->lmax().val() * k[mu]->lmax().val());
+		    r += ( jmu / i[mu]->lmax().val() ) * klf[mu].val();
+		    l += ( jmu % i[mu]->lmax().val() ) * ilf[mu].val();
+		}, NoF{} );
+	    }
+	    else {
+		const SizeT NI = mI->dim();
+		for(SizeT mu = 0; mu != NI; ++mu){
+		    const SizeT jmu = (mLex / (ilf[mu].val()*klf[mu].val())) %
+			(i[mu]->lmax().val() * k[mu]->lmax().val());
+		    r += ( jmu / i[mu]->lmax().val() ) * klf[mu].val();
+		    l += ( jmu % i[mu]->lmax().val() ) * ilf[mu].val();
+		}
+	    }
+	    *mI = l;
+	    *mK = r;
+	    IB::mPos = mK->pos() * mI->pmax().val() + mI->pos();
+	    return *this;
+	}
+
+	template <class IndexI, class IndexK>
+	RIndex<IndexI,IndexK>& RIndex<IndexI,IndexK>::operator++()
+	{
+	    *this = lex() + 1; // room for optimization
+	    return *this;
+	}
+
+	template <class IndexI, class IndexK>
+	RIndex<IndexI,IndexK>& RIndex<IndexI,IndexK>::operator--()
+	{
+	    *this = lex() - 1; // room for optimization
+	    return *this;
+	}
+
+	template <class IndexI, class IndexK>
+	RIndex<IndexI,IndexK> RIndex<IndexI,IndexK>::operator+(Int n) const
+	{
+	    RIndex<IndexI,IndexK> o(*this);
+	    return o += n;
+	}
+
+	template <class IndexI, class IndexK>
+	RIndex<IndexI,IndexK> RIndex<IndexI,IndexK>::operator-(Int n) const
+	{
+	    RIndex<IndexI,IndexK> o(*this);
+	    return o -= n;
+	}
+
+	template <class IndexI, class IndexK>
+	SizeT RIndex<IndexI,IndexK>::operator-(const RIndex& i) const
+	{
+	    return lex() - i.lex();
+	}
+
+	template <class IndexI, class IndexK>
+	RIndex<IndexI,IndexK>& RIndex<IndexI,IndexK>::operator+=(Int n)
+	{
+	    *this = lex() + n;
+	    return *this;
+	}
+
+	template <class IndexI, class IndexK>
+	RIndex<IndexI,IndexK>& RIndex<IndexI,IndexK>::operator-=(Int n)
+	{
+	    *this = lex() - n;
+	    return *this;
+	}
+
+	template <class IndexI, class IndexK>
+	SizeT RIndex<IndexI,IndexK>::lex() const
+	{
+	    return mLex;
+	}
+
+	template <class IndexI, class IndexK>
+	constexpr decltype(auto) RIndex<IndexI,IndexK>::pmax() const
+	{
+	    return UPos(mK->pmax().val() * mI->pmax().val());
+	}
+
+	template <class IndexI, class IndexK>
+	constexpr decltype(auto) RIndex<IndexI,IndexK>::lmax() const
+	{
+	    return UPos(mK->lmax().val() * mI->lmax().val());
+	}
+
+	template <class IndexI, class IndexK>
+	IndexId<0> RIndex<IndexI,IndexK>::id() const
+	{
+	    return IndexId<0>(this->ptrId());
+	}
+
+	template <class IndexI, class IndexK>
+	typename RIndex<IndexI,IndexK>::MetaType RIndex<IndexI,IndexK>::operator*() const
+	{
+	    return meta();
+	}
+	
+	template <class IndexI, class IndexK>
+	constexpr SizeT RIndex<IndexI,IndexK>::dim() const
+	{
+	    return mI->dim();
+	}
+
+	template <class IndexI, class IndexK>
+	Sptr<typename RIndex<IndexI,IndexK>::RangeType> RIndex<IndexI,IndexK>::range() const
+	{
+	    return mRange;
+	}
+
+	template <class IndexI, class IndexK>
+	template <SizeT I>
+	decltype(auto) RIndex<IndexI,IndexK>::stepSize(const IndexId<I>& id) const
+	{
+	    const auto own = [&]() {
+		if constexpr(I != 0){ return SPos<0> {}; }
+		else { return UPos(id == this->id() ? 1 : 0); }
+	    };
+	    return mI->stepSize(id) + own();
+	    //return getRankStepSize(id);
+	}
+
+	template <class IndexI, class IndexK>
+	String RIndex<IndexI,IndexK>::stringMeta() const
+	{
+	    const SizeT r = mK->lex();
+	    String o;
+	    if(r == getRankNumber()){
+		o = mI->stringMeta();
+	    }
+	    bcast(o, r);
+	    return o;
+	}
+
+	template <class IndexI, class IndexK>
+	typename RIndex<IndexI,IndexK>::MetaType RIndex<IndexI,IndexK>::meta() const
+	{
+	    MetaType o;
+	    if constexpr(Typemap<MetaType>::exists){
+		const SizeT r = mK->lex();
+		if(r == getRankNumber()){
+		    o = mI->meta();
+		}
+		bcast(o, r);
+	    }
+	    else {
+		// TODO: In general don't allow RIndices if broadcast for MetaType
+		// does not exitst (once DType broadcast is implemented)!!!
+		CXZ_ERROR("no broadcast implementation for given meta type ("
+			  << typeid(MetaType).name() << ") available");
+	    }
+	    return o;
+	}
+
+	template <class IndexI, class IndexK>
+	RIndex<IndexI,IndexK>& RIndex<IndexI,IndexK>::at(const MetaType& metaPos)
+	{
+	    mI->at(metaPos);
+	    const size_t lex = mI->lex();
+	    Vector<size_t> lexs(mK->lmax().val());
+	    MPI_Allgather(&lex, 1, MPI_UNSIGNED_LONG, lexs.data(), 1, MPI_UNSIGNED_LONG,
+			  MPI_COMM_WORLD);
+	    SizeT root = 0;
+	    for(; root != lexs.size(); ++root) {
+		if(lexs[root] != mI->lmax().val()){
+		    break;
+		}
+	    }
+	    if(root == lexs.size()){ // metaPos not in rrange
+		*this = lmax().val();
+		VCHECK(toString(metaPos));
+		assert(0);
+	    }
+	    else {
+		*mK = root;
+		*mI = lexs[root];
+		(*this)();
+	    }
+	    return *this;
+	}
+
+	template <class IndexI, class IndexK>
+	RangePtr RIndex<IndexI,IndexK>::prange(const RIndex<IndexI,IndexK>& last) const
+	{
+	    CXZ_ERROR("not implemented");
+	    return nullptr;
+	}
+
+	template <class IndexI, class IndexK>
+	auto RIndex<IndexI,IndexK>::deepFormat() const
+	{
+	    return concat( mul(mK->deepFormat(), mI->lmax().val() ), mI->deepFormat() );
+	}
+
+	template <class IndexI, class IndexK>
+	auto RIndex<IndexI,IndexK>::deepMax() const
+	{
+	    return concat( mK->deepMax(), mI->deepMax() );
+	}
+
+	template <class IndexI, class IndexK>
+	RIndex<IndexI,IndexK>& RIndex<IndexI,IndexK>::reformat(const Vector<SizeT>& f, const Vector<SizeT>& s)
+	{
+	    CXZ_ERROR("not implemented");
+	    return *this;
+	}
+
+	template <class IndexI, class IndexK>
+	template <class Xpr, class F>
+	constexpr decltype(auto) RIndex<IndexI,IndexK>::ifor(const Xpr& xpr, F&& f) const
+	{
+	    return mI->ifor(xpr, std::forward<F>(f));
+	}
+
+	template <class IndexI, class IndexK>
+	bool RIndex<IndexI,IndexK>::formatIsTrivial() const
+	{
+	    return mI->formatIsTrivial();
+	}
+
+	template <class IndexI, class IndexK>
+	decltype(auto) RIndex<IndexI,IndexK>::xpr(const Sptr<RIndex<IndexI,IndexK>>& _this) const
+	{
+	    return _this->local()->xpr( _this->local() );
+	}
+
+	template <class IndexI, class IndexK>
+	RIndex<IndexI,IndexK>& RIndex<IndexI,IndexK>::operator()(const Sptr<IndexI>& i)
+	{
+	    mI = i;
+	    return (*this)();
+	}
+	
+	template <class IndexI, class IndexK>
+	RIndex<IndexI,IndexK>& RIndex<IndexI,IndexK>::operator()()
+	{
+	    if(mI->lex() >= mI->lmax().val()){
+		mLex = lmax().val();
+	    }
+	    if constexpr(has_static_sub<IndexI>::value){
+		constexpr SizeT NI = index_dim<IndexI>::value;
+		mLex = iter<0,NI>
+		    ([&](auto i) {
+			return mK->pack()[i]->lex() * mK->lexFormat()[i].val() *
+			    mI->lexFormat()[i].val() * mI->pack()[i]->lmax().val() +
+			    mI->pack()[i]->lex() * mI->lexFormat()[i].val() *
+			    mK->lexFormat()[i].val();
+		    }, [](const auto&... e) { return (e + ...); });
+	    }
+	    else if constexpr( has_static_sub<IndexK>::value){
+		constexpr SizeT NI = index_dim<IndexK>::value;
+		mLex = iter<0,NI>
+		    ([&](auto i) {
+			return mK->pack()[i]->lex() * mK->lexFormat()[i].val() *
+			    mI->lexFormat()[i].val() * mI->pack()[i]->lmax().val() +
+			    mI->pack()[i]->lex() * mI->lexFormat()[i].val() *
+			    mK->lexFormat()[i].val();
+		    }, [](const auto&... e) { return (e + ...); });
+	    }
+	    else {
+		const SizeT NI = mI->dim();
+		mLex = 0;
+		for(SizeT i = 0; i != NI; ++i){
+		    mLex += mK->pack()[i]->lex() * mK->lexFormat()[i].val() *
+			mI->lexFormat()[i].val() * mI->pack()[i]->lmax().val() +
+			mI->pack()[i]->lex() * mI->lexFormat()[i].val() *
+			mK->lexFormat()[i].val();
+		}
+	    }
+	    IB::mPos = mK->pos() * mI->pmax().val() + mI->pos();
+	    return *this;
+	}
+
+	template <class IndexI, class IndexK>
+	SizeT RIndex<IndexI,IndexK>::rank() const
+	{
+	    return mK->lex();
+	}
+
+	template <class IndexI, class IndexK>
+	Sptr<IndexI> RIndex<IndexI,IndexK>::local() const
+	{
+	    return mI;
+	}
+	
+	template <class IndexI, class IndexK>
+	Sptr<IndexK> RIndex<IndexI,IndexK>::rankI() const
+	{
+	    return mK;
+	}
+
+	template <class IndexI, class IndexK, class I1>
+	decltype(auto) operator*(const Sptr<RIndex<IndexI,IndexK>>& a, const Sptr<I1>& b)
+	{
+	    return iptrMul(a, b);
+	}
+
+	/*=====================+
+	 |    RRangeFactory    |
+	 +=====================*/
+
+	template <class RangeI, class RangeK>
+	RRangeFactory<RangeI,RangeK>::RRangeFactory(const Sptr<RangeI>& ri,
+						    const Sptr<RangeK>& rk):
+	    mRI(ri),
+	    mRK(rk)
+	{
+	    int s = 0;
+	    MPI_Comm_size(MPI_COMM_WORLD, &s);
+	    CXZ_ASSERT(rk->size() == static_cast<SizeT>(s),
+		       "geometry rank size ( = " << rk->size()
+		       << ") does not match number of ranks ( = " << s <<  ")");
+	    if constexpr(has_static_sub<typename RangeI::IndexType>::value and
+			 has_static_sub<typename RangeK::IndexType>::value) {
+		constexpr SizeT NRI = RangeI::NR;
+		constexpr SizeT NRK = RangeK::NR;
+		static_assert(NRI == NRK, "ranges have to be of same dimension");
+	    }
+	    else {
+		CXZ_ASSERT(ri->dim() == rk->dim(), "ranges have to be of same dimension, got "
+			   << ri->dim() << " and " << rk->dim());
+	    }
+	}
+
+	template <class RangeI, class RangeK>
+	void RRangeFactory<RangeI,RangeK>::make()
+	{
+	    Vector<Uuid> key = { mRI->id(), mRK->id() };
+	    const auto& info = typeid(RRange<RangeI,RangeK>);
+	    mProd = this->fromCreated(info, key);
+	    if(mProd == nullptr) {
+		mProd = std::shared_ptr<RRange<RangeI,RangeK>>
+		    ( new RRange<RangeI,RangeK>(mRI, mRK) );
+		this->addToCreated(info, key, mProd);
+	    }
+	}
+
+	/*==============+
+	 |    RRange    |
+	 +==============*/
+
+	template <class RangeI, class RangeK>
+	RangePtr RRange<RangeI,RangeK>::sub(SizeT num) const
+	{
+	    if(num == 0){
+		return mGeom;
+	    }
+	    else if(num == 1){
+		return mLocal;
+	    }
+	    else {
+		CXZ_ERROR("out of range: " << num);
+		return nullptr;
+	    }
+	}
+	
+	template <class RangeI, class RangeK>
+	MArray<RangePtr> RRange<RangeI,RangeK>::sub() const
+	{
+	    RangePtr sr = SRangeFactory<SizeT,2>( Arr<SizeT,2> { 0, 1 } ).create();
+	    return MArray<RangePtr>( sr, Vector<RangePtr>( { mGeom, mLocal } ) );
+	}
+	
+	template <class RangeI, class RangeK>
+	SizeT RRange<RangeI,RangeK>::size() const
+	{
+	    return mGeom->size() * mLocal->size();
+	}
+	
+	template <class RangeI, class RangeK>
+	SizeT RRange<RangeI,RangeK>::dim() const
+	{
+	    return 2;
+	}
+	
+	template <class RangeI, class RangeK>
+	String RRange<RangeI,RangeK>::stringMeta(SizeT pos) const
+	{
+	    return (this->begin()+pos).stringMeta();
+	}
+	
+	template <class RangeI, class RangeK>
+	const TypeInfo& RRange<RangeI,RangeK>::type() const
+	{
+	    return typeid(RRange<RangeI,RangeK>);
+	}
+	
+	template <class RangeI, class RangeK>
+	const TypeInfo& RRange<RangeI,RangeK>::metaType() const
+	{
+	    return typeid(MetaType);
+	}
+	
+	template <class RangeI, class RangeK>
+	RangePtr RRange<RangeI,RangeK>::extend(const RangePtr& r) const
+	{
+	    CXZ_ERROR("not implemented");
+	    return nullptr;
+	}
+
+	template <class RangeI, class RangeK>
+	Sptr<RangeI> RRange<RangeI,RangeK>::local() const
+	{
+	    return mLocal;
+	}
+
+	template <class RangeI, class RangeK>
+	Sptr<RangeK> RRange<RangeI,RangeK>::geom() const
+	{
+	    return mGeom;
+	}
+	    
+	template <class RangeI, class RangeK>
+	const typename RRange<RangeI,RangeK>::MetaType RRange<RangeI,RangeK>::get(SizeT pos) const
+	{
+	    return (this->begin()+pos)->meta();
+	}
+
+	template <class RangeI, class RangeK>
+	SizeT RRange<RangeI,RangeK>::getMeta(const MetaType& metaPos) const
+	{
+	    auto i = this->begin();
+	    return i.at(metaPos).pos();
+	}
+
+	template <class RangeI, class RangeK>
+	SizeT RRange<RangeI,RangeK>::getRank(SizeT pos) const
+	{
+	    return (this->begin()+pos)->rank();
+	}
+	
+	/*=========================+
+	 |    RRange (protected)   |
+	 +=========================*/
+
+	template <class RangeI, class RangeK>
+	RRange<RangeI,RangeK>::RRange(const Sptr<RangeI>& loc, const Sptr<RangeK>& geom) :
+	    mLocal(loc),
+	    mGeom(geom)
+	{}
+
+	template <class RangeI, class RangeK>
+	Vector<Uuid> RRange<RangeI,RangeK>::key() const
+	{
+	    return Vector<Uuid> { mLocal->id(), mGeom->id() };
+	}
+
+	/*============================+
+	 |    non-member functions    |
+	 +============================*/
+
+    } // namespace mpi
+
+    template <class RangeI, class RangeK>
+    Sptr<mpi::RRange<RangeI,RangeK>> RangeCast<mpi::RRange<RangeI,RangeK>>::func(const RangePtr& r)
+    {
+	CXZ_ASSERT(r->dim() == 2, "expected RRange");
+	Sptr<RangeI> loc = rangeCast<RangeI>(r->sub(1));
+	Sptr<RangeK> geom = rangeCast<RangeK>(r->sub(0));
+	return std::dynamic_pointer_cast<mpi::RRange<RangeI,RangeK>>
+	    ( mpi::RRangeFactory( loc, geom ).create() );
+    }
+	
+} // namespace CNORXZ
+
+#endif

src/opt/mpi/include/rrange.h

@@ -0,0 +1,306 @@
+// -*- C++ -*-
+/**
+   
+   @file opt/mpi/include/rrange.h
+   @brief RRange and RIndex declaration.
+   
+   Copyright (c) 2024 Christian Zimmermann. All rights reserved.
+   Mail: chizeta@f3l.de
+   
+**/
+
+#ifndef __cxz_mpi_rrange_h__
+#define __cxz_mpi_rrange_h__
+
+#include "mpi_base.h"
+
+namespace CNORXZ
+{
+    namespace mpi
+    {
+
+	/** ****
+	    Specific index for RRange.
+	    Every call on an instance of this class has to be done on each rank.
+	    @tparam IndexI Index type used to indicate local position.
+	    @tparam IndexK Index type used to indicate the rank.
+	*/
+	template <class IndexI, class IndexK>
+	class RIndex : public IndexInterface<RIndex<IndexI,IndexK>,typename IndexI::MetaType>
+	{
+	public:
+	    typedef IndexInterface<RIndex<IndexI,IndexK>,typename IndexI::MetaType> IB;
+	    typedef typename IndexI::MetaType MetaType;
+	    typedef RRange<typename IndexI::RangeType,typename IndexK::RangeType> RangeType;
+
+	    INDEX_RANDOM_ACCESS_ITERATOR_DEFS(MetaType);
+
+	    /** Default constructor. */
+	    RIndex() = default;
+
+	    /** Move constructor (default). */
+	    RIndex(RIndex&& in) = default;
+
+	    /** Move assignment (default). */
+	    RIndex& operator=(RIndex&& in) = default;
+
+	    /** Copy constructor (no default, copy local index instance). */
+	    RIndex(const RIndex& in);
+
+	    /** Copy assignment (no default, copy local index instance). */
+	    RIndex& operator=(const RIndex& in);
+
+	    /** Construct from global range and format.
+		@param global Pointer to global range (RRange).
+		@param lexpos Start position.
+	     */
+	    RIndex(const RangePtr& global, SizeT lexpos = 0);
+
+	    /** Construct from local index and rank index.
+		@param i Local index.
+		@param k Rank index.
+	     */
+	    RIndex(const Sptr<IndexI>& i, const Sptr<IndexK>& k);
+	    
+	    /** @copydoc IndexInterface::operator=(SizeT) */
+	    RIndex& operator=(SizeT pos);
+
+	    /** @copydoc IndexInterface::operator++() */
+	    RIndex& operator++();
+
+	    /** @copydoc IndexInterface::operator--() */
+	    RIndex& operator--();
+
+	    /** @copydoc IndexInterface::operator+() */
+	    RIndex operator+(Int n) const;
+
+	    /** @copydoc IndexInterface::operator-() */
+	    RIndex operator-(Int n) const;
+
+	    /** @copydoc IndexInterface::operator-(CIndex) */
+	    SizeT operator-(const RIndex& i) const;
+
+	    /** @copydoc IndexInterface::operator+=() */
+	    RIndex& operator+=(Int n);
+
+	    /** @copydoc IndexInterface::operator-=() */
+	    RIndex& operator-=(Int n);
+
+	    /** @copydoc IndexInterface::lex() */
+	    SizeT lex() const;
+
+	    /** @copydoc IndexInterface::pmax() */
+	    constexpr decltype(auto) pmax() const;
+
+	    /** @copydoc IndexInterface::lmax() */
+	    constexpr decltype(auto) lmax() const;
+
+	    /** @copydoc IndexInterface::id() */
+	    IndexId<0> id() const;
+
+	    /** @copydoc IndexInterface::operator*() */
+	    MetaType operator*() const;
+	
+	    /** @copydoc IndexInterface::dim() */
+	    constexpr SizeT dim() const;
+
+	    /** @copydoc IndexInterface::range() */
+	    Sptr<RangeType> range() const;
+
+	    /** @copydoc IndexInterface::stepSize() */
+	    template <SizeT I>
+	    decltype(auto) stepSize(const IndexId<I>& id) const;
+
+	    /** @copydoc IndexInterface::stringMeta() */
+	    String stringMeta() const;
+
+	    /** @copydoc IndexInterface::meta() */
+	    MetaType meta() const;
+
+	    /** @copydoc IndexInterface::at() */
+	    RIndex& at(const MetaType& metaPos);
+
+	    /** @copydoc IndexInterface::prange() */
+	    RangePtr prange(const RIndex<IndexI,IndexK>& last) const;
+
+	    /** @copydoc IndexInterface::deepFormat() */
+	    auto deepFormat() const;
+
+	    /** @copydoc IndexInterface::deepMax() */
+	    auto deepMax() const;
+
+	    /** @copydoc IndexInterface::reformat() */
+	    RIndex& reformat(const Vector<SizeT>& f, const Vector<SizeT>& s);
+
+	    /** @copydoc IndexInterface::ifor() */
+	    template <class Xpr, class F>
+	    constexpr decltype(auto) ifor(const Xpr& xpr, F&& f) const;
+
+	    /** @copydoc IndexInterface::formatIsTrivial() */
+	    bool formatIsTrivial() const;
+
+	    /** @copydoc IndexInterface::xpr() */
+	    decltype(auto) xpr(const Sptr<RIndex<IndexI,IndexK>>& _this) const;
+
+	    /** Replace local index instance and update index position correspondingly.
+		@param i New index instances.
+	     */
+	    RIndex& operator()(const Sptr<IndexI>& i);
+
+	    /** Update index position according to the sub-indices. */
+	    RIndex& operator()();
+	    
+	    /** Get the current rank. */
+	    SizeT rank() const;
+
+	    /** Get the local index on THIS rank. */
+	    Sptr<IndexI> local() const;
+
+	    /** Get index indicating the current rank this index points to. */
+	    Sptr<IndexK> rankI() const;
+	    
+	private:
+	    SizeT mLex = 0;
+	    Sptr<RangeType> mRange; /**< RRange. */
+	    Sptr<IndexI> mI; /**< Index on the local range of the THIS rank. */
+	    Sptr<IndexK> mK; /**< Multi-index indicating the current rank. */
+
+	};
+
+	template <class IndexI, class IndexK>
+	constexpr decltype(auto) rindex(const Sptr<IndexI>& i, const Sptr<IndexK>& k)
+	{
+	    return RIndex<IndexI,IndexK>(i,k);
+	}
+
+	template <class IndexI, class IndexK>
+	constexpr decltype(auto) rindexPtr(const Sptr<IndexI>& i, const Sptr<IndexK>& k)
+	{
+	    return std::make_shared<RIndex<IndexI,IndexK>>(i,k);
+	}
+
+	/** Make index pack of a RIndex and another index.
+	    @param a pointer to RIndex.
+	    @param b pointer to another index.
+	*/
+	template <class IndexI, class IndexK, class I1>
+	decltype(auto) operator*(const Sptr<RIndex<IndexI,IndexK>>& a, const Sptr<I1>& b);
+
+	// Traits!!!
+	template <class I>
+	struct is_rank_index
+	{
+	    static constexpr bool value = false;
+	};
+
+	template <class IndexI, class IndexK>
+	struct is_rank_index<RIndex<IndexI,IndexK>>
+	{
+	    static constexpr bool value = true;
+	};
+	
+	/** ****
+	    Specific factory for RRange.
+	    @tparam RangeI Local range type.
+	    @tparam RangeK Geometry range type.
+	 */
+	template <class RangeI, class RangeK>
+	class RRangeFactory : public RangeFactoryBase
+	{
+	public:
+	    /** Construct and setup factory.
+		@param ri Local range.
+		@param rk Geometry range.
+	     */
+	    RRangeFactory(const Sptr<RangeI>& ri, const Sptr<RangeK>& rk);
+	    
+	private:
+	    RRangeFactory() = default;
+	    virtual void make() override final;
+
+	    Sptr<RangeI> mRI;
+	    Sptr<RangeK> mRK;
+	};	
+
+	/** ****
+	    Range-Wrapper for ranges that are distributed on MPI ranks.
+	    @tparam RangeI Local range type.
+	    @tparam RangeK Geometry range type.
+	*/
+	template <class RangeI, class RangeK>
+	class RRange : public RangeInterface<RRange<RangeI,RangeK>>
+	{
+	public:
+	    typedef RangeBase RB;
+	    typedef RIndex<typename RangeI::IndexType,typename RangeK::IndexType> IndexType;
+	    typedef typename RangeI::MetaType MetaType;
+
+	    friend RRangeFactory<RangeI,RangeK>;
+
+	    virtual RangePtr sub(SizeT num) const override final;
+	    virtual MArray<RangePtr> sub() const override final;
+	    virtual SizeT size() const override final;
+	    virtual SizeT dim() const override final;
+	    virtual String stringMeta(SizeT pos) const override final;
+	    virtual const TypeInfo& type() const override final;
+	    virtual const TypeInfo& metaType() const override final;
+	    virtual RangePtr extend(const RangePtr& r) const override final;
+
+	    /** Get local range. */
+	    Sptr<RangeI> local() const;
+
+	    /** Get range of the rank geometry. */
+	    Sptr<RangeK> geom() const;
+	    
+	    /** Get meta data for given lexicographic position.
+		@param pos Lexicographic position.
+	    */
+	    const MetaType get(SizeT pos) const;
+
+	    /** Get lexicographic position according to the given meta data value.
+		@param metaPos Meta data value.
+	    */
+	    SizeT getMeta(const MetaType& metaPos) const;
+
+	    /** Get rank from lexicographic meta data position.
+		@param pos Lexicographic meta data position.
+	    */
+	    SizeT getRank(SizeT pos) const;
+
+	protected:
+
+	    /** Dafault constructor */
+	    RRange() = default;
+
+	    RRange(const RRange& in) = delete;
+	    RRange& operator=(const RRange& in) = delete;
+
+	    /** Construct from local range and geometry.
+		@param loc Local range.
+		@param geom Rank geometry range.
+	    */
+	    RRange(const Sptr<RangeI>& loc, const Sptr<RangeK>& geom);
+	
+	    Sptr<RangeI> mLocal; /**< Local range of THIS rank. */
+	    Sptr<RangeK> mGeom; /**< Rank geometry range. */
+	
+	    virtual Vector<Uuid> key() const override final;
+	};
+
+	/** Create RRange from global range and given rank geometry.
+	    @param global Global range.
+	    @param geom Rank geometry.
+	*/
+	RangePtr rrange(const RangePtr& global, const RangePtr& geom);
+
+    } // namespace mpi
+
+    template <class RangeI, class RangeK>
+    struct RangeCast<mpi::RRange<RangeI,RangeK>>
+    {
+	static Sptr<mpi::RRange<RangeI,RangeK>> func(const RangePtr& r);
+    };
+	
+} // namespace CNORXZ
+
+#endif

src/opt/mpi/include/typemap.h

@@ -0,0 +1,87 @@
+
+#ifndef __cxz_mpi_typemap_h__
+#define __cxz_mpi_typemap_h__
+
+#include "mpi.h"
+
+namespace CNORXZ
+{
+    namespace mpi
+    {
+	template <typename T>
+	struct Typemap
+	{
+	    static constexpr bool exists = false;
+	    static MPI_Datatype value() { return MPI_BYTE; }
+	};
+
+	template <>
+	struct Typemap<int>
+	{
+	    static constexpr bool exists = true;
+	    static MPI_Datatype value() { return MPI_INT; }
+	};
+
+	template <>
+	struct Typemap<unsigned>
+	{
+	    static constexpr bool exists = true;
+	    static MPI_Datatype value() { return MPI_UNSIGNED; }
+	};
+
+	template <>
+	struct Typemap<long int>
+	{
+	    static constexpr bool exists = true;
+	    static MPI_Datatype value() { return MPI_LONG; }
+	};
+
+	template <>
+	struct Typemap<unsigned long>
+	{
+	    static constexpr bool exists = true;
+	    static MPI_Datatype value() { return MPI_UNSIGNED_LONG; }
+	};
+
+	template <>
+	struct Typemap<double>
+	{
+	    static constexpr bool exists = true;
+	    static MPI_Datatype value() { return MPI_DOUBLE; }
+	};
+
+	template <>
+	struct Typemap<float>
+	{
+	    static constexpr bool exists = true;
+	    static MPI_Datatype value() { return MPI_FLOAT; }
+	};
+
+	template <typename T, SizeT N>
+	struct Typemap<Arr<T,N>>
+	{
+	    static constexpr bool exists = Typemap<T>::exists;
+	    static MPI_Datatype value() { return Typemap<T>::value(); }
+	};
+
+	template <typename T>
+	struct Typemap<Vector<T>>
+	{
+	    static constexpr bool exists = Typemap<T>::exists;
+	    static MPI_Datatype value() { return Typemap<T>::value(); }
+	};
+
+	template <typename... Ts>
+	struct Typemap<Tuple<Ts...>>
+	{
+	    static constexpr bool exists = ( Typemap<Ts>::exists and ... );
+	    static MPI_Datatype value() { return MPI_BYTE; }
+	};
+	
+	
+	// further !!!
+	
+    } // namespace mpi
+} // namespace CNORXZ
+
+#endif

src/opt/mpi/lib/CMakeLists.txt

@@ -0,0 +1,18 @@
+
+set(libcnorxzmpi_a_SOURCES
+  ${CMAKE_CURRENT_SOURCE_DIR}/mpi_base.cc
+  ${CMAKE_CURRENT_SOURCE_DIR}/rrange.cc
+  )
+
+add_library(cnorxzmpi_obj OBJECT
+  ${libcnorxzmpi_a_SOURCES}
+  )
+
+set_target_properties(cnorxzmpi_obj PROPERTIES POSITION_INDEPENDENT_CODE TRUE)
+
+add_library(cnorxzmpi SHARED
+  $<TARGET_OBJECTS:cnorxzmpi_obj>
+  )
+set_target_properties(cnorxzmpi PROPERTIES POSITION_INDEPENDENT_CODE TRUE)
+
+install(TARGETS cnorxzmpi LIBRARY DESTINATION ${INSTALL_PATH}/lib)

src/opt/mpi/lib/mpi_base.cc

@@ -0,0 +1,38 @@
+// -*- C++ -*-
+
+/**
+   
+   @file
+   @brief
+
+   Copyright (c) 2024 Christian Zimmermann. All rights reserved.
+   Mail: chizeta@f3l.de
+   
+**/
+
+#include "mpi_base.h"
+
+namespace CNORXZ
+{
+    namespace mpi
+    {
+
+	SizeT getRankNumber()
+	{
+	    int s;
+	    const int ret = MPI_Comm_rank(MPI_COMM_WORLD, &s);
+	    CXZ_ASSERT(ret == MPI_SUCCESS, "got mpi error = " << ret);
+	    return static_cast<SizeT>(s);
+	}
+
+	SizeT getNumRanks()
+	{
+	    int s;
+	    const int ret = MPI_Comm_size(MPI_COMM_WORLD, &s);
+	    CXZ_ASSERT(ret == MPI_SUCCESS, "got mpi error = " << ret);
+	    return static_cast<SizeT>(s);
+	}
+	
+    } // namespace mpi
+} // namespace CNORXZ
+

src/opt/mpi/lib/rrange.cc

@@ -0,0 +1,80 @@
+// -*- C++ -*-
+/**
+   
+   @file opt/mpi/lib/rrange.h
+   @brief RRange and RIndex implementation.
+   
+   Copyright (c) 2024 Christian Zimmermann. All rights reserved.
+   Mail: chizeta@f3l.de
+   
+**/
+
+#include "rrange.h"
+#include "rrange.cc.h"
+
+namespace CNORXZ
+{
+    namespace mpi
+    {
+
+	/*===========================+
+	 |    non-member functions   |
+	 +===========================*/
+
+	MArray<RangePtr> rsplit(const RangePtr& global, const RangePtr& geom)
+	{
+	    if(geom->dim() > 1){
+		CXZ_ASSERT(global->dim() == geom->dim(),
+			   "global range has to be of same dimension as geometry range, got "
+			   << global->dim() << " and " << geom->dim());
+		MArray<RangePtr> o(geom);
+		YIndex k(geom);
+		YIndex ke(geom, geom->size());
+		auto mu = std::make_shared<CIndex>(geom->sub().range());
+		MArray<MArray<RangePtr>> rs(geom->sub().range());
+		rs(mu) = operation( [&](const SizeT i) { return rsplit(global->sub(i), geom->sub(i)); } , xpr(mu) );
+		Vector<RangePtr> elem(geom->dim());
+		for(; k != ke; ++k){
+		    mu->ifor( operation( [&](const SizeT i, const MArray<RangePtr>& el){
+			elem[i] = el[DIndex(k.pack()[i])];
+		    }, xpr(mu), rs(mu) ), NoF {} )();
+		    o[k] = YRangeFactory(elem).create();
+		}
+		return o;
+	    }
+	    else {
+		CXZ_ASSERT(global->size() % geom->size() == 0,
+			   "global range has to be dividible by geometry range, got "
+			   << global->size() << " and " << geom->size());
+		const SizeT n = global->size() / geom->size();
+		auto k = std::make_shared<CIndex>(geom);
+		auto jb = global->begin();
+		auto je = global->begin();
+		MArray<RangePtr> o(geom);
+		o(k) = operation( [&](const SizeT x){
+		    jb = n*x; je = n*(x+1)-1; return jb.prange(je);
+		} , xpr(k) );
+		return o;
+	    }
+	}
+	
+	RangePtr rrange(const RangePtr& global, const RangePtr& geom)
+	{
+	    const MArray<RangePtr> rs = rsplit(global, geom);
+	    RangePtr o = nullptr;
+	    for(auto ri = rs.begin(); ri != rs.end(); ++ri){
+		if(ri.lex() == getRankNumber()){
+		    o = *ri;
+		}
+	    }
+	    assert(o);
+	    auto loc = rangeCast<YRange>(o);
+	    auto geo = rangeCast<YRange>(geom);
+	    RRangeFactory<YRange,YRange> xx(loc, geo);
+	    return RRangeFactory<YRange,YRange>(loc, geo).create();
+	}
+	
+	
+    } // namespace mpi
+} // namespace CNORXZ
+	

src/opt/mpi/tests/CMakeLists.txt

@@ -0,0 +1,27 @@
+
+add_definitions(-DTEST_NUMBER_FILE="${CMAKE_SOURCE_DIR}/src/tests/numbers.txt")
+include_directories(${CMAKE_SOURCE_DIR}/src/tests)
+
+add_executable(mpirrutest rrange_unit_test.cc)
+add_dependencies(mpirrutest cnorxz cnorxzmpi test_lib)
+target_link_libraries(mpirrutest ${GTEST_BOTH_LIBRARIES} ${CMAKE_THREAD_LIBS_INIT} ${MPI_LIBS} cnorxz cnorxzmpi test_lib)
+set(MPI_TEST_COMMAND mpirun -n 4 mpirrutest)
+add_test(NAME mpirrutest COMMAND ${MPI_TEST_COMMAND})
+
+add_executable(mpirautest rarray_unit_test.cc)
+add_dependencies(mpirautest cnorxz cnorxzmpi test_lib)
+target_link_libraries(mpirautest ${GTEST_BOTH_LIBRARIES} ${CMAKE_THREAD_LIBS_INIT} ${MPI_LIBS} cnorxz cnorxzmpi test_lib)
+set(MPI_TEST_COMMAND mpirun -n 4 mpirautest)
+add_test(NAME mpirautest COMMAND ${MPI_TEST_COMMAND})
+
+add_executable(mpisbutest setbuf_unit_test.cc)
+add_dependencies(mpisbutest cnorxz cnorxzmpi test_lib)
+target_link_libraries(mpisbutest ${GTEST_BOTH_LIBRARIES} ${CMAKE_THREAD_LIBS_INIT} ${MPI_LIBS} cnorxz cnorxzmpi test_lib)
+set(MPI_TEST_COMMAND mpirun -n 4 mpisbutest)
+add_test(NAME mpisbutest COMMAND ${MPI_TEST_COMMAND})
+
+add_executable(mpiroptest roperation_unit_test.cc)
+add_dependencies(mpiroptest cnorxz cnorxzmpi test_lib)
+target_link_libraries(mpiroptest ${GTEST_BOTH_LIBRARIES} ${CMAKE_THREAD_LIBS_INIT} ${MPI_LIBS} cnorxz cnorxzmpi test_lib)
+set(MPI_TEST_COMMAND mpirun -n 4 mpiroptest)
+add_test(NAME mpiroptest COMMAND ${MPI_TEST_COMMAND})

src/opt/mpi/tests/mpi_env.h

@@ -0,0 +1,47 @@
+// -*- C++ -*-
+/**
+   
+   @file opt/mpi/tests/mpi_env.h
+   @brief MPI Environment.
+
+   Copyright (c) 2024 Christian Zimmermann. All rights reserved.
+   Mail: chizeta@f3l.de
+   
+**/
+
+#ifndef __cxz_mpi_mpi_env_h__
+#define __cxz_mpi_mpi_env_h__
+
+#include <cstdlib>
+#include "gtest/gtest.h"
+
+namespace CNORXZ
+{
+    namespace mpi
+    {
+	class MPIEnv : public ::testing::Environment
+	{
+	public:
+
+	    MPIEnv(int argc, char** argv) : mArgc(argc), mArgv(argv) {}
+	
+	    virtual ~MPIEnv() override {}
+
+	    virtual void SetUp() override
+	    {
+		MPI_Init(&mArgc, &mArgv);
+	    }
+
+	    virtual void TearDown() override
+	    {
+		MPI_Finalize();
+	    }
+
+	protected:
+	    int mArgc;
+	    char** mArgv;
+	};
+    }
+}
+
+#endif

src/opt/mpi/tests/rarray_unit_test.cc

@@ -0,0 +1,108 @@
+// -*- C++ -*-
+/**
+   
+   @file opt/mpi/tests/rarray_unit_test.cc
+   @brief RArray unit tests.
+
+   Copyright (c) 2024 Christian Zimmermann. All rights reserved.
+   Mail: chizeta@f3l.de
+   
+**/
+
+#include <cstdlib>
+#include <iostream>
+
+#include "gtest/gtest.h"
+
+#include "cnorxz.h"
+#include "cnorxz_mpi.h"
+#include "test_numbers.h"
+#include "rrange.cc.h"
+#include "rarray.cc.h"
+#include "mpi_env.h"
+
+namespace
+{
+    using namespace CNORXZ;
+    using Test::Numbers;
+    using namespace CNORXZ::mpi;
+
+    class RCArray_Test : public ::testing::Test
+    {
+    protected:
+
+	RCArray_Test()
+	{
+	    CXZ_ASSERT(getNumRanks() == 4, "exptected 4 ranks");
+	    Vector<Int> xs(12);
+	    Vector<Int> ts(16);
+	    for(SizeT i = 0; i != xs.size(); ++i){
+		const Int x = static_cast<Int>(i) - static_cast<Int>(xs.size()/2);
+		xs[i] = x;
+	    }
+	    for(SizeT i = 0; i != ts.size(); ++i){
+		const Int t = static_cast<Int>(i) - static_cast<Int>(ts.size()/2);
+		ts[i] = t;
+	    }
+	    mXRange = URangeFactory<Int>(xs).create();
+	    mTRange = URangeFactory<Int>(ts).create();
+	    Vector<RangePtr> rs { mTRange, mXRange, mXRange, mXRange };
+	    mGRange = YRangeFactory(rs).create();
+	    RangePtr g1 = CRangeFactory(1).create();
+	    RangePtr g2 = CRangeFactory(2).create();
+	    Vector<RangePtr> gs { g2, g1, g1, g2 };
+	    mGeom = YRangeFactory(gs).create();
+	    mRRange = rrange(mGRange, mGeom);
+	    const SizeT size = ts.size()*xs.size()*xs.size()*xs.size();
+	    const SizeT locsize = size/4;
+	    Vector<Double> vec = Numbers::get(0,locsize+10);
+	    Vector<Double> data(locsize);
+	    mData.resize(size);
+	    const SizeT myrank = getRankNumber();
+	    for(SizeT i = 0; i != locsize; ++i){
+		assert(i < data.size());
+		data[i] = vec[i] * vec[i+myrank] / vec[i+2*myrank];
+		assert(i + locsize*myrank < mData.size());
+		mData[i + locsize*myrank] = data[i];
+	    }
+	    MPI_Barrier(MPI_COMM_WORLD);
+	    for(SizeT r = 0; r != 4; ++r){
+		MPI_Bcast(mData.data() + locsize*r, locsize, MPI_DOUBLE, r, MPI_COMM_WORLD);
+		MPI_Barrier(MPI_COMM_WORLD);
+	    }
+	    mLoc = MArray<Double>( mRRange->sub(1), data);
+	    mA = RCArray<Double>(mLoc, mGeom);
+	}
+
+	RangePtr mXRange;
+	RangePtr mTRange;
+	RangePtr mGRange;
+	RangePtr mGeom;
+	RangePtr mRRange;
+        MArray<Double> mLoc;
+	RCArray<Double> mA;
+	Vector<Double> mData;
+    };
+
+    TEST_F(RCArray_Test, Basics)
+    {
+	EXPECT_EQ(mA.size(), mRRange->size());
+    }
+
+    TEST_F(RCArray_Test, GlobalIterate)
+    {
+	const SizeT size = mRRange->sub(1)->size();
+	auto e = mA.end();
+	for(auto i = mA.begin(); i != e; ++i){
+	    const Double x = *i;
+	    EXPECT_EQ(x, mData[i.rank()*size + i.local()->pos()]);
+	}
+    }
+}
+
+int main(int argc, char** argv)
+{
+    ::testing::InitGoogleTest(&argc, argv);
+    ::testing::AddGlobalTestEnvironment( new MPIEnv(argc, argv) );
+    return RUN_ALL_TESTS();
+}

src/opt/mpi/tests/roperation_unit_test.cc

@@ -0,0 +1,239 @@
+// -*- C++ -*-
+/**
+   
+   @file opt/mpi/tests/roperation_unit_test.cc
+   @brief Rank operation unit tests.
+
+   Copyright (c) 2024 Christian Zimmermann. All rights reserved.
+   Mail: chizeta@f3l.de
+   
+**/
+
+#include <cstdlib>
+#include <iostream>
+
+#include "gtest/gtest.h"
+
+#include "cnorxz.h"
+#include "cnorxz_mpi.h"
+#include "test_numbers.h"
+#include "rrange.cc.h"
+#include "rarray.cc.h"
+#include "mpi_env.h"
+
+namespace
+{
+    using namespace CNORXZ;
+    using Test::Numbers;
+    using namespace CNORXZ::mpi;
+
+    class ROp_Test : public ::testing::Test
+    {
+    protected:
+
+	ROp_Test()
+	{
+	    CXZ_ASSERT(getNumRanks() == 4, "exptected 4 ranks");
+	    const SizeT myrank = getRankNumber();
+	    T = 16;
+	    L = 12;
+	    mTRange = CRangeFactory(T).create();
+	    mLRange = CRangeFactory(L).create();
+	    mSpRange = SRangeFactory(Arr<SizeT,4>{0,1,2,3}).create();
+	    mXRange = mTRange*mLRange*mLRange*mLRange;
+	    
+	    RangePtr g1 = CRangeFactory(1).create();
+	    RangePtr g2 = CRangeFactory(2).create();
+	    mGeom = g2*g1*g1*g2*g1;
+	    mGeom2 = g2*g1*g1*g2;
+	    
+	    mRXRange = rrange(mXRange, g2*g1*g1*g2);
+	    RangePtr scr = mSpRange*mSpRange;
+	    const Vector<Double> vec = Numbers::get(0,mRXRange->sub(1)->size()+2);
+	    RangePtr ltr = mRXRange->sub(1)->sub(0);
+	    RangePtr ll1r = mRXRange->sub(1)->sub(1);
+	    RangePtr ll2r = mRXRange->sub(1)->sub(2);
+	    RangePtr ll3r = mRXRange->sub(1)->sub(3);
+	    mMRange = ltr*ll1r*ll2r*ll3r*scr;
+	    Vector<Double> data(mMRange->size());
+	    Vector<Double> data2(mMRange->size());
+	    for(SizeT i = 0; i != mRXRange->sub(1)->size(); ++i){
+		for(SizeT j = 0; j != scr->size(); ++j){
+		    const SizeT k = i*scr->size() + j; 
+		    data[k] = vec[i] * static_cast<Double>(j+2);
+		    data2[k] = vec[i] / static_cast<Double>(j+2);
+		    if(k > 0){ 
+			assert(data[k] != data[k-1]);
+		    }
+		}
+	    }
+	    mM1 = RCArray<Double>( MArray<Double>(mMRange, data), mGeom );
+	    mM2 = RCArray<Double>( MArray<Double>(mMRange, data2), mGeom );
+	    mAll1 = Vector<Double>(data.size() * getNumRanks());
+	    mAll2 = Vector<Double>(data2.size() * getNumRanks());
+	    typedef RIndex<MIndex<CIndex,CIndex,CIndex,CIndex>,MIndex<CIndex,CIndex,CIndex,CIndex>> RI;
+	    const SizeT scrs = scr->size();
+	    auto rix = RI(mRXRange);
+	    assert(rix.lmax().val() == 27648);
+	    assert(scrs == 16);
+	    for(auto ri = RI(mRXRange); ri.lex() != ri.lmax().val(); ++ri){
+		Double* buf = mAll1.data() + scrs*ri.lex();
+		Double* buf2 = mAll2.data() + scrs*ri.lex();
+		if(ri.rank() == myrank){
+		    std::memcpy(buf, data.data()+ri.local()->lex()*scrs, scrs*sizeof(Double));
+		    std::memcpy(buf2, data2.data()+ri.local()->lex()*scrs, scrs*sizeof(Double));
+		}
+		MPI_Bcast(buf, scrs, MPI_DOUBLE, ri.rank(), MPI_COMM_WORLD);
+		MPI_Bcast(buf2, scrs, MPI_DOUBLE, ri.rank(), MPI_COMM_WORLD);
+	    }
+	}
+
+	SizeT T;
+	SizeT L;
+	RangePtr mTRange;
+	RangePtr mLRange;
+	RangePtr mXRange; // global T*L*L*L
+	RangePtr mRXRange;
+	RangePtr mSpRange;
+	RangePtr mMRange;
+	RangePtr mGeom;
+	RangePtr mGeom2;
+
+	RCArray<Double> mM1;
+	RCArray<Double> mM2;
+	Vector<Double> mAll1;
+	Vector<Double> mAll2;
+    };
+
+    TEST_F(ROp_Test, Check)
+    {
+	EXPECT_EQ(mM1.size(), mM2.size());
+    }
+
+    TEST_F(ROp_Test, Difference)
+    {
+	RArray<Double> res( MArray<Double>(mM1.range()->sub(1)), mGeom );
+	Vector<Double> comp( mXRange->size()*mSpRange->size()*mSpRange->size() );
+	EXPECT_EQ(res.size(), comp.size());
+
+	typedef UIndex<SizeT> UI;
+	
+	auto xp = std::make_shared<RIndex<MIndex<UI,UI,UI,UI>,MIndex<UI,UI,UI,UI>>>(mRXRange);
+	auto xm = std::make_shared<RIndex<MIndex<UI,UI,UI,UI>,MIndex<UI,UI,UI,UI>>>(mRXRange);
+	auto x = std::make_shared<RIndex<MIndex<UI,UI,UI,UI>,MIndex<UI,UI,UI,UI>>>(mRXRange);
+	auto A = std::make_shared<SIndex<SizeT,4>>(mSpRange);
+	auto B = std::make_shared<SIndex<SizeT,4>>(mSpRange);
+	auto AB = mindexPtr(A*B);
+	
+	Sptr<Vector<SizeT>> imap1;
+	imap1 = setupMap(xp, x, [&](const auto& vec) {
+	    return std::make_tuple((std::get<0>(vec)+1)%T, (std::get<1>(vec)+1)%L,
+				   (std::get<2>(vec)+1)%L, (std::get<3>(vec)+1)%L); } );
+	Vector<bool> req(xp->range()->size(), false);
+	for(const auto& r: *imap1){
+	    req[r] = true;
+	}
+	mM1.load(xp, AB, req);
+	res.rop(x*A*B) = mapXpr(xp,x,imap1, mM1(xp*A*B) - mM1(x*A*B) );
+
+	for(SizeT x0 = 0; x0 != T; ++x0) { 
+	for(SizeT x1 = 0; x1 != L; ++x1)
+	for(SizeT x2 = 0; x2 != L; ++x2)
+	for(SizeT x3 = 0; x3 != L; ++x3)
+	for(SizeT A = 0; A != 4; ++A)
+	for(SizeT B = 0; B != 4; ++B) {
+	    const SizeT xi = x0*L*L*L + x1*L*L + x2*L + x3;
+	    const SizeT x0p = (x0+1)%T;
+	    const SizeT x1p = (x1+1)%L;
+	    const SizeT x2p = (x2+1)%L;
+	    const SizeT x3p = (x3+1)%L;
+	    const SizeT xpi = x0p*L*L*L + x1p*L*L + x2p*L + x3p;
+	    const SizeT pi = xpi*4*4 + A*4 + B;
+	    const SizeT ri = xi*4*4 + A*4 + B;
+	    comp[ri] = mAll1[pi] - mAll1[ri];
+	}}
+	
+	for(auto i = res.begin(); i.lex() != i.lmax().val(); ++i){
+	    const auto a1 = *i;
+	    const auto a2 = comp[i.lex()];
+	    EXPECT_EQ(a1, a2);
+	}
+    }
+
+    TEST_F(ROp_Test, Contract)
+    {
+	Vector<Double> comp(mRXRange->size());
+	RArray<Double> res( MArray<Double>(mRXRange->sub(1)), mGeom2 );
+	EXPECT_EQ(res.size(), comp.size());
+	typedef UIndex<SizeT> UI;
+	auto xl = std::make_shared<RIndex<MIndex<UI,UI,UI,UI>,MIndex<UI,UI,UI,UI>>>(mRXRange);
+	auto y = std::make_shared<RIndex<MIndex<UI,UI,UI,UI>,MIndex<UI,UI,UI,UI>>>(mRXRange);
+	auto xy = std::make_shared<RIndex<MIndex<UI,UI,UI,UI>,MIndex<UI,UI,UI,UI>>>(mRXRange);
+	auto x = std::make_shared<RIndex<MIndex<UI,UI,UI,UI>,MIndex<UI,UI,UI,UI>>>(mRXRange);
+	auto A = std::make_shared<SIndex<SizeT,4>>(mSpRange);
+	auto B = std::make_shared<SIndex<SizeT,4>>(mSpRange);
+
+	Sptr<Vector<SizeT>> imap1;
+	Sptr<Vector<SizeT>> imap2;
+	auto AB = mindexPtr(A*B);
+	auto BA = mindexPtr(B*A);
+	auto myx = mindexPtr(y*xl);
+	
+	// block 1:
+	//imap1 = setupMap(x, xl, [](const auto& vec) { return vec; } );
+	//imap2 = setupMap(xy, myx, [](const auto& vec) { return std::get<0>(vec); } );
+	//mM1.load(xl, x, AB, imap1);
+	//mM2.load(myx, xy, BA, imap2);
+	//res(y).a(mindexPtr(y*xl), [](auto& t, const auto& s) { t += s; },
+	//	 mapXpr( x, xl, imap1, mapXpr( xy, myx, imap2, (mM1(x*A*B) * mM2(xy*B*A)).c(AB) ) ) );
+	/*
+	// block 2:
+	mM1.load(,A*B*a*b);
+	mM2.load(,A*B*a*b);
+	res(y) += (mM1(x*A*B*a*b) * mM2(xy*B*A*b*a)).c(mindexPtr(x*A*B));
+	// block 3:
+	mM1.load(,A*B*a*b);
+	mM2.load(,A*B*a*b);
+	res(y) += (mM1(x*A*B*a*b) * mM2(xy*B*A*b*a)).c(x*A*B*a*b);
+	// block 4:
+	mM1.load(,A*B*a*b);
+	mM2.load(,A*B*a*b);
+	res(y) += (mM1(x*A*B*a*b) * mM2(xy*B*A*b*a)).c(x*A*B*a*b);
+	// comparison loop
+	for(SizeT x0 = 0; x0 != T; ++x0) { VCHECK(x0);
+	for(SizeT x1 = 0; x1 != L; ++x1)
+	for(SizeT x2 = 0; x2 != L; ++x2)
+	for(SizeT x3 = 0; x3 != L; ++x3)
+	for(SizeT y0 = 0; y0 != T; ++y0)
+	for(SizeT y1 = 0; y1 != L; ++y1)
+	for(SizeT y2 = 0; y2 != L; ++y2)
+	for(SizeT y3 = 0; y3 != L; ++y3)
+	for(SizeT A = 0; A != 4; ++A)
+	for(SizeT B = 0; B != 4; ++B) {
+	    const SizeT yi = y0*L*L*L + y1*L*L + y2*L + y3;
+	    const SizeT xi = x0*L*L*L + x1*L*L + x2*L + x3;
+	    const SizeT xy0 = (x0 + y0) % T;
+	    const SizeT xy1 = (x1 + y1) % L;
+	    const SizeT xy2 = (x2 + y2) % L;
+	    const SizeT xy3 = (x3 + y3) % L;
+	    const SizeT xyi = xy0*L*L*L + xy1*L*L + xy2*L + xy3;
+	    const SizeT i1 = xi*4*4 + A*4 + B;
+	    const SizeT i2 = xyi*4*4 + B*4 + A;
+	    comp[yi] += mAll1[i1] * mAll2[i2];
+	}}
+	VCHECK(comp[123]);
+	*/
+	/*
+	for(auto i = res.begin(); i.lex() != i.lmax().val(); ++i){
+	    EXPECT_EQ(*i, comp[i.lex()]);
+	}
+	*/
+    }
+}
+
+int main(int argc, char** argv)
+{
+    ::testing::InitGoogleTest(&argc, argv);
+    ::testing::AddGlobalTestEnvironment( new MPIEnv(argc, argv) );
+    return RUN_ALL_TESTS();
+}

src/opt/mpi/tests/rrange_unit_test.cc

@@ -0,0 +1,133 @@
+// -*- C++ -*-
+/**
+   
+   @file opt/mpi/tests/rrange_unit_test.cc
+   @brief RRange unit tests.
+
+   Copyright (c) 2024 Christian Zimmermann. All rights reserved.
+   Mail: chizeta@f3l.de
+   
+**/
+
+#include <cstdlib>
+#include <iostream>
+
+#include "gtest/gtest.h"
+
+#include "cnorxz.h"
+#include "cnorxz_mpi.h"
+#include "test_numbers.h"
+#include "rrange.cc.h"
+#include "mpi_env.h"
+
+namespace
+{
+    using namespace CNORXZ;
+    using Test::Numbers;
+    using namespace CNORXZ::mpi;
+    
+    class RRange_Test : public ::testing::Test
+    {
+    protected:
+
+	RRange_Test()
+	{
+	    CXZ_ASSERT(getNumRanks() == 4, "exptected 4 ranks");
+	    Vector<Int> xs(12);
+	    Vector<Int> ts(16);
+	    for(SizeT i = 0; i != xs.size(); ++i){
+		const Int x = static_cast<Int>(i) - static_cast<Int>(xs.size()/2);
+		xs[i] = x;
+	    }
+	    for(SizeT i = 0; i != ts.size(); ++i){
+		const Int t = static_cast<Int>(i) - static_cast<Int>(ts.size()/2);
+		ts[i] = t;
+	    }
+	    mXRange = URangeFactory<Int>(xs).create();
+	    mTRange = URangeFactory<Int>(ts).create();
+	    Vector<RangePtr> rs { mTRange, mXRange, mXRange, mXRange };
+	    mGRange = YRangeFactory(rs).create();
+	    RangePtr g1 = CRangeFactory(1).create();
+	    RangePtr g2 = CRangeFactory(2).create();
+	    Vector<RangePtr> gs { g2, g1, g1, g2 };
+	    mGeom = YRangeFactory(gs).create();
+	    mRRange = rrange(mGRange, mGeom);
+	}
+
+	RangePtr mXRange;
+	RangePtr mTRange;
+	RangePtr mGRange;
+	RangePtr mGeom;
+	RangePtr mRRange;
+    };
+
+    TEST_F(RRange_Test, Basics)
+    {
+	EXPECT_EQ(mRRange->size(), mGRange->size());
+    }
+
+    TEST_F(RRange_Test, Local)
+    {
+	typedef UIndex<Int> UI;
+	MIndex<UI,UI,UI,UI> mi(mRRange->sub(1));
+	MIndex<CIndex,CIndex,CIndex,CIndex> ri(mRRange->sub(0));
+	EXPECT_EQ(mi.lmax().val(), mGRange->size()/mRRange->sub(0)->size());
+	EXPECT_EQ(ri.lmax().val(), getNumRanks());
+	const SizeT rank = getRankNumber();
+	ri = rank;
+	iter<0,4>( [&](auto i) {
+	    UI xi = *mi.pack()[i];
+	    const SizeT max = xi.lmax().val();
+	    const Int s = mGRange->sub(i)->size()/2;
+	    const SizeT rx = ri.pack()[i]->lex();
+	    for(; xi.lex() != max; ++xi){
+		EXPECT_EQ(*xi, static_cast<Int>(xi.lex()) - s + static_cast<Int>(max*rx));
+	    }
+	}, NoF {} );
+    }
+    
+    TEST_F(RRange_Test, Global)
+    {
+	RIndex<YIndex,YIndex> rgi(mRRange);
+	YIndex gi(mGRange);
+	for(SizeT c = 0; gi.lex() != gi.lmax().val(); ++gi, ++rgi, ++c) {
+	    const String s1 = gi.stringMeta();
+	    const String s2 = rgi.stringMeta();
+	    EXPECT_EQ(rgi.lex(), c);
+	    EXPECT_EQ(s1 ,s2);
+	}
+    }
+
+    TEST_F(RRange_Test, Global2)
+    {
+	typedef UIndex<Int> UI;
+	typedef MIndex<UI,UI,UI,UI> LocI;
+	typedef MIndex<CIndex,CIndex,CIndex,CIndex> RankI;
+	RIndex<LocI,RankI> rgi(mRRange);
+	LocI gi(mGRange);
+	for(SizeT c = 0; gi.lex() != gi.lmax().val(); ++gi, ++rgi, ++c) {
+	    const String s1 = gi.stringMeta();
+	    const String s2 = rgi.stringMeta();
+	    EXPECT_EQ(rgi.lex(), c);
+	    EXPECT_EQ(s1 ,s2);
+	}
+	rgi = 0;
+	auto j = rgi;
+	for(auto i = rgi; i.lex() != i.lmax().val(); ++i){
+	    j.at( i.meta() );
+	    if(getRankNumber() == 0){
+		EXPECT_EQ(i.lex(), j.lex());
+		EXPECT_EQ(i.rank(), j.rank());
+		EXPECT_EQ(i.local()->lex(), j.local()->lex());
+	    }
+	}
+    }
+    
+}
+
+int main(int argc, char** argv)
+{
+    ::testing::InitGoogleTest(&argc, argv);
+    ::testing::AddGlobalTestEnvironment( new MPIEnv(argc, argv) );
+    return RUN_ALL_TESTS();
+}

src/opt/mpi/tests/setbuf_unit_test.cc

@@ -0,0 +1,160 @@
+// -*- C++ -*-
+/**
+   
+   @file opt/mpi/tests/setbuf_unit_test.cc
+   @brief Setbuf unit tests.
+
+   Copyright (c) 2024 Christian Zimmermann. All rights reserved.
+   Mail: chizeta@f3l.de
+   
+**/
+
+#include <cstdlib>
+#include <iostream>
+
+#include "cnorxz.h"
+#include "cnorxz_mpi.h"
+#include "test_numbers.h"
+#include "rrange.cc.h"
+#include "mpi_env.h"
+
+namespace
+{
+    using namespace CNORXZ;
+    using Test::Numbers;
+    using namespace CNORXZ::mpi;
+
+    class Setbuf_Test : public ::testing::Test
+    {
+    protected:
+	
+	Setbuf_Test()
+	{
+	    CXZ_ASSERT(getNumRanks() == 4, "exptected 4 ranks");
+	    Vector<Int> xs(L);
+	    Vector<Int> ts(T);
+	    for(SizeT i = 0; i != xs.size(); ++i){
+		const Int x = static_cast<Int>(i) - static_cast<Int>(xs.size()/2);
+		xs[i] = x;
+	    }
+	    for(SizeT i = 0; i != ts.size(); ++i){
+		const Int t = static_cast<Int>(i) - static_cast<Int>(ts.size()/2);
+		ts[i] = t;
+	    }
+	    mSRange = CRangeFactory(4).create();
+	    mXRange = URangeFactory<Int>(xs).create();
+	    mTRange = URangeFactory<Int>(ts).create();
+	    Vector<RangePtr> rs { mTRange, mXRange, mXRange, mXRange };
+	    mGRange = YRangeFactory(rs).create();
+	    RangePtr g1 = CRangeFactory(1).create();
+	    RangePtr g2 = CRangeFactory(2).create();
+	    Vector<RangePtr> gs { g2, g1, g1, g2 };
+	    mGeom = YRangeFactory(gs).create();
+	    mRRange = rrange(mGRange, mGeom);
+	}
+
+	SizeT T = 16;
+	SizeT L = 12;
+	RangePtr mSRange;
+	RangePtr mXRange;
+	RangePtr mTRange;
+	RangePtr mGRange;
+	RangePtr mGeom;
+	RangePtr mRRange;
+    };
+
+
+
+    TEST_F(Setbuf_Test, run)
+    {
+	const SizeT myrank = getRankNumber();
+    
+	typedef UIndex<Int> UI;
+	typedef MIndex<UI,UI,UI,UI> LocI;
+	typedef MIndex<CIndex,CIndex,CIndex,CIndex> RankI;
+	auto rgi = std::make_shared<RIndex<LocI,RankI>>(mRRange);
+	auto rgj = std::make_shared<RIndex<LocI,RankI>>(mRRange);
+	auto rgk = std::make_shared<RIndex<LocI,RankI>>(mRRange);
+	LocI gi(mGRange);
+	LocI gj(mGRange);
+	auto ri = std::make_shared<RankI>(mGeom);
+	constexpr auto C0 = CSizeT<0> {};
+	constexpr auto C2 = CSizeT<2> {};
+	constexpr auto C3 = CSizeT<3> {};
+
+	const SizeT LSize = mRRange->sub(1)->size();
+	const SizeT blocks = mSRange->size();
+
+	MArray<Double> data(mRRange->sub(1)*mSRange);
+	for(SizeT i = 0; i != data.size(); ++i){
+	    data.data()[i] = static_cast<Double>(LSize*myrank*blocks+i);
+	}
+
+	*rgj = 0;
+	while(rgj->rank() != 1){
+	    ++*rgj;
+	}
+	*rgj->local() = 0;
+
+	Vector<Double> buf;
+	Vector<const Double*> map(mRRange->size(),nullptr);
+
+	auto shift = [&](const auto& x){
+	    auto o = x;
+	    std::get<0>(o) += 1;
+	    if(std::get<0>(o) >= static_cast<int>(T)/2) { std::get<0>(o) -= T; }
+	    std::get<2>(o) += 1;
+	    if(std::get<2>(o) >= static_cast<int>(L)/2) { std::get<2>(o) -= L; }
+	    std::get<3>(o) += 1;
+	    if(std::get<3>(o) >= static_cast<int>(L)/2) { std::get<3>(o) -= L; }
+	    return o;
+	};
+
+	Vector<bool> req(mRRange->size(), false);
+	for(auto ii = *rgi; ii.lex() != ii.lmax().val(); ++ii){
+	    auto jj = ii;
+	    jj.at( shift( ii.meta() ) );
+	    req[jj.pos()] = true;
+	}
+	setupBuffer(rgi, req, data, buf, map, mSRange->size());
+
+	EXPECT_EQ(mRRange->sub(1)->size(), 16*12*12*12/4);
+	const SizeT locsz = rgj->local()->lmax().val();
+	const SizeT myrankoff = myrank*locsz;
+	const SizeT mapsize = map.size();
+	// Fourth loop: Check:
+	for(*rgi = 0, gi = 0; rgi->lex() != rgi->lmax().val(); ++*rgi, ++gi){
+	    gj = gi.lex();
+	    *gj.pack()[C0] = (gj.pack()[C0]->lex() + 1) % gj.pack()[C0]->lmax().val();
+	    *gj.pack()[C2] = (gj.pack()[C2]->lex() + 1) % gj.pack()[C2]->lmax().val();
+	    *gj.pack()[C3] = (gj.pack()[C3]->lex() + 1) % gj.pack()[C3]->lmax().val();
+	    gj();
+	    *rgj = gj.lex();
+
+	    if(rgi->rank() == myrank){
+		const SizeT mpidx = (rgj->pos() - myrankoff + mapsize) % mapsize;
+		EXPECT_TRUE(map.data()[mpidx] != nullptr);
+		
+		const Double vn = *map[mpidx]/blocks;
+		const SizeT xp = static_cast<SizeT>(vn);
+		const SizeT orank = xp / mRRange->sub(1)->size();
+		if(myrank == 0){
+		    std::cout << " pos = " << rgj->pos() << " , val = " << *map[mpidx]
+			      << " , val_norm = " << vn << " , origin rank =  "
+			      << orank << std::endl;
+		}
+		EXPECT_EQ(orank, rgj->rank());
+		EXPECT_EQ(vn, rgj->pos());
+	    }
+	}
+	MPI_Barrier(MPI_COMM_WORLD);
+    }
+    
+}
+
+int main(int argc, char** argv)
+{
+    ::testing::InitGoogleTest(&argc, argv);
+    ::testing::AddGlobalTestEnvironment( new MPIEnv(argc, argv) );
+    return RUN_ALL_TESTS();
+}

src/tests/operation_unit_test.cc

@@ -155,6 +155,120 @@ namespace
 	COpRoot<double,MCCI1x> mOR1i2j;
     };
 
+    class Op_Test : public ::testing::Test
+    {
+    protected:
+
+	Op_Test()
+	{
+	    mXr = CRangeFactory(7).create();
+	    mAr = CRangeFactory(4).create();
+	    mIr = CRangeFactory(5).create();
+
+	    SizeT off = 15;
+
+	    const SizeT m1s = mAr->size() * mXr->size() * mAr->size();
+	    mM1 = Numbers::get(off, m1s);
+	    off += m1s;
+	    const SizeT m2s = mXr->size() * mXr->size() * mIr->size();
+	    mM2 = Numbers::get(off, m2s);
+	    off += m2s;
+	    const SizeT m3s = mXr->size() * mAr->size() * mAr->size() * mIr->size();
+	    mM3 = Numbers::get(off, m3s);
+	    off += m3s;
+	}
+
+	RangePtr mXr;
+	RangePtr mAr;
+	RangePtr mIr;
+	Vector<Double> mM1; // a,x,a
+	Vector<Double> mM2; // x,x,i
+	Vector<Double> mM3; // x,a,a,i
+    };
+
+    TEST_F(Op_Test, Assign)
+    {
+	Vector<Double> r1(mXr->size()*mXr->size()*mIr->size(),0);
+	Vector<Double> r2(mXr->size()*mIr->size(),0);
+
+	auto x = std::make_shared<CIndex>(mXr);
+	auto y = std::make_shared<CIndex>(mXr);
+	auto i = std::make_shared<SIndex<SizeT,5>>(mIr);
+	
+	oproot(r1.data(), mindexPtr(y*x*i)) = coproot(mM2.data(), mindexPtr(y*x*i));
+	oproot(r2.data(), mindexPtr(x*i)) = coproot(mM2.data(), mindexPtr(x*x*i));
+
+	for(SizeT ii = 0; ii != r1.size(); ++ii){
+	    EXPECT_EQ(r1[ii], mM2[ii]);
+	}
+
+	for(SizeT xi = 0; xi != mXr->size(); ++xi){
+	    for(SizeT ii = 0; ii != mIr->size(); ++ii){
+		const SizeT r2i = xi*mIr->size() + ii;
+		const SizeT m2i = xi*mIr->size()*(1 + mXr->size()) + ii;
+		EXPECT_EQ(r2[r2i], mM2[m2i]);
+	    }
+	}
+    }
+
+    TEST_F(Op_Test, Multiply)
+    {
+	Vector<Double> r1(mIr->size()*mXr->size()*mAr->size(),0);
+
+	auto x = std::make_shared<CIndex>(mXr);
+	auto a = std::make_shared<SIndex<SizeT,4>>(mAr);
+	auto i = std::make_shared<SIndex<SizeT,5>>(mIr);
+
+	oproot(r1.data(), mindexPtr(i*x*a)) = coproot(mM2.data(), mindexPtr(x*x*i))
+	    * coproot(mM3.data(), mindexPtr(x*a*a*i)) * coproot(mM1.data(), mindexPtr(a*x*a));
+
+	for(SizeT ii = 0; ii != mIr->size(); ++ii){
+	    for(SizeT xi = 0; xi != mXr->size(); ++xi){
+		for(SizeT ai = 0; ai != mAr->size(); ++ai){
+		    const SizeT r1i = (ii*mXr->size() + xi)*mAr->size() + ai;
+		    const SizeT m2i = (xi*mXr->size() + xi)*mIr->size() + ii;
+		    const SizeT m3i = ((xi*mAr->size() + ai)*mAr->size() + ai)*mIr->size() + ii;
+		    const SizeT m1i = (ai*mXr->size() + xi)*mAr->size() + ai;
+		    EXPECT_EQ(r1[r1i], mM2[m2i]*mM3[m3i]*mM1[m1i]);
+		}
+	    }
+	}
+    }
+
+    TEST_F(Op_Test, Contract)
+    {
+	Vector<Double> r1(mIr->size()*mXr->size(),0);
+
+	auto x = std::make_shared<CIndex>(mXr);
+	auto y = std::make_shared<CIndex>(mXr);
+	auto a = std::make_shared<SIndex<SizeT,4>>(mAr);
+	auto b = std::make_shared<SIndex<SizeT,4>>(mAr);
+	auto i = std::make_shared<SIndex<SizeT,5>>(mIr);
+
+	oproot(r1.data(), mindexPtr(i*x)) =
+	    (coproot(mM2.data(), mindexPtr(y*x*i)) * coproot(mM3.data(), mindexPtr(x*a*b*i)) *
+	     coproot(mM1.data(), mindexPtr(b*y*a))).c(mindexPtr(a*b*y));
+
+	for(SizeT ii = 0; ii != mIr->size(); ++ii){
+	    for(SizeT xi = 0; xi != mXr->size(); ++xi){
+		Double rx = 0;
+		const SizeT r1i = ii*mXr->size() + xi;
+		for(SizeT ai = 0; ai != mAr->size(); ++ai){
+		    for(SizeT bi = 0; bi != mAr->size(); ++bi){
+			for(SizeT yi = 0; yi != mXr->size(); ++yi){
+			    const SizeT m2i = (yi*mXr->size() + xi)*mIr->size() + ii;
+			    const SizeT m3i = ((xi*mAr->size() + ai)*mAr->size() + bi) *
+				mIr->size() + ii;
+			    const SizeT m1i = (bi*mXr->size() + yi)*mAr->size() + ai;
+			    rx += mM2.at(m2i) * mM3.at(m3i) * mM1.at(m1i);
+			}
+		    }
+		}
+		EXPECT_FLOAT_EQ(r1[r1i], rx);
+	    }
+	}
+    }
+    
     TEST_F(OpCont_CR_Test, Basics)
     {
 	EXPECT_EQ(mOp2.data(), mData1.data());

src/tests/range_unit_test.cc

@@ -51,6 +51,20 @@ namespace
 	RangePtr ur;
     };
 
+    class WR_Test : public ::testing::Test
+    {
+    protected:
+
+	WR_Test()
+	{
+	    mMeta = { "These", "are", "test", "strings", "foo", "bar", "baz" };
+	    wr = WRangeFactory<String>(mMeta).create();
+	}
+
+	Vector<String> mMeta;
+	RangePtr wr;
+    };
+
     class MR_Test : public ::testing::Test
     {
     protected:
@@ -211,6 +225,38 @@ namespace
 	EXPECT_EQ(prx->orig(),ur);
     }
 
+    TEST_F(WR_Test, Basics)
+    {
+	auto wrx = std::dynamic_pointer_cast<WRange<String>>(wr);
+	EXPECT_FALSE(wrx == nullptr);
+	EXPECT_EQ(wr->size(), mMeta.size());
+	EXPECT_EQ(wrx->size(), mMeta.size());
+	
+	EXPECT_TRUE(wrx->begin() != wrx->end());
+	EXPECT_FALSE(wrx->begin() == wrx->end());
+	EXPECT_EQ(wrx->begin().pos(), 0u);
+	EXPECT_EQ(wrx->end().pos(), mMeta.size());
+
+	EXPECT_TRUE(wr->begin() != wr->end());
+	EXPECT_FALSE(wr->begin() == wr->end());
+	EXPECT_EQ(wr->begin().pos(), 0u);
+	EXPECT_EQ(wr->end().pos(), mMeta.size());
+
+	SizeT cnt = 0;
+	auto endxi = wr->end();
+	for(auto xi = wr->begin(); xi != endxi; ++xi){
+	    EXPECT_EQ(xi.pos(), cnt);
+	    EXPECT_TRUE(*xi == DType(mMeta[cnt]));
+	    EXPECT_EQ((*xi).str(), mMeta[cnt]);
+	    ++cnt;
+	}
+
+	SizeT cnt2 = 0;
+	for(const auto& x: *wrx){
+	    EXPECT_EQ(x, mMeta[cnt2++]);
+	}
+    }
+
     TEST_F(MR_Test, Basics2d)
     {
 	auto mrx = std::dynamic_pointer_cast<MRange<CRange,URange<String>>>(mr);

src/tests/xpr_unit_test.cc

@@ -17,12 +17,6 @@ namespace
 {
     using namespace CNORXZ;
 
-    template <class PosT1, class PosT2>
-    constexpr auto mkMPos(const PosT1& a, const PosT2& b)
-    {
-	return MPos<PosT1,PosT2>(a,b);
-    }
-
     class Pos_Test : public ::testing::Test
     {
     protected:
@@ -253,7 +247,7 @@ namespace
     TEST_F(For_Test, SFor)
     {
 	auto loop = mkSFor<sSize>(IndexId<0>(mId1), TestXpr1( IndexId<0>(mId1) ),
-				  [](const auto& a, const auto& b) { return a + b; });
+				  [](auto& a, const auto& b) { return a += b; });
 
 	const UPos rs = loop.rootSteps(IndexId<0>(mId1));
 	EXPECT_EQ(rs.val(), 1u);