Метод Симпсона для интеграции реальных функций с CUDA

Основная проблема здесь в том, что вы не понимаете свой собственный алгоритм.

Ваша функция integralSimpsonMethod() разработана таким образом, что каждый поток производит выборку не менее 3 точек квадратуры на подинтервал в интегральной области. Следовательно, если вы выберете n так, чтобы оно было менее чем в четыре раза больше количества потоков в вызове ядра, неизбежно, что каждый подинтервал будет перекрываться, и результирующий интеграл будет неверным. Вы должны убедиться, что код проверяет и масштабирует количество потоков или n, чтобы они не создавали перекрытия при вычислении интеграла.

Если вы делаете это для чего-то другого, кроме самоназидания, я рекомендую вам посмотреть составную версию правила Симпсона. Это гораздо лучше подходит для параллельной реализации и будет значительно более производительным при правильной реализации.

Я бы предложил подход к интеграции Симпсона с использованием CUDA Thrust. В основном вам нужно пять шагов:

1. Сгенерируйте квадратурные веса Симпсона;
2. Создайте точки выборки функции;
3. Сгенерируйте значения функции;
4. Вычислить поэлементное произведение между квадратурными весами и значениями функции;
5. Суммируйте вышеуказанные произведения.

Шаг №1 требует создания массива с многократно повторяющимися элементами, а именно 1 4 2 4 2 4 ... 1 для случая Симпсона. Этого можно добиться, позаимствовав подход Роберта Кровеллы в библиотека толчков cuda повторяет вектор несколько раз< /а>.

Шаг № 2 можно выполнить, используя подход couting_iterators и заимствования talonmies в Назначение и использование counting_iterators в библиотеке CUDA Thrust.

Шаг № 3 — это приложение thrust::transform.

Шаги №4 и №5 можно выполнить вместе с помощью thrust::inner_product.

Этот подход также можно использовать, когда интерес представляют другие правила квадратурной интеграции.

Вот код

#include <thrust/iterator/counting_iterator.h>
#include <thrust/iterator/transform_iterator.h>
#include <thrust/iterator/permutation_iterator.h>
#include <thrust/iterator/counting_iterator.h>
#include <thrust/iterator/constant_iterator.h>
#include <thrust/inner_product.h>
#include <thrust/functional.h>

#include <thrust/fill.h>
#include <thrust/device_vector.h>
#include <thrust/host_vector.h>

// for printing
#include <thrust/copy.h>
#include <ostream>

#define STRIDE 2
#define N 100

#define pi_f  3.14159265358979f                 // Greek pi in single precision

struct sin_functor
{
    __host__ __device__
    float operator()(float x) const
    {
        return sin(2.f*pi_f*x);
    }
};

template <typename Iterator>
class strided_range
{
    public:

    typedef typename thrust::iterator_difference<Iterator>::type difference_type;

    struct stride_functor : public thrust::unary_function<difference_type,difference_type>
    {
        difference_type stride;

        stride_functor(difference_type stride)
            : stride(stride) {}

        __host__ __device__
        difference_type operator()(const difference_type& i) const
        {
            return stride * i;
        }
    };

    typedef typename thrust::counting_iterator<difference_type>                   CountingIterator;
    typedef typename thrust::transform_iterator<stride_functor, CountingIterator> TransformIterator;
    typedef typename thrust::permutation_iterator<Iterator,TransformIterator>     PermutationIterator;

    // type of the strided_range iterator
    typedef PermutationIterator iterator;

    // construct strided_range for the range [first,last)
    strided_range(Iterator first, Iterator last, difference_type stride)
    : first(first), last(last), stride(stride) {}

    iterator begin(void) const
    {
        return PermutationIterator(first, TransformIterator(CountingIterator(0), stride_functor(stride)));
    }

    iterator end(void) const
    {
        return begin() + ((last - first) + (stride - 1)) / stride;
    }

    protected:
        Iterator first;
        Iterator last;
        difference_type stride;
};

int main(void)
{
    // --- Generate the integration coefficients
    thrust::host_vector<float> h_coefficients(STRIDE);
    h_coefficients[0] = 4.f;
    h_coefficients[1] = 2.f;

    thrust::device_vector<float> d_coefficients(N);

    typedef thrust::device_vector<float>::iterator Iterator;
    strided_range<Iterator> pos1(d_coefficients.begin()+1, d_coefficients.end()-2, STRIDE);
    strided_range<Iterator> pos2(d_coefficients.begin()+2, d_coefficients.end()-1, STRIDE);

    thrust::fill(pos1.begin(), pos1.end(), h_coefficients[0]);
    thrust::fill(pos2.begin(), pos2.end(), h_coefficients[1]);

    d_coefficients[0]       = 1.f;
    d_coefficients[N-1]     = 1.f;

    // print the generated d_coefficients
    std::cout << "d_coefficients: ";
    thrust::copy(d_coefficients.begin(), d_coefficients.end(), std::ostream_iterator<float>(std::cout, " "));  std::cout << std::endl;

    // --- Generate sampling points
    float a     = 0.f;
    float b     = .5f;

    float Dx    = (b-a)/(float)(N-1);

    thrust::device_vector<float> d_x(N);

    thrust::transform(thrust::make_counting_iterator(a/Dx),
        thrust::make_counting_iterator((b+1.f)/Dx),
        thrust::make_constant_iterator(Dx),
        d_x.begin(),
        thrust::multiplies<float>());

    // --- Calculate function values
    thrust::device_vector<float> d_y(N);

    thrust::transform(d_x.begin(), d_x.end(), d_y.begin(), sin_functor());

    // --- Calculate integral
    float integral = (Dx/3.f) * thrust::inner_product(d_y.begin(), d_y.begin() + N, d_coefficients.begin(), 0.0f);

    printf("The integral is = %f\n", integral);

    getchar();

    return 0;
}