diff --git a/src/include/ranges/mrange.cc.h b/src/include/ranges/mrange.cc.h
index 6d03ff2..34785c1 100644
--- a/src/include/ranges/mrange.cc.h
+++ b/src/include/ranges/mrange.cc.h
@@ -631,6 +631,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)
     {
diff --git a/src/include/ranges/mrange.h b/src/include/ranges/mrange.h
index ec8f37d..6d59cb9 100644
--- a/src/include/ranges/mrange.h
+++ b/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.
      */
diff --git a/src/include/xpr/pos_type.cc.h b/src/include/xpr/pos_type.cc.h
index 7fead30..98c6409 100644
--- a/src/include/xpr/pos_type.cc.h
+++ b/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
     {
@@ -214,6 +215,7 @@ namespace CNORXZ
     template <class PosT1>
     constexpr UPos FPos::operator()(const PosT1& a) const
     {
+	assert(0);
 	return UPos(mExt * mMap[a.val()]);
     }
 
@@ -234,6 +236,11 @@ namespace CNORXZ
 	return val();
     }
 
+    constexpr SizeT FPos::max() const
+    {
+	return mMax;
+    }
+    
     /*===========+
      |   SFPos   |
      +===========*/
@@ -320,6 +327,12 @@ namespace CNORXZ
 	return val();
     }
 
+    template <SizeT N, SizeT... Ms>
+    constexpr SizeT SFPos<N,Ms...>::max() const
+    {
+	return sizeof...(Ms); // CHECK!!!
+    }
+    
     /*==========+
      |   MPos   |
      +==========*/
diff --git a/src/include/xpr/pos_type.h b/src/include/xpr/pos_type.h
index 09f300e..6d524d9 100644
--- a/src/include/xpr/pos_type.h
+++ b/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,6 +120,10 @@ 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 <SizeT N, SizeT... Ms>
@@ -151,6 +157,8 @@ namespace CNORXZ
 	explicit constexpr operator FPos() const;
 
 	explicit constexpr operator SizeT() const;
+
+	constexpr SizeT max() const;
     };
 
     template <class BPosT, class NPosT>
diff --git a/src/opt/mpi/include/rrange.cc.h b/src/opt/mpi/include/rrange.cc.h
index 986cea1..e0111ed 100644
--- a/src/opt/mpi/include/rrange.cc.h
+++ b/src/opt/mpi/include/rrange.cc.h
@@ -53,6 +53,15 @@ namespace CNORXZ
 	    *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)
 	{
@@ -196,7 +205,7 @@ namespace CNORXZ
 	template <SizeT I>
 	decltype(auto) RIndex<IndexI,IndexK>::stepSize(const IndexId<I>& id) const
 	{
-	    return mK->stepSize(id) * mI->lmax() + mI->stepSize(id);
+	    return getRankStepSize(id) * mI->pmax() + mI->stepSize(id);
 	}
 
 	template <class IndexI, class IndexK>
@@ -283,8 +292,7 @@ namespace CNORXZ
 	template <class Xpr, class F>
 	constexpr decltype(auto) RIndex<IndexI,IndexK>::ifor(const Xpr& xpr, F&& f) const
 	{
-	    CXZ_ERROR("not implemented");
-	    return xpr;
+	    return mI->ifor(xpr, std::forward<F>(f));
 	}
 
 	template <class IndexI, class IndexK>
@@ -299,6 +307,13 @@ namespace CNORXZ
 	    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()()
 	{
diff --git a/src/opt/mpi/include/rrange.h b/src/opt/mpi/include/rrange.h
index 1179642..7c2c92f 100644
--- a/src/opt/mpi/include/rrange.h
+++ b/src/opt/mpi/include/rrange.h
@@ -56,6 +56,12 @@ namespace CNORXZ
 	     */
 	    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);
 
@@ -136,6 +142,11 @@ namespace CNORXZ
 	    /** @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()();
 	    
@@ -150,8 +161,50 @@ namespace CNORXZ
 	    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. */
+
+	    mutable std::map<PtrId,Vector<SizeT>> mRankMap;
+
+	    template <SizeT I>
+	    decltype(auto) getRankStepSize(const IndexId<I>& id) const
+	    {
+		auto ss = mI->stepSize(id);
+		FPos x;
+		VCHECK(typeid(ss).name());
+		VCHECK(typeid(x).name());
+		assert(0);
+		if constexpr(std::is_same<decltype(ss),FPos>::value){
+		    assert(0);
+		    if(mRankMap.count(id.id()) != 0){
+			return FPos(ss.val(), mRankMap[id.id()].data());
+		    }
+		    else {
+			Vector<SizeT> mp(ss.max());
+			for(SizeT i = 0; i != mp.size(); ++i){
+			    //max2 = num ranks in this dir (preliminary solution)!!!
+			    mp[i] = ( ss.map()[i] / ss.max() ) % ss.max2(); 
+			}
+			mRankMap[id.id()] = mp;
+			return FPos(ss.val(), mp.data());
+		    }
+		}
+		else {
+		    return SPos<0> {};
+		}
+	    }
 	};
 
