more operation tests + fix SFor (compatibility with For) + more RIndex utils

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     |
      +=================*/

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...>>
     {

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/opt/mpi/include/rarray.cc.h

@@ -152,7 +152,7 @@ namespace CNORXZ
 	    typedef typename std::remove_reference<decltype(*pack[CSizeT<0>{}])>::type I0;
 	    if constexpr(is_rank_index<I0>::value){
 		// preliminary:
-		CXZ_ASSERT(this->formatIsTrivial(),
+		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)>
@@ -235,15 +235,15 @@ namespace CNORXZ
 	template <typename T>
 	const Vector<const T*>& RCArray<T>::buffermap() const
 	{
-	    return mBuf;
+	    return mMap;
 	}
 	
 	template <typename T>
-	template <class TarI, class RTarI, class SrcI, class RSrcI, class I>
+	template <class LoopI, class SrcI, class RSrcI, class I>
-	void RCArray<T>::load(const Sptr<RIndex<TarI,RTarI>>& lpi, const Sptr<RIndex<SrcI,RSrcI>>& ai,
+	void RCArray<T>::load(const Sptr<LoopI>& lpi, const Sptr<RIndex<SrcI,RSrcI>>& ai,
 			      const Sptr<I>& i, const Sptr<Vector<SizeT>>& imap) const
 	{
-	    mA->checkFormatCompatibility(ai->local()*i);
+	    mA->checkFormatCompatibility(mindex(ai->local()*i));
 	    const SizeT blocks = i->pmax().val();
 	    setupBuffer(lpi, ai, imap, *mA, mBuf, mMap, blocks);
 	}

src/opt/mpi/include/rarray.h

@@ -130,8 +130,8 @@ namespace CNORXZ
 		@param i Index for non-ranked dimensions.
 		@param imap Index position map.
 	     */
-	    template <class TarI, class RTarI, class SrcI, class RSrcI, class I>
+	    template <class LoopI, class SrcI, class RSrcI, class I>
-	    void load(const Sptr<RIndex<TarI,RTarI>>& lpi, const Sptr<RIndex<SrcI,RSrcI>>& ai,
+	    void load(const Sptr<LoopI>& lpi, const Sptr<RIndex<SrcI,RSrcI>>& ai,
 		      const Sptr<I>& i, const Sptr<Vector<SizeT>>& imap) const;
 
 	    /** Load all data from other ranks that is accessed by f(i).

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

@@ -11,14 +11,14 @@ namespace CNORXZ
 	 +================*/
     
 	template <typename T, class RIndexT, class IndexT>
-	constexpr CROpRoot<T,RIndexT,IndexT>::CROpRoot(const RArray<T>& a, const Sptr<RIndexT>& ri,
+	constexpr CROpRoot<T,RIndexT,IndexT>::CROpRoot(const RCArray<T>& a, const Sptr<RIndexT>& ri,
 						       const Sptr<IndexT>& li) :
-	    mData(a.map()),
+	    mData(a.buffermap().data()),
 	    mRIndex(ri),
 	    mIndex(li)
 	{
-	    CXZ_ASSERT(a.map().size() == ri->lmax().val(),
+	    CXZ_ASSERT(a.buffermap().size() == ri->lmax().val(),
-		       "data map not properly initialized: map size = " << a.map.size()
+		       "data map not properly initialized: map size = " << a.buffermap().size()
 		       << ", rank index range size = " << ri->lmax().val());
 	}
 
@@ -43,7 +43,7 @@ namespace CNORXZ
 	}
 
 	template <typename T, class RIndexT, class IndexT>
-	constexpr decltype(auto) croproot(const RArray<T>& a, const Sptr<RIndexT>& ri,
+	constexpr decltype(auto) croproot(const RCArray<T>& a, const Sptr<RIndexT>& ri,
 					  const Sptr<IndexT>& li)
 	{
 	    return CROpRoot<T,RIndexT,IndexT>(a,ri,li);
@@ -56,12 +56,12 @@ namespace CNORXZ
 	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.map()),
+	    mData(a.buffermap().data()),
 	    mRIndex(ri),
 	    mIndex(li)
 	{
-	    CXZ_ASSERT(a.map().size() == ri->lmax().val(),
+	    CXZ_ASSERT(a.buffermap().size() == ri->lmax().val(),
-		       "data map not properly initialized: map size = " << a.map.size()
+		       "data map not properly initialized: map size = " << a.buffermap().size()
 		       << ", rank index range size = " << ri->lmax().val());
 	}
 

src/opt/mpi/include/rop_types.h

@@ -25,7 +25,7 @@ namespace CNORXZ
 	    typedef COpInterface<CROpRoot<T,RIndexT,IndexT>> OI;
 
 	    constexpr CROpRoot() = default;
-	    constexpr CROpRoot(const RArray<T>& a, const Sptr<RIndexT>& ri,
+	    constexpr CROpRoot(const RCArray<T>& a, const Sptr<RIndexT>& ri,
 			       const Sptr<IndexT>& li);
 	    
 	    template <class PosT>
@@ -38,13 +38,13 @@ namespace CNORXZ
 
 	private:
 
-	    const T** mData;
+	    const T* const* mData;
 	    Sptr<RIndexT> mRIndex;
 	    Sptr<IndexT> mIndex;
 	};
 
 	template <typename T, class RIndexT, class IndexT>
-	constexpr decltype(auto) croproot(const RArray<T>& a, const Sptr<RIndexT>& ri,
+	constexpr decltype(auto) croproot(const RCArray<T>& a, const Sptr<RIndexT>& ri,
 					  const Sptr<IndexT>& li);
 	
 	template <typename T, class RIndexT, class IndexT>

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

@@ -382,6 +382,12 @@ namespace CNORXZ
 	    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    |
 	 +=====================*/

src/opt/mpi/include/rrange.h

@@ -179,6 +179,13 @@ namespace CNORXZ
 	    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

src/opt/mpi/tests/roperation_unit_test.cc

@@ -108,19 +108,25 @@ namespace
 	RArray<Double> res( MArray<Double>(mRXRange->sub(1)), mGeom2 );
 	EXPECT_EQ(res.size(), comp.size());
 	typedef UIndex<SizeT> UI;
-	RIndex<MIndex<UI,UI,UI,UI>,MIndex<UI,UI,UI,UI>> x(mRXRange);
+	auto x = std::make_shared<RIndex<MIndex<UI,UI,UI,UI>,MIndex<UI,UI,UI,UI>>>(mRXRange);
-	RIndex<MIndex<UI,UI,UI,UI>,MIndex<UI,UI,UI,UI>> y(mRXRange);
+	auto y = std::make_shared<RIndex<MIndex<UI,UI,UI,UI>,MIndex<UI,UI,UI,UI>>>(mRXRange);
-	SIndex<SizeT,4> A(mSpRange);
+	auto xy = std::make_shared<RIndex<MIndex<UI,UI,UI,UI>,MIndex<UI,UI,UI,UI>>>(mRXRange);
-	SIndex<SizeT,4> B(mSpRange);
+	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;
+	
 	// block 1:
-	mM1.load(,A*B*a*b);
+	//mM1.load(mindexPtr(y*x), x, mindexPtr(A*B), imap1);
-	mM2.load(,A*B*a*b);
+	//mM2.load(mindexPtr(y*x), xy, mindexPtr(A*B), imap2);
-	res(y) += (mM1(x*A*B*a*b) * mM2(xy*B*A*b*a)).c(x*A*B*a*b);
+	/*
+	res(y) = ( mM1(x*A*B) * mM2(xy*B*A)).c(mindexPtr(x*A*B));
 	// 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(x*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);

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/xpr_unit_test.cc

@@ -247,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);