Un example d'application en cuda de l'algorithme de scan pour le calcul des intégrales

Télécharger le projet
Commenter

Description

Utilisation de l'algorithme de scan pour le calcul de l'aire sous la courbe d'une fonction f ===> intégrale.
L'algorithme de scan est le plus basique j'ajouterais plus tard des versions plus performants.
C'est un code fait rapidement qui fonctionne correctement.
Dans cet example je calcule l'intégrale de x² entre 0..1 avec un pas de 0.001, mais vous pouvez modifier ces paramètres et prendre une autre fonction.

Source / Exemple : 


#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;
}
Télécharger le projet

Codes Sources

A voir également

S'inscrire maintenant

Vous n'êtes pas encore membre ?

inscrivez-vous, c'est gratuit et ça prend moins d'une minute !

Les membres obtiennent plus de réponses que les utilisateurs anonymes.

Le fait d'être membre vous permet d'avoir un suivi détaillé de vos demandes et codes sources.

Le fait d'être membre vous permet d'avoir des options supplémentaires.