Per confermare la potenza del calcolo su scheda grafica mediante CUDA proponiamo un programma di calcolo del prodotto scalare di due vettori casuali di lunghezza data che lo esegue sia su GPU che su CPU. I risultati sono sorprendenti, soprattutto aumentando il numero di punti. Ovviamente una scheda grafica potente avrà prestazioni maggiori di un processore medio, ma un processore molto potente può eseguire calcoli più velocemente di una scheda grafica con poca memoria allocata. In ogni caso, per vettori molto lunghi la supremazia della scheda grafica risulta schiacciante.


// PRODOTTO SCALARE TRA DUE VETTORI OTTIMIZZATO
// implementare lettura vettori da file
// confronto tra CPU e GPU

#include < iostream >
using namespace std;
#include < cstdlib >

#define imin(a,b) (a
const int N = 10000000;

const int threadsPerBlock = 1024; // grandezze per l'ottimizzazione del problema (alternativa al 128*128 di prima)
const int blocksPerGrid = imin(32, (N+threadsPerBlock-1)/threadsPerBlock); // trucco per ottimizzare il numero di blocchi
//const int blocksPerGrid = 1024;

/////////////////////////////////////////////////////////////////////////////////////////////////////

__global__ void scalar_prod( float *a,  float *b,  float *c ){

  __shared__ float cache[threadsPerBlock];   // shared: interna e comune ai threads di un singolo blocco
  int tid = threadIdx.x + blockDim.x*blockIdx.x;
  int cacheIndex = threadIdx.x; // indice del thread
  float temp = 0;

  while(tid < N){
    temp += a[tid]*b[tid]; // incremento dell'indice della componente sommata (salto di blocco in blocco!)
        tid += blockDim.x*gridDim.x;
  };

  // il numero di elementi da valutare può essere maggiore del numero di thread * numero di blocchi
  // la somma (temp) corre sui diversi gruppi di blocchi combinando le componenti associate 
  // allo stesso valore di tid (modulo il numero totale di threads)

  // scopo del prodotto scalare è di restituire un numero; parte di questo numero viene calcolato sulla GPU
  // la somma parziale appena calcolata viene messa nella memoria shared del blocco in corrispondenza del valore
  // di threadIdx.x
  cache[cacheIndex] = temp; // 512 (threads) somme parziali per blocco, ogni blocco ha la sua cache[]
 // sono privati all'interno del blocco: non si mischiano tra di loro!
  __syncthreads(); // aspetto che tutti i thread del blocco abbiano finito

        // valuto la somma parziale dei risultati dei vari threads di un blocco
  int i=blockDim.x/2;
  while(i!=0){
    if (cacheIndex < i) cache[cacheIndex] += cache[cacheIndex+i];
        // questa somma combina i valori della metà superiore del blocco di threads con quelli della 
        // meta' inferiore e li associa a indici della meta' inferiore
    __syncthreads();  // attendo che tutti i threads abbiano effettuato la somma
    // tutti gli elementi rilevanti hanno indici da zero a meta' del blocco nella seconda iterazione si
    // sommano gli elementi del secondo quarto con quelli del primo quarto del blocco, e così via
    i/=2; // poichè i è una variabile intera, nel caso 1/2 = 0.5 viene restituito 0
  };

  if (cacheIndex ==0) c[blockIdx.x] = cache[0];
} // il risultato finale e' la somma parziale di un blocco

/////////////////////////////////////////////////////////////////////////////////////////////////////
// prodotto scalare su cpu (confronto dei tempi)
void sp_cpu_check(float *a, float *b, float &check){
int i;
  check=0;
  for (i=0; i
    check += a[i]*b[i];
  }
}

/////////////////////////////////////////////////////////////////////////////////////////////////////

int main( void ){

  float *a, *b, *partial_c;
  float *dev_a, *dev_b, *dev_partial_c;
  float  sp_gpu = 0, sp_cpu = 0;

  a = (float *)malloc(N*sizeof(float));
  b = (float *)malloc(N*sizeof(float));
  partial_c = (float *)malloc(blocksPerGrid*sizeof(float));

  // riempimento casuale di vettori (escluso dal conteggio dei tempi)
  for( int i=0; i
    a[i] = drand48();
    b[i] = drand48();
  };

  // alloco la memoria sulla scheda
  cudaMalloc( (void **)&dev_a, N*sizeof(float) );
  cudaMalloc( (void **)&dev_b, N*sizeof(float) );
  cudaMalloc( (void **)&dev_partial_c, blocksPerGrid*sizeof(float) );

  // faccio partire il cronometro della gpu
  cudaEvent_t start1,stop1;
  float time1;
  cudaEventCreate(&start1);
  cudaEventCreate(&stop1);
  cudaEventRecord(start1,0);
  
  cudaMemcpy( dev_a, a, N*sizeof(float), cudaMemcpyHostToDevice );
  cudaMemcpy( dev_b, b, N*sizeof(float), cudaMemcpyHostToDevice );

  scalar_prod<<>>(dev_a, dev_b, dev_partial_c);

  cudaMemcpy( partial_c, dev_partial_c, blocksPerGrid*sizeof(float), cudaMemcpyDeviceToHost );

  // infine le somme parziali dei vari blocchi vengono combinate 
  for( int i=0; i
      sp_gpu += partial_c[i];
  };
cout << "scalar product gpu = " << sp_gpu << "\n";

  // fermo il cronometro della gpu
  cudaEventRecord(stop1,0);
  cudaEventSynchronize(stop1);
  cudaEventElapsedTime(&time1, start1, stop1);
  cout << "time gpu: " << time1 << "\n";
  
  // delloco dalla gpu
  cudaEventDestroy(start1);
  cudaEventDestroy(stop1);
  cudaFree(dev_a);
  cudaFree(dev_b);
  cudaFree(dev_partial_c);

  ///////////////////////////////////////////
  // rifaccio il calcolo con la cpu
  cudaEvent_t start2,stop2;
  float time2;
  cudaEventCreate(&start2);
  cudaEventCreate(&stop2);
  cudaEventRecord(start2,0);

  //check del prodotto scalare, calcolato sulla cpu
  sp_cpu_check(a,b,sp_cpu);
cout << "scalar product cpu = " << sp_cpu << "\n";

  // fermo il cronometro della cpu
  cudaEventRecord(stop2,0);
  cudaEventSynchronize(stop2);
  cudaEventElapsedTime(&time2, start2, stop2);
cout << "time cpu: " << time2 << "\n";

  // dealloco dalla cpu
  cudaEventDestroy(start2);
  cudaEventDestroy(stop2);
  free(a);
  free(b);
  free(partial_c);

  return 0;
}