Параллельные алгоритмы вычислительной алгебры. Разделение переменных

Содержание

Слайд 2

Часть 5: Разделение переменных Разделение переменных для оператора Лапласа Вычислительный алгоритм OMP и MPI реализация

Часть 5: Разделение переменных

Разделение переменных для оператора Лапласа
Вычислительный алгоритм
OMP и MPI

реализация
Слайд 3

Сеточная аппроксимация оператора Лапласа Кусочно-линейная аппроксимация, квадратная сетка

Сеточная аппроксимация оператора Лапласа

Кусочно-линейная аппроксимация, квадратная сетка

Слайд 4

Сеточная аппроксимация оператора Лапласа Как найти собственные вектора и собственные значения матрицы А?

Сеточная аппроксимация оператора Лапласа

Как найти собственные вектора и собственные значения матрицы

А?
Слайд 5

Базисные функции оператора Лапласа Пусть - собственный вектор матрицы A. Тогда

Базисные функции оператора Лапласа

Пусть - собственный вектор матрицы A.
Тогда оказывается:
Следовательно для

нашей матрицы A:
Слайд 6

Базисные функции оператора Лапласа Разложим теперь решение задачи и правую часть

Базисные функции оператора Лапласа

Разложим теперь решение задачи и правую часть по

базисным функциям матрицы A
Подставим полученные выражения в изначальное уравнение получим и учтя собственные значения матрицы A:
Слайд 7

Вычислительный алгоритм (двойной ряд) - Тригонометрическое разложение, трудоемкость С*N*log(N)

Вычислительный алгоритм (двойной ряд)

- Тригонометрическое разложение, трудоемкость С*N*log(N)

Слайд 8

Вычислительный алгоритм (двойной ряд) - Тригонометрическое разложение, трудоемкость С*N*log(N)

Вычислительный алгоритм (двойной ряд)

- Тригонометрическое разложение, трудоемкость С*N*log(N)

Слайд 9

Вычислительный алгоритм (однократный ряд) - Серия из трехдиагональных матриц, трудоемкость решения 6*N

Вычислительный алгоритм (однократный ряд)

- Серия из трехдиагональных матриц, трудоемкость решения 6*N

Слайд 10

Вычислительный алгоритм (однократный ряд)

Вычислительный алгоритм (однократный ряд)

Слайд 11

Вычислительный алгоритм (однократный ряд) Разложение по синусам горизонтальных компонент. Присутствует в библиотеках Intel MKL, Netlib

Вычислительный алгоритм (однократный ряд)

Разложение по синусам горизонтальных компонент. Присутствует в библиотеках

Intel MKL, Netlib
Слайд 12

Вычислительный алгоритм (однократный ряд) Решение трехдиагональных систем (каждая система разная!!! к диагонали добавляютя различные собственные числа)

Вычислительный алгоритм (однократный ряд)

Решение трехдиагональных систем (каждая система разная!!! к диагонали

добавляютя различные собственные числа)
Слайд 13

Вычислительный алгоритм (однократный ряд) Обратное разложение по синусам горизонтальных компонент. Присутствует в тех же библиотеках

Вычислительный алгоритм (однократный ряд)

Обратное разложение по синусам горизонтальных компонент. Присутствует в

тех же библиотеках
Слайд 14

Вычислительный алгоритм (однократный ряд) subroutine laplas_solution(f(*,*),N) …… Do i=1,N Call forward_trig_transform(f(i,*),….)

Вычислительный алгоритм (однократный ряд)

subroutine laplas_solution(f(*,*),N)
……
Do i=1,N
Call forward_trig_transform(f(i,*),….)
enddo
Do j=1,N
Call three_diagonal_lu(f(*,j),lambda(j),…)
enddo
Do

i=1,N
Call backward_trig_transform(f(i,*),….)
enddo
…..
End subroutine
Слайд 15

Вычислительный алгоритм (однократный ряд) (OMP) subroutine laplas_solution(f(*,*),N) …… C$OMP PARALLEL DO

Вычислительный алгоритм (однократный ряд) (OMP)

subroutine laplas_solution(f(*,*),N)
……
C$OMP PARALLEL DO
Do i=1,N
Call

forward_trig_transform(f(i,*),….)
enddo
C$OMP END PARALLEL DO
C$OMP PARALLEL DO
Do j=1,N
Call three_diagonal_lu(f(*,j),lambda(j),…)
enddo
C$OMP END PARALLEL DO
C$OMP PARALLEL DO
Do i=1,N
Call backward_trig_transform(f(i,*),….)
enddo
C$OMP END PARALLEL DO
…..
End subroutine
Слайд 16

Вычислительный алгоритм (однократный ряд) ( MPI) Для MPI идеалогии такое неподходит,

Вычислительный алгоритм (однократный ряд) ( MPI)

Для MPI идеалогии такое неподходит, так

как данные распределены по различным процессам, следовательно невозможно иметь либо LU на одном процессе (если данные распределены по строкам), либо разложение по синусам (если данные распределены по столбцам), либо оба (если данные распределены “шахматкой”).
Слайд 17

Вариант 1 (транспонирование, MPI) Пусть данные распределены по строкам Proc 0 Proc 1 Proc 2

Вариант 1 (транспонирование, MPI)

Пусть данные распределены по строкам

Proc 0

Proc 1

Proc 2

Слайд 18

Вариант 1 (транспонирование, MPI) Так как каждая строка целиком лежит на

Вариант 1 (транспонирование, MPI)

Так как каждая строка целиком лежит на одном

процессе, то каждый процесс независимо может сделать разложение по синусам

