update

include/amscppnarray/amscppnarray.hpp

@@ -15,8 +15,8 @@ namespace narray
 {
 
 // Forward declaration of the operator- template function
-template<typename T> class narray;
+//template<typename T> class narray;
-template<typename T> narray<T> operator-(const narray<T> &n);
+//template<typename T> narray<T> operator-(const narray<T> &n);
 
 typedef int64_t narray_size_t;
 
@@ -48,6 +48,7 @@ template<typename T> class narray
     ~narray();
 
     int resize(const narray_size_t newlength);
+    const narray_size_t size() const;
 
     T& operator[](const narray_size_t ind);
     const T& operator[](const narray_size_t ind) const;
@@ -71,14 +72,35 @@ template<typename T> class narray
     narray<T> operator*(const T val) const;
     narray<T> operator/(const T val) const;
     template<class U> friend narray<U> operator-(const narray<U> &n);
+    narray<T> applyfn(T (*fnptr)(T)) const;
+
+    //min,max
+    T min() const;
+    T max() const;
+
     
     //subarray operations
+    narray<T> subarray(const narray_size_t ind1, const narray_size_t ind2) const;
 
     //permutations and sorting
+
+
+
 };
 
+//Initializers
+template<typename T> narray<T> arange(const T low, const T high, const T incr);
+template<typename T> narray<T> zeros(const narray_size_t N);
+template<typename T> narray<T> ones(const narray_size_t N);
+
+narray<float> linspacef(const float low, const float high, narray_size_t N);
+narray<double> linspace(const double low, const double high, narray_size_t N);
 
 
+/////////
+//Tests//
+/////////
+
 void test_narray1();
 void test_narray2();
 void test_narray3();
@@ -88,5 +110,9 @@ void test_narray3();
 
 #include <amscppnarray/amscppnarray_impl.hpp>
 #include <amscppnarray/amscppnarray_imploper.hpp>
+#include <amscppnarray/amscppnarray_implsubarray.hpp>
+#include <amscppnarray/amscppnarray_implinits.hpp>
+
+#include <amscppnarray/amscppnarray_random.hpp>
 
 #endif

include/amscppnarray/amscppnarray_impl.hpp

@@ -297,6 +297,10 @@ template<typename T> void narray<T>::setall(const T& val)
     return;
 }
 
+template<typename T> const narray_size_t narray<T>::size() const
+{
+    return length;
+}
 
 }; //end namespace narray
 }; //end namespace ams

include/amscppnarray/amscppnarray_implinits.hpp

@@ -0,0 +1,108 @@
+#ifndef __AMSCPPNARRAY_IMPLINITS_HPP__
+#define __AMSCPPNARRAY_IMPLINITS_HPP__
+
+namespace ams
+{
+namespace narray
+{
+
+    template<typename T> void narray_arange_tf(
+        narray<T> *in,
+        const T low, const T high, const T incr,
+        int threadnum, int nthreads
+    )
+    {
+        int I0,I1,Is,I;
+        Is = in->length/nthreads;
+        I0 = Is*threadnum;
+        I1 = (threadnum>=(nthreads-1))? in->length: Is*(threadnum+1);
+        for(I=I0;I<I1;I++)
+        {
+            in->data[I] = low + incr*I;
+        }
+    }
+
+template<typename T> narray<T> arange(const T low, const T high, const T incr)
+{
+    narray<T> ret;
+    narray_size_t I;
+    int J;
+    int nthreads;
+    std::vector<std::thread*> threads;
+
+    narray_size_t N = (narray_size_t)((high-low)/incr);
+    if(N<=0)
+    {
+        N = 0;
+    }
+    ret.resize(N);
+
+    if(ret.length==0)
+    {
+        return ret;
+    }
+
+    if(ret.length<narray_thread_sz)
+    {
+        for(I=0;I<ret.length;I++)
+        {
+            ret.data[I] = low + incr*I;
+        }
+    }
+    else
+    {
+        //threaded operation
+        nthreads = std::thread::hardware_concurrency();
+        nthreads = (nthreads<1) ? 1 : nthreads;
+        nthreads = (nthreads>narray_max_threads) ? narray_max_threads : nthreads;
+        //if(nthreads<1) nthreads=1;
+        //if(nthreads>narray_max_threads) nthreads = narray_max_threads;
+        threads.resize(nthreads);
+        for(J=0;J<nthreads;J++)
+        {
+            threads[J] = new(std::nothrow) std::thread(
+                &narray_arange_tf<T>,
+                &ret,
+                low,high,incr,
+                J,nthreads
+            );
+        }
+        for(J=0;J<nthreads;J++)
+        {
+            if(threads[J]!=NULL)
+            { 
+                threads[J]->join(); 
+                delete threads[J]; 
+                threads[J]= NULL;
+            }
+        }
+    }
+
+
+
+    return ret;
+}
+
+template<typename T> narray<T> zeros(const narray_size_t N)
+{
+    narray<T> ret;
+    ret.resize(N);
+    ret.setall(0);
+
+    return ret;
+}
+
+template<typename T> narray<T> ones(const narray_size_t N)
+{
+    narray<T> ret;
+    ret.resize(N);
+    ret.setall(1);
+
+    return ret;
+}
+
+
+};
+};
+
+#endif

