无法在 CUDA 向量乘法中获得最大值

Question

我将 pB 中的每一行多次播放到 pA 中的每一行，并将最大值放入 pC。 问题是：在内部循环中，唯一的最后一行受体被视为“最大值”。结果右栏完全错误。

void TestCalcDotMax_2x5x3()
{
    const size_t m = 2;  // nReceptors
    const size_t k = 5;  // nSources
    const size_t n = 3;  // nChemicals

float pA[m * k] =   { 1, 2, 3, 4, 5
                    , 2, 4, 6, 8, 2};

float pB[k * n] =   { 9, 8, 7, 6, 5
                    , 4, 3, 2, 1, 9
                    , 8, 7, 6, 5, 4 };

float expected[k * n] = { 18, 32, 42, 48, 25
                        , 8, 12, 12,  8, 45
                        ,16, 28, 36, 40, 20 };

float pC[k * n] =   { 18, 32, 42, 48, 10
                    , 8, 12, 12,  8, 18
                    ,16, 28, 36, 40,  8 };

int rst = ::CalcDotMax( pA, pB, m, k, n, pC );

    CPPUNIT_ASSERT_EQUAL_MESSAGE( "passed processing",  0, rst ); 
}
// pDevB and pDevC nave the same size
__global__ void KernelDotMax( const float* pDevA, const float* pDevB, const size_t m, const size_t k, float* pDevC ) 
{ 
    int i = blockDim.x * blockIdx.x + threadIdx.x;

    if( i < m )
    {
        for( size_t j = 0; j < k; j++ )
        {
            const float value = pDevA[ i * k + j ] * pDevB[j];

            if( value > pDevC[j] )
            {
                pDevC[j] = value;
            }
        }
    }
}    

__host__ int CalcDotMax( const float* pA, const float* pB, int m, int k, int n, float* pC, pfnMsg fnMsg )
{
    int nbrCtas = m;
    int threadsPerCta = 64;

    if( nbrCtas >= 32 ) 
    {
         nbrCtas = 32;
         threadsPerCta = 64;
    }
    float* pDevA = nullptr;
    float* pDevB = nullptr;
    float* pDevC = nullptr;

    cudaError_t code = ::cudaMalloc( (void**)&pDevA, m * k * sizeof(float) );

    code = ::cudaMalloc( (void**)&pDevB, k * n * sizeof(float) );
    code = ::cudaMalloc( (void**)&pDevC, k * n * sizeof(float) );

    code = ::cudaMemcpy( pDevA, pA, m * k * sizeof(float), cudaMemcpyHostToDevice); 
    code = ::cudaMemcpy( pDevB, pB, k * n * sizeof(float), cudaMemcpyHostToDevice); 
    code = ::cudaMemcpy( pDevC, pC, k * n * sizeof(float), cudaMemcpyHostToDevice); 

    for( size_t index = 0; index < n * k; index += k )
    {
        KernelDotMax<<<nbrCtas,threadsPerCta>>>( pDevA, &pDevB[index], m, k, &pDevC[index] );
    }
    code = ::cudaMemcpy( pC, pDevC, k * n * sizeof(float), cudaMemcpyDeviceToHost);
    code = ::cudaFree( pDevA );
    code = ::cudaFree( pDevB );
    code = ::cudaFree( pDevC );

    return 0;
}

Answer 1

抱歉，我错过了您编辑代码的某个时间点。

您遇到的问题是竞争条件。在失败的情况下，您将启动 2 个块。您的算法设计使得每个块都在同一组输出元素（在pdevC）上运行。因此，由于两个块可以同时执行，因此两个块可以同时写入相同的输出元素。这是一个碰撞，有两种方法可以避免它：

1. 重新设计您的算法，以不同的方式划分工作 块。而不是每个块检查 all （或同一组）输出元素 针对一组特定的输入，只有每个块 负责输出元素的一部分，但检查 针对所有输入。这是一个常见的代码重构 转换顺序/串行时完成的操作 算法，一种并行运行的算法。
2. 使用atomic operations 防止发生冲突。如果您的算法只有少量此类冲突，则使用原子可能很方便且成本不高。但是，当算法对每个输出元素使用原子（可能多次，如本例所示）时，尝试按照上述方法 1 重构代码可能会更好（以获得更高的性能）。

下面是我说明第二种方法的一些代码（因为它更容易编写）。没有原子函数可以在 float 上提供 atomicMax 操作，所以我按照原子函数文档中给出的模板制作了自己的函数，用于使用 atomicCAS 创建任意原子操作。这就是atomicMaxf。

如果您选择使用第一种方法（推荐），我会指出循环调用内核对于您的算法可能不是必需的。我将制作一个新内核，为每个输出点分配一个线程，然后在内核中的一个循环（或嵌套循环）中计算各种输入点上所有必要的最大操作。由于每个线程都在写入一个且只有一个唯一的输出点，因此线程之间不可能发生写入冲突。

无论如何，这段代码应该提供正确的结果：

#include <stdio.h>

__device__ float atomicMaxf(float* address, float val)
{
    int *address_as_int =(int*)address;
    int old = *address_as_int, assumed;
    while (val > __int_as_float(old)) {
        assumed = old;
        old = atomicCAS(address_as_int, assumed,
                        __float_as_int(val));
        }
    return __int_as_float(old);
}

