Нижегородский государственный университет им. Н.И.Лобачевского Факультет Вычислительной математики и кибернетики Лабораторная работа 2: Вычисление определенного интеграла Параллельные численные методы Козинов Е.А., Мееров И.Б., Сысоев А.В. Кафедра математического обеспечения ЭВМ

Содержание Введение Цель работы Тестовая инфраструктура Вычисление определенного интеграла Программная реализация Задания для самостоятельной работы Литература Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 2

1. Приступаем к работе Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 3

Введение В лабораторной работе рассматривается известная математическая задача – вычисление значений определенных интегралов. –Рассматривается случай, когда задача не решается аналитически. –Для численного интегрирования используется метод прямоугольников. В работе изучаются подходы к распараллеливанию метода прямоугольников, идеи по алгоритмической оптимизации, приводящие к уменьшению времени вычислений, вопросы параллельной отладки. В работе широко используется пакет программных инструментов Intel Parallel Studio (компилятор, отладчик, профилировщик). Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 4

Цель работы Изучение принципов написания высокопроизводительных реализаций алгоритмов с использованием современных компилирующих и отладочных средств Изучение общей схемы алгоритма численного расчета определенного интеграла, обсуждение программной реализации и возможных подходов к ее оптимизации по скорости. Написание нескольких программных реализаций численного расчета определенного интеграла и сравнение их производительности. Демонстрация использования инструментов пакета Intel Parallel Studio в процессе реализации и оптимизации программного кода. Изучение подходов, позволяющих увеличить производительность программных реализаций алгоритмов. Сравнение на данном учебном примере нескольких подходов к распараллеливанию. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 5

Тестовая инфраструктура Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 6 Процессор2 четырехъядерных процессора Intel Xeon E5520 (2.27 GHz) Память16 Gb Операционная система Microsoft Windows 7 Среда разработкиMicrosoft Visual Studio 2008 Компилятор, профилировщик, отладчик Intel Parallel Studio SP1

2. Интегрирование по методу прямоугольников Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 7

Численное интегрирование по методу прямоугольников и трапеций Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 8

Численное интегрирование по методу прямоугольников и трапеций Для случая деления отрезка интегрирования на равные части и вычисления функции в центре отрезков: –где N – количество отрезков интегрирования, а h = (b – a) / N. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 9

Численное интегрирование по методу прямоугольников и трапеций Погрешность вычислений квадратурных формул методов трапеций и прямоугольников: –где Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 10

Численное интегрирование по методу прямоугольников для двумерной функции –где N – количество участков интегрирования по оси x, M – количество участков интегрирования по оси y, h 1 = (b 1 – a 1 ) / N и h 2 = (b 2 – a 2 ) / M. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 11

Тестовая функция В качестве учебного примера рассмотрим тестовую функцию вида: Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 12

3. Программная реализация Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 13

Задание #1. Последовательная версия. Создание нового проекта Новый проект 01_Reference, консольное приложение. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 14

Задание #1. Последовательная версия. Базовая реализация (1) Создайте файл main_s.cpp, в который будет помещена реализация. Основная функция должна иметь структуру: int main () { int i; // переменная цикла double time; // время проведенного эксперимента double res; // значение вычисленного интеграла double min_time; // минимальное время работы // реализации алгоритма double max_time; // максимальное время работы // реализации алгоритма double avg_time; // среднее время работы // реализации алгоритма int numbExp = 10; // количество запусков программы … Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 15

Задание #1. Последовательная версия. Базовая реализация (2) // первый запуск min_time = max_time = avg_time = experiment(&res); // оставшиеся запуски for(i = 0; i < numbExp - 1; i ++) { time = experiment(&res); avg_time += time; if(max_time < time) max_time = time; if(min_time > time) min_time = time; } // вывод результатов эксперимента printf("integral value : %lf; \n", res); printf("execution time : %lf; %lf; %lf \n", avg_time / numbExp, min_time, max_time); return 0; } Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 16

