OMP_ValArray.h

#ifndef OMP_ValArray_H
#define OMP_ValArray_H
 
//base class header
#include "CORE_Object.h"
 
//exception header
#include "CORE_Exception.h"
 
//temmplate base class for array  header
#include "CORE_ValArray.h"
 
 
template<class T>
class OMP_ValArray : public CORE_ValArray<T,OMP_ValArray<T>> {
    
private:
    //This class type
    typedef OMP_ValArray<T>  Self;
    
public:
 
    
    
    OMP_ValArray() : CORE_ValArray<T,Self>() {
    }
 
 
    explicit OMP_ValArray(const tIndex& n) : CORE_ValArray<T,Self>(n) {
    }
 
    
    
   
 
    // DESTRUCTORS
    virtual ~OMP_ValArray() {
    }
 
 
 
     
    //MEMORY
    //=====
public:
    virtual  tMemSize getMemorySize() const override {
        return sizeof(*this)+this->getContentsMemorySize();
    }
 
    
    // ----------------
    // New constructors
    // ----------------
    static inline CORE_UniquePointer<Self> New() {
        return CORE_UniquePointer<Self>(new Self(),
                                        CORE_Object::Delete());
    };
    
public:
 
    
    //accessor operators
    //====================
    
 
    //accessor iterator
    //=================
 
    
    //accessor methods
    //=================
 
    //assignement operator
    //====================
   
public:
    
    inline  Self& operator=(const T& v)  {
        initialize(v);
        return *this;
    }
    inline Self& operator=(const std::initializer_list<T>& values)  {
        Copy(values.size(),values,this->getArray());
        return *this;
    }
    inline Self& operator=(std::initializer_list<T>&& values)  {
        Copy(values.size(),values,this->getArray());
        return *this;
    }
    
    template<size_t N>
    inline  Self& operator=(const std::array<T,N>& values)  {
        Copy(N,values,this->getArray());
        return *this;
        
    }
    
    inline  Self& operator=(const std::valarray<T>& values)  {
        Copy(values.size(),values,this->getArray());
        return *this;
    }
    
    inline  Self& operator=(const std::vector<T>& values)  {
        Copy(values.size(),values,this->getArray());
        return *this;
    }
    inline  Self& operator=(const Self& values)  {
        Copy(values,this->getArray());
        return *this;
    }
    inline  Self& operator=(Self&& values)  {
        Copy(values,this->getArray());
        return *this;
    }
 
    template<typename Q,class I1>
    inline  Self& operator=(const CORE_Array<Q,I1>& values)  {
        Copy(values,this->getArray());
        return *this;
    }
    template<typename Q,class I1>
    inline  Self& operator=(const CORE_Array<Q,I1>&& values)  {
        Copy(values,this->getArray());
        return *this;
    }
    
   
    
    
    //copy methods
    //============
public:
       
   
    template<typename Q,class I>
    inline void copy(const CORE_Array<Q,I>& X) {
        Copy(X,this->getArray());
    }
 
    template<typename Q,class I>
    inline void copy(CORE_Array<Q,I>&& X) {
        Copy(X,this->getArray());
    }
     
    template<typename Q,class I>
    inline void copy(const CORE_ValArray<Q,I>& X) {
        Copy(X.getArray(),this->getArray());
    }
   
    template<typename Q,class I>
    inline void copy(CORE_ValArray<Q,I>&& X)  {
        Copy(X.getArray(),this->getArray());
    }
    
   
    
    inline void copy(const Self& cpy) {
        Copy(cpy.getArray(),this->getArray());
    }
    
    inline void copy(Self&& cpy) {
       Copy(cpy.getArray(),this->getArray());
    }
 
  
 
    template<typename Q>
    inline void copy(const tIndex& n , const Q* Vs) {
        Copy(n,Vs,this->getArray());
    }
    
    inline void copy(const tIndex& n,const std::initializer_list<T>& Vs) {
        tIndex p=functions_numeric::min(n,Vs.size());
        Copy(p,Vs,this->getArray());
    }
    inline void copy(const tIndex& n,std::initializer_list<T>&& Vs) {
        tIndex p=functions_numeric::min(n,Vs.size());
        Copy(p,Vs,this->getArray());
    }
    template<typename Q,size_t N>
    inline void copy(const tIndex& n, const std::array<Q,N>& Vs) {
        tIndex p=functions_numeric::min(n,N);
        Copy(p,Vs.getData(),this->getArray());
    }
    
