Для иллюстрации возможности технологии CUDA рассмотрим для начала решение задачи Эйлера, написанное на языке программирования C с использованием только возможностей центрального процессора компьютера. Подробно эту задачу мы обсуждали в одном из предыдущих параграфов, поэтому просто приведем исходный текст программы, реализующий численный метод ее решения.

#include <stdio.h>
#include <math.h>
#include <sys/types.h>
#include <time.h>
#include <stdlib.h>
 
#define xmax 49152
#define ymax 1000
 
//Подпрограмма задания начальных значений
void initial(float *f[ymax]) {
	int x,y;
    //Задаем начальные значения массива
	for(x=0;x<xmax;x++) {
		for(y=0;y<ymax;y++) {
			f[x][y]=0.0;
	}}
}
 
//Подпрограмма задания граничных условий
void boundary(float *f[ymax]) {
	int x,y;
	//Задаем граничные значения массива
	for(y=0;y<ymax;y++) {
		f[0][y]=0.0;
		f[xmax-1][y]=0.0;
	}
	for(x=0;x<xmax;x++) {
		f[x][0]=sin((x*1.0)/(xmax/2));
		f[x][ymax-1]=0.0;
	}
	f[xmax/2][ymax-1]=-5.0;
}
 
//Подпрограмма итерационного процесса
float iter(float *a[ymax], float *b[ymax]) {
	float dt = 0.01;
	float dx = 0.5;
	float dy = 0.5;
	float dff, df;
	int x,y;
 
	df=0.0;
	for(x=1;x<xmax-1;x++) {
	for(y=1;y<ymax-1;y++) {
		dff = dt * (
			((a[x+1][y]-(2*a[x][y])+a[x-1][y])/(dx*dx))
			+
			((a[x][y+1]-(2*a[x][y])+a[x][y-1])/(dy*dy))
			);
		b[x][y]=a[x][y]+dff;
		if( df < fabs(dff) ) { df=fabs(dff); }
 
	}}
	return df;
}
 
void main() {
	printf("Start... [%dx%d]\n",xmax,ymax);
	float df = 0.0; int n=0,i;
	float elapsed;
	clock_t start, end;
 
	float **f0 = malloc(xmax * sizeof(float *));
	float **f1 = malloc(xmax * sizeof(float *));
	for(i=0;i<xmax;i++) f0[i] = malloc(ymax * sizeof(float));
	for(i=0;i<xmax;i++) f1[i] = malloc(ymax * sizeof(float));
 
	initial(f0); initial(f1);
	boundary(f0); boundary(f1);
 
	printf("Begin calculations...\n");
 
	start=clock();
	do {
		if ( ((n/2)*2) == n ) {
			df = iter(f0,f1);
		}else{
			df = iter(f1,f0);
		}
		n++;
		//printf(" [%d] Diff: %f\n",n,df);
	} while ( df > 0.01 );
	end=clock();
	elapsed=((float) (end - start)) / CLOCKS_PER_SEC;
 
	for(i=0;i<xmax;i++) free((void *) f0[i]);
	for(i=0;i<xmax;i++) free((void *) f1[i]);
	free((void *) f0); free((void *) f1);
 
	printf("Calculation time: %2lf sec.\n",elapsed);
	printf("Iteration count: %d\n",n);
	printf("Averege Iteration time: %.2lf sec.\n",elapsed/n);
}
 
Загрузить исходный код программы

Вычисления мы будем проводить на матрице размера 49152x1000 элементов. Как и раньше особого смысла ни в начальных, ни в граничных условиях задачи нет, поскольку нас интересует в данном случае не физический смысл, а сам процесс счета, точнее его скоростные характеристики.

Поскольку графический процессор видеокарты по сути является мультипроцессором, каждое ядро которого способно выполнять несколько потоков команд обновременно, то логично будет разнести внутренний цикл итерации по оси Y разностной сетки(строка программы 45) в отдельные параллельно исполняемые блоки, получив значительный прирост быстродействия. Этим мы и займемся в следующем параграфе.