增加数组大小时 CUDA nbody 平铺计算错误代码 77

Question

我无法在 CUDA 代码中解决此问题。

我基本上是在计算 gems3 中的 nbody 问题，并增加数组大小。

粒子在__global__ void ParticleAmplification()内核中的特定数组dev_Ionisation中创建，并通过主机函数DynamicAlloc()添加到全局数组中。

在这一个中，空位置被删除，新的位置被放在新数组的末尾。 由于我抛出的线程多于可用粒子，因此我有一个转义变量以避免浪费时间检查是否存在粒子。

块和瓦片的数量是动态分配的，设备数组是通过以下方式重新分配的：

checkCudaErrors( cudaFree( dev_vector ) );
checkCudaErrors( cudaMalloc( (void**)&dev_vector, N * sizeof(ParticleProperties) ) );

然后经过一些步骤，通常当粒子数增加到 28000 左右时，内核就会崩溃。它给了我错误代码 77，这似乎归因于 __device__ float3 computeBodyAccel 函数中外部共享变量大小 extern __shared__ float3 sharedPos[] 的错误大小 (cudaDeviceSynchronize() error code 77: cudaErrorIllegalAddress)。但是，那个似乎总是以相同的大小正确传递给内核：

size_t sharedMemSize = ThreadsInit * sizeof(float3);

integrateBodies<<<blocksInit, ThreadsInit, sharedMemSize>>>( dev_vectorIonisation, dt, numTiles, nbodyTemp );

当使用一个固定但很大的数组时，一切都很好。

我做错了什么？共享内存是否会因未释放的内存而变满？

这是完整的可编译代码：

// ----- Libraries C/C++ ----- //

#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <iomanip>
#include <string>
#include <math.h>
#include <time.h>
#include <dirent.h>

// ----- Libraries CUDA ----- //

#include "cuda.h"
#include <helper_cuda.h>

#include "curand_kernel.h"

// ----- Global variables ----- //

#define El_DIM 512

#define imin(a,b) (a<b?a:b)

using namespace std;

__constant__ float softening_ = 1.0e-12;    // softening factor for nbody interaction

__device__ __managed__ int NewParticles = 0;
__device__ __managed__ int TotalProcesses = 0;
__device__ __managed__ bool Ampl = false;

const int ThreadsInit = 512;
const int blocksPerGrid = (int)( ( El_DIM * El_DIM + ThreadsInit -1 ) / ThreadsInit );


struct ParticleProperties{
    float3 Position, Velocity, Force;
};


__device__ void initVector( ParticleProperties *dev_Vect, int index ){

    dev_Vect[index].Position.x = -1.0;
    dev_Vect[index].Position.y = -1.0;
    dev_Vect[index].Position.z = -1.0;

    dev_Vect[index].Velocity.x = 0.0;
    dev_Vect[index].Velocity.y = 0.0;
    dev_Vect[index].Velocity.z = 0.0;

    dev_Vect[index].Force.x = 0.0;
    dev_Vect[index].Force.y = 0.0;
    dev_Vect[index].Force.z = 0.0;
}

__device__ void SetVector( ParticleProperties *dev_Vect, float3 position, float4 v, int index ){

    dev_Vect[index].Position.x = position.x;
    dev_Vect[index].Position.y = position.y;
    dev_Vect[index].Position.z = position.z;

    dev_Vect[index].Velocity.x = v.x;
    dev_Vect[index].Velocity.y = v.y;
    dev_Vect[index].Velocity.z = v.z;

    dev_Vect[index].Force.x = 0.0;
    dev_Vect[index].Force.y = 0.0;
    dev_Vect[index].Force.z = 0.0;  
}