    template<typename Q>
    inline void copy(const tIndex& n, const std::valarray<Q>& Vs) {
        tIndex p=functions_numeric::min(n,Vs.size());
        Copy(p,Vs,this->getArray());
    }
    template<typename Q>
    inline void copy(const tIndex& n, const std::vector<Q>& Vs) {
        tIndex p=functions_numeric::min(n,Vs.size());
        Copy(p,Vs,this->getArray());
    }
   
 
 
    //static Copy
    //----------
    template<typename  Q,class I>
    inline static  void  Copy(const CORE_Array<Q,I>& X,std::valarray<T>& R) {
     
        //resize R
        tIndex p=X.getSize();
        if (R.size()!=p) R.resize(p);
        
        //get begin iterator of R
        auto iR=std::ranges::begin(R);
        
        //get const begin iterator on X
        auto iX=X.cbegin();
        
        OMP_PARALLEL_SHARED((p,iX,iR)) {//begin parallel section
            // thread id
            tInteger threadId=OMP_GET_THREAD_ID();
 
            // threads number
            tInteger nThreads=OMP_GET_THREADS_NUMBER();
            
            //start index
            tIndex start=threadId*p/nThreads;
            
            //end index
            tIndex end=(threadId+1)*p/nThreads;
 
            //iterator on R values at start
            auto iRi=iR;
            iRi+=start;
            
            //iterator on cpy values at start
            auto iXi=iX;
            iXi+=start;
 
            //iterator on cpy values at end
            auto  eXi=iX;
            eXi+=end;
            
            //copy the values of cpyA into aValue for element in [start,end[
            std::for_each(iXi,eXi,[&iRi](const auto& Xi) {(*iRi)=Xi;iRi++;});
            
        }//end parallel section
    }
 
 
    
    template<typename Q>
    inline static  void  Copy(const std::valarray<Q>& X,std::valarray<T>& R) {
        Copy(X.size(),X,R);
    }
    template<typename Q>
    inline static  void  Copy(const std::valarray<T>& X,std::valarray<T>& R) {
        R=X;
    }
    
    template<typename Q>
    inline static  void Copy(const tIndex& n,const std::valarray<Q>& X,std::valarray<T>& R) {
        Copy(n,&X[0],R);
    }
 
