Added a few missing kernels.
authorTito Dal Canton <tito.canton@ligo.org>
Fri, 25 Jan 2013 17:08:56 +0000 (18:08 +0100)
committerTito Dal Canton <tito.canton@ligo.org>
Fri, 25 Jan 2013 17:08:56 +0000 (18:08 +0100)
Fixed bug in Array::var().
Improved cl_complex_* classes.

common/include/array/GWArray.hpp
common/include/types/GWTypes.hpp
common/src/GWArrayGeneric.cl

index 430dcab..e175e5d 100644 (file)
@@ -337,6 +337,10 @@ template <typename T> cl_int Array <T>::Init() {
          o third is the right hand side argument type
     */
 
+#define STORE_KERNEL(op, t1, t2, name) \
+    ocl.kernels[op][t1][t2] = cl::Kernel(ocl.program, name, &err); \
+    ocl.checkErr(err, "cl::Kernel(" name ")");
+
     // COPY CONSTRUCTOR KERNELS WHEN CASTING IS NECESSARY
 
     K = gwarrayAssign;
@@ -419,6 +423,10 @@ template <typename T> cl_int Array <T>::Init() {
                                 = cl::Kernel(ocl.program, "floatAddEqual", &err);
     ocl.kernels[gwarrayAddEq][gwcl_double][gwcl_double]
                                 = cl::Kernel(ocl.program, "doubleAddEqual", &err);
+    ocl.kernels[gwarrayAddEq][cF][cF]
+                                = cl::Kernel(ocl.program, "complexFloatAddEqual", &err);
+    ocl.kernels[gwarrayAddEq][cD][cD]
+                                = cl::Kernel(ocl.program, "complexDoubleAddEqual", &err);
 
     // (Array + Number ) ADDITION KERNELS
     ocl.kernels[gwarrayAddnum][gwcl_int][gwcl_int]
@@ -469,7 +477,7 @@ template <typename T> cl_int Array <T>::Init() {
                                 = cl::Kernel(ocl.program, "doubleSubEqual", &err);
     ocl.checkErr(err,  "cl::Kernel(doubleSubEqual)");
 
-    // (Array -= Num ) SUBSTRACTION KERNELS
+    // (Array - Num ) SUBSTRACTION KERNELS
     ocl.kernels[gwarraySubnum][gwcl_int][gwcl_int]
                                 = cl::Kernel(ocl.program, "intSubNum", &err);
     ocl.kernels[gwarraySubnum][gwcl_uint][gwcl_uint]
@@ -478,37 +486,28 @@ 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)");
+    ocl.kernels[gwarraySubnum][cF][cF]
+                                = cl::Kernel(ocl.program, "complexFloatSubNum", &err);
+    ocl.kernels[gwarraySubnum][cD][cD]
+                                = cl::Kernel(ocl.program, "complexDoubleSubNum", &err);
+
+    // (Array -= Num ) SUBSTRACTION KERNELS
+
+    STORE_KERNEL(gwarraySubEqnum, gwcl_float, gwcl_float, "floatSubEqNum")
+    STORE_KERNEL(gwarraySubEqnum, gwcl_double, gwcl_double, "doubleSubEqNum")
+    STORE_KERNEL(gwarraySubEqnum, gwcl_complex_float, gwcl_complex_float, "complexFloatSubEqNum")
+    STORE_KERNEL(gwarraySubEqnum, gwcl_complex_double, gwcl_complex_double, "complexDoubleSubEqNum")
 
     // (Array * Array ) MULTIPLICATION KERNELS
 
-    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);
+    STORE_KERNEL(gwarrayMul, gwcl_int, gwcl_int, "intMul")
+    STORE_KERNEL(gwarrayMul, gwcl_uint, gwcl_uint, "uintMul")
+    STORE_KERNEL(gwarrayMul, gwcl_float, gwcl_float, "floatMul")
+    STORE_KERNEL(gwarrayMul, gwcl_double, gwcl_double, "doubleMul")
+    STORE_KERNEL(gwarrayMul, gwcl_complex_int, gwcl_complex_int, "complexIntMul")
+    STORE_KERNEL(gwarrayMul, gwcl_complex_uint, gwcl_complex_uint, "complexUintMul")
+    STORE_KERNEL(gwarrayMul, gwcl_complex_float, gwcl_complex_float, "complexFloatMul")
+    STORE_KERNEL(gwarrayMul, gwcl_complex_double, gwcl_complex_double, "complexDoubleMul")
 
     // (Array *= Array ) MULTIPLICATION KERNELS
     ocl.kernels[gwarrayMulEq][gwcl_int][gwcl_int]
@@ -1093,39 +1092,54 @@ template <typename T>  T  Array<T>::mean(void)
 
 template <typename T>  T  Array<T>::var(void)
 {
-    Array delta(*this);
+    Array<T> tmp(length, m_commif, m_logger, m_runtime);
     
-    delta -= mean();
-    delta.sqr();
-    return delta.mean();
+    tmp = *this - mean();
+    tmp.sqr();
+    return tmp.mean();
 }
 
 
-#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)");\
+#define _M_GWARRAY_0ARG_METHODS(mymessage, myoperation) \
+{ \
+    T mtype = 0; \
+    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, mymessage " cl::Kernel::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("array::conj()", gwarrayConj);
+
+template <typename T>  void Array<T>::sqr()
+    _M_GWARRAY_0ARG_METHODS("array::sqr()", gwarraySqr);
+
+template <typename T>  void Array<T>::sqrt()
+    _M_GWARRAY_0ARG_METHODS("array::sqrt()", gwarraySqrt);
+
+template <typename T>  void Array<T>::sin()
+    _M_GWARRAY_0ARG_METHODS("array::sin()", gwarraySin);
+
+template <typename T>  void Array<T>::cos()
+    _M_GWARRAY_0ARG_METHODS("array::cos()", gwarrayCos);
+
+template <typename T>  void Array<T>::ln()
+    _M_GWARRAY_0ARG_METHODS("array::ln()", gwarrayLn);
+
+template <typename T>  void Array<T>::abs()
+    _M_GWARRAY_0ARG_METHODS("array::abs()", gwarrayAbs);
 
-template <typename T>  void Array<T>::conj()  _M_GWARRAY_0ARG_METHODS("array::conj()",gwarrayConj);
-template <typename T>  void Array<T>::sqr()   _M_GWARRAY_0ARG_METHODS("array::sqr()",gwarraySqr);
-template <typename T>  void Array<T>::sqrt()  _M_GWARRAY_0ARG_METHODS("array::sqrt()",gwarraySqrt);
-template <typename T>  void Array<T>::sin()   _M_GWARRAY_0ARG_METHODS("array::sin()",gwarraySin);
-template <typename T>  void Array<T>::cos()   _M_GWARRAY_0ARG_METHODS("array::cos()",gwarrayCos);
-template <typename T>  void Array<T>::ln()    _M_GWARRAY_0ARG_METHODS("array::ln()",gwarrayLn);
-template <typename T>  void Array<T>::abs()   _M_GWARRAY_0ARG_METHODS("array::abs()",gwarrayAbs);
-template <typename T>  void Array<T>::invert()_M_GWARRAY_0ARG_METHODS("array::invert()",gwarrayInvert);
+template <typename T>  void Array<T>::invert()
+    _M_GWARRAY_0ARG_METHODS("array::invert()", gwarrayInvert);
 
 
 template <typename T>  void Array<T>::pow(T exponent) {
index 113d77e..ee1f053 100644 (file)
@@ -147,8 +147,8 @@ typedef union cl_complex_int
         return ret;
     }
 
-    cl_int real() {return s[0];}
-    cl_int imag() {return s[1];}
+    cl_int real() const { return s[0]; }
+    cl_int imag() const { return s[1]; }
 
     #if defined( __GNUC__) && ! defined( __STRICT_ANSI__ )
    __extension__ struct{ cl_int  x, y; };
@@ -292,8 +292,8 @@ typedef union cl_complex_uint
         return ret;
     }
 
-    cl_uint real() {return s[0];}
-    cl_uint imag() {return s[1];}
+    cl_uint real() const { return s[0]; }
+    cl_uint imag() const { return s[1]; }
 
     #if defined( __GNUC__) && ! defined( __STRICT_ANSI__ )
    __extension__ struct{ cl_uint  x, y; };
@@ -320,11 +320,11 @@ typedef union cl_complex_float
 {
     cl_float  CL_ALIGNED(8) s[2];
 
-    cl_complex_float ()          { s[0] = (float) 0;   s[1] = 0; }
-    cl_complex_float (int src)   { s[0] = (float) src; s[1] = 0; }
-    cl_complex_float (uint src)  { s[0] = (float) src; s[1] = 0; }
-    cl_complex_float (float src) { s[0] = (float) src; s[1] = 0; }
-    cl_complex_float (double src){ s[0] = (float) src; s[1] = 0; }
+    cl_complex_float() { s[0] = 0; s[1] = 0; }
+    cl_complex_float(int re, int im=0) { s[0] = (float)re; s[1] = (float)im; }
+    cl_complex_float(uint re, uint im=0) { s[0] = (float)re; s[1] = (float)im; }
+    cl_complex_float(float re, float im=0) { s[0] = re; s[1] = im; }
+    cl_complex_float(double re, double im=0) { s[0] = (float)re; s[1] = (float)im; }
 
     cl_complex_float & operator  = (const cl_complex_float & rhs) {
         s[0] = rhs.s0;
@@ -348,6 +348,16 @@ typedef union cl_complex_float
         s[0] = (float) rhs; s[1] = 0; return *this;
     };
 
+    bool operator == (const cl_complex_float & rhs)
+    {
+        return (s[0] == rhs.s[0]) && (s[1] == rhs.s[1]);
+    }
+
+    bool operator != (const cl_complex_float & rhs)
+    {
+        return (s[0] != rhs.s[0]) || (s[1] != rhs.s[1]);
+    }
+
 
     cl_complex_float & operator  += (const cl_complex_float & rhs) {
         s[0] += rhs.s[0];
@@ -371,7 +381,7 @@ typedef union cl_complex_float
         return *this;
     };
 
-    cl_complex_float & operator  *= (const cl_complex_float & rhs) {
+    cl_complex_float & operator *= (const cl_complex_float & rhs) {
         cl_complex_float ret = *this;
         ret.s[0] = s[0] * rhs.s[0] - s[1] * rhs.s[1];
         ret.s[1] = s[1] * rhs.s[0] - s[0] * rhs.s[1];
@@ -379,13 +389,13 @@ typedef union cl_complex_float
         return *this;
     };
 
-    cl_complex_float & operator  *= (const cl_float & rhs) {
+    cl_complex_float & operator *= (const cl_float & rhs) {
         s[0] = s[0] * rhs;
         s[1] = s[1] * rhs;
         return *this;
     };
 
-    cl_complex_float & operator  /= (const cl_complex_float & rhs) {
+    cl_complex_float & operator /= (const cl_complex_float & rhs) {
         cl_float nominator   = (rhs.s[0] * rhs.s[0] + rhs.s[1] * rhs.s[1]);
         cl_complex_float ret = *this;
         ret.s[0] = (s[0] * rhs.s[0] + s[1] * rhs.s[1]) / nominator;
@@ -394,7 +404,7 @@ typedef union cl_complex_float
         return *this;
     };
 
-    cl_complex_float & operator  /= (const cl_float & rhs) {
+    cl_complex_float & operator /= (const cl_float & rhs) {
         s[0] /= rhs;
         s[1] /= rhs;
         return *this;
@@ -448,22 +458,21 @@ typedef union cl_complex_float
         return ret;
     }
 
-    cl_float abs() {
-        return sqrt(s[0] * s[0] + s[1] * s[1]);
-    }
+    cl_float abs() const { return sqrt(s[0] * s[0] + s[1] * s[1]); }
 
-    cl_float real() {return s[0];}
-    cl_float imag() {return s[1];}
+    cl_float real() const { return s[0]; }
 
-    #if defined( __GNUC__) && ! defined( __STRICT_ANSI__ )
-   __extension__ struct{ cl_float  x, y; };
-   __extension__ struct{ cl_float  s0, s1; };
-   __extension__ struct{ cl_float  lo, hi; };
-    #endif
+    cl_float imag() const { return s[1]; }
 
-    #if defined( __CL_FLOAT2__)
-    __cl_float2     v2;
-    #endif
+#if defined( __GNUC__) && ! defined( __STRICT_ANSI__ )
+    __extension__ struct { cl_float x, y; };
+    __extension__ struct { cl_float s0, s1; };
+    __extension__ struct { cl_float lo, hi; };
+#endif
+
+#if defined( __CL_FLOAT2__)
+    __cl_float2 v2;
+#endif
 } cl_complex_float;
 
 // Overloaded functions for cl_complex_float
@@ -485,11 +494,11 @@ typedef union cl_complex_double
 {
     cl_double  CL_ALIGNED(16) s[2];
 
-    cl_complex_double ()          { s[0] = (double) 0;   s[1] = 0; }
-    cl_complex_double (int src)   { s[0] = (double) src; s[1] = 0; }
-    cl_complex_double (uint src)  { s[0] = (double) src; s[1] = 0; }
-    cl_complex_double (float src) { s[0] = (double) src; s[1] = 0; }
-    cl_complex_double (double src){ s[0] = (double) src; s[1] = 0; }
+    cl_complex_double() { s[0] = (double) 0; s[1] = 0; }
+    cl_complex_double(int re, int im=0) { s[0] = (double)re; s[1] = (double)im; }
+    cl_complex_double(uint re, uint im=0) { s[0] = (double)re; s[1] = (double)im; }
+    cl_complex_double(float re, float im=0) { s[0] = (double)re; s[1] = (double)im; }
+    cl_complex_double(double re, double im=0) { s[0] = re; s[1] = im; }
 
     cl_complex_double & operator  = (const cl_complex_double & rhs) {
         s[0] = rhs.s0;
@@ -602,12 +611,10 @@ typedef union cl_complex_double
         return ret;
     }
 
-    cl_double abs() {
-        return sqrt(s[0] * s[0] + s[1] * s[1]);
-    }
+    cl_double abs() const { return sqrt(s[0] * s[0] + s[1] * s[1]); }
 
-    cl_double real() {return s[0];}
-    cl_double imag() {return s[1];}
+    cl_double real() const { return s[0]; }
+    cl_double imag() const { return s[1]; }
 
     #if defined( __GNUC__) && ! defined( __STRICT_ANSI__ )
    __extension__ struct{ cl_double  x, y; };
@@ -639,12 +646,32 @@ inline std::istream & operator>> (std::istream &in,  cl_complex_int  &data) { ch
 inline std::ostream & operator<< (std::ostream &out, const cl_complex_uint &data) { out << data.s0 << " " << data.s1; return out;};
 inline std::istream & operator>> (std::istream &in,  cl_complex_uint  &data) { char delim; in  >> data.s0 >> delim >> data.s1; return in;};
 
-inline std::ostream & operator<< (std::ostream &out, const cl_complex_float &data) { out << data.s0 << " " << data.s1; return out;};
-inline std::istream & operator>> (std::istream &in,  cl_complex_float  &data) { char delim; in  >> data.s0 >> delim >> data.s1; return in;};
+inline std::ostream & operator<< (std::ostream &out, const cl_complex_float &data)
+{
+    out.precision(10);
+    out << data.s0 << " " << data.s1;
+    return out;
+}
 
-inline std::ostream & operator<< (std::ostream &out, const cl_complex_double &data) { out << data.s0 << " " << data.s1; return out;};
-inline std::istream & operator>> (std::istream &in,  cl_complex_double  &data) { char delim; in  >> data.s0 >> delim >> data.s1; return in;};
+inline std::istream & operator>> (std::istream &in, cl_complex_float &data)
+{
+    char delim;
+    in  >> data.s0 >> delim >> data.s1;
+    return in;
+}
 
+inline std::ostream & operator<< (std::ostream &out, const cl_complex_double &data)
+{
+    out.precision(10);
+    out << data.s0 << " " << data.s1;
+    return out;
+}
+
+inline std::istream & operator>> (std::istream &in, cl_complex_double &data)
+{
+    char delim; in  >> data.s0 >> delim >> data.s1;
+    return in;
+}
 
 #endif // _GWTYPES_H
 
index 25acceb..317dfd9 100644 (file)
@@ -164,6 +164,22 @@ __kernel void intAddEqual(__global int * a, __global int * b, int batch_size) _M
 __kernel void uintAddEqual(__global uint * a, __global uint * b, int batch_size) _M_ADDEQUAL_KERNEL\r
 __kernel void floatAddEqual(__global float * a, __global float * b, int batch_size) _M_ADDEQUAL_KERNEL\r
 __kernel void doubleAddEqual(__global double * a, __global double * b, int batch_size) _M_ADDEQUAL_KERNEL\r
+__kernel void complexFloatAddEqual(__global float *a, __global float *b, int batch_size)\r
+{\r
+    int gid = get_global_id(0) % batch_size;\r
+    int rhs_idx = gid % batch_size;\r
+    \r
+    a[2 * gid] += b[2 * rhs_idx];\r
+    a[2 * gid + 1] += b[2 * rhs_idx + 1];\r
+}\r
+__kernel void complexDoubleAddEqual(__global float *a, __global float *b, int batch_size)\r
+{\r
+    int gid = get_global_id(0) % batch_size;\r
+    int rhs_idx = gid % batch_size;    \r
+    \r
+    a[2 * gid] += b[2 * rhs_idx];\r
+    a[2 * gid + 1] += b[2 * rhs_idx + 1];\r
+}\r
 \r
 \r
 /* Kernels for vector addition with a number*/\r
@@ -218,6 +234,20 @@ __kernel void intSubNum(__global int * a, int b, __global int * c)  _M_SUBNUM_KE
 __kernel void uintSubNum(__global unsigned int * a, unsigned int b, __global unsigned int * c)  _M_SUBNUM_KERNEL \r
 __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
+__kernel void complexFloatSubNum(__global float *a, float2 b, __global float *c)\r
+{\r
+    int gid = get_global_id(0) * 2;\r
+\r
+    c[gid] = a[gid] - b.s0;\r
+    c[gid + 1] = a[gid + 1] - b.s1;\r
+}\r
+__kernel void complexDoubleSubNum(__global double *a, double2 b, __global double *c)\r
+{\r
+    int gid = get_global_id(0) * 2;    \r
+    \r
+    c[gid] = a[gid] - b.s0;\r
+    c[gid + 1] = a[gid + 1] - b.s1;\r
+}\r
 \r
 /* Kernels for compound vector subtraction with a number */\r
 \r
@@ -528,12 +558,12 @@ __kernel void floatSqr(__global float * a)  _M_SQR_REAL_KERNEL
 __kernel void doubleSqr(__global double * a)  _M_SQR_REAL_KERNEL\r
 \r
 \r
-#define _M_SQR_COMPLEX_KERNEL                          \\r
-{                                              \\r
-               int gid  = get_global_id(0);            \\r
-        a[2 * gid]   = a[2 * gid] * a[2 * gid]  \\r
-                       - a[2 * gid + 1] * a[2 * gid + 1];\\r
-       a[2 * gid + 1] = 2 * a[2 * gid] * a[2 * gid + 1];\\r
+#define _M_SQR_COMPLEX_KERNEL \\r
+{ \\r
+    int gid  = get_global_id(0); \\r
+    a[2 * gid]     = a[2 * gid] * a[2 * gid] \\r
+                   - a[2 * gid + 1] * a[2 * gid + 1]; \\r
+    a[2 * gid + 1] = 2 * a[2 * gid] * a[2 * gid + 1]; \\r
 }\r
 \r
 __kernel void complexIntSqr(__global int * a)  _M_SQR_COMPLEX_KERNEL\r