Matlab & CUDA

Scuola Estiva di Calcolo Avanzato 2008

Piero Lanucara, Fabio BonaccorsoCASPUR

Via dei Tizii,6b – 00185 Roma – Italyhttp://www.caspur.it

E-mail:[email protected]: [email protected]

La flessibilita’ di MatLab e la potenza delle GPU

http://www.caspur.it/

mailto:[email protected]

mailto:[email protected]


Matlab e Cuda: 2

1. La GPU

2. La piattaforma CUDA

3. Esecuzione di programmi CUDA

4. Integrazione con Matlab attraverso i Mex file

5. Jacket

6. Esempi

Sommario


Matlab e Cuda: 3

• La GPU e' un processore dedicato a calcoli basilariper la grafica– Incorpora hardware specifico per velocizzare

operazioni matematiche tipiche in grafica

• La GPU “scarica” la CPU per questi calcoli, raggiungendo un'efficienza non uguagliabile in software– Le vecchie GPU acceleravano operazioni 2d, oggi si

occupano di tutto il processo anche nel 3d

– Hanno capacita' di calcolo intorno ai 1000 GFlop/s

La GPU: Graphic Processing Unit


Matlab e Cuda: 4

• Performance crescenti nel tempo

GPU vs CPU: prestazioni


Matlab e Cuda: 5

• Performance del sistema di memoria

GPU vs CPU: prestazioni

La “banda” dimemoria e' la misuradella velocita' con cui la CPU (GPU) legge(o scrive) i dati in memoria

Una banda maggiorepermette di trasferirepiu' dati per unità ditempo


Matlab e Cuda: 6

GPU vs CPU: architettura interna

Nella GPU la maggior parte dell’hardware e’ utilizzato per le molte unita' di calcolo

• Architettura “dedicata” (specialpurpose)

Nella CPU molto hardware e’ dedicato al controllo del flusso di esecuzione e alla cache

• Architettura flessibile (generalpurpose)


Matlab e Cuda: 7

• Una CPU (oggi) si sfrutta al massimomassimizzando l'utilizzo della cache– La cache “avvicina” i dati alla CPU

– La memoria ha maggiore latenza e minore banda

– Essendo general-purpose e' adatta a tutti glialgoritmi

• La GPU nasce per il calcolo data-parallel– Stesso codice su operandi diversi (vettori, matrici)

– Idonea per algoritmi con alta intensita' numerica

GPU vs CPU


Matlab e Cuda: 8

• CUDA e’ un’architettura hardware e software per sfruttare la potenza delle moderne GPU per ilcalcolo (scientifico)• Acronimo di Compute Unified Device Architecture

• Supera i problemi dei precedenti approcci GP-GPU– Non usa un’API grafica (OpenGL, Direct3D)

• E' uno standard proprietario di NVIDIA• E' multipiattaforma (C con estensioni + runtime)

– Linux (32 e 64 bit)– Windows (32 e 64 bit)

CUDA: Cos'e'


Matlab e Cuda: 9

CUDA: modello di esecuzione

• Il programma inizia sulla CPU (Host)

• La CPU esegue una parte di codice– Trasferimenti CPU<->GPU

• Il codice (sulla CPU) invoca un “kernel” eseguito dalla GPU (Device)– In parallelo da piu’ “threads”

• La CPU esegue altro codice– Trasferimenti CPU<->GPU

• ...

• Il programma termina sulla CPU


Matlab e Cuda: 10

Una GPU e' composta da:

• Uno o piu' multiprocessori

– Il numero dipende dal modello della GPU

• Memoria (device memory)

– Globale

– Costante

– Texture

CUDA: implementazione Hardware


Matlab e Cuda: 11

• Ogni multiprocessore ha:

– Una sola unita’ per l’esecuzione del codice (istruzioni)

– Memoria condivisa

– Uno o piu’ processori

• Ogni processore consiste di:

– Una sola unita’ di calcolo

– Registri

– Memoria locale (opzionale)

CUDA: implementazione Hardware


Matlab e Cuda: 12

• Thread– Flusso di esecuzione indipendente

• I thread eseguono lo stesso codice su operandi diversi