    template<typename Q>
    inline static  void Copy(const tIndex& n,const std::vector<Q>& X,std::valarray<T>& R) {
        //resize R
        tIndex p=functions_numeric::min(n,X.size());
        if (R.size()!=p) R.resize(p);
        
        //get begin iterator of R
        auto iR=std::ranges::begin(R);
        
        //get const begin iterator on X
        auto iX=X.cbegin();
        
        OMP_PARALLEL_SHARED((p,iX,iR)) {//begin parallel section
            // thread id
            tInteger threadId=OMP_GET_THREAD_ID();
 
            // threads number
            tInteger nThreads=OMP_GET_THREADS_NUMBER();
            
            //start index
            tIndex start=threadId*p/nThreads;
            
            //end index
            tIndex end=(threadId+1)*p/nThreads;
 
            //iterator on R values at start
            auto iRi=iR;
            iRi+=start;
            
            //iterator on cpy values at start
            auto iXi=iX;
            iXi+=start;
 
            //iterator on cpy values at end
            auto  eXi=iX;
            eXi+=end;
            
            //copy the values of cpyA into aValue for element in [start,end[
            std::for_each(iXi,eXi,[&iRi](const auto& Xi) {(*iRi)=Xi;iRi++;});
            
        }//end parallel section
    }
 
    
    inline static  void Copy(const tIndex& n,const std::initializer_list<T>& X,std::valarray<T>& R) {
        //resize R
        tIndex p=functions_numeric::min(n,X.size());
        if (R.size()!=p) R.resize(p);
        
        //get begin iterator of R
        auto iR=std::ranges::begin(R);
        
        //get const begin iterator on X
        auto iX=X.begin();
        
        OMP_PARALLEL_SHARED((p,iX,iR)) {//begin parallel section
            // thread id
            tInteger threadId=OMP_GET_THREAD_ID();
 
            // threads number
            tInteger nThreads=OMP_GET_THREADS_NUMBER();
            
            //start index
            tIndex start=threadId*p/nThreads;
            
            //end index
            tIndex end=(threadId+1)*p/nThreads;
 
            //iterator on R values at start
            auto iRi=iR;
            iRi+=start;
            
            //iterator on cpy values at start
            auto iXi=iX;
            iXi+=start;
 
            //iterator on cpy values at end
            auto  eXi=iX;
            eXi+=end;
            
            //copy the values of cpyA into aValue for element in [start,end[
            std::for_each(iXi,eXi,[&iRi](const auto& Xi) {(*iRi)=Xi;iRi++;});
            
        }//end parallel section
    }
    inline static  void Copy(const tIndex& n,std::initializer_list<T>&& X,std::valarray<T>& R) {
        //resize R
        tIndex p=functions_numeric::min(n,X.size());
        if (R.size()!=p) R.resize(p);
        
        //get begin iterator of R
        auto iR=std::ranges::begin(R);
        
        //get const begin iterator on X
        auto iX=X.begin();
        
        OMP_PARALLEL_SHARED((p,iX,iR)) {//begin parallel section
            // thread id
            tInteger threadId=OMP_GET_THREAD_ID();
 
            // threads number
            tInteger nThreads=OMP_GET_THREADS_NUMBER();
            
            //start index
            tIndex start=threadId*p/nThreads;
            
            //end index
            tIndex end=(threadId+1)*p/nThreads;
 
            //iterator on R values at start
            auto iRi=iR;
            iRi+=start;
            
            //iterator on cpy values at start
            auto iXi=iX;
            iXi+=start;
 
            //iterator on cpy values at end
            auto  eXi=iX;
            eXi+=end;
            
            //copy the values of cpyA into aValue for element in [start,end[
            std::for_each(iXi,eXi,[&iRi](const auto& Xi) {(*iRi)=Xi;iRi++;});
            
        }//end parallel section
    }
 
    
    template<typename Q>
    requires functions_type::isArithmeticType<Q>
    inline static  void Copy(const tIndex& n,const Q* X,std::valarray<T>& R) {
        R.setSize(n);
        
        //get begin iterator of R
        auto iR=std::ranges::begin(R);
       
        OMP_PARALLEL_SHARED((n,X,iR)) {//begin parallel section
            // thread id
            tInteger threadId=OMP_GET_THREAD_ID();
 
            // threads number
            tInteger nThreads=OMP_GET_THREADS_NUMBER();
            
            //start index
            tIndex start=threadId*n/nThreads;
            
            //end index
            tIndex end=(threadId+1)*n/nThreads;
 
            //iterator on R values at start
            auto iRi=iR;
            iRi+=start;
            
            //iterator on R values at end
            auto iRe=iR;
            iRe+=end;
            
            //iterator on X at start
            const Q* iX=X+start;
 
            
            //copy the values of cpyA into aValue for element in [start,end[
            std::for_each(iRi,iRe,[&iX](const auto& Ri) {Ri=(*iX);iX++;});
            
        }//end parallel section
        
    }
    
   
   
    //initializer
    //===========
  
    inline void initialize(const T& x) {
        Initialize(x,this->getArray());
    }
    inline void Initialize(const T& x,std::valarray<T>& values) {
        auto F=[&x](const auto &Ri){return x;};
        this->transform(F);
        
    }
    inline void uniformRandomize(const T& min,const T& max) {
        UniformRandomize(min,max,this->getArray());
    }
   
    
    inline static void UniformRandomize(const T& min,const T& max,std::valarray<T>& values){
        T alpha=(max-min)/RAND_MAX;
        
        //init the random data
        auto iR=std::ranges::begin(values);
 
        //number of elements
        tIndex n=values.size();
        
        OMP_PARALLEL_SHARED((alpha,n,min,iR)) {//begin parallel section
            
            // threads number
            tInteger nThreads=OMP_GET_THREADS_NUMBER();
            
            // thread id
            tInteger threadId=OMP_GET_THREAD_ID();
 
            tIndex start=0;
            
           //set random seed by process
            tULLInt seed=0;
            for(start=0;start<=threadId;start++) {
                seed=std::rand();
            }
            std::srand(seed);
            
            //start index
            start=threadId*n/nThreads;
            
            //end index
            tIndex end=(threadId+1)*n/nThreads;
 
            
 
            
            std::transform(iR+start,iR+end,
                           iR+start,[&alpha,&min](const auto& Ri) {
                                        return alpha*std::rand()+min;
                                    });
 
            //V3
            // auto iRi=iR+start;
            // auto iRe=iR+end;
            // while (iRi!=iRe) {
            //     (*iRi)=alpha*F(drand)+min;
            //     iRi++;
            // }
            
        }//end parallel section
        
    }
    
 
    //compount assignement operator
    //==============================
    