include/amscppnarray/amscppnarray_imploper.hpp

@@ -223,7 +223,7 @@ namespace narray
     
         if(this->length<narray_thread_sz)
         {
-            for(I=0;I<this->length;I++) ret.data[I] = this->data[I] + other.data[I];
+            for(I=0;I<this->length;I++) ret.data[I] = this->data[I] - other.data[I];
         }
         else
         {
@@ -730,6 +730,246 @@ namespace narray
         return ret;
     }
 
+    template<typename T> void narray_applyfn_tf( 
+        const narray<T> *a,
+        T (*fnptr)(T),
+        narray<T> *out,
+        int threadnum, int nthreads
+    )
+    {
+        int I0,I1,Is,I;
+    
+        Is = a->length/nthreads;
+        I0 = Is*threadnum;
+        I1 = (threadnum>=(nthreads-1))? a->length: Is*(threadnum+1);
+        for(I=I0;I<I1;I++)
+        {
+            out->data[I] = fnptr(a->data[I]);
+        }
+    
+        return;
+    }
+
+    template<typename T> narray<T> narray<T>::applyfn(T (*fnptr)(T)) const
+    {
+        narray<T> ret;
+        narray_size_t I;
+        int J;
+        int nthreads;
+        std::vector<std::thread*> threads;
+    
+        ret.resize(this->length);
+    
+        if(this->length<narray_thread_sz)
+        {
+            for(I=0;I<this->length;I++) ret.data[I] = fnptr(this->data[I]);
+        }
+        else
+        {
+            //threaded operation
+            nthreads = std::thread::hardware_concurrency();
+            nthreads = (nthreads<1) ? 1 : nthreads;
+            nthreads = (nthreads>narray_max_threads) ? narray_max_threads : nthreads;
+            //if(nthreads<1) nthreads=1;
+            //if(nthreads>narray_max_threads) nthreads = narray_max_threads;
+            threads.resize(nthreads);
+            for(J=0;J<nthreads;J++)
+            {
+                threads[J] = new(std::nothrow) std::thread(
+                    &narray_applyfn_tf<T>,
+                    this,
+                    fnptr,
+                    &ret,
+                    J,nthreads
+                );
+            }
+            for(J=0;J<nthreads;J++)
+            {
+                if(threads[J]!=NULL)
+                { 
+                    threads[J]->join(); 
+                    delete threads[J]; 
+                    threads[J]= NULL;
+                }
+            }
+        }
+        return ret;
+    }
+
+    template<typename T> void narray_min_tf( 
+        const narray<T> *in,
+        narray<T> *out,
+        int threadnum, int nthreads
+    )
+    {
+        int I0,I1,Is,I;
+    
+        Is = in->length/nthreads;
+        I0 = Is*threadnum;
+        I1 = (threadnum>=(nthreads-1))? in->length: Is*(threadnum+1);
+        for(I=I0;I<I1;I++)
+        {
+            if(in->data[I]<out->data[threadnum])
+                out->data[threadnum] = in->data[I];
+        }
+    
+        return;
+    }
+
+    template<typename T> T narray<T>::min() const
+    {
+        T ret = T();
+        narray<T> minarr;
+        narray_size_t I;
+        int J;
+        int nthreads;
+        std::vector<std::thread*> threads;
+    
+        if(this->length<=0)
+        {
+            return ret;
+        }
+    
+        if(this->length<narray_thread_sz)
+        {
+            ret = this->data[0];
+            for(I=0;I<this->length;I++)
+            {
+                if(this->data[I]<ret)
+                {
+                    ret = this->data[I];
+                }
+            }
+        }
+        else
+        {
+            //threaded operation
+            
+            nthreads = std::thread::hardware_concurrency();
+            nthreads = (nthreads<1) ? 1 : nthreads;
+            nthreads = (nthreads>narray_max_threads) ? narray_max_threads : nthreads;
+            //if(nthreads<1) nthreads=1;
+            //if(nthreads>narray_max_threads) nthreads = narray_max_threads;
+            threads.resize(nthreads);
+            minarr.resize(nthreads);
+            minarr.setall(this->data[0]);
+            for(J=0;J<nthreads;J++)
+            {
+                threads[J] = new(std::nothrow) std::thread(
+                    &narray_min_tf<T>,
+                    this,
+                    &minarr,
+                    J,nthreads
+                );
+            }
+            for(J=0;J<nthreads;J++)
+            {
+                if(threads[J]!=NULL)
+                { 
+                    threads[J]->join(); 
+                    delete threads[J]; 
+                    threads[J]= NULL;
+                }
+            }
+            ret = this->data[0];
+            for(J=0;J<nthreads;J++)
+            {
+                if(minarr[J]<ret)
+                {
+                    ret = minarr[J];
+                }
+            }
+        }
+        return ret;
+    }
+
+    template<typename T> void narray_max_tf( 
+        const narray<T> *in,
+        narray<T> *out,
+        int threadnum, int nthreads
+    )
+    {
+        int I0,I1,Is,I;
+    
+        Is = in->length/nthreads;
+        I0 = Is*threadnum;
+        I1 = (threadnum>=(nthreads-1))? in->length: Is*(threadnum+1);
+        for(I=I0;I<I1;I++)
+        {
+            if(in->data[I]>out->data[threadnum])
+                out->data[threadnum] = in->data[I];
+        }
+    
+        return;
+    }
+
+    template<typename T> T narray<T>::max() const
+    {
+        T ret = T();
+        narray<T> maxarr;
+        narray_size_t I;
+        int J;
+        int nthreads;
+        std::vector<std::thread*> threads;
+    
+        if(this->length<=0)
+        {
+            return ret;
+        }
+    
+        if(this->length<narray_thread_sz)
+        {
+            ret = this->data[0];
+            for(I=0;I<this->length;I++)
+            {
+                if(this->data[I]>ret)
+                {
+                    ret = this->data[I];
+                }
+            }
+        }
+        else
+        {
+            //threaded operation
+            
+            nthreads = std::thread::hardware_concurrency();
+            nthreads = (nthreads<1) ? 1 : nthreads;
+            nthreads = (nthreads>narray_max_threads) ? narray_max_threads : nthreads;
+            //if(nthreads<1) nthreads=1;
+            //if(nthreads>narray_max_threads) nthreads = narray_max_threads;
+            threads.resize(nthreads);
+            maxarr.resize(nthreads);
+            maxarr.setall(this->data[0]);
+            for(J=0;J<nthreads;J++)
+            {
+                threads[J] = new(std::nothrow) std::thread(
+                    &narray_max_tf<T>,
+                    this,
+                    &maxarr,
+                    J,nthreads
+                );
+            }
+            for(J=0;J<nthreads;J++)
+            {
+                if(threads[J]!=NULL)
+                { 
+                    threads[J]->join(); 
+                    delete threads[J]; 
+                    threads[J]= NULL;
+                }
+            }
+            ret = this->data[0];
+            for(J=0;J<nthreads;J++)
+            {
+                if(maxarr[J]>ret)
+                {
+                    ret = maxarr[J];
+                }
+            }
+        }
+        return ret;
+    }
+
 }; //end namespace narray
 }; //end namespace ams
 