__device__ float3 bodyBodyInteraction( float3 fi, float3 bi, float3 bj ){

    float3 r;

    // r_ij  [4 FLOPS]
    r.x = ( bj.x - bi.x );
    r.y = ( bj.y - bi.y );
    r.z = ( bj.z - bi.z );
    r.z = 0.0;

    // distSqr = dot(r_ij, r_ij) + EPS^2  [7 FLOPS]
    float distSqr = r.x * r.x + ( r.y * r.y + ( r.z * r.z + softening_ * softening_ ) );

    // invDistCube =1/distSqr^(3/2)  [4 FLOPS (2 mul, 1 sqrt, 1 inv)]
    float invDist = rsqrt(distSqr);
    float invDistCube =  invDist * invDist * invDist;

    // s = m_j * invDistCube [2 FLOP]
    float s = invDistCube;
    // a_i =  a_i + s * r_ij [6 FLOPS]
    fi.x += r.x * s;
    fi.y += r.y * s;
    fi.z += r.z * s;

    return fi;
}


__device__ float3 computeBodyAccel( float3 force, float3 bodyPos, ParticleProperties * positions, const int numTiles, const int nbody ){

    extern __shared__ float3 sharedPos[];

    int computedNbody = 0;

    for( int tile = 0; tile < numTiles; tile++ ){

        sharedPos[threadIdx.x] = positions[tile * blockDim.x + threadIdx.x].Position;

        __syncthreads();

        // This is the "tile_calculation" from the GPUG3 article.

#pragma unroll 128

        for( unsigned int counter = 0; counter < blockDim.x; counter++ ){
            force = bodyBodyInteraction(force, bodyPos, sharedPos[counter]);
            computedNbody++;
            if( computedNbody == nbody ) break;
        }
        __syncthreads();
    }

    return force;
}

__global__ void integrateBodies( ParticleProperties * __restrict__ dev_vector, float deltaTime, int numTiles, int nbody ){

    int index = blockIdx.x * blockDim.x + threadIdx.x;

    float3 position = {0.0, 0.0, 0.0};
    float3 force = {0.0, 0.0, 0.0};

    if( index < nbody ){
        position = dev_vector[index].Position;

        force = computeBodyAccel( force, position, dev_vector, numTiles, nbody );

        // store new force
        dev_vector[index].Position = position;

        dev_vector[index].Force = force;
    }
}


__global__ void IntegrationKernel( ParticleProperties * __restrict__ dev_vector, const float deltaT, const int nbody ){

    int tid = blockIdx.x * blockDim.x + threadIdx.x;

    float3 dvel;
    float3 velocity;

    if( tid < nbody ){
        // integrate
        dvel.x = dev_vector[tid].Force.x * deltaT * 0.5;
        dvel.y = dev_vector[tid].Force.y * deltaT * 0.5;
        dvel.z = dev_vector[tid].Force.z * deltaT * 0.5;

        velocity.x = dev_vector[tid].Velocity.x + dvel.x;
        velocity.y = dev_vector[tid].Velocity.y + dvel.y;
        velocity.z = dev_vector[tid].Velocity.z + dvel.z;

        dev_vector[tid].Position.x += velocity.x * deltaT;
        dev_vector[tid].Position.y += velocity.y * deltaT;
        dev_vector[tid].Position.z += velocity.z * deltaT;

        dev_vector[tid].Velocity.x = velocity.x + dvel.x;
        dev_vector[tid].Velocity.y = velocity.y + dvel.y;
        dev_vector[tid].Velocity.z = velocity.z + dvel.z;
    }
}


__global__ void ParticleAmplification( curandState *state, ParticleProperties * __restrict__ dev_vectorIonisation, 
                                        ParticleProperties * __restrict__ dev_Ionisation, 
                                        const float dt, int numbodies ){

    int tid = threadIdx.x + blockIdx.x * blockDim.x;

    int LocalProcesses = 0;

    float3 position = {0.0, 0.0, 0.0};

    float4 v_new = {0.0, 0.0, 0.0, 0.0};

    float prob = 0.0;

    if( TotalProcesses >= El_DIM * El_DIM - 1 ) Ampl = false;


    if( tid < numbodies ){
        position.x = dev_vectorIonisation[tid].Position.x;
        position.y = dev_vectorIonisation[tid].Position.y;
        position.z = dev_vectorIonisation[tid].Position.z;

        prob = curand_uniform( &state[tid] );

        if( Ampl ){
            if( prob < 1.e-3 ){
                atomicAdd( &TotalProcesses, 1 );
                LocalProcesses = atomicAdd( &NewParticles, 1 );

                v_new.x = 0.0;
                v_new.y = 0.0;
                v_new.z = 0.0;          

                SetVector( dev_Ionisation, position, v_new, LocalProcesses );
            }
        }
    }
}


