Введение

Данная статья кратко описывает распараллеливание расчетов на вычислительных мощностях CPU и GPU. Перед тем как перейти к описанию самих алгоритмов, ознакомлю вас с поставленной задачей.

Необходимо смоделировать систему решения задач методом конечных разностей. С математической точки зрения это выглядит следующим образом. Дана некоторая конечная сетка:



Неизвестные значения сетки находятся по следующей формуле методом конечных разностей:

Распараллеливание на CPU

Для параллельных вычислений на CPU используются технологию Parallel, которая аналогична OpenMP. Parallel – внутренняя технология, используемая в языке C#, предоставляющая поддержку параллельных циклов и областей.

Параллельные вычисления на Parallel:

/* n*m размерность конечной сетки
T начальные значения конечной сетки
eps допустимая погрешность */
/* Так как процесс выполняется параллельно, то необходимо для хранения 
 ошибок и новых значений узлов использовать двумерные массивы, так 
 как выполняются вычисления без блокировок. Это тратит больше
 пространства памяти, но позволяет обойтись без блокировок.*/
private void Parallelization(int n, int m, float[,] T, float eps)
{
	int time; //Секундомер
	bool flag = false; //Условие завершения
	int interetion = 0; //Количество итераций
	float epsilint; // Наибольшая погрешность в конечной сетке
	float[,] count_eps = new float[n,m]; // Погрешность узла
	float[,] T_new = new float[n, m]; // Новые значения неизвестных эл-тов сетки
	time = Environment.TickCount; // Записывает время начала итераций
	do
	{
		epsilint = eps;
		Parallel.For(1, n-1, i =>
			{
			Parallel.For(1, m - 1, j =>
				{
					//Находим новое значение неизвестной
					T_new[i, j] = (T[i - 1, j] + T[i + 1, j] + T[i, j - 1] + T[i, j + 1]) / 4;
					//Вычисляем значение разницу между старым и новым значением
					count_eps[i, j] = Math.Abs(T_new[i, j] - T[i, j]);
					 //Проверяем погрешность полученного значения
					if (count_eps[i,j] > epsilint)
					{
						epsilint = count_eps[i,j];
					}
					T[i, j] = T_new[i, j];
				});
			});
		interetion++;
	}while(epsilint > eps || epsilint != eps); //повторяем пока погрешность не удовлетворяет условиям
	time = Environment.TickCount - time; //Записываем время окончания итераций
	Output(n, m, time, interetion, "OpenMP Parallezetion"); //Вывод
}

Распараллеливание на GPU

Для параллельных вычислений на GPU используется технология CUDA. CUDA – это архитектура параллельных вычислений от NVIDIA, позволяющая существенно увеличить вычислительную производительность благодаря использованию GPU.

Параллельные вычисления на CUDA:

/*Упрощаем работу с CUDA,
вылавливаем ошибки при выполнении команд*/
#define CUDA_DEBUG
#ifdef CUDA_DEBUG
#define CUDA_CHECK_ERROR(err)           if (err != cudaSuccess) {          printf("Cuda error: %s\n", cudaGetErrorString(err));    printf("Error in file: %s, line: %i\n", __FILE__, __LINE__);  }                 #else
#define CUDA_CHECK_ERROR(err)
#endif

/*Параллельные вычисления на GPU*/
__global__ void VectorAdd(float* inputMatrix, float* outputMatrix, int n, int m)
{
	int i = threadIdx.x + blockIdx.x * blockDim.x; //Индексация по столбцам
	int j = threadIdx.y  + blockIdx.y * blockDim.y; //Индексация по строкам
	if(i < n -1 && i > 0)
	{
		if(  j < m - 1 &&  j > 0)
			//Находим новое значение неизвестной
			outputMatrix[i * n + j] = (inputMatrix[(i - 1) * n + j ] + inputMatrix[(i + 1) * n + j] + inputMatrix[i * n + (j - 1)] + inputMatrix[i * n + (j + 1)])/4;
	}
}
/* n*m размерность конечной сетки
T начальные значения конечной сетки
eps допустимая погрешность */
/* В GPU все двумерные массивы передаются в виде одномерной строки,
тем самым заранее двумерный массив T был преобразован в одномерный*/
void OpenCL_Parallezetion(int n, int m, float *T, float eps)
{
	int matrixsize = n * m; // Размер памяти для CPU
	int byteSize = matrixsize * sizeof(float); // Размер памяти под GPU
	time_t start, end; // Время начала и конца итераций
	float time; // Время просчетов
	float* T_new = new float[matrixsize]; // Новые значения неизвестных эл-тов сетки
	float *cuda_T_in; // Старые значения неизвестных эл-тов сетки на GPU
	float *cuda_T_out; // Новые значения неизвестных эл-тов сетки на GPU
	CUDA_CHECK_ERROR(cudaMalloc((void**)&cuda_T_in, byteSize)); // Выделение памяти массива на GPU
	CUDA_CHECK_ERROR(cudaMalloc((void**)&cuda_T_out, byteSize));
	float epsilint; // Наибольшая погрешность в конечной сетке
	float count_eps; // Погрешность узла
	int interetaion = 0; //Количество итераций
	start = clock(); // Записывает время начала итераций
	dim3 gridsize = dim3(n,m,1); // Диапазон индексов (x,y,z) для GPU
	do{
		epsilint = eps;
		CUDA_CHECK_ERROR(cudaMemcpy(cuda_T_in, T, byteSize, cudaMemcpyHostToDevice)); // Копируем в память GPU
		VectorAdd<<< gridsize, m >>>(cuda_T_in, cuda_T_out, n, m);
		CUDA_CHECK_ERROR(cudaMemcpy(T_new, cuda_T_out, byteSize, cudaMemcpyDeviceToHost)); // Извлекаем из памяти GPU 
		for(int i = 1; i < n -1; i++)
		{
			for(int j = 1; j < m -1; j++)
			{
				//Вычисляем значение разницу между старым и новым значением
				count_eps = T_new[i* n + j ] - T[i* n + j];
				//Проверяем погрешность полученного значения
				if(count_eps > epsilint)
				{
					epsilint = count_eps;
				}
				T[i * n + j] = T_new[i * n + j];
			}
		}
		interetaion++;
	}while(epsilint > eps || epsilint != eps);//повторяем пока погрешность не удовлетворяет условиям
	end = clock(); //Записываем время окончания итераций
	time = (end - start); //Высчитываем время
	/*Освобождаем память массивов из памяти
	CPU и GPU*/
	free(T);
	free(T_new);
	cudaFree(cuda_T_in);
	cudaFree(cuda_T_out);
	Output(n, m, time, interetaion); //Вывод
}