无穷大导致双重操作答案

【问题标题】：infinity as result in double operation无穷大导致双重操作
【发布时间】：2013-09-16 13:35:42
【问题描述】：

我会理解为什么结果是无穷大的。我写了下面的代码，我总是收到 inf 作为结果。我的代码有任何精度问题吗？

#include <stdio.h>
#include <stdlib.h>

#include "cuda.h"
#include "curand_kernel.h"

#define NDIM 30
#define NPAR 5

#define DIMPAR NDIM*NPAR

__device__ double uniform(int index){
    return (double) 0.767341;
}


__global__ void iteracao(double *pos){

    int thread = threadIdx.x + blockDim.x * blockIdx.x;
    double tvel;
    int i = 0;

    double l, r, t;

    if(thread < DIMPAR){
        do{
            t = (double) uniform(thread);
            l = (double) 2.05 * t * ( pos[thread] );
            r = (double) 2.05 * t * ( pos[thread] );
            tvel = (double) l+t+r;
            pos[thread] =  tvel;
            i++;
        }while(i < 10000);
    }

}


int main(int argc, char *argv[])
{

    double *d_pos,    *h_pos;


    h_pos = (double *) malloc(sizeof( double ) * DIMPAR);


    cudaMalloc((void**)&d_pos, DIMPAR   * sizeof( double ));


    int i, j, k, numthreadsperblock, numblocks;

    numthreadsperblock = 512;
    numblocks = (DIMPAR / numthreadsperblock) + ((DIMPAR % numthreadsperblock)?1:0);
    //
    printf("numthreadsperblock: %i;; numblocks:%i\n", numthreadsperblock, numblocks);

    cudaMemset(d_pos,  0.767341, DIMPAR   * sizeof( double ));
    iteracao<<<numblocks,numthreadsperblock>>>(d_pos);
    cudaMemcpy(h_pos, d_pos, DIMPAR * sizeof( double ), cudaMemcpyDeviceToHost);

    printf("\n");
    for(i = 0; i < NPAR; i++){
        for(j = i*NDIM, k = j; j < (k+30); j++){
            printf("%f,", h_pos[j]);
        }
        printf("***\n\n");
    }

    system("PAUSE");
    return 0;
}

输出总是这样：

numthreadsperblock: 512;;块数：1

inf,inf,inf,inf,inf,inf,inf,inf,inf,inf,inf,inf,inf,inf,inf,inf,inf,inf,inf,inf,inf,inf,inf,inf ,inf,inf,inf,inf,inf,inf,*

【问题讨论】：

您很可能将 d_pos 设置为垃圾。 cudaMemset 对字节值进行操作，并且您给出的浮点数很可能仅导致零字节。

标签： cuda double-precision

【解决方案1】：

你有两个问题。第一个是 @Anycorn 在 cmets 中描述的。 cudaMemset，就像memset 需要一个字节值并设置字节位置。您不能使用它来初始化float 值。

第二个是你的内核有一个循环，它在每个 pos 数组元素上运行 10000 次。实际上，您正在寻找复杂表达式的 10000 阶乘。由于该表达式始终是肯定的，因此您的答案会爆炸。您的内核很可能没有正确编写。它没有做你想让它做的事情。即使您解决了第一个问题并将pos 正确初始化为零，您的计算仍然会失败。

你正在执行的算术是：

pos[idx] =  0.767341 + (3.1460981 * pos[idx]);

对于每个idx，您正在执行上述操作 10000 次。即使初始 pos[idx] 值为零，在循环的第二次迭代中，它也会开始以几何方式起飞。

【讨论】：

【解决方案2】：

您以错误的方式初始化d_pos。 cudaMemset() 只能逐字节设置内存。见cudaMemset() doc 了解更多详情。

要按照您的意图初始化数组，您可以使用 Thrust 作为一种快捷方式。

thrust::fill(
    thrust::device_pointer_cast(d_pos),
    thrust::device_pointer_cast(d_pos) + DIMPAR,
    0.767341);

【讨论】：