GPU = Graphics Processing Unit. Σε αντιδιαστολή με την CPU (Central Processing Unit), δηλαδή τον κεντρικό επεξεργαστή ο οποίος συχνά αποκαλείται και host, η GPU είναι μια κάρτα γραφικών, η οποία συχνά αποκαλείται και device. Η GPU διαθέτει μεγάλο πλήθος επεξεργαστικών μονάδων -συνήθως αρκετές δεκάδες ή εκατοντάδες- που καλούνται πυρήνες (cores). Οι πυρήνες είναι διαθέσιμοι στον προγραμματιστή ώστε μέσω μιας κατάλληλης βιβλιοθήκης μπορούν να χρησιμοποιηθούν για παράλληλη εκτέλεση προγραμμάτων.

Η πλατφόρμα CUDA (Compute Unified Device Architecture) είναι ένα προγραμματιστικό μοντέλο της εταιρείας nVidia το οποίο υποστηρίζει την εκτέλεση παράλληλων προγραμμάτων χρησιμοποιώντας τους πυρήνες (cores) μιας ή περισσοτέρων καρτών GPU. Η CUDA επιτρέπει στον προγραμματιστή να δημιουργήσει εκατοντάδες ή και χιλιάδες νήματα τα οποία εκτελούνται στους πυρήνες μιας ή περισσοτέρων GPU. Με τον τρόπο αυτό μπορούν να υλοποιηθούν αλγόριθμοι μαζικής παράλληλης επεξεργασίας. Σειριακοί αλγόριθμοι ή αλγόριθμοι με μικρό παραλληλισμό δε συμφέρει να υλοποιούνται με CUDA.

Από την ιστοσελίδα CUDA downloads θα χρειαστεί να κατεβούν και να εγκατασταθούν τρία βασικά πακέτα:

• το CUDA Toolkit
• οι CUDA drivers
• το CUDA SDK (Software Development Kit)

Η τρέχουσα έκδοση CUDA 5.5 υποστηρίζει τα εξής λειτουργικά συστήματα:

Υπάρχει πληθώρα καρτών GPU με υποστήριξη CUDA. Δείτε πληροφορίες σχετικά με τις συμβατές κάρτες γραφικών στο https://developer.nvidia.com/cuda-gpus.

Στο πακέτο της CUDA περιλαμβάνεται ένα προγραμματιστικό περιβάλλον που επιτρέπει στον χρήστη να αναπτύξει εφαρμογές GPU χρησιμοποιώντας μια επέκταση της γλώσσας C. Ωστόσο υποστηρίζονται και άλλες γλώσσες όπως FORTRAN, DirectCompute, OpenACC ενώ υπάρχουν επίσης JavaΠχ. jcuda και PythonΠχ. pycuda WrappersWrapper: ένα αντικείμενο (object) που προσφέρει μια συγκεκριμένη λειτουργικότητα.

Τα τελευταία χρόνια υπήρξαν ραγδαίες εξελίξεις στο αρχιτεκτονικό μοντέλο των καρτών GPU. Αν και υπήρχαν ήδη από το 1999 ο προγραμματισμός τους με σκοπό την υλοποίηση αλγορίθμων γενικού σκοπού (δηλ. όχι μόνο graphics) ήταν εξαιρετικά δύσκολος. Το μοντέλο ορόσημο που γέννησε την ιδέα των καρτών GPU για προγραμματισμό γενικού σκοπού ήταν η κάρτα GeForce 8800 που βασίστηκε στην αρχιτεκτονική G80 της nVidia. Συγκεκριμένα η αρχιτεκτονική G80

• ήταν η πρώτη αρχιτεκτονική GPU που υποστήριζε τη γλώσσα C
• διέθετε για πρώτη φορά ένα ενιαίο επεξεργαστή για υπολογισμούς graphics (vertex, geometry, pixel) και πράξεις γενικού σκοπού.
• εισήγαγε για πρώτη φορά το προγραμματιστικό μοντέλο Single-Instruction Multiple-Thread (SIMT) όπου πολλοί κλώνοι του ίδιου προγράμματος εκτελούνται ως διαφορετικά νήματα με διαφορετικά id.
• εισήγαγε την κοινή μνήμη (shared memory) και τον συγχρονισμό μέσω φραγμάτων (barrier synchronization) για την επικοινωνία μεταξύ νημάτων

Αν και η αρχιτεκτονική G80 βελτιώθηκε με το λανσάρισμα της G200 το 2008, η δεύτερη γενιά αρχιτεκτονικής GPU θεωρείται η Αρχιτεκτονική Fermi. Το παρακάτω σχήμα δείχνει το αρχιτεκτονικό διάγραμμα ενός chip GPU που ακολουθεί την Αρχιτεκτονική Fermi.

Στο chip αυτό υπάρχουν

• 16 Streaming Processors (SM) με 32 πυρήνες (cores) ο καθένας οπότε συνολικά υπάρχουν 512 πυρήνες. Οι SM είναι τοποθετημένοι γύρω από μια κοινή μνήμη cache 2ου επιπέδου μεγέθους 768 KB.
• 6 τράπεζες μνήμης 64bit που υποστηρίζουν μέχρι 6GB μνήμης τύπου GDDR5 DRAM.
• Ένα host-interface που συνδέει την GPU με την CPU μέσω διαύλου PCI Express.
• Ο δρομολογητής GigaThread engine ο οποίος κατανέμει τα blocks των νημάτων στους SM.

