1

2D/3D非線形拡散用の最初のCUDAコードを書いています。2Dケースは正常に機能していますが、3Dケースで苦労しています。基本的に、有限差分計算の段階でゼロを取得しています。驚くべきことに、「deltaN」（以下のコードを参照）は正しい答えを示していますが、他のものは機能していません（答えはゼロ）。256x256x256ボリュームを処理しようとしています。何か提案はありますか？ありがとう！

 #define BLKXSIZE 8
 #define BLKYSIZE 8
 #define BLKZSIZE 8
 #define idivup(a, b) ( ((a)%(b) != 0) ? (a)/(b)+1 : (a)/(b) )

void AnisotropDiff_GPU(double* A, double* B, int N, int M, int Z, double sigma, int iter, double tau, int type)
{ 
// Nonlinear Diffusion in 3D 
double *Ad;     

dim3 dimBlock(BLKXSIZE, BLKYSIZE, BLKZSIZE);           
dim3 dimGrid(idivup(N,BLKXSIZE), idivup(M,BLKYSIZE), idivup(Z,BLKYSIZE));    

cudaMalloc((void**)&Ad,N*M*Z*sizeof(double));  
cudaMemcpy(Ad,A,N*M*Z*sizeof(double),cudaMemcpyHostToDevice);  

int n = 1;
while (n <= iter) {    
anis_diff3D<<<dimGrid,dimBlock>>>(Ad, N, M, Z, sigma, iter, tau, type);  
n++;}
cudaMemcpy(B,Ad,N*M*Z*sizeof(double),cudaMemcpyDeviceToHost);
cudaFree(Ad);
}