    inline OMP_ValArray<T>& operator+=(const T& v) {
        auto F=[&v](const auto &x){return x+v;};
        this->transform(F);
        return (*this);
    }
    inline OMP_ValArray<T>& operator-=(const T& v) {
        auto F=[&v](const auto& x){return x-v;};
        transform(F);
        return (*this);
    }
    
    inline OMP_ValArray<T>& operator*=(const T& v) {
        auto F=[&v](const auto &x){return x*v;};
        this->transform(F);
        return (*this);
       
    }
    inline OMP_ValArray<T>& operator/=(const T& v) {
        if (fabs(v)<std::numeric_limits<T>::epsilon()) {
            throw CORE_Exception("core",
                                 "OMP_ValArray::/=("+std::to_string(v)+")",
                                 "division by zero");
        }
        auto F=[&v](const auto &x){return x/v;};
        this->transform(F);
        return (*this);
    }
    
   
    //integer type compond operators
    //===============================
    inline Self& operator%=(const T& v) requires functions_type::isIntegerType<T>{
        auto F=[&v](const auto &x){return x%v;};
        this->transform(F);
        return (*this);
    }
    
    inline OMP_ValArray<T>& operator&=(const T& v) requires functions_type::isIntegerType<T>{
        auto F=[&v](const auto &x){return x&v;};
        this->transform(F);
        return (*this);
    }
    
    inline OMP_ValArray<T>& operator|=(const T& v) requires functions_type::isIntegerType<T> {
        auto F=[&v](const auto &x){return x|v;};
        this->transform(F);
        return (*this);
    }
 
    inline OMP_ValArray<T>& operator^=(const T& v) requires functions_type::isIntegerType<T> {
        auto F=[&v](const auto &x){return x^v;};
        this->transform(F);
        return (*this);
    }
    inline OMP_ValArray<T>& operator<<=(const T& v) requires functions_type::isIntegerType<T> {
        auto F=[&v](const auto &x){return x<<v;};
        this->transform(F);
        return (*this);
    }
    inline OMP_ValArray<T>& operator>>=(const T& v) requires functions_type::isIntegerType<T> {
        auto F=[&v](const auto &x){return x>>v;};
        this->transform(F);
        return (*this);
    }
 
    //class compount assignement operator
    //-----------------------------------
    template<typename Q,class I>
    inline Self& operator+=(const CORE_Array<Q,I>& X) {
        tIndex p=functions_numeric::min(X.getSize(),this->getSize());
        auto iX=X.cbegin();
        auto iV=this->begin();
        OMP_PARALLEL_SHARED((p,iX,iV)) {//begin parallel section
            // thread id
            tInteger threadId=OMP_GET_THREAD_ID();
            
            // threads number
            tInteger nThreads=OMP_GET_THREADS_NUMBER();
            
            //start index
            tIndex start=threadId*p/nThreads;
            
            //end index
            tIndex end=(threadId+1)*p/nThreads;
 
            std::transform(iX+start,iX+end,
                           iV+start,
                           iV+start,[](const auto& Xi,const auto &Vi) {return Vi+Xi;});
            
            
        }//end parallel section
        return *this;
    }
 
   
    template<typename Q,class I>
    inline Self& operator-=(const CORE_Array<Q,I>& X) {
        tIndex p=functions_numeric::min(X.getSize(),this->getSize());
        auto iX=X.cbegin();
        auto iV=this->begin();
        OMP_PARALLEL_SHARED((p,iX,iV)) {//begin parallel section
            // thread id
            tInteger threadId=OMP_GET_THREAD_ID();
            
            // threads number
            tInteger nThreads=OMP_GET_THREADS_NUMBER();
            
            //start index
            tIndex start=threadId*p/nThreads;
            
            //end index
            tIndex end=(threadId+1)*p/nThreads;
 
            std::transform(iX+start,iX+end,
                           iV+start,
                           iV+start,[](const auto& Xi,const auto &Vi) {return Vi-Xi;});
            
            
        }//end parallel section
        return (*this);
    }
    