Задание #1. Последовательная версия. Базовая реализация (3) Реализуйте функцию однократного вычисления интеграла experiment() double experiment(double *res) { double stime, ftime; // время начала и конца расчета double a1 = 0.0 ; // левая граница интегрирования по x double b1 = 16.0; // правая граница интегрирования по x double a2 = 0.0 ; // левая граница интегрирования по y double b2 = 16.0; // правая граница интегрирования по y double h = 0.001; // шаг интегрирования stime = omp_get_wtime( ); integral(a1, b1, a2, b2, h, res); // вычисление интеграла ftime = omp_get_wtime( ); return (ftime - stime); } Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 17

Задание #1. Последовательная версия. Базовая реализация (4) Реализуйте функцию integral() для вычисления приближенного значения интеграла тестовой функции. void integral(const double a1, const double b1, const double a2, const double b2, const double h, double *res) { int i, j, n1, n2; double sum; // локальная переменная // для подсчета интеграла double x; // координата точки сетки по оси x double y; // координата точки сетки по оси y // количество точек сетки интегрирования // n1 - по координате x // n2 - по координате y n1 = (int)((b1 - a1) / h); n2 = (int)((b2 - a2) / h); … Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 18

Задание #1. Последовательная версия. Базовая реализация (5) Реализуйте функцию integral() для вычисления приближенного значения интеграла тестовой функции. sum = 0.0; for(i = 0; i < n1; i++) { for(j = 0; j < n2; j++) { // вычисление координат точек x = a1 + i * h + h / 2; y = a2 + j * h + h / 2; // вычисление интеграла sum += ((exp(sin(x * PI)*cos(y * PI)) + 1) / ((b1 - a1) * (b2 - a2))) * h * h; } *res = sum; } Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 19

Задание #1. Последовательная версия. Запуск приложения Соберите получившийся код (команда BuildBuild Solution) и запустите его на выполнение. Сравните результат вычисления интеграла с представленным на слайде Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 20

Задание #2. Последовательная версия. Применение компилятора Intel C/C++ Compiler На данном этапе лабораторной работы предлагается использовать оптимизирующий компилятор Intel Windows C/C++ Compiler, входящий в состав продукта Intel Parallel Studio (пакет Intel Parallel Composer). –В окне Solution Explorer выберите проект 01_Reference и выполните команду контекстного меню Intel Parallel ComposerUse Intel C++…. –Перестройте решение и запустите полученный исполняемый файл. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 21

Сравнение времени выполнения алгоритма численного интегрирования с использованием разных компиляторов. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 22 Задание #2. Последовательная версия. Сравнение эффективности компиляторов

Задание #3. Параллельная версия. Схема распараллеливания по координате x Используя средства OpenMP, реализуйте первую параллельную версию подсчета интеграла с разделением сетки интегрирования по столбцам (по координате x). Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 23

Задание #3. Параллельная версия. Создание нового проекта Создайте в рамках решения 02_Integral новый проект с названием 02_Integral_col. В проекте создайте файл main_col.cpp и скопируйте в него код из файла main_s.cpp проекта 01_Reference. Выполните переход к использованию компилятора Intel C++. Настройте в свойствах проекта использование OpenMP. –В дереве Configuration Properties перейдите к разделу C/C++Language и в поле OpenMP Support справа выберите вариант: Generate Parallel Code (/openmp, equiv. to /Qopenmp). Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 24

Задание #3. Параллельная версия. Программная реализация Согласно схеме разбиения данных, необходимо распараллеливать внешний цикл функции integral(). –Примените директиву OpenMP omp parallel for, позволяющую создать потоки и распределить итерации цикла между ними. //разделение точек сетки интегрирования по оси x #pragma omp parallel for for(i = 0; i < n1; i++) { for(j = 0; j < n2; j++) { //вычисление координат точки x = a1 + i*h + h / 2; y = a2 + j*h + h / 2; //вычисление интеграла sum += ((exp(sin(x * PI) * cos(y * PI)) + 1) / ((b1 - a1) * (b2 - a2))) * h * h; } Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 25