– Ogni thread ha un ID univoco (threadIdx)

• Permette di selezionare i dati su cui opera il thread

• Il blocco e’ un gruppo di thread– Il blocco e' assegnato ad un unico multiprocessore

– Se il blocco contiene piu’ thread del numero di processori, i thread verranno eseguiti in time-slice

– Puo' essere un insieme lineare, 2D o 3D di thread

– Ha una sua memoria condivisa

• Fino al termine dell'esecuzione del blocco di threads da parte del multiprocessore

CUDA: Definizioni


Matlab e Cuda: 13

• La griglia CUDA e’ un gruppo di blocchi– Puo' essere un insieme lineare, 2D o 3D di blocchi

– Permette di mappare il problema in maniera naturale sui differenti blocchi

– Ogni blocco ha un ID univoco nella griglia (blockIdx)

• CUDA e' progettato per lavorare efficientemente con un alto numero di threads

– La priorita’ e’ far calcolare tutti i processori

– Bisogna trovare la dimensione ottimale del blocco

CUDA: Definizioni


Matlab e Cuda: 14

• Condivisa (Shared) Accesso veloce– Visibile da tutti i threads di uno stesso blocco

• Globale Accesso costoso (senza cache)

– Memoria visibile (lettura e/o scrittura) da tutti i threads

• Costante Accesso medio (cache)

• Texture Accesso medio (cache)– Sfrutta unita' hw per operazioni di filtraggio

CUDA: Memoria


Matlab e Cuda: 15

• Per compilare il codice CUDA si usa il compilatore nvcc

• Supporta le estensioni CUDA al C

• Compila file sorgenti con sintassi C++

– Piu' restrittivo del “C” sul controllo dei tipi delle variabili

– La sintassi per la programmazione ad oggetti non e' (ad oggi) supportata

CUDA: il compilatore nvcc


Matlab e Cuda: 16

• Il codice CUDA permette di definire un “kernel”

• Il “kernel” viene eseguito dalla griglia CUDA– Opera sui dati nelle memorie della GPU (globale, locale)

• Bisogna portare i dati dalla CPU alla GPU e viceversa– Funzioni dell'API di CUDA

• Alcune funzioni possono eseguire in modo asincrono

– La CPU puo' lavorare mentre la GPU esegue la chiamata

CUDA: modello di programmazione


Matlab e Cuda: 17

• Per definire un kernel CUDA si usa la parola chiave__global__ prefisso alla funzione.

• Esempio:__global__ void vecAdd(float* A, float* B){

int i = threadIdx.x;

B[i] = A[i] + B[i];

}

• In questo modo, nvcc compila il codice per la GPU

• Le variabili su cui opera il kernel sono sulla GPU– Distinte da quelle sulla CPU

– Vanno trasferite prima di chiamare il kernel

CUDA: estensioni al “C”


Matlab e Cuda: 18

• Se in un kernel CUDA si vuole chiamare un’altra funzione, questa va dichiarata con la parola chiave __device__

• Esempio:__device__ float scalarAdd(float a, float b){

return a + b;

}

__global__ void vecAdd(float* A, float* B){

int i = threadIdx.x;

B[i] = scalarAdd(A[i], B[i]);

}



Matlab e Cuda: 19

• Per invocare il kernel nel codice si usa l'operatore<<< dimGrid, dimBlock >>>

• Esempio:vecAdd <<<1 , N >>> (a_dev, b_dev);

La griglia e' da 1 blocco da N threads• Per dimensionare le griglie e i blocchi CUDA c'e' un

tipo nuovo, dim3 che specifica un vettore 1D, 2D o 3D

• dim3 dimGrid1(10); // Griglia di 10 blocchi• dim3 dimGrid2(20,10); // Griglia di 20x10 blocchi• dim3 dimBlock1(3,3,3); // Blocco di 3x3x3 thread



Matlab e Cuda: 20

• Funzioni per allocare/deallocare memoria sulla GPU:cudaError_t cudaMalloc(void** devPtr, size_t count)

cudaError_t cudaFree(void* devPtr)

• Funzioni per trasferire dati CPU <-> GPU:cudaError_t cudaMemcpy( void* dst, const void* src,

size_t count, enum cudaMemcpyKind kind)