__global__ void initCurand( curandState *state, unsigned long seed ){
    int tid = threadIdx.x + blockIdx.x * blockDim.x;

    curand_init(seed, tid, 0, &state[tid]);
}


__global__ void initProcessIoni( ParticleProperties *dev_Vect ){
    int x = threadIdx.x + blockIdx.x * blockDim.x;
    initVector( dev_Vect, x );
}


__global__ void Enumerate_Nbody( ParticleProperties *dev_Vect, int *N, int PrevNbody ){
    int tid = threadIdx.x + blockIdx.x * blockDim.x;

    int gid = blockIdx.x;

    extern __shared__ int cache[];

    if( tid == 0 ) 
        *N = 0;

    if( threadIdx.x == 0 )  
        cache[gid] = 0;

    __syncthreads();

    while( tid < PrevNbody ){
        if( dev_Vect[tid].Position.x > -1.0 )
            atomicAdd( &(cache[gid]), 1 );
        tid += blockDim.x * gridDim.x;
    }

    __syncthreads();

    if( threadIdx.x == 0 )
        atomicAdd( N, cache[gid] );
}



void DynamicAlloc( ParticleProperties **DynamicVector, const ParticleProperties *StaticVector, const int n, int nbody, const int max ){

    ParticleProperties *h_vectorIonisation = new ParticleProperties [nbody];

    ParticleProperties *VectTemporary = new ParticleProperties [n]; 


    checkCudaErrors( cudaMemcpy( VectTemporary, *DynamicVector, n * sizeof(ParticleProperties), cudaMemcpyDeviceToHost ) );

    checkCudaErrors( cudaFree( *DynamicVector ) );

    int i = 0;
    int j = 0;

        for( i = 0; i < n; i++ ){
            if( VectTemporary[i].Position.x > -1.0 ){
                h_vectorIonisation[j] = VectTemporary[i];
                j++;
            }
        }

    delete [] VectTemporary;

    if( NewParticles != 0 ){
        ParticleProperties *StaticVectTemporary = new ParticleProperties [max]; 
        checkCudaErrors( cudaMemcpy( StaticVectTemporary, StaticVector, max * sizeof(ParticleProperties), cudaMemcpyDeviceToHost ) );
        int k = 0;

#pragma unroll 32       
        for( i = 0; i < max; i++ ){
            if( StaticVectTemporary[i].Position.x > -1.0 ){
                h_vectorIonisation[j + k] = StaticVectTemporary[i];
                k++;
            }
        }

        delete [] StaticVectTemporary;
    }

    if( nbody > 0 ){
        checkCudaErrors( cudaMalloc( (void**)DynamicVector, nbody * sizeof(ParticleProperties) ) );

        checkCudaErrors( cudaMemcpy( *DynamicVector, h_vectorIonisation, nbody * sizeof(ParticleProperties), cudaMemcpyHostToDevice ) );
    }

    delete [] h_vectorIonisation;
}