Η ιεραρχία των threads της CUDA (grid → blocks of threads → threads) αντιστοιχεί στην ιεραρχία των επεξεργαστών της GPU (GPU → SMs → cores):

Η Αρχιτεκτονική Kepler είναι η τρίτη και πιο πρόσφατη γενιά GPU.
Σύγκριση αρχιτεκτονικών Fermi, Kepler:

O Streaming Multiprocessor εκτελεί ένα ή περισσότερα block νημάτων. Αποτελείται από πολλούς πυρήνες (cores) όπως φαίνεται στο παρακάτω σχήμα. Το σχήμα δείχνει το αρχιτεκτονικό διάγραμμα ενός SM της αρχιτεκτονικής Fermi.

Ένας SM της αρχιτεκτονικής Fermi περιέχει:

• 32 πυρήνες (cores). Κάθε πυρήνας είναι ένας πολύ απλός επεξεργαστής ο οποίος περιέχει μια πλήρως pipelined Integer Arithmetic Logic Unit (ALU) και μια Floating Point Unit (FPU). Το διάγραμμα ενός πυρήνα δίνεται στο παρακάτω σχήμα:
  
  
• 16 μονάδες Load/Store
• 4 Special Function Units (SFU) οι οποίες εκτελούν πολύπλοκες μαθηματικές συναρτήσεις όπως ημίτονο (sine), συνημίτονο (cosine), αντίστροφο (1/x) και τετραγωνική ρίζα.
• 32768 Καταχωρητές (Registers) μήκους 32 bit.
• Μια Κοινή μνήμη μεγέθους 64 KB
• Μια κοινή Data Cache επιπέδου 1 και μια Instruction Cache
• Ένα κύκλωμα δικτύου για διασύνδεση με τους υπόλοιπους SM
• Δύο warp schedulers και δύο instruction dispatch units.

Ένας SM δρομολογεί νήματα σε groups που καλούνται warps. Το πλήθος των νημάτων σε ένα warp σχετίζεται άμεσα με το πλήθος των πυρήνων που διαθέτει ο SM. Στο συγκεκριμένο παράδειγμα κάθε warp έχει 32 νήματα. Υπάρχουν 2 warp-schedulers και δύο μονάδες ανάθεσης εντολών (instruction dispatch units). Δηλαδή, ο Fermi SM εκτελεί ταυτόχρονα 2 εντολές προερχόμενες από 2 διαφορετικά warps. Μια εντολή από ένα warp μπορεί, ανάλογα με το είδος της, να ανατεθεί

• είτε σε ένα πυρήνα από μια 16-άδα πυρήνων
• είτε σε μια από τις 16 Load/Store units
• είτε σε ένα από τα 4 Special Function Units

Ένα kernel προκαλεί την εκτέλεση ενός Πλέγματος (Grid)

• Ένα Πλέγμα αποτελείται από Ν Blocks νημάτων CUDA (Block₀, ..., Block_N-1).
  • Ένα Block αποτελείται από Τ Νήματα CUDA (Νήμα₀, ..., Νήμα_T-1).

Έτσι συνολικά υπάρχουν ΝΤ Νήματα.

Τα νήματα μέσα σε ένα block μπορούν να συνεργαστούν

• είτε μοιραζόμενα την κοινή μνήμη του block (shared memory)

• είτε κάνοντας συγχρονισμό εκτέλεσης χρησιμοποιώντας την εντολή __synchthreads(). Ένα νήμα που εκτελεί την εντολή αυτή σταματάει και περιμένει όλα τα νήματα του ιδίου block να την εκτελέσουν επίσης για να συνεχίσει. (Προσοχή στα if-then-else: κίνδυνος αδιεξόδου!!)

Νήματα σε διαφορετικά blocks δεν μπορούν να συνεργαστούν ούτε να συγχρονιστούν με την __synchthreads(). Μπορούν ωστόσο να προσπελάσουν την καθολική μνήμη (global memory) που είναι κοινή για όλα τα νήματα ανεξαρτήτως block. Η CPU βλέπει επίσης την καθολική μνήμη την οποία μπορεί να γράψει ή να διαβάσει χρησιμοποιώντας την εντολή cudaMemcpy().

Όλα τα νήματα εκτελούν τον ίδιο κώδικα, με τη διαφορά ότι έχουν ξεχωριστή ταυτότητα (id). Η ταυτότητα ενός νήματος καθορίζεται από 6 αριθμούς:

• την ταυτότητα του block όπου βρίσκεται το νήμα: blockIdx.x, blockIdx.y, blockIdx.z
• την ταυτότητα του νήματος μέσα στο block: threadIdx.x, threadIdx.y, threadIdx.z

