OMP_PtrField.h

#ifndef OMP_PtrField_H
#define OMP_PtrField_H
 
#include "CORE_Field.h"
 
//corresponding array header
#include "OMP_PtrArray.h"
 
template <typename T,typename K,K D>
class OMP_PtrField : public CORE_Field<T,K,D,OMP_PtrArray<T>,OMP_PtrField<T,K,D>> {
    
private:
    
    //This class type
    typedef OMP_PtrField<T,K,D>  Self;
    
    
    
    // CONSTRUCTORS
public:
    OMP_PtrField() : CORE_Field<T,K,D,OMP_PtrArray<T>,Self>() {
    }
    
   
    
    // DESTRUCTORS
    virtual ~OMP_PtrField() {
    }
    
    
private:
    
    
    
    
    //MEMORY
    //=====
public:
    virtual  tMemSize getMemorySize() const override {
        return sizeof(*this)+this->getContentsMemorySize();
    }
    
    //allocation methods
    //===================
    
public:
    static inline CORE_UniquePointer<Self> New() {
        return CORE_UniquePointer<Self>(new Self(),CORE_Object::Delete());
    }
    
   
    
    
    //accessor operators
    //====================
    
 
    //accessor iterator
    //=================
    
    
    //accessor methods
    //=================
    
    
    //assignment operators
    //====================
    inline  Self& operator=(const T& v)  {
        initialize(v);
        return *this;
    }
    inline Self& operator=(const std::initializer_list<T>& values)  {
        if (values.size()==D) {
            initialize(values);
        } else {
            copy(values);
        }
        return *this;
    }
    inline Self& operator=(std::initializer_list<T>&& values)  {
        if (values.size()==D) {
            initialize(values);
        } else {
            copy(values);
        }
        return *this;
    }
    
   
    template<size_t N,typename Q>
    inline  Self& operator=(const std::array<Q,N>& values)  {
        initialize(values);
        return *this;
        
    }
    
    template<typename Q>
    inline  Self& operator=(const std::valarray<Q>& values)  {
        copy(values);
        return *this;
    }
    
    template<typename Q>
    inline  Self& operator=(const std::vector<Q>& values)  {
        copy(values);
        return *this;
    }
 
    inline  Self& operator=(const Self& values)  {
        copy(values);
        return *this;
    }
    inline  Self& operator=(Self&& values)  {
        copy(values);
        return *this;
    }
 
