работа Белорусский Государственный Университет Факультет радиофизики и электроники Кафедра информатики Разработка и анализ эффективности программы нахождения интегралов методом Монте-Карло на кластере MPI Выполнил: студент 4 группы 3 курса Строев В.Г. Научный руководитель: доцент Шпаковский Г.И. г. Минск, 2004 ОГЛАВЛЕНИЕ ВВЕДЕНИЕ ……………………………………………………………... 3 Раздел 1. Многопроцессорные системы. c тандарт MPI ………………………………………... 4 – 11 1.1. Организация и эффективность параллельных ЭВМ………. 4 – 6 1.2. Разновидности многопроцессорных систем……………...... 6 – 8 1.3. Техническая реализация параллельных ЭВМ……….……... 8 – 9 1.4. Системы программирования для многопроцессорных систем. Стандарт MPI………………………………………... 9 – 11 Раздел 2. Метод Монте-Карло ……………………………... 12 – 15 2.1. Применение метода статистического моделирования…….. 12 – 13 2.3. Способы получения случайных чисел……………………… 13 – 14 2.4. Идея и суть метода Монте-Карло........................................... 14 – 15 Раздел 3. Алгоритмы вычисления интегралов методом Монте-Карло ……………………........ 16 – 27 3.1. Определенный интеграл……………………………………... 16 – 18 3.2. Параллельные алгоритмы вычисления...…………………… 18 – 25 3.3. Эксперимент. Анализ эффективности некоторых из реализованных алгоритмов………………………………….. 25 – 27 Заключение ………………………………………………………… 28 Литература ………………………………………………………….. 29 ПриложениЯ ………………………………………………………... 30 – 41 ВВЕДЕНИЕ Классическое понимание алгоритма включает в себя в неявном виде утверждение о том, что процедура, которая может быть названа алгоритмом, имеет последовательный пошаговый характер. Поскольку человек, который решает задачу, и ЭВМ, построенная по классической схеме, выполняют программы вычислений последовательно, то и в информатике на первых порах рассматривались только последовательные методы вычислений. Однако развитие архитектуры ЭВМ привело к появлению структур, в которых вычислительный процесс может протекать по нескольким ветвям параллельно. Параллелизм возник естественным образом, как средство увеличения производительности вычислительных машин. Идея параллелизма была технически реализована в многопроцессорных системах , однородных средах, локальных вычислительных сетях и в ряде других аппаратных решений. Появление технических средств, способных организовать параллельное протекание вычислительного процесса, поставило две задачи. Одна из них – создание специальных программных средств для записи параллельных процессов в форме, которую может использовать операционная система вычислительной машины, системы или сети, чтобы фактически организовать параллельный процесс. Появление этой задачи породило новое направление в программировании – параллельное программирование . Другая задача была связана с созданием специальных методов решения задач, ориентированных на параллельную реализацию, поиском численных методов, которые допускали бы распараллеливание вычислительных процедур. Прежде всего, специалистов привлекли методы, опирающиеся на идею независимых статистических испытаний. Их еще называют методами Монте-Карло . Главное преимущество применения метода Монте-Карло в параллельных вычислениях состоит в том, что вычислительные процессы не зависят друг от друга, обеспечивая наибольший выигрыш во времени по сравнению с другими методами [5]. Цель данной работы – нахождение алгоритмов вычисления интегралов методом Монте-Карло, без обменов сообщениями между процессами на этапе накопления статистического материала. Раздел 1. Многопроцессорные системы. Стандарт MPI 1.1. Организация и эффективность параллельных ЭВМ Под параллельной ЭВМ понимают ЭВМ, состоящую из множества связанных определенным образом вычислительных блоков, которые способны функционировать совместно и одновременно выполнять множество арифметико-логических операций, принадлежащих одной задаче (программе) [8]. Классифицировать параллельные ЭВМ можно по различным критериям: виду соединения процессоров, способу функционирования процессорного поля, области применения и т.п. Однако наиболее употребительная классификация параллельных ЭВМ была предложена Флинном [8, 9] и отражает форму реализуемого ЭВМ параллелизма. Основными понятиями классификации являются «поток команд» и «поток данных». Под потоком команд упрощенно понимают последовательность команд одной программы. Поток данных – это последовательность данных, обрабатываемых одной программой. Согласно классификации Флинна имеется четыре больших класса ЭВМ: 1) ОКОД – О диночный поток К оманд – О диночный поток Д анных (SISD – S ingle I nstruction – S ingle D ata). Это последовательные ЭВМ, в которых выполняется единственная программа, т.е. имеется только один счетчик команд и одно арифметико-логическое устройство; 2) ОКМД – О диночный поток К оманд – М ножественный поток Д анных (SIMD – S ingle I nstruction – M ultiple D ata). В таких ЭВМ выполняется единственная программа, но каждая ее команда обрабатывает много чисел. Это соответствует векторной форме параллелизма; 3) МКОД – М ножественный поток К оманд – О диночный поток Д анных (MISD – M ultiple I nstruction – S ingle D ata). Подразумевается, что в данном классе несколько команд одновременно работает с одним элементом данных, однако эта позиция классификации Флинна на практике не нашла применения; 4) МКМД – М ножественный поток К оманд – М ножественный поток Д анных (MIMD – M ultiple I nstruction – M ultiple D ata). В таких ЭВМ одновременно и независимо друг от друга выполняется несколько программных ветвей, в определенные промежутки времени обменивающихся данными. Такие системы обычно называют многопроцессорными . Далее, в работе будут рассматриваться только многопроцессорные ЭВМ. Вопрос об эффективности параллельных ЭВМ возникает на разных стадиях исследования и разработки ЭВМ. Следует различать эффективность параллельных ЭВМ в процессе их функционирования и эффективность параллельных алгоритмов [8]. В зависимости от стадии разработки полезными оказываются различные характеристики эффективности ЭВМ. Одной из главных является ускорение r параллельной системы, которое определяется выражением: (1.1.1) где T 1 – время решения задачи на однопроцессорной системе, а Tn – время решения той же задачи на n -процессорной системе. Если обозначить W = W ск + W пр , где W – общее число операций в задаче, W пр – число операций, которые можно выполнять параллельно, а Wск – число скалярных (нераспараллеливаемых) операций; t – время выполнения одной операции, то закон Амдала будет иметь вид [8, 9]: (1.1.2) где a = W ск / W – удельный вес скалярных операций. Закон Амдала определяет принципиально важные для параллельных вычислений положения: 1. Ускорение вычислений зависит как от потенциального параллелизма задачи (величина 1-a ), так и от параметров аппаратуры (число процессоров n ). 2. Предельное ускорение вычислений определяется свойствами задачи и не может превосходить 1/a при любом числе процессоров, то есть максимальное ускорение определяется потенциальным параллелизмом задачи и весьма чувствительно к изменению величины a . Тем не менее закон Амдала не учитывает временных затрат на обмен сообщениями между процессорами. А ведь во многих случаях эти потери могут не только снизить ускорение вычислений, но и замедлить вычисления по сравнению с однопроцессорным вариантом . Поэтому, необходимо ввести величины W с – количество передач данных и tc – время одной передачи данных. Тогда модернизированный закон Амдала, называемый сетевым законом Амдала , будет выглядеть следующим образом: (1.1.3) где с = Wc ·tc / W ·t = cA ·cT – коэффициент сетевой деградации вычислений определяет объем вычислений, приходящийся на одну передачу данных (по затратам времени). При этом cA определяет алгоритмическую составляющую коэффициента деградации, обусловленную свойствами алгоритма, а cT – техническую составляющую, которая зависит от соотношения технического быстродействия процессора и аппаратуры сети. Из (1.1.3) видно, что для повышения скорости вычислений следует воздействовать на оба составляющие коэффициента деградации, которые зависят от множества факторов таких, как алгоритм задачи, размер данных, технические характеристики оборудования и пр. 1.2. Разновидности многопроцессорных систем В соответствии с классификацией Флинна к многопроцессорным ЭВМ относятся ЭВМ с множественным потоком команд и множественным потоком данных (МКМД-ЭВМ или MIMD-ЭВМ). В основе МКМД-ЭВМ лежит традиционная последовательная организация программы, расширенная добавлением специальных средств указания независимых, параллельно исполняемых фрагментов . Параллельно-последовательная программа привычна для пользователя и позволяет относительно просто собирать параллельную программу из обычных последовательных программ. МКМД-ЭВМ имеют две разновидности: ЭВМ с разделяемой (общей) и распределенной (индивидуальной) памятью. Структура этих ЭВМ представлена на рис. 1.2.1 [8]. Главное различие между МКМД-ЭВМ с общей и индивидуальной (локальной, распределенной) памятью состоит в характере адресной системы. a b Рис.1.2.1. Структура ЭВМ с разделяемой (а ) и индивидуальной (b) памятью. Здесь: П – процессор, ИП - индивидуальная память. В машинах с общей памятью адресное пространство всех процессоров является единым, следовательно, если в программах нескольких процессоров встречается одна и та же переменная VAL , то эти процессоры будут обращаться в одну физическую ячейку общей памяти. Это вызывает как положительные, так и отрицательные последствия [8, 9]: 1. Не нужно физически перемещать данные между взаимодействующими программами, которые параллельно выполняются в разных процессорах. Это исключает затраты времени на межпроцессорный обмен. 2. Поскольку одновременное обращение нескольких процессоров к общим данным может привести к получению неверных результатов, необходимо ввести систему синхронизации параллельных процессов (например, семафоры), что усложняет механизмы операционной системы. 3. Так как при выполнении каждой команды процессорам необходимо обращаться в общую память, то требования к пропускной способности коммутатора этой памяти чрезвычайно высоки, что ограничивает число процессоров в системах величиной 10...20. МКМД-ЭВМ с индивидуальной памятью получили большое распространение ввиду относительной простоты их архитектуры. В таких ЭВМ межпроцессорный обмен осуществляется с помощью передачи сообщений . Каждый процессор имеет независимое адресное пространство, и наличие одной и той же переменной VAL в программах разных процессоров приводит к обращению в физически разные ячейки индивидуальной памяти этих процессоров. Это вызывает физическое перемещение данных (обмен сообщениями) между взаимодействующими программами в разных процессорах, однако, поскольку основная часть обращений производится каждым процессором в собственную память, то требования к коммутатору ослабляются, и число процессоров в системах с распределенной памятью и скоростным коммутатором может достигать нескольких сотен и даже тысяч. 1.3. Техническая реализация параллельных ЭВМ Современные многопроцессорные системы делятся на три группы: симметричные мультипроцессоры (SMP); системы с массовым параллелизмом (МРР); кластеры [9]. Симметричные мультипроцессоры SMP (S ymmetric M ulti P rocessors) используют принцип разделяемой памяти. В этом случае система состоит из нескольких однородных процессоров и массива общей памяти (обычно из нескольких независимых блоков). Все процессоры имеют доступ к любой ячейке памяти с одинаковой скоростью. Процессоры подключены к памяти с помощью общей шины или коммутатора. Аппаратно поддерживается когерентность кэш-памяти. Наличие общей памяти сильно упрощает взаимодействие процессоров между собой, однако накладывает сильные ограничения на их число – не более 32 в реальных системах. Вся система работает под управлением единой ОС. Системы с массовым параллелизмом MPP (M assively P arallel P rocessing) содержат множество процессоров c индивидуальной памятью в каждом из них (прямой доступ к памяти других узлов невозможен), которые связаны через некоторую коммуникационную среду (высокоскоростная сеть, коммутатор и т.д.). Как правило, системы MPP благодаря специализированной высокоскоростной системе обмена, обеспечивают наивысшее быстродействие (и наивысшую стоимость). Кластерные системы - более дешевый вариант MPP-систем, поскольку они также используют принцип передачи сообщений, но строятся из компонентов высокой степени готовности. Базовым элементом кластера является локальная сеть. В качестве компьютеров могут выступать, например, рабочие станции. Оказалось, что на многих классах задач и при достаточном числе узлов такие системы дают производительность, сравнимую с суперкомпьютерной. Первым кластером на рабочих станциях был Beowulf [10]. Проект Beowulf начался летом 1994 года сборкой в научно-космическом центре NASA 16-процессорного кластера на Ethernet-кабеле. С тех пор кластеры на рабочих станциях обычно называют Beowulf-кластерами. Любой Beowulf-кластер состоит из отдельных машин (узлов) и объединяющей их сети (коммутатора). Кроме ОС, необходимо установить и настроить сетевые драйверы, компиляторы, ПО поддержки параллельного программирования и распределения вычислительной нагрузки. В качестве узлов обычно используются однопроцессорные ПЭВМ с быстродействием 1 ГГц и выше или SMP-серверы с небольшим числом процессоров (обычно 2-4). Для получения хорошей производительности межпроцессорных обменов, как правило, используют полнодуплексную сеть Fast Ethernet с пропускной способностью 100 Mбит/с. При этом для уменьшения числа коллизий устанавливают несколько "параллельных" сегментов Ethernet, или соединяют узлы кластера через коммутатор (switch). В качестве операционных систем обычно используют OC Linux или Windows NT и ее варианты. В качестве языков программирования используются С, С++ и старшие версии языка Fortran. Наиболее распространенным интерфейсом параллельного программирования в модели передачи сообщений является MPI (M essage P assing I nterface). 1.4. Системы программирования для многопроцессорных систем. Стандарт MPI Основными средствами программирования для многопроцеccорных систем являются две библиотеки, оформленные как стандарты: библиотека OpenMP для систем с общей памятью (для SMP-систем) и библиотека MPI для систем с индивидуальной памятью. Библиотека OpenMP является стандартом для программирования на масштабируемых SMP-системах. В стандарт входят описания набора директив компилятора, переменных среды и процедур. За счет идеи "инкрементального распараллеливания" OpenMP идеально подходит для разработчиков, желающих быстро распараллелить свои вычислительные программы с большими параллельными циклами. Разработчик не создает новую параллельную программу, а просто добавляет в текст последовательной программы OpenMP-директивы. Предполагается, что программа OpenMP на однопроцессорной платформе может быть использована в качестве последовательной программы, т.е. нет необходимости поддерживать последовательную и параллельную версии. Директивы OpenMP просто игнорируются последовательным компилятором, а для вызова процедур OpenMP могут быть поставлены заглушки, текст которых приведен в спецификациях. В OpenMP любой процесс состоит из нескольких нитей управления, которые имеют общее адресное пространство, но разные потоки команд и раздельные стеки. В простейшем случае, процесс состоит из одной нити. Стандарт MPI. Система программирования MPI предназначена для ЭВМ с индивидуальной памятью, то есть для многопроцессорных систем с обменом сообщениями. MPI имеет следующие особенности: 1. MPI - библиотека функций, а не язык. Она определяет состав, имена, вызовы функций и результаты их работы. Программы, которые пишутся на языках FORTRAN, C и C++, компилируются обычными компиляторами и связаны с MPI-библиотекой. 2. MPI - описание функций, а не реализация. Все поставщики параллельных компьютерных систем предлагают реализации MPI для своих машин бесплатно, реализации общего назначения также могут быть получены из Internet. Правильная MPI-программа должна выполняться на всех реализациях без изменения. В модели передачи сообщений процессы, выполняющиеся параллельно, имеют раздельные адресные пространства. Обмен происходит, когда часть адресного пространства одного процесса скопирована в адресное пространство другого процесса. Эта операция совместная и возможна только, когда первый процесс выполняет операцию передачи сообщения, а второй процесс - операцию его получения. Процессы в MPI принадлежат группам . Если группа содержит n процессов, то процессы нумеруются внутри группы номерами, которые являются целыми числами от 0 до n- l. Имеется начальная группа, которой принадлежат все процессы в реализации MPI. Контекст есть свойство коммуникаторов, которое позволяет разделять пространство обмена. Сообщение, посланное в одном контексте, не может быть получено в другом контексте. Контексты не являются явными объектами MPI, они проявляются только как часть реализации коммуникаторов. Понятия контекста и группы объединены в едином объекте, называемом коммуникатором . Таким образом, отправитель или получатель, определенные в операции посылки или получения, всегда обращается к номеру процесса в группе, идентифицированной данным коммуникатором. Раздел 2. Метод Монте-Карло 2.1. Применение метода статистического моделирования На практике существует множество операций[1] , построение аналитических моделей которых весьма затруднительно, а то и вообще невозможно, например, вследствие сложности рассматриваемой операции, содержащей множество случайных и неопределенных неконтролируемых факторов. В таких случаях применяется метод математического моделирования, известный под названием метода статистических испытаний или метода Монте-Карло [6]. Применение этого метода неразрывно связано с параллельными ЭВМ, позволяющими производить большие расчеты за сравнительно короткие промежутки времени. Метод статистических испытаний основан на построении модели операции для просчета на ЭВМ. Каждый раз в результате просчета на ЭВМ получаются результаты для конкретных исходных данных. Вследствие изменения исходных данных, контролируемых и неконтролируемых факторов накапливается большое количество статистического материала, который обрабатывается затем в соответствии с методами математической статистики. Полученный в результате накопления и обработки материал используется для выработки решений по оптимальному управлению операцией. Одним из недостатков метода Монте-Карло является его громоздкость и трудоемкость. Накапливаемое большое количество различной информации делает ее трудно обозримой. Поэтому, хотя метод Монте-Карло и является весьма мощным аппаратом в руках исследователя операции, его все же нельзя противопоставлять аналитическим методам. Прежде чем использовать метод Монте-Карло, полезно построить хотя бы грубую, приближенную аналитическую модель операции с целью выявления основных факторов, оказывающих влияние на исход операции. Преимущество метода статистических испытаний заключается в его универсальности, т.е. он может быть использован для исследования практически любых процессов или явлений, протекающих в различных областях человеческой деятельности и поддающихся количественному описанию. Метод статистического моделирования позволяет получать различные характеристики и зависимости между параметрами операции, что особенно важно, когда исследуемая операция является новой и находится в стадии проектирования. Метод статистического моделирования находит самое широкое применение в экономических, физических, математических исследованиях, а также в других областях естествознания. Так, в экономике метод Монте-Карло используется для изучения различных внутренних взаимодействий, протекающих в исследуемой операции; изучения свойств операции и выявления законов, оказывающих влияние на ее функционирование. Решение вопросов, связанных с прогнозированием экономики, оптимального планирования и управления рассматриваемой экономической операции [6]. В физике метод Монте-Карло применялся с целью нахождения интегралов от сложных уравнений при разработке первых ядерных бомб (интегралы квантовой механики). С помощью формирования больших выборок случайных чисел из, например, нескольких распределений, интегралы этих (сложных) распределений могут быть аппроксимированы из (сгенерированных) данных. Также метод Монте-Карло широко используется в задачах ядерной физики, задачах гидромеханики, газодинамики и пр. В математике метод статистических испытаний применяется для интегрирования, решения интегральных и дифференциальных уравнений, минимизации функций, обращения матриц, решения задач Дирихле и др. Кроме всего, метод Монте-Карло может применяться для решения логических задач [4]. 2.2. Способы получения случайных чисел При статистическом моделировании любых операций задача моделирования исходных данных – одна из важнейших. Как правило, потоки событий, моделирующие входную информацию, являются случайными, подчиняющимися некоторым известным законам распределения. Величина называется случайной , если она принимает те или иные значения в зависимости от появления или не появления некоторого случайного события [2, 4]. Имеется ряд таблиц [1, 4, 6], которые могут быть использованы для получения случайных чисел с заданным законом распределения. Но при работе на ЭВМ применение случайных чисел, записанных в таблицах, весьма неудобно, так как вначале эти числа необходимо ввести в ЭВМ, а потом хранить в памяти, занимая весьма значительную ее часть. Поэтому в настоящее время используются различные приемы, позволяющие получать (генерировать) случайные числа, необходимые по ходу реализации модели на ЭВМ, без хранения их предварительно в памяти компьютера. Применение физических генераторов случайных чисел основано на использовании случайных сигналов (шумов, т.е. случайным образом меняющегося напряжения u ( t ) ), получаемых от некоторых электронных ламп. Но и этот способ, как и табличный не лишен недостатков. Один из них состоит в том, что расчеты на ЭВМ в силу различных причин приходится выполнять несколько раз. А с помощью этого способа генерировать каждый раз одни и те же случайные числа невозможно, если их по ходу расчетов на ЭВМ не запоминать. Поэтому при реализации модели на ЭВМ используются специальные аналитические зависимости, позволяющие получать случайные числа, обладающие заданными свойствами. Такие числа называются псевдослучайными . 2.3. Идея метода Монте-Карло Метод Монте-Карло – метод, в котором используется искусственная реализация вероятностных законов. При применении метода статистических испытаний каждый раз проводится одна реализация операции (один просчет модели на ЭВМ) при случайных исходных данных . Многократный просчет операции по ее модели позволяет накопить данные, на основании которых уже можно судить об исследуемой операции. Сущность метода Монте-Карло состоит в следующем: требуется найти значение а некоторой изучаемой величины. Для этого выбирают такую случайную величину Х, математическое ожидание которой равно а: М(Х) = а. Практически же поступают так: производят n испытаний, в результате которых получают n возможных значений Х ; вычисляют их среднее арифметическое: (2.3.1) и принимают в качестве оценки (приближённого значения) a * искомого числа a : (2.3.2) Теория этого метода указывает, как наиболее целесообразно выбрать случайную величину Х , как найти её возможные значения. В частности, разрабатываются способы уменьшения дисперсии используемых случайных величин, в результате чего уменьшается ошибка, допускаемая при замене искомого математического ожидания а его оценкой а * . Раздел 3. Алгоритмы вычисления интегралов методом Монте-Карло 3.1. Определенный интеграл Пусть функция f ( x ) определена на конечном промежутке [a , b ]. Разобьем отрезок [a , b ] каким-либо образом на N частей (не обязательно равных) точками деления [7]: a = x 0 < x 1 < x 2 < x 3 < … < xN -1 < xN = b . (3.1.1) Число назовем мелкостью разбиения [ a , b ] . Выберем произвольные точки ξi из интервала [xi , xi +1 ] и составим интегральную сумму (3.1.2) функции f ( x ) , отвечающую разбиению (3.1.1) и выбору точек ξi . Если последовательность интегральных сумм IN при δN → 0 имеет конечный предел I , не зависящий ни от способа разбиения отрезка [a , b ], ни от выбора точек ξi , то этот предел называют определенным интегралом функции f ( x ) в промежутке от a до b : (3.1.3) Пусть, теперь, требуется приближенно вычислить интеграл (3.1.3). Предположим, что каким-то образом удается получить N случайных, попарно независимых точек ξ0 , …, ξN -1 , одинаково распределенных на интервале L = [a , b ] с функцией распределения F ( ξ) . Введем обозначение: (3.1.4) Случайные величины si попарно независимы и одинаково распределены: (3.1.5) где M – оператор математического ожидания величины si [2]. Дисперсия D случайной величины si [2]: (3.1.6) Определим некоторую величину SN : (3.1.7) Из указанных выше свойств величин si следует, что [1]: (3.1.8) (3.1.9) Таким образом, получив N случайных чисел ξi и вычислив в них значения функции f , можем определить величину IN , которая при N больших даст, приближенно, интеграл I . Для случайных чисел с равномерным законом распределения получаем: (3.1.10) Формулы (3.1.7) и (3.1.10) определяют первый способ вычисления интегралов методом статистических испытаний или методом Монте-Карло. Назовем его методом средних значений. Задав уровень надежности η , с которым желательно получить приближенное значение интеграла, можем применить неравенство Чебышева [1, 2] для определения условия прекращения итерационного процесса, причем это неравенство выполняется с вероятностью 1 – η: (3.1.11) Для законности применения неравенства Чебышева нужно выполнение предположения о независимости распределения любых точек ξi . Так как нам не известна функция распределения F ( ξ) , на практике будем использовать следующее условие прекращения итерационного процесса (дальше – условие выхода ): если модуль разности значений оценок интеграла (3.1.3) на предыдущем и текущем шагах меньше либо равен заданной величины ошибки ε, прекратить выполнение процедуры : (3.1.12) Кроме вариации метода Монте-Карло, описанного в предыдущем пункте, для вычисления однократного интеграла можно определить еще один способ – метод площадей . Поясним метод площадей на примере. Пусть дана некоторая, в общем случае знакопеременная, функция f ( x ) , пределы интегрирования a и b (рис. 3.1.1). Интеграл (3.1.1) будет представлять собой, в этом случае, общую площадь заштрихованных на рисунке 3.1.1 частей. Обозначим площадь прямоугольника, в который вписана кривая f ( x ) , – S , искомую площадь – σ. Каким-либо образом получим N случайных точек x и N случайных точек y , лежащих внутри прямоугольника. Рис. 3.1.1 Если в заштрихованную область попало n пар точек, то искомая площадь будет выражаться формулой: (3.1.13) 3.2. Параллельные алгоритмы вычисления В случае, последовательной ЭВМ для применение метода Монте-Карло необходимо наличие генератора случайных числе, вычислителя, анализатора. Все эти элементы действуют линейно, по цепочке: генератор → вычислитель → анализатор: Если же мы обладаем параллельной системой, например, – кластером, то промежуточные расчеты в методе Монте-Карло могут осуществляться на разных узлах этого кластера, а окончательный результат обрабатываться на некоторой главной машине, так называемой корневой. Такая схема изображена на рис. 3.2.1: Согласно рис. 3.2.1, генератор случайных чисел один для системы и, проходя сверху вниз по вычислителям, выдает каждому из них сгенерированное случайное число. Однако, из-за того, что информация должна быть передана между компьютерами, снижается производительность параллельной системы. Передаваться информация будет средствами, доступными в настоящее время. Скорость же передачи определяется для технологии Fast Ethernet до 12 Мбайт/с, более современные технологии SCI или Myrinet от 80 Мбайт/с [9]. Кроме значительной разницы в пропускной способности, последние высокоскоростные технологии имеют значительно меньшую латентность (5-10 мкс против 150-300 мкс Fast Ethernet). Основной недостаток этих технологий – стоимость, превышающая, порой, стоимость вычислительных элементов. Рис. 3.2.1. Упрощенная схема параллельного алгоритма метода Монте-Карло Чтобы свести обмен сообщениями к минимуму можно каждому вычислителю предоставить свой генератор случайных (псевдослучайных) чисел. Кроме того, генератором и вычислителем может быть один и тот же процесс. Кроме того, в одном, главном процессе могут быть объединены генератор, вычислитель и анализатор. Именно такая схема и была реализована в работе. Графически она изображена на рис. 3.2.2: В работе были разработаны и реализованы три параллельных алгоритма вычисления однократных интегралов методом Монте-Карло на кластере MPI, основанные на модели рис. 3.2.2. Все три алгоритма базируются на использовании трех функций MPI [9]: 1. MPI _ BARRIER ( comm ) : входной параметр comm (тип – дескриптор) является коммуникатором группы, в которой выполняется коллективная операция. Эта функция – функция барьерной синхронизации, блокирует вызывающий ее процесс, пока все процессы группы не вызовут функцию. В каждом процессе управление возвращается только тогда, когда все процессы в группе вызовут MPI _ BARRIER . Синхронизацией называется установление правильной очередности процессов. Необходимость в синхронизации вызывается либо распределением общего ресурса между процессами, либо логической зависимостью процессов друг от друга. 2. MPI _ REDUCE ( sendbuf , recvbuf , count , datatype , op , root , comm ) : входной параметр sendbuf (тип – альтернатива) содержит адрес буфера посылки; recvbuf (тип – альтернатива) содержит адрес принимающего буфера (используется только корневым процессом); count (тип – целое) количество элементов в посылающем буфере; datatype (тип – дескриптор) тип данных буфера посылки; op (тип – дескриптор) содержит операцию редукции; root (тип – целое) номер главного процесса; comm (тип – дескриптор) коммуникатор группы, в которой выполняется коллективная операция. Функция MPI _ REDUCE объединяет элементы входного буфера каждого процесса, используя операцию op , и возвращает объединенное значение в выходной буфер процесса с номером root . 3. MPI _ BCAST ( buffer , count , datatype , root , comm ) : двусторонний (входной/выходной) параметр buffer (тип – альтернатива) содержит адрес начала буфера посылки/приема; входной параметр count (тип – целое) содержит количество записей в буфере; входной параметр datatype (тип – дескриптор) описывает тип данных в буфере; входной параметр root (тип – целое) содержит номер корневого процесса; входной параметр comm (тип – дескриптор) является коммуникатором группы, в которой выполняется коллективная операция. Функция широковещательной передачи MPI _ BCAST посылает сообщение из корневого процесса всем процессам группы, включая себя. Она вызывается всеми процессами группы с одинаковыми аргументами для comm , root . В момент возврата управления содержимое корневого буфера обмена будет уже скопировано во все процессы. Интерфейсы функций (Integral 1 DMK _ par (…) , Integral 1 DMK _ par 1(…) , Integral 1 DMK _ par 2(…) ), реализующих все три алгоритма приведены в приложении 3, а их реализации – в приложении 4. Рассмотрим все три алгоритма. 1. Алгоритм, в основе которого лежит скалярный параллелизм, т.е. распараллеливание цикла , реализован в функции Integral 1 DMK _ par (…) . В этом алгоритме генерация чисел осуществляется с некоторым шагом равным отношению интервала интегрирования к количеству генерируемых точек. Это позволяет, по мере уменьшения шага, добиваться равномерного распределения чисел по оси координат. Как следствие – увеличение сходимости метода, а также повышение точности вычислений. В этом алгоритме каждый процесс выполняет следующие основные шаги: - Получение количества запущенных процессов, получение процессом своего номера. - Начало внешнего цикла. - Инициализация переменных, содержащих значение интеграла на двух соседних итерациях (для проверки условия выхода). - Инициализация генератора псевдослучайных чисел номером процесса и текущим временем. - Вычисление шага, с которым ведется генерация чисел. - Параллельный внутренний цикл, в котором вычисляется случайная точка, значение функции в этой точке и накапливается частичная сумма. В цикле каждый процесс совершает, в среднем , totpt / totproc итераций (totproc – количество запущенных процессов, totpt – количество генерируемых точек). В таблице 3.2.1 показан пример распределения итераций цикла по процессам для 6-ти процессорной группы (totproc = 6) и totpt = 20: Таблица 3.2.1. Процессы 0 1 2 3 4 5 Итерации 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 - - - - - Синхронизация всех процессов. - Объединение частичных сумм всех процессов. - Вычисление корневым процессом значения интеграла. - Синхронизация процессов. - Широковещание всем процессам из корневого процесса значения интеграла. - Увеличение числа генерируемых точек для следующего прохода внешнего цикла. - Проверка условия выхода. Если ложь – прервать внешний цикл, если истина – продолжить внешний цикл с увеличенных количеством генерируемых точек. - Возвращение функцией полученного значения интеграла. 2. Алгоритм, в основе которого лежит идея сужения интервала интегрирования для каждого процесса (процессорный интервал ), реализован в функции Integral 1 DMK _ par 1(…) . В этом алгоритме, как и в предыдущем, генерация чисел осуществляется с некоторым шагом равным отношению интервала интегрирования L к количеству генерируемых точек. Здесь каждый процесс вычисляет интеграл на своем участке (частичный интеграл ), а затем эти значения объединяются операцией редукции. Это позволяет, как и в предыдущем алгоритме, по мере уменьшения шага, добиваться равномерного распределения чисел по оси координат, а также точнее изучить поведение функции на узком интервале. Как следствие – увеличение сходимости метода, а также повышение точности вычислений. В этом алгоритме каждый процесс выполняет следующие основные шаги: - Получение количества запущенных процессов, получение процессом своего номера. - Вычисление переменных, определяющих рабочие интервалы каждого процесса. - Начало внешнего цикла. - Инициализация переменных, содержащих значение интеграла на двух соседних итерациях (для проверки условия выхода). - Инициализация генератора псевдослучайных чисел номером процесса и текущим временем. - Вычисление шага, с которым ведется генерация чисел. - Внутренний цикл, в котором вычисляется случайная точка, значение функции в этой точке и накапливается частичная сумма. - Вычисление процессом частичного интеграла на узком промежутке. - Увеличение числа генерируемых точек для следующего прохода внешнего цикла. - Синхронизация всех процессов. - Объединение частичных интегралов всех процессов. - Синхронизация процессов. - Широковещание из корневого процесса интеграла на L всем процессам. - Проверка условия выхода. Если ложь – прервать внешний цикл, если истина – продолжить внешний цикл с увеличенных количеством генерируемых точек. - Возвращение функцией полученного значения интеграла. 3. Алгоритм, в основе которого лежит идея сужения интервала интегрирования для каждого процесса и вычисление полного интеграла на этом участке , реализован в функции Integral 1 DMK _ par 2(…) . В этом алгоритме, как и в предыдущих, генерация чисел осуществляется с некоторым шагом равным отношению интервала интегрирования L к количеству генерируемых точек. Здесь каждый процесс вычисляет интеграл на своем участке до тех пор, пока не получит значения с заданной точностью, а затем эти значения объединяются операцией редукции после завершения внешнего цикла. Это позволяет, как и в предыдущих алгоритмах по мере уменьшения шага, добиваться равномерного распределения чисел по оси координат, а также точнее изучить поведение функции на узком интервале. Кроме того, обмен сообщениями осуществляется только один раз. Как следствие – увеличение сходимости метода, повышение точности вычислений и уменьшение времени вычислений на некоторых классах функций. В этом алгоритме каждый процесс выполняет следующие основные шаги: - Получение количества запущенных процессов, получение процессом своего номера. - Вычисление переменных, определяющих рабочие интервалы каждого процесса. - Начало внешнего цикла. - Инициализация переменных, содержащих значение интеграла на двух соседних итерациях (для проверки условия выхода). - Инициализация генератора псевдослучайных чисел номером процесса и текущим временем. - Вычисление шага, с которым ведется генерация чисел. - Внутренний цикл, в котором вычисляется случайная точка, значение функции в этой точке и накапливается частичная сумма. - Вычисление процессом частичного интеграла на узком промежутке. - Увеличение числа генерируемых точек для следующего прохода внешнего цикла. - Проверка условия выхода. Если ложь – прервать внешний цикл, если истина – продолжить внешний цикл с увеличенных количеством генерируемых точек. - Синхронизация процессов. - Редукция полных интегралов на всех процессорных интервалах. - Возвращение функцией полученного значения интеграла. Временные диаграммы синхронных вычислительных процессов с независимыми ветвями приведены на рис. 3.2.3. Еще раз подчеркнем, что для увеличения сходимости метода существенную роль играет, состав последовательности псевдослучайных чисел или, иными словами, ее качество. Рис. 3.2.3 3.3. Эксперимент. Анализ эффективности некоторых из реализованных алгоритмов Были проведены исследования эффективности разработанных алгоритмов на кластере MPI, с компонентами: процессор – Pentium 1 Ггц, ОЗУ – 128 Мб, сеть – Fast Ethernet. В экспериментах тестировалась подынтегральная функция . Пределы интегрирования составляли a = 0, b = 1e6. По результатам тестирования были получены кривые 3.3.1 и 3.3.2. В таблице 3.3.1 приведены результаты тестирования алгоритмов на различных функциях. Сравнивались результаты вычисления на 1 и 5 процессорных системах Таблица 3.3.1 Функция a b t, с (1 проц.) t, с (5 проц.) 2 4 0,001784 0,007947 sinx/x 0,0000001 1 0,000555 0,002694 x 3 -100 100 75,140435 7,556625 0 3 39,392135 15,766564 0 100000 6,454987 2,638202 100 -0,0010000 1000 0,004833 0,00334 0 1000000 131,795787 20,773000 По результатам экспериментов были сделаны следующие выводы: 1. Алгоритм 1 работает наиболее эффективно со сложными функциями в широком диапазоне L . В то время как остальные заметно уступают ему в работе с такими функциями. 2. Для небольших функций в небольших пределах L 3-ий алгоритм работают быстрее других, т.к. содержит только одну операцию обмена данными. Для того же класса функций 1-ый алгоритм проигрывает, т.к. время, затраченное на обмен данными оказывается больше времени вычислений. Таким образом, для вычислений интегралов могут быть использованы все три алгоритма – в зависимости от класса функции и пределов интегрирования. ЗАКЛЮЧЕНИЕ В заключении хотим отметить следующее. 1. Полученные алгоритмы легко обобщаются на случай многомерных пространств. Это осуществляется путем генерации в функции чисел, отвечающих за нужную координату. Так, в приложениях 3, 4 приведены тексты функций, реализующих вычисление двойных интегралов на последовательной и параллельной ЭВМ. 2. Для вычисления n -кратного интеграла необходимо знать величину n -1-кратного интеграла. Поэтому при расчете интегралов высокой кратности необходимо рассчитывать распределение вычислений по процессам особенно тщательно. Для решения этой проблемы в MPI предусмотрены виртуальная топология, группы, коммуникаторы. 3. Вследствие низкой сходимости метода Монте-Карло целесообразно применять его только при расчетах интегралов кратности 2 и выше. ЛИТЕРАТУРА 1. Бахвалов Н.С. Численные методы. – М.: «Наука», 1975. – 631 с. 2. Грэхем Р., Кнут Д., Паташник О. Конкретная математика. – М.: «Мир», 1998. – 703 с. 3. Демидович Б.П. Сборник задач и упражнений по математическому анализу. – М.: «Наука», 1972. – 544 с. 4. Демидович Б.П., Марон И.А. Основы вычислительной математики. – М.: «Наука», 1966. – 664 с. 5. Информатика. Энциклопедический словарь для начинающих. – М.: «Педагогика-пресс», 1994. – 349 с. 6. Костевич Л.С., Лапко А.А. Теория игр. Исследование операций. – Мн.: «Вышэйшая школа», 1982. – 210 с. 7. Куст Ю., Юмагулов М. Математика. Основы математического анализа. – М.: «Айрис Пресс», 1999. – 270 с. 8. Шпаковский Г.И. Организация параллельных ЭВМ и суперскалярных процессоров. – Мн.: БГУ, 1996. – 283 с. 9. Шпаковский Г.И., Серикова Н.В. Программирование для многопроцессорных систем в стандарте MPI. – Мн.: БГУ, 2002. – 324 с. 10. Что такое Beowulf? (http://acdwp.narod.ru/la/podr/05020010.htm). ПРИЛОЖЕНИЕ 1 // Файл DataType.h // Содержит определенные и используемые в работе типы данных // 20.05.04 г. // Целые знаковые и беззнаковые typedef unsigned int CODE_ERR; // Тип для возврата результата-ошибки typedef unsigned int UINT; // Тип для счетчиков typedef unsigned long int ULONG; // Тип для счетчиков // Вещественные знаковые и беззнаковые typedef float PTCOUNT; // Хранит кол-во точек typedef long double LIMIT; // Тип для границ диапазона typedef long double FUNC_TYPE; /* Тип результата подынтегральной функции */ typedef long double ARG_TYPE; // Тип для аргументов функций typedef unsigned long double UFUNC_TYPE; // Тип для площади typedef unsigned long double ERR_T; // Тип для точности // Указатели на функции typedef FUNC_TYPE (*SURF) (ARG_TYPE, ARG_TYPE); /* Функция двух аргументов */ typedef FUNC_TYPE (*CURVLINE) (ARG_TYPE); // Функция одного аргумента. ПРИЛОЖЕНИЕ 2 // __MyHeader__.h // Инлайн реализации вспомогательных функций // Последние изменения 20 мая 2004 г. /////////////////////////////////////////////////////////////////////////////////////////////////////////// // Включения # include <time.h> #include <stdio.h> #include <stdlib.h> #include "DataType.h" /////////////////////////////////////////////////////////////////////////////////////////////// inline CODE_ERR TotTimePrint (clock_t start, clock_t stop, FILE *fp, char mode = 'f') /* Инлайн функция печатает в файл (по умолчанию) или на дисплей (mode == 'd ') время в минутах и секундах между отметками start и stop , полученными с помощью функции clock() . При этом, если время меньше 1 секунды функция выводит, что общее время выполнения чего-либо – < 1 sec . Возвращает код ошибки: "0" – если ошибок нет, "−1" – если значения start и/или stop – ошибочные или указан несуществующий режим */ { if (start < 0 || stop < 0) return −1; // Подсчитываем время выполнения clock_t t = (stop-start) / CLK_TCK; switch (mode) { case 'f': // Печать в файл if (t < 1) fprintf (fp, "Total time < 1 sec\n"); else fprintf (fp, "Total time = %d sec (%d min %d sec)\n", t, t/60, (t - (t/60)*60)); break ; case 'd': // Печать на экран if (t < 1) printf ("Total time < 1 sec\n"); else printf ("Total time = %d sec (%d min %d sec)\n", t, t/60, (t - (t/60)*60)); break ; default : return −1; break ; } return 0; } inline long double random (LIMIT left, LIMIT right) /*Инлайн функция возвращает случайное число в интервале от left до right */ { return ((long double )rand()) * (right-left) / RAND_MAX + left; } ПРИЛОЖЕНИЕ 3 // Файл IntegralMK.h /* Описания функций для подсчета интегралов методом Монте-Карло */ // 4 мая 2004 г. – 25 мая 2004 г. /////////////////////////////////////////////////////////////////////////////////////////////// // Включения # include <mpi.h> #include "DataType.h" #include "__MyHeader__.h" /////////////////////////////////////////////////////////////////////////////////////////////// // Константы # define EPS_ERR 1 // Неккоректная величина абсолютной ошибки # define TOTPT_MAX 1e30 // Максимально допустимое число точек /////////////////////////////////////////////////////////////////////////////////////////////// // Поверхность FUNC_TYPE Surf (ARG_TYPE x, ARG_TYPE y); // Границы у -правильной области, т.е. bot ( x ) и top ( x ) FUNC_TYPE bot (ARG_TYPE arg); // Нижняя граница FUNC_TYPE top (ARG_TYPE arg); // Верхняя граница // Подынтегральная функция для 1-кратного интеграла FUNC_TYPE f (ARG_TYPE arg); /////////////////////////////////////////////////////////////////////////////////////////////// // Интерфейсы функций FUNC_TYPE Integral1DSF (CURVLINE f, LIMIT left, LIMIT right, ERR_T eps = 1e-4); // Рассчет 1-кратного интеграла по формуле парабол (Симпсона) /* Использовалась для сравнения со значениями, полученными методом Монте-Карло */ // П Р О Т Е С Т И Р О В А Н А ! ! ! – 22.05.04 г. FUNC_TYPE Integral1DMK (CURVLINE f, LIMIT left, LIMIT right, ERR_T eps = 1e-3, PTCOUNT totpt = 100, ULONG mult = 2); // Вычисление 1-кратного интеграла методом Монте-Карло // Вариант со средним значением функции // П Р О Т Е С Т И Р О В А Н А ! ! ! – 25.05.04 г. FUNC_TYPE Integral1DMK_par (CURVLINE f, LIMIT left, LIMIT right, ERR_T eps = 1e-3, PTCOUNT totpt = 100, ULONG mult = 2); // Вычисление 1-кратного интеграла методом Монте-Карло на кластере MPI // Вариант со средним значением функции // П Р О Т Е С Т И Р О В А Н А // НА КЛАСТЕРЕ ! ! ! – 26.05.04 г. FUNC_TYPE Integral1DMK_par1 (CURVLINE f, LIMIT left, LIMIT right, ERR_T eps = 1e-3, PTCOUNT totpt = 100, ULONG mult = 2); /* Вычисление 1-кратного интеграла методом Монте-Карло на кластере MPI. Вариант со средним значением функции. Каждый процесс работает в более узком диапазоне, чем первоначальный */ // П Р О Т Е С Т И Р О В А Н А // НА КЛАСТЕРЕ ! ! ! – 26.05.04 г. FUNC_TYPE Integral1DMK_par2 (CURVLINE f, LIMIT left, LIMIT right, ERR_T eps = 1e-3, PTCOUNT totpt = 100, ULONG mult = 2); /* Вычисление 1-кратного интеграла методом Монте-Карло на кластере MPI. Вариант со средним значением функции. Каждый процесс работает в более узком диапазоне, чем начальный. Операция редукции проводится однажды, когда каждый процесс подсчитал интеграл на своем узком участке */ // П Р О Т Е С Т И Р О В А Н А // НА КЛАСТЕРЕ ! ! ! – 26.05.04 г. UFUNC_TYPE Square (CURVLINE f, LIMIT left, LIMIT right, ERR_T eps = 1e-4); /* Вычисление площади под кривой f ( x ) . Используется формула Симпсона */ // П Р О Т Е С Т И Р О В А Н А ! ! ! – 24.05.04 г. inline UFUNC_TYPE RegionY (CURVLINE bottom, CURVLINE top, LIMIT left, LIMIT right, unsigned char var, ERR_T eps = 1e-4) /* Вычисление площади между кривыми bottom ( x ) и top ( x ) . Кривые сшиваются в точках left и right */ // П Р О Т Е С Т И Р О В А Н А ! ! ! – 25.05.04 г. { // var определяет случай вычисления площади области switch (var) { // Кривые top и bottom лежат в положительной полуоси case 1: return Square(top, left, right, eps) – Square(bottom, left, right, eps); break ; // Кривые top и bottom лежат в отрицательной полуоси case 2: return Square(bottom, left, right, eps) – Square(top, left, right, eps); break ; // Кривая top лежит в положительной полуоси // кривая bottom лежит в отрицательной полуоси case 3: return Square(top, left, right, eps) + Square(bottom, left, right, eps); break ; default : return 0; } } FUNC_TYPE Integral2DMK (SURF surf, LIMIT left, LIMIT right, CURVLINE bottom, CURVLINE top, unsigned char var, ERR_T eps = 1e-3, PTCOUNT totpt = 100, ULONG mult = 2); /* Рассчет двойного интеграла для у -правильной области. Вариант со средним значением функции */ FUNC_TYPE Integral2DMK_par (SURF surf, LIMIT left, LIMIT right, CURVLINE bottom, CURVLINE top, unsigned char var, ERR_T eps = 1e-3, PTCOUNT totpt = 100, ULONG mult = 2); /* Рассчет двойного интеграла для у-правильной области на кластере MPI. Вариант со средним значением функции */ ПРИЛОЖЕНИЕ 4 // Файл IntegralMK.h // Реализация метода Монте-Карло для подсчета интегралов // Начало – 4 мая 2004 г. /////////////////////////////////////////////////////////////////////////////////////////////// #include <math.h> #include <limits.h> #include "IntegralMK.h" /////////////////////////////////////////////////////////////////////////////////////////////// FUNC_TYPE Integral1DSF (CURVLINE f, LIMIT left, LIMIT right, ERR_T eps) // Реализация формулы Симпсона для однократных интегралов /* Использовалась для сравнения со значениями, полученными методом Монте-Карло */ // П Р О Т Е С Т И Р О В А Н А ! ! ! – 22.05.04 г. { if (eps <= 0) exit (EPS_ERR); FUNC_TYPE _s, s = 0.0; FUNC_TYPE t0 = f(left) + f(right); ULONG n = 1024; do { _s = s; ARG_TYPE h = (right - left) / (2*n); FUNC_TYPE t1 = 0, t2 = 0; // t1 - сумма нечетных, t2 – сумма четных for (ULONG j = 1; j <= 2*n-1; j++) (j % 2 == 0) ? t2 += f(left+j*h) : t1 += f(left+j*h); s = (t0 + 4*t1 + 2*t2) *(h/3); n *= 2; // mult == 2 } while (fabsl(s – _s) > eps && n <= ULONG_MAX); return s; } /////////////////////////////////////////////////////////////////////////////////////////////// FUNC_TYPE Integral1DMK (CURVLINE f, LIMIT left, LIMIT right, ERR_T eps, PTCOUNT totpt, ULONG mult) // Вычисление однократного интеграла методом Монте-Карло. // Вариант со средним значением функции // П Р О Т Е С Т И Р О В А Н А ! ! ! – 25.05.04 г. { // Проверка корректности введенных данных if (eps <= 0 || mult < 1 || totpt < 10 || totpt >= TOTPT_MAX) exit(1); FUNC_TYPE _s = 0.0, s; ARG_TYPE h, left_lim; clock_t start = clock(), stop; // Запускаем таймер do { s = _s; // Запоминаем предыдущее значение интеграла _s = 0.0; // Обнуляем сумму // Инициализируем генератор случайных чисел временем srand((UINT)(time(NULL))); do { // Обеспечиваем малый шаг h = (right - left) / totpt; if (h > 1) totpt *= mult; } while (h > 1); if (!h) exit(1); left_lim = left; // Накапливаем сумму while ((h > 0 && left_lim + h <= right) || (h < 0 && left_lim + h >= right)) { _s += f(random(left_lim, left_lim+h)); left_lim += h; } // Вычисляем интеграл по формуле (3.1.10) _s = (_s/totpt)*(right-left); // Вывод на экран промежуточных данных printf ("\nValue = %G (total point = %G)\n", _s, totpt); // Интеграл и кол-во точек printf ("Current error = %G\n", fabsl(s - _s)); // Абсолютная ошибка totpt *= mult; } while (fabsl(s - _s) > eps && totpt <= TOTPT_MAX); stop = clock(); // Останавливаем таймер TotTimePrint(start, stop, NULL, 'd'); // Печатаем время выполнения цикла return _s; } /////////////////////////////////////////////////////////////////////////////////////////////// FUNC_TYPE Integral1DMK_par (CURVLINE f, LIMIT left, LIMIT right, ERR_T eps, PTCOUNT totpt, ULONG mult) /* Вычисление однократного интеграла методом Монте-Карло на кластере MPI. Вариант со средним значением функции */ // П Р О Т Е С Т И Р О В А Н А // НА КЛАСТЕРЕ ! ! ! – 26.05.04 г. { // Проверка корректности введенных значений if (eps <= 0 || mult < 1 || totpt < 10 || totpt >= TOTPT_MAX) exit(1); int totproc, self, root = 0; MPI_Comm comw = MPI_COMM_WORLD; MPI_Comm_size(comw, &totproc); // Получение кол-ва процессов MPI_Comm_rank(comw, &self); // Получение процессом своего номера FUNC_TYPE _s = 0.0, _s_root = 0.0, s; ARG_TYPE h; //srand(self); //MPI_Barrier(comw); //double mpi_start = MPI_Wtime(); do { s = _s; _s = 0.0; // Инициализация генератора случайных чисел номером процесса self srand((UINT)(time(NULL) >> self)); do { // Обеспечиваем малый шаг h = (right - left) / totpt; if (h > 1) totpt *= mult; } while (h > 1); if (!h) exit(1); /* Распараллеливание цикла накапливания суммы _s значений подынтегральной функции f в случайных точках, лежащих в пределах [left, left + i*h] */ ULONG i = self + 1; while ((h > 0 && left + i*h <= right) || (h < 0 && left + i*h >= right)) { // Накапливаем сумму _s += f(random(left + (i-1)*h, left + i*h)); i += totproc; } // Синхронизация всех процессов, т.к процессы выполняют // разное кол-во итераций в цикле MPI_Barrier(comw); // Суммирование _s во всех процессах и сохранение полученного // значения в корневом процессе в переменной _s_root MPI_Reduce(&_s, &_s_root, 1, MPI_LONG_DOUBLE, MPI_SUM, root, comw); if (!self) { // В корневом процессе: // Рассчет интеграла по формуле (3.1.10) _s = (_s_root/totpt)*(right-left); // Вывод на экран корневого процесса промежуточных данных printf("\nValue = %G (total point = %G)\n", _s, totpt); printf ("Current error = %G\n", fabsl(s - _s)); } // Синхронизация всех процессов MPI_Barrier(comw); /* Широковещание полученного значения интеграла из корнево го процесса всем остальным, включая самого себя */ MPI_Bcast(&_s, 1, MPI_LONG_DOUBLE, root, comw); /* Если ошибка вычислений больше допустимой, увеличение числа разыгрываемых точек в mult раз */ totpt *= mult; } while (fabsl(s – _s) > eps && totpt <= TOTPT_MAX); //double mpi_stop = MPI_Wtime(); //printf("Total time = %f\n", mpi_stop-mpi_start); return _s; } /////////////////////////////////////////////////////////////////////////////////////////////// FUNC_TYPE Integral1DMK_par1 (CURVLINE f, LIMIT left, LIMIT right, ERR_T eps, PTCOUNT totpt, ULONG mult) /* Вычисление однократного интеграла методом Монте-Карло на кластере MPI Вариант со средним значением функции. Каждый процесс считает интеграл на уменьшенном интервале*/ // П Р О Т Е С Т И Р О В А Н А // НА КЛАСТЕРЕ ! ! ! – 26.05.04 г. { // Проверка корректности введенных значений if (eps <= 0 || mult < 1 || totpt < 10 || totpt >= TOTPT_MAX) exit(1); int totproc, self, root = 0; MPI_Comm comw = MPI_COMM_WORLD; MPI_Comm_size(comw, &totproc); // Получение кол-ва процессов MPI_Comm_rank(comw, &self); // Получение процессом своего номера FUNC_TYPE _s = 0.0, _s_root = 0.0, s; ARG_TYPE h, H, left_lim; // Для каждого из процессов определяем свой интервал интегрирования H = (ARG_TYPE)(right - left) / totproc; left += self*H; right = left + H; printf("%f\t%f\n", left, right); //srand(self); //MPI_Barrier(comw); //double mpi_start = MPI_Wtime(); do { s = _s; // Запоминаем предыдущее значение интеграла на промежутке [left, right] _s = 0.0; // Обнуляем сумму // Инициализируем генератор случайных чисел временем srand((UINT)(time(NULL) >> self)); do { // Обеспечиваем малый шаг h = (right - left) / totpt; if (h > 1) totpt *= mult; } while (h > 1); if (!h) exit(1); left_lim = left; // Накапливаем сумму while ((h > 0 && left_lim + h <= right) || (h < 0 && left_lim + h >= right)) { _s += f(random(left_lim, left_lim+h)); left_lim += h; } // Вычисляем интеграл "через" среднее арифметическое функции _s = (_s/totpt)*(right-left); totpt *= mult; MPI_Barrier(comw); // Cинхронизируем процессы // Суммируем значения интеграла, полученные на малых промежутках MPI_Reduce(&_s, &_s_root, 1, MPI_LONG_DOUBLE, MPI_SUM, root, comw); if (!self) { _s = _s_root; // Вывод на экран промежуточных данных printf("\nValue = %G (total point = %G)\n", _s, totpt); // Интеграл и кол-во точек printf ("Current error = %G\n", fabsl(s - _s)); // Абсолютная ошибка } MPI_Barrier(comw); // Cинхронизируем процессы MPI_Bcast(&_s, 1, MPI_LONG_DOUBLE, root, comw); } while (fabsl(s – _s) > eps && totpt <= TOTPT_MAX); //double mpi_stop = MPI_Wtime(); //printf("Total time = %f\n", mpi_stop-mpi_start); return _s; } /////////////////////////////////////////////////////////////////////////////////////////////////////////// FUNC_TYPE Integral1DMK_par2 (CURVLINE f, LIMIT left, LIMIT right, ERR_T eps, PTCOUNT totpt, ULONG mult) /* Вычисление однократного интеграла методом Монте-Карло на кластере MPI. Вариант со средним значением функции. Каждый процесс вычисляет интеграл на уменьшенном интервале */ // П Р О Т Е С Т И Р О В А Н А // НА КЛАСТЕРЕ ! ! ! - 26.05.04 г. { // Проверка корректности введенных значений if (eps <= 0 || mult < 1 || totpt < 10 || totpt >= TOTPT_MAX) exit(1); int totproc, self, root = 0; MPI_Comm comw = MPI_COMM_WORLD; MPI_Comm_size(comw, &totproc); // Получение кол-ва процессов MPI_Comm_rank(comw, &self); // Получение процессом своего номера FUNC_TYPE _s = 0.0, _s_root = 0.0, s; ARG_TYPE h, H, left_lim; // Для каждого из процессов определяем свой интервал интегрирования H = (ARG_TYPE)(right - left)/totproc; left += self*H; right = left + H; //srand(self); //MPI_Barrier(comw); //double mpi_start = MPI_Wtime(); do { s = _s; // Запоминаем предыдущее значение интеграла на промежутке [left, right] _s = 0.0; // Обнуляем сумму // Инициализируем генератор случайных чисел временем srand((UINT)(time(NULL) >> self)); do { // Обеспечиваем малый шаг h = (right - left) / totpt; if (h > 1) totpt *= mult; } while (h > 1); if (!h) exit(1); left_lim = left; // Накапливаем сумму while ((h > 0 && left_lim + h <= right) || (h < 0 && left_lim + h >= right)) { _s += f(random(left_lim, left_lim+h)); left_lim += h; } // Вычисляем интеграл по формуле (3.1.10) _s = (_s/totpt)*(right-left); totpt *= mult; } while (fabsl(s – _s) > eps && totpt <= TOTPT_MAX); MPI_Barrier(comw); MPI_Reduce(&_s, &_s_root, 1, MPI_LONG_DOUBLE, MPI_SUM, root, comw); //double mpi_stop = MPI_Wtime(); //printf("Total time = %f\n", mpi_stop-mpi_start); return _s_root; } ///////////////////////////////////////////////////////////////////////////////////////////// UFUNC_TYPE Square (CURVLINE f, LIMIT left, LIMIT right, ERR_T eps) // Вычисление площади под кривой f(x) по формуле Симпсона // П Р О Т Е С Т И Р О В А Н А ! ! ! – 24.05.04 г. { if (eps <= 0) exit (EPS_ERR); ARG_TYPE _s, s = 0.0; ARG_TYPE t0 = fabsl(f(left)) + fabsl(f(right)); ULONG n=1024; do { _s = s; ARG_TYPE h = (right - left) / (2*n); ARG_TYPE t1 = 0, t2 = 0; // t1 – сумма нечетных, t2 – сумма четных for (ULONG j = 1; j <= 2*n-1; j++) j % 2 == 0 ? t2 += fabsl(f(left+j*h)) : t1 += fabsl(f(left+j*h)); s = (t0 + 4*t1 + 2*t2) * (fabsl(h)/3); n *= 2; // mult == 2 } while (fabsl(s-_s) > eps && n <= ULONG_MAX); return s; } /////////////////////////////////////////////////////////////////////////////////////////////// FUNC_TYPE Integral2DMK (SURF surf, LIMIT left, LIMIT right, CURVLINE bottom, CURVLINE top, unsigned char var, ERR_T eps, PTCOUNT totpt, ULONG mult) // Рассчет двойного интеграла для у -правильной области { // Проверка корректности введенных значений if (eps <= 0 || mult < 1 || totpt < 10 || totpt >= TOTPT_MAX) exit(1); FUNC_TYPE V = 0.0, _s = 0.0, s; clock_t start = clock(), stop; // Запускаем таймер do { s = _s; _s = 0.0; // Инициализируем генератор случайных чисел временем srand((UINT)(time(NULL))); ARG_TYPE x, y; for (ULONG i = 1; i <= totpt; i++) { x = random(left, right); y = random(bottom(x), top(x)); _s += surf(x, y); } _s = _s/totpt; // Среднее значение функции // Вывод на экран промежуточных данных printf("\nValue = %G (total point = %G)\n", _s, totpt); printf ("Current error = %G\n", fabsl(s - _s)); totpt *= mult; } while (fabsl(s – _s) > eps && totpt <= TOTPT_MAX); // Вычисляем площадь области интегрирования FUNC_TYPE field = RegionY(bottom, top, left, right, var, eps); printf("region = %f\n", field); // Вычисляем двойной интеграл V = _s*field; stop = clock(); TotTimePrint(start, stop, NULL, 'd'); return V; } /////////////////////////////////////////////////////////////////////////////////////////////// FUNC_TYPE Integral2DMK_par (SURF surf, LIMIT left, LIMIT right, CURVLINE bottom, CURVLINE top, unsigned char var, ERR_T eps, PTCOUNT totpt, ULONG mult) // Рассчет двойного интеграла для у-правильной области { // Проверка корректности введенных значений if (eps <= 0 || mult < 1 || totpt < 10 || totpt >= TOTPT_MAX) exit(1); FUNC_TYPE V = 0.0, _s = 0.0, _s_root = 0.0, s; int totproc, self, root = 0; MPI_Comm comw = MPI_COMM_WORLD; MPI_Comm_size(comw, &totproc); // Получение кол-ва процессов MPI_Comm_rank(comw, &self); // Получение процессом своего номера do { s = _s; _s = 0.0; // Инициализируем генератор случайных чисел временем srand((UINT)(time(NULL) >> self)); ARG_TYPE x, y; for (ULONG i = 1; i <= totpt; I += totproc) { x = random(left, right); y = random(bottom(x), top(x)); _s += surf(x, y); } MPI_Barrier(comw); MPI_Reduce(&_s, &_s_root, 1, MPI_LONG_DOUBLE, MPI_SUM, root, comw); if (!self) { _s = _s_root/totpt; // Среднее значение функции // Вывод на экран промежуточных данных printf("\nValue = %G (total point = %G)\n", _s, totpt); printf ("Current error = %G\n", fabsl(s - _s)); } totpt *= mult; MPI_Barrier(comw); MPI_Bcast (&_s, 1, MPI_LONG_DOUBLE, root, comw); } while (fabsl(s – _s) > eps && totpt <= TOTPT_MAX); // Вычисляем площадь области интегрирования FUNC_TYPE field = RegionY(bottom, top, left, right, var, eps); printf("region = %f\n", field); // Вычисляем двойной интеграл V = _s*field; return V; } [1] Операция понимается в самом широком смысле этого слова: операция интегрирования, операция решения СЛАУ и т.д.