Ο λόγος που χρησιμοποιούνται τρισδιάστες (3-D) ταυτότητες για τα blocks και για τα νήματα είναι ότι αυτό βολεύει συχνά όταν γράφουμε αλγορίθμους για επεξεργασία εικόνων (2-D προβλήματα) ή όγκων (3-D προβλήματα). Η ταυτότητα ενός νήματος ή η ταυτότητα του block στο οποίο εμπεριέχεται αυτό, μπορεί να χρησιμοποιηθεί για να γίνει προσπέλαση σε διαφορετικές θέσεις μνήμης και για να παίρνοται διαφορετικές αποφάσεις if-then-else σε σχέση με τα υπόλοιπα νήματα. Έτσι διαφοροποιούνται οι πράξεις που εκτελεί κάθε νήμα. Προφανώς, σε διαφορετική περίπτωση, αν όλα τα νήματα εκτελούσαν τις ίδιες πράξεις στα ίδια δεδομένα δεν θα υπήρχε παράλληλη επεξεργασία...
Το υλικό (hardware) δρομολογεί τα blocks σε οποιονδήποτε streaming multiprocessor (αν υπάρχουν περισσότεροι από 1 διαθέσιμοι).

Η βιβλιοθήκη CUBLASCUDA Basic Linear Algebra Subprograms αποτελείται από παράλληλες συναρτήσεις CUDA που υλοποιούν βασικές ρουτίνες γραμμικής άλγεβρας, όπως ανάγνωση και εγγραφή διανυσμάτων και πινάκων, πολλαπλασιασμό διανυσμάτων με πίνακες, πολλαπλασιασμό πινάκων με πίνακες, κλπ. Πλήρης τεκμηρίωση της βιβλιοθήκης καθώς και παραδείγματα κώδικα δίνονται εδώ. Η CUBLAS προσφέρεται ήδη ενσωματωμένη στο NVIDIA CUDA SDK οπότε δεν απαιτείται κάποιο επί πλέον βήμα εγκατάστασης. Ο πλήρης κατάλογος των συναρτήσεων της βιβλιοθήκης είναι ο παρακάτω. Χρησιμοποιείται η συντόμευση <t> που αφορά τον τύπο της εκάστοτε συνάρτησης. Οι δυνατές συντομεύσεις είναι:

Έτσι, γα παράδειγμα, η συνάρτηση cublasIdamax() βρίσκει τον δείκτη του στοιχείου με την μεγαλύτερη απόλυτη τιμή σε ένα διάνυσμα τύπου double.

• CUBLAS Datatypes Reference
  • cublasHandle_t
  • cublasStatus_t
  • cublasOperation_t
  • cublasFillMode_t
  • cublasDiagType_t
  • cublasSideMode_t
  • cublasPointerMode_t
  • cublasAtomicsMode_t
• CUBLAS Helper Function Reference
  • cublasCreate()
  • cublasDestroy()
  • cublasGetVersion()
  • cublasSetStream()
  • cublasGetStream()
  • cublasGetPointerMode()
  • cublasSetPointerMode()
  • cublasSetVector()
  • cublasGetVector()
  • cublasSetMatrix()
  • cublasGetMatrix()
  • cublasSetVectorAsync()
  • cublasGetVectorAsync()
  • cublasSetMatrixAsync()
  • cublasGetMatrixAsync()
  • cublasSetAtomicsMode()
  • cublasGetAtomicsMode()
• CUBLAS Level-1 Function Reference
  • cublasI<t>amax()
  • cublasI<t>amin()
  • cublas<t>asum()
  • cublas<t>axpy()
  • cublas<t>copy()
  • cublas<t>dot()
  • cublas<t>nrm2()
  • cublas<t>rot()
  • cublas<t>rotg()
  • cublas<t>rotm()
  • cublas<t>rotmg()
  • cublas<t>scal()
  • cublas<t>swap()
• CUBLAS Level-2 Function Reference
  • cublas<t>gbmv()
  • cublas<t>gemv()
  • cublas<t>ger()
  • cublas<t>sbmv()
  • cublas<t>spmv()
  • cublas<t>spr()
  • cublas<t>spr2()
  • cublas<t>symv()
  • cublas<t>syr()
  • cublas<t>syr2()
  • cublas<t>tbmv()
  • cublas<t>tbsv()
  • cublas<t>tpmv()
  • cublas<t>tpsv()
  • cublas<t>trmv()
  • cublas<t>trsv()
  • cublas<t>hemv()
  • cublas<t>hbmv()
  • cublas<t>hpmv()
  • cublas<t>her()
  • cublas<t>her2()
  • cublas<t>hpr()
  • cublas<t>hpr2()
• CUBLAS Level-3 Function Reference
  • cublas<t>gemm()
  • cublas<t>gemmBatched()
  • cublas<t>symm()
  • cublas<t>syrk()
  • cublas<t>syr2k()
  • cublas<t>syrkx()
  • cublas<t>trmm()
  • cublas<t>trsm()
  • cublas<t>trsmBatched()
  • cublas<t>hemm()
  • cublas<t>herk()
  • cublas<t>her2k()
  • cublas<t>herkx()
• BLAS-like Extension
  • cublas<t>geam()
  • cublas<t>dgmm()
  • cublas<t>getrfBatched()
  • cublas<t>getriBatched()

Η συγγραφή προγραμμάτων CUDA σε Windows γίνεται με χρήση της πλατφόρμας Visual Studio 2005, Visual Studio 2008 ή Visual Studio 2010 σε γλώσσα C/C++. H CUDA χρησιμοποιεί μια επέκταση της γλώσσας C προσθέτοντας επί πλέον σύμβολα και λέξεις κλειδιά (keywords) έτσι ώστε να γίνεται εκτέλεση προγραμμάτων στους επεξεργαστές GPU. Το NVIDIA GPU Computing Toolkit (CUDA Toolkit) προσφέρει τον μεταφραστή nvcc ο οποίος είναι απαραίτητος για την μετάφραση κώδικα CUDA.