int main( int argc_, char **argv_ ){    

    cudaDeviceReset();  

    cudaDeviceProp prop;

    checkCudaErrors( cudaGetDeviceProperties( &prop, 0 ) );

    int Newers = 256;

    int nbody = 1;

    Ampl = true;

    int *dev_nbody;
    checkCudaErrors( cudaMalloc( (void**)&dev_nbody, sizeof(int) ) );
    checkCudaErrors( cudaMemcpy( dev_nbody, &nbody, sizeof(int), cudaMemcpyHostToDevice ) );

    float dt = 0.5e-13;

    float3 pos;
    pos.x = 1.0 / 2.0 * 1.0e-3;
    pos.y = 1.0 / 2.0 * 1.0e-3;
    pos.z = 1.0 / 2.0 * 1.0e-3;

    float3 speed;
    speed.x = 0.0;
    speed.y = 0.0;
    speed.z = 0.0;

    ParticleProperties *dev_vectorIonisation;
    checkCudaErrors( cudaMalloc( (void**)&dev_vectorIonisation, nbody * sizeof(ParticleProperties) ) );

    ParticleProperties *host_vectorIonisation = new ParticleProperties [nbody];

    clog << "Particles array initialisation...";

    for( int i = 0; i < nbody; i++ ){
        host_vectorIonisation[i].Position.x = drand48() * 1.0e-6 + pos.x;
        host_vectorIonisation[i].Position.y = drand48() * 1.0e-6 + pos.y;
        host_vectorIonisation[i].Position.z = 0.0;

        host_vectorIonisation[i].Velocity.x = speed.x;
        host_vectorIonisation[i].Velocity.y = speed.y;
        host_vectorIonisation[i].Velocity.z = speed.z;

        host_vectorIonisation[i].Force.x = 0.0;
        host_vectorIonisation[i].Force.y = 0.0;
        host_vectorIonisation[i].Force.z = 0.0;
    }

    checkCudaErrors( cudaMemcpy( dev_vectorIonisation, host_vectorIonisation, nbody * sizeof(ParticleProperties), cudaMemcpyHostToDevice ) );

    delete [] host_vectorIonisation;

    clog << "Done" << endl;

    ParticleProperties *dev_Ionisation;
    checkCudaErrors( cudaMalloc( (void**)&dev_Ionisation, Newers * sizeof(ParticleProperties) ) );  

    curandState *RndState;  
    checkCudaErrors( cudaMalloc( (void**)&RndState, El_DIM * El_DIM * sizeof(curandState) ) );

    unsigned long seed = 1773;  

    clog << "cuRand array initialisation...";

    initCurand<<<blocksPerGrid, ThreadsInit>>>( RndState, seed );

    initProcessIoni<<<1, Newers>>>( dev_Ionisation );

    clog << "Done" << endl;

    clog << "Propagation of " << nbody << " primary particle(s)." << endl;

    int ProcessTemp = 0; 

    int nbodyTemp = nbody;

    int blocksInit = (nbody + ThreadsInit - 1) / ThreadsInit;
    int numTiles = (nbody + ThreadsInit - 1) / ThreadsInit;

    size_t sharedMemSize = ThreadsInit * sizeof(float3);

    char buffer[64];

    setvbuf(stdout, buffer, _IOFBF, sizeof(buffer));


    while( nbody > 0 ){
        integrateBodies<<<blocksInit, ThreadsInit, sharedMemSize>>>( dev_vectorIonisation, dt, numTiles, nbodyTemp );

        IntegrationKernel<<<blocksInit, ThreadsInit>>>( dev_vectorIonisation, dt, nbodyTemp );

        ParticleAmplification<<<blocksInit, ThreadsInit>>>( RndState, dev_vectorIonisation, dev_Ionisation, dt, nbodyTemp );
        checkCudaErrors( cudaDeviceSynchronize() );

        Enumerate_Nbody<<<blocksInit, ThreadsInit, blocksInit * sizeof(int)>>>( dev_vectorIonisation, dev_nbody, nbodyTemp );
        checkCudaErrors( cudaDeviceSynchronize() );
        getLastCudaError("Kernel enumerate bodies execution failed");

        checkCudaErrors( cudaMemcpy( &nbody, dev_nbody, sizeof(int), cudaMemcpyDeviceToHost ) );

        nbody += NewParticles;

        if( NewParticles > ProcessTemp ) ProcessTemp = NewParticles;

        if( nbody != nbodyTemp ){
            DynamicAlloc( &dev_vectorIonisation, dev_Ionisation, nbodyTemp, nbody, Newers );

            numTiles = blocksInit = ( nbody + ThreadsInit - 1) / ThreadsInit;

            if( NewParticles != 0 ){
                initProcessIoni<<<1, Newers>>>( dev_Ionisation );
                checkCudaErrors( cudaDeviceSynchronize() );
            }

            nbodyTemp = nbody;

            NewParticles = 0;
            checkCudaErrors( cudaDeviceSynchronize() );
        }
        printf("\r nbodies: %d", nbodyTemp);
    }

    checkCudaErrors( cudaFree( dev_Ionisation ) );
}

