cuda - CUDA を使用した要素ごとのベクトル乗算

Question

CUDA で基本的なカーネルを構築して、2 つの複素ベクトルの要素ごとのベクトル間の乗算を行いました。カーネル コードは下に挿入されます ( multiplyElementwise)。それは正常に動作しますが、他の一見単純な操作 (ベクトルのスケーリングなど) が CUBLAS や CULA などのライブラリで最適化されていることに気付いたので、私のコードをライブラリ呼び出しで置き換えることができるかどうか疑問に思っていました。驚いたことに、CUBLAS にも CULA にもこのオプションはありません。ベクトルの 1 つを対角行列ベクトル積の対角にすることで偽装しようとしましたが、結果は非常に遅くなりました。

multiplyElementwiseFast最後の手段として、2 つのベクトルを共有メモリにロードし、そこから作業することで、このコードを自分で最適化しようとしましたが (以下を参照)、元のコードよりも遅くなりました。

だから私の質問：

1. 要素ごとのベクトル間の乗算を行うライブラリはありますか?
2. そうでない場合、コードを高速化できますか ( multiplyElementwise)?

どんな助けでも大歓迎です！

__global__ void multiplyElementwise(cufftComplex* f0, cufftComplex* f1, int size)
{
    const int i = blockIdx.x*blockDim.x + threadIdx.x;
    if (i < 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 < 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(&(Re[thId]), a[thId]*c[thId]);        
        if (rem4 == 1)
            atomicAdd(&(Re[thId]), -b[thId]*d[thId]);
        if (rem4 == 2)
            atomicAdd(&(Im[thId]), b[thId]*c[thId]);
        if (rem4 == 3)
            atomicAdd(&(Im[thId]), a[thId]*d[thId]);
        __syncthreads();

        if (rem4 == 0)
            f0[i4].x = Re[thId];
        if (rem4 == 1)
            f0[i4].y = Im[thId];
    }
}