Για να γράψουμε ένα πρόγραμμα C/C++ χρησιμοποιώντας τις δυνατότητες της βιβλιοθήκης CUDA ακολουθούμε τα εξής βήματα:

1. Από το Visual Studio 2010 δημιουργούμε νέο project: New Project → Visual C++ → Win32 Console Application. Επιλέγουμε Win32 Console Application ακόμη και σε έκδοση Windows 64bit
2. Δίνουμε όνομα στο project, πχ. cuFirst. Έτσι στον φάκελο Source Files δημιουργούνται αυτόματα τα αρχεία cuFirst.cpp και stdafx.cpp
3. Μετονομάζουμε το cuFirst.cpp σε cuFirst.cu. Επίσης αφαιρούμε το αρχείο stdafx.cpp από το project: [δεξί κλικ] Exclude From Project
4. Στο παράθυρο Solution Explorer επιλέγουμε το project cuFirst [δεξί κλικ] Build Customizations... τσεκάρουμε την επιλογή CUDA 4.2
5. Ανοίγουμε το μενού Properties του project: στο παράθυρο Solution Explorer επιλέγουμε το project cuFirst [δεξί κλικ] Properties. Στο μενού αυτό κάνουμε τις παρακάτω ρυθμίσεις:
   • VC++ Directories → Include Directories: προσθήκη του $(DXSDK_DIR)\include
   • CUDA C/C++ → Common → Additional Include Directories: προσθήκη των
     ./
     C:/ProgramData/NVIDIA Corporation/NVIDIA GPU Computing SDK 4.2/C/common/inc
     C:/ProgramData/NVIDIA Corporation/NVIDIA GPU Computing SDK 4.2/shared/inc
     
   • CUDA C/C++ → Device → Code Generation: προσθήκη των
     compute_10,sm_10
     compute_20,sm_20
     compute_30,sm_30
   • CUDA C/C++ → Host → Runtime Library: επιλογή του Multi-Threaded Debug (/MTd)
   • Linker → General → Additional Library Directories: προσθήκη των
     $(CudaToolkitLibDir)
     C:/NVIDIA/CUDA/CUDASamples/common/lib/$(PlatformName)
   • Linker → Input → Additional Dependencies: προσθήκη των cuda.lib και cublas.lib (εφόσον χρησιμοποιηθεί η CUBLAS).

Γράφουμε τον κώδικά μας στο αρχείο cuFirst.cu. Ο παρακάτω κώδικας προσθέτει δύο διανύσματα a, b μήκους N και βάζει το αποτέλεσμα στο διάνυσμα c.

	// cuDokimi.cu : Παράλληλη πρόσθεση διανυσμάτων σε GPU με χρήση CUDA
	//
	#include <stdlib.h>
	#include <stdio.h>
	#include <cuda.h>

	/*
	 * Συνάρτηση πυρήνα (kernel). Εκτελείται στους επεξεργαστές της GPU
	 */
	__global__ void add_vectors(float *a, float *b, float *c) {
		int i = blockIdx.x;
		
		c[i] = a[i] + b[i];
	}


	int main(void)
	{
		int count, i;
		
		// Βρες πόσες GPU διαθέτει το σύστημα (συνήθως 1)
		cudaGetDeviceCount(&count);
		printf("There are %d GPU devices in your system\n", count);

		int N = 10;	// Μήκος διανυσμάτων

		/*
		 * Δημιούργησε τα διανύσματα "a", "b", "c" στον host (CPU)
		 */
		float *a = (float *)malloc(N*sizeof(float));
		float *b = (float *)malloc(N*sizeof(float));
		float *c = (float *)malloc(N*sizeof(float));
		// Δώσε αρχικές τιμές στα "a", "b"
		for (i=0; i<N; i++) {
			a[i] = i - 0.5;
			b[i] = i*i - 3;
		}

		/*
		 * Δημιούργησε τα διανύσματα "a_dev", "b_dev", "c_dev" στο device (GPU)
		 */
		float *a_dev, *b_dev, *c_dev;
		cudaMalloc((void **)&a_dev, N*sizeof(float));
		cudaMalloc((void **)&b_dev, N*sizeof(float));
		cudaMalloc((void **)&c_dev, N*sizeof(float));
		
		/*
		 * Αντιγραφή των διανυσμάτων "a", "b" από τον host (CPU) στο device (GPU)
		 * δηλ. a --> a_dev,  b --> d_dev
		 */
		cudaMemcpy(a_dev, a, N*sizeof(float), cudaMemcpyHostToDevice);
		cudaMemcpy(b_dev, b, N*sizeof(float), cudaMemcpyHostToDevice);
		
		// Παράλληλη πρόσθεση c_dev[i] = a_dev[i] + b_dev[i]
		add_vectors<<< N, 1 >>>(a_dev, b_dev, c_dev);

		// Αντιγραφή του αποτελέσματος από το device (GPU) στον host (CPU)
		cudaMemcpy(c, c_dev, N*sizeof(float), cudaMemcpyDeviceToHost);
		
		/*
		 * Απελευθέρωση της μνήμης που δεσμεύεται από τα a_dev, b_dev, c_dev
		 */
		cudaFree(a_dev);
		cudaFree(b_dev);
		cudaFree(c_dev);

		// Εκτύπωση αποτελέσματος από τον host (CPU)
		printf("\n\n\nVector Addition Result:\n");
		for (i=0; i<N; i++) {
			printf("a[%d] : %0.2f \t+\t", i, a[i]);
			printf("b[%d] : %0.2f \t=\t", i, b[i]);
			printf("c[%d] : %0.2f\n", i, c[i]);
		}
		
		return 0;
	}
	