+	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);
+	}
+
 	// Traits!!!
 	template <class I>
 	struct is_rank_index
diff --git a/src/opt/mpi/tests/rindex_exp.cc b/src/opt/mpi/tests/rindex_exp.cc
index 736baaf..93c4c27 100644
--- a/src/opt/mpi/tests/rindex_exp.cc
+++ b/src/opt/mpi/tests/rindex_exp.cc
@@ -52,14 +52,163 @@ namespace
 	RangePtr mGeom;
 	RangePtr mRRange;
     };
+
+    template <class Index, class IndexR>
+    class IndexMap : public Index
+    {
+    public:
+	typedef typename Index::IB IB;
+	typedef typename Index::RangeType RangeType;
+
+	DEFAULT_MEMBERS(IndexMap);
+
+	IndexMap(const Sptr<Index>& i, const Vector<SizeT>& map, const Sptr<IndexR>& r) :
+	    Index(*i), mI(i), mR(r), mMap(map) {}
+
+	constexpr decltype(auto) id() const
+	{
+	    return mI->id();
+	}
+
+	template <SizeT I>
+	decltype(auto) stepSize(const IndexId<I>& id) const
+	{
+	    // TODO: new Pos Type!!!
+	    //assert(0);
+	    return FPos( mI->stepSize(id).val(), mMap.data(), mI->lmax().val(), mR->lmax().val() );
+	}
+
+    private:
+	Sptr<Index> mI;
+	Sptr<IndexR> mR;
+	Vector<SizeT> mMap;
+    };
+
+    template <class Index, class IndexR>
+    constexpr decltype(auto) shift(const Sptr<Index>& i, const Sptr<IndexR>& r)
+    {
+	// only one-dim indices!!!
+	auto j = *i;
+	Vector<SizeT> mp(j.lmax().val());
+	const SizeT L = i->lmax().val()*r->lmax().val();
+	for(j = 0; j.lex() != j.lmax().val(); ++j){
+	    mp[j.lex()] = ( j.lex() + 1 + L) % L;
+	}
+	return std::make_shared<IndexMap<Index,IndexR>>(i, mp, r);
+    }
+    
+    template <class Index>
+    class PosOp : public COpInterface<PosOp<Index>>
+    {
+    public:
+	typedef COpInterface<PosOp<Index>> OI;
+
+	constexpr PosOp() = default;
+
+	constexpr PosOp(const Sptr<Index>& i, SizeT block) :
+	    mMyrank(getRankNumber()), mI(i), mBlock(block) {}
+
+	template <class PosT>
+	constexpr decltype(auto) operator()(const PosT& pos) const
+	{
+	    return static_cast<SizeT>(pos)+mMyrank*mBlock;
+	}
+
+	constexpr decltype(auto) operator()() const
+	{
+	    return static_cast<SizeT>(mMyrank*mBlock);
+	}
+
+	template <SizeT I>
+	constexpr decltype(auto) rootSteps(const IndexId<I>& id) const
+	{
+	    return mI->stepSize(id);
+	}
+
+    private:
+	SizeT mMyrank;
+	Sptr<Index> mI;
+	SizeT mBlock;
+    };
+
+    template <class Index>
+    constexpr decltype(auto) posop(const Sptr<Index>& i, SizeT block)
+    {
+	return PosOp<Index>(i, block);
+    }
 }
 
