Update Spherical kernel interface to use the new data organisation

 - Divide the FSphericalCell inner data layout into two sub-types:
   multipole_t and local_expansion_t. As the multipole and the local
   expansion layout is the same, they share a common implementation.
   Two class attribute are accessible through the `getMultipoledata` and
   `getLocalExpansionData` methods.

 - Change FMM operators signature to take advantage of the new
   layout. The operators only take as parameters the data they may need
   instead of the whole cell. For instance (simplified):

   void M2M(CellClass* parent, CellClass** children);

   becomes

   void M2M(multipole_t*  parent_m,    symbolic_data_t* parent_s,
            multipole_t** children_ms, symbolic_data_t* children_ss);

Src/Kernels/Spherical/FAbstractSphericalKernel.hpp

@@ -144,10 +144,14 @@ public:
     }
 
     /** P2M with a cell and all its particles */
-    void P2M(CellClass* const inPole, const ContainerClass* const inParticles) override {
-        FComplex<FReal>* FRestrict const cellMultiPole = inPole->getMultipole();
+    template<class SymbolicData>
+    void P2M(typename CellClass::multipole_t* const LeafMultipole,
+             const SymbolicData* const LeafSymb,
+             const ContainerClass* const inParticles)
+    {
+        FComplex<FReal>* FRestrict const cellMultiPole = LeafMultipole->get();
         // Copying the position is faster than using cell position
-        const FPoint<FReal> polePosition = getLeafCenter(inPole->getCoordinate());
+        const FPoint<FReal> polePosition = getLeafCenter(LeafSymb->getCoordinate());
         // For all particles in the leaf box
         const FReal*const physicalValues = inParticles->getPhysicalValues();
         const FReal*const positionsX = inParticles->getPositions()[0];
@@ -162,38 +166,58 @@ public:
     }
 
     /** M2M with a cell and all its child */
-    void M2M(CellClass* const FRestrict inPole, const CellClass *const FRestrict *const FRestrict inChild, const int inLevel) override {
-        FComplex<FReal>* FRestrict const multipole_exp_target = inPole->getMultipole();
+    template<class SymbolicData>
+    void M2M(typename CellClass::multipole_t* const FRestrict ParentMultipole,
+             const SymbolicData* const ParentSymb,
+             const typename CellClass::multipole_t*const FRestrict *const FRestrict ChildMultipoles,
+             const SymbolicData* const /*ChildSymbs*/[])
+    {
+        int level = static_cast<int>(ParentSymb->getLevel());
+        FComplex<FReal>* FRestrict const multipole_exp_target = ParentMultipole->get();
         // iter on each child and process M2M
-        const FComplex<FReal>* FRestrict const preM2MTransitionsAtLevel = preM2MTransitions[inLevel];
+        const FComplex<FReal>* FRestrict const preM2MTransitionsAtLevel
+            = preM2MTransitions[level];
         for(int idxChild = 0 ; idxChild < 8 ; ++idxChild){
-            if(inChild[idxChild]){
-                multipoleToMultipole(multipole_exp_target, inChild[idxChild]->getMultipole(), &preM2MTransitionsAtLevel[idxChild * harmonic.getExpSize()]);
+            if(ChildMultipoles[idxChild]){
+                multipoleToMultipole(multipole_exp_target, ChildMultipoles[idxChild]->get(), &preM2MTransitionsAtLevel[idxChild * harmonic.getExpSize()]);
             }
         }
     }
 
     /** M2L with a cell and all the existing neighbors */
-    virtual void M2L(CellClass* const FRestrict inLocal, const CellClass* inInteractions[],
-                     const int neighborPositions[], const int inSize, const int inLevel) override
-        = 0;
+    // virtual void M2L(CellClass* const FRestrict inLocal, const CellClass* inInteractions[],
+    //                  const int neighborPositions[], const int inSize, const int inLevel) override
+    //     = 0;
 
     /** L2L with a cell and all its child */
-    void L2L(const CellClass* const FRestrict pole, CellClass* FRestrict *const FRestrict child, const int inLevel) override {
-        // iter on each child and process L2L
-        const FComplex<FReal>* FRestrict const preL2LTransitionsAtLevel = preL2LTransitions[inLevel];
+    template<class SymbolicData>
+    void L2L(const typename CellClass::local_expansion_t * const FRestrict ParentExpansion,
+             const SymbolicData* const ParentSymb,
+             typename CellClass::local_expansion_t * FRestrict * const FRestrict ChildExpansions,
+             const SymbolicData* const /*ChildSymbs*/[])
+    {
+        int level = static_cast<int>(ParentSymb->getLevel());
+        // iter on each ChildrenExpansions and process L2L
+        const FComplex<FReal>* FRestrict const preL2LTransitionsAtLevel
+            = preL2LTransitions[level];
         for(int idxChild = 0 ; idxChild < 8 ; ++idxChild){
-            if(child[idxChild]){
-                localToLocal(child[idxChild]->getLocal(), pole->getLocal(), &preL2LTransitionsAtLevel[idxChild * harmonic.getExpSize()]);
+            if(ChildExpansions[idxChild]){
+                localToLocal(ChildExpansions[idxChild]->get(),
+                             ParentExpansion->get(),
+                             &preL2LTransitionsAtLevel[idxChild * harmonic.getExpSize()]);
             }
         }
     }
 
     /** L2P with a cell and all its particles */
-    void L2P(const CellClass* const local, ContainerClass* const inParticles)override {
-        const FComplex<FReal>* const cellLocal = local->getLocal();
+    template<class SymbolicData>
+    void L2P(const typename CellClass::local_expansion_t* const local,
+             const SymbolicData* const LeafSymb,
+             ContainerClass* const inParticles)
+    {
+        const FComplex<FReal>* const cellLocal = local->get();
         // Copying the position is faster than using cell position
-        const FPoint<FReal> localPosition = getLeafCenter(local->getCoordinate());
+        const FPoint<FReal> localPosition = getLeafCenter(LeafSymb->getCoordinate());
         // For all particles in the leaf box
         const FReal*const physicalValues = inParticles->getPhysicalValues();
         const FReal*const positionsX = inParticles->getPositions()[0];

Src/Kernels/Spherical/FSphericalBlasKernel.hpp

@@ -181,15 +181,24 @@ public:
 	}
 
 	/** M2L with a cell and all the existing neighbors */
-    void M2L(CellClass* const FRestrict inLocal, const CellClass* distantNeighbors[],
-             const int neighborPositions[], const int inSize, const int inLevel)  override {
-		// For all neighbors compute M2L
+    template<class SymbolicData>
+    void M2L(typename CellClass::local_expansion_t * const FRestrict TargetExpansion,
+             const SymbolicData* const TargetSymb,
+             const typename CellClass::multipole_t * const FRestrict SourceMultipoles[],
+             const SymbolicData* const FRestrict /*SourceSymbs*/[],
+             const int neighborPositions[],
+             const int inSize)
+    {
+        int level = static_cast<int>(TargetSymb->getLevel());
+        // For all neighbors compute M2L
         for(int idxExistingNeigh = 0 ; idxExistingNeigh < inSize ; ++idxExistingNeigh){
             const int idxNeigh = neighborPositions[idxExistingNeigh];
-            const FComplex<FReal>* const transitionVector = preM2LTransitions[inLevel][idxNeigh];
-            multipoleToLocal(inLocal->getLocal(), distantNeighbors[idxExistingNeigh]->getMultipole(), transitionVector);
+            const FComplex<FReal>* const transitionVector = preM2LTransitions[level][idxNeigh];
+            multipoleToLocal(TargetExpansion->get(),
+                             SourceMultipoles[idxExistingNeigh]->get(),
+                             transitionVector);
         }
-	}
+    }
 
 
 	/** preExpNExp
@@ -240,22 +249,26 @@ public:
 	 *       }
 	 *   }
 	 */
-    void multipoleToLocal(FComplex<FReal>*const FRestrict local_exp, const FComplex<FReal>* const FRestrict multipole_exp_src,
-            const FComplex<FReal>* const FRestrict M2L_Outer_transfer){
+    void multipoleToLocal(FComplex<FReal>*const FRestrict local_exp,
+                          const FComplex<FReal>* const FRestrict multipole_exp_src,
+                          const FComplex<FReal>* const FRestrict M2L_Outer_transfer)
+    {
 		// Copy original vector and compute exp2nexp
-        FMemUtils::copyall<FComplex<FReal>>(temporaryMultiSource, multipole_exp_src, CellClass::GetPoleSize());
-		// Get a computable vector
-		preExpNExp(temporaryMultiSource);
-
-		const FReal one[2] = {1.0 , 0.0};
-
-		FBlas::c_gemva(
-				FF_MATRIX_ROW_DIM,
-				FF_MATRIX_COLUMN_DIM,
-				one,
-                FComplex<FReal>::ToFReal(M2L_Outer_transfer),
-                FComplex<FReal>::ToFReal(temporaryMultiSource),
-                FComplex<FReal>::ToFReal(local_exp));
+        FMemUtils::copyall<FComplex<FReal>>(temporaryMultiSource,
+                                            multipole_exp_src,
+                                            CellClass::multipole_t::getSize());
+        // Get a computable vector
+        preExpNExp(temporaryMultiSource);
+
+        const FReal one[2] = {1.0 , 0.0};
+
+        FBlas::c_gemva(
+            FF_MATRIX_ROW_DIM,
+            FF_MATRIX_COLUMN_DIM,
+            one,
+            FComplex<FReal>::ToFReal(M2L_Outer_transfer),
+            FComplex<FReal>::ToFReal(temporaryMultiSource),
+            FComplex<FReal>::ToFReal(local_exp));
 
 // #ifdef DEBUG_SPHERICAL_M2L
 

Src/Kernels/Spherical/FSphericalBlockBlasKernel.hpp

@@ -175,15 +175,23 @@ public:
     }
 
     /** M2L with a cell and all the existing neighbors */
-    void M2L(CellClass* const FRestrict inLocal, const CellClass* distantNeighbors[],
-             const int neighborPositions[], const int inSize, const int inLevel)  override {
+    template<class SymbolicData>
+    void M2L(typename CellClass::local_expansion_t * const FRestrict TargetExpansion,
+             const SymbolicData* const TargetSymb,
+             const typename CellClass::multipole_t * const FRestrict SourceMultipoles[],
+             const SymbolicData* const FRestrict /*SourceSymbs*/[],
+             const int neighborPositions[],
+             const int inSize)
+    {
         // For all neighbors compute M2L
         for(int idxExistingNeigh = 0 ; idxExistingNeigh < inSize ; ++idxExistingNeigh){
             const int idxNeigh = neighborPositions[idxExistingNeigh];
-            interactions[idxNeigh].push(ComputationPair(distantNeighbors[idxExistingNeigh]->getMultipole(), inLocal->getLocal()));
+            interactions[idxNeigh].push(
+                ComputationPair(SourceMultipoles[idxExistingNeigh]->get(),
+                                TargetExpansion->get()));
 
             if( interactions[idxNeigh].getSize() == BlockSize){
-                multipoleToLocal( idxNeigh, inLevel);
+                multipoleToLocal( idxNeigh, static_cast<int>(TargetSymb->getLevel()));
             }
         }
     }

Src/Kernels/Spherical/FSphericalCell.hpp

@@ -4,13 +4,13 @@
 // This software is a computer program whose purpose is to compute the FMM.
 //
 // This software is governed by the CeCILL-C and LGPL licenses and
-// abiding by the rules of distribution of free software.  
-// 
+// abiding by the rules of distribution of free software.
+//
 // This program is distributed in the hope that it will be useful,
 // but WITHOUT ANY WARRANTY; without even the implied warranty of
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 // GNU General Public and CeCILL-C Licenses for more details.
-// "http://www.cecill.info". 
+// "http://www.cecill.info".
 // "http://www.gnu.org/licenses".
 // ===================================================================================
 #ifndef FSPHERICALCELL_HPP
@@ -30,12 +30,81 @@ template <class FReal>
 class FSphericalCell : public FBasicCell, public FAbstractSendable {
 protected:
     static int DevP;
-    static int LocalSize;
-    static int PoleSize;
     static bool UseBlas;
 
-    FComplex<FReal>* multipole_exp; //< For multipole extenssion
-    FComplex<FReal>* local_exp;     //< For local extenssion
+    template<class PoleTag>
+    struct exp_impl {
+        using cell_t = FSphericalCell;
+        static int Size;
+
+        static int getSize() noexcept {
+            return Size;
+        }
+
+        FComplex<FReal>* exp;
+
+        exp_impl() : exp{nullptr} {
+            this->exp = new FComplex<FReal>[Size];
+        }
+
+        exp_impl(const exp_impl& other) : exp_impl() {
+            *this = other;
+        }
+
+        exp_impl(exp_impl&& other) : exp(nullptr) {
+            *this = std::move(other);
+        }
+
+        exp_impl& operator=(const exp_impl& other) {
+            FMemUtils::copyall(exp, other.exp, Size);
+            return *this;
+        }
+
+        exp_impl& operator=(exp_impl&& other) {
+            using std::swap;
+            swap(this->exp, other.exp);
+            return *this;
+        }
+
+        ~exp_impl() {
+            delete[] this->exp;
+        }
+
+        FComplex<FReal>* get() noexcept {
+            return exp;
+        }
+        const FComplex<FReal>* get() const noexcept {
+            return exp;
+        }
+
+        void reset() noexcept {
+            for(int idx = 0; idx < Size; ++idx) {
+                exp[idx].setRealImag(FReal(0.0), FReal(0.0));
+            }
+        }
+
+        template<class BufferWriterClass>
+        void serialize(BufferWriterClass& buffer) const {
+            buffer.write(exp, Size);
+        }
+        template<class BufferReaderClass>
+        void deserialize(BufferReaderClass& buffer) {
+            buffer.fillArray(exp, Size);
+        }
+    };
+
+
+public:
+
+    using multipole_t = exp_impl<struct multipole_tag>;
+    using local_expansion_t = exp_impl<struct local_expansion_tag>;
+
+
+
+protected:
+
+    multipole_t m_data {};
+    local_expansion_t l_data {};
 
 public:
     static void Init(const int inDevP, const bool inUseBlas = false){
@@ -43,136 +112,106 @@ public:
         const int ExpP  = int((inDevP+1) * (inDevP+2) * 0.5);
         const int NExpP = (inDevP+1) * (inDevP+1);
 
-        LocalSize = ExpP;
+        local_expansion_t::Size = ExpP;
         if(inUseBlas) {
-            PoleSize = NExpP;
+            multipole_t::Size = NExpP;
         }
         else{
-            PoleSize = ExpP;
+            multipole_t::Size = ExpP;
         }
     }
 
-    static int GetLocalSize(){
-        return LocalSize;
-    }
-
-    static int GetPoleSize(){
-        return PoleSize;
-    }
 
     /** Default constructor */
-    FSphericalCell()
-        : multipole_exp(nullptr), local_exp(nullptr){
-        multipole_exp = new FComplex<FReal>[PoleSize];
-        local_exp = new FComplex<FReal>[LocalSize];
-    }
-
-    /** Constructor */
-    FSphericalCell(const FSphericalCell<FReal>& other)
-        : multipole_exp(nullptr), local_exp(nullptr){
-        multipole_exp = new FComplex<FReal>[PoleSize];
-        local_exp = new FComplex<FReal>[LocalSize];
-        (*this) = other;
-    }
-
-    /** Default destructor */
-    virtual ~FSphericalCell(){
-        delete[] multipole_exp;
-        delete[] local_exp;
-    }
-
+    FSphericalCell() : m_data(), l_data() {}
     /** Copy constructor */
-    FSphericalCell<FReal>& operator=(const FSphericalCell<FReal>& other) {
-        FMemUtils::copyall(multipole_exp, other.multipole_exp, PoleSize);
-        FMemUtils::copyall(local_exp, other.local_exp, LocalSize);
-        return *this;
-    }
+    FSphericalCell(const FSphericalCell& other) = default;
+    /** Move constructor */
+    FSphericalCell(FSphericalCell&& other) = default;
+    /** Copy operator */
+    FSphericalCell& operator=(const FSphericalCell&) = default;
+    /** Move operator */
+    FSphericalCell& operator=(FSphericalCell&&) = default;
+    /** Default destructor */
+    virtual ~FSphericalCell(){}
 
-    /** Get Multipole */
-    const FComplex<FReal>* getMultipole() const {
-        return multipole_exp;
+    multipole_t& getMultipoleData() noexcept {
+        return m_data;
     }
-    /** Get Local */
-    const FComplex<FReal>* getLocal() const {
-        return local_exp;
+    const multipole_t& getMultipoleData() const noexcept {
+        return m_data;
     }
 
-    /** Get Multipole */
-    FComplex<FReal>* getMultipole() {
-        return multipole_exp;
+    local_expansion_t& getLocalExpansionData() noexcept {
+        return l_data;
     }
-    /** Get Local */
-    FComplex<FReal>* getLocal() {
-        return local_exp;
+    const local_expansion_t& getLocalExpansionData() const noexcept {
+        return l_data;
     }
 
+
     /** Make it like the begining */
     void resetToInitialState(){
-        for(int idx = 0 ; idx < PoleSize ; ++idx){
-            multipole_exp[idx].setRealImag(FReal(0.0), FReal(0.0));
-        }
-        for(int idx = 0 ; idx < LocalSize ; ++idx){
-            local_exp[idx].setRealImag(FReal(0.0), FReal(0.0));
-        }
+        m_data.reset();
+        l_data.reset();
     }
 
     ///////////////////////////////////////////////////////
     // to extend FAbstractSendable
     ///////////////////////////////////////////////////////
     template <class BufferWriterClass>
-    void serializeUp(BufferWriterClass& buffer) const{
-        buffer.write(multipole_exp, PoleSize);
+    void serializeUp(BufferWriterClass& buffer) const {
+        m_data.serialize(buffer);
     }
     template <class BufferReaderClass>
     void deserializeUp(BufferReaderClass& buffer){
-        buffer.fillArray(multipole_exp, PoleSize);
+        m_data.deserialize(buffer);
     }
 
     template <class BufferWriterClass>
-    void serializeDown(BufferWriterClass& buffer) const{
-        buffer.write(local_exp, LocalSize);
+    void serializeDown(BufferWriterClass& buffer) const {
+        l_data.serialize(buffer);
     }
     template <class BufferReaderClass>
     void deserializeDown(BufferReaderClass& buffer){
-        buffer.fillArray(local_exp, LocalSize);
+        l_data.deserialize(buffer);
     }
 
     ///////////////////////////////////////////////////////
     // to extend Serializable
     ///////////////////////////////////////////////////////
     template <class BufferWriterClass>
-    void save(BufferWriterClass& buffer) const{
+    void save(BufferWriterClass& buffer) const {
         FBasicCell::save(buffer);
-        buffer.write(multipole_exp, PoleSize);
-        buffer.write(local_exp, LocalSize);
+        m_data.serialize(buffer);
+        l_data.serialize(buffer);
     }
     template <class BufferReaderClass>
     void restore(BufferReaderClass& buffer){
         FBasicCell::restore(buffer);
-        buffer.fillArray(multipole_exp, PoleSize);
-        buffer.fillArray(local_exp, LocalSize);
+        m_data.deserialize(buffer);
+        l_data.deserialize(buffer);
     }
 
     FSize getSavedSize() const {
-        return ((FSize) sizeof(FComplex<FReal>)) * (PoleSize+LocalSize)
+        return ((FSize) sizeof(FComplex<FReal>)) * (multipole_t::Size+local_expansion_t::Size)
                 + FBasicCell::getSavedSize();
     }
 
     FSize getSavedSizeUp() const {
-        return ((FSize) sizeof(FComplex<FReal>)) * (PoleSize);
+        return ((FSize) sizeof(FComplex<FReal>)) * (multipole_t::Size);
     }
 
     FSize getSavedSizeDown() const {
-        return ((FSize) sizeof(FComplex<FReal>)) * (LocalSize);
+        return ((FSize) sizeof(FComplex<FReal>)) * (local_expansion_t::Size);
     }
 };
 
 template <class FReal>
 int FSphericalCell<FReal>::DevP(-1);
-template <class FReal>
-int FSphericalCell<FReal>::LocalSize(-1);
-template <class FReal>
-int FSphericalCell<FReal>::PoleSize(-1);
+
+template <class FReal> template<class Tag>
+int FSphericalCell<FReal>::exp_impl<Tag>::Size(-1);
 
 
 /**
@@ -204,5 +243,3 @@ public:
 
 
 #endif //FSPHERICALCELL_HPP
-
-

Src/Kernels/Spherical/FSphericalKernel.hpp

@@ -95,15 +95,20 @@ public:
     }
 
     /** M2L with a cell and all the existing neighbors */
-    void M2L(CellClass* const FRestrict inLocal, const CellClass* distantNeighbors[],
-             const int neighborPositions[], const int inSize, const int inLevel)
-        override
+    template<class SymbolicData>
+    void M2L(typename CellClass::local_expansion_t * const FRestrict TargetExpansion,
+             const SymbolicData* const TargetSymb,
+             const typename CellClass::multipole_t * const FRestrict SourceMultipoles[],
+             const SymbolicData* const FRestrict /*SourceSymbs*/[],
+             const int neighborPositions[],
+             const int inSize)
     {
         // For all neighbors compute M2L
+        int level = static_cast<int>(TargetSymb->getLevel());
         for(int idxExistingNeigh = 0 ; idxExistingNeigh < inSize ; ++idxExistingNeigh){
             const int idxNeigh = neighborPositions[idxExistingNeigh];
-            const FComplex<FReal>* const transitionVector = &preM2LTransitions[inLevel][idxNeigh * devM2lP];
-            multipoleToLocal(inLocal->getLocal(), distantNeighbors[idxExistingNeigh]->getMultipole(), transitionVector);
+            const FComplex<FReal>* const transitionVector = &preM2LTransitions[level][idxNeigh * devM2lP];
+            multipoleToLocal(TargetExpansion->get(), SourceMultipoles[idxExistingNeigh]->get(), transitionVector);
         }
     }
 

Src/Kernels/Spherical/FSphericalRotationKernel.hpp

@@ -419,13 +419,22 @@ public:
     }
 
     /** M2L with a cell and all the existing neighbors */
-    void M2L(CellClass* const FRestrict inLocal, const CellClass* distantNeighbors[],
-             const int neighborPositions[], const int inSize, const int inLevel)  override {
+    template<class SymbolicData>
+    void M2L(typename CellClass::local_expansion_t * const FRestrict TargetExpansion,
+             const SymbolicData* const TargetSymb,
+             const typename CellClass::multipole_t * const FRestrict SourceMultipoles[],
+             const SymbolicData* const FRestrict /*SourceSymbs*/[],
+             const int neighborPositions[],
+             const int inSize)
+    {
+        const int level = static_cast<int>(TargetSymb->getLevel());
         // For all neighbors compute M2L
         for(int idxExistingNeigh = 0 ; idxExistingNeigh < inSize ; ++idxExistingNeigh){
             const int idxNeigh = neighborPositions[idxExistingNeigh];
-            const RotationM2LTransfer& transitionVector = preM2LTransitions[inLevel][idxNeigh];
-            multipoleToLocal(inLocal->getLocal(), distantNeighbors[idxExistingNeigh]->getMultipole(), transitionVector);
+            const RotationM2LTransfer& transitionVector = preM2LTransitions[level][idxNeigh];
+            multipoleToLocal(TargetExpansion->get(),
+                             SourceMultipoles[idxExistingNeigh]->get(),
+                             transitionVector);
         }
     }