WIP: MPI: RAIndex + RCArray + corresp tests

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

@@ -0,0 +1,119 @@
+// -*- 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__
+
+namespace CNOXRZ
+{
+    namespace mpi
+    {
+	template <typename T>
+	RAIndex<T>::RAIndex(const T* loc, const RangePtr& range, SizeT lexpos) :
+	    RIndex<YIndex,YIndex>(range, lexpos),
+	    mLoc(loc),
+	    mCur(i.rank()),
+	    mThisRank(getRankNumber())
+	{
+	    setBufferSize();
+	    setBuffer();
+	}
+	    
+
+	template <typename T>
+	RAIndex<T>::RAIndex(const T* loc, const RIndex<YIndex,YIndex>& i) :
+	    RIndex<YIndex,YIndex>(i)
+	    mLoc(loc),
+	    mCur(i.rank()),
+	    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(rank() != mThisRank){
+		if(mCur != rank()){
+		    setBuffer();
+		}
+		return mBuf[local()->pos() % mBufSize];
+	    }
+	    else {
+		mLoc[local()->pos()];
+	    }
+	}
+
+	template <typename T>
+	const T* RAIndex<T>::operator->() const
+	{
+	    if(rank() != mThisRank){
+		if(mCur != rank()){
+		    setBuffer();
+		}
+		return mBuf + 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()
+	{
+	    if(mBuf.size() != mBufSize){
+		mBuf.resize(mBufSize);
+	    }
+	    // A Bcast alternative with const pointer to source would be better...
+	    std::memcpy(mBuf.data(), mLoc + local()->pos() / mBufSize, mBufSize*sizeof(T));
+	    MPI_Bcast(mBuf.data(), mBufSize*sizeof(T), MPI_BYTE, static_cast<int>(rank()),
+		      MPI_COMM_WORLD );
+	}
+	
+    } // namespace mpi
+} // namespace CNOXRZ

src/opt/mpi/include/raindex.h

@@ -16,35 +16,49 @@
 
 namespace CNORXZ
 {
-
+    namespace mpi
-    /** ****
-	Index for multi-ranked array.
-	This is meant to be a global index. For local iterations
-	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);
+	    Index for multi-ranked array.
-	RAIndex(const T* loc, const RIndex<YIndex,YIndex>& i);
+	    This is meant to be a global index, its position must be equal on each rank,
-	RAIndex(const T* loc, const RAIndex<T>& i);
+	    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;
 
-	RAIndex operator+(Int n) const;
+	    DEFAULT_MEMBERS(RAIndex);
-	RAIndex operator-(Int n) const;
+	    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);
 
-	const T& operator*() const;
+	    RAIndex operator+(Int n) const;
-	const T& operator->() const;
+	    RAIndex operator-(Int n) const;
 
-    private:
+	    const T& operator*() const;
-	const T* mLoc = nullptr;
+	    const T* operator->() const;
-	Vector<T> mBuf; // used if iterating over content on different rank 
-    };
 
-}
+	    void setBufferSize();
+	    
+	private:
+
+	    void setBuffer();
+	    
+	    const T* mLoc = nullptr;
+	    Vector<T> mBuf; // used if iterating over content on different rank
+	    SizeT mCur; // current rank in the buffer
+	    SizeT mBufSize;
+	    SizeT mThisRank;
+	};
+
+    } // namespace mpi
+} // namespace CNOXRZ
 
 #endif

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

@@ -20,7 +20,8 @@ namespace CNORXZ
 	template <typename T>
 	RCArray<T>::RCArray(const Sptr<CArrayBase<T>> a, const RangePtr& geom) :
 	    mA(a),
-	    mGeom(geom)
+	    mGeom(geom),
+	    mGlobal(RRangeFactory(a->range(),mGeom).create())
 	{}
 
 	template <typename T>
@@ -53,22 +54,30 @@ namespace CNORXZ
 	template <class... Indices>
 	T RCArray<T>::operator[](const SPack<Indices...>& pack) const
 	{
+	    CXZ_ERROR("not implemented");
+	    return T();
 	}
 
 	template <typename T>
 	template <class... Indices>
 	T RCArray<T>::at(const SPack<Indices...>& pack) const
 	{
+	    CXZ_ERROR("not implemented");
+	    return T();
 	}
 
 	template <typename T>
 	T RCArray<T>::operator[](const DPack& pack) const
 	{
+	    CXZ_ERROR("not implemented");
+	    return T();
 	}
 
 	template <typename T>
 	T RCArray<T>::at(const DPack& pack) const
 	{
+	    CXZ_ERROR("not implemented");
+	    return T();
 	}
 
 	template <typename T>
@@ -76,23 +85,31 @@ namespace CNORXZ
 	Sptr<CArrayBase<T>> RCArray<T>::sl(const IndexInterface<I,M>& begin,
 					   const IndexInterface<I,M>& end) const
 	{
+	    CXZ_ERROR("not implemented");
+	    return nullptr;
 	}
 
 	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
 	{
+	    CXZ_ERROR("not implemented");
+	    return COpRoot<T,Index>();
 	}
 
 	template <typename T>
 	inline decltype(auto) RCArray<T>::operator()(const DPack& pack) const
 	{
+	    CXZ_ERROR("not implemented");
+	    return COpRoot<T,Index>();
 	}
 
 	template <typename T>
@@ -110,27 +127,31 @@ namespace CNORXZ
 	template <typename T>
 	RangePtr RCArray<T>::range() const
 	{
-	    return RRangeFactory(mA->range(),mGeom).create();
+	    return mGlobal;
 	}
 
 	template <typename T>
 	const_iterator RCArray<T>::begin() const
 	{
+	    return const_iterator(mA.data(), mGlobal);
 	}
 
 	template <typename T>
 	const_iterator RCArray<T>::end() const
 	{
+	    return const_iterator(mA.data(), mGlobal, mGlobal->size());
 	}
 
 	template <typename T>
 	const_iterator RCArray<T>::cbegin() const
 	{
+	    return const_iterator(mA.data(), mGlobal);
 	}
 
 	template <typename T>
 	const_iterator RCArray<T>::cend() const
 	{
+	    return const_iterator(mA.data(), mGlobal, mGlobal->size());
 	}
 
 	template <typename T>
@@ -145,6 +166,11 @@ namespace CNORXZ
 	    return *mA;
 	}
 
+	template <typename T>
+	RangePtr RCArray<T>::geom() const
+	{
+	    return mGeom;
+	}
 	
     } // namespace mpi
 } // namespace CNORXZ

src/opt/mpi/include/rarray.h

@@ -103,9 +103,13 @@ namespace CNORXZ
 	    /** Get local array object. */
 	    const CArrayBase<T>& local() const;
 
+	    /** Get rank geometry. */
+	    RangePtr geom() const;
+	    
 	private:
 	    ObjHandle<CArrayBase<T>> mA;
 	    RangePtr mGeom;
+	    RangePtr mGlobal;
 	};
 
 	/** ****

src/opt/mpi/tests/rarray_unit_test.cc

@@ -0,0 +1,123 @@
+// -*- 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"
+
+namespace
+{
+    using namespace CNORXZ;
+    using Test::Numbers;
+    using namespace CNORXZ::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;
+    };
+
+
+    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();
+	    Vector<Double> vec = Numbers::get(0,size/4+10);
+	    Vector<Double> data(size);
+	    Vector<Double> mData(size*4);
+	    const SizeT myrank = getRankNumber();
+	    for(SizeT i = 0; i != size; ++i){
+		data[i] = vec[i] * vec[i+myrank] / vec[i+2*myrank];
+		mData[i + size*myrank] = data[i];
+		MPI_Bcast(mData.data() + size*i, size, MPI_DOUBLE, i, MPI_COMM_WORLD);
+	    }
+	    mLoc = std::make_shared<MArray<Double>>( mRRange->sub(1), data);
+	    mA = RCArray<Double>(mLoc, mGeom);
+	}
+
+	RangePtr mXRange;
+	RangePtr mTRange;
+	RangePtr mGRange;
+	RangePtr mGeom;
+	RangePtr mRRange;
+	Sptr<MArray<Double>> mLoc;
+	RCArray mA;
+	Vector<Double> mData;
+    };
+
+    TEST_F(RArray_Test, Basics)
+    {
+	EXPECT_EQ(mA.size(), mRRange->size());
+    }
+    
+    TEST_F(RArray_Test, GlobalIterate)
+    {
+	const SizeT size = mRRange->sub(1)->size();
+	auto e = mA.end();
+	for(auto i = mA.begin(); i != e; ++i){
+	    EXPECT_EQ(*i, 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();
+}