Ora un programma ottimizzato per il calcolo del prodotto scalare tra vettori.


// PRODOTTO SCALARE TRA DUE VETTORI
// implementare lettura vettori da file

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

#define NBLOCKS 1024
#define THREADS 1024

const long int N = 100000;
//const int blocksPerGrid = imin( 32, (N+threadsPerBlock-1)/threadsPerBlock );

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

__global__ void scalar_prod( long int *a,  long int *b,  long int *c){

// la componente trattata viene individuata sommando un offset dovuto
  // al fatto di trovarsi in un certo blocco e quindi pari all'indice
    // del blocco per il numero di threads contenuti in un blocco
    // a cui si aggiunge l'indice del trhread corrente in quel blocco
    int tid = threadIdx.x + blockDim.x*blockIdx.x;

    // tid continua a venire incrementato, fino a quando non si eccede il numero massimo N di componenti

    while(tid < N){
      c[tid]=a[tid]*b[tid]; // salto un numero di threads pari a tutti quelli permessi lungo x    
      tid += blockDim.x*gridDim.x;
    }
    // questo while è un loop che mi permette di applicare il kernel ripetutamente (itero le 
    // 128*128 = 16.000 circa istanze = gridDim.x)
} 

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

int main( void ){

  long int a[N], b[N], c[N];
  long int *dev_a, *dev_b, *dev_c;
  long int  scalarprod = 0;

  // riempimento dei vettori secondo questa scelta
  for(int i=0; i
    a[i] = 2*i;
    b[i] = -i;
  };

  // alloco memoria sulla scheda
  cudaMalloc( (void **)&dev_a, N*sizeof(long int) );
  cudaMalloc( (void **)&dev_b, N*sizeof(long int) );
  cudaMalloc( (void **)&dev_c, N*sizeof(long int) );

  // faccio partire il cronometro (non conto il tempo di creazione dei vettori)
  cudaEvent_t start,stop;
  cudaEventCreate(&start);
  cudaEventCreate(&stop);
  cudaEventRecord(start,0);

  cudaMemcpy( dev_a, a, N*sizeof(long int), cudaMemcpyHostToDevice );
  cudaMemcpy( dev_b, b, N*sizeof(long int), cudaMemcpyHostToDevice );

// scalar_prod<<<(N+127)/128,128>>>(dev_a, dev_b, dev_c);
   scalar_prod<<>>(dev_a, dev_b, dev_c);

  cudaMemcpy( c, dev_c, N*sizeof(long int), cudaMemcpyDeviceToHost ); // passo all'host
  // può ora essere grande come la memoria ram della scheda grafica!!

   for( int i=0; i
      cout << "scalar product = " << scalarprod << "\n";

  // fermo il cronometro (chiedo che la gpu abbia finito prima di restituire il val di stop
  cudaEventRecord(stop,0);
  cudaEventSynchronize(stop);
  float elapsed;
  cudaEventElapsedTime(&elapsed, start, stop);
  cout << "time: " << elapsed/1000. << endl;

  // dealloco
  cudaEventDestroy(start);
  cudaEventDestroy(stop);
  cudaFree(dev_a);
  cudaFree(dev_b);
  cudaFree(dev_c);

  return 0;
}