Bessel function on CUDA to within some epsilon

Question

I try to execute Bessel functions( J0(x) as example) using CUDA. Heres the formula

I try to get result to within some epsilon value. So here's the code

    __device__ void Bessel_j0(int totalBlocks, int totalThreads, float z, float epsilon, float* result){
            int n = 1;
            *result = 0;

            bool epsilonFlag = true;

            int idx_start;
            int idx_end;

            while(epsilonFlag == true){ 
                initThreadBounds(&idx_start, &idx_end, n, totalBlocks, totalThreads);
                float a_k;
                for (int k = idx_start; k < idx_end; k++) {
                    a_k = m_power((-0.25 * z * z), k)/(m_factorial(k) * m_factorial(k)); 
                    *result += a_k;
                }
                if(a_k < epsilon){
                        epsilonFlag = false;
                }
                n++;
            }
        }

__global__ void J0(int totalBlocks, int totalThreads,  float x, float* result){
        float res = 0;

        Bessel_j0(totalBlocks, totalThreads, 10, 0.01, &res);
        result[(blockIdx.x*totalThreads + threadIdx.x)] = res;
}

__host__ void J0test(){

    const int blocksNum = 32;
    const int threadNum = 32;

    float   *device_resultf; //для устройства
    float   host_resultf[threadNum*blocksNum]   ={0};


    cudaMalloc((void**) &device_resultf, sizeof(float)*threadNum*blocksNum);

    J0<<<blocksNum, threadNum>>>(blocksNum, threadNum, 10, device_resultf); 
    cudaThreadSynchronize();

    cudaMemcpy(host_resultf, device_resultf, sizeof(float)*threadNum*blocksNum, cudaMemcpyDeviceToHost);

    float sum = 0;

    for (int i = 0; i != blocksNum*threadNum; ++i) {
        sum += host_resultf[i];
        printf ("result in %i cell = %f \n", i, host_resultf[i]);
    }
    printf ("Bessel res = %f \n", sum);
    cudaFree(device_resultf);
}
int main(int argc, char* argv[])
{
    J0test();   
}

When I run it appears black screen and Windows says that nVidia driver didn't respond and it recovered it. And in console output there are only zeros in host_resultf array. What's wrong? How can I execute properly functions to within some epsilon?

Answer 1

It is unlikely, but may be your kernel execution hit the allowable kernel execution time limit. Your code doesn't show an upper limit for iteration numbers. It can happen that epsilon is never reached and your kernel keeps executing beyond the time limits. This site can help.

In all cases, I will add an upper limit to the epsilon loop, never leave a code to run without limit on the iterations number.