title: Element-by-element vector multiplication with CUDA
body: <p>I have build a rudimentary kernel in CUDA to do an <em>elementwise</em> vector-vector multiplication of two complex vectors. The kernel code is inserted below (<code>multiplyElementwise</code>). It works fine, but since I noticed that other seemingly straightforward operations (like scaling a vector) are optimized in libraries like CUBLAS or CULA, I was wondering if it is possible to replace my code by a library call? To my surprise, neither CUBLAS nor CULA have this option, I tried to fake it by making one of the vectors the diagonal of a diagonal matrix-vector product, but the result was really slow.   </p>
  
  <p>As a matter of last resort I tried to optimize this code myself (see <code>multiplyElementwiseFast</code> below) by loading the two vectors in shared memory and then work from there, but that was slower than my original code.</p>
  
  <p>So my questions:</p>
  
  <ol>
  <li>Is there library that does elementwise vector-vector multiplications?</li>
  <li>If not, can I accelerate my code (<code>multiplyElementwise</code>)?</li>
  </ol>
  
  <p>Any help would be greatly appreciated!</p>
  
  <pre><code>__global__ void multiplyElementwise(cufftComplex* f0, cufftComplex* f1, int size)
  {
      const int i = blockIdx.x*blockDim.x + threadIdx.x;
      if (i &lt; size)
      {
          float a, b, c, d;
          a = f0[i].x; 
          b = f0[i].y;
          c = f1[i].x; 
          d = f1[i].y;
          float k;
          k = a * (c + d);
          d =  d * (a + b);
          c =  c * (b - a);
          f0[i].x = k - d;
          f0[i].y = k + c;
      }
  }
  
  __global__ void multiplyElementwiseFast(cufftComplex* f0, cufftComplex* f1, int size)
  {
      const int i = blockIdx.x*blockDim.x + threadIdx.x;
      if (i &lt; 4*size)
      {
          const int N = 256;
          const int thId = threadIdx.x / 4;
          const int rem4 = threadIdx.x % 4;
          const int i4 = i / 4;
  
          __shared__ float a[N];
          __shared__ float b[N];
          __shared__ float c[N];
          __shared__ float d[N];
          __shared__ float Re[N];
          __shared__ float Im[N];
  
          if (rem4 == 0)
          {
              a[thId] = f0[i4].x;
              Re[thId] = 0.f;
          }
          if (rem4 == 1)
          {
              b[thId] = f0[i4].y;
              Im[thId] = 0.f;
          }
          if (rem4 == 2)
              c[thId] = f1[i4].x;
          if (rem4 == 0)
              d[thId] = f1[i4].y;
          __syncthreads();
  
          if (rem4 == 0)
              atomicAdd(&amp;(Re[thId]), a[thId]*c[thId]);        
          if (rem4 == 1)
              atomicAdd(&amp;(Re[thId]), -b[thId]*d[thId]);
          if (rem4 == 2)
              atomicAdd(&amp;(Im[thId]), b[thId]*c[thId]);
          if (rem4 == 3)
              atomicAdd(&amp;(Im[thId]), a[thId]*d[thId]);
          __syncthreads();
  
          if (rem4 == 0)
              f0[i4].x = Re[thId];
          if (rem4 == 1)
              f0[i4].y = Im[thId];
      }
  }        
  </code></pre>