Αφού κάνουμε Build Project το εκτελέσιμο αρχείο cuFirst.exe βρίσκεται στον φάκελο Τα έγγραφά μου > Visual Studio 2010 > Projects > cuDokimi > Debug. Η εφαρμογή μας δεν είναι παραθυρική. Εκτελείται σε Command Window (Γραμμή Εντολών).

Το πρόγραμμα προσθέτει δύο διανύσματα a και b και γράφει το αποτέλεσμα στο διάνυσμα c. Όλα διανύσματα έχουν μήκος N. Τα αθροίσματα των στοιχείων εκτελούνται σε ξεχωριστά νήματα στους πυρήνες της GPU.
 

__global__ Δήλωση ότι η συνάρτηση θα εκτελεστεί από τους πυρήνες (cores) της GPU void add_vectors Συνάρτηση πυρήνα (kernel). ( float *a, float *b, float *c Οι μεταβλητές αναφέρονται στην μνήμη της συσκευής GPU και όχι στην κεντρική μνήμη της CPU ) {
     int i = blockIdx.x; Δείχνει την θέση x του block στο 2-Δ πλέγμα των blocks στο οποίο ανήκει το συγκεκριμένο νήμα
    
     c[i] = a[i] + b[i];
}

Tο πρόθεμα __global__ δηλώνει ότι η συνάρτηση add_vectors() θα εκτελεστεί στους πυρήνες (cores) της GPU. Γι' αυτό το λόγο η συνάρτηση λέγεται και συνάρτηση πυρήνα (kernel function ή απλά kernel). Προφανώς το πρόθεμα __global__ δεν ανήκει στη standard C, επομένως αν αυτό το πρόγραμμα μεταφραστεί από έναν κλασσικό C-compiler θα λάβουμε μήνυμα σφάλματος. Το πρόγραμμά μας γράφεται σε ένα αρχείο με κατάληξη .cu και θα μεταφραστεί από τον μεταφραστή nvcc που προσφέρει η NVIDIA. Ο μεταφραστής αυτός αποτελεί υπερσύνολο του κλασικού μεταφραστή C γνωρίζοντας επί πλέον τις έξτρα λέξεις κλειδιά και τις εντολές που απαιτούνται από την CUDA. Έτσι γνωρίζει την λέξη κλειδί __global__ και θα παράξει τον κατάλληλο κώδικα μηχανής ώστε η συνάρτηση να τρέξει στους πυρήνες της GPU. Σημειώνεται ότι ο nvcc θα μπορούσε επίσης να χρησιμοποιηθεί σαν ένας κλασικός μεταφραστής C ακόμη και σε προγράμματα που δεν χρησιμοποιούν CUDA.
 

