OMP_PtrArray.h

#ifndef OMP_PtrArray_H
#define OMP_PtrArray_H
 
//base class header
#include "CORE_Object.h"
 
//exception header
#include "CORE_Exception.h"
 
//temmplate base class for array  header
#include "CORE_PtrArray.h"
 
 
//openmp header
#include "openMP.h"
 
template<class T>
class OMP_PtrArray : public CORE_PtrArray<T,OMP_PtrArray<T>> {
 
private:
    //This class type
    typedef OMP_PtrArray<T>  Self;
    
public:
 
    
    
    OMP_PtrArray() : CORE_PtrArray<T,Self>() {
    }
 
 
    explicit OMP_PtrArray(const tIndex& n)  {
        this->setSize(n);
    }
 
    
    
   
 
    // DESTRUCTORS
    virtual ~OMP_PtrArray() {
    }
 
 
 
     
    //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=(std::initializer_list<T>&& values)  {
        copy(values);
        return *this;
    }
    inline Self& operator=(const std::initializer_list<T>& values)  {
        copy(values);
        return *this;
    }
    
    template<size_t N>
    inline  Self& operator=(const std::array<T,N>& values)  {
        copy(values);
        return *this;
        
    }
    
    inline  Self& operator=(const std::valarray<T>& values)  {
        copy(values);
        return *this;
    }
    
    inline  Self& operator=(const std::vector<T>& 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 I1>
    inline  Self& operator=(const CORE_Array<Q,I1>& values)  {
        copy(values);
        return *this;
    }
    template<typename Q,class I1>
    inline  Self& operator=(const CORE_Array<Q,I1>&& values)  {
        copy(values);
        return *this;
    }
    
   
    
    
    //copy methods
    //============
public:
       
      
    template<typename Q,class I>
    inline void copy(const CORE_Array<Q,I>& cpy) {
        this->setSize(cpy.getSize());
        Copy(cpy.getSize(),cpy.getValues(),this->getValues());
    }
   
    template<typename Q,class I>
    inline void copy(CORE_Array<Q,I>&& cpy) {
        this->setSize(cpy.getSize());
        Copy(cpy.getSize(),cpy.getValues(),this->getValues());
    }
   
    template<class I>
    inline void copy(const Self& cpy) {
        this->setSize(cpy.getSize());
        Copy(cpy.getSize(),cpy.getValues(),this->getValues());
    }
   
    template<class I>
    inline void copy(Self&& cpy) {
        this->setSize(cpy.getSize());
        Copy(cpy.getSize(),cpy.getValues(),this->getValues());
    }
 
  
    template<typename Q>
    inline void copy(const tIndex& n , const Q* Vs) {
        this->setSize(n);
        Copy(n,Vs,this->getValues());
    }
    
    inline void copy(const tIndex& n,const std::initializer_list<T>& Vs) {
        tIndex p=functions_numeric::min(n,Vs.size());
        this->setSize(p);
        Copy(p,Vs,this->getValues());
    }
    inline void copy(const tIndex& n,std::initializer_list<T>&& Vs) {
        tIndex p=functions_numeric::min(n,Vs.size());
        this->setSize(p);
        Copy(p,Vs,this->getValues());
    }
    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);
        this->setSize(p);
        Copy(p,Vs.getData(),this->getValues());
    }
    
    template<typename Q>
    inline void copy(const tIndex& n,const std::valarray<Q>& Vs) {
        tIndex p=functions_numeric::min(n,Vs.size());
        this->setSize(p); 
        Copy(p,&Vs[0],this->getValues());
    }
    template<typename Q>
    inline void copy(const tIndex&n , const std::vector<Q>& Vs) {
        tIndex p=functions_numeric::min(n,Vs.size());
        this->setSize(p);
        Copy(Vs.size(),Vs,this->getValues());
    }
   
 
 
    //static Copy
    //----------
   
    template<class Q>
    inline static  void Copy(const tIndex& n,const Q* X,T* R) {
        
     
        OMP_PARALLEL_SHARED((n,X,R)) {//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
            T* iR=R;
            iR+=start;
           
            //iterator on cpy values at start
            const Q* iX=X;
            iX+=start;
 
            if (sizeof(T)==sizeof(Q)) {
                memcpy(iR,iX,sizeof(T)*(end-start));
            } else {
                //iterator on R values at end
                const T* iRe=R;
                iRe+=end;
                
                //copy the values of cpyA into aValue for element in [start,end[
                while (iR!=iRe) {
                    (*iR)=(*iX);
                    iR++;
                    iX++;
                }
            }
            
        }//end parallel section
    }
    
    template<typename Q>
    inline static  void Copy(const tIndex& n,const std::vector<Q>& X,T *R) {
     
   
        //values of X
        auto vX=X.cbegin();
        
        OMP_PARALLEL_SHARED((n,vX,R)) {//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
            T* iR=R;
            iR+=start;
            
            //iterator on R values at end
            const T* iRe=R;
            iRe+=end;
            
            //iterator on cpy values at start
            auto iX=vX;
            iX+=start;
 
            
            while (iR!=iRe) {
                (*iR)=(*iX);
                iR++;
                iX++;
            }
            
        }//end parallel section
    }
 
    
   
    inline static  void Copy(const tIndex& n,const std::initializer_list<T>& X,T *R) {
     
   
        //get values of X
        auto vX=X.begin();
        
        OMP_PARALLEL_SHARED((n,vX,R)) {//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
            T* iR=R;
            iR+=start;
            
            //iterator on R values at end
            const T* iRe=R;
            iRe+=end;
            
            //iterator on X at start
            auto iX=vX;
            iX+=start;
 
            
            while (iR!=iRe) {
                (*iR)=(*iX);
                iR++;
                iX++;
            }
            
        }//end parallel section
    }
 
   
    inline static  void Copy(const tIndex& n,std::initializer_list<T>&& X,T *R) {
     
   
        //get values of X
        auto vX=X.begin();
        
        OMP_PARALLEL_SHARED((n,vX,R)) {//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
            T* iR=R;
            iR+=start;
            
            //iterator on R values at end
            const T* iRe=R;
            iRe+=end;
            
            //iterator on X at start
            auto iX=vX;
            iX+=start;
 
            
            while (iR!=iRe) {
                (*iR)=(*iX);
                iR++;
                iX++;
            }
            
        }//end parallel sections
    }
 
   
   
    //initializers
    //============
    inline void initialize(const T& x) {
        Initialize(x,this->getSize(),this->getValues());
    }
    
    inline static  void Initialize(const T& v, const tIndex& n,T* R) {
        
        if (fabs(v)<std::numeric_limits<T>::epsilon()) {
 
            //set to 0 by memset method
            
            OMP_PARALLEL_SHARED((n,R)) {//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;
 
                if (end>start) memset(&R[start],0,(end-start)*sizeof(T));
                
            }
        
        } else {
            
            OMP_PARALLEL_SHARED((n,v,R)) {//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
                T* iR=R;
                iR+=start;
                
                //iterator on R values at end
                const T* iRe=R;
                iRe+=end;
                
                
                
                
                //copy the values of cpyA into aValue for element in [start,end[
                while (iR!=iRe) {
                    (*iR)=v;
                    iR++;
                }
                
            }//end parallel section
        }
    }
    
    
    inline static void UniformRandomize(const T& min,const T& max,
                                        const tIndex& n,T* values) {
 
        //x=d*(max-min)+min
        T alpha=(max-min)/RAND_MAX;
 
      
        //init the random data
        OMP_PARALLEL_SHARED((alpha,n,min,values)) {//begin parallel section
            
            // threads number
            tInteger nThreads=OMP_GET_THREADS_NUMBER();
            
            // thread id
            tInteger threadId=OMP_GET_THREAD_ID();
 
             
            //start index
            tIndex start=threadId*n/nThreads;
            
            //end index
            tIndex end=(threadId+1)*n/nThreads;
 
            //iterator on start values
            T* iV=values;iV+=start;
            //iterator on end values
            const T* iVe=values;iVe+=end;
 
            //set random seed by process
            tULLInt seed;
            for(start=0;start<=threadId;start++) {
                seed=std::rand();
            }
            std::srand(seed);
          
            
            //copy the values of vs  into A for element in [start,end[
            while (iV!=iVe) {
                (*iV)=alpha*std::rand()+min;
                iV++;
            }
            
        }//end parallel section
    }
   
 
    inline void uniformRandomize(const T& min,const T& max) {
        UniformRandomize(min,max,this->getSize(),this->getValues());
    }
 
   
 
    //compount assignement operator
    //==============================
    
    inline Self& operator+=(const T& v) {
        Add(v,this->getSize(),this->getValues());
        return (*this);
    }
    inline Self& operator-=(const T& v) {
        Add(-v,this->getSize(),this->getValues());
        return (*this);
    }
    
    inline Self& operator*=(const T& v) {
        Multiply(v,this->getSize(),this->getValues());
        return (*this);
       
    }
    inline Self& operator/=(const T& v) {
        if (fabs(v)<std::numeric_limits<T>::epsilon()) {
            throw CORE_Exception("core",
                                 "OMP_PtrArray::/=("+std::to_string(v)+")",
                                 "division by zero");
        }
        Multiply(1./v,this->getSize(),this->getValues());
        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 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 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 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 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 Self& 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>& v) {
        tIndex p=functions_numeric::min(v.getSize(),this->getSize());
        Add(p,v.getValues(),this->getValues());
        return *this;
    }
 
   
    template<typename Q,class I>
    inline Self& operator-=(const CORE_Array<Q,I>& v) {
        tIndex p=functions_numeric::min(v.getSize(),this->getSize());
        Sub(p,v.getValues(),this->getValues());
        return (*this);
    }
    template<typename Q,class I>
    inline Self& operator*=(const CORE_Array<Q,I>& v) {
        tIndex p=functions_numeric::min(v.getSize(),this->getSize());
        Multiply(p,v.getValues(),this->getValues());
        return (*this);
    }
    template<typename Q,class I>
    inline Self& operator/=(const CORE_Array<Q,I>& v) {
        tIndex p=functions_numeric::min(v.getSize(),this->getSize());
        Divide(p,v.getValues(),this->getValues());
        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_PtrArray::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
    }
    
    template<typename Q>
    static void Add(const Q& v,const tIndex& n,T* values) {
       
        OMP_PARALLEL_SHARED((n,values,v)) {//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 start index
            T* iV=values;iV+=start;
 
            //iterator on end index
            T* iVe=values;iVe+=end;
          
            //loop on index
            while(iV!=iVe) {
                (*iV)+=v;
                iV++;
            }
            
        }//end parallel section
    }
 
    template<typename Q>
    inline static void Add(const tIndex & n,const Q* Y,T* X) {
        
        OMP_PARALLEL_SHARED((n,X,Y)) {//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 start index of X
            T* iX=X;iX+=start;
            
            //iterator on start index of Y
            const Q* iY=Y;iY+=start;
            
            //iterator on end index of X
            T* iXe=X;iXe+=end;
          
            //loop on index
            while(iX!=iXe) {
                (*iX)+=(*iY);
                iY++;
                iX++;
            }
            
        }//end parallel section
        
    }
    
    template<typename Q>
    inline static void Sub(const tIndex & n,const Q* Y,T* X) {
        
        OMP_PARALLEL_SHARED((n,X,Y)) {//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 start index of X
            T* iX=X;iX+=start;
            
            //iterator on start index of Y
            const Q* iY=Y;iY+=start;
            
            //iterator on end index of X
            T* iXe=X;iXe+=end;
          
            //loop on index
            while(iX!=iXe) {
                (*iX)-=(*iY);
                iY++;
                iX++;
            }
            
        }//end parallel section
        
    }
 
    
 
    template<typename Q>
    static void Multiply(const Q& v,const tIndex& n,T* values) {
        
        OMP_PARALLEL_SHARED((n,values,v)) {//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 start index
            T* iV=values;iV+=start;
            
            //iterator on end index
            T* iVe=values;iVe+=end;
            
            //loop on index
            while(iV!=iVe) {
                (*iV)*=v;
                iV++;
            }
            
        }//end parallel section
    }
 
    
    template<typename Q>
    inline static void Multiply(const tIndex & n,const Q* Y,T* X) {
        
        OMP_PARALLEL_SHARED((n,X,Y)) {//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 start index of X
            T* iX=X;iX+=start;
            
            //iterator on start index of Y
            const Q* iY=Y;iY+=start;
            
            //iterator on end index of X
            T* iXe=X;iXe+=end;
            
            //loop on index
            while(iX!=iXe) {
                (*iX)*=(*iY);
                iY++;
                iX++;
            }
            
        }//end parallel section
        
    }
    template<typename Q>
    inline static void Divide(const tIndex & n,const Q* Y,T* X) {
 
        Q eps=std::numeric_limits<Q>::epsilon();
        
        OMP_PARALLEL_SHARED((eps,n,X,Y)) {//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 start index of X
            T* iX=X;iX+=start;
            
            //iterator on start index of Y
            const Q* iY=Y;iY+=start;
            
            //iterator on end index of X
            T* iXe=X;iXe+=end;
          
            //loop on index
            while(iX!=iXe) {
                if (fabs(*iY)<eps) throw CORE_Exception("openmp",
                                                        "OMP_PtrArray::Divide()",
                                                        "division by zero");
                (*iX)/=(*iY);
                iY++;
                iX++;
            }
            
        }//end parallel section
        
    }
 
    inline void normalize() {
        Normalize(this->getSize(),this->getValues());
    }
 
    inline static T norm2(const tIndex& n,const T *values)  {
        //initialization of norm
        T norm2=0;
        
        OMP_PARALLEL_SHARED_REDUCTION((n,values),(+:norm2)) {//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;;
            
            const T* iV =values;iV+=start;//iterator on start index
 
            const T* eiV=values;eiV+=end;//end iterator on end index
            
            norm2=0;
            while (iV!=eiV) {
                norm2+=(*iV)*(*iV);
                iV++;
            }
                
        }
        return norm2;
    }
    inline static void Normalize(const tIndex& n,T* values) {
       
        //inverse of norm of V
        T inorm=norm2(n,values);
        if (fabs(inorm)>std::numeric_limits<T>::epsilon()) {//norm is not null
            //inverse of norm
            inorm=1./sqrt(inorm);
            //loop on 
            OMP_PARALLEL_SHARED((n,values,inorm)) {//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;;
                T* iV =values;iV+=start;//iterator on values
                const T* iVe=values;iVe+=end;//iterator on end values
                while (iV!=iVe) {
                    (*iV)*=inorm;
                    iV++;
                }
               
            }    
            
        }
    }
    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());
        AXPY(p,alpha,X.getValues(),beta,this->getValues());
    }
 
    template<typename Q>
    inline void AXPY(const tIndex& n,const Q& alpha, const Q* X,const T& beta, T* Y) {
       
        OMP_PARALLEL_SHARED((n,alpha,beta,X,Y)) {//begin parallel section
               
             // thread id
             tInteger threadId=OMP_GET_THREAD_ID();
             
             // threads number
             tInteger nThreads=OMP_GET_THREADS_NUMBER();
             
             //start index of the values within the thread
             tIndex start=threadId*n/nThreads;
             
             //end index of the values within the thread
             tIndex end=(threadId+1)*n/nThreads;
 
             //start iterator on X
             const Q *iX=X+start;
             
             //start iterator on Y
             T *iY=Y+start;
             
             //end iterator on Y
             const T  *iYe=Y+end;
             
             while (iY!=iYe) {
                 (*iY)=beta*(*iY)+alpha*(*iX);
                 iX++;
                 iY++;
             }
         }
    }
    
    
    //Data consistency
    //=================
    inline  tBoolean isNANContained() const {
        return IsNANContained(this->getSize(),this->getValues());
    }
    inline  static tBoolean IsNANContained(const tIndex& n,const T* values)  {
        
        
        tBoolean isnan=false;//0
        OMP_PARALLEL_SHARED_REDUCTION((n,values),(max:isnan)) {//begin parallel section
            
            // thread id
            tInteger threadId=OMP_GET_THREAD_ID();
            
            // threads number
            tInteger nThreads=OMP_GET_THREADS_NUMBER();
            
            //start index
            const T* iV=values+threadId*n/nThreads;
            
            //end index
            const T* iVe=values+(threadId+1)*n/nThreads;
            
            
            //test if there is a nan value within  [start,end[
            while (iV!=iVe) {
                if (std::isnan(*iV)) break;
                iV++;
            }
            
            isnan=(iV!=iVe);
            
        }//end parallel section
        
        return isnan;
        
    }
    
   
   
 
    
    
    //algebric method
    //===================
 
    inline void  sum(T& s) const {
        s=Sum(this->getSize(),this->getValues());
    }
 
     
    inline static T Sum(const tIndex& n,const T* values)  {
        
       
       
        T s=0;//0
        OMP_PARALLEL_SHARED_REDUCTION((n,values),(+:s)) {//begin parallel section
            
            // thread id
            tInteger threadId=OMP_GET_THREAD_ID();
            // threads number
            tInteger nThreads=OMP_GET_THREADS_NUMBER();
      
            //start index
            const T* iV=values+threadId*n/nThreads;
            
            //end index
            const T* iVe=values+(threadId+1)*n/nThreads;
            
            
            //test if there is a nan value within  [start,end[
            s=0;
            while (iV!=iVe){s+=(*iV);iV++;}
                         
            
        }//end parallel section
        
        return s;
    }
    inline void prod(T& p) const {
        p=Prod(this->getSize(),this->getValues());
    }
      
    inline static T Prod(const tIndex& n,const T* values)  {
 
        T p=1;
        OMP_PARALLEL_SHARED_REDUCTION((n,values),(*:p)) {//begin parallel section
            
            // thread id
            tInteger threadId=OMP_GET_THREAD_ID();
            // threads number
            tInteger nThreads=OMP_GET_THREADS_NUMBER();
      
            //start index
            const T* iV=values+threadId*n/nThreads;
            
            //end index
            const T* iVe=values+(threadId+1)*n/nThreads;
            
            
            //test if there is a nan value within  [start,end[
            p=1;
            while (iV!=iVe){p*=(*iV);iV++;}
                         
            
        }//end parallel section
        
        return p;
       
                 
    }
    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());
        s=ScalarProduct(p,X.getValues(),this->getValues());
        return s;
    }
 
    
 
    template<typename Q>
    inline static T ScalarProduct(const tIndex& n,
                                  const Q* X,
                                  const T* Y) {
 
        //s=sum_i  Xi Yi
        T s=0;
        
        OMP_PARALLEL_SHARED_REDUCTION((n,Y,X),(+: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 of the loop
            auto iX=X;iX+=start;
            auto iY=Y;iY+=start;
            
            
            //end iterator of the loop
            auto iXe=X;iXe+=end;
            s=0;
            while (iX!=iXe) {  
                s+=(*iX)*(*iY);
                iY++;
                iX++;
            }
        }//end parallel loop
        
        return s;
    }
      
  
    
    inline tReal l2Norm2() const {
        return L2Distance2(0,(tReal*)null,this->getSize(),this->getValues());
    }
 
    template<typename Q,class I>
    inline tReal l2Distance2(const CORE_Array<Q,I>& X)  const  {
        return L2Distance2(X.getSize(),X.getValues(),this->getSize(),this->getValues());
    }
 
     
    
    template<typename Q>
    inline static tReal  L2Distance2(const tIndex& nX,const Q* Xs,
                                     const tIndex& nY,const T* Ys) {
    
        
        //d2=\sum_i |Xi-Yi|^2
        tReal d2=0;
 
        //size of common loop
        tIndex n=functions_numeric::min(nX,nY);
    
       
        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;
 
            //iterator on start
            const T *iY=Ys+start;
            const Q *iX=Xs+start;
            //iterator on end
            const Q  *iXe=Xs+end;
            //temporary value
            tReal tmp;
 
            d2=0;
            while (iX!=iXe) {
                tmp=(*iX)-(*iY);
                d2+=tmp*tmp;
                iX++;
                iY++;
            }//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=nY-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;
                
                //iterator on start
                const T *iY=Ys+n+start;
                //iterator on end
                const T *iYe=Ys+n+end;
                while (iY!=iYe) {
                    d2+=(*iY)*(*iY);
                    iY++;
                }
               
            }//end parallel section
        }//and loop on Y
 
        return d2;
    }
 
     
    inline tReal linfNorm(tIndex& imax) const {
        return LinfDistance(0,(tReal*)null,this->getSize(),this->getValues(),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->getSize(),this->getValues(),imax);
    }
    
    
    
    template<typename Q>
    inline static tReal LinfDistance(const tIndex& nX,const Q* Xs,
                                     const tIndex& nY,const T* Ys,
                                     tIndex& imax)  {
        
       
        //common size
        tIndex n=functions_numeric::min(nX,nY);
        
 
        //S=(index,value)
        std::pair<tIndex,tReal> S(0,0);
 
        //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 T *iY=Ys+start;
            
            //begin iterator on start
            const Q *iX=Xs+start;
            
            //start index 
            tIndex i;
            //local variable
            T norm2;
            for (i=start;i<end;i++) {
                norm2=(*iY)-(*iX);
                norm2*=norm2;
                if (S.second<norm2) {
                    S.second=norm2;
                    S.first=i;
                }
                iY++;
                iX++;
                
            }
        }//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=nY-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;
              
                
                //begin iterator on start
                const T *iY=Ys+n+start;  
                //begin iterator on end
                const T *iYe=Ys+n+end;
                
                //start index 
                tIndex i=n+start;
                //local variable
                tReal norm2;
                while (iY!=iYe) {
                    norm2=(*iY)*(*iY);
                    if (S.second<norm2) {
                        S.second=norm2;
                        S.first=i;
                    }
                    iY++;
                    i++;
                }
            }//end parallel
            
        }//end loop on X
 
 
        //return value
        imax=S.first;
        return sqrt(S.second);
        
    }
    
  
 
    //ordered methods
    //===============
 
    inline void min(T& m) const  requires functions_type::isOrderedType<T> {
        m=Min(this->getSize(),this->getValues());
    }
    inline static T Min(const tIndex& n,const T* values)  requires functions_type::isOrderedType<T> {
 
 
        //default min value set to the first value
        T p=(n>0)?(*values):0;
        OMP_PARALLEL_SHARED_REDUCTION((n,values),(min:p)) {//begin parallel section
            
            // thread id
            tInteger threadId=OMP_GET_THREAD_ID();
            // threads number
            tInteger nThreads=OMP_GET_THREADS_NUMBER();
      
            //start index
            const T* iV=values+threadId*n/nThreads;
            
            //end index
            const T* iVe=values+(threadId+1)*n/nThreads;
            
            while (iV!=iVe){p=(p<(*iV))?p:(*iV);iV++;}
                         
            
        }//end parallel section
        
        return p;
       
                 
    }
 
    inline void max(T& m) const requires functions_type::isOrderedType<T> {
        m=Max(this->getSize(),this->getValues());
    }
    inline static T Max(const tIndex& n,const T* values)  requires functions_type::isOrderedType<T> {
 
 
        //defaut max values set to the first values 
        T p=(n>0)?(*values):0;
        
        OMP_PARALLEL_SHARED_REDUCTION((n,values),(max:p)) {//begin parallel section
            
            // thread id
            tInteger threadId=OMP_GET_THREAD_ID();
            // threads number
            tInteger nThreads=OMP_GET_THREADS_NUMBER();
      
            //start index
            const T* iV=values+threadId*n/nThreads;
            
            //end index
            const T* iVe=values+(threadId+1)*n/nThreads;
            
            while (iV!=iVe){p=(p>(*iV))?p:(*iV);iV++;}
                         
            
        }//end parallel section
        
        return p;
       
                 
    }
 
    
    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_PtrArray<tReal> OMP_RealPtrArray;
 
typedef  OMP_PtrArray<tFloat> OMP_FloatPtrArray;
typedef  OMP_PtrArray<tDouble> OMP_DoublePtrArray;
typedef  OMP_PtrArray<tLDouble> OMP_LDoublePtrArray;
 
 
typedef  OMP_PtrArray<tIndex> OMP_IndexPtrArray;
 
typedef  OMP_PtrArray<tRelativeInteger> OMP_RelativeIntegerPtrArray;
typedef  OMP_PtrArray<tInteger> OMP_IntegerPtrArray;
 
typedef  OMP_PtrArray<tBoolean> OMP_BooleanPtrArray;
 
//default field
#ifndef DEFAULT_OMP_ARRAY
#define DEFAULT_OMP_ARRAY
typedef OMP_RealPtrArray OMP_RealArray ;
#endif
 
 
#endif