UN EXAMPLE D'APPLICATION EN CUDA DE L'ALGORITHME DE SCAN POUR LE CALCUL DES INTÉGRALES

#include <cuda.h>
#include <cuda_runtime.h>
#include <stdlib.h>
#include <iostream>


__device__ double f( double const x )
{
	double res = x*x;
	return res;
}

__device__ double Trap( double const h, double const a, double const b)
{
  double Base = f(b);
  double base = f(a);
  double Aire = h*( Base + base )/2;
  return Aire;
}


template<int BLOCKSIZE>
__global__ void scan( double* vect, double* sum, int n)
{
  int  i = threadIdx.x + blockDim.x*blockIdx.x;
  int tx = threadIdx.x;
  __shared__ double shared_mem[BLOCKSIZE]; 

  if( tx < BLOCKSIZE && i < n )
    shared_mem[threadIdx.x] = vect[i];
  __syncthreads();
  for(int step = 1 ; step < BLOCKSIZE; step *= 2)
    {
      int index = tx + step;
      if( tx%2 == 0 && i < n && index  < BLOCKSIZE )
	{
	  shared_mem[tx] += shared_mem[index];	
	}
      __syncthreads();
    }
  if( tx == 0)
    sum[blockIdx.x] = shared_mem[tx];
  __syncthreads();
}


__global__ void InitVect(double* vect, double a, double b, double h, int n)
{
  int i = threadIdx.x + blockDim.x*blockIdx.x;
  double x = a + i*h;
  if( i < n )
    vect[i] = Trap(h, x, x + h);
}


double CalculIntegrale( double a, double b, int n)
{
  double h = (b - a )/n;
  double* vect;
  double* sum_h, *sum_d;
  double res = 0;
  int size = ceil( (float)n/256);
  sum_h = new double[size];
  cudaMalloc( (void**)&vect, n*sizeof(double));
  if( cudaGetLastError() != cudaSuccess )
    std::cout << "Malloc 1  "<<cudaGetErrorString(cudaGetLastError()) << std::endl;
  cudaMalloc( (void**)&sum_d, size*sizeof(double));
  if( cudaGetLastError() != cudaSuccess )
    std::cout << "Malloc 2 "<<cudaGetErrorString(cudaGetLastError()) << std::endl;

  dim3 grid(ceil( (float)n/256),1,1 );
  dim3 block( 256,1,1);
  InitVect<<<grid, block>>>(vect, a, b, h, n);
  cudaThreadSynchronize();
  if( cudaGetLastError() != cudaSuccess )
    std::cout << "Vect Init "<<cudaGetErrorString(cudaGetLastError()) << std::endl;
  scan<256><<<grid, block >>>(vect, sum_d,n);
  cudaThreadSynchronize();
  if( cudaGetLastError() != cudaSuccess )
    std::cout << "Scan "<<cudaGetErrorString(cudaGetLastError()) << std::endl;
  
  cudaMemcpy( sum_h, sum_d, size*sizeof(double), cudaMemcpyDeviceToHost);
  if( cudaGetLastError() != cudaSuccess )
    std::cout << "Memcpy 1  "<<cudaGetErrorString(cudaGetLastError()) << std::endl;
  

  for(int i = 0; i < size; ++i)//Boucle pour finir la somme 
    res += sum_h[i];

  cudaFree( vect );
  cudaFree( sum_d );
  delete[] sum_h;
  return res;
}

int main(int argc, char* argv[])
{
  
  double integral = CalculIntegrale(0.,1.,1000);
  std::cout << integral << std::endl;
  return 0;
}