    template<typename Q,class I>
    inline Self& operator*=(const CORE_Array<Q,I>& X) {
        tIndex p=functions_numeric::min(X.getSize(),this->getSize());
        auto iX=X.cbegin();
        auto iV=this->begin();
        OMP_PARALLEL_SHARED((p,iX,iV)) {//begin parallel section
            // thread id
            tInteger threadId=OMP_GET_THREAD_ID();
            
            // threads number
            tInteger nThreads=OMP_GET_THREADS_NUMBER();
            
            //start index
            tIndex start=threadId*p/nThreads;
            
            //end index
            tIndex end=(threadId+1)*p/nThreads;
 
            std::transform(iX+start,iX+end,
                           iV+start,
                           iV+start,[](const auto& Xi,const auto &Vi) {return Vi*Xi;});
            
            
        }//end parallel section
        return (*this);
    }
    template<typename Q,class I>
    inline Self& operator/=(const CORE_Array<Q,I>& X) {
        tIndex p=functions_numeric::min(X.getSize(),this->getSize());
        Q eps=std::numeric_limits<Q>::epsilon();
        auto iX=X.cbegin();
        auto iV=this->begin();
        OMP_PARALLEL_SHARED((p,eps,iX,iV)) {//begin parallel section
            // thread id
            tInteger threadId=OMP_GET_THREAD_ID();
            
            // threads number
            tInteger nThreads=OMP_GET_THREADS_NUMBER();
            
            //start index
            tIndex start=threadId*p/nThreads;
            
            //end index
            tIndex end=(threadId+1)*p/nThreads;
 
            std::transform(iX+start,iX+end,
                           iV+start,
                           iV+start,[&eps](const auto& Xi,const auto &Vi) {
                                        if (fabs(Xi)<eps) throw CORE_Exception("openmp",
                                                                               "OMP_ValArray/=",
                                                                               "division by 0");
                                        return Vi/Xi;});
            
            
        }//end parallel section
        return (*this);
    }
 
    
    //transform methods
    //=================
    template<typename LambdaFct>
    inline void transform(LambdaFct&& F) {
        Transform(F,(*this));
    }
    
    template<typename LambdaFct>
    inline void transform(LambdaFct&& F,
                          const CORE_Array<T,Self>& X) {
        Transform(F,X,(*this));
    }
    template<typename LambdaFct>
    inline void transform(LambdaFct&& F,
                          const CORE_Array<T,Self>& X,
                          const CORE_Array<T,Self>& Y) {
        Transform(F,X,Y,*this);
    }
 
    
    template<typename LambdaFct>
    static inline void Transform(LambdaFct&& F,
                                 Self& X) {
        
        //size of all the array
        tIndex n=X.getSize();
       
      
        
        //get the cont begin iterator of X
        auto iX=X.cbegin();
        
         
        //get the begin iterator of R
        auto iR=X.begin();
        
      
        
        //init the random data
        OMP_PARALLEL_SHARED((n,iX,iR,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,
                           iR+start,F);
            
            
        }//end parallel section
    }
     
   
    
    template<typename LambdaFct>
    static inline void Transform(LambdaFct&& F,
                                 const CORE_Array<T,Self>& X,
                                 Self& R) {
        
        //size of all the array
        tIndex n=X.getSize();
       
        //resize the result
        R.setSize(n);
 
       
        
        
        //get the cont begin iterator of X
        auto iX=X.cbegin();
        
         
        //get the begin iterator of R
        auto iR=R.begin();
        
      
        
        //init the random data
        OMP_PARALLEL_SHARED((n,iX,iR,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,
                           iR+start,F);
            
            
        }//end parallel section
    }
       
   
    template<typename LambdaFct>
    static inline void Transform(LambdaFct&& F,
                                 const CORE_Array<T,Self>& X,
                                 const CORE_Array<T,Self>& Y,
                                 Self& R) {
 
        //siez of all the array
        tIndex n=X.getSize();
        
        if (n!=Y.getSize())
            throw CORE_Exception("opemp",
                                 "OMP_ValArray::Transform(F,X,Y,R)",
                                 "X ("+std::to_string(X.getSize())+") & Y ("+std::to_string(Y.getSize())+") have not same size");
        
        R.setSize(n);
 
       
         
        //get the const begin iterator of X
        auto iX=X.cbegin();
        
        //get the const begin iterator of Y
        auto iY=Y.cbegin();
        
         
        //get the begin iterator of R
        auto iR=R.begin();
      
        
       
      
        
        //init the random data
        OMP_PARALLEL_SHARED((n,iX,iY,iR,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,
                           iY+start,
                           iR+start,F);
            
            
        }//end parallel section
    }
 
