OMP_ValField.h

#ifndef OMP_ValField_H
#define OMP_ValField_H
 
#include "CORE_Field.h"
 
//corresponding array header
#include "OMP_ValArray.h"
 
template <typename T,typename K,K D>
class OMP_ValField : public CORE_Field<T,K,D,OMP_ValArray<T>,OMP_ValField<T,K,D>> {
    
private:
    
    //This class type
    typedef OMP_ValField<T,K,D>  Self;
    
    
    // CONSTRUCTORS
public:
    OMP_ValField() : CORE_Field<T,K,D,OMP_ValArray<T>,Self>() {
    }
    
   
    
    // DESTRUCTORS
    virtual ~OMP_ValField() {
    }
    
    
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)  {
        copy(v);
        return *this;
    }
    inline Self& operator=(const std::initializer_list<T>& values)  {
        copy(values);
        return *this;
    }
    inline Self& operator=(std::initializer_list<T>&& values)  {
        copy(values);
        return *this;
    }
    
   
    template<size_t N,typename Q>
    inline  Self& operator=(const std::array<Q,N>& values)  {
        copy(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
    //====================
    
 
    template<typename Q>
    inline void initialize(const std::array<Q,D>& a) {
 
        //begin iterator of a
        auto iA=a.cbegin();
 
        //begin iterator on this values
        auto iV=this->sbegin();
        
        //number of elements
        tInteger n=this->getElementsNumber();
        
        OMP_PARALLEL_SHARED((n,iA,iV)) {
 
            // thread id
            tInteger threadId=OMP_GET_THREAD_ID();
            //threads number
            tInteger nThreads=OMP_GET_THREADS_NUMBER();
 
            //start index
            tInteger start=threadId*n/nThreads;
 
            //end index
            tInteger end=(threadId+1)*n/nThreads;
            
            //compute Vi:=Vi/norm
            auto F=[&iA](auto& Vi) {           
                       //compute Vi:=Vi/norm
                       auto Vid=&Vi;
                       std::for_each_n(iA,D,[&Vid](const auto& Aid){(*Vid)=Aid;Vid++;});
                       
                       //return the value
                       return Vi;
                   };
            
            //apply the function F
            std::transform(iV+start,iV+end,
                           iV+start,F);
        }
    }
    
    
    
 
    
    
    //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 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();
            
            tInteger start=threadId*n/nThreads;
            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 the const begin iterator of X
        auto iV=this->sbegin();
        
     
        //init the random data
        OMP_PARALLEL_SHARED((n,iV,eps)) {//begin parallel section
            
            // thread id
            tInteger threadId=OMP_GET_THREAD_ID();
            //threads number
            tInteger nThreads=OMP_GET_THREADS_NUMBER();
            //start index
            tInteger start=threadId*n/nThreads;
            //end index
            tInteger end=(threadId+1)*n/nThreads;
 
            T norm=0;
            //compute Vi:=Vi/norm
            auto F=[&norm,&eps](auto& Vi) {
                       //compute Norm2
                       norm=0;
                       std::for_each_n(&Vi,D,[&norm](const auto& Vid){
                                                 norm+=Vid*Vid;
                                             });                        
                       
                       //compute invers of norm
                       if (norm>eps) {
                           norm=1./sqrt(norm);
                       }
                       
                       //compute Vi:=Vi/norm
                       std::for_each_n(&Vi,D,[&norm](auto& Vid){Vid*=norm;});
                       
                       //return the value
                       return Vi;
                   };
            std::transform(iV+start,iV+end,
                           iV+start,F);
            
        }//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
        X.setSize(n);
         
        //get the const begin iterator of this
        auto iV=this->csbegin();
 
        //get the  begin iterator of  X
        auto iX=X.begin();
        
     
        //init the random data
        OMP_PARALLEL_SHARED((n,iX,iV)) {//begin parallel section
            
            
            // thread id
            tInteger threadId=OMP_GET_THREAD_ID();
            //threads number
            tInteger nThreads=OMP_GET_THREADS_NUMBER();
            //start index
            tInteger start=threadId*n/nThreads;
            //end index
            tInteger end=(threadId+1)*n/nThreads;
 
            T norm2;
            std::transform(iV+start,iV+end,
                           iX+start,[&norm2](const auto& Vi) {
                                        //compute sum_d=0^D-1  Vid^2
                                        norm2=0;
                                        std::for_each_n(&Vi,D,[&norm2](const auto& Vid) {norm2+=Vid*Vid;});
                                        return norm2;
                                    });
            
            
        }//end parallel section
    }//end 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();
 
        //s=sum_i  Xi Yi
        s=0;
 
        if (weights.size()==0) {
            return CORE_Field<T,K,D,OMP_PtrArray<T>,Self>::scalarProduct(X,s);
        }
 
        //increment of W
        tBoolean incW=(weights.size()>1);
        //values of W
        const tReal *Ws=&weights[0];
 
        //values of this
        const T* Ys=this->getValues();
        
        //values of X
        const Q* Xs=X.getValues();
 
        //increment on X
        tBoolean incXs=(X.getElementsNumber()>1);
        
        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;
            
            //end index oft the loop
            tIndex end=(threadId+1)*n/nThreads;;
            
            //start iterator on X
            const Q* iX=Xs;iX+=incXs*start*D;
            const Q* iXd=null;//iterator on d-coordinate
            //start iterator on Y
            const T* iY=Ys;iY+=start*D;
            //end  iterator on Y
            const T* iYe=Ys;iYe+=end*D;
            //end iterator on Yfor element
            const T* iYd=null;
 
            //begin weight iterator
            const Q* iW=Ws;iW+=start*incW;
            
            //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;
    }
      
  
    
    
};
 
 
 
#ifndef DIMENSION_TYPE
#define  DIMENSION_TYPE
typedef  tUCInt tDimension;//d in [0,256[
#endif
 
typedef  OMP_ValField<tReal,tDimension,3> OMP_Real3DValField;
 
 
 
#endif