Форум программистов [Powered by Invision Power Board]

Версия для печати темы
Нажмите сюда для просмотра этой темы в оригинальном формате

Форум программистов > C/C++: Для новичков > CUDA. Вычисление интеграла

Автор: ioManip 8.11.2015, 14:49

Здравствуйте! Есть такая проблема, что при вычислении определенного интеграла на отрезке [0;1] (

Код

f(x) = 4 / (1 + x * x)

), при разном количестве прямоугольников, результат не всегда является верным.
Например при n = 10000, выдает 3,14159 и т.д
а при n = 1000000, выдает 0,20384.

Я не пойму почему так происходит.
Я использую GTS 250 1Gb with CC 1.1 и Ubuntu 14.04 x64

Вот код
main.cpp

Код


#include <iostream>
#include <iomanip>

#include "kernel.h"

#include <cuda_runtime.h>

int main()
{
   const long stepCount = 10000;

   float *hArray = new float[stepCount];
   float h = 1.f / static_cast<float>(stepCount);

   for (size_t i = 0; i < stepCount; i++)
   {
      hArray[i] = (i + 0.5f) * h;
   }

   float *dArray = nullptr;

   cudaMalloc((void**)&dArray, stepCount * sizeof(float));

   cudaEvent_t start, stop;
   cudaEventCreate(&start);
   cudaEventCreate(&stop);

   cudaMemcpy(dArray, hArray, stepCount * sizeof(float), cudaMemcpyKind::cudaMemcpyHostToDevice);
   cudaEventRecord(start);
   calcIntegral(dArray, stepCount);
   cudaEventRecord(stop);

   cudaEventSynchronize(stop);

   cudaMemcpy(hArray, dArray, stepCount * sizeof(float), cudaMemcpyKind::cudaMemcpyDeviceToHost);

   float sum = 0.f;
   for (size_t i = 0; i < stepCount; i++)
   {
      sum += hArray[i];
   }

   float ms = 1.f;
   cudaEventElapsedTime(&ms, start, stop);
   std::cout << "Time: " << ms << " ms" << std::endl;
   std::cout << "Result: " << std::setprecision(7) << std::fixed << sum / stepCount << std::endl;

   delete[] hArray;

   cudaFree(dArray);

   return 0;
}

kernel.h

Код


#ifndef KERNEL_H
#define KERNEL_H

void calcIntegral(float *array, const long stepCount);

#endif // KERNEL_H

kernel.cu

Код


#include "kernel.h"

#include <cuda_runtime.h>
#include <device_launch_parameters.h>


__device__ float f(float x)
{
   return 4.f / (1.f + x * x);
}

__global__ void calcIntegralGPU(float *array, const long stepCount)
{
   int idx = blockIdx.x * blockDim.x + threadIdx.x;

   float x = array[idx];
   if (idx < stepCount)
   {
      array[idx] = f(x);
   }
}

void calcIntegral(float *array, const long stepCount)
{
   int threadPerBlocks = 256;
   int blocksPerGrid = 256;

   calcIntegralGPU<<< blocksPerGrid , threadPerBlocks >>>(array, stepCount);
}

Автор: math64 9.11.2015, 09:42

Скорее всего, при сложении кучи маленьких чисел, накапливается ошибка округления.
Будьте аккуратнее!
(a+b) + c != a + (b + c)

Автор: tzirechnoy 9.11.2015, 21:05

Цитата

if (idx < stepCount) {

заменить на

Код

 while (idx < stepCount) {
...
 idx += blockDim.x*gridDim.x;
}

А то у вас только максимум первые 65536 значений вычисляются (blockDim * gridDim). Потому и, кстати, правильно всё получается если меньшэ 65536 шагов.

Ну, и по поводу точности тожэ правильное замечание -- если складывать от бОльшых значений к мЕньшым -- то за 1000000 итэрацый можэт, при некотом невезении, до 1/17 от результирующего значения ошыбка накопиться. В общем, чтобы не париться, складывайте хотя бы в тип double. К тому жэ, сложэние тожэ можно распараллелить -- ну, в смысле, в рамках одного потока провести на GPU. Чтобы в итоге подсложыть только до 64k чисел. Но делайте через double всё равно.

Автор: ioManip 11.11.2015, 14:58

Цитата(tzirechnoy @ 9.11.2015, 21:05

)

Но делайте через double всё равно.

У меня аппаратно не поддерживается данный тип. Карта старая gts250 на g92. На сколько я знаю, если даже double и написать, то все равно к float приводить будет.

Цитата(tzirechnoy @ 9.11.2015, 21:05

)

заменить на

Ок, попробую

Автор: tzirechnoy 11.11.2015, 22:49

Так ты ведь не на карте суммируешь.

Автор: ioManip 13.11.2015, 13:02

tzirechnoy, Да, все верно, просто бегло прочитал и не сразу понял о чем речь.

Попробовал использовать Ваши изменения и вот что вышло

Просто поменял код ядра
n = 10000 --> result: 3.1415908
n = 1000000 --> result: 3.9995644

Поменял код ядра и тип данных на double
n = 10000 --> result: 0.4370185
n = 1000000 --> result: 1.1864045