    inline T norm2() const {
        tIndex n=this->getSize();
        auto iV=this->cbegin();
        
        T norm2=0;
        OMP_PARALLEL_SHARED_REDUCTION((n,iV),(+:norm2)) {//begin parallel section
            // thread id
            tInteger threadId=OMP_GET_THREAD_ID();
            
            // threads number
            tInteger nThreads=OMP_GET_THREADS_NUMBER();
            
            //start index
            tIndex start=threadId*n/nThreads;
            
            //end index
            tIndex end=(threadId+1)*n/nThreads;
            
            
            norm2=0;
            std::for_each(iV+start,iV+end,[&norm2](const auto &Vi){norm2+=Vi*Vi;});
            
            
            
        }//end parallel section
        
        return norm2;
    }
    
    inline void normalize() {
        auto iV=this->begin();
        T inorm=norm2();
 
        if (inorm>std::numeric_limits<T>::epsilon()) {
            inorm=1./inorm;
            tIndex n=this->getSize();
            auto iV=this->cbegin();
            OMP_PARALLEL_SHARED((n,inorm,iV)) {//begin parallel section
                // thread id
                tInteger threadId=OMP_GET_THREAD_ID();
                
                // threads number
                tInteger nThreads=OMP_GET_THREADS_NUMBER();
                
                //start index
                tIndex start=threadId*n/nThreads;
                
                //end index
                tIndex end=(threadId+1)*n/nThreads;
                
                std::transform(iV+start,iV+end,iV+start,
                               iV+start,[&inorm](const auto& Vi){return Vi*inorm;});
                
                
                
            }//end parallel section
            
        }
    }
    
    template<typename Q,class I>
    inline void axpy(const Q& alpha, const CORE_Array<Q,I>& X,const T& beta) {
        tIndex p=functions_numeric::min(X.getSize(),this->getSize());
        
        auto iX=X.cbegin();
        auto iY=this->begin();
        OMP_PARALLEL_SHARED((p,alpha,beta,iX,iY)) {//begin parallel section
               
            // thread id
            tInteger threadId=OMP_GET_THREAD_ID();
            // threads number
            tInteger nThreads=OMP_GET_THREADS_NUMBER();
             
            //start index of the loop
            tIndex start=threadId*p/nThreads;
             
            //end index of the loop
            tIndex end=(threadId+1)*p/nThreads;
 
            
            std::transform(iX+start,iX+end,
                           iY+start,
                           iY+start,[&alpha,&beta](const auto& Xi,const auto& Yi) {return beta*Yi+alpha*Xi;});
        }
    }
    
    //Data consistency
    //=================
    inline  tBoolean isNANContained() const {
        tBoolean isnan=false;//0
        auto iV=this->cbegin();
        tIndex n=this->getSize();
        OMP_PARALLEL_SHARED_REDUCTION((n,iV),(max:isnan)) {//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;
 
            isnan=(std::find_if(iV+start,iV+end,[](const auto& Vi) {return std::isnan(Vi);})!=iV+end);
            
            
        }//end parallel section
        
        return isnan;
    }
   
   
 
    
    
    //algebric method
    //===================
 
   
      
    template<typename Q,class I>
    inline T& scalarProduct(const CORE_Array<Q,I>& X,T& s)  const  {
        tIndex p=functions_numeric::min(this->getSize(),X.getSize());
        auto iV=this->cbegin();
        auto iX=X.cbegin();
        s=0;
        OMP_PARALLEL_SHARED_REDUCTION((p,iX,iV),(+:s)) {//begin parallel section
            
            // thread id
            tInteger threadId=OMP_GET_THREAD_ID();
            
            // threads number
            tInteger nThreads=OMP_GET_THREADS_NUMBER();
            
            //start index
            tInteger start=threadId*p/nThreads;
            
            //end index
            tInteger end=(threadId+1)*p/nThreads;
 
            //start iterator on X
            auto iXi=iX+start;
            
            s=0;
            std::for_each(iV+start,iV+end,[&s,&iXi](const auto& Vi) {s+=Vi*(*iXi);iXi++;});
            
            
        }//end parallel section
        return s;
    }
    
    inline void sum(T& s) const {
        auto iV=this->cbegin();
        tIndex n=this->getSize();
        s=0;
        OMP_PARALLEL_SHARED_REDUCTION((n,iV),(+:s)) {//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;
 
           
            s=0;
            std::for_each(iV+start,iV+end,[&s](const auto& Vi) {s+=Vi;});
            
            
        }//end parallel section
        
    }
 
    
    inline void prod(T& p) const {
        auto iV=this->cbegin();
        tIndex n=this->getSize();
        p=1;
        OMP_PARALLEL_SHARED_REDUCTION((n,iV),(*:p)) {//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;
 
           
            p=1;
            std::for_each(iV+start,iV+end,[&p](const auto& Vi) {p*=Vi;});
            
            
        }//end parallel section
        
        
    }
 
     
    inline tReal l2Norm2() const {
        return L2Distance2(0,null,this->getArray());
    }
    