    template<typename Q,class S1,class I1>
    inline  Self& operator=(const CORE_Field<Q,K,D,S1,I1>& values)  {
        copy(values);
        return *this;
    }
    template<typename Q,class S1,class I1>
    inline  Self& operator=(CORE_Field<Q,K,D,S1,I1>&& values)  {
        copy(values);
        return *this;
    }
    
   
     //initializers methods
    //====================
    inline void initialize(const T& a) {
        this->getStorage().initialize(a);
    }
    template<typename Q>
    inline void initialize(const std::initializer_list<Q>&& a) {
        //values of A
        auto vA =std::ranges::cbegin(a);
        
        //end iterator of A
        const auto iAe=std::ranges::cend(a);
        
        //begin iterator on this values
        auto values=this->getValues();
        
        //number of elements
        tInteger n=this->getElementsNumber();
        
        OMP_PARALLEL_SHARED((n,vA,iAe,values)) {
 
            // thread id
            tInteger threadId=OMP_GET_THREAD_ID();
            
            //threads number
            tInteger nThreads=OMP_GET_THREADS_NUMBER();
 
            //start index of element
            tInteger start=threadId*n/nThreads;
 
            //end index of element
            tInteger end=(threadId+1)*n/nThreads;
 
            //start iterator on values
            T* iV=values+D*start;
            
            //end iterator on values
            const T* iVe=values+D*end;
            
            //iterator on A
            auto iA=vA;
            
            while (iV!=iVe) {//loop on i
 
                iA=vA;
                while (iA!=iAe) {//loop on d
 
                    //Vi[d]=A[d];
                    (*iV)=(*iA);
                    
                    iA++;
                    iV++;
                }
 
              
            }
 
          
        }
        
    }
    template<typename Q>
    inline void initialize(const std::initializer_list<Q>& a) {
        //values of a
        auto vA =std::ranges::cbegin(a);
        
        //end iterator of A
        const auto iAe=std::ranges::cend(a);
        
        //begin iterator on this values
        auto values=this->getValues();
        
        //number of elements
        tInteger n=this->getElementsNumber();
        
        OMP_PARALLEL_SHARED((n,vA,iAe,values)) {
 
            // thread id
            tInteger threadId=OMP_GET_THREAD_ID();
            
            //threads number
            tInteger nThreads=OMP_GET_THREADS_NUMBER();
 
            //start index of element
            tInteger start=threadId*n/nThreads;
 
            //end index of element
            tInteger end=(threadId+1)*n/nThreads;
 
            //start iterator on values
            T* iV=values+D*start;
            
            //end iterator on values
            const T* iVe=values+D*end;
            
            //iterator on A
            auto iA=vA;
            
            while (iV!=iVe) {//loop on i
 
                iA=vA;
                while (iA!=iAe) {//loop on d
 
                    //Vi[d]=A[d];
                    (*iV)=(*iA);
                    
                    iA++;
                    iV++;
                }
 
              
            }
 
          
        }
        
    }
    template<typename Q>
    inline void initialize(const std::array<Q,D>& a) {
 
        //values of A
        const Q* vA =a.data();
        
        //end iterator of A
        const Q* iAe=vA+D;
        
        //begin iterator on this values
        auto values=this->getValues();
        
        //number of elements
        tInteger n=this->getElementsNumber();
        
        OMP_PARALLEL_SHARED((n,vA,iAe,values)) {
 
            // thread id
            tInteger threadId=OMP_GET_THREAD_ID();
            
            //threads number
            tInteger nThreads=OMP_GET_THREADS_NUMBER();
 
            //start index of element
            tInteger start=threadId*n/nThreads;
 
            //end index of element
            tInteger end=(threadId+1)*n/nThreads;
 
            //start iterator on values
            T* iV=values+D*start;
            
            //end iterator on values
            const T* iVe=values+D*end;
            
            //iterator on A
            const Q* iA=null;
            
            while (iV!=iVe) {//loop on i
 
                iA=vA;
                while (iA!=iAe) {//loop on d
 
                    //Vi[d]=A[d];
                    (*iV)=(*iA);
                    
                    iA++;
                    iV++;
                }
            }
 
          
        }
    }
    
    
    
 
    
    
    //compound assignement operators (+=,-=,*=,/=,^=...)
    //===================================================
    //arithmetic type compound operators
    //====================================
    inline Self& operator+=(const T& v) {
        this->getStorage()+=v;
        return (*this);
    }
    inline Self& operator-=(const T& v) {
        this->getStorage()-=v;
        return (*this);
      
    }
    
    inline Self& operator*=(const T& v) {
        this->getStorage()*=v;
        return (*this);
        
    }
    inline Self& operator/=(const T& v) {
        this->getStorage()/=v;
        return (*this);
    }
 
    //integer type compond operators
    //===============================
    inline Self& operator%=(const T& v) requires functions_type::isIntegerType<T>{
        this->getStorage()%=v;
        return (*this);
    }
    
    inline Self& operator&=(const T& v) requires functions_type::isIntegerType<T>{
        this->getStorage()&=v;
        return (*this);
    }
    
    inline Self& operator|=(const T& v) requires functions_type::isIntegerType<T>{
        this->getStorage()|=v;
        return (*this);
    }
    
    inline Self& operator^=(const T& v) requires functions_type::isIntegerType<T> {
         this->getStorage()^=v;
        return (*this);
    }
    inline  Self& operator<<=(const T& v) requires functions_type::isIntegerType<T> {
        this->getStorage()<<=v;
        return (*this);
    }
    
    inline Self& operator>>=(const T& v) requires functions_type::isIntegerType<T> {
        this->getStorage()>>=v;
        return (*this);
    }
    
    //class type compound operatos
    //=============================
    