ここに有限差分を計算するための部分があります

 __global__ void anis_diff3D(double* A, int N, int M, int Z, double sigma, int iter, double tau, int type)
{

 int xIndex = blockDim.x * blockIdx.x + threadIdx.x;
 int yIndex = blockDim.y * blockIdx.y + threadIdx.y;
 int zIndex = blockDim.z * blockIdx.z + threadIdx.z;

 if ( (xIndex < N) && (yIndex < M) && (zIndex < Z) )
 {
    int index_out = xIndex + M*yIndex + N*M*zIndex;

    double deltaN=0, deltaS=0, deltaW=0, deltaE=0, deltaU=0, deltaD=0;
    double cN, cS, cW, cE, cU, cD;

    int indexN = (xIndex-1) + M*yIndex + N*M*zIndex;
    int indexS = (xIndex+1) + M*yIndex + N*M*zIndex;
    int indexW = xIndex + M*(yIndex-1) + N*M*zIndex;
    int indexE = xIndex + M*(yIndex+1) + N*M*zIndex;
    int indexU = xIndex + M*yIndex + N*M*(zIndex-1);
    int indexD = xIndex + M*yIndex + N*M*(zIndex+1);


    if (xIndex>1)
        deltaN = A[indexN]-A[index_out];
    if (xIndex<N)
        deltaS = A[indexS]-A[index_out];    
    if (yIndex>1)
        deltaW = A[indexW]-A[index_out];    
    if (yIndex<M)
        deltaE = A[indexE]-A[index_out];
    if (zIndex>1)
        deltaU = A[indexU]-A[index_out];    
    if (zIndex<Z)
        deltaD = A[indexD]-A[index_out];

   A[index_out] = deltaN ; // works for deltaN but not for deltaS, deltaW... . 

 }

助けてくれてありがとう！

4

1 に答える 1

2

カーネルで計算しようとしている値の一部に対して、範囲外のインデックスが作成されています。

次のように最後のカーネル行を使用してコードをコンパイルする場合：

A[index_out] = deltaS ;

次に、それを実行するとcuda-memcheck、cuda-memcheckは範囲外のアクセスを報告します。

========= Invalid __global__ read of size 8
=========     at 0x000000b0 in anis_diff3D(double*, int, int, int, double, int, double, int)
=========     by thread (7,7,7) in block (31,31,31)
=========     Address 0x408100000 is out of bounds

では、何が起こっているのでしょうか。インデックスの計算を見てみましょう。

int indexS = (xIndex+1) + M*yIndex + N*M*zIndex;

グリッドの最後のスレッド（ブロック（31、31、31）のスレッド（7,7,7））の場合、このインデックス計算は、次の行のメモリ配列の終わりを超えてインデックスを作成します。A

    deltaS = A[indexS]-A[index_out];

物事を適切に実行するには、これらの境界条件を処理する必要があります。

その間、エラーチェックを行っていたとしたら、カーネルはエラーをスローしていたでしょう。カーネルの最後に格納するために選択した値に応じて、コンパイラが他の計算を最適化し、カーネルが正しく実行されているように見える場合があることに注意してください（たとえば、deltaSではなくdeltaNを格納する場合）。エラーチェックを含むコードの例を次に示します。

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

#define cudaCheckErrors(msg) \
    do { \
        cudaError_t __err = cudaGetLastError(); \
        if (__err != cudaSuccess) { \
            fprintf(stderr, "Fatal error: %s (%s at %s:%d)\n", \
                msg, cudaGetErrorString(__err), \
                __FILE__, __LINE__); \
            fprintf(stderr, "*** FAILED - ABORTING\n"); \
            exit(1); \
        } \
    } while (0)

 #define BLKXSIZE 8
 #define BLKYSIZE 8
 #define BLKZSIZE 8
 #define idivup(a, b) ( ((a)%(b) != 0) ? (a)/(b)+1 : (a)/(b) )
 #define sizeX 256
 #define sizeY 256
 #define sizeZ 256
 #define sizeT (sizeX*sizeY*sizeZ)


 __global__ void anis_diff3D(double* A, int N, int M, int Z, double sigma, int iter, double tau, int type)
{

 int xIndex = blockDim.x * blockIdx.x + threadIdx.x;
 int yIndex = blockDim.y * blockIdx.y + threadIdx.y;
 int zIndex = blockDim.z * blockIdx.z + threadIdx.z;

 if ( (xIndex < N) && (yIndex < M) && (zIndex < Z) )
 {
    int index_out = xIndex + M*yIndex + N*M*zIndex;

    double deltaN=0, deltaS=0, deltaW=0, deltaE=0, deltaU=0, deltaD=0;
    double cN, cS, cW, cE, cU, cD;

    int indexN = (xIndex-1) + M*yIndex + N*M*zIndex;
    if (indexN > ((N*M*Z)-1)) indexN = (N*M*Z) -1;
    if (indexN < 0) indexN = 0;
    int indexS = (xIndex+1) + M*yIndex + N*M*zIndex;
    if (indexS > ((N*M*Z)-1)) indexS = (N*M*Z) -1;
    if (indexS < 0) indexS = 0;
    int indexW = xIndex + M*(yIndex-1) + N*M*zIndex;
    if (indexW > ((N*M*Z)-1)) indexW = (N*M*Z) -1;
    if (indexW < 0) indexW = 0;
    int indexE = xIndex + M*(yIndex+1) + N*M*zIndex;
    if (indexE > ((N*M*Z)-1)) indexE = (N*M*Z) -1;
    if (indexE < 0) indexE = 0;
    int indexU = xIndex + M*yIndex + N*M*(zIndex-1);
    if (indexU > ((N*M*Z)-1)) indexU = (N*M*Z) -1;
    if (indexU < 0) indexU = 0;
    int indexD = xIndex + M*yIndex + N*M*(zIndex+1);
    if (indexD > ((N*M*Z)-1)) indexD = (N*M*Z) -1;
    if (indexD < 0) indexD = 0;    

    if (xIndex>1)
        deltaN = A[indexN]-A[index_out];
    if (xIndex<N)
        deltaS = A[indexS]-A[index_out];
    if (yIndex>1)
        deltaW = A[indexW]-A[index_out];
    if (yIndex<M)
        deltaE = A[indexE]-A[index_out];
    if (zIndex>1)
        deltaU = A[indexU]-A[index_out];
    if (zIndex<Z)
        deltaD = A[indexD]-A[index_out];

   A[index_out] = deltaS; // works for deltaN but not for deltaS, deltaW... .

 }
}
void AnisotropDiff_GPU(double* A, double* B, int N, int M, int Z, double sigma, int iter, double tau, int type)
{
  // Nonlinear Diffusion in 3D
  double *Ad;

  dim3 dimBlock(BLKXSIZE, BLKYSIZE, BLKZSIZE);
  dim3 dimGrid(idivup(N,BLKXSIZE), idivup(M,BLKYSIZE), idivup(Z,BLKYSIZE));

  cudaMalloc((void**)&Ad,N*M*Z*sizeof(double));
  cudaCheckErrors("cm1");
  cudaMemcpy(Ad,A,N*M*Z*sizeof(double),cudaMemcpyHostToDevice);
  cudaCheckErrors("cc1");
  int n = 1;
  while (n <= iter) {
    anis_diff3D<<<dimGrid,dimBlock>>>(Ad, N, M, Z, sigma, iter, tau, type);
    n++;
    cudaDeviceSynchronize();
    cudaCheckErrors("kernel");
  }
  cudaMemcpy(B,Ad,N*M*Z*sizeof(double),cudaMemcpyDeviceToHost);
  cudaCheckErrors("cc2");
  cudaFree(Ad);
}

int main(){

  double *A;
  A = (double *)malloc(sizeT * sizeof(double));
  if (A == 0) {printf("malloc fail\n"); return 1;}

  for (int i=0; i< sizeT; i++)
    A[i] = (double)(rand()/(double)RAND_MAX);

  printf("data:\n");
  for (int i = 0; i < 8; i++)
    printf("A[%d] = %f\n", i, A[i]);

  AnisotropDiff_GPU(A, A, sizeX, sizeY, sizeZ, 0.5f, 3, 0.5, 3);
  printf("results:\n");
  for (int i = 0; i < 8; i++)
    printf("A[%d] = %f\n", i, A[i]);

  return 0;
}

編集 上記のコードを編集して、インデックス計算を定義された配列の境界にクランプしました。これにより、範囲外のアクセスを防ぐことができます。アルゴリズムの観点からそれが賢明かどうかはわかりません。

于 2013-03-02T17:00:49.973 に答える