kind puo’ essere cudaMemcpyHostToDevice oppurecudaMemcpyDeviceToHost

CUDA: principali funzioni dell’API


Matlab e Cuda: 21

__global__ void vecAdd(float* A, float* B)

{

int x = threadIdx.x;

int y = threadIdx.y * blockDim.x;

int blSz = blockDim.x * blockDim.y;

int bl = (blockIdx.y * gridDim.x + blockIdx.x) * blSz;

int i = bl + ( x + y );

if (i < N)

B[i] = alpha*A[i] + B[i];

}

int main()

{

dim3 dimBlock(16,8);

int grid_x = 16*4;

int bloc_x = dimBlock.x;

int bloc_y = dimBlock.y;

int grid_y = (N + grid_x*bloc_x*bloc_y - 1) / (grid_x * bloc_x*bloc_y);

dim3 dimGrid (grid_x, grid_y);

float A[N], B[N];

VecInit(A,B, N);

const unsigned int mem_size = sizeof(float) * N;

void *a_dev, *b_dev;

cudaMalloc( (void**) &a_dev, mem_size); // Alloco memoria per A sulla GPU

cudaMalloc( (void**) &b_dev, mem_size); // Alloco memoria per A sulla GPU

cudaMemcpy( a_dev, A, mem_size, cudaMemcpyHostToDevice); // Copio A sulla GPU

cudaMemcpy( b_dev, B, mem_size, cudaMemcpyHostToDevice); // Copio B sulla GPU

vecAdd <<<dimGrid,dimBlock>>> (a_dev, b_dev);

cudaMemcpy( B, c_dev, mem_size, cudaMemcpyDeviceToHost); // Copio B dalla GPU

PrintVec(B, N);

}

CUDA: esempio di kernel


Matlab e Cuda: 22

• L’accesso alla memoria globale della GPU e’ costoso• Il blocco puo’ utilizzare la memoria shared che e’ piu’ veloce• La parola chiave da utilizzare e’ __shared__• Esempio:

__shared__ float As[BLOCK_SIZE];

• Tutti i processi di uno stesso blocco (fino alla completaesecuzione del kernel) possono leggere/scrivere la memoriacondivisa

• L’accesso da parte dei threads di un blocco deve esserecoordinato dal programmatore– Si utilizza la funzione __syncthreads() che rende visibile a tutti i

thread una modifica alla memoria condivisa

CUDA: memoria shared


Matlab e Cuda: 23

Esempio:__global__ void vecIncr(float* A)

{

__shared__ float A_blocco[32];

int id = threadIdx.x; // Ogni thread opera su un indice differente

A_blocco[id] = A[id]; // Ogni thread riempie una parte della memoria shared

__syncthreads(); // Ora i cambiamenti sono visibili a tutto il blocco

A_blocco[32-id] ++; // Ogni thread aggiorna una posizione

__syncthreads(); // Tutti vedranno le modifiche

A[id] = A_blocco[id]; // Ogni thread copia dalla memoria shared alla

// memoria globale, altrimenti i cambiamenti verranno

// perduti

}

CUDA: memoria shared


Matlab e Cuda: 24

• Per verificare il corretto comportamento del programma CUDA, si puo’ compilarlo utilizzandol’opzione –deviceemu di nvcc

• Il codice girera’ sulla CPU e sara’ possibile usare un debugger per eseguire le istruzioni un passo allavolta, esaminare variabili, ecc.

• Il flusso dell’algoritmo sara’ diverso– La CPU non eseguira’ tutti i thread in parallelo

• In questa modalita’ nel kernel si possono usarefunzioni di stampa, apertura file, ... non disponibilinell’esecuzione sulla GPU

CUDA: modalita’ emulazione


Matlab e Cuda: 25

• La piattaforma CUDA viene fornita insieme a due librerie, che permettono di utilizzareimplementazioni efficienti di algoritmi di base per ilcalcolo scientifico– CUBLAS

– CUFFT

• Permettono un’astrazione maggiore rispetto allaprogrammazione con i kernel CUDA

• Nella SDK ci sono numerosi esempi pronti per l’usoo da prendere a modello

