vector_axpy.cu

Aller à la documentation de ce fichier.

// includes, system
#include <iostream>
#include <stdlib.h>

// Utilities and system includes
#include <cuda_runtime.h>
#include <curand.h>

#include <helper_functions.h>
#include <helper_cuda.h>

#include <curand.h>

__device__ double axpy(double a,double x,double y){
  return  a*x+y;
}

__global__ void vector_axpy_kernel(const double a,const double *d_X, const double *d_Y, double *d_Z, int N){
  const int      tid = blockDim.x * blockIdx.x + threadIdx.x;
  const int      Incl= blockDim.x * gridDim.x; 

  for (int Pos=tid;Pos<N ;Pos+=Incl)
    d_Z[Pos]= a*d_X[Pos]+d_Y[Pos];
}

void vector_axpy(const double a,const double *X, const double *Y, double *Z, int N){
  for (int i=0;i<N;i++)
    Z[i]=a*X[i]+Y[i];
}

double errorinf(const double *X,const double *Y,int N){
  double norm=fabs(X[0]-Y[0]);
  for (int i=1;i<N;i++){
    double s=fabs(X[i]-Y[i]);
    if (s>norm) norm=s;
  }
  return norm;
}

void printArray(const double *X,int N,int n){
  int nd=N-n,np=n;
  if (n>N){
    np=N;nd=N+1;
  }
  std::cout.precision(16);  
  for (int i=0;i<np;i++)
    std::cout << " [" << i << "]: " << X[i] << std::endl;
  if (nd<N)
    std::cout << " ...\n";
  for (int i=nd;i<N;i++)
    std::cout << " [" << i << "]: " << X[i] << std::endl;
  
}

int main(int argc, char** argv){
  double *d_X,*d_Y,*d_Z;         // device variable  (on GPU memory)
  double *h_X,*h_Y,*h_Z,*h_Zgpu; // host   variable  (on CPU memory)
  cudaError_t err;
  int N=1<<24;
  std::cout << "Start ...\n";
  
  // 1) Arrays allocations on <host> 
  h_X = new double[N]; 
  h_Y = new double[N]; 
  h_Z = new double[N];
  h_Zgpu = new double[N];
  if ((h_X == NULL)||(h_Y == NULL)||(h_Z == NULL)||(h_Zgpu == NULL)){
    fprintf(stderr,"Allocation error on CPU\n");
    exit(EXIT_FAILURE);
  }
  
  // 2) Arrays allocations on <device> 
  checkCudaErrors( cudaMalloc((void**) &d_X, N * sizeof(double)) );
  cudaMalloc((void**) &d_Y, N * sizeof(double));
  err = cudaGetLastError();
  if( cudaSuccess != err) { // d_Y allocation failed
    fprintf(stderr, "%s(%i) : CUDA Malloc error : (%d) %s.\n",
                __FILE__,__LINE__, (int)err, cudaGetErrorString( err ) );
    checkCudaErrors(cudaFree(d_X));  // Free <device> array d_X
    exit(EXIT_FAILURE);
  } 
  cudaMalloc((void**) &d_Z, N * sizeof(double));
  err = cudaGetLastError();
  if( cudaSuccess != err) { // d_Z allocation failed
    fprintf(stderr, "%s(%i) : CUDA Malloc error : (%d) %s.\n",
                __FILE__,__LINE__, (int)err, cudaGetErrorString( err ) );
    checkCudaErrors(cudaFree(d_X)); // Free <device> array d_X
    checkCudaErrors(cudaFree(d_Y)); // Free <device> array d_X
    exit(EXIT_FAILURE);
  }
  
  // 3) Set pseudo random generator on <host> using CURAND library
  curandGenerator_t prngCPU;
  checkCudaErrors(curandCreateGeneratorHost(&prngCPU, CURAND_RNG_PSEUDO_MTGP32));
  checkCudaErrors(curandSetPseudoRandomGeneratorSeed(prngCPU, 777));
  
  // 4) Generate uniformly distributed random numbers in double precision. Values are between 0.0 and 1.0, 
  // where 0.0 is excluded and 1.0 is included.
  checkCudaErrors(curandGenerateUniformDouble(prngCPU, h_X, N));
  checkCudaErrors(curandGenerateUniformDouble(prngCPU, h_Y, N));
  
  // 5) Commpute Z <- a*X+Y on <host>
  vector_axpy(2.,h_X,h_Y,h_Z,N);
  std::cout << "h_Z=\n";
  printArray(h_Z,N,3);
  
  // 6) Copy <host> arrays h_X and h_Y respectively in <device> arrays d_X and d_Y 
  checkCudaErrors( cudaMemcpy(d_X, h_X, N * sizeof(double), cudaMemcpyHostToDevice) );
  checkCudaErrors( cudaMemcpy(d_Y, h_Y, N * sizeof(double), cudaMemcpyHostToDevice) );
  
  // 7) Commpute Z <- a*X+Y on <device> using 256 threads and 512 blocks
  /* dim3 dimBlock(256, 1, 1);
     dim3 dimGrid(512, 1, 1);
     vector_axpy_kernel<<< dimGrid, dimBlock>>>(2.,d_X,d_Y,d_Z,N);*/
  vector_axpy_kernel<<< 512, 256>>>(2.,d_X,d_Y,d_Z,N);
  
  // 8) Copy <device> array d_Z on <host> array h_Zgpu
  checkCudaErrors( cudaMemcpy(h_Zgpu, d_Z, N * sizeof(double), cudaMemcpyDeviceToHost) );
  
  std::cout << "\nh_Zgpu=\n";
  printArray(h_Zgpu,N,3);
  
  // 9) Print max|h_Z - h_Zgpu|
  std::cout.precision(16);
  std::cout << "Error : " << errorinf(h_Zgpu,h_Z,N);
  
  // 10) Free <device> arrays
  checkCudaErrors(cudaFree(d_X));
  checkCudaErrors(cudaFree(d_Y));
  checkCudaErrors(cudaFree(d_Z));
  
  // 11) Free <host> arrays
  delete [] h_X; delete [] h_Y; delete [] h_Z; delete [] h_Zgpu;
  std::cout << "\n...Stop\n";
}