Задание #3. Параллельная версия. Запуск приложения Соберите проект 02_Integral_col и запустите на исполнение несколько раз. –Убедитесь, что на многоядерной/многопроцессорной системе результат работы функции подсчета интеграла существенно отличается от последовательного и, кроме того, меняется от запуска к запуску. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 26

Задание #3. Параллельная версия. Поиск ошибок с Intel Parallel Inspector (1) В чем же дело? Ответить на этот вопрос может помочь инструмент Intel Parallel Inspector из пакета Intel Parallel Studio. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 27

Задание #3. Параллельная версия. Поиск ошибок с Intel Parallel Inspector (2) Схема поиска ошибок многопоточности с использованием Intel Parallel Inspector: –собрать Debug версия программы; –выбрать в IPI режим работы Threading errors; –выбрать режим анализа Analysis level: Extreme; –нажать Run Analysis. Найдите ошибки многопоточности в проекте 02_Integral_col Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 28

Задание #3. Параллельная версия. Поиск ошибок с Intel Parallel Inspector (3) «Гонки» данных: Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 29

Задание #3. Параллельная версия. Исправление «гонок» данных (1) «Гонки» данных связаны с доступом к следующим переменным: –x и y – координаты точки, где необходимо вычислить функцию; –j – переменная внутреннего цикла; –sum – переменная для накопления суммы интеграла. Первые три переменные необходимо сделать локальными для потоков, а по переменной sum организовать редукцию данных с операцией суммирования. –Для локализации данных в директиве parallel можно использовать параметр private( ). –Для редукции данных в директиве omp for используется параметр reduction( : ). Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 30

Задание #3. Параллельная версия. Исправление «гонок» данных (2) Исправленная версия кода: // Сетка интегрирования делится по столбцам // разделение точек сетки интегрирования по оси x #pragma omp parallel for private (x, y, j) / reduction(+: sum) for(i = 0; i < n1; i++) { for(j = 0; j < n2; j++) { //вычисление координат точки x = a1 + i*h + h / 2; y = a2 + j*h + h / 2; //вычисление интеграла sum += ((exp(sin(x * PI) * cos(y * PI)) + 1) / ((b1 - a1) * (b2 - a2))) * h * h; } Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 31

Задание #3. Параллельная версия. Повторный поиск ошибок с IPI Проанализируйте полученный код с помощью Intel Parallel Inspector. –Результат анализа должен показать, что в полученной версии кода гонок данных нет. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 32

Задание #3. Параллельная версия. Сравнение результатов вычислений Результаты последовательной и параллельной версии не совпали. Всегда ли это говорит об ошибке? Результаты последовательной и параллельной версии совпали. Всегда ли это говорит об отсутствии ошибок? Инструмент показывает наличие ошибок. Всегда ли это правда? Инструмент показывает отсутствие ошибок. Всегда ли это правда? Алгоритмически порядок вычислений в посл. и паралл. версии одинаков. Всегда ли это будет так в машинной реализации? Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 33

Задание #3. Параллельная версия. Сравнение времени работы Сравнение времени численного интегрирования для последовательной и параллельной реализации Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 34

Задание #4. Параллельная версия. Схема распараллеливания по координате y Сетку интегрирования можно разделять не только по столбцам, но и по строкам. –Схематичное изображение метода распараллеливания по строкам представлено на слайде. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 35

Задание #4. Параллельная версия. Создание нового проекта Создайте в рамках решения 02_Integral новый проект с названием 02_Integral_row. Выполните переход к использованию компилятора Intel C++ и в свойствах проекта включите поддержку использования OpenMP. В рамках созданного проекта реализуйте параллельный алгоритм численного интегрирования с разделением сетки по строкам. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 36