Librerie CUDA


Matlab e Cuda: 26

• NVIDIA fornisce un'implementazione della libreriaBLAS (Basic Linear Susbsystem Subroutine), checomprende funzioni tipiche dell'algebra lineare:

– Operazioni tra vettori (BLAS1)

– Operazioni tra matrici e vettori (BLAS2)

– Operazioni tra matrici (BLAS3)

• La libreria e' in C

– Vengono forniti i wrapper per il Fortran che sioccupano anche dei traferimenti CPU<->GPU

CUDA: libreria BLAS


Matlab e Cuda: 27

• NVIDIA fornisce un'implementazione di una libreria per calcolare le fft (fast fourier transforms), modellata sulla libreria open-source FFTW (v2):

– Trasformate 1D, 2D e 3D

• La libreria e' in C

– Vengono forniti i wrapper per il Fortran

CUDA: libreria FFT


Matlab e Cuda: 28

• Dopo aver compilato il programma CUDA, si lancia normalmente

• L'utilizzo della GPU puo' durare al piu' 5 sec

– Altrimenti “congela” il monitor

– Si mette una scheda video secondaria

• Il runtime CUDA deve essere nel path corrente

CUDA stand-alone


Matlab e Cuda: 29

• Calcolare prodotto di due matrici con CUBLAS

– Demo presente nell'SDK CUDA

• Basta chiamare cublas_sgemm invece di sgemm!

– Esempio con performance su AMD Opteron 2.8Ghz

CUDA: esempio Blas3


Matlab e Cuda: 30

do n=128,2560,32

! Initialize the matrices A,B and C

A = 1.

B = 2.

C = 3.

! With the prescribed inputs, each element of the C matrix

! should be equal to c_right

c_right= 2*n + 3

call cpu_time(time_start)

call cublas_SGEMM ('n','n',n,n,n,alpha,A,n,B,n,beta,C,n)

call cpu_time(time_end)

gflops = 1.e-9*(2*n*n*n - n*n)/(time_end-time_start)

print "(i5,1x,a,1x,f8.4,2x,a,f12.4)", n, &

" time =",time_end-time_start, GFLOPS=", gflops

end do

CUDA: esempio Blas3


Matlab e Cuda: 31

CUDA: prestazioni Blas3

• CUDA ver 1.1 • CUDA ver 2.0 (beta)


Matlab e Cuda: 32

• Una reduction e’ un’operazione che opera su un insieme di dati e restituisce un valore

– Minimo, massimo, media, ecc. di un vettore

• Nell'SDK c’e’ il codice sorgente di implementazioni CUDA

• Esempio avanzato da studiare e utilizzare!

CUDA: reduction


Matlab e Cuda: 33

• Utilizzando Matlab, l'implementazione CUDA si occupera’ soltanto del calcolo vero e proprio

– Le matrici, vettori, ecc. vengono inizializzati nell'ambiente Matlab

– Il calcolo viene lanciato sulla GPU attraverso il MEX

• Nel MEX dobbiamo trasferire i dati verso la GPU e poi indietro negli array predisposti da Matlab

• Il tutto sempre ricordando che si lavora in singola precisione floating point (per ora)

CUDA con Mex


Matlab e Cuda: 34

CUDA+Mex: moltiplicazione di matrici

#include "mex.h"

#include "mclcppclass.h"

void mexFunction(int nlhs,

mxArray * plhs[], int nrhs,

const mxArray * prhs[])

