Я пишу программу численного интегрирования, которая реализует правило трапеций с адаптивным размером шага. Не вдаваясь в подробности, алгоритм использует рекурсию для вычисления интеграла закодированной математической функции в заданном интервале с заданным относительным допуском. Я упростил код для публикации, но сохранил все основные моменты, поэтому некоторые части могут показаться ненужными или слишком сложными. Вот:
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <cmath>
#include <iostream>
#include <iomanip>
class Integral {
public:
double value; // the value of the integral
__device__ __host__ Integral() : value(0) {};
__device__ __host__ Integral& operator+=(Integral &I);
};
__device__ Integral trapezoid(double a, double b, double tolerance, double fa, double fb);
__device__ __host__ double f(double x); // the integrand function
const int BLOCKS = 1;
const int THREADS = 1;
__global__ void controller(Integral *blockIntegrals, double a, double b, double tolerance) {
extern __shared__ Integral threadIntegrals[]; // an array of thread-local integrals
double fa = f(a), fb = f(b);
threadIntegrals[threadIdx.x] += trapezoid(a, b, tolerance, fa, fb);
blockIntegrals[blockIdx.x] += threadIntegrals[0];
}
int main() {
// *************** Input parameters ***************
double a = 1, b = 800; // integration bounds
double tolerance = 1e-7;
// ************************************************
cudaError cudaStatus;
Integral blockIntegrals[BLOCKS]; // an array of total integrals computed by each block
Integral *devBlockIntegrals;
cudaStatus = cudaMalloc((void**)&devBlockIntegrals, BLOCKS * sizeof(Integral));
if (cudaStatus != cudaSuccess)
std::cout << "cudaMalloc failed!\n";
double estimate = 0; // a rough 10-point estimate of the whole integral
double h = (b - a) / 10;
for (int i = 0; i < 10; i++)
estimate += f(a + i*h);
estimate *= h;
tolerance *= estimate; // compute relative tolerance
controller<<<BLOCKS, THREADS, THREADS*sizeof(Integral)>>>(devBlockIntegrals, a, b, tolerance);
cudaStatus = cudaGetLastError();
if (cudaStatus != cudaSuccess)
std::cout << "addKernel launch failed: " << cudaGetErrorString(cudaStatus) << "\n";
cudaStatus = cudaMemcpy(blockIntegrals, devBlockIntegrals, BLOCKS * sizeof(Integral), cudaMemcpyDeviceToHost);
if (cudaStatus != cudaSuccess)
std::cout << "cudaMemcpy failed: " << cudaGetErrorString(cudaStatus) << "\n";
Integral result; // final result
for (int i = 0; i < BLOCKS; i++) // final reduction that sums the results of all blocks
result += blockIntegrals[i];
std::cout << "Integral = " << std::setprecision(15) << result.value;
cudaFree(devBlockIntegrals);
getchar();
return 0;
}
__device__ double f(double x) {
return log(x);
}
__device__ Integral trapezoid(double a, double b, double tolerance, double fa, double fb) {
double h = b - a; // compute the new step
double I1 = h*(fa + fb) / 2; // compute the first integral
double m = (a + b) / 2; // find the middle point
double fm = f(m); // function value at the middle point
h = h / 2; // make step two times smaller
double I2 = h*(0.5*fa + fm + 0.5*fb); // compute the second integral
Integral I;
if (abs(I2 - I1) <= tolerance) { // if tolerance is satisfied
I.value = I2;
}
else { // if tolerance is not satisfied
if (tolerance > 1e-15) // check that we are not requiring too high precision
tolerance /= 2; // request higher precision in every half
I += trapezoid(a, m, tolerance, fa, fm); // integrate the first half [a m]
I += trapezoid(m, b, tolerance, fm, fb); // integrate the second half [m b]
}
return I;
}
__device__ Integral& Integral::operator+=(Integral &I) {
this->value += I.value;
return *this;
}
Для простоты я использую здесь только один поток. Теперь, если я запускаю этот код, я получаю сообщение "Сбой cudaMemcpy: обнаружен незаконный доступ к памяти". Когда я запускаю «cuda-memcheck», я получаю эту ошибку:
========= Invalid __local__ write of size 4
========= at 0x00000b18 in C:/Users/User/Desktop/Integrator Stack/Integrator_GPU/kernel.cu:73:controller(Integral*, double, double, double)
========= by thread (0,0,0) in block (0,0,0)
========= Address 0x00fff8ac is out of bounds
Он говорит, что проблема связана со строкой 73, которая просто
double m = (a + b) / 2;
Может быть, в этот момент у меня заканчивается память?
Если я уменьшу интервал интегрирования, изменив правую границу с b = 800
на b = 700
в main
, программа будет работать нормально и даст правильный результат. Почему я получаю ошибку недопустимого доступа к памяти при простом создании новой переменной?
Кроме того, у меня есть идентичная версия этой программы для процессора, и она работает без нареканий, так что алгоритм расчета, скорее всего, правильный.
+=
неверна. И если вы думаете, что я скачиваю что-то анонимное по ссылке из дропбокса, вы очень ошибаетесь..... - person talonmies   schedule 21.05.2018const
к параметру функции перегрузки, чтобы заставить его скомпилироваться. Как уже указывалось, у вас, вероятно, не хватает места в стеке. Ограничение по умолчанию в моем тесте составляло 1024 байта на поток. Я понятия не имею, сколько рекурсии вы делаете или сколько места в стеке на самом деле требуется, но когда я добавилcudaSetLimit(cudaLimitStackSize, 32768ULL);
в начало вашего кода, он заработал для меня без каких-либо ошибок во время выполнения. Я не говорю, что результаты правильные; похоже, вы используете неинициализированную общую память, и в вашем коде могут быть другие ошибки. - person Robert Crovella   schedule 22.05.2018cudaDeviceSetLimit
действительно заставляет этот конкретный код работать правильно (и вычисляет правильное значение), но я получаю ту же ошибку, когда увеличиваю интервал интегрирования. Кроме того, я получаю предупреждение о размере стека в выводе компиляции. Итак, проблема теперь выявлена, думаю, лучше переписать алгоритм без рекурсии. Спасибо! - person Yashman   schedule 22.05.2018