Задание #4. Параллельная версия. Программная реализация Параллельная реализация с разделением сетки интегрирования по строкам (по координате y): //Сетка интегрирования делится по строкам for(i = 0; i < n1; i++) { //разделение точек сетки интегрирования по оси y #pragma omp parallel for private (x, y) / reduction(+: sum) for(j = 0; j < n2; j++) { //вычисление координат точки y = a2 + j*h + h / 2; x = a1 + i*h + h / 2; //вычисление интеграла sum += ((exp(sin(x * PI) * cos(y * PI)) + 1) / ((b1 - a1) * (b2 - a2))) * h * h; } Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 37

Задание #4. Параллельная версия. Запуск приложения Пересоберите получившийся код (команда BuildRebuild Solution) и запустите его на выполнение. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 38

Задание #4. Параллельная версия. Сравнение с предыдущими реализациями Сравнение времени численного интегрирования для последовательной и параллельных реализаций Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 39

Задание #5. Параллельная версия. Блочная схема распараллеливания Еще один популярный способ геометрического разделения данных – разделения данных на блоки. –Схематичное изображение параллельного алгоритма численного интегрирования с блочным разделением данных представлено на слайде. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 40

Задание #4. Параллельная версия. Создание нового проекта Создайте в рамках решения 02_Integral новый проект с названием 02_Integral_block. Выполните переход к использованию компилятора Intel C++ и в свойствах проекта включим поддержку использования OpenMP. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 41

Задание #4. Параллельная версия. Реализация блочной схемы Блочная схема легко реализуется с помощью использования вложенного параллелизма в OpenMP. –Чтобы использовать вложенный параллелизм, необходимо сначала включить его поддержку с помощью вызова библиотечной функции: // включаем возможность использования вложенного // параллелизма omp_set_nested(true); Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 42

Задание #4. Параллельная версия. Программная реализация Реализуйте параллельный алгоритм численного интегрирования с разделением сетки интегрирования на блоки. omp_set_nested(true); #pragma omp parallel for for (i = 0; i < n1; i++) { //разделение точек сетки интегрирования по оси y #pragma omp parallel for private (x, y) reduction(+: sum) for(j = 0; j < n2; j++) { //вычисление координат точки x = a1 + i * h + h / 2; y = a2 + j * h + h / 2; //вычисление интеграла sum += ((exp(sin(x * PI) * cos(y * PI)) + 1) / ((b1 - a1) * (b2 - a2))) * h * h; } Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 43

Задание #4. Параллельная версия. Запуск приложения Пересоберите получившийся код (команда BuildRebuild Solution) и запустите его на выполнение. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 44

Задание #4. Параллельная версия. Сравнение с предыдущими реализациями Сравнение времени численного интегрирования для последовательной и параллельных реализаций Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 45

Задание #5. Последовательная версия. Поиск участков кода для оптимизации (1) Как еще можно увеличить производительность программы? –В общем случае, чтобы ответить на данный вопрос, необходимо найти участки кода, которые нуждаются в оптимизации и занимают основное время вычислений. –Для поиска таких мест в пакете Intel Parallel Studio имеется инструмент Intel Parallel Amplifier. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 46

Задание #5. Последовательная версия. Поиск участков кода для оптимизации (2) Для профилировки необходимо выбрать первый пункт – поиск «горячих точек». –Будем профилировать исходную последовательную версию алгоритма. Запустите Intel Parallel Amplifier. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 47

Задание #5. Последовательная версия. Использование предварительных вычислений Из результатов анализа видно, что основное время в работе программы занимает вычисление математических функций. Значения математических функций sin и cos вычисляются много раз. –При этом точки, в которых вычисляются эти функции, очень часто повторяются. Данные значения можно вычислить предварительно. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 48

Задание #5. Последовательная версия. Создание нового проекта Создайте в рамках решения 02_Integral новый проект с названием 03_AlgOptV1. Выполните переход к использованию компилятора Intel C++. Добавьте в проект файл main_alg1.cpp. Скопируйте в него весь код из main_s.cpp кроме функции вычисления интеграла. Далее необходимо реализовать функцию integral() использующую предварительные вычисления. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 49