int main(void) Εκτελείται στην CPU
{
     int count, i;
     // Βρες πόσες GPU διαθέτει το σύστημα (συνήθως 1)
     cudaGetDeviceCount(&count); Βρίσκει πόσες κάρτες GPU διαθέτει ο υπολογιστής.
     printf("There are %d GPU devices in your system\n", count);

Ο κώδικας της συνάρτησης main() εκτελείται στον host (CPU). Κατ'αρχήν ρωτάμε πόσες συσκευές GPU διαθέτει ο υπολογιστής μας χρησιμοποιώντας την συνάρτηση cudaGetDeviceCount()Αν και η συνάρτηση δίνει πληροφορίες σχετικα με την GPU εκτελείται στην CPU και όχι στη GPU..
 

     int N = 10; // Μήκος διανυσμάτων
    
     /*
     * Δημιούργησε τα διανύσματα "a", "b", "c" στον host (CPU)
     */
     float *a = (float *)malloc(N*sizeof(float));
     float *b = (float *)malloc(N*sizeof(float));
     float *c = (float *)malloc(N*sizeof(float));
     // Δώσε αρχικές τιμές στα "a", "b"
     for (i=0; i<N; i++) {
         a[i] = i - 0.5;
         b[i] = i*i - 3;
     }

Δημιουργούμε στη μνήμη της CPU τα διανύσματα a[], b[] και c[] μεγέθους N*sizeof(float). Γίνεται δέσμευση της κατάλληλης ποσότητας μνήμης με χρήση της standard συνάρτησης της C, malloc(). Οι πίνακες a[], b[] αρχικοποιούνται σύμφωνα με τον τύπο a[i] = i-0.5 και b[i] = i^2-3.
 

     /*
     * Δημιούργησε τα διανύσματα "a_dev", "b_dev", "c_dev" στο device (GPU)
     */
     float *a_dev, *b_dev, *c_dev;
     cudaMalloc((void **)&a_dev, N*sizeof(float));
     cudaMalloc((void **)&b_dev, N*sizeof(float));
     cudaMalloc((void **)&c_dev, N*sizeof(float));

Δημιουργούμε στη μνήμη της συσκευής GPU αντίστοιχους πίνακες a_dev[], b_dev[] και c_dev[] , χρησιμοποιώντας τη συνάρτηση cudaMalloc. Αυτό είναι απαραίτητο διότι η συνάρτηση πυρήνα add_vectors() μπορεί να χρησιμοποιήσει μόνο τη μνήμη της συσκευής.
 

     /*
     * Αντιγραφή των διανυσμάτων "a", "b" από τον host (CPU) στο device (GPU)
     * δηλ. a --> a_dev, b --> d_dev
     */
     cudaMemcpy(a_dev, a, N*sizeof(float), cudaMemcpyHostToDevice);
     cudaMemcpy(b_dev, b, N*sizeof(float), cudaMemcpyHostToDevice);

Αντιγράφουμε τους πίνακες a, b από τη μνήμη της CPU στους αντίστοιχους πίνακες a_dev, b_dev που βρίσκονται στη μνήμη της GPU χρησιμοποιώντας τη συνάρτηση cudaMemcpy(). Είμαστε έτοιμοι για την παράλληλη εκτέλεση της πρόσθεσης των διανυσμάτων a_dev, b_dev.
 

     // Παράλληλη πρόσθεση c_dev[i] = a_dev[i] + b_dev[i]
     add_vectors<<< N, 1 >>> Εκτέλεση της συνάρτησης πυρήνα add_vectors() στην GPU σε Ν blocks νημάτων με 1 νήμα ανά block (a_dev, b_dev, c_dev);

Γενικά, μια συνάρτηση πυρήνα (kernel) εκτελείται παράλληλα σε πολλαπλά αντίγραφα. Σε κάθε αντίγραφο αντιστοιχεί και ένα διαφορετικό νήμα (thread) οπότε δημιουργείται ένα πλήθος νημάτων ίσο με το πλήθος των αντιγράφων που καλούμε να εκτελέσουμε. Αυτό το μοντέλο παράλληλου προγραμματισμού είναι γνωστό ως Single Program Multiple Threads (SPMT). Τα νήματα (threads) οργανώνονται σε 3-Δ blocks νημάτων και τα blocks οργανώνονται σε ένα 3-Δ πλέγμα (grid) από blocks. Κάθε νήμα διακρίνεται από τα υπόλοιπα χάρη

• (α) στις συντεταγμένες του νήματος μέσα στο block, που δίνονται από τριάδα των ακεραίων threadIdx.x, threadIdx.y, threadIdx.z, και
• (β) στις συντεταγμένες του block μέσα στο πλέγμα blockIdx.x, blockIdx.y, blockIdx.z.

Για παράδειγμα, στον παρακάτω κώδικα εκτελείται η συνάρτηση πυρήνα f() σε ένα πλέγμα διαστάσεων [4 blocks] x [6 blocks] x [3 blocks], όπου κάθε block αποτελείται από ένα 3-Δ array νημάτων με διαστάσεις: [7 νήματα] x [2 νήματα] x [4 νήματα]:

Συνολικά λοιπόν δημιουργούνται 4 x 6 x 3 x 7 x 2 x 4 = 4032 νήματα. Κάθε νήμα διακρίνεται από τα υπόλοιπα χάρη στην εξάδα δεικτών [blockIdx.x, blockIdx.y, blockIdx.z, threadIdx.x, threadIdx.y, threadIdx.z]. Κάθε δείκτης μεταβάλλεται από την τιμή 0 έως την αντίστοιχη διάσταση - 1. Στο συγκεκριμένο παράδειγμα, το εύρος τιμών των διαφόρων δεικτών είναι:

Η κλήση της συνάρτησης πυρήνα: add_vectors<<< N, 1 >>>(a_dev, b_dev, c_dev); ισοδυναμεί με τον παρακάτω κώδικα:

Δημιουργείται συνεπώς ένα 1-Δ array νημάτων μήκους N=10. Μόνο ο δείκτης blockIdx.x μεταβάλλεται μεταξύ των τιμών 0,...,9. Όλοι οι άλλοι δείκτες είναι καρφωμένοι στο μηδέν.

Με τον τρόπο αυτό εκτελούνται παράλληλα οι 10 αθροίσεις. Η δρομολόγηση των νημάτων και η ανάθεσή τους στους κατάλληλους πυρήνες της GPU γίνεται από την CUDA χωρίς τον έλεγχο του χρήστη.
 

     // Αντιγραφή του αποτελέσματος από το device (GPU) στον host (CPU)
     cudaMemcpy(c, c_dev, N*sizeof(float), cudaMemcpyDeviceToHost);
    
     /*
     * Απελευθέρωση της μνήμης που δεσμεύεται από τα a_dev, b_dev, c_dev
     */
     cudaFree(a_dev);
     cudaFree(b_dev);
     cudaFree(c_dev);

Αντιγράφουμε το αποτέλεσμα από το array c_dev[] που βρίσκεται στη μνήμη της GPU στο array c[] που βρίσκεται στη μνήμη της CPU. Κατόπιν απελευθερώνουμε την μνήμη της GPU από τους πίνακες a_dev[], b_dev[], c_dev[] που δεν χρειάζονται πλέον.
 

     // Εκτύπωση αποτελέσματος από τον host (CPU)
     printf("\n\n\nVector Addition Result:\n");
     for (i=0; i<N; i++) {
         printf("a[%d] : %0.2f \t+\t", i, a[i]);
         printf("b[%d] : %0.2f \t=\t", i, b[i]);
         printf("c[%d] : %0.2f\n", i, c[i]);
     }
    
     return 0;
}

