М. Н. Геворкян, А. В. Королькова, Д. С. Кулябов Параллельное программирование
жүктеу/скачать 0,58 Mb. Pdf көрінісі
|
параллель
do l = 1,m,1 do i = 1,k,1 c(l,i) = a(l,1)*b(1,i) do j=2,n,1 c(l,i) = c(l,i) + a(l,j)*b(j,i) end do end do end do На первый взгляд, никак нельзя уменьшить объём работы, необходи- мый для перемножения двух матриц. Однако исследователям не удалось доказать минимальность и в результате были найдены другие алгоритмы, умножающие матрицы более эффективно. Рассмотрим один из таких алго- ритмов.
1.7.2. Умножение матриц по алгоритму Винограда Одним из алгоритмов умножения матриц является алгоритм Винограда (Shmuel Winograd). Он даёт небольшой выигрыш в скорости вычисления произведения, но может быть эффективно распараллелен. Основная идея алгоритма весьма проста и заключается в преобразовании суммы
∑
a j i b l j следующим образом: c l i =
∑
(a 2j
+ b l 2j )(a 2j i + b l 2j −1 ) | {z } (1) − m/2 ∑
a 2j −1 i a 2j i | {z } (2)
− m/2 ∑
b l 2j −1 b l 2j | {z
(3) . Количество операций сложения и умножения приведено в следующей таблице: Слагаемое Сложений Умножений (1)
(2)
n(m/2 − 1) nm/2 (3)
k(m/2 − 1) km/2 28 Глава 1. Fortran В слагаемом (2) участвуют только элементы матрицы A, а в слагаемом (3) — только элементы матрицы B. Для каждой матрицы A и B, участвую- щих в умножении, эти слагаемые можно вычислить заранее и сохранить в память (особенно это актуально для больших матриц). При необходимости нахождения произведения A
результату прибавить (2) и (3), вычисленные ранее. Ниже приведён код, реализующий алгоритм Винограда [4]: d = n/2
! Вычисление rowFactor для матрицы A do i = 1,m,1 rowFactor(i) = A(i,1)*A(i,2) do j = 2,d,1 rowFactor(i) = rowFactor(i) + A(i,2*j-1)*A(i,2*j) end do end do ! Вычисление columnFactor для матрицы B do i = 1,k,1 columnFactor(i) = B(1,i)*B(2,i) do j = 2,d,1 columnFactor(i) = columnFactor(i) + B(2*j-1,i)*B(2*j,i)
C(i,j) = -rowFactor(i) - columnFactor(j) do l = 1,d,1 C(i,j) = C(i,j) + (A(i,2*l-1)+B(2*l,j))* (A(i,2*l)+B(2*l-1,j))
C(i,j) = C(i,j) + A(i,n)*B(n,j) end do end do end if Геворкян М. Н., Королькова А. В., Кулябов Д. С. Параллельное программирование 29
Задания для лабораторной работы 1.8.1. Задание №1 Написать программу, решающую квадратное уравнение с корнями лю- бого вида.
Написать программу, вычисляющую числа Фибоначчи с помощью ре- курсивной функции.
Написать подпрограмму, распечатывающую переданную ей в качестве аргумента матрицу (двухмерный массив) с любым числом строк и столб- цов. Для этого воспользуйтесь приёмом, указанным в конце раздела 1.6. Подпрограмму требуется оформить как модуль.
Написать программу перемножения матриц стандартным способом и с помощью алгоритма Винограда. Сравнить производительность этих алго- ритмов и встроенной функции matmul для больших матриц. Текущее время можно получить с помощью процедуры ctime(t)
. 1.9. Содержание отчёта Отчёт должен включать: 1) титульный лист; 2) формулировку цели работы; 3) описание результатов выполнения задания:
графики в зависимости от задания); 4) выводы по каждому заданию лабораторной работы.
1. Какие встроенные типы данных поддерживает Фортран? Как описать вещественную переменную двойной точности? 2. Почему комплексная переменная двойной точности занимает в памяти 16 байт?
30 Глава 1. Fortran 3. Для чего применяются дополнительные атрибуты при описании пере- менной?
4. Как с помощью встроенных функций вычислить дробную часть веще- ственного числа? 5. Как правильно проверить, равны ли две вещественные переменные? 6. Какой оператор прерывает виток цикла и переходит к следующей ите- рации? 7. Чем подпрограмма отличается от функции? 8. Какова структура программы на фортране? 9. В чём преимущество модулей? 10. Необходимо описать 5 массивов размерности (1:100,-100:1). Как сделать это наиболее кратко? 11. Что такое экстент, размер и ранг массива? 12. Что такое сечение и вырезка массива? 13. Какие способы задания значений элементов массива вы знаете? 14. Как называется встроенная функция, вычисляющая определитель мат- рицы (двухмерного массива)? 15. Что такое динамический массив и как он задаётся? 16. Опишите основные операторы для ввода-вывода данных. Как записать данные в файл? 31 2. Библиотека LAPACK 2.1. Назначение библиотеки LAPACK Библиотека LAPACK (Linear Algebra PACKage) служит для решения за- дач линейной алгебры (решение систем линейных алгебраических уравне- ний (СЛАУ), различные матричные разложения). Подпрограммы LAPACK написаны на языке Fortran
, хорошо оптимизированы и для повышения быстродействия используют BLAS (Basic Linear Algebra Subprograms) — подпрограммы, оптимизированные под конкретные архитектуры. Изначально процедуры LAPACK были написаны под Fortran 77 , од-
нако к версии 3.2 (2008) были переписаны под Fortran 90 . Процедуры LAPACK делятся на следующие виды: – процедуры-драйверы (driver routines) — для решения конкретных задач линейной алгебры (решение СЛАУ, нахождение собственных векторов и т.д.);
дачи (различные разложения матриц); – вспомогательные процедуры (auxiliary routines). 2.2. Процедуры LAPACK При компиляции программы, использующей процедуры LAPACK, сле- дует использовать одну из следующих команд: gfortran -L /usr/local/lib -llapack -lblas src.f90 -o prog gfortran
Названия всех процедур закодированы по одинаковой схеме: pmmaaa ,
p задаёт тип чисел (вещественный или комплексный) матрицы и может принимать следующие значения: p Расшифровка s действительный тип одинарной точности real(kind=4) d действительный тип двойной точности real(kind=8) c комплексный тип одинарной точности complex(kind=8) z комплексный тип двойной точности complex(kind=16) Буквы
mm задают тип матрицы, что позволяет использовать более оптими- зированные процедуры для работе с ними. Наиболее часто используемыми типами являются: 32 Глава 2. Библиотека LAPACK mm Расшифровка GE матрица общего вида DI диагональная матрица OR ортогональная матрица (вещественная) UN унитарная матрица (комплексная) GT трёхдиагональная матрица SY симметричная матрица Наконец aa задают конкретный алгоритм, например sv — решение СЛАУ, или trf
— вычисление LUP
-разложения матрицы. Познакомимся теперь с некоторыми конкретными процедурами. Будем рассматривать процедуры для действительных чисел одинарной точности. Однако все рассматриваемые процедуры имеют реализации и для действи- тельных чисел двойной точности, и для комплексных чисел одинарной и двойной точности. 2.3. Процедура sgesv решения СЛАУ Рассмотрим применение процедуры для решения систем линейных ал- гебраических уравнений общего вида:
11
1 + a 12
2 +
13 x 3 + . . .+ a 1n x n =
1
21
1 +
22 x 2 + a 23
3 +
a 2n x n =
2
31
1 +
32 x 2 + a 33
3 +
a 3n x n =
3
.. . .. . .. . . ..
.. . .. . a m1 x 1 + a m2 x 2 + a m3 x 3 + . . .+ a mn x n =
m . Исходными данными, передаваемыми подпрограмме, являются матрица A и вектор b: A =
a 11
12
13
a 1,n −1 a 1n a 21
22
23
a 2,n −1 a 1n a 31
32
33
a 3,n −1 a 1n .. .
. .. . . .. .. . .. .
m −1,1 a m −1,2 a m −1,3 . . . a m −1,n−1 a m −1,n a m1 a m2 a m3 . . . a m,n −1 a mn , b = b 1
2
3 .. . b m −1 b m . Компактно система запишется следующим образом: Ax = b. Как уже говорилось выше, процедура, решающая СЛАУ с квадратной Геворкян М. Н., Королькова А. В., Кулябов Д. С. Параллельное программирование 33 матрицей A[n × n] общего вида и действительными коэффициентами оди- нарной точности называется sgesv и имеет вид: subroutine sgesv(n, nb, A, lda, pivot, B, ldb, info) integer, intent(in) :: n integer, intent(in) :: nb real, dimension(lda,*), intent(inout) :: A integer, intent(in) :: lda integer, dimension(*), intent(out) :: pivot, real, dimension(ldb,*), intent(inout) :: B integer, intent(in) :: ldb integer, intent(out) :: info Расшифруем математический смысл параметров. Параметр Расшифровка n размерность системы уравнений (матрицы A) nb число векторов b A матрица A lda главная размерность массива A (
) pivot
вектор, задающий матрицу перестановок B массив, содержащий вектор(ы) b ldb число правых частей (векторов b) info индикатор успешности выполнения процедуры Следует отметить, что обычно lda = n
и ldb = 1
. Более подробного разъяснения требует pivot . Пусть задана матрица
11
12
13
21
22
23
31
32
33
Известно, что элементарное преобразование, заключающееся в пере- становке первой и второй строки матрицы, можно задать с помощью эле- ментарной матрицы вида 0 1 0 1 0 0 0 0 1 , так как
a 11
12
13
21
22
23
31
32
33 0 1 0 1 0 0 0 0 1 =
21
22
23
11
12
13
31
12
13
Матрицы, задающие перестановку строк или столбцов, называются мат- рицами перестановки (permutation matrix). Они представляют собой матри- цу, каждая строка которой содержит лишь один ненулевой элемент, равный 34 Глава 2. Библиотека LAPACK единице. Поэтому для экономии памяти компьютера данная матрица зада- ётся в виде вектора (одномерного массива). При этом i-й элемент массива равен номеру столбца, где в i-й строке стоит единица. В нашем случае 0 1 0 1 0 0 0 0 1 , p =
1 3 . Ниже приведён пример использования процедуры sgesv .
program lapack_ex1 implicit none !----------------------------------------- real, allocatable :: A(:,:), b(:) integer, allocatable :: pivot(:) integer :: n, i, j, ok ! переменная для задания форматирования character(len=20) :: fmt !----------------------------------------- write(*,*), "Введите размерность СЛАУ" read(*,*) n allocate(A(1:n,1:n)) allocate(b(1:n)) allocate(pivot(1:n)) ! Заполняем числами от 0 до 1 call random_number(A) call random_number(b) !! write(fmt,*), n write(*,*), "----- Матрица A -----" write (*,"("//adjustl(fmt)//"(f8.3), t1)") & ((A(i,j), j = 1,n), i = 1,n) !! call SGESV(n, 1, A, n, pivot, b, n, ok) !! write(*,*), "----- Матрица A -----" write (*,"("//adjustl(fmt)//"(f8.3), t1)") & ((A(i,j), j = 1,n), i = 1,n) write(*,*), "--- Вектор перестановок ---" write (*,"("//adjustl(fmt)//"(i3), t1)") & (pivot(i), i = 1,n) write(*,*), "--- Решение системы ---" write (*,"("//adjustl(fmt)//"(f8.3), t1)") & (b(i), i = 1,n) end program lapack_ex1 Геворкян М. Н., Королькова А. В., Кулябов Д. С. Параллельное программирование 35
Процедура sgetrf вычисления LUP-разложения При работе с матрицами большую роль играют различные матричные разложения. В LAPACK не реализованы, например, такие операции, как на- хождение определителя матрицы и нахождение обратной матрицы. Одна- ко вычислить определитель можно, воспользовавшись LUP-разложением, а найти обратную матрицу — с помощью SVD-разложения. Рассмотрим вначале LUP-разложение и соответствующую процедуру LAPACK.
Любую невырожденную матрицу A можно представить в виде PA = LU, где L — нижнедиагональная матрица, U — верхнедиагональная матрица, P — матрица перестановок (элементарная матрица). Алгоритм LUP-разложения, реализованный процедурой sgetrf
, может обрабатывать любые невырожденные матрицы и при этом обладает высо- кой устойчивостью:
Многие параметры совпадают с параметрами процедуры sgesv :
Расшифровка m число строк матрицы A n число столбцов матрицы A A матрица A lda главная размерность массива A (
) pivot
вектор, задающий матрицу перестановок info
индикатор успешности выполнения процедуры (0 — успех) Результат вычислений (коэффициенты матриц L и U) записывается в массив A
такой записи информация о ней не теряется: A =
u 11
12
13
u 1,n −1 u 1n l 21
22
23
u 2,n −1 u 1n l 31
32
33
u 3,n −1 u 1n .. .
. .. . . .. .. . .. .
m −1,1 l m −1,2 l m −1,3 . . . u m −1,n−1 u m −1,n l m1 l m2 l m3 . . . l m,n −1 u mn . 36 Глава 2. Библиотека LAPACK Матрица P восстанавливается по вектору перестановок p, коэффициен- ты которого записываются в массив pivot .
sgesv .
жүктеу/скачать 0,58 Mb. Достарыңызбен бөлісу:
|