{

float *A, *B, *C;

float one = 1.0, zero = 0.0;

int m, n, p;

char *chn = "N";

/* Retrieve the input data */

A = (float*)mxGetPr(prhs[0]);

B = (float*)mxGetPr(prhs[1]);

/* Find the dimension of the

data */

m = mxGetM(prhs[0]);

p = mxGetN(prhs[0]);

n = mxGetN(prhs[1]);

/* Error checking */

if (p != mxGetM(prhs[1]))

mexErrMsgTxt(“Dim mismatch");

/* Create an mxArray for the

output data */

mwSize dims[2] = {m,n};

mxClassID tipo = mxSINGLE_CLASS;

plhs[0] = mxCreateNumericArray(2,

dims, tipo, mxREAL);

/* Create a pointer to the

output data */

C = mxGetPr(plhs[0]);

/* Pass all arguments to Fortran

by reference */

cublas_sgemm(chn, chn, &m, &n, &p,

&one, A, &m, B, &p, &zero, C, &m);

}


Matlab e Cuda: 35

CUDA con Mex: prestazioni

• Matlab su Opteron • La versione sulla GPU raggiunge i 100 Gflop/s

• In questo caso incide l'overhead del mex file, rispetto alla versione CUDA stand-alone


Matlab e Cuda: 36

• L'utilizzo della GPU da un programma Matlab richiede la scrittura di un Mex file e dell'implementazione della routine in CUDA

• Jacket permette invece un uso trasparente della GPU

– Si lavora a livello di vettori e matrici

– Crea i mex file necessari

– Lancia il calcolo sulla GPU

Matlab -> Jacket -> CUDA:alto livello e buone prestazioni


Matlab e Cuda: 37

• Jacket e’ un plugin per Matlab prodotto dalla AccelerEyes(ancora in beta–version).

• Jacket, intervenendo ad un livello di astrazione notevolmente superiore ad altri strumenti analoghi, rimpiazza le strutture dati (originarie) di MATLAB, normalmente sulla CPU, in strutture dati sulla GPU.

• Tutte le operazioni saranno effettuate sulla GPU invece che sulla CPU (in modalità data–parallel).

• Alla fine del calcolo sarà possibile riavere i dati sulla CPU utilizzando le funzionalità già presenti in MATLAB.

Jacket e Matlab - 1


Matlab e Cuda: 38

Vediamo velocemente l’uso di Jacket con un primo esempio.

Per utilizzare Jacket dare (dalla linea comandi di MATLAB):>> addpath(’/scratch/jacket/engine’);

Questo comando ci consente di utilizzare i files MEX presenti nella directory e corrispondenti a tutte le funzioni incluse in Jacket.

Verifichiamo il nostro ambiente di lavoro:>> ginfoJacket v0.4 (build ST67) data memory: 0mb used, 1491mb freeGPU0 Quadro FX 5600, 1350 MHz, 1535mb VRAM, Capability 1.0

Jacket e Matlab - 2


Matlab e Cuda: 39

>> X=gsingle(magic(3));>> Y=gsingle(ones(3));>> A=X*YA =15 15 1515 15 1515 15 15>> whosName Size Bytes Class AttributesA 3x3 496 gsingleX 3x3 496 gsingleY 3x3 496 gsingle

Jacket e Matlab - 3


Matlab e Cuda: 40

• L’output del comando whos mostra come la matrice A abbia lo stessotipo delle matrici di input. Il calcolo (prodotto matriciale) è statoeffettuato, quindi, sulla GPU.

• La funzione gsingle consente il cast di una matrice MATLAB in una matrice in singola precisione sulla GPU:

>> A=gsingle(B); %push B to the GPU from the CPU

• Per tornare indietro sulla CPU dobbiamo utilizzare la funzione built–inMATLAB double (o single):

>> B=single(A);>> whosName Size Bytes Class AttributesB 3x3 36 single

Jacket e Matlab - 4


Matlab e Cuda: 41

Le tecniche per vettorizzare un codice MATLAB mostrate in precedenza sonodi aiuto per avere un codice efficiente sulla GPU. Lo stesso codice avràbeneficio immediato dall’utilizzo di Jacket.Codice vettorizzato MATLAB:>> tic;A=(ones(5000)+eye(5000));B=sqrt(A);max(max(B)),tocans = 1.4142Elapsed time is 1.487942 seconds.

Le funzioni gones e geye creano due strutture dati sulla GPU analoghe allefunzioni built–in MATLAB.Codice Jacket:>> tic;A=(gones(5000)+geye(5000));B=sqrt(A);max(max(B)),tocans = 1.4142Elapsed time is 0.029784 seconds.

Jacket e Matlab - 5


Matlab e Cuda: 42

• Non c’e speranza di migliorare le prestazioni di un codice cheesegue una serie di calcoli su variabili scalari, come quello sotto:

a = 2 . 5 ; b = 3 . 3 3 ; c = 6 . 2 3 ;for i = 1:1000000

a = b + c ;b = a − c ;c = b * a ;i f ( a > b )

a = b − c ;e l s e

b = b + 1 2 . 3 ;end

End

• I codici vanno, pertanto, vettorizzati.

Jacket e Matlab - 6


Matlab e Cuda: 43

Il prossimo esempio e’ la versione Jacket dello script benchdgemm utilizzato per testare le prestazioni del prodotomatriciale built--in in MATLAB.Il codice utilizzato e’:addpath('/scratch/jacket/engine‘)% Warm-up GPUa = gsingle(rand(200));b = a;c=a*b; gforce(c);% Begin benchicont=1;for n=128:32:2560

a = gsingle(rand(n)); b = a;tic; c=a*b; gforce(c); ttt=toc;jacket(icont) = (2*n^3-n^2)/(ttt*10^9);icont = icont+1;

end

v=128:32:2560;plot(v,jacket);xlabel('Matrix dimension')ylabel('GFLOP/s‘)

Jacket e Matlab - 7


Matlab e Cuda: 44

Osservazioni sullo script:•Nel modello di programmazione di CUDA e’ possibile che il calcolo sulla GPUsia dilazionato di modo che questo venga effettuato solo quando il risultato e’necessario.

•La funzione gforce viene utilizzata per forzare il calcolo sulla GPU. In generalepuo’ essere utilizzata per ottimizzare un codice Jacket.•E’ buona norma “riscaldare" la GPU chiamando una prima volta la parte dicodice di cui vogliamo misurare le prestazioni. Questo consentira’ diprecaricare le librerie per poter prendere il massimo delle prestazioni dalcalcolo.

Mandiamo in esecuzione il codice:>> addpath('/scratch/jacket/engine');>> benchdgemm_jacket

Jacket e Matlab - 8


Matlab e Cuda: 45

Come si vede le prestazioni sono eccellenti, con cambiamentiminimi rispetto alla versione Matlab

Jacket e Matlab - 9


Matlab e Cuda: 46

Riprendiamo l'esempio di Laboratorio relativo all'algoritmo di bilateral image denoising filter e partendo dalla versione vettorizzata mostriamo il codice perGPU:B = bilateral_jacket(A,sd,sr,R);% The bilateral image denoising filter, Jacket versionA=gsingle(A); % Jacket #1B = gzeros(size(A)); % Jacket #2zsum=gzeros(size(A)); % Jacket #3ilast=size(A,1);jlast=size(A,2);i = 1:ilast;j = 1:jlast;for m = -R:R

ipm=i+m;ip=ipm(ipm>0 & ipm<=ilast);im=i(ipm>0 & ipm<=ilast);for n = -R:R

jpn=j+n;jp=jpn(jpn>0 & jpn<=jlast);jm=j(jpn>0 & jpn<=jlast);z = exp(-(A(ip,jp)-A(im,jm)).^2/(2*sd^2)) * exp(-(m^2 + n^2)/(2*sr^2));B(im,jm) = B(im,jm) + z.*A(ip,jp);zsum(im,jm)=zsum(im,jm)+z;

endendB=B./zsum;gforce(B); % Jacket #2

Laboratorio: bilateral image denoisingfilter con Jacket – 1


Matlab e Cuda: 47

Mandiamo in esecuzione il nuovo codice:>> A=single(rgb2gray(imread('susi.jpg')));>> sd=10;sr=10;R=3;>> addpath('/scratch/jacket/engine');

>> tic;Bj=bilateral_jacket(A,sd,sr,R);tocElapsed time is 0.102315 seconds.

>> Bj=double(Bj);>> norm(Bv-Bj,2)ans = 0.0023

Il confronto rispetto alla versione vettorizzata per CPU mostra comunque un buonguadagno anche per questo codice (di un fattore 6) e una leggera degradazione dellaqualita’ dell'output che non sembra pregiudicare la bonta’ del risultato.

Laboratorio: bilateral image denoisingfilter con Jacket – 2

Matlab & CUDA

Documents

Transcript of Matlab & CUDA