Задание #5. Последовательная версия. Программная реализация (1) Объявите все требуемые переменные. void integral(const double a1, const double b1, const double a2, const double b2, const double h, double *res) { int i, j, n1, n2; double sum;// локальная переменная для подсчета интеграла double x; // координата точки сетки по оси x double y; // координата точки сетки по оси y double *sinx; // значение sin(x * pi) double *cosy; // значение cos(y * pi) //количество точек сетки интегрирования n1 = (int)((b1 - a1) / h); n2 = (int)((b2 - a2) / h); … Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 50

Задание #5. Последовательная версия. Программная реализация (2) Предварительно вычислите значения функции sin, используемые в последующих расчетах. // вычисление значений sin(x * pi) sinx = new double [n1]; for(i = 0; i < n1; i++){ x = a1 + i*h + h / 2; sinx[i] = sin(x * PI); } Аналогичным образом вычислите значения функции cos. // вычисление значений cos(y * pi) cosy = new double [n2]; for(j = 0; j < n2; j++){ y = a2 + j*h + h / 2; cosy[j] = cos(y * PI); } Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 51

Задание #5. Последовательная версия. Программная реализация (3) Последним шагом вычислите значение интеграла с использованием результатов предвычислений. //вычисление интеграла sum = 0.0; for(i = 0; i < n1; i++) { for(j = 0; j < n2; j++) { //вычисление интеграла //(значение sin и cos уже подсчитаны) sum += ((exp(sinx[i] * cosy[j]) + 1) / ((b1 - a1) * (b2 - a2))) * h * h; } *res = sum; delete [] sinx; delete [] cosy; } Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 52

Задание #5. Последовательная версия. Запуск приложения Пересоберите получившийся код (команда BuildRebuild Solution) и запустите его на выполнение. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 53

Задание #5. Последовательная версия. Повторный анализ «горячих» точек Проанализируйте полученный код на предмет нахождения «горячих» точек. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 54

Задание #5. Последовательная версия. Сравнение с предыдущими реализациями Сравнение времени численного интегрирования последовательных и параллельных реализаций Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 55

Задание #5. Последовательная версия. Анализ разработанного решения Недостаток предыдущей реализации заключается в том, что при увеличении области интегрирования резко возрастают расходы на память для хранения подсчитанных значений. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 56

Задание #6. Последовательная версия. Использование буферизации Необходимо разбить область вычислений на прямоугольники равных размеров, в точках которых будут выполняться предварительные вычисления. Для каждого прямоугольника подсчитать интеграл и полученные значения суммировать. Размер буфера обычно подбирают таким образом, чтобы количество точек сетки интегрирования делилось на него без остатка. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 57

Задание #6. Последовательная версия. Создание нового проекта Создайте в рамках решения 02_Integral новый проект с названием 04_Buf. Выполните переход к использованию компилятора Intel C++. Создайте файл main_Buf.cpp и скопируйте в него код из файла main_alg1.cpp проекта 01_Reference. Модифицируйте код согласно предложенным улучшениям. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 58

Задание #6. Последовательная версия. Программная реализация (1) Объявите необходимые переменные. //размер буфера #define BLOCK_SIZE 2000 void integral(const double a1, const double b1, const double a2, const double b2, const double h, double *res) { int i, j, ii, jj, n1, n2, nb1, nb2; double sum;// локальная переменная для подсчета // интеграла double x; // координата точки сетки по оси x double y; // координата точки сетки по оси y double *sinx; // значение sin(x * pi) double *cosy; // значение cos(y * pi) Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 59

