Le operazioni tra le matrici si può dire che siano il punto forte di CUDA: la scheda grafica infatti è appositamente costruita per gestire tutti i pixel dello schermo contemporaneamente, quindi una matrice ad esempio 1280x800. Si capisce che un calcolo effettuato in serie con una potente CPU può essere svolto molto più rapidamente in parallelo tramite una GPU, e anche se ogni processore della scheda grafica presenta una potenza di calcolo nettamente inferiore rispetto al processore principale del computer, il calcolo del prodotto degli elementi di due matrici punto per punto è solo una banale operazione di moltiplicazione, svolta quindi alla medesima velocità.
Cercheremo quindi in questo programma di far gestire alla scheda grafica tutta la matrice contemporaneamente (o almeno a blocchi), di modo da effettuare tutti i calcoli elemento per elemento in un solo singolo colpo di clock.


#include < iostream > 
#include < cstdlib >
#include < cmath >
#define EPS 1e-6
#include < ctime >
using namespace std;
#define MAX 30
#define THREADS 32
typedef float dato; // così­ modifico solo questo float/double se voglio cambiare scheda

// funzione per controllare la coerenza dei risultati cpu-gpu
void check(dato *a, dato *b,int N) {

for(int i = 0; i < N; ++i)
for(int j = 0 ; j < N ; ++j )
if(fabs(a[i*N+j]-b[i*N+j])>EPS) {
cerr << "Errore in posizione (" << i << "," << j <<"), a: "\
<< a[i*N+j] << ", check: " << b[i*N+j] << endl; 
return;
}
cerr << "OK" << endl;
return;
};

//Versione cpu del prodotto matrice matrice
void matrix_product_cpu(dato *a,dato *b, dato *c,int N){

dato temp;
for(int i = 0; i < N; ++i) 
for(int j = 0; j < N; ++j) {
temp = 0; 
for(int k = 0; k < N;++k) 
temp += a[i*N+k]*b[k*N+j];
c[i*N+j] = temp;
}
return;
};

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

// Prodotto matriciale: versione ottimizzata. 
// Si considerano matrici quadrate di ordine N  multiplo della dimensione del blocco di threads THREADS (quadrato).
// E' possibile ottimizzare il calcolo di un gruppo di righe per un gruppo di colonne
// specificamente di un blocco di elementi di A di dimensione THREADSxN e un blocco di B di dimensione Nx THREADS. 
// Per ridurre il numero di accessi alla memoria è conveniente sfruttare la memoria condivisa della scheda grafica.
// La memoria condivisa è presente in quantità limitata e quindi diventa necessario ridurre il problema in tanti sottoproblemi.
// Piu' precisamente, il sottoproblema elementare è costituito da un prodotto di due matrici di dimensione THREADSxTHREADS
// memorizzate sulla memoria veloce (shared) della GPU

__global__ void matrix_product( dato* A, dato* B, dato* C, int width) {

    // la mia posizione è data da due coordinate
    int ib = blockIdx.y;
    int jb = blockIdx.x;        
    int it = threadIdx.y;
    int jt = threadIdx.x;
        
    // Indice della prima sottomatrice di A elaborata dal blocco
    // width è un multiplo intero di THREADS (larghezza matrice completata)
    // aBegin  include un certo numero ib  di gruppi di blocchi rettangolari THREADSxwidth
    int aBegin = width*THREADS*ib; // posiz in alto a sx dell'ib-esimo blocco 
    
    // Indice dell'ultima sottomatrice di A elaborata dal blocco
    int aEnd   = aBegin + width - 1; // condizione che mi dica qual è l'ultimo elemento della riga per farci i conti
    
    // numero di colonne tra una sottomatrice e la successiva
    int aStep  = THREADS;
    
    // Indice della prima sottomatrice di B elaborata dal blocco
    // bBegin include un certo numero jb di blocchi di colonne, blocchi larghi THREADS 
    int bBegin = THREADS*jb;
         
    // numero di elementi tra una sottomatrice e la successiva    
    int bStep  = THREADS * width;
     
    // Csub è usata come variabile in cui memorizzare il valore dell'elemento di C
    // calcolato dal thread
    // Viene aggiornato ripetutamente nel ciclo for seguente
    dato Csub = 0;    
     
    // Iterazione sulle sottomatrici in cui viene suddiviso il calcolo degli elementi del blocco 
// a <= aEnd controllo solo sulle a, tanto è quadrata
// questo ciclo for mi fa saltare da una sottomatrice di a alla successiva, idem per b
    for (int a = aBegin, b = bBegin; a <= aEnd; a += aStep, b += bStep) {

        // Dichiarazione della variabile in cui salvare la sottomatrice di A in esame
__shared__ dato As[THREADS][THREADS];
        __shared__ dato Bs[THREADS][THREADS];

        // Carico gli elementi di ciascuna sottomatrice in memoria condivisa: ogni thread del 
// blocco carica un elemento!      
As[it][jt] = A[a +  width*it + jt]; 
Bs[it][jt] = B[b +  width*it + jt];

        // Sincronizzo per assicurarmi che ogni thread del blocco abbia caricato gli elementi.
__syncthreads();

  // Calcolo i contributi agli elementi di matrice di C dovute alle sottomatrici in esame        
for (int k = 0; k < THREADS ; ++k ) Csub += As[it][k]*Bs[k][jt];
// l'elemento C[it][jt] viene aggiornato in nblocks passaggi, pari al numero di iterazioni
// del ciclo for
  
        // Sincronizzo per assicurarmi che il calcolo precedente sia terminato prima di caricare nuove
// sottomatrici
__syncthreads();
    }

    // Scrivo i risultati in C. Ogni thread elabora un elemento di C. 
    int c = width*THREADS*ib + THREADS*jb; // conto quante righe ci sono prima + di quante mi sposto         
    C[c +  width*it + jt] = Csub;    
};