Proc 0

Proc 1

Proc 2

Слайд 19

Вариант 1 (транспонирование, MPI) Обычный алгоритм прогонки не может быть реализован,

Вариант 1 (транспонирование, MPI)

Обычный алгоритм прогонки не может быть реализован, так

как компоненты правых частей лежат на разных процессах

Proc 0

Proc 1

Proc 2

Слайд 20

Вариант 1 (транспонирование, MPI) Данные можно транспонировать между процессами. ВАЖНО! Не

Вариант 1 (транспонирование, MPI)

Данные можно транспонировать между процессами.
ВАЖНО! Не забыть про

порядок нумерации на процессах

Proc 0

Proc 1

Proc 2

Proc 0

Proc 1

Proc 2

MPI_Alltoallv

Слайд 21

Вариант 1 (транспонирование, MPI) После транспонирования данных правые части для прогонок

Вариант 1 (транспонирование, MPI)

После транспонирования данных правые части для прогонок лежат

на одном процессе, следовательно может быть выполнен обычный алгоритм

Proc 0

Proc 1

Proc 2

Слайд 22

Вариант 1 (транспонирование, MPI) Данные необходимо вернуть в изначальное расположения для

Вариант 1 (транспонирование, MPI)

Данные необходимо вернуть в изначальное расположения для релизации

обратного разложения по синусам

Proc 0

Proc 1

Proc 2

Proc 0

Proc 1

Proc 2

MPI_Alltoallv

Слайд 23

Вариант 1 (транспонирование, MPI) Теперь каждая строка опять целиком лежит на

Вариант 1 (транспонирование, MPI)

Теперь каждая строка опять целиком лежит на одном

процессе, и каждый процесс независимо может сделать обратное разложение по синусам

Proc 0

Proc 1

Proc 2

Слайд 24

Вариант 1 (транспонирование, MPI) subroutine laplas_solution(f(*,*),N) …… Call MPI_Init(); Do i=ny_first,ny_last

Вариант 1 (транспонирование, MPI)

subroutine laplas_solution(f(*,*),N)
……
Call MPI_Init();
Do i=ny_first,ny_last
Call forward_trig_transform(f(i-ny_first+1,*),….)
enddo
call transpose _y_to_x(f,f_work,…)
call

MPI_All_to_allv(f_work,f,…,MPI_COMM_WORLD)
Do j=nx_first,nx_last
Call three_diagonal_lu(f(*,j-nx_first+1),lambda(j),…)
enddo
call MPI_All_to_allv(f,f_work,…,MPI_COMM_WORLD)
call transpose _x_to_y(f_work,f,…)
Do i=ny_first,ny_last
Call backward_trig_transform(f(i-ny_first+1,*),….)
Enddo
…..
Call MPI_Finalize();
End subroutine
Слайд 25

Вариант 2 (параллельная прогонка, MPI) Так как обмены данных может быть

Вариант 2 (параллельная прогонка, MPI)

Так как обмены данных может быть достаточно

затратной операций, предлагается аглоритм параллельной прогонки (прогонка для правой части, распределенной между процессами)

Proc 0

Proc 1

Proc 2

Слайд 26

Вариант 2 (параллельная прогонка, MPI) Разобьем компоненты вектора решения (X)i и

Вариант 2 (параллельная прогонка, MPI)

Разобьем компоненты вектора решения (X)i и правой

части (F)I между процессами таким образом, чтоб на каждом процессоре лежала часть вектора с компонентами [k,k+M], каждая компонента хранилась только на одном процессе.
Слайд 27

Вариант 2 (параллельная прогонка, MPI) Теорема: Если вектор решения записать в

Вариант 2 (параллельная прогонка, MPI)

Теорема: Если вектор решения записать в следующем

виде:
,где max_proc – число процессов, то решение задачи с трехдиагональной матрицей записывается в следующем виде:
Слайд 28

Вариант 2 (параллельная прогонка, MPI) Теорема: Если вектор решения записать в

Вариант 2 (параллельная прогонка, MPI)

Теорема: Если вектор решения записать в следующем

виде:
,где max_proc – число процессов, то решение задачи с трехдиагональной матрицей записывается в следующем виде:
Относительно Wk записывается трехдиагональная система уравнений, с коэффициентами, зависящими от a,b,c,P,Q.
Слайд 29

Вариант 2 (параллельная прогонка, MPI) Алгоритм: Решаем на каждом процессора 3

Вариант 2 (параллельная прогонка, MPI)

Алгоритм:
Решаем на каждом процессора 3 системы уравнений

методом прогонки с одной и той же матрицей, разными правыми частями
Рассылаем с каждого процесса малый объем данных на один процесс для того, чтоб собрать трехдиагональную систему уравнений размерностью=число процессов
Полученные компоненты решения рассылаем на каждый процесс (по два значения на процесс) и собираем на них требуемые компоненты решения
Слайд 30

Вариант 1 (транспонирование, MPI) subroutine laplas_solution(f(*,*),N) …… Call MPI_Init(); Do i=ny_first,ny_last

Вариант 1 (транспонирование, MPI)

subroutine laplas_solution(f(*,*),N)
……
Call MPI_Init();
Do i=ny_first,ny_last
Call forward_trig_transform(f(i-ny_first+1,*),….)
enddo
Do j=1,nx
Call

MPI_three_diagonal_lu(f(*,j),lambda(j),…)
enddo
Do i=ny_first,ny_last
Call backward_trig_transform(f(i-ny_first+1,*),….)
Enddo
…..
Call MPI_Finalize();
End subroutine