common/array: added mean, variance and in-place subtraction of constant.
authorTito Dal Canton <tito.canton@ligo.org>
Wed, 23 Jan 2013 16:40:29 +0000 (17:40 +0100)
committerTito Dal Canton <tito.canton@ligo.org>
Wed, 23 Jan 2013 16:40:29 +0000 (17:40 +0100)
common/array: improved some OpenCL error checking.
cbc/findchirp: print SNR and chi2 stats for debugging.

cbc/findchirp/src/findchirp.cpp
common/include/array/GWArray.hpp
common/include/array/GWArray.hpp.old [deleted file]
common/include/array/GWArrayGeneric.hpp
common/src/GWArrayGeneric.cl

index 6c894d7..3e7af3c 100644 (file)
@@ -139,7 +139,7 @@ int main(int ac, char* av[]) {
     std::vector<GW::FrequencySeries<cReal> > strain_fft_segments;
 
     // Copying frame data to a complex strain time series
-    *strain =  *frame_data;
+    *strain = *frame_data;
 
     // Applying the dynamic exponent
     scaling_factor = pow(10, dynamic_exp);
@@ -249,8 +249,6 @@ int main(int ac, char* av[]) {
     std::vector <GW::FrequencySeries <cReal> *> out_series;
     out_series.push_back(h);
 
-    cReal * buff = new cReal[data_size];
-
     myTimer.Stop(t_precond);
 
     // Here the real filtering starts
@@ -311,6 +309,14 @@ int main(int ac, char* av[]) {
         myTimer.Continue(t_cluster);
         mySMCluster.FindClusters(cluster_args, &myIntrParams, snr, chi2, &myTriggerList);
         myTimer.Stop(t_cluster);
+
+        // print stats for debugging
+        snr->abs();
+        chi2->abs();
+        myLogger.Log(GW::DEBUG, __func__, "SNR mean %.2f variance %.2f", snr->mean().real(), snr->var().real());
+        if (number_of_chi2_bands > 0) {
+            myLogger.Log(GW::DEBUG, __func__, "Chi2 mean %.2f variance %.2f", chi2->mean().real(), chi2->var().real());
+        }
     }
 
     // Outputing the TriggerList
index ade3f6e..3da87e1 100644 (file)
@@ -32,7 +32,7 @@ using namespace std;
     err |= ocl.kernels[myop][idx.s[1]][idx.s[2]].setArg(1, par1); ocl.checkErr(err,  "cl::addKernel::setArg(1)");\
     err |= ocl.kernels[myop][idx.s[1]][idx.s[2]].setArg(2, par2); ocl.checkErr(err,  "cl::addKernel::setArg(2)");\
     err |= m_runtime->appQueues[0].enqueueNDRangeKernel(ocl.kernels[myop][idx.s[1]][idx.s[2]], cl::NullRange, globalSize, 256, NULL, &kernelEvent); \
-    ocl.checkErr(err,  "mymessage");\
+    ocl.checkErr(err, mymessage);\
     m_runtime->appQueues[0].finish();\
     return result;
 
@@ -63,7 +63,7 @@ using namespace std;
     err |= ocl.kernels[myop][idx.s[1]][idx.s[2]].setArg(2, batch_size); ocl.checkErr(err,  "cl::addKernel::setArg(3)"); \
     err |= m_runtime->appQueues[0].enqueueNDRangeKernel(ocl.kernels[myop][idx.s[1]][idx.s[2]], \
                                           cl::NullRange, globalSize, 256, NULL, &kernelEvent); \
-    ocl.checkErr(err,  "mymessage");\
+    ocl.checkErr(err, mymessage);\
     m_runtime->appQueues[0].finish();
 
 #define _M_GWARRAY_COMPOUND_ARRAY_NUM_ARITHMETIC_KERNELCALL(mymessage,par0,par1) \
@@ -73,12 +73,17 @@ using namespace std;
     cl::Event kernelEvent; \
        cl_uint globalSize = this->length; \
        cl_uint batch_size = 1;\
-       if ((idx.s[1] > 3) && (myop < 4)) { globalSize = globalSize * 2;}; \
-    err  = ocl.kernels[myop][idx.s[1]][idx.s[2]].setArg(0, par0); ocl.checkErr(err,  "cl::addKernel::setArg(0)"); \
-    err |= ocl.kernels[myop][idx.s[1]][idx.s[2]].setArg(1, par1); ocl.checkErr(err,  "cl::addKernel::setArg(1)"); \
-    err |= m_runtime->appQueues[0].enqueueNDRangeKernel(ocl.kernels[myop][idx.s[1]][idx.s[2]], \
-                                          cl::NullRange, globalSize, 256, NULL, &kernelEvent); \
-    ocl.checkErr(err,  "mymessage");\
+       if ((idx.s[1] > 3) && (myop < 4)) { \
+        globalSize = globalSize * 2; \
+    }; \
+    err = ocl.kernels[myop][idx.s[1]][idx.s[2]].setArg(0, par0); \
+    ocl.checkErr(err, "cl::addKernel::setArg(0)"); \
+    err = ocl.kernels[myop][idx.s[1]][idx.s[2]].setArg(1, par1); \
+    ocl.checkErr(err, "cl::addKernel::setArg(1)"); \
+    err = m_runtime->appQueues[0].enqueueNDRangeKernel( \
+        ocl.kernels[myop][idx.s[1]][idx.s[2]], \
+        cl::NullRange, globalSize, 256, NULL, &kernelEvent); \
+    ocl.checkErr(err, mymessage);\
     m_runtime->appQueues[0].finish();
 
 /******************************************************************************/
@@ -160,6 +165,7 @@ template <typename T> class Array : public ArrayBase
     Array <T>   operator *  (const T rhs) const;
 
     Array <T> & operator -= (const Array <T> & rhs);
+    Array <T> & operator -= (const T rhs);
     Array <T>   operator -  (const Array <T>) const;
     Array <T>   operator -  (const T rhs) const;
 
@@ -183,6 +189,8 @@ template <typename T> class Array : public ArrayBase
     void            invert();
     T               max();
        T               sum();
+    T               mean();
+    T               var();
 
        // Transfer methods
     cl_int          Copy(cl_int src_off, cl_int dest_off, cl_int size, Array <T> * dest);
@@ -219,6 +227,7 @@ enum arrayKernelEnum {
         gwarrayMulEqnum = 11,
         gwarrayDivnum   = 12,
         gwarrayDivEqnum = 13,
+        gwarraySubEqnum = 14,
         gwarrayFill     = 20,
         gwarrayConj     = 21,
         gwarraySqr      = 22,
@@ -454,8 +463,10 @@ template <typename T> cl_int Array <T>::Init() {
                                 = cl::Kernel(ocl.program, "uintSubEqual", &err);
     ocl.kernels[gwarraySubEq][gwcl_complex_float][gwcl_complex_float]
                                 = cl::Kernel(ocl.program, "floatSubEqual", &err);
+    ocl.checkErr(err,  "cl::Kernel(floatSubEqual)");
     ocl.kernels[gwarraySubEq][gwcl_complex_double][gwcl_complex_double]
                                 = cl::Kernel(ocl.program, "doubleSubEqual", &err);
+    ocl.checkErr(err,  "cl::Kernel(doubleSubEqual)");
 
     // (Array -= Num ) SUBSTRACTION KERNELS
     ocl.kernels[gwarraySubnum][gwcl_int][gwcl_int]
@@ -466,6 +477,18 @@ template <typename T> cl_int Array <T>::Init() {
                                 = cl::Kernel(ocl.program, "floatSubNum", &err);
     ocl.kernels[gwarraySubnum][gwcl_double][gwcl_double]
                                 = cl::Kernel(ocl.program, "doubleSubNum", &err);
+    ocl.kernels[gwarraySubEqnum][gwcl_float][gwcl_float]
+                                = cl::Kernel(ocl.program, "floatSubEqNum", &err);
+    ocl.checkErr(err, "cl::Kernel(floatSubEqNum)");
+    ocl.kernels[gwarraySubEqnum][gwcl_double][gwcl_double]
+                                = cl::Kernel(ocl.program, "doubleSubEqNum", &err);
+    ocl.checkErr(err, "cl::Kernel(doubleSubEqNum)");
+    ocl.kernels[gwarraySubEqnum][cF][cF]
+                                = cl::Kernel(ocl.program, "complexFloatSubEqNum", &err);
+    ocl.checkErr(err, "cl::Kernel(complexFloatSubEqNum)");
+    ocl.kernels[gwarraySubEqnum][cD][cD]
+                                = cl::Kernel(ocl.program, "complexDoubleSubEqNum", &err);
+    ocl.checkErr(err, "cl::Kernel(complexDoubleSubEqNum)");
 
     // (Array * Array ) MULTIPLICATION KERNELS
 
@@ -881,10 +904,16 @@ template <typename T> Array <T>  Array <T>::operator + (const T rhs) const {
 
 template <typename T> Array <T> & Array <T>::operator -= (const Array<T> &rhs) {
        cl_uint myop = gwarraySubEq;
-        _M_GWARRAY_COMPOUND_ARRAY_ARRAY_ARITHMETIC_KERNELCALL("Addeq kernel", this->data, rhs.data);
+        _M_GWARRAY_COMPOUND_ARRAY_ARRAY_ARITHMETIC_KERNELCALL("-= kernel", this->data, rhs.data);
        return *this;
 };
 
+template <typename T> Array <T> & Array <T>::operator -= (const T rhs) {
+    cl_uint myop = gwarraySubEqnum;
+        _M_GWARRAY_COMPOUND_ARRAY_NUM_ARITHMETIC_KERNELCALL("-= kernel", this->data, rhs);
+    return *this;
+};
+
 template <typename T> Array <T>  Array <T>::operator - (const Array<T> rhs) const {
        cl_uint myop = gwarraySub;
         M_GWARRAY_ARRAY_ARRAY_ARITHMETIC_KERNELCALL("Substraction kernel",this->data,rhs.data, result.data);
@@ -908,9 +937,9 @@ template <typename T> Array <T>  Array <T>::operator * (const Array<T> rhs) cons
         M_GWARRAY_ARRAY_ARRAY_ARITHMETIC_KERNELCALL("Multiplication kernel",this->data, rhs.data, result.data);
 };
 
-template <typename T> Array <T>  & Array <T>::operator *= (const T rhs) {
+template <typename T> Array <T> & Array <T>::operator *= (const T rhs) {
     cl_uint myop = gwarrayMulEqnum;
-         _M_GWARRAY_COMPOUND_ARRAY_NUM_ARITHMETIC_KERNELCALL("Multiplication kernel", this->data, rhs);
+        _M_GWARRAY_COMPOUND_ARRAY_NUM_ARITHMETIC_KERNELCALL("Multiplication kernel", this->data, rhs);
     return *this;
 };
 
@@ -1055,6 +1084,20 @@ template <typename T>  T  Array<T>::max()  {
     return mymax;
 }
 
+template <typename T>  T  Array<T>::mean(void)
+{
+    return sum() / length;
+}
+
+template <typename T>  T  Array<T>::var(void)
+{
+    Array delta(*this);
+    
+    delta -= mean();
+    delta.sqr();
+    return delta.mean();
+}
+
 
 #define _M_GWARRAY_0ARG_METHODS(mymessage,myoperation) {    \
     T mtype;                                                               \
diff --git a/common/include/array/GWArray.hpp.old b/common/include/array/GWArray.hpp.old
deleted file mode 100644 (file)
index aeb6a38..0000000
+++ /dev/null
@@ -1,826 +0,0 @@
-/*! \file
- *
- *   \brief Definition of the Array class and some of its derivatives
- *   \author Gergely Debreczeni, Karsten Wiesner
- *
- *   \version 1.0
- *   \date 2012m_runtime
- *
- */
-
-#ifndef _GWArray_H
-#define _GWArray_H
-
-#include <fstream>
-#include <iostream>
-#include <string>
-#include "array/GWArrayBase.hpp"
-#include "GWOpenCl.hpp"
-
-using namespace std;
-
-#define M_GWARRAY_ARRAY_ARRAY_ARITHMETIC_KERNELCALL(mymessage,par0,par1,par2) \
-       T mtype; \
-               cl_int err = 0; \
-       cl_uint4 idx = ocl.KernelSelector(mtype,mtype); \
-        cl::Event kernelEvent; \
-        Array <T> result(length, m_commif, m_logger, m_runtime); \
-       cl_uint globalSize = length ; \
-       if ((idx.s[1] > 3) && (myop < 4)) { globalSize = globalSize * 2;} \
-        err  = ocl.kernels[myop][idx.s[1]][idx.s[2]].setArg(0, par0); ocl.checkErr(err,  "cl::addKernel::setArg(0)");\
-        err |= ocl.kernels[myop][idx.s[1]][idx.s[2]].setArg(1, par1); ocl.checkErr(err,  "cl::addKernel::setArg(1)");\
-        err |= ocl.kernels[myop][idx.s[1]][idx.s[2]].setArg(2, par2); ocl.checkErr(err,  "cl::addKernel::setArg(2)");\
-        err |= m_runtime->appQueues[0].enqueueNDRangeKernel(ocl.kernels[myop][idx.s[1]][idx.s[2]], cl::NullRange, globalSize, 256, NULL, &kernelEvent); \
-        ocl.checkErr(err,  "mymessage");\
-        m_runtime->appQueues[0].finish();\
-        return result;
-
-/// Macro for compound array - array arithmetic operators
-/**
- This is a macro for compound arithmetic operators involving only two argument. A const
- reference on the right-hand-side (rhs) and the operand itself. Returns reference
- in order to allow operator chaining. Allows for non equal length operations. The
- length of the left-hand-side argument must be an integer multiple that of the right-hand-side
- one. This feature yields c.c. 30% in performance in case of long arrays and one order of
- magnitude speedup for short arrays, compared to a loop solution.
-*/
-
-#define _M_GWARRAY_COMPOUND_ARRAY_ARRAY_ARITHMETIC_KERNELCALL(mymessage,par0,par1) \
-       T mtype; \
-               cl_int err = 0; \
-       cl_uint4 idx = ocl.KernelSelector(mtype,mtype); \
-        cl::Event kernelEvent; \
-       cl_uint globalSize = this->length; \
-       cl_uint batch_size = par1.length; \
-       if ((idx.s[1] > 3) && (myop < 4)) { globalSize = globalSize * 2;}; \
-        err  = ocl.kernels[myop][idx.s[1]][idx.s[2]].setArg(0, par0); ocl.checkErr(err,  "cl::addKernel::setArg(0)"); \
-        err |= ocl.kernels[myop][idx.s[1]][idx.s[2]].setArg(1, par1.data); ocl.checkErr(err,  "cl::addKernel::setArg(1)"); \
-        err |= ocl.kernels[myop][idx.s[1]][idx.s[2]].setArg(2, batch_size); ocl.checkErr(err,  "cl::addKernel::setArg(3)"); \
-        err |= m_runtime->appQueues[0].enqueueNDRangeKernel(ocl.kernels[myop][idx.s[1]][idx.s[2]], \
-                                          cl::NullRange, globalSize, 256, NULL, &kernelEvent); \
-        ocl.checkErr(err,  "mymessage");\
-        m_runtime->appQueues[0].finish();
-
-#define _M_GWARRAY_COMPOUND_ARRAY_NUM_ARITHMETIC_KERNELCALL(mymessage,par0,par1) \
-       T mtype; \
-               cl_int err = 0; \
-       cl_uint4 idx = ocl.KernelSelector(mtype,mtype); \
-        cl::Event kernelEvent; \
-       cl_uint globalSize = this->length; \
-       cl_uint batch_size = 1;\
-       if ((idx.s[1] > 3) && (myop < 4)) { globalSize = globalSize * 2;}; \
-        err  = ocl.kernels[myop][idx.s[1]][idx.s[2]].setArg(0, par0); ocl.checkErr(err,  "cl::addKernel::setArg(0)"); \
-        err |= ocl.kernels[myop][idx.s[1]][idx.s[2]].setArg(1, par1); ocl.checkErr(err,  "cl::addKernel::setArg(1)"); \
-        err |= m_runtime->appQueues[0].enqueueNDRangeKernel(ocl.kernels[myop][idx.s[1]][idx.s[2]], \
-                                          cl::NullRange, globalSize, 256, NULL, &kernelEvent); \
-        ocl.checkErr(err,  "mymessage");\
-        m_runtime->appQueues[0].finish();
-
-
-/*! \namespace GW
- *
- * \brief All the objects of GWTools are defined in this namespace
- *
- */
-namespace GW {
-
-/// Array
-/**
- This is the basic object of GWTools. An Array stored on the GPU. Various
- member function are defined. TimeSeries and FrequencySeries are derived
- from this class. It is a templates class. It can be instantiated as
- int, uint, float, double.
-*/
-
-template <typename T> class Array : public ArrayBase
-{
-    private:
-
-       static OpenCLObject ocl;
-
-        cl_int Create(const size_t lngth);
-        cl_int Fillin(const cl_int value, size_t min, size_t max);
-        cl_int Fillin(const cl_uint value, size_t min, size_t max);
-        cl_int Fillin(const cl_float value, size_t min, size_t max);
-        cl_int Fillin(const cl_double value, size_t min, size_t max);
-        cl_int Fillin(const cl_complex_int value, size_t min, size_t max);
-        cl_int Fillin(const cl_complex_uint value, size_t min, size_t max);
-        cl_int Fillin(const cl_complex_float value, size_t min, size_t max);
-        cl_int Fillin(const cl_complex_double value, size_t min, size_t max);
-
-    public:
-/*
-    cl::Buffer data;
-       size_t length;
-       size_t elementsize;
-*/
-       ///static cl_int Init();
-              cl_int Init();
-
-    Array(const size_t,
-          CLIBase* commif,
-          LoggerBase* logger,
-          Runtime* runtime);
-
-    Array(const T * ptr, const size_t,
-          CLIBase* commif,
-          LoggerBase* logger,
-          Runtime* runtime);
-
-    Array(const T  value, const size_t,
-          CLIBase* commif,
-          LoggerBase* logger,
-          Runtime* runtime);
-
-    ~Array();
-    Array <T> & operator =  (const Array <T> &rhs);
-    Array <T> & operator =  (const T &rhs);
-
-    Array <T> & operator += (const Array <T> &rhs);
-    Array <T>   operator +  (const Array <T>) const;
-    Array <T>   operator +  (const T rhs ) const;
-
-    Array <T> & operator *= (const Array <T> & rhs);
-    Array <T>   operator *  (const Array <T>) const;
-    Array <T>   operator *  (const T rhs) const;
-
-    Array <T> & operator -= (const Array <T> & rhs);
-    Array <T>   operator -  (const Array <T>) const;
-    Array <T>   operator -  (const T rhs) const;
-
-    Array <T> & operator /= (const Array <T> & rhs);
-    Array <T> & operator /= (const T rhs);
-    Array <T>   operator /  (const Array <T>) const;
-    Array <T>   operator /  (const T rhs) const;
-
-
-    std::vector<Array<T> >     Segment(cl_uint segment_length, cl_uint overlap = 0);
-       T                       Sum();
-
-       // Mathematical methods
-       void                    conj();
-       void                    sqr();
-       void                    sin();
-       void                    cos();
-       void                    ln();
-       void                    pow(T exponent);
-       void            abs();
-
-       // Transfer methods
-       cl_int                  Copy(cl_int src_off, cl_int dest_off, cl_int size, Array <T> * dest);
-    cl_int          Fill(T value, size_t min, size_t max);
-    cl_int             Fetch(T * ptr);
-    cl_int             Dispatch(const T * ptr);
-
-       // Other methods
-       cl_int                  ReadFromFile(const char * filename, const char * format);
-       cl_int                  WriteToFile(const char * filename, const char * format);
-
-};
-
-enum arrayKernelEnum {
-        gwarrayAdd      = 0,
-        gwarrayAddEq    = 1,
-        gwarraySub      = 2,
-        gwarraySubEq    = 3,
-        gwarrayMul      = 4,
-        gwarrayMulEq    = 5,
-        gwarrayDiv      = 6,
-        gwarrayDivEq    = 7,
-        gwarrayAddnum   = 8,
-        gwarraySubnum   = 9,
-        gwarrayMulnum   = 10,
-        gwarrayDivnum   = 11,
-        gwarrayDivEqnum = 12,
-        gwarrayFill     = 20,
-        gwarrayConj     = 21,
-        gwarraySqr      = 22,
-        gwarraySin      = 23,
-        gwarrayCos      = 24,
-        gwarrayLn       = 25,
-        gwarrayLog      = 26,
-        gwarrayPow      = 27,
-        gwarrayAbs      = 28,
-        gwarraySumFirst = 40,
-        gwarraySumFinal = 41,
-        gwarraySum      = 42,
-};
-
-
-
-template <typename T> OpenCLObject Array<T>::ocl;
-
-
-/// Initializing and building kernels.
-/**
-  The Init() method builds and assigns the kernels for the first instantiation of
-  the template class for a given template parameter.
-*/
-
-template <typename T> cl_int Array <T>::Init() {
-
-       cl_int err;
-       if (ocl.isbuilt ==  false) {
-        string kernelFileName(m_commif->get_home());
-        kernelFileName += "/kernels/GWArray.cl";
-        ocl.Build(kernelFileName.c_str(), m_logger, m_runtime);
-        ocl.isbuilt = true;
-    }
-
-        ocl.kernels[gwarrayAdd][gwcl_int][gwcl_int] = cl::Kernel(ocl.program, "intAdd", &err);
-        ocl.kernels[gwarrayAdd][gwcl_uint][gwcl_uint] = cl::Kernel(ocl.program, "uintAdd", &err);
-        ocl.kernels[gwarrayAdd][gwcl_float][gwcl_float] = cl::Kernel(ocl.program, "floatAdd", &err);
-        ocl.kernels[gwarrayAdd][gwcl_double][gwcl_double] = cl::Kernel(ocl.program, "doubleAdd", &err);
-        ocl.kernels[gwarrayAdd][gwcl_complex_int][gwcl_complex_int] = cl::Kernel(ocl.program, "intAdd", &err);
-        ocl.kernels[gwarrayAdd][gwcl_complex_uint][gwcl_complex_uint] = cl::Kernel(ocl.program, "uintAdd", &err);
-        ocl.kernels[gwarrayAdd][gwcl_complex_float][gwcl_complex_float] = cl::Kernel(ocl.program, "floatAdd", &err);
-        ocl.kernels[gwarrayAdd][gwcl_complex_double][gwcl_complex_double] = cl::Kernel(ocl.program, "doubleAdd", &err);
-
-        ocl.kernels[gwarrayAddEq][gwcl_int][gwcl_int] = cl::Kernel(ocl.program, "intAddEqual", &err);
-        ocl.kernels[gwarrayAddEq][gwcl_uint][gwcl_uint] = cl::Kernel(ocl.program, "uintAddEqual", &err);
-        ocl.kernels[gwarrayAddEq][gwcl_float][gwcl_float] = cl::Kernel(ocl.program, "floatAddEqual", &err);
-        ocl.kernels[gwarrayAddEq][gwcl_double][gwcl_double] = cl::Kernel(ocl.program, "doubleAddEqual", &err);
-/*
-        ocl.kernels[gwarrayAddEq][gwcl_complex_int][gwcl_complex_int] = cl::Kernel(ocl.program, "complexntAddEqual", &err);
-        ocl.kernels[gwarrayAddEq][gwcl_complex_uint][gwcl_complex_uint] = cl::Kernel(ocl.program, "uintAddEqual", &err);
-        ocl.kernels[gwarrayAddEq][gwcl_complex_float][gwcl_complex_float] = cl::Kernel(ocl.program, "floatAddEqual", &err);
-        ocl.kernels[gwarrayAddEq][gwcl_complex_double][gwcl_complex_double] = cl::Kernel(ocl.program, "doubleAddEqual", &err);
-*/
-        ocl.kernels[gwarrayAddnum][gwcl_int][gwcl_int] = cl::Kernel(ocl.program, "intAddNum", &err);
-        ocl.kernels[gwarrayAddnum][gwcl_uint][gwcl_uint] = cl::Kernel(ocl.program, "uintAddNum", &err);
-        ocl.kernels[gwarrayAddnum][gwcl_float][gwcl_float] = cl::Kernel(ocl.program, "floatAddNum", &err);
-        ocl.kernels[gwarrayAddnum][gwcl_double][gwcl_double] = cl::Kernel(ocl.program, "doubleAddNum", &err);
-
-
-
-
-        ocl.kernels[gwarraySub][gwcl_int][gwcl_int] = cl::Kernel(ocl.program, "intSub", &err);
-        ocl.kernels[gwarraySub][gwcl_uint][gwcl_uint] = cl::Kernel(ocl.program, "uintSub", &err);
-        ocl.kernels[gwarraySub][gwcl_float][gwcl_float] = cl::Kernel(ocl.program, "floatSub", &err);
-        ocl.kernels[gwarraySub][gwcl_double][gwcl_double] = cl::Kernel(ocl.program, "doubleSub", &err);
-        ocl.kernels[gwarraySub][gwcl_complex_int][gwcl_complex_int] = cl::Kernel(ocl.program, "intSub", &err);
-        ocl.kernels[gwarraySub][gwcl_complex_uint][gwcl_complex_uint] = cl::Kernel(ocl.program, "uintSub", &err);
-        ocl.kernels[gwarraySub][gwcl_complex_float][gwcl_complex_float] = cl::Kernel(ocl.program, "floatSub", &err);
-        ocl.kernels[gwarraySub][gwcl_complex_double][gwcl_complex_double] = cl::Kernel(ocl.program, "doubleSub", &err);
-
-        ocl.kernels[gwarraySubEq][gwcl_int][gwcl_int] = cl::Kernel(ocl.program, "intSubEqual", &err);
-        ocl.kernels[gwarraySubEq][gwcl_uint][gwcl_uint] = cl::Kernel(ocl.program, "uintSubEqual", &err);
-        ocl.kernels[gwarraySubEq][gwcl_float][gwcl_float] = cl::Kernel(ocl.program, "floatSubEqual", &err);
-        ocl.kernels[gwarraySubEq][gwcl_double][gwcl_double] = cl::Kernel(ocl.program, "doubleSubEqual", &err);
-        ocl.kernels[gwarraySubEq][gwcl_complex_int][gwcl_complex_int] = cl::Kernel(ocl.program, "intSubEqual", &err);
-        ocl.kernels[gwarraySubEq][gwcl_complex_uint][gwcl_complex_uint] = cl::Kernel(ocl.program, "uintSubEqual", &err);
-        ocl.kernels[gwarraySubEq][gwcl_complex_float][gwcl_complex_float] = cl::Kernel(ocl.program, "floatSubEqual", &err);
-        ocl.kernels[gwarraySubEq][gwcl_complex_double][gwcl_complex_double] = cl::Kernel(ocl.program, "doubleSubEqual", &err);
-
-        ocl.kernels[gwarraySubnum][gwcl_int][gwcl_int] = cl::Kernel(ocl.program, "intSubNum", &err);
-        ocl.kernels[gwarraySubnum][gwcl_uint][gwcl_uint] = cl::Kernel(ocl.program, "uintSubNum", &err);
-        ocl.kernels[gwarraySubnum][gwcl_float][gwcl_float] = cl::Kernel(ocl.program, "floatSubNum", &err);
-        ocl.kernels[gwarraySubnum][gwcl_double][gwcl_double] = cl::Kernel(ocl.program, "doubleSubNum", &err);
-
-
-        ocl.kernels[gwarrayMul][gwcl_int][gwcl_int] = cl::Kernel(ocl.program, "intMul", &err);
-        ocl.kernels[gwarrayMul][gwcl_uint][gwcl_uint] = cl::Kernel(ocl.program, "uintMul", &err);
-        ocl.kernels[gwarrayMul][gwcl_float][gwcl_float] = cl::Kernel(ocl.program, "floatMul", &err);
-        ocl.kernels[gwarrayMul][gwcl_double][gwcl_double] = cl::Kernel(ocl.program, "doubleMul", &err);
-        ocl.kernels[gwarrayMul][gwcl_complex_int][gwcl_complex_int] = cl::Kernel(ocl.program, "complexIntMul", &err);
-        ocl.kernels[gwarrayMul][gwcl_complex_uint][gwcl_complex_uint] = cl::Kernel(ocl.program, "complexUintMul", &err);
-        ocl.kernels[gwarrayMul][gwcl_complex_float][gwcl_complex_float] = cl::Kernel(ocl.program, "complexFloatMul", &err);
-        ocl.kernels[gwarrayMul][gwcl_complex_double][gwcl_complex_double] = cl::Kernel(ocl.program, "complexDoubleMul", &err);
-
-        ocl.kernels[gwarrayMulEq][gwcl_int][gwcl_int] = cl::Kernel(ocl.program, "intMulEqual", &err);
-        ocl.kernels[gwarrayMulEq][gwcl_uint][gwcl_uint] = cl::Kernel(ocl.program, "uintMulEqual", &err);
-        ocl.kernels[gwarrayMulEq][gwcl_float][gwcl_float] = cl::Kernel(ocl.program, "floatMulEqual", &err);
-        ocl.kernels[gwarrayMulEq][gwcl_double][gwcl_double] = cl::Kernel(ocl.program, "doubleMulEqual", &err);
-        ocl.kernels[gwarrayMulEq][gwcl_complex_int][gwcl_complex_int] = cl::Kernel(ocl.program, "complexIntMulEqual", &err);
-        ocl.kernels[gwarrayMulEq][gwcl_complex_uint][gwcl_complex_uint] = cl::Kernel(ocl.program, "complexUintMulEqual", &err);
-        ocl.kernels[gwarrayMulEq][gwcl_complex_float][gwcl_complex_float] = cl::Kernel(ocl.program, "complexFloatMulEqual", &err);
-        ocl.kernels[gwarrayMulEq][gwcl_complex_double][gwcl_complex_double] = cl::Kernel(ocl.program, "complexDoubleMulEqual", &err);
-
-        ocl.kernels[gwarrayMulnum][gwcl_int][gwcl_int] = cl::Kernel(ocl.program, "intMulNum", &err);
-        ocl.kernels[gwarrayMulnum][gwcl_uint][gwcl_uint] = cl::Kernel(ocl.program, "uintMulNum", &err);
-        ocl.kernels[gwarrayMulnum][gwcl_float][gwcl_float] = cl::Kernel(ocl.program, "floatMulNum", &err);
-        ocl.kernels[gwarrayMulnum][gwcl_double][gwcl_double] = cl::Kernel(ocl.program, "doubleMulNum", &err);
-
-        ocl.kernels[gwarrayDiv][gwcl_int][gwcl_int] = cl::Kernel(ocl.program, "intDiv", &err);
-        ocl.kernels[gwarrayDiv][gwcl_uint][gwcl_uint] = cl::Kernel(ocl.program, "uintDiv", &err);
-        ocl.kernels[gwarrayDiv][gwcl_float][gwcl_float] = cl::Kernel(ocl.program, "floatDiv", &err);
-        ocl.kernels[gwarrayDiv][gwcl_double][gwcl_double] = cl::Kernel(ocl.program, "doubleDiv", &err);
-        ocl.kernels[gwarrayDiv][gwcl_complex_int][gwcl_complex_int] = cl::Kernel(ocl.program, "complexIntDiv", &err);
-        ocl.kernels[gwarrayDiv][gwcl_complex_uint][gwcl_complex_uint] = cl::Kernel(ocl.program, "complexUintDiv", &err);
-        ocl.kernels[gwarrayDiv][gwcl_complex_float][gwcl_complex_float] = cl::Kernel(ocl.program, "complexFloatDiv", &err);
-        ocl.kernels[gwarrayDiv][gwcl_complex_double][gwcl_complex_double] = cl::Kernel(ocl.program, "complexDoubleDiv", &err);
-
-        ocl.kernels[gwarrayDivEq][gwcl_int][gwcl_int] = cl::Kernel(ocl.program, "intDivEqual", &err);
-        ocl.checkErr(err,  "cl::Kernel(intDivEq)");
-        ocl.kernels[gwarrayDivEq][gwcl_uint][gwcl_uint] = cl::Kernel(ocl.program, "uintDivEqual", &err);
-        ocl.checkErr(err,  "cl::Kernel(uintDivEq)");
-        ocl.kernels[gwarrayDivEq][gwcl_float][gwcl_float] = cl::Kernel(ocl.program, "floatDivEqual", &err);
-        ocl.checkErr(err,  "cl::Kernel(floatDivEq)");
-        ocl.kernels[gwarrayDivEq][gwcl_double][gwcl_double] = cl::Kernel(ocl.program, "doubleDivEqual", &err);
-        ocl.checkErr(err,  "cl::Kernel(doubleDivEq)");
-        ocl.kernels[gwarrayDivEq][gwcl_complex_int][gwcl_complex_int] = cl::Kernel(ocl.program, "complexIntDivEqual", &err);
-        ocl.checkErr(err,  "cl::Kernel(complexIntDivEq)");
-        ocl.kernels[gwarrayDivEq][gwcl_complex_uint][gwcl_complex_uint] = cl::Kernel(ocl.program, "complexUintDivEqual", &err);
-        ocl.checkErr(err,  "cl::Kernel(complexUintDivEq)");
-        ocl.kernels[gwarrayDivEq][gwcl_complex_float][gwcl_complex_float] = cl::Kernel(ocl.program, "complexFloatDivEqual", &err);
-        ocl.checkErr(err,  "cl::Kernel(complexFloatDivEq)");
-        ocl.kernels[gwarrayDivEq][gwcl_complex_double][gwcl_complex_double] = cl::Kernel(ocl.program, "complexDoubleDivEqual", &err);
-        ocl.checkErr(err,  "cl::Kernel(complexDoubleDivEq)");
-
-        ocl.kernels[gwarrayDivnum][gwcl_int][gwcl_int] = cl::Kernel(ocl.program, "intDivNum", &err);
-        ocl.kernels[gwarrayDivnum][gwcl_uint][gwcl_uint] = cl::Kernel(ocl.program, "uintDivNum", &err);
-        ocl.kernels[gwarrayDivnum][gwcl_float][gwcl_float] = cl::Kernel(ocl.program, "floatDivNum", &err);
-        ocl.kernels[gwarrayDivnum][gwcl_double][gwcl_double] = cl::Kernel(ocl.program, "doubleDivNum", &err);
-
-        ocl.kernels[gwarrayDivEqnum][gwcl_complex_float][gwcl_complex_float] = cl::Kernel(ocl.program, "complexFloatDivEqNum", &err);
-        ocl.kernels[gwarrayDivEqnum][gwcl_complex_double][gwcl_complex_double] = cl::Kernel(ocl.program, "complexDoubleDivEqNum", &err);
-
-        ocl.kernels[gwarrayConj][gwcl_complex_int][gwcl_complex_int] = cl::Kernel(ocl.program, "complexIntConj", &err);
-        ocl.checkErr(err,  "cl::Kernel(complexIntConj)");
-        ocl.kernels[gwarrayConj][gwcl_complex_uint][gwcl_complex_uint] = cl::Kernel(ocl.program, "complexUintConj", &err);
-        ocl.checkErr(err,  "cl::Kernel(complexUintConj)");
-        ocl.kernels[gwarrayConj][gwcl_complex_float][gwcl_complex_float] = cl::Kernel(ocl.program, "complexFloatConj", &err);
-        ocl.checkErr(err,  "cl::Kernel(complexFloatConj)");
-        ocl.kernels[gwarrayConj][gwcl_complex_double][gwcl_complex_double] = cl::Kernel(ocl.program, "complexDoubleConj", &err);
-        ocl.checkErr(err,  "cl::Kernel(complexDoubleConj)");
-
-       ocl.kernels[gwarraySumFirst][gwcl_int][gwcl_int] = cl::Kernel(ocl.program, "intSumFirst", &err);
-        ocl.checkErr(err,  "cl::Kernel(intSumFirst)");
-        ocl.kernels[gwarraySumFirst][gwcl_uint][gwcl_uint] = cl::Kernel(ocl.program, "uintSumFirst", &err);
-        ocl.checkErr(err,  "cl::Kernel(uintSumFirst)");
-        ocl.kernels[gwarraySumFirst][gwcl_float][gwcl_float] = cl::Kernel(ocl.program, "floatSumFirst", &err);
-        ocl.checkErr(err,  "cl::Kernel(floatSumFirst)");
-        ocl.kernels[gwarraySumFirst][gwcl_double][gwcl_double] = cl::Kernel(ocl.program, "doubleSumFirst", &err);
-        ocl.checkErr(err,  "cl::Kernel(doubleSumFirst)");
-       ocl.kernels[gwarraySumFinal][gwcl_int][gwcl_int] = cl::Kernel(ocl.program, "intSumFinal", &err);
-        ocl.checkErr(err,  "cl::Kernel(intSumFinal)");
-        ocl.kernels[gwarraySumFinal][gwcl_uint][gwcl_uint] = cl::Kernel(ocl.program, "uintSumFinal", &err);
-        ocl.checkErr(err,  "cl::Kernel(uintSumFinal)");
-        ocl.kernels[gwarraySumFinal][gwcl_float][gwcl_float] = cl::Kernel(ocl.program, "floatSumFinal", &err);
-        ocl.checkErr(err,  "cl::Kernel(floatSumFinal)");
-        ocl.kernels[gwarraySumFinal][gwcl_double][gwcl_double] = cl::Kernel(ocl.program, "doubleSumFinal", &err);
-        ocl.checkErr(err,  "cl::Kernel(doubleSumFinal)");
-
-        ocl.kernels[gwarraySum][gwcl_int][gwcl_int] = cl::Kernel(ocl.program, "intSum", &err);
-        ocl.checkErr(err,  "cl::Kernel(intSum)");
-        ocl.kernels[gwarraySum][gwcl_uint][gwcl_uint] = cl::Kernel(ocl.program, "uintSum", &err);
-        ocl.checkErr(err,  "cl::Kernel(uintSum)");
-        ocl.kernels[gwarraySum][gwcl_float][gwcl_float] = cl::Kernel(ocl.program, "floatSum", &err);
-        ocl.checkErr(err,  "cl::Kernel(floatSum)");
-        ocl.kernels[gwarraySum][gwcl_double][gwcl_double] = cl::Kernel(ocl.program, "doubleSum", &err);
-        ocl.checkErr(err,  "cl::Kernel(doubleSum)");
-        ocl.kernels[gwarraySum][gwcl_complex_int][gwcl_complex_int] = cl::Kernel(ocl.program, "complexIntSum", &err);
-        ocl.checkErr(err,  "cl::Kernel(complexIntSum)");
-        ocl.kernels[gwarraySum][gwcl_complex_uint][gwcl_complex_uint] = cl::Kernel(ocl.program, "complexUintSum", &err);
-        ocl.checkErr(err,  "cl::Kernel(complexUintSum)");
-        ocl.kernels[gwarraySum][gwcl_complex_float][gwcl_complex_float] = cl::Kernel(ocl.program, "complexFloatSum", &err);
-        ocl.checkErr(err,  "cl::Kernel(complexFloatSum)");
-        ocl.kernels[gwarraySum][gwcl_complex_double][gwcl_complex_double] = cl::Kernel(ocl.program, "complexDoubleSum", &err);
-        ocl.checkErr(err,  "cl::Kernel(complexDoubleSum)");
-
-        ocl.kernels[gwarraySqr][gwcl_int][gwcl_int] = cl::Kernel(ocl.program, "intSqr", &err);
-        ocl.checkErr(err,  "cl::Kernel(intSqr)");
-        ocl.kernels[gwarraySqr][gwcl_uint][gwcl_uint] = cl::Kernel(ocl.program, "uintSqr", &err);
-        ocl.checkErr(err,  "cl::Kernel(uintSqr)");
-        ocl.kernels[gwarraySqr][gwcl_float][gwcl_float] = cl::Kernel(ocl.program, "floatSqr", &err);
-        ocl.checkErr(err,  "cl::Kernel(floatSqr)");
-        ocl.kernels[gwarraySqr][gwcl_double][gwcl_double] = cl::Kernel(ocl.program, "doubleSqr", &err);
-        ocl.checkErr(err,  "cl::Kernel(doubleSqr)");
-        ocl.kernels[gwarraySqr][gwcl_complex_int][gwcl_complex_int] = cl::Kernel(ocl.program, "complexIntSqr", &err);
-        ocl.checkErr(err,  "cl::Kernel(complexIntSqr)");
-        ocl.kernels[gwarraySqr][gwcl_complex_uint][gwcl_complex_uint] = cl::Kernel(ocl.program, "complexUintSqr", &err);
-        ocl.checkErr(err,  "cl::Kernel(complexUintSqr)");
-        ocl.kernels[gwarraySqr][gwcl_complex_float][gwcl_complex_float] = cl::Kernel(ocl.program, "complexFloatSqr", &err);
-        ocl.checkErr(err,  "cl::Kernel(complexFloarSqr)");
-        ocl.kernels[gwarraySqr][gwcl_complex_double][gwcl_complex_double] = cl::Kernel(ocl.program, "complexDoubleSqr", &err);
-        ocl.checkErr(err,  "cl::Kernel(complexDoubleSqr)");
-
-
-
-        ocl.kernels[gwarraySin][gwcl_float][gwcl_float] = cl::Kernel(ocl.program, "floatSin", &err);
-        ocl.checkErr(err,  "cl::Kernel(floatSin)");
-        ocl.kernels[gwarraySin][gwcl_double][gwcl_double] = cl::Kernel(ocl.program, "doubleSin", &err);
-        ocl.checkErr(err,  "cl::Kernel(doubleSin)");
-
-        ocl.kernels[gwarrayCos][gwcl_float][gwcl_float] = cl::Kernel(ocl.program, "floatCos", &err);
-        ocl.checkErr(err,  "cl::Kernel(floatCos)");
-        ocl.kernels[gwarrayCos][gwcl_double][gwcl_double] = cl::Kernel(ocl.program, "doubleCos", &err);
-        ocl.checkErr(err,  "cl::Kernel(doubleCos)");
-
-        ocl.kernels[gwarrayLn][gwcl_float][gwcl_float] = cl::Kernel(ocl.program, "floatLn", &err);
-        ocl.checkErr(err,  "cl::Kernel(floatLn)");
-        ocl.kernels[gwarrayLn][gwcl_double][gwcl_double] = cl::Kernel(ocl.program, "doubleLn", &err);
-        ocl.checkErr(err,  "cl::Kernel(doubleLn)");
-
-        ocl.kernels[gwarrayAbs][gwcl_complex_float][gwcl_complex_float] = cl::Kernel(ocl.program, "complexFloatAbs", &err);
-        ocl.checkErr(err,  "cl::Kernel(complexFloatAbs)");
-        ocl.kernels[gwarrayAbs][gwcl_complex_double][gwcl_complex_double] = cl::Kernel(ocl.program, "complexDoubleAbs", &err);
-        ocl.checkErr(err,  "cl::Kernel(complexDoubleAbs)");
-
-        ocl.kernels[gwarrayFill][gwcl_int][gwcl_int] = cl::Kernel(ocl.program, "intFill", &err);
-        ocl.checkErr(err,  "cl::Kernel(intFill)");
-        ocl.kernels[gwarrayFill][gwcl_uint][gwcl_uint] = cl::Kernel(ocl.program, "uintFill", &err);
-        ocl.checkErr(err,  "cl::Kernel(uintFill)");
-        ocl.kernels[gwarrayFill][gwcl_float][gwcl_float] = cl::Kernel(ocl.program, "floatFill", &err);
-        ocl.checkErr(err,  "cl::Kernel(floatFill)");
-        ocl.kernels[gwarrayFill][gwcl_double][gwcl_double] = cl::Kernel(ocl.program, "doubleFill", &err);
-        ocl.checkErr(err,  "cl::Kernel(doubleFill)");
-
-        ocl.kernels[gwarrayFill][gwcl_complex_int][gwcl_complex_int] = cl::Kernel(ocl.program, "complexIntFill", &err);
-        ocl.checkErr(err,  "cl::Kernel(complexIntFill)");
-        ocl.kernels[gwarrayFill][gwcl_complex_uint][gwcl_complex_uint] = cl::Kernel(ocl.program, "complexUintFill", &err);
-        ocl.checkErr(err,  "cl::Kernel(complexUintFill)");
-        ocl.kernels[gwarrayFill][gwcl_complex_float][gwcl_complex_float] = cl::Kernel(ocl.program, "complexFloatFill", &err);
-        ocl.checkErr(err,  "cl::Kernel(complexFloatFill)");
-        ocl.kernels[gwarrayFill][gwcl_complex_double][gwcl_complex_double] = cl::Kernel(ocl.program, "complexDoubleFill", &err);
-        ocl.checkErr(err,  "cl::Kernel(complexDoubleFill)");
-
-        ocl.isbuilt = true;
-
-        return (cl_int) 0;
-};
-
-
-template <typename T> Array <T> & Array <T>::operator = (const Array<T> &rhs) {
-               cl_int err = 0;
-        cl::Event cpyEvent;
-        if (this != &rhs) {
-        err = m_runtime->appQueues[0].enqueueCopyBuffer(rhs.data, data,
-                      (size_t) 0, (size_t) 0, (size_t) length * elementsize, NULL, &cpyEvent);
-        ocl.checkErr(err,  "cl:Copy Array");
-        m_runtime->appQueues[0].finish();
-       }
-       return *this;
-};
-
-template <typename T> Array <T> & Array <T>::operator = (const T &rhs) {
-       Fillin(rhs, 0, length - 1);
-       return *this;
-};
-
-
-template <typename T> Array <T>  Array <T>::operator + (const Array<T> rhs) const {
-       cl_uint myop = gwarrayAdd;
-        M_GWARRAY_ARRAY_ARRAY_ARITHMETIC_KERNELCALL("Addition kernel",this->data,rhs.data, result.data);
-};
-
-template <typename T> Array <T> & Array <T>::operator += (const Array<T> &rhs) {
-       cl_uint myop = gwarrayAddEq;
-        _M_GWARRAY_COMPOUND_ARRAY_ARRAY_ARITHMETIC_KERNELCALL("Addeq kernel", this->data, rhs);
-       return *this;
-};
-
-template <typename T> Array <T>  Array <T>::operator + (const T rhs) const {
-       cl_uint myop = gwarrayAddnum;
-        M_GWARRAY_ARRAY_ARRAY_ARITHMETIC_KERNELCALL("Addition kernel", this->data, rhs, result.data);
-};
-
-
-/* Substraction operators */
-
-template <typename T> Array <T> & Array <T>::operator -= (const Array<T> &rhs) {
-       cl_uint myop = gwarraySubEq;
-        _M_GWARRAY_COMPOUND_ARRAY_ARRAY_ARITHMETIC_KERNELCALL("Addeq kernel", this->data, rhs.data);
-       return *this;
-};
-
-template <typename T> Array <T>  Array <T>::operator - (const Array<T> rhs) const {
-       cl_uint myop = gwarraySub;
-        M_GWARRAY_ARRAY_ARRAY_ARITHMETIC_KERNELCALL("Substraction kernel",this->data,rhs.data, result.data);
-};
-
-template <typename T> Array <T>  Array <T>::operator - (const T rhs) const {
-       cl_uint myop = gwarraySubnum;
-        M_GWARRAY_ARRAY_ARRAY_ARITHMETIC_KERNELCALL("Substraction kernel", this->data, rhs, result.data);
-};
-
-/* Multiplication operators */
-
-template <typename T> Array <T> & Array <T>::operator *= (const Array<T> &rhs) {
-       cl_uint myop = gwarrayMulEq;
-        _M_GWARRAY_COMPOUND_ARRAY_ARRAY_ARITHMETIC_KERNELCALL("Multiply equal kernel", this->data, rhs);
-       return *this;
-};
-
-template <typename T> Array <T>  Array <T>::operator * (const Array<T> rhs) const {
-       cl_uint myop = gwarrayMul;
-        M_GWARRAY_ARRAY_ARRAY_ARITHMETIC_KERNELCALL("Multiplication kernel",this->data, rhs.data, result.data);
-};
-
-template <typename T> Array <T>  Array <T>::operator * (const T rhs) const {
-       cl_uint myop = gwarrayMulnum;
-        M_GWARRAY_ARRAY_ARRAY_ARITHMETIC_KERNELCALL("Multiplication kernel", this->data, rhs, result.data);
-};
-
-/* Division operators */
-
-template <typename T> Array <T> & Array <T>::operator /= (const Array<T> &rhs) {
-       cl_uint myop = gwarrayDivEq;
-        _M_GWARRAY_COMPOUND_ARRAY_ARRAY_ARITHMETIC_KERNELCALL("/= kernel", this->data, rhs);
-       return *this;
-};
-
-template <typename T> Array <T>  Array <T>::operator / (const Array<T> rhs) const {
-       cl_uint myop = gwarrayDiv;
-        M_GWARRAY_ARRAY_ARRAY_ARITHMETIC_KERNELCALL("Division kernel",this->data,rhs.data, result.data);
-};
-
-template <typename T> Array <T>  & Array <T>::operator /= (const T rhs) {
-       cl_uint myop = gwarrayDivEqnum;
-         _M_GWARRAY_COMPOUND_ARRAY_NUM_ARITHMETIC_KERNELCALL("Division kernel", this->data, rhs);
-       return *this;
-};
-
-template <typename T> Array <T>  Array <T>::operator / (const T rhs) const {
-       cl_uint myop = gwarrayDivnum;
-        M_GWARRAY_ARRAY_ARRAY_ARITHMETIC_KERNELCALL("Division kernel", this->data, rhs, result.data);
-};
-
-
-template <typename T> std::vector<Array<T> > Array<T>::Segment(cl_uint segment_length, cl_uint overlap) {
-       cl_int  err          = 0;
-       cl_uint num_segments = (length - segment_length) / ( segment_length - overlap ) + 1;
-       std::vector<Array<T> > result;
-       for (int i = 0; i < num_segments; i++) {
-               result.push_back(Array<T>(segment_length, m_commif, m_logger, m_runtime));
-               Copy(i * (segment_length - overlap), 0, segment_length, &result[i]);
-       }
-       return result;
-} 
-
-
-
-/* The SUM method */
-
-#define        _M_CPUSUM_REAL    {for (int i = 0; i < length; i++) {*sum += buff[i];}}
-#define        _M_CPUSUM_COMPLEX {for (int i = 0; i < length; i++) {sum->s[0] += (buff[i]).s[0]; sum->s[1] += (buff[i]).s[1];}}
-
-template <typename T> void CPUSum(int * buff, int * sum, uint length)  _M_CPUSUM_REAL
-template <typename T> void CPUSum(uint * buff, uint * sum, uint length)  _M_CPUSUM_REAL
-template <typename T> void CPUSum(float * buff, float * sum, uint length)  _M_CPUSUM_REAL
-template <typename T> void CPUSum(double * buff, double * sum, uint length)  _M_CPUSUM_REAL
-template <typename T> void CPUSum(cl_complex_int * buff, cl_complex_int * sum, uint length)  _M_CPUSUM_COMPLEX
-template <typename T> void CPUSum(cl_complex_uint * buff, cl_complex_uint * sum, uint length)  _M_CPUSUM_COMPLEX
-template <typename T> void CPUSum(cl_complex_float * buff, cl_complex_float * sum, uint length)  _M_CPUSUM_COMPLEX
-template <typename T> void CPUSum(cl_complex_double * buff, cl_complex_double * sum, uint length)  _M_CPUSUM_COMPLEX
-
-template <typename T>  T  Array<T>::Sum()  {
-
-       T mtype;                        
-       T sum;
-       cl_int  err     = 0;            
-       cl_uint myop    = gwarraySum;   
-       cl_uint4 toCall = ocl.KernelSelector(mtype,mtype);
-       std::vector<cl::Event> kernelEvents(1);
-                                       
-       // Creating a new buffer where the summation will happen
-        cl::Buffer tmpGPUBuffer = cl::Buffer(m_runtime->appContexts[0], CL_MEM_READ_WRITE, length * elementsize, NULL, &err);
-        ocl.checkErr(err,  "cl::Buffer(tmpGPUBuffer)");
-                                               
-       // Copy the array into the new buffer
-       err = m_runtime->appQueues[0].enqueueCopyBuffer(this->data, tmpGPUBuffer, (size_t) 0, (size_t) 0, (size_t) length * elementsize, NULL, kernelEvents.data());
-        ocl.checkErr(err,  "cl:Copy Array");
-                                       
-       // Setting the kernel argument 
-       err = ocl.kernels[myop][toCall.s[1]][toCall.s[1]].setArg(0, tmpGPUBuffer);
-       ocl.checkErr(err,  "cl::sumKernel::setArg(0)");
-                                       
-       // Start the summation          
-       uint minlength = 512;           
-       uint i = (uint) length / 2;     
-       while (i * 2 > minlength ) {    
-               err = m_runtime->appQueues[0].enqueueNDRangeKernel(ocl.kernels[myop][toCall.s[1]][toCall.s[1]], cl::NullRange, i, 256, NULL, kernelEvents.data());
-               ocl.checkErr(err,  "cl::CommandQueue::enqueueNDRangeKernel(SumKernel)");
-               i = i / 2;
-       }                               
-                                       
-       // Now finish the rest on the CPU
-       T * tmpCPUBuffer  = (T  * ) malloc(sizeof(T) * minlength);
-       err = m_runtime->appQueues[0].enqueueReadBuffer(tmpGPUBuffer, CL_TRUE, 0, minlength * elementsize, tmpCPUBuffer, NULL,  kernelEvents.data());
-        ocl.checkErr(err,  "ReadBuffer");
-               
-       CPUSum <T> (tmpCPUBuffer, &sum, minlength);
-       return sum;
-}
-
-
-
-#define _M_GWARRAY_0ARG_METHODS(mymessage,myoperation) { \
-        T mtype;                                        \
-       cl_uint myop    = myoperation;                   \
-        cl_int  err     = 0;                             \
-        cl_uint4 toCall = ocl.KernelSelector(mtype,mtype);   \
-        std::vector<cl::Event> kernelEvents(1);          \
-        cl::NDRange globalWorkSize = this->length;       \
-        cl::NDRange localWorkSize  = 256;                \
-        err = ocl.kernels[myop][toCall.s[1]][toCall.s[1]].setArg(0, this->data);       \
-        ocl.checkErr(err,  "cl::ConjKernel::setArg(0)");               \
-        err = m_runtime->appQueues[0].enqueueNDRangeKernel(ocl.kernels[myop][toCall.s[1]][toCall.s[1]], cl::NullRange, globalWorkSize, localWorkSize, NULL, kernelEvents.data());      \
-        ocl.checkErr(err,  "cl::CommandQueue::enqueueNDRangeKernel(mymessage)");} 
-
-
-template <typename T>  void Array<T>::conj()  _M_GWARRAY_0ARG_METHODS("Conjugation of array",gwarrayConj);
-template <typename T>  void Array<T>::sqr()   _M_GWARRAY_0ARG_METHODS("Conjugation of array",gwarraySqr);
-template <typename T>  void Array<T>::sin()   _M_GWARRAY_0ARG_METHODS("Conjugation of array",gwarraySin);
-template <typename T>  void Array<T>::cos()   _M_GWARRAY_0ARG_METHODS("Conjugation of array",gwarrayCos);
-template <typename T>  void Array<T>::ln()    _M_GWARRAY_0ARG_METHODS("Conjugation of array",gwarrayLn);
-template <typename T>  void Array<T>::abs()   _M_GWARRAY_0ARG_METHODS("Conjugation of array",gwarrayAbs);
-
-
-template <typename T>  void Array<T>::pow(T exponent) {
-       cl_uint myop    = gwarrayPow;
-       T mtype;
-       cl_int  err     = 0;
-       cl_uint4 toCall = ocl.KernelSelector(mtype,mtype);
-       std::vector<cl::Event> kernelEvents(1);
-
-       cl::NDRange globalWorkSize = this->length;
-       cl::NDRange localWorkSize  = 256;
-
-       err = ocl.kernels[myop][toCall.s[1]][toCall.s[1]].setArg(0, this->data); 
-       ocl.checkErr(err,  "cl::PowKernel::setArg(0)");
-       err = ocl.kernels[myop][toCall.s[1]][toCall.s[1]].setArg(1, &exponent); 
-       ocl.checkErr(err,  "cl::PowKernel::setArg(1)");
-
-        err = m_runtime->appQueues[0].enqueueNDRangeKernel(ocl.kernels[myop][toCall.s[1]][toCall.s[1]], cl::NullRange, globalWorkSize, localWorkSize, NULL, kernelEvents.data());
-        ocl.checkErr(err,  "cl::CommandQueue::enqueueNDRangeKernel(PowKernel)");
-}
-
-
-template <typename T> cl_int Array <T>::Copy(cl_int src_offset, cl_int dest_offset, cl_int size, Array <T> * dest) {
-       cl_int err = 0;
-        cl::Event cpyEvent;
-        if (this != dest) {
-               err = m_runtime->appQueues[0].enqueueCopyBuffer(this->data, dest->data,
-                      src_offset * this->elementsize, dest_offset * this->elementsize , size * this->elementsize, NULL, &cpyEvent);
-               ocl.checkErr(err,  "cl:Copy Array");
-               m_runtime->appQueues[0].finish();
-        }
-        return err;
-};
-
-template <typename T> cl_int Array <T>::Fill(T value, size_t min, size_t max) {
-    Fillin(value, min, max);
-    return 0;
-}
-
-template <typename T> cl_int Array <T>::Fetch(T * ptr) {
-       cl_int err = 0;
-        cl::Event cpyEvent;
-        err = m_runtime->appQueues[0].enqueueReadBuffer(data, CL_TRUE, 0, length*elementsize, ptr, NULL, &cpyEvent);
-        ocl.checkErr(err,  "Array::Fetch");
-        m_runtime->appQueues[0].finish();
-       return err;
-}
-
-
-template <typename T> cl_int Array <T>::Dispatch(const T * ptr) {
-       cl_int err = 0;
-        cl::Event cpyEvent;
-        err = m_runtime->appQueues[0].enqueueWriteBuffer(data, CL_TRUE, 0, length*elementsize, ptr, NULL, &cpyEvent);
-        ocl.checkErr(err,  "Array::Dispatch");
-        m_runtime->appQueues[0].finish();
-       return err;
-}
-
-template <typename T> cl_int Array <T>:: Create(const size_t lngth) {
-       cl_int err  = 0;
-       length      = lngth;
-    elementsize = sizeof(T);
-    data = cl::Buffer(m_runtime->appContexts[0], CL_MEM_READ_WRITE, length * elementsize, NULL, &err);
-    ocl.checkErr(err,  "Array::Create");
-       return err;
-}
-
-template <typename T> cl_int Array <T>::WriteToFile(const char * filename, const char * format) {
-       T * myBuff = (T  *) malloc(sizeof(T) * length); 
-       Fetch(myBuff);
-       if (format == "binary")  {
-               ofstream myFile (filename, ios::out | ios::binary | ios::trunc);
-               myFile.write( ( char *) myBuff, length * elementsize);
-               myFile.close();
-       } else if  (format == "ascii") {
-               ofstream myFile (filename, ios::out | ios::trunc);
-               for (int i = 0; i < length; i++) {
-                       myFile << myBuff[i] << endl;
-               }
-               myFile.close();
-       }
-};
-
-template <typename T> cl_int Array <T>::ReadFromFile(const char * filename, const char * format) {
-       T * myBuff = (T  * ) malloc(sizeof(T) * length); 
-       if (format == "binary") {               
-               ifstream myFile (filename, ios::in | ios::binary );
-               myFile.read( ( char *) myBuff, length * elementsize);
-               myFile.close();
-       } else if (format == "ascii") {
-               ifstream myFile (filename, ios::in );
-                for (int i = 0; i < length; i++) {
-                        myFile >> myBuff[i];
-                }
-                myFile.close();
-       }
-       Dispatch(myBuff);
-};
-
-#define _M_REAL_FILL_KERNEL(mytype)  {                                      \
-       cl_int    err = 0;                                                      \
-       mytype    mtype;                                                        \
-    cl::Event kernelEvent;                                                  \
-    cl_uint   gsize = (max - min) / 2;                                      \
-       cl_uint4  idx   = ocl.KernelSelector(mtype,mtype);                      \
-    ocl.kernels[gwarrayFill][idx.s[1]][idx.s[2]].setArg(0, this->data);     \
-    ocl.kernels[gwarrayFill][idx.s[1]][idx.s[2]].setArg(1, value);          \
-    ocl.kernels[gwarrayFill][idx.s[1]][idx.s[2]].setArg(2, min);            \
-    err = m_runtime->appQueues[0].enqueueNDRangeKernel                      \
-    (ocl.kernels[gwarrayFill][idx.s[1]][idx.s[2]], cl::NullRange,           \
-       gsize, cl::NullRange,  NULL, &kernelEvent);                             \
-    ocl.checkErr(err,  "fillKernel");                                       \
-    m_runtime->appQueues[0].finish();                                       \
-       return err; }                                                           \
-
-
-#define _M_COMPLEX_FILL_KERNEL(mytype)  {                               \
-       cl_int err = 0;                                                     \
-       mytype mtype;                                                       \
-    cl::Event kernelEvent;                                              \
-    cl_uint gsize = max - min;                                          \
-       cl_uint4 idx  = ocl.KernelSelector(mtype,mtype);                    \
-    ocl.kernels[gwarrayFill][idx.s[1]][idx.s[2]].setArg(0, this->data); \
-    ocl.kernels[gwarrayFill][idx.s[1]][idx.s[2]].setArg(1, value.s[0]); \
-    ocl.kernels[gwarrayFill][idx.s[1]][idx.s[2]].setArg(2, value.s[1]); \
-    ocl.kernels[gwarrayFill][idx.s[1]][idx.s[2]].setArg(3, min);        \
-    err = m_runtime->appQueues[0].enqueueNDRangeKernel                  \
-    (ocl.kernels[gwarrayFill][idx.s[1]][idx.s[2]], cl::NullRange,       \
-       gsize, cl::NullRange,  NULL, &kernelEvent);                         \
-    ocl.checkErr(err,  "fillKernel");                                   \
-    m_runtime->appQueues[0].finish();                                   \
-       return err; }                                                       \
-
-template <typename T> cl_int Array <T>:: 
-  Fillin(const cl_int value, size_t min, size_t max) _M_REAL_FILL_KERNEL(cl_int)
-template <typename T> cl_int Array <T>:: 
-  Fillin(const cl_uint value, size_t min, size_t max) _M_REAL_FILL_KERNEL(cl_uint)
-template <typename T> cl_int Array <T>:: 
-  Fillin(const cl_float value, size_t min, size_t max) _M_REAL_FILL_KERNEL(cl_float)
-template <typename T> cl_int Array <T>:: 
-  Fillin(const cl_double value, size_t min, size_t max) _M_REAL_FILL_KERNEL(cl_double)
-template <typename T> cl_int Array <T>:: 
-  Fillin(const cl_complex_int value, size_t min, size_t max) _M_COMPLEX_FILL_KERNEL(cl_complex_int)
-template <typename T> cl_int Array <T>:: 
- Fillin(const cl_complex_uint value, size_t min, size_t max) _M_COMPLEX_FILL_KERNEL(cl_complex_uint)
-template <typename T> cl_int Array <T>:: 
- Fillin(const cl_complex_float value, size_t min, size_t max) _M_COMPLEX_FILL_KERNEL(cl_complex_float)
-template <typename T> cl_int Array <T>:: 
- Fillin(const cl_complex_double value, size_t min, size_t max) _M_COMPLEX_FILL_KERNEL(cl_complex_double)
-
-// Constructors/Destructors:
-    
-template <typename T> Array <T>:: Array(const size_t lngth, CLIBase* commif, 
-                                                            LoggerBase* logger,
-                                                            Runtime* runtime) :
-    ArrayBase(commif, logger, runtime)
-{
-       Init();
-       length = lngth;
-       elementsize = sizeof(T);
-    Create(length);
-};
-
-template <typename T> Array <T>:: Array(const T * ptr, const size_t lngth, 
-                                        CLIBase* commif, 
-                                        LoggerBase* logger,
-                                        Runtime* runtime) :
-    ArrayBase(commif, logger, runtime)
-{
-       Init();
-       length = lngth;
-       elementsize = sizeof(T);
-
-    Create(length);
-    Dispatch(ptr);
-};
-
-template <typename T> Array <T>:: Array(const T value, const size_t lngth, 
-                                        CLIBase* commif, 
-                                        LoggerBase* logger,
-                                        Runtime* runtime) :
-    ArrayBase(commif, logger, runtime)
-{
-       Init();
-       length = lngth;
-       elementsize = sizeof(T);
-
-    Create(length);
-       Fillin(value,0, length - 1);
-};
-
-template <typename T> Array <T>:: ~Array() 
-{};
-
-
-
-} // End of GW namespace
-
-#endif
index 0c15728..9c28571 100644 (file)
@@ -395,6 +395,18 @@ template <typename T> cl_int Array <T>::Init() {
                                 = cl::Kernel(ocl.program, "floatSubNum", &err);
     ocl.kernels[gwarraySubnum][gwcl_double][gwcl_double]
                                 = cl::Kernel(ocl.program, "doubleSubNum", &err);
+    ocl.kernels[gwarraySubEqnum][gwcl_float][gwcl_float]
+                                = cl::Kernel(ocl.program, "floatSubEqualNum", &err);
+    ocl.checkErr(err, "cl::Kernel(floatSubEqualNum)");
+    ocl.kernels[gwarraySubEqnum][gwcl_double][gwcl_double]
+                                = cl::Kernel(ocl.program, "doubleSubEqualNum", &err);
+    ocl.checkErr(err, "cl::Kernel(doubleSubEqualNum)");
+    ocl.kernels[gwarraySubEqnum][gwcl_complex_float][gwcl_complex_float]
+                                = cl::Kernel(ocl.program, "complexFloatSubEqualNum", &err);
+    ocl.checkErr(err, "cl::Kernel(complexFloatSubEqualNum)");
+    ocl.kernels[gwarraySubEqnum][gwcl_complex_double][gwcl_complex_double]
+                                = cl::Kernel(ocl.program, "complexDoubleSubEqualNum", &err);
+    ocl.checkErr(err, "cl::Kernel(complexDoubleSubEqualNum)");
 
     // (Array * Array ) MULTIPLICATION KERNELS
 
index 738c63f..25acceb 100644 (file)
@@ -193,7 +193,7 @@ __kernel void uintSub(__global unsigned int * a, __global unsigned int * b, __gl
 __kernel void floatSub(__global float * a, __global float * b, __global float * c)  _M_SUB_KERNEL\r
 __kernel void doubleSub(__global double * a, __global double * b, __global double * c)  _M_SUB_KERNEL\r
 \r
-/* Kernels for compbound vector substraction (-=) */\r
+/* Kernels for compound vector substraction (-=) */\r
 \r
 #define _M_SUBEQUAL_KERNEL              \\r
 {                                       \\r
@@ -219,6 +219,37 @@ __kernel void uintSubNum(__global unsigned int * a, unsigned int b, __global uns
 __kernel void floatSubNum(__global float * a, float b, __global float * c)  _M_SUBNUM_KERNEL \r
 __kernel void doubleSubNum(__global double * a, double b, __global double * c)  _M_SUBNUM_KERNEL \r
 \r
+/* Kernels for compound vector subtraction with a number */\r
+\r
+__kernel void floatSubEqNum(__global float *a, float b)\r
+{\r
+    int gid = get_global_id(0);\r
+\r
+    a[gid] -= b;\r
+}\r
+\r
+__kernel void doubleSubEqNum(__global double *a, double b)\r
+{\r
+    int gid = get_global_id(0);\r
+\r
+    a[gid] -= b;\r
+}\r
+\r
+__kernel void complexFloatSubEqNum(__global float *a, float2 b)\r
+{\r
+    int gid = get_global_id(0) * 2;\r
+\r
+    a[gid] -= b.s0;\r
+    a[gid + 1] -= b.s1;\r
+}\r
+\r
+__kernel void complexDoubleSubEqNum(__global double *a, double2 b)\r
+{\r
+    int gid = get_global_id(0) * 2;\r
+\r
+    a[gid] -= b.s0;\r
+    a[gid + 1] -= b.s1;\r
+}\r
 \r
 /*  MULTIPLICATION */\r
 \r