20150420_msu08_Integr (Лекции)

МГУ ВМКЛекция 8Параллельные алгоритмычисленного интегрированияЯкобовский Михаил Владимировичпроф., д.ф.-м.н.Институт прикладной математикиим. М.В.Келдыша РАН, МоскваМосква, 2015 г.Параллельные алгоритмы: Параллельные алгоритмы численногоинтегрирования © Якобовский М.В.1 из 40Постановка задачиВычислить с точностью  значение определенного интегралаBJ  A, B    f x dxAПусть на отрезке [A,B] задана равномерная сетка, содержащая n+1 узел:xi  A B Ai, i  0,, nnx0=A x1xn1=BТогда, согласно методу трапеций, можно численно найтиопределенный интеграл от функции на отрезке [A,B]:J n  A, B  B  A  f x0  n 1f xn    f xi  n  22 i 1Будем полагать значение J найденным с точностью ,если выполнено условие:J n1  J n 2   J n 2n2  n1x0=A x1Москва, 2015 г.Параллельные алгоритмы: Параллельные алгоритмы численногоинтегрирования © Якобовский М.В.xn2=B2 из 40Последовательный алгоритм интегрированияIntTrap01(A,B){n=1J2n=(f(A)+f(B))(B-A)/2do{Jn= J2nx0=An=2ns=f(A)+f(B)for(i=1;i<n;i++)s+=2f(A+(B-A)i/n);J2n=s(B-A)/n;}while(|J2n- Jn| J2n)return J2n}Москва, 2015 г.x1x2=BНедостатки:- в некоторых точках значениеподынтегральной функции вычисляетсяболее одного раза- на всем интервале интегрированияиспользуется равномерная сетка, тогда какчисло узлов сетки на единицу длины наразных участках интервалаинтегрирования, необходимое длядостижения заданной точности, зависит отвида функцииПараллельные алгоритмы: Параллельные алгоритмы численногоинтегрирования © Якобовский М.В.3 из 40Пример функции1000011f x   2 sin 2  , 0  A  1x x8000BB11 1 21J  A, B    2 sin 2  dx   sin x2x 4  x  A xA600040001 B A2 2 J  A, B    2 sin   sin  4  ABB A 20000.020.030.040.05Результаты вычисления интеграла на разных отрезках [A,B]A0.000010.00010.0010.010.1110Москва, 2015 г.B0.00010.0010.010.1110100Npoints1 553 5681 726 123360 07579 973105 108396412289903831845653149331eps real-2.77E-111.90E-102.05E-11-2.22E-128.67E-11-6.00E-11-6.30E-11time, c434.55470.9974.1216.4421.420.0940.096Параллельные алгоритмы: Параллельные алгоритмы численногоинтегрирования © Якобовский М.В.4 из 40Адаптивный алгоритмmain(){J= IntTrap03(}A, B, f(A),f(B))IntTrap03(A,B,fA,fB){J=0C=(A+B)/2fC=f(C)Преимущества:- нет повторных вычисленийфункции- малое число вычисленийна гладких участкахsAB=(fA+fB)*(B-A)/2sAC=(fA+fC)*(C-A)/2sCB=(fC+fB)*(B-C)/2sACB=sAC+sCBx0=Ax1=Cx2=BНедостаток:координаты концов отрезковхранятся в программном стекепроцесса и фактическинедоступны программистуМосква, 2015 г.if(|sAB-sACB||sACB|)J=IntTrap03(A,C,fA,fC)+IntTrap03(C,B,fC,fB)elseJ=sACBreturn J}Параллельные алгоритмы: Параллельные алгоритмы численногоинтегрирования © Якобовский М.В.5 из 40Метод локального стекаIntTrap04(A,B){J=0C=(A+B)/2fA=f(A)fB=f(B)sAC=(fA+fC)*(C-A)/2sCB=(fC+fB)*(B-C)/2sAB=(fA+fB)*(B-A)/2sACB=sAC+sCBwhile(1){Обработка отрезка}return J}if(|sAB-sACB||sACB|){PUT_INTO_STACK( A, C, fA, fC, sAC)A=CfA=fCsAB=sCB}else{J+=sACBACfC=f(C)BABif(STACK_IS_FREE)breakGET_FROM_STACK( A, B, fA, fB, sAB)}Москва, 2015 г.Параллельные алгоритмы: Параллельные алгоритмы численногоинтегрирования © Якобовский М.В.6 из 40Процедуры и данные обеспечивающие стек// данные, описывающие стек// указатель вершины стекаsp=0// массив структур в которых// хранятся отложенные заданияstruct{A,B,fA,fB,s}stk[1000]// макроопределения доступа к стеку#define STACK_IS_FREE (sp==0)#define PUT_INTO_STACK(A,B,fA,fB,s){stk[sp].A=Astk[sp].B=Bstk[sp].fA=fAstk[sp].fB=fBstk[sp].s=ssp++}#define GET_FROM_STACK(A,B,fA,fB,s){sp-A=stk[sp].AB=stk[sp].BfA=stk[sp].fAfB=stk[sp].fBs=stk[sp].s}Москва, 2015 г.Параллельные алгоритмы: Параллельные алгоритмы численногоинтегрирования © Якобовский М.В.7 из 40К вопросу о времени выполненияТестирование показало, что при расчете спомощью алгоритма локального стекаIntTrap04 время работы было меньше,примерно на 5%, чем при использованииIntTrap03Примем алгоритм IntTrap04 за «наилучший»последовательныйМосква, 2015 г.Параллельные алгоритмы: Параллельные алгоритмы численногоинтегрирования © Якобовский М.В.8 из 40Параллельный алгоритм интегрированияМетод геометрического параллелизма? Метод коллективного решения??Москва, 2015 г.Параллельные алгоритмы: Параллельные алгоритмы численногоинтегрирования © Якобовский М.В.9 из 40Метод геометрического параллелизмаmain(){…for(i=0;i<p;i++)StartParallelProcess( IntTrap04, A+(B-A)*i/p, A+(B-A)*(i+1)/p, &(s[i]) )WaitAllParallelProcessJ=0for(i=0;i<p;i++)J+=s[i]}Москва, 2015 г.Недостаток:Значительный дисбалансзагрузки процессоровПараллельные алгоритмы: Параллельные алгоритмы численногоинтегрирования © Якобовский М.В.10 из 40Расчет интеграла на разных отрезках10000Недостаток:Большой дисбаланс загрузкипроцессоров800012 1 f x   2 sin  x x600040002000A=1e5p,0.02C0.030.040.05Bинтервал 1интервал2время1, свремя2, с(B-C)/ (C-A)101001 00010 000100 000Москва, 2015 г.[1e-5, 0.10000900000][1e-5, 0.01000990000][1e-5, 0.00100999000][1e-5, 0.00010999900][1e-5, 0.00001999990][0.10000900000, 1][0.01000990000, 1][0.00100999000, 1][0.00010999900, 1][0.00001999990, 1]37.67937.27436.98934.06418.869Параллельные алгоритмы: Параллельные алгоритмы численногоинтегрирования © Якобовский М.В.0.0040.0370.3693.36418.82211 из 40Метод коллективного решенияmain(){// Порождение p параллельных процессов,// каждый из которых выполняет процедуру slavefor(k=0;k<p;k++)StartParallel(slave #k)i=0 // число переданных для обработки интервалов// n – число отрезков интегрированияfor(k=0;k<p;k++){ // Передача концов отрезков интегрированияSend(slave k, A+(B-A)*i/n, A+(B-A)*(i+1)/n)i++}// J - значение интеграла на всем интервале [A,B]J=0n=?Недостаток:- Либо большой дисбаланс загрузкипроцессоров- Либо большой объем лишнихвычисленийМосква, 2015 г.Пока есть отрезки, не переданные для отработки,следует дождаться сообщения от любого из процессов slave,вычислившего частичную сумму на переданном ему отрезке,Получить значение этой суммы, прибавить к общему значениюИнтеграла и передать освободившемуся процессу очереднойотрезокwhile(i<n){k = Recv(slave ANY, s)J+=sSend(slave k, A+(B-A)*i/n, A+(B-A)*(i+1)/n)i++}// Получить результаты вычислений переданных отрезков// и прибавить их к общей суммеfor(k=0;k<p;k++){Recv(slave any, s)J+=s}}slave(){подчиненный процесс, вычисляющийзначение интеграла на отрезке[a,b]while(1){}Параллельные алгоритмы: Параллельные алгоритмы численногоинтегрирования © Якобовский М.В.Recv(main,a,b)s=IntTrap04(a,b)Send(main,s)}12 из 40ВыводПрактически непригодны для решенияпоставленной задачи методыгеометрического параллелизма(статическая балансировка)иколлективного решения(динамическая балансировка)Москва, 2015 г.Параллельные алгоритмы: Параллельные алгоритмы численногоинтегрирования © Якобовский М.В.13 из 40Метод глобального стекаВычислительные системы с общей памятью Динамическая балансировка загрузки Отсутствие централизованного управленияИнициализация данныхВычисленияПорождение pпараллельных процессовПроцесс 1Процесс 2Процесс 3Процесс pОжидание завершениявыполнения запущенныхпроцессовЗавершение работыМосква, 2015 г.Параллельные алгоритмы: Параллельные алгоритмы численногоинтегрирования © Якобовский М.В.14 из 40Стеки алгоритмаГлобальный стекЛокальные стекиIntTrap04pМосква, 2015 г.IntTrap04pIntTrap04p IntTrap04pПараллельные алгоритмы: Параллельные алгоритмы численногоинтегрирования © Якобовский М.В.15 из 40Стеки алгоритмаГлобальный стекЛокальные стекиТолько стек НикакогопроцессанетIntTrap04pМосква, 2015 г.IntTrap04pIntTrap04p IntTrap04pПараллельные алгоритмы: Параллельные алгоритмы численногоинтегрирования © Якобовский М.В.16 из 40Стеки алгоритмаГлобальный стекЛокальные стекиIntTrap04pМосква, 2015 г.IntTrap04pIntTrap04p IntTrap04pПараллельные алгоритмы: Параллельные алгоритмы численногоинтегрирования © Якобовский М.В.17 из 40Идея алгоритмаПусть есть доступный всем параллельным процессам список отрезков интегрирования,организованный в виде стека.