// costruisco una cornice di zeri per ignorare ciò che eccede la dimensione N*N della matrice originaria
__global__ void completa_matrice(dato *a, int N) {

int i = blockIdx.y*blockDim.y + threadIdx.y;
int j = blockIdx.x*blockDim.x + threadIdx.x;
int offset = gridDim.x*blockDim.x*i + j;
        
        // se sono fuori dalla matrice N*N pongo tutto a zero!
if(i >= N || j >= N) a[offset] = 0.;

};

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

void matrix(dato *a, dato *b, dato *c, int N) {

dato *a_dev,*b_dev,*c_dev;
int nblock = (N+THREADS-1)/THREADS; // se N = nblock è già un multiplo intero di threads questa operazione non fa nulla
int width = nblock*THREADS; // calcolo larghezza matrice completata

cout << "width=" << width << "  dim matrix=" << N << "\n";
        
        // alloco la memoria su cuda per le due matrici di input e per quella di output
cudaMalloc((void**)&a_dev, width*width*sizeof(dato));
cudaMalloc((void**)&b_dev, width*width*sizeof(dato));
cudaMalloc((void**)&c_dev, width*width*sizeof(dato));

        // faccio partire il cronometro
cudaEvent_t start,stop;
cudaEventCreate(&start);
cudaEventCreate(&stop);
cudaEventRecord(start,0);

// copio da host a device: copio sulla scheda le matrici create nel main
        // copia bidimensionale: voglio accederci con doppia indicizzazione
cudaMemcpy2D(a_dev, width*sizeof(dato), a, N*sizeof(dato), N*sizeof(dato), N, cudaMemcpyHostToDevice);
cudaMemcpy2D(b_dev, width*sizeof(dato), b, N*sizeof(dato), N*sizeof(dato), N, cudaMemcpyHostToDevice);
// a ha N righe di larghezza N*sizeof(dato)
// a_dev è la destinazione sulla scheda di W*W elementi (width*sizeof(dato))
             // la matrice di partenza ha N*N elementi, ma la matrice di arrivo ha W*W elementi: alloco W*W elementi

        // lancio il kernel cuda: blocchi lungo x e lungo y, e all´interno threads lungo x e lungo y
dim3 c_threads(THREADS,THREADS);
dim3 c_blocks(nblock,nblock);

dim3 threads(THREADS,THREADS);
dim3 blocks(nblock,nblock);

// ho messo in a_dev e in b_dev una cornicetta di zeri
if(N % THREADS != 0) {
completa_matrice<<>>(a_dev,N);
       completa_matrice<<>>(b_dev,N);
}
// blocks contiene nblock nblock, threads contiene THREADS THREADS  
matrix_product<<>>(a_dev,b_dev,c_dev,width);

// copio il risultato in c (sulla cpu)
cudaMemcpy2D(c, N*sizeof(dato), c_dev, width*sizeof(dato), N*sizeof(dato), N, cudaMemcpyDeviceToHost);
cudaEventRecord(stop,0);
cudaEventSynchronize(stop);
float elapsed;
cudaEventElapsedTime(&elapsed,start,stop);
cout << N << " " << elapsed/1000. << endl;

// dealloco
cudaEventDestroy(start);
cudaEventDestroy(stop);
cudaFree(a_dev);
cudaFree(b_dev);
cudaFree(c_dev);
};

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

int main(int argc, char **argv) {

cudaSetDevice(0);
const int N = atoi(argv[1]); // N è la dimensione della matrice
// metto le matrici in vettori di N*N elementi
dato *a_host = new dato[N*N];
dato *b_host = new dato[N*N];
dato *c_host = new dato[N*N];
dato *c_check = new dato[N*N];

srand48(time(NULL));

// riempio matrici a e b con i numeri random
for(int i = 0; i < N ;++i)
for(int j = 0; j < N;++j) {
int index = i*N+j;
a_host[index] = drand48()*MAX; // numeri tra 0 e 30 (define MAX)
b_host[index] = drand48()*MAX;
}

float elapsed_cpu;
matrix(a_host,b_host,c_host,N); // esegue la maggior parte delle operazioni, definita sopra

clock_t start,stop;
start = clock();
matrix_product_cpu(a_host,b_host,c_check,N);
stop = clock();
elapsed_cpu = (float)(stop - start)/CLOCKS_PER_SEC;
cout << elapsed_cpu << endl;

check(c_host,c_check,N);
delete[](a_host);
delete[](b_host);
delete[](c_host);
delete[](c_check);

return 0;
};