使用CUDA进行逐元素向量乘法

Question

我已经用CUDA构建了一个简单的内核，以执行两个复向量的逐元素向量乘法。内核代码如下（multiplyElementwise）。它能正常工作，但是我注意到像缩放向量这样的其他看似简单的操作在像CUBLAS或CULA这样的库中进行了优化，所以我想知道是否可能用库调用替换我的代码？令我惊讶的是，无论是CUBLAS还是CULA都没有这个选项，我试图通过使其中一个向量成为对角线矩阵-向量乘积的对角线来模拟它，但结果非常慢。

作为最后一招，我尝试通过在共享内存中加载两个向量，然后从那里工作来优化这段代码（请参见下面的multiplyElementwiseFast），但这比我的原始代码更慢。

所以我的问题是：

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];
    }
}        

- WVDB

---

"逐元素向量-向量乘法" 你是指点积吗？ - BenC

@Benc... 不是的。对于真正的向量，点乘积是元素逐个相乘的和。 - sgarizvi

@sgar91：如果他正在“乘法”复数，他可能实际上想要计算一个称为内积/点积的共轭双线性形式（参见此处）。我只是想确保他的意图。 - BenC

Answer 1

如果您想要实现复数的简单逐元素乘积，您似乎在multiplyElementwise内核中执行了一些额外步骤，从而增加了寄存器使用量。您想要计算的是：

f0[i].x = a*c - b*d;
f0[i].y = a*d + b*c;

由于 (a + ib)*(c + id) = (a*c - b*d) + i(a*d + b*c)，通过使用改进的复数乘法，您可以将1次乘法转换为3次加法和一些额外的寄存器。这是否可以证明有效可能取决于您使用的硬件。例如，如果您的硬件支持FMA（融合乘加），那么这种优化可能不是有效的。您应该考虑阅读此文档：“精度和性能：NVIDIA GPU的浮点和IEEE 754兼容性”，该文档还解决了浮点精度问题。

尽管如此，您应该考虑使用Thrust。该库提供了许多高级工具，可在主机和设备向量上操作。您可以在此处查看长列表的示例：https://github.com/thrust/thrust/tree/master/examples。这将使您的生活变得更加轻松。

更新的代码

在您的情况下，您可以使用此示例并将其调整为以下内容：

#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
#include <time.h>

struct ElementWiseProductBasic : public thrust::binary_function<float2,float2,float2>
{
    __host__ __device__
    float2 operator()(const float2& v1, const float2& v2) const
    {
        float2 res;
        res.x = v1.x * v2.x - v1.y * v2.y;
        res.y = v1.x * v2.y + v1.y * v2.x;
        return res;
    }
};

/**
 * See: http://www.embedded.com/design/embedded/4007256/Digital-Signal-Processing-Tricks--Fast-multiplication-of-complex-numbers%5D
 */
struct ElementWiseProductModified : public thrust::binary_function<float2,float2,float2>
{
    __host__ __device__
    float2 operator()(const float2& v1, const float2& v2) const
    {
        float2 res;
        float a, b, c, d, k;
        a = v1.x;
        b = v1.y;
        c = v2.x;
        d = v2.y;
        k = a * (c + d);
        d =  d * (a + b);
        c =  c * (b - a);
        res.x = k -d;
        res.y = k + c;
        return res;
    }
};

int get_random_int(int min, int max)
{
    return min + (rand() % (int)(max - min + 1));
}

thrust::host_vector<float2> init_vector(const size_t N)
{
    thrust::host_vector<float2> temp(N);
    for(size_t i = 0; i < N; i++)
    {
        temp[i].x = get_random_int(0, 10);
        temp[i].y = get_random_int(0, 10);
    }
    return temp;
}

int main(void)
{
    const size_t N = 100000;
    const bool compute_basic_product    = true;
    const bool compute_modified_product = true;

    srand(time(NULL));

    thrust::host_vector<float2>   h_A = init_vector(N);
    thrust::host_vector<float2>   h_B = init_vector(N);
    thrust::device_vector<float2> d_A = h_A;
    thrust::device_vector<float2> d_B = h_B;

    thrust::host_vector<float2> h_result(N);
    thrust::host_vector<float2> h_result_modified(N);

    if (compute_basic_product)
    {
        thrust::device_vector<float2> d_result(N);

        thrust::transform(d_A.begin(), d_A.end(),
                          d_B.begin(), d_result.begin(),
                          ElementWiseProductBasic());
        h_result = d_result;
    }

    if (compute_modified_product)
    {
        thrust::device_vector<float2> d_result_modified(N);

        thrust::transform(d_A.begin(), d_A.end(),
                          d_B.begin(), d_result_modified.begin(),
                          ElementWiseProductModified());
        h_result_modified = d_result_modified;
    }

    std::cout << std::fixed;
    for (size_t i = 0; i < 4; i++)
    {
        float2 a = h_A[i];
        float2 b = h_B[i];

        std::cout << "(" << a.x << "," << a.y << ")";
        std::cout << " * ";
        std::cout << "(" << b.x << "," << b.y << ")";

        if (compute_basic_product)
        {
            float2 prod = h_result[i];
            std::cout << " = ";
            std::cout << "(" << prod.x << "," << prod.y << ")";
        }

        if (compute_modified_product)
        {
            float2 prod_modified = h_result_modified[i];
            std::cout << " = ";
            std::cout << "(" << prod_modified.x << "," << prod_modified.y << ")";
        }
        std::cout << std::endl;
    }   

    return 0;
}

这个命令将返回以下结果：

(6.000000,5.000000)  * (0.000000,1.000000)  = (-5.000000,6.000000)
(3.000000,2.000000)  * (0.000000,4.000000)  = (-8.000000,12.000000)
(2.000000,10.000000) * (10.000000,4.000000) = (-20.000000,108.000000)
(4.000000,8.000000)  * (10.000000,9.000000) = (-32.000000,116.000000)

您可以比较两种不同的乘法策略的时间，并选择适合您的硬件的最佳策略。

- BenC

---

非常感谢，我会研究一下Thrust，看看能否让它为我工作。复杂产品被计算在这样一个“反向”的方式是因为我在这里读到[http://www.embedded.com/design/embedded/4007256/Digital-Signal-Processing-Tricks--Fast-multiplication-of-complex-numbers]，处理器用加法比乘法更容易处理，所以可能会更快一些。不过仍然有点奇怪的是，向量的缩放或两个向量的求和在许多库中都很容易找到，但逐元素乘法却没有...无论如何，感谢您的回答！ - WVDB

@WVDB：我明白了，我想这值得比较，但这种优化真的取决于硬件，并且您需要记住，如果您正在使用编译器优化，则编译器可能会重新排序操作。您可以使用nvprof / nvvp甚至生成的PTX代码来比较计时，以便更好地了解实际情况。 - BenC

@WVDB：我更新了我的答案，使其可以与NVIDIA的性能分析工具进行比较。 - BenC

@ BenC 目前我正在尝试使用 Nsight 对我的代码进行分析。再次感谢您的实用建议！ - WVDB

Answer 2

你可以使用cublasZdgmm。

cublasStatus_t cublasZdgmm(cublasHandle_t handle, cublasSideMode_t mode,
                      int m, int n,
                      const cuDoubleComplex *A, int lda,
                      const cuDoubleComplex *x, int incx,
                      cuDoubleComplex *C, int ldc)

- Phoz