    template<class Q,class S,class I>
    inline Self& operator+=(const CORE_Field<Q,K,D,S,I>& X) {
        this->getStorage()+=X.getStorage();
        return (*this);
    }
    template<class Q,class S,class I>
    inline Self& operator-=(const CORE_Field<Q,K,D,S,I>& X) {
        this->getStorage()-=X.getStorage();
        return (*this);
    }
    template<class Q,class S,class I>
    inline Self& operator*=(const CORE_Field<Q,K,D,S,I>& X) {
        //T=F(T,V)
        this->getStorage()*=X.getStorage();
        return (*this);
    }
    template<class Q,class S,class I>
    inline Self& operator/=(const CORE_Field<Q,K,D,S,I>& X) {
        this->getStorage()/=X.getStorage();
        return (*this);
    }
    
 
    //transform methods
    //=================
    
   
 
    template<typename LambdaFct>
    static inline void ElementsTransform(LambdaFct&& F,
                                         Self& X) {
 
        //number of elements
        tIndex n=X.getElementsNumber();
       
         
        //get the const begin stride iterator of X
        auto iX=X.sbegin();
        
     
        //init the random data
        OMP_PARALLEL_SHARED((n,iX,F)) {//begin parallel section
            
            // thread id
            tInteger threadId=OMP_GET_THREAD_ID();
            //threads number
            tInteger nThreads=OMP_GET_THREADS_NUMBER();
 
            //start index of element
            tInteger start=threadId*n/nThreads;
            
            //end index of element
            tInteger end=(threadId+1)*n/nThreads;
 
 
            std::transform(iX+start,iX+end,
                           iX+start,F);
            
            
        }//end parallel section
        
        
    }
    
  
    template<typename LambdaFct>
    inline void elementsTransform(LambdaFct&& F) {
        this->ElementsTransform(F,(*this));
    }
    
   
   
    inline void normalize()  {
      
        //error to detect null
        T eps=std::numeric_limits<T>::epsilon();
 
        
        //number of elements
        tIndex n=this->getElementsNumber();
       
         
        //get values of this
        T*  values=this->getValues();
        
     
        //init the random data
        OMP_PARALLEL_SHARED((n,values,eps)) {//begin parallel section
            
            // thread id
            tInteger threadId=OMP_GET_THREAD_ID();
            //threads number
            tInteger nThreads=OMP_GET_THREADS_NUMBER();
            //start index of element
            tInteger start=threadId*n/nThreads;
            //end index of element
            tInteger end=(threadId+1)*n/nThreads;
 
            //iterator of index of the values of the start element 
            T* iV=values+D*start;
            
            //iterator of index of the values of the end element 
            const T* iVe=values+D*end;
            
            //begin iterator on dimension ie iVd=iV
            const T* iVd=null;
            
            //end iterator on dimension
            const T* iVde=null;
            
            //norm2
            T norm2=0;
            while (iV!=iVe)  {//loop on element i
 
 
                //iterator on end of element i
                iVde=iV+D;
 
                //compute the norm at element iV
                iVd=iV;
                norm2=0;
                while (iVd!=iVde) {
                    norm2+=(*iVd)*(*iVd);
                    iVd++;
                }
                //normlize
                if (norm2>eps) {
                    norm2=1./sqrt(norm2);
                    while (iV!=iVde) {
                        (*iV)*=norm2;
                        iV++;
                    }
                } else {
                    //go to next element
                    iV+=D;
                }
            }
            
        }//end parallel section
        
    }
    
   
    //Data consistency
    //=================
   
 
    //transform methods
    //==================
  
 
    //algebric methods
    //================
    
   
     