Задание #6. Последовательная версия. Программная реализация (2) Подсчитайте необходимые размеры сеток и проверьте тот факт, что выбран буфер правильного размера //количество точек сетки интегрирования n1 = (int)((b1 - a1) / h); n2 = (int)((b2 - a2) / h); //правильность размера блока assert((n1 % BLOCK_SIZE) == 0); assert((n2 % BLOCK_SIZE) == 0); //Вычисление количества точек в блоке // nb1 - по координате x // nb2 - по координате y nb1 = n1 / BLOCK_SIZE; nb2 = n2 / BLOCK_SIZE; Выделите необходимую память. sinx = new double [BLOCK_SIZE]; cosy = new double [BLOCK_SIZE]; Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 60

Задание #6. Последовательная версия. Программная реализация (3) Далее пройдите в цикле по всем блокам и вычислите значение интеграла. //проход по всем блокам for(ii = 0; ii < n1; ii += BLOCK_SIZE) { // вычисление значений sin(x * pi) for(i = 0; i < BLOCK_SIZE; i++) { x = a1 + i*h + ii + h / 2; sinx[i] = sin(x * PI); } /* for(i = 0; i < BLOCK_SIZE; i++) */ for(jj = 0; jj < n2; jj += BLOCK_SIZE) { // вычисление значений cos(y * pi) for(j = 0; j < BLOCK_SIZE; j++) { y = a2 + j*h + jj + h / 2; cosy[j] = cos(y * PI); } /* for(j = 0; j < BLOCK_SIZE; j++) */ //вычисление интеграла for(i = 0; i < BLOCK_SIZE; i++) { for(j = 0; j < BLOCK_SIZE; j++) { //вычисление интеграла sum += ((exp(sinx[i] * cosy[j]) + 1) / ((b1 - a1) * (b2 - a2))) * h * h; } // for(j = 0; j < BLOCK_SIZE; j++) } // for(i = 0; i < BLOCK_SIZE; i++) } // for(jj = 0; jj < n2; jj += BLOCK_SIZE) } // for(ii = 0; ii < n1; ii += BLOCK_SIZE) Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 61

Задание #6. Последовательная версия. Запуск приложения Пересоберите получившийся код (команда BuildRebuild Solution) и запустите его на выполнение. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 62

Задание #6. Последовательная версия. Сравнение с предыдущими реализациями Сравнение времени численного интегрирования последовательных и параллельных реализаций Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 63

Задание #7. Последовательная версия. Алгоритмическая оптимизация Предпоследняя оптимизация основана на математических особенностях вычисляемого интеграла, а именно на том факте, что функции sin и cos периодические. –Предположим, что количество точек в сетке интегрирования кратно размеру подобранного буфера. –В этом случае легко видеть, что итерации цикла по блокам будут вычислять одно и то же значение. Можно также заметить, что: Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 64

Задание #7. Последовательная версия. Создание нового проекта Реализуйте алгоритмическое улучшение в проекте 05_AlgOptV2, с использованием техник, рассмотренных в предыдущих разделах. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 65

Задание #7. Последовательная версия. Программная реализация (1) //количество точек в периоде интегрирования npi = (int)(2.0 / h); //правильность размера блока assert((n1 % npi) == 0); assert((n2 % npi) == 0); // вычисление значений sin(x * pi) for(i = 0; i < npi; i++) { x = a1 + i*h + h / 2; sinx[i] = sin(x * PI); } // вычисление значений cos(y * pi) for(j = 0; j < npi; j++) { y = a2 + j*h + h / 2; cosy[j] = cos(y * PI); } Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 66

//вычисление интеграла sum = 0.0; for(i = 0; i < npi; i++) { for(j = 0; j < npi; j++) { //вычисление интеграла sum += (exp(sinx[i] * cosy[j]) / ((b1 - a1) * (b2 - a2))) * h * h; } *res = sum * (n1 / npi) * (n2 / npi) + 1; Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 67 Задание #7. Последовательная версия. Программная реализация (2)

Пересоберите получившийся код (команда BuildRebuild Solution) и запустите его на выполнение. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 68 Задание #7. Последовательная версия. Запуск приложения