include/amscppnarray/amscppnarray_implsubarray.hpp

@@ -0,0 +1,44 @@
+#ifndef __AMSCPPNARRAY_IMPLSUBARRAY_HPP__
+#define __AMSCPPNARRAY_IMPLSUBARRAY_HPP__
+
+namespace ams
+{
+namespace narray
+{
+
+    template<typename T> narray<T> narray<T>::subarray(const narray_size_t ind1, const narray_size_t ind2) const
+    {
+        int res;
+        narray<T> ret;
+        T defval = T();
+        narray_size_t I,J;
+
+        int nl = ind2-ind1;
+
+        if(nl<=0)
+        {
+            ret.resize(0);
+            return ret;
+        }
+
+        ret.resize(nl);
+        for(I=0;I<nl;I++)
+        {
+            J = I + ind1;
+            if(J<0 || J>=this->length)
+            {
+                ret.data[I] = defval;
+            }
+            else
+            {
+                ret.data[I] = this->data[J];
+            }
+        }
+        
+        return ret;
+    }
+
+};
+};
+
+#endif

include/amscppnarray/amscppnarray_random.hpp

@@ -0,0 +1,12 @@
+#ifndef __AMSCPPNARRAY_RANDOM_HPP__
+#define __AMSCPPNARRAY_RANDOM_HPP__
+
+namespace ams
+{
+namespace narray
+{
+
+};
+};
+
+#endif