    template<class I>
    inline void mod2(CORE_Array<T,I>& X) const {
 
           
        //number of elements
        tIndex n=this->getElementsNumber();
 
        //resize x to the number of elements
        X.setSize(n);
        
        //get values of this
        const T* values=this->getValues();
 
        //get the  values of  X
        T* vX=X.getValues();
        
     
        //init the random data
        OMP_PARALLEL_SHARED((n,vX,values)) {//begin parallel section
            
            
            // thread id
            tInteger threadId=OMP_GET_THREAD_ID();
            //threads number
            tInteger nThreads=OMP_GET_THREADS_NUMBER();
            //start element index
            tInteger start=threadId*n/nThreads;
            //end element index
            tInteger end=(threadId+1)*n/nThreads;
            
            // iterator on values of element  start
            const T* iV=values+start*D;
 
            //iterator on values of element  send
            const T* iVe=values+end*D;
 
            //end iterator on element i
            const T* iVde=null;
 
            //start iterator on X
            T* iX=vX+start;
 
            while (iV!=iVe) {//loop on element
 
                (*iX)=0;
                iVde=iV+D;//end iterator on dimension
                while (iV!=iVde) {//loopon dimension
                    (*iX)+=(*iV)*(*iV);
                    iV++;
                }
                
                iX++;
                
            }//end loop on element
 
        }//end parallel section
       
            
        
    }//end norm2(X) method 
 
    template<class Q,class S1,class I1>
    inline T& scalarProduct(const CORE_Field<Q,K,D,S1,I1>& X,T& s)  const  {
        return this->getStorage().scalarProduct(X.getStorage(),s);
    }
    template<class Q,class S1,class I1>
    inline T& scalarProduct(const std::valarray<Q>& weights,const CORE_Field<Q,K,D,S1,I1>& X,T& s)  const  {
        
        //number of elements
        tIndex n=this->getElementsNumber();
 
         //increment on X
        tBoolean incXs=(X.getElementsNumber()!=1);
 
        
        //s=sum_i  Xi Yi
        s=0;
 
        if (weights.size()==0) {
            return CORE_Field<T,K,D,OMP_PtrArray<T>,Self>::scalarProduct(X,s);
        }
        if ((weights.size()==1)&& (incXs)) {
            CORE_Field<T,K,D,OMP_PtrArray<T>,Self>::scalarProduct(X,s);
            s*=weights[0];
            return s;
        }
 
       
        //values of W
        const tReal *Ws=&weights[0];
        tBoolean incW=(weights.size()!=1);
 
        //values of this
        const T* Ys=this->getValues();
        
        //values of X
        const Q* Xs=X.getValues();
 
       
        
        OMP_PARALLEL_SHARED_REDUCTION((n,incW,Ws,Ys,Xs,incXs),
                                      (+:s)) {//begin parallel section
 
             
            // thread id
            tInteger threadId=OMP_GET_THREAD_ID();
            // threads number
            tInteger nThreads=OMP_GET_THREADS_NUMBER();
 
            //start idnex of the loop
            tIndex start=threadId*n/nThreads;
 
            
            //begin weight iterator
            const Q* iW=Ws;
            iW+=incW*start;
 
          
            
            start*=D;
 
            //start iterator on X
            const Q* iX=Xs;
            iX+=incXs*start;
            const Q* iXd=null;//iterator on d-coordinate
                 
                 
            //start iterator on Y
            const T* iY=Ys;
            iY+=start;
            
            //end index oft the loop
            tIndex end=(threadId+1)*n/nThreads;;
            end*=D;
            
            
            //end  iterator on Y
            const T* iYe=Ys;
            iYe+=end;
            
            //end iterator on Y for element
            const T* iYd=null;
 
            
            //s=wi <X^i.Y^i>_d
            s=0;
            while (iY!=iYe) {//loop on elements
 
                iYd=iY+D;
                iXd=iX;
                while(iY!=iYd) {//loop on <X^i,Y^i>
                    s+=(*iW)*(*iXd)*(*iY);
                    iY++;
                    iXd++;
                }//end loop on <X^i,Y^i>
 
                //weight of next elemnt
                iW+=incW;
                iX+=incXs*D;
            }//end loop on element
            
        }//end parallel loop
        
        return s;
    }
      
  
    
    
    
};
 
 
 
 
//defualt dimension
#ifndef DIMENSION_TYPE
#define  DIMENSION_TYPE
typedef  tUCInt tDimension;//d in [0,256[
#endif
 
typedef  OMP_PtrField<tReal,tDimension,3> OMP_Real3DPtrField;
 
 
 
//default field
#ifndef DEFAULT_OMP_FIELD
#define DEFAULT_OMP_FIELD
typedef OMP_Real3DPtrField OMP_Real3DField ;
#endif
 
#endif