这是在计算能力为 3.5 的 GTX Titan Black 上执行的

Answer 1

你的问题从这行代码开始（main）：

numTiles = blocksInit = ( nbody + ThreadsInit - 1) / ThreadsInit;

这会创建足够的图块以完全覆盖nbody 的大小，但并非每个图块都完全填充了主体。

然后问题实际上在您的 computeBodyAccel 例程中表现出来，从 integrateBodies 调用：

for( int tile = 0; tile < numTiles; tile++ ){

    sharedPos[threadIdx.x] = positions[tile * blockDim.x + threadIdx.x].Position;

您对positions 的索引没有任何保护，并且您假设每个磁贴对于threadIdx.x 的每个值都有一个有效的positions 条目。但事实并非如此，通过使用-lineinfo 编译您的代码并使用cuda-memcheck 运行它可以发现问题的第一个表现。在这种情况下，由于cuda-memcheck 提供了严格的内存保护，您的代码（对我而言）在大约 500 个主体处失败，而不是 28000 个。具体的失败是在最后一行代码指示的大小为 4 的无效全局读取更多。 （因此，这不是与共享内存的写入相关的索引问题。）从根本上说，问题是tile*blockDim.s + threadIdx.x 可以超过nbody，并且您在读取时索引越界的positions。 （使用-lineinfo 来识别失败的特定内核代码行已涵盖here）

对computeBodyAccel 例程的以下限制检查修改允许我将您的代码运行到最多约 262,000 个主体，在那里它停止增加（由于Ampl 限制为El_DIM*El_DIM）并停留在那里：

__device__ float3 computeBodyAccel( float3 force, float3 bodyPos, ParticleProperties * positions, const int numTiles, const int nbody ){

    extern __shared__ float3 sharedPos[];

    int computedNbody = 0;

    for( int tile = 0; tile < numTiles; tile++ ){
        if ((tile*blockDim.x + threadIdx.x) < nbody)

          sharedPos[threadIdx.x] = positions[tile * blockDim.x + threadIdx.x].Position;

        __syncthreads();

        // This is the "tile_calculation" from the GPUG3 article.

        int limit = blockDim.x;
        if (tile = (numTiles - 1)) limit -= (numTiles*blockDim.x)-nbody;
#pragma unroll 128

        for( unsigned int counter = 0; counter < limit; counter++ ){
            force = bodyBodyInteraction(force, bodyPos, sharedPos[counter]);
            computedNbody++;
            if( computedNbody == nbody ) break;
        }
        __syncthreads();
    }

    return force;
}

您的代码中似乎还存在其他问题，即使通过上述修复它似乎可以运行。如果您使用cuda-memcheck 并使用here 描述的-lineinfo 方法运行代码（使用上面的“修复”），您会发现（由于更严格的内存范围检查cuda-memcheck 强制执行）当数字物体的数量变大，最终你在ParticleAmplification 尝试创建一个新粒子并在最后调用SetVector 时遇到另一个内存访问错误。看来您在此行之间存在竞争条件：

if( TotalProcesses >= El_DIM * El_DIM - 1 ) Ampl = false;

以及以下可能同时增加TotalProcesses 和LocalProcesses 的行：

            atomicAdd( &TotalProcesses, 1 );
            LocalProcesses = atomicAdd( &NewParticles, 1 );

由于您有许多并行运行的线程，这种限制检查是没有用的。您需要通过检查 atomicAdd 操作的实际返回值并查看它们是否超出限制来更仔细地管理新粒子的建立。