    template<typename Q,class I>
    inline tReal l2Distance2(const CORE_Array<Q,I>& X)  const  {
        return L2Distance2(X.getSize(),X.getValues(),this->getArray());
    }
 
     
    
    template<typename Q>
    inline static tReal  L2Distance2(const tIndex& nX,const Q* Xs,
                                     const std::valarray<T>& Y) {
    
        
        //d2=\sum_i |Xi-Yi|^2
        tReal d2=0;
 
        //size of common loop
        tIndex n=functions_numeric::min(nX,Y.size());
      
        auto Ys=std::cbegin(Y);
        
        OMP_PARALLEL_SHARED_REDUCTION((n,Xs,Ys),(+:d2)) {
 
            //third id
            tInteger threadId=OMP_GET_THREAD_ID();
            // number of threads 
            tInteger nThreads=OMP_GET_THREADS_NUMBER();
            
            //start index
            tIndex start=threadId*n/nThreads;
            //end index
            tIndex end=(threadId+1)*n/nThreads;
 
            
            
            const Q *iX=Xs+start;
            
            //temporary value
            tReal tmp;
 
            d2=0;
            std::for_each(Ys+start,Ys+end,
                          [&d2,&iX,&tmp](const auto& Yi) {
                              tmp=(*iX)-Yi;
                              d2+=tmp*tmp;
                              iX++;
                          }//end loop on X & Y
                );
        }//end parallel section
        
        tIndex p=nX-n;
        if (p>0) {
            OMP_PARALLEL_SHARED_REDUCTION((p,n,Xs),(+:d2)) {
                
                //third id
                tInteger threadId=OMP_GET_THREAD_ID();
                // number of threads 
                tInteger nThreads=OMP_GET_THREADS_NUMBER();
                
                //start index
                tIndex start=threadId*p/nThreads;
                //end index
                tIndex end=(threadId+1)*p/nThreads;
                
                //iterator on start
                const Q *iX=Xs+n+start;
                //iterator on end
                const Q *iXe=Xs+n+end;
               
                while (iX!=iXe) {
                    d2+=(*iX)*(*iX);
                    iX++;
                }
            }//end parallel section
        }//and loop on X
        
        
          
        p=Y.size()-n;
        if (p>0) {
            OMP_PARALLEL_SHARED_REDUCTION((p,n,Ys),(+:d2)) {
                
                //third id
                tInteger threadId=OMP_GET_THREAD_ID();
                // number of threads 
                tInteger nThreads=OMP_GET_THREADS_NUMBER();
                
                //start index
                tIndex start=threadId*p/nThreads;
                //end index
                tIndex end=(threadId+1)*p/nThreads;
                
                std::for_each(Ys+n+start,Ys+n+end,
                              [&d2](const auto& Yi) {d2+=Yi*Yi;});
               
            }//end parallel section
        }//and loop on Y
        
        return d2;
    }
 
     
    inline tReal linfNorm(tIndex& imax) const {
        return LinfDistance(0,(T*)null,this->getArray(),imax);
    }
 
    template<typename Q,class I>
    inline tReal linfDistance(const CORE_Array<Q,I>& X,tIndex& imax)  const  {
        return LinfDistance(X.getSize(),X.getValues(),this->getArray(),imax);
    }
    
    
    