Η CPU τυπώνει τις πράξεις που γίνονται καθώς και το αποτέλεσμα έτσι ώστε να ελέγξουμε την ορθότητα του προγράμματος.
 

Ναι υπάρχουν περιορισμοί. Ο πιο απλός τρόπος αποτύπωσης των περιορισμών είναι η ανάγνωση των ιδιοτήτων της συσκευής μέσω της συνάρτησης cudaGetDeviceProperties. Η συνάρτηση αυτή επιστρέφει ένα struct τύπου cudaDeviceProp που περιλαμβάνει όλες τις πληροφορίες σχετικά με την κάρτα GPU. Για παράδειγμα, το παρακάτω πρόγραμμα κάνει χρήση αυτής της συνάρτησης και εμφανίζει στην οθόνη όλες τις πληροφορίες που επιστρέφονται στο struct prop.

	#include <stdlib.h>
	#include <stdio.h>
	#include <cuda.h>

	int main(void)
	{
		int count, i;
		cudaDeviceProp prop;

		cudaGetDeviceCount(&count);

		printf("There are %d GPU devices in your system\n", count);

		for (i=0; i<count; i++) {
			printf("----------\nDevice #%d description\n", i);
			// Get device properties
			cudaGetDeviceProperties(&prop, i);
			// Display device properties
			printf(" Device name: %s\n",
				prop.name);
			printf(" Amount of global memory on device: %d Bytes\n",
				prop.totalGlobalMem);
			printf(" Maximum amount of shared memory a single block may use: %d Bytes\n",
				prop.sharedMemPerBlock);
			printf(" Number of 32-bit registers available per block %d\n",
				prop.regsPerBlock);
			printf(" Number of threads in a warp: %d\n",
				prop.warpSize);
			printf(" Maximum pitch allowed for memory copies: %d Bytes\n",
				prop.memPitch);
			printf(" Maximum number of threads that a block may contain: %d\n",
				prop.maxThreadsPerBlock);
			printf(" Maximum number of threads allowed along each block dimension: %d x %d x %d\n",
				prop.maxThreadsDim[0], prop.maxThreadsDim[1], prop.maxThreadsDim[2]);
			printf(" Number of blocks allowed along each grid dimension: %d x %d x %d\n",
				prop.maxGridSize[0], prop.maxGridSize[1], prop.maxGridSize[2]);
			printf(" Amount of available constant memory: %d Bytes\n",
				prop.totalConstMem);
			printf(" Major revision of the device's compute capability: %d\n", 
				prop.major);
			printf(" Minor revision of the device's compute capability: %d\n",
				prop.minor);
			printf(" Device's requirement for texture alignment: %d\n",
				prop.textureAlignment);
			printf(" Device can simultaneously perform a cudaMemcpy() and kernel execution: %s\n",
				(prop.deviceOverlap ? "yes" : "no"));
			printf(" Number of multiprocessors on device: %d\n",
				prop.multiProcessorCount);
			printf(" Runtime limit for kernels executed on this device: %s\n",
				(prop.kernelExecTimeoutEnabled ? "yes" : "no"));
			printf(" Device is an integrated GPU (part of the chipset and not a discrete GPU): %s\n",
				(prop.integrated ? "yes" : "no"));
			printf(" Device can map host memory into the CUDA device address space: %s\n",
				(prop.canMapHostMemory ? "yes" : "no"));
			printf(" Device's computing mode: default, exclusive, or prohibited: %d\n",
				prop.computeMode);
			printf(" Maximum size supported for 1D textures: %d\n",
				prop.maxTexture1D);
			printf(" Maximum dimensions supported for 2D textures: %d x %d\n",
				prop.maxTexture2D[0], prop.maxTexture2D[1]);
			printf(" Maximum dimensions supported for 3D textures: %d x %d x %d\n",
				prop.maxTexture3D[0], prop.maxTexture3D[1], prop.maxTexture3D[2]);
			/*
			printf(" Maximum dimensions supported for 2D texture arrays: %d x %d x %d\n",
				prop.maxTexture2DArray[0], prop.maxTexture2DArray[1], prop.maxTexture2DArray[2]);
			*/
			printf(" Device supports executing multiple kernels within the same context simultaneously: %s\n",
				(prop.concurrentKernels ? "yes" : "no"));
		}

		return 0;
	}
	