// pDevB and pDevC have the same size
__global__ void KernelDotMax( const float* pDevA, const float* pDevB, const size_t m, const size_t k, float* pDevC )
{
    int i = blockDim.x * blockIdx.x + threadIdx.x;

    if( i < m )
    {
        for( size_t j = 0; j < k; j++ )
        {
            const float value = pDevA[ i * k + j ] * pDevB[j];

            atomicMaxf(pDevC+j, value);
//            if( value > pDevC[j] )
//            {
//                pDevC[j] = value;
//            }
        }
    }
}


__host__ int CalcDotMax( const float* pA, const float* pB, int m, int k, int n, float* pC )
{
    int nbrCtas = m;
    int threadsPerCta = 64;

    if( nbrCtas >= 32 )
    {
         nbrCtas = 32;
         threadsPerCta = 64;
    }
    float* pDevA = NULL;
    float* pDevB = NULL;
    float* pDevC = NULL;

    cudaError_t code = ::cudaMalloc( (void**)&pDevA, m * k * sizeof(float) );

    code = ::cudaMalloc( (void**)&pDevB, k * n * sizeof(float) );
    code = ::cudaMalloc( (void**)&pDevC, k * n * sizeof(float) );

    code = ::cudaMemcpy( pDevA, pA, m * k * sizeof(float), cudaMemcpyHostToDevice);
    code = ::cudaMemcpy( pDevB, pB, k * n * sizeof(float), cudaMemcpyHostToDevice);
    code = ::cudaMemcpy( pDevC, pC, k * n * sizeof(float), cudaMemcpyHostToDevice);

    for( size_t index = 0; index < n * k; index += k )
    {
        KernelDotMax<<<nbrCtas,threadsPerCta>>>( pDevA, &pDevB[index], m, k, &pDevC[index] );
    }
    code = ::cudaMemcpy( pC, pDevC, k * n * sizeof(float), cudaMemcpyDeviceToHost);
    code = ::cudaFree( pDevA );
    code = ::cudaFree( pDevB );
    code = ::cudaFree( pDevC );

    return 0;
}


void TestCalcDotMax_2x5x3()
{
    const size_t m = 2;  // nReceptors
    const size_t k = 5;  // nSources
    const size_t n = 3;  // nChemicals

  float pA[m * k] =   { 1.0f, 2.0f, 3.0f, 4.0f, 5.0f
                    , 2.0f, 4.0f, 6.0f, 8.0f, 2.0f};

  float pB[k * n] =   { 9.0f, 8.0f, 7.0f, 6.0f, 5.0f
                    , 4.0f, 3.0f, 2.0f, 1.0f, 9.0f
                    , 8.0f, 7.0f, 6.0f, 5.0f, 4.0f };

  float expected[k * n] = { 18.0f, 32.0f, 42.0f, 48.0f, 25.0f
                        , 8.0f, 12.0f, 12.0f,  8.0f, 45.0f
                        ,16.0f, 28.0f, 36.0f, 40.0f, 20.0f };

  float pC[k * n] =   { 0.0f, 0.0f, 0.0f, 0.0f, 0.0f
                    , 0.0f, 0.0f, 0.0f, 0.0f, 0.0f
                    , 0.0f, 0.0f, 0.0f, 0.0f, 0.0f };

  int rst = ::CalcDotMax( pA, pB, m, k, n, pC );

  printf("passed processing: %d \n",  rst );
  for (int i=0; i<(k*n); i++)
    if (pC[i] != expected[i]) printf("mismatch at %d, should be: %f was: %f\n", i, expected[i], pC[i]);
}

int main(){


  TestCalcDotMax_2x5x3();
  return 0;

}

Answer 2

非常感谢 - 现在可以使用了。在比较时是否可以保留迭代索引[idx]？像这样：

struct ValIndex_t
{
    float  value;
    int    index;
};

__device__ float atomicMaxPare( float* address, float val, int* index, int idx )
{
    int *address_as_int = reinterpret_cast<int*>( address->value ); // assume that float has size of integer 32 bit
    int old = *address_as_int, assumed;

    while( val > ::__int_as_float(old) ) 
    {
        assumed = old;
        old     = ::atomicCAS( address_as_int, assumed, ::__float_as_int(val) );
        *index  = idx;
    }
    return ::__int_as_float(old);
}

__global__ void CudaPareDotMax( float* pDevA, const float* pDevB, ValIndex_t* pDevC, const size_t m, const size_t k, const size_t n ) 
{ 
    int idx = blockDim.x * blockIdx.x + threadIdx.x;

    if( idx < m )
    {
        for( size_t row = 0; row < n; row++ )
        {
            for( size_t col = 0; col < k; col++ )
            {
                const size_t slice = col + row * k;
                const size_t index = slice + k * n * idx;

                pDevA[index] *= pDevB[ col + k * idx ];

                float& prvalue = (pDevC +  slice )->value;
                int&   prindex = (pDevC + slice  )->index;

                ::atomicMaxPare( &prvalue, pDevA[ index ], &prindex, idx );
            }
        }
    }
}

或者我必须使用另一个原子函数进行交换？在价值变为最大值的那一刻，不太明白如何加入它。再次感谢