src/amscppnarray/amscppnarray.cpp

@@ -74,6 +74,7 @@ void test_narray2()
     
 }
 
+
 void test_narray3()
 {
     int I;
@@ -83,17 +84,30 @@ void test_narray3()
     b = {0,0.25,0.5,0.675,0.8675};
     c.resize(5); c.setall(2.0f);
 
-    d = -a*(b+c);
+    d = -a*(b+c*2+2);
     for(I=0;I<5;I++)
     {
         printf("d[%d]=%1.3f\n",I,d[I]);
 
     }
 
-    printf("d==-a*(b+c):%d\n",d==(-a*(b+c)));
+    printf("d==-a*(b+c*2+2):%d\n",d==(-a*(b+c*2+2)));
     d[4] = 0.0f;
-    printf("d==-a*(b+c):%d\n",d==(-a*(b+c)));
+    printf("d==-a*(b+c*2+2):%d\n",d==(-a*(b+c*2+2)));
+    d = d.applyfn(::sinf);
+
+    for(I=0;I<5;I++)
+    {
+        printf("d[%d]=%1.3f\n",I,d[I]);
+
+    }
 
+    d = d.subarray(0,1);
+    for(I=0;I<d.length;I++)
+    {
+        printf("d[%d]=%1.3f\n",I,d[I]);
+
+    }
 }
 
 

src/amscppnarray/amscppnarray_inits.cpp