Сравнение времени численного интегрирования последовательных реализаций Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 69 Задание #7. Последовательная версия. Сравнение с предыдущей реализацией

Задание #8. Параллельная версия оптимизированного алгоритма Реализуйте параллельную версию алгоритмически улучшенной реализации в проекте 06_AlgOptV2par, используя идеи распараллеливания предыдущих версий. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 70

Задание #8. Параллельная версия оптимизированного алгоритма. Реализация (1) В начале, объявите необходимые переменные. void integral(const double a1, const double b1, const double a2, const double b2, const double h, double *res) { int i, j, n1, n2, npi; double sum;// локальная переменная для подсчета интеграла double x; // координата точки сетки по оси x double y; // координата точки сетки по оси y double *sinx; // значение sin(x * pi) double *cosy; // значение cos(y * pi) Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 71

Задание #8. Параллельная версия оптимизированного алгоритма. Реализация (2) Подсчитайте размеры сетки интегрирования и выделите необходимую память. //количество точек сетки интегрирования // n1 - по координате x // n2 - по координате y n1 = (int)((b1 - a1) / h); n2 = (int)((b2 - a2) / h); //количество точек в периоде интегрирования npi = (int)(2.0 / h); //правильность размера блока assert((n1 % npi) == 0); assert((n2 % npi) == 0); sinx = new double [npi]; cosy = new double [npi]; Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 72

Задание #8. Параллельная версия оптимизированного алгоритма. Реализация (3) Параллельно вычислите необходимые значения функций sin и cos. // вычисление значений sin(x * pi) #pragma omp parallel for private(x) for(i = 0; i < npi; i++) { x = a1 + i*h + h / 2; sinx[i] = sin(x * PI); } // вычисление значений cos(y * pi) #pragma omp parallel for private(y) for(j = 0; j < npi; j++) { y = a2 + j*h + h / 2; cosy[j] = cos(y * PI); } Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 73

Задание #8. Параллельная версия оптимизированного алгоритма. Реализация (4) Далее распараллельте основные вычислительные циклы и освободите выделенную память //вычисление интеграла sum = 0.0; #pragma omp parallel for private (x, y, j) / reduction(+: sum) for(i = 0; i < npi; i++) { for(j = 0; j < npi; j++) { //вычисление интеграла sum += (exp(sinx[i] * cosy[j]) / ((b1 - a1) * (b2 - a2))) * h * h; } *res = sum * (n1 / npi) * (n2 / npi) + 1; delete [] sinx; delete [] cosy; } Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 74

Задание #8. Параллельная версия. Запуск приложения Пересоберите получившийся код (команда BuildRebuild Solution) и запустите его на выполнение. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 75

Задание #8. Параллельная версия. Сравнение с последовательной версией Сравнение времени численного интегрирования последовательной и параллельной реализаций Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 76

Задания для самостоятельной работы… Реализуйте алгоритмы, использующие буфер, на случай, когда размер буфера не кратен количеству узлов сетки интегрирования. Предложите дальнейшие алгоритмические методы сокращения времени вычисления последней версии реализации алгоритма численного интегрирования. Рассмотреть возможность использования Intel MKL для векторного вычисления математических функций (функции VML). Рассмотреть другие квадратурные формулы. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 77

Использованные источники информации Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. – Бином. Лаборатория знаний, – 640c. Гергель В.П. Теория и практика параллельных вычислений. – Бином. Лаборатория знаний, – 424c. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 78

Авторский коллектив Козинов Евгений Александрович ассистент кафедры Математического обеспечения ЭВМ факультета ВМК ННГУ. Мееров Иосиф Борисович, к.т.н., доцент, зам. зав. кафедры Математического обеспечения ЭВМ факультета ВМК ННГУ. Сысоев Александр Владимирович, ассистент кафедры Математического обеспечения ЭВМ факультета ВМК ННГУ. Н.Новгород, 2010 г. Параллельные численные методы. Лабораторная работа 2: «Вычисление определенного интеграла» 79