-
-int main(int argc, char** argv)
+void run2(const Env& env)
 {
-    MPI_Init(&argc, &argv);
+    const SizeT myrank = getRankNumber();
+    const SizeT Nranks = getNumRanks();
+    
+    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>>(env.mRRange);
+    auto rgj = std::make_shared<RIndex<LocI,RankI>>(env.mRRange);
+    auto rgk = std::make_shared<RIndex<LocI,RankI>>(env.mRRange);
+    LocI gi(env.mGRange);
+    LocI gj(env.mGRange);
+    auto ri = std::make_shared<RankI>(env.mGeom);
+    constexpr auto C0 = CSizeT<0> {};
+    constexpr auto C1 = CSizeT<1> {};
+    constexpr auto C2 = CSizeT<2> {};
+    constexpr auto C3 = CSizeT<3> {};
 
-    Env env;
+    const SizeT LSize = env.mRRange->sub(1)->size();
+    
+    const SizeT blocks = env.mSRange->size();
+    Vector<Double> data(LSize*blocks);
+    Vector<Double> buf;
+    Vector<Double*> map(env.mRRange->size(),nullptr);
+    for(SizeT i = 0; i != data.size(); ++i){
+	data[i] = static_cast<Double>(LSize*myrank*blocks+i);
+    }
+    Vector<Vector<SizeT>> cnt(Nranks);
+    for(auto& c: cnt){
+	c.resize(Nranks);
+    }
+    Vector<Vector<Double>> sendbuf(Nranks);
+    for(auto& sb: sendbuf){
+	sb.reserve(data.size());
+    }
+
+    auto srgi = rindexPtr( mindexPtr(
+	shift(rgi->local()->pack()[C0], ri->pack()[C0]),
+	rgi->local()->pack()[C1],
+	shift(rgi->local()->pack()[C2], ri->pack()[C2]),
+	shift(rgi->local()->pack()[C3], ri->pack()[C3])
+    ), ri );
+    VCHECK(srgi->lmax().val());
+    
+    *rgj = 0;
+    while(rgj->rank() != 1){
+	++*rgj;
+    }
+    *rgj->local() = 0;
+    rgi->ifor( operation
+	       ( [&](SizeT p) {
+		   gj = rgj->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();
+		   *rgk = gj.lex();
+		   if(myrank == 1){
+		       std::cout << p << " " << rgk->pos() << " " << rgj->pos() << std::endl;
+		   }
+		   ++*rgj->local();
+		   (*rgj)();
+	       } , posop(srgi, rgi->local()->pmax().val()) ) ,
+	       NoF {} )();
+    
+    MPI_Barrier(MPI_COMM_WORLD);
+}
+
+void run1(const Env& env)
+{
     const SizeT myrank = getRankNumber();
     const SizeT Nranks = getNumRanks();
     
@@ -92,6 +241,7 @@ int main(int argc, char** argv)
     for(auto& sb: sendbuf){
 	sb.reserve(data.size());
     }
+
     // First loop: setup send buffer
     for(rgi = 0, gi = 0; rgi.lex() != rgi.lmax().val(); ++rgi, ++gi){
 	gj = gi.lex();
@@ -183,6 +333,16 @@ int main(int argc, char** argv)
 	    assert(vn == rgj.pos());
 	}
     }
+    MPI_Barrier(MPI_COMM_WORLD);
+}
+
+int main(int argc, char** argv)
+{
+    MPI_Init(&argc, &argv);
+
+    Env env;
+    run1(env);
+    run2(env);
     
     MPI_Finalize();
     return 0;