Назовем его глобальным стеком. Пусть у каждого процесса есть свой, доступный только этому процессу локальный стек Перед запуском параллельных процессов в глобальный стек помещаетсяединственная запись (в дальнейшем "отрезок"):– координаты концов отрезка интегрирования,– значения функции на концах,– приближенное значение интеграла на этом отрезке. Тогда, предлагается следующая схема алгоритма, выполняемого каждым изпараллельных процессов:Пока в глобальном стеке есть отрезки:- взять один отрезок из глобального стека- выполнить алгоритм локального стека, но,если при записи в локальный стек в нем есть несколько отрезков, а в глобальном стекеотрезки отсутствуют, то:- переместить часть отрезков из локального стека в глобальный стек.Москва, 2015 г.Параллельные алгоритмы: Параллельные алгоритмы численногоинтегрирования © Якобовский М.В.18 из 40Вопросыкакую часть отрезков следует перемещать излокального стека в глобальный стек?- в какой момент интеграл вычислен?- что должен делать процесс у которого пустлокальный стек, если глобальный стек тоже пуст?-- должен ли процесс закончить работу, если в еголокальном и в глобальном стеке отрезков нет?Москва, 2015 г.Параллельные алгоритмы: Параллельные алгоритмы численногоинтегрирования © Якобовский М.В.19 из 40Схема Интегрирующего процессаIntTrap04p(){// цикла обработки стека отрезковwhile(sdat.ntask>0){// чтение одного интервала из списка интерваловsdat.ntask-- // указатель глобального стекаGET_FROM_GLOBAL_STACK[sdat.ntask](a,b,fa,fb,sab)ИНТЕГРИРОВАНИЕ ОДНОГО ОТРЕЗКА}sdat.s_all = sdat.s_all + s}Москва, 2015 г.Параллельные алгоритмы: Параллельные алгоритмы численногоинтегрирования © Якобовский М.В.20 из 40Правильное определение общей суммыmain(){.Sem_init(sdat.sem_sum,1) //доступ к глобальной сумме открыт.}IntTrap04p(){…// Начало критической секции сложения частичных суммsem_wait(sdat.sem_sum)sdat.s_all = sdat.s_all + ssem_post(sdat.sem_sum)// Конец критической секции сложения частичных сумм}Москва, 2015 г.Параллельные алгоритмы: Параллельные алгоритмы численногоинтегрирования © Якобовский М.В.21 из 40Схема Интегрирующего процессаIntTrap04p(){// цикла обработки стека отрезковwhile(sdat.ntask>0){// чтение одного интервала из списка интерваловsdat.ntask-- // указатель глобального стекаGET_FROM_GLOBAL_STACK[sdat.ntask](a,b,fa,fb,sab)ИНТЕГРИРОВАНИЕ ОДНОГО ОТРЕЗКА}sem_wait(sdat.sem_sum)sdat.s_all = sdat.s_all + ssem_post(sdat.sem_sum)}Москва, 2015 г.Параллельные алгоритмы: Параллельные алгоритмы численногоинтегрирования © Якобовский М.В.22 из 40Схема интегрирования отрезкаwhile(1) // интегрирование одного отрезка{c=(a+b)/2;fc=f(c)sac=(fa+fc)*(c-a)/2scb=(fc+fb)*(b-c)/2sacb=sac+scbif(!BreakCond(sacb,sab)) // Точность на части отрезка достигнута{s+=sacbif(sp==0) break; // локальный стек пуст, выходsp--; GET_FROM_LOCAL_STACK[sp]( a, b, fa, fb, sab)}else{PUT_INTO_LOCAL_STACK[sp]( a, c, fa, fc, sac); sp++a=cfa=fcsab=scb}if((sp>SPK) && (!sdat.ntask)) // перемещение части локального стека в общий список интервалов{while((sp>1) && (sdat.ntask<sdat.maxtask)){sp--; GET_FROM_LOCAL_STACK[sp](a,b,fa,fb,sab)PUT_INTO_GLOBAL_STACK[sdat.ntask](a,b,fa,fb,sab); sdat.ntask++}}}Москва, 2015 г.Параллельные алгоритмы: Параллельные алгоритмы численногоинтегрирования © Якобовский М.В.23 из 40Схема интегрирования отрезка// интегрирование одного отрезкаwhile(1){ИнициализацияТочность на части отрезка достигнута?Добавлять отрезки в глобальный стек?}Москва, 2015 г.c=(a+b)/2;fc=f(c)sac=(fa+fc)*(c-a)/2scb=(fc+fb)*(b-c)/2sacb=sac+scbПараллельные алгоритмы: Параллельные алгоритмы численногоинтегрирования © Якобовский М.В.24 из 40Схема интегрирования отрезкаif(!BreakCond(sacb,sab)){ // Точность на части отрезка достигнутаs+=sacbif(sp==0) break; // локальный стек пуст, выход// интегрирование одного отрезкаwhile(1){ИнициализацияТочность на части отрезка достигнута?Добавлять отрезки в глобальный стек?}sp--;GET_FROM_LOCAL_STACK[sp]( a, b, fa, fb, sab)}else{PUT_INTO_LOCAL_STACK[sp]( a, c, fa, fc, sac);sp++a=cfa=fcsab=scb}Москва, 2015 г.Параллельные алгоритмы: Параллельные алгоритмы численногоинтегрирования © Якобовский М.В.25 из 40Схема интегрирования отрезка// интегрирование одного отрезкаwhile(1){ИнициализацияТочность на части отрезка достигнута?Добавлять отрезки в глобальный стек?}if((sp>SPK) && (!sdat.ntask)) // перемещение частилокального стека в общий список интервалов{while((sp>1) && (sdat.ntask<sdat.maxtask)){sp--;GET_FROM_LOCAL_STACK[sp](a,b,fa,fb,sab)PUT_INTO_GLOBAL_STACK[sdat.ntask](a,b,fa,fb,sab);sdat.ntask++}}Москва, 2015 г.Параллельные алгоритмы: Параллельные алгоритмы численногоинтегрирования © Якобовский М.В.26 из 40Преждевременное окончание работы процессаwhile(1){// Начало критической секции чтения из глобального// стека очередного интервала интегрированияsem_wait(sdat.sem_list)if(sdat.ntask0){sem_post(sdat.sem_list) // разрешить другим процессам// доступ к глобальному стекуbreak}sdat.ntask-- // указатель глобального стекаGET_FROM_GLOBAL_STACK[sdat.ntask](a,b,fa,fb,sab)sem_post(sdat.sem_list)// Конец критической секции чтения из глобального// стека очередного интервала интегрирования…}Москва, 2015 г.Параллельные алгоритмы: Параллельные алгоритмы численногоинтегрирования © Якобовский М.В.27 из 40Преждевременный выходПлохое условие выхода из цикла обработкистека интервалов Интегрирующие процессы не должнызаканчивать работу до тех пор, пока все отрезкиинтервала интегрирования не будут полностьюобработаны Преждевременное завершение работы приведетк получению верного ответа, но за большеевремяМосква, 2015 г.Параллельные алгоритмы: Параллельные алгоритмы численногоинтегрирования © Якобовский М.В.28 из 40Если глобальный и локальный стеки пустыОтрезок интегрирования может находиться в нескольких состояниях:- находится в глобальном стеке интервалов;- обрабатывается некоторым интегрирующим процессом;- находится в локальном стеке интервалов некоторого процесса;- полностью обработан: известно значение интеграла на этом отрезке и оноприбавлено к локальной частичной сумме соответствующего процесса."Время жизни" отрезка, после того, как некоторый процесс начал его обработку,относительно невелико - отрезок разбивается на две части и перестаетсуществовать, породив два новых отрезка.