2

CUDAでSimpsonの方法による統合をコーディングしようとしています。

これがシンプソンの法則の公式です

ここに画像の説明を入力

どこでx_k = a + k*h。

これが私のコードです

    __device__ void initThreadBounds(int *n_start, int *n_end, int n, 
                                        int totalBlocks, int blockWidth)
    {
        int threadId = blockWidth * blockIdx.x + threadIdx.x;
        int nextThreadId = threadId + 1;

        int threads = blockWidth * totalBlocks;

        *n_start = (threadId * n)/ threads;
        *n_end =  (nextThreadId * n)/ threads;
    }

    __device__ float reg_func (float x)
    {
        return x;
    }

    typedef float (*p_func) (float);

    __device__ p_func integrale_f = reg_func;

    __device__ void integralSimpsonMethod(int totalBlocks, int totalThreads, 
                    double a, double b, int n, float p_function(float), float* result)
    {
        *result = 0;

        float h = (b - a)/n; 
        //*result = p_function(a)+p_function(a + h * n);
        //parallel
        int idx_start;
        int idx_end;
        initThreadBounds(&idx_start, &idx_end, n-1, totalBlocks, totalThreads);
        //parallel_ends
        for (int i = idx_start; i < idx_end; i+=2) {
            *result +=  ( p_function(a + h*(i-1)) + 
                          4 * p_function(a + h*(i)) + 
                          p_function(a + h*(i+1)) ) * h/3;

        }   
    } 


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

        integralSimpsonMethod(totalBlocks, totalThreads, 0, 10, 1000, integrale_f, &res);
        result[(blockIdx.x*totalThreads + threadIdx.x)] = res;

        //printf ("Simpson method\n");
    }


    __host__ void inttest()
    {

        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);
            integralSimpson<<<blocksNum, threadNum>>>(blocksNum, threadNum, 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 ("sum = %f \n", sum);
        cudaFree(device_resultf);
    }

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


   inttest();


    int i;
    scanf ("%d",&i);

}

問題は、が100000nよりも小さい場合に正しく動作しないことです。010~99n = 100000~50

どうしたの？

4

2 に答える 2

5

ここでの基本的な問題は、独自のアルゴリズムを理解していないことです。

関数integralSimpsonMethod()は、各スレッドが積分領域のサブ間隔ごとに少なくとも 3 つの直交点をサンプリングするように設計されています。したがって、カーネル呼び出しのスレッド数の 4 倍未満になるように n を選択すると、各サブ間隔がオーバーラップし、結果の積分が不正確になることは避けられません。コードがスレッド数または n をチェックしてスケーリングし、積分の計算時にオーバーラップが発生しないようにする必要があります。

自己啓発以外の目的でこれを行っている場合は、シンプソンのルールの複合バージョンを調べることをお勧めします. これは並列実装に非常に適しており、正しく実装するとパフォーマンスが大幅に向上します。

于 2013-04-22T10:11:05.480 に答える