Число π (пи) является одной из важнейших математических констант, часто применяемых в математике, физике и инженерии. Существует множество методов для приближенного вычисления значения π. Одним из таких методов является численное интегрирование, где значение π можно рассчитать как площадь под определенной кривой. Метод прямоугольников позволяет получить приближенное значение интеграла путем деления области под кривой на малые прямоугольники и суммирования их площадей.
С развитием вычислительной техники и методов параллельных вычислений стало возможным значительно ускорить такие расчеты. В данной работе будет рассмотрено последовательное и несколько параллельных решений для вычисления числа π, используя методы и библиотеки для параллельных вычислений: OpenMP, MPI и CUDA.
Постановка задачи
Цель данной работы – реализовать программы на языке C для вычисления числа π методом прямоугольников с использованием численного интегрирования. Для этого требуется:
Разработать последовательную программу вычисления числа π методом прямоугольников.
Реализовать параллельный вариант программы с использованием OpenMP для многопоточной обработки на процессоре.
Реализовать параллельный вариант программы с использованием MPI, распределяя вычисления по нескольким узлам или процессам.
Реализовать параллельный вариант программы на языке CUDA C для выполнения вычислений на графическом процессоре.
Провести сравнительный анализ времени выполнения между последовательной и параллельными реализациями.
Реализация вычислений.
Последовательный вариант
Листинг программы:
#include <stdio.h>
#include <time.h>
int main() {
int n = 1000000000;
double pi, h, sum = 0.0;
h = 1.0 / (double)n;
clock_t start = clock(); // Старт замера времени
for (int i = 1; i <= n; i++) {
double x = h * ((double)i - 0.5);
sum += 4.0 / (1.0 + x * x);
}
pi = h * sum;
clock_t end = clock(); // Конец замера времени
double time_spent = (double)(end - start) / CLOCKS_PER_SEC;
printf("Sequential pi is approximately %.16f\n", pi);
printf("Time taken: %.6f seconds\n", time_spent);
return 0;
}
Результат работы:
Рис.1. Последовательный вариант.
Первыми 16 числами после 3 числа пи является: 3.1415926535897932
Как видно вычисления после 12 знака теряют точность из-за особенностей метода вычислений, в целом точность терпимая, а вот время выполнения оставляет желать лучшего.
Параллельный вариант с использованием OpenMP
Листинг программы:
#include <stdio.h>
#include <omp.h>
int main() {
int n = 1000000000;
double pi, h, sum = 0.0;
h = 1.0 / (double)n;
double start = omp_get_wtime(); // Старт замера времени
omp_set_num_threads(1);
#pragma omp parallel for reduction(+:sum)
for (int i = 1; i <= n; i++) {
double x = h * ((double)i - 0.5);
sum += 4.0 / (1.0 + x * x);
}
pi = h * sum;
double end = omp_get_wtime(); // Конец замера времени
printf("Parallel (OpenMP) pi is approximately %.16f\n", pi);
printf("Time taken: %.6f seconds\n", end - start);
return 0;
}
Результат работы:
Рис.2. Вычисления используя openMP.
Программа была запущена с 8, 4 и 2 потоками соответственно, скорость вычислений с OpenMP используя даже 1 поток сильно быстрее чем последовательный вариант, вероятно данная библиотека использует какие-либо улучшения для столь быстрых вычислений. Стоит еще упомянуть, что чем больше потоков, тем точнее получался результат.
Параллельный вариант используя MPI
Листинг программы:
#include <stdio.h>
#include <mpi.h>
int main(int argc, char* argv[]) {
int n = 1000000000;
int rank, size, i;
double pi, h, sum = 0.0, local_sum = 0.0;
MPI_Init(&argc, &argv); // Инициализация MPI
MPI_Comm_rank(MPI_COMM_WORLD, &rank); // Получение номера процесса
MPI_Comm_size(MPI_COMM_WORLD, &size); // Получение общего числа процессов
h = 1.0 / (double)n;
int local_n = n / size; // Количество итераций на процесс
int start = rank * local_n + 1;
int end = start + local_n;
double start_time = MPI_Wtime(); // Старт замера времени
for (i = start; i < end; i++) {
double x = h * ((double)i - 0.5);
local_sum += 4.0 / (1.0 + x * x);
}
MPI_Reduce(&local_sum, &sum, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
if (rank == 0) {
pi = h * sum;
double end_time = MPI_Wtime(); // Конец замера времени
printf("Parallel (MPI) pi is approximately %.16f\n", pi);
printf("Time taken: %.6f seconds\n", end_time - start_time);
}
MPI_Finalize(); // Завершение MPI
return 0;
}
Результат выполнения:
Рис.3. Параллельный вариант с MPI
Скорость выполнения растет от количества потоков, однако точность вычислений никак не зависит от количества потоков, при каждом прогоне программы получается близкий, но никак не точный результат, обусловленный методом исчислений.
Параллельный вариант используя CUDA
#include <stdio.h>
#include <cuda_runtime.h>
__global__ void calculate_pi(double *partial_sums, double h, int n) {
extern __shared__ double sdata[];
int tid = threadIdx.x;
int i = blockIdx.x * blockDim.x + threadIdx.x + 1;
double x = h * ((double)i - 0.5);
// Вычисляем частичное значение для каждого потока
sdata[tid] = (i <= n) ? 4.0 / (1.0 + x * x) : 0.0;
__syncthreads();
// Параллельная редукция для суммирования значений в блоке
for (int s = blockDim.x / 2; s > 0; s >>= 1) {
if (tid < s) {
sdata[tid] += sdata[tid + s];
}
__syncthreads();
}
// Записываем результат блока в массив частичных сумм
if (tid == 0) {
partial_sums[blockIdx.x] = sdata[0];
}
}
int main() {
int n = 1000000000; // Количество шагов
double h = 1.0 / (double)n;
// Вычисление числа блоков и потоков
int blockSize = 1024;
int numBlocks = (n + blockSize - 1) / blockSize;
// Выделяем память для хранения частичных сумм
double *d_partial_sums;
cudaMalloc((void**)&d_partial_sums, numBlocks * sizeof(double));
// Создаем события для замера времени
cudaEvent_t start, stop;
cudaEventCreate(&start);
cudaEventCreate(&stop);
// Запуск события start
cudaEventRecord(start);
// Запускаем вычисление на GPU
calculate_pi<<<numBlocks, blockSize, blockSize * sizeof(double)>>>(d_partial_sums, h, n);
cudaDeviceSynchronize();
// Запуск события stop
cudaEventRecord(stop);
cudaEventSynchronize(stop);
// Копируем частичные суммы обратно на хост
double *h_partial_sums = (double*)malloc(numBlocks * sizeof(double));
cudaMemcpy(h_partial_sums, d_partial_sums, numBlocks * sizeof(double), cudaMemcpyDeviceToHost);
cudaFree(d_partial_sums);
// Суммируем частичные результаты на CPU
double pi = 0.0;
for (int i = 0; i < numBlocks; i++) {
pi += h_partial_sums[i];
}
pi *= h; // Окончательное значение π
free(h_partial_sums);
// Подсчет времени
float milliseconds = 0;
cudaEventElapsedTime(&milliseconds, start, stop);
// Вывод результата и времени
printf("Parallel (CUDA) pi is approximately %.16f\n", pi);
printf("Time taken: %.6f seconds\n", milliseconds / 1000.0);
// Очистка событий
cudaEventDestroy(start);
cudaEventDestroy(stop);
return 0;
}
Результат выполнения:
Рис.4. Выполненная Cuda программа.
Вариант достаточно быстр и точен, во всяком случае программа написанная на CUDA выдавала на 1 правильный знак больше, чем MPI, а так же выполнилась быстрее чем OpenMP, вероятно можно и сильно быстрее, но при попытке увеличения кол-ва блоков программа перестала выдавать результат.