@@ -0,0 +1,173 @@
+#include <amscppnarray/amscppnarray.hpp>
+
+namespace ams
+{
+namespace narray
+{
+
+void narray_linspacef_tf(
+    narray<float> *in,
+    const float low, const float high, const narray_size_t N,
+    int threadnum, int nthreads
+)
+{
+    int I0,I1,Is,I;
+    float g;
+
+    Is = in->length/nthreads;
+    I0 = Is*threadnum;
+    I1 = (threadnum>=(nthreads-1))? in->length: Is*(threadnum+1);
+    for(I=I0;I<I1;I++)
+    {
+        g = ((float)I)/((float)(N-1));
+        in->data[I] = low*(1.0f-g) + high*g;
+    }
+}
+
+
+narray<float> linspacef(const float low, const float high, narray_size_t N)
+{
+    narray<float> ret;
+    narray_size_t I;
+    int J;
+    int nthreads;
+    std::vector<std::thread*> threads;
+    float g;
+
+    if(N<=0)
+    {
+        N = 0;
+    }
+    ret.resize(N);
+
+    if(ret.length==0)
+    {
+        return ret;
+    }
+
+    if(ret.length<narray_thread_sz)
+    {
+        for(I=0;I<ret.length;I++)
+        {
+            g = ((float)I)/((float)(N-1));
+            ret.data[I] = low*(1.0f-g) + high*g;
+        }
+    }
+    else
+    {
+        //threaded operation
+        nthreads = std::thread::hardware_concurrency();
+        nthreads = (nthreads<1) ? 1 : nthreads;
+        nthreads = (nthreads>narray_max_threads) ? narray_max_threads : nthreads;
+        //if(nthreads<1) nthreads=1;
+        //if(nthreads>narray_max_threads) nthreads = narray_max_threads;
+        threads.resize(nthreads);
+        for(J=0;J<nthreads;J++)
+        {
+            threads[J] = new(std::nothrow) std::thread(
+                &narray_linspacef_tf,
+                &ret,
+                low,high,N,
+                J,nthreads
+            );
+        }
+        for(J=0;J<nthreads;J++)
+        {
+            if(threads[J]!=NULL)
+            { 
+                threads[J]->join(); 
+                delete threads[J]; 
+                threads[J]= NULL;
+            }
+        }
+    }
+
+
+
+    return ret;
+}
+
+void narray_linspace_tf(
+    narray<double> *in,
+    const double low, const double high, const narray_size_t N,
+    int threadnum, int nthreads
+)
+{
+    int I0,I1,Is,I;
+    double g;
+
+    Is = in->length/nthreads;
+    I0 = Is*threadnum;
+    I1 = (threadnum>=(nthreads-1))? in->length: Is*(threadnum+1);
+    for(I=I0;I<I1;I++)
+    {
+        g = ((double)I)/((double)(N-1));
+        in->data[I] = low*(1.0-g) + high*g;
+    }
+}
+
+narray<double> linspace(const double low, const double high, narray_size_t N)
+{
+    narray<double> ret;
+    narray_size_t I;
+    int J;
+    int nthreads;
+    std::vector<std::thread*> threads;
+    double g;
+
+    if(N<=0)
+    {
+        N = 0;
+    }
+    ret.resize(N);
+
+    if(ret.length==0)
+    {
+        return ret;
+    }
+
+    if(ret.length<narray_thread_sz)
+    {
+        for(I=0;I<ret.length;I++)
+        {
+            g = ((double)I)/((double)(N-1));
+            ret.data[I] = low*(1.0-g) + high*g;
+        }
+    }
+    else
+    {
+        //threaded operation
+        nthreads = std::thread::hardware_concurrency();
+        nthreads = (nthreads<1) ? 1 : nthreads;
+        nthreads = (nthreads>narray_max_threads) ? narray_max_threads : nthreads;
+        //if(nthreads<1) nthreads=1;
+        //if(nthreads>narray_max_threads) nthreads = narray_max_threads;
+        threads.resize(nthreads);
+        for(J=0;J<nthreads;J++)
+        {
+            threads[J] = new(std::nothrow) std::thread(
+                &narray_linspace_tf,
+                &ret,
+                low,high,N,
+                J,nthreads
+            );
+        }
+        for(J=0;J<nthreads;J++)
+        {
+            if(threads[J]!=NULL)
+            { 
+                threads[J]->join(); 
+                delete threads[J]; 
+                threads[J]= NULL;
+            }
+        }
+    }
+
+
+
+    return ret;
+}
+
+
+}; //end namespaces
+};

src/amscppnarray/amscppnarray_random.cpp

@@ -0,0 +1,11 @@
+#include <amscppnarray/amscppnarray.hpp>
+
+namespace ams
+{
+namespace narray
+{
+
+
+    
+}; //end namespaces
+};

src/main.cpp

@@ -6,7 +6,7 @@ int main(int argc, char* argv[])
     printf("ams c++ numeric array library tests.\n");
     //ams::narray::test_narray1();
     //ams::narray::test_narray2();
-    ams::narray::test_narray3();
+    //ams::narray::test_narray3();
 
     return ret;
 }