    template<typename Q>
    inline static tReal LinfDistance(const tIndex& nX,const Q* Xs,
                                     const std::valarray<T>& Y,
                                     tIndex& imax)  {
        
        
 
        //common size
        tIndex n=functions_numeric::min(nX,Y.size());
        
        
        std::pair<tIndex,T> S(0,0);
 
        auto Ys=std::cbegin(Y);
        
        //common loop
        OMP_PARALLEL_SHARED_REDUCTION((n,Ys,Xs), (argrmax:S)) {
            // thread id
            tInteger threadId=OMP_GET_THREAD_ID();
            // threads number
            tInteger nThreads=OMP_GET_THREADS_NUMBER();
            
            //start index
            tIndex start=threadId*n/nThreads;
            //end index
            tIndex  end=(threadId+1)*n/nThreads;
            
           
            
            //begin iterator on start
            const Q *iX=Xs+start;
            
            //start index 
            tIndex i=start;
            //local variable
            tReal norm2;
            std::for_each(Ys+start,Ys+end,
                          [&i,&norm2,&iX,&S](const auto& Yi) {
                              norm2=Yi-(*iX);
                              norm2*=norm2;
                              if (S.second<norm2) {
                                  S.second=norm2;
                                  S.first=i;
                              }
                              iX++;
                              i++;
                          });
        }//end parallel
        
        tIndex p=nX-n;
        if (p>0) {
            OMP_PARALLEL_SHARED_REDUCTION((n,p,Xs), (argrmax:S)) {
                // thread id
                tInteger threadId=OMP_GET_THREAD_ID();
                // threads number
                tInteger nThreads=OMP_GET_THREADS_NUMBER();
                
                //start index
                tIndex start=threadId*p/nThreads;
                //end index
                tIndex  end=(threadId+1)*p/nThreads;
              
                
                //begin iterator on start
                const Q*iX=Xs+n+start;  
                //begin iterator on end
                const Q *iXe=Xs+n+end;
                
                //start index 
                tIndex i=n+start;
                //local variable
                tReal norm2;
                while (iX!=iXe) {
                    norm2=(*iX)*(*iX);
                    if (S.second<norm2) {
                        S.second=norm2;
                        S.first=i;
                    }
                    iX++;
                    i++;
                }
            }//end parallel
            
        }//end loop on X
        p=Y.size()-n;
        if (p>0) {
            OMP_PARALLEL_SHARED_REDUCTION((n,p,Ys), (argrmax:S)) {
                // thread id
                tInteger threadId=OMP_GET_THREAD_ID();
                // threads number
                tInteger nThreads=OMP_GET_THREADS_NUMBER();
                
                //start index
                tIndex start=threadId*p/nThreads;
                //end index
                tIndex  end=(threadId+1)*p/nThreads;
 
                //start index
                tIndex i=n+start;
                //module of Y
                tReal mod2;
                std::for_each(Ys+n+start,Ys+n+end,
                              [&i,&S,&mod2](const auto& Yi) {
                                  mod2=Yi*Yi;
                                  i++;
                                  if (S.second<mod2) {
                                      S.second=mod2;
                                      S.first=i;
                                  }});
                
            }//end parallel
            
        }//end loop on Y
 
        
        imax=S.first;
        return sqrt(S.second);
        
    }
    
  
 
  
 
 
    //ordered methods
    //===============
 
    inline void min(T& minV) const  requires functions_type::isOrderedType<T> {
        auto iV=this->cbegin();
        tIndex n=this->getSize();
        minV=0;
        if (n==0) return;
        minV=(*iV);
        OMP_PARALLEL_SHARED_REDUCTION((n,iV),(min:minV)) {//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;
 
            if (end>start) {
                auto iVs=iV+start;
                minV=*(iVs);
                std::for_each(iVs,iV+end,[&minV](const auto& Vi) {minV=(minV<Vi)?minV:Vi;});
            }
            
        }//end parallel section
        
    }
    
 
    inline void max(T& maxV) const requires functions_type::isOrderedType<T> {
        auto iV=this->cbegin();
        tIndex n=this->getSize();
        maxV=0;
        if (n==0) return;
        maxV=(*iV);
        OMP_PARALLEL_SHARED_REDUCTION((n,iV),
                                      (max:maxV)) {//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;
 
            if (end>start) {
                auto iVs=iV+start;
                maxV=*(iVs);
                std::for_each(iVs,iV+end,[&maxV](const auto& Vi) {maxV=(maxV>Vi)?maxV:Vi;});
            }
            
        }//end parallel section
        
    }
   
 
    
    inline void directionalSort(const tUChar& order) requires functions_type::isOrderedType<T>{
        switch(order) {
        case 'i':
            std::sort(this->begin(),this->end(),std::less<T>());
            break;
        case 'd':
            std::sort(this->begin(),this->end(),std::greater<T>());
            break;
        }
        
    }
    
};
 
 
 
typedef  OMP_ValArray<tReal> OMP_RealValArray;
 
typedef  OMP_ValArray<tFloat> OMP_FloatValArray;
typedef  OMP_ValArray<tDouble> OMP_DoubleValArray;
typedef  OMP_ValArray<tLDouble> OMP_LDoubleValArray;
 
 
typedef  OMP_ValArray<tIndex> OMP_IndexValArray;
 
typedef  OMP_ValArray<tRelativeInteger> OMP_RelativeIntegerValArray;
typedef  OMP_ValArray<tInteger> OMP_IntegerValArray;
 
typedef  OMP_ValArray<tBoolean> OMP_BooleanValArray;
 
 
#endif