Ανάλογα με την κάρτα GPU που διαθέτουμε, μια πιθανή έξοδος του παραπάνω προγράμματος είναι η παρακάτω. Όπως βλέπουμε υπάρχει όριο διαστάσεων 65535 x 65535 x 65535 για το πλέγμα των blocks και όριο 1024 x 1024 x 64 για το array των νημάτων μέσα σε κάθε block. Παρατηρούμε επίσης ότι υπάρχει όριο 1024 νημάτων για κάθε block. Με άλλα λόγια, ένα array νημάτων 1024 x 1 x 1 για κάθε block είναι αποδεκτό, αλλά τα arrays νημάτων με διαστάσεις 1024 x 2 x 1, 512 x 4 x 1, 512 x 2 x 2, 1024 x 1 x 2, δεν είναι αποδεκτά διότι το πλήθος των νημάτων υπερβαίνει το άνω όριο, 1024. Επίσης, για παράδειγμα, δεν είναι αποδεκτό το array νημάτων 1 x 2 x 128 διότι η τρίτη διάσταση υπερβαίνει το όριο 64.

	There are 1 GPU devices in your system
	----------
	Device #0 description
	 Device name: GeForce GT 520M
	 Amount of global memory on device: 1073545216 Bytes
	 Maximum amount of shared memory a single block may use: 49152 Bytes
	 Number of 32-bit registers available per block 32768
	 Number of threads in a warp: 32
	 Maximum pitch allowed for memory copies: 2147483647 Bytes
	 Maximum number of threads that a block may contain: 1024
	 Maximum number of threads allowed along each block dimension: 1024 x 1024 x 64
	 Number of blocks allowed along each grid dimension: 65535 x 65535 x 65535
	 Amount of available constant memory: 65536 Bytes
	 Major revision of the device's compute capability: 2
	 Minor revision of the device's compute capability: 1
	 Device's requirement for texture alignment: 512
	 Device can simultaneously perform a cudaMemcpy() and kernel execution: yes
	 Number of multiprocessors on device: 2
	 Runtime limit for kernels executed on this device: yes
	 Device is an integrated GPU (part of the chipset and not a discrete GPU): no
	 Device can map host memory into the CUDA device address space: yes
	 Device's computing mode: default, exclusive, or prohibited: 0
	 Maximum size supported for 1D textures: 65536
	 Maximum dimensions supported for 2D textures: 65536 x 65535
	 Maximum dimensions supported for 3D textures: 2048 x 2048 x 2048
	 Device supports executing multiple kernels within the same context simultaneously: yes
	

Το πρόθεμα __global__ δηλώνει ότι η συνάρτηση που ακολουθεί θα κληθεί μεν από την CPU αλλά θα εκτελεστεί από τους πυρήνες της GPU. Οι σχετικές λέξεις κλειδιά είναι οι παρακάτω

Είναι δυνατόν να χρησιμοποιηθούν μαζί τα δύο προθέματα __device__ και __host__ σε μια δήλωση συνάρτησης, πχ.
__device__ __host__ void f() { ... }
Στην περίπτωση αυτή ο μεταφραστής δημιουργεί δύο εκδόσεις της συνάρτησης f(): η μια θα καλείται μόνο από την CPU και θα εκτελείται στην CPU ενώ η δεύτερη θα καλείται από την GPU και θα εκτελείται στην GPU.

Σε μια συνάρτηση πυρήνα μπορούν να γίνουν οποιεσδήποτε πράξεις μεταξύ μεταβλητών, ανάγνωση και εγγραφή από τη μνήμη της συσκευής, ενώ μπορούν επίσης να χρησιμοποιηθούν όλα τα είδη βρόχων και γενικά όλες οι εντολές της C. Υπάρχουν ωστόσο κάποιοι περιορισμοί:

• Δεν μπορούν να χρησιμοποιηθούν συναρτήσεις που χρησιμοποιούν συσκευές I/O (οθόνη, πληκτρολόγιο, κάρτα δικτύου, κλπ). Για παράδειγμα, δεν μπορεί να χρησιμοποιηθεί η εντολή printf() η οποία γράφει στην οθόνη, ούτε η εντολή scanf() η οποία διαβάζει από το πληκτρολόγιο.
• Oι εντολές σε μια συνάρτηση πυρήνα δεν μπορούν να προσπελάσουν απ' ευθείας την κύρια μνήμη της CPU (host). Αν θέλουμε να προσπελάσουμε τα δεδομένα της κύριας μνήμης πρέπει πρώτα να τα μεταφέρουμε στη μνήμη της GPU χρησιμοποιώντας την εντολή cudaMemcpy().
• Δεν υποστηρίζεται ο τύπος δεδομένων double. Σε μια συνάρτηση πυρήνα, οι μεταβλητές αυτού του τύπου μετατρέπονται αυτομάτως κατά την μετάφραση σε float.

Οι παραπάνω περιορισμοί ισχύουν για συναρτήσεις τύπου __device__ και __global__ οι οποίες εκτελούνται στην GPU. Προφανώς δεν ισχύουν για συναρτήσεις τύπου __host__ οι οποίες εκτελούνται στην CPU όπως οι κλασικές συναρτήσεις C / C++.

Υπάρχουν διάφορες πηγές πληροφόρησης σχετικά με την CUDA.

• Δείτε το CUDA Toolkit Documentation στο CUDA Developer Zone
• Δείτε πληροφορίες σχετικά με τις συμβατές κάρτες γραφικών https://developer.nvidia.com/cuda-gpus
• Getting Started with CUDA (Univ. Gratz)
• Supercomputing blog - CUDA tutorials