program lapack_ex2 implicit none

Геворкян М. Н., Королькова А. В., Кулябов Д. С. Параллельное программирование

37

real, dimension(ldvl,*), intent(out) :: Vl,

integer, intent(in) :: ldvl,

real, dimension(ldvr,*), intent(out) :: Vr,

integer, intent(in) :: ldvr,

real, dimension(*), intent(out) :: work,

integer, intent(in) :: lwork,

integer, intent(out) :: info

Поясним смысл передаваемых параметров:

–

jobvl

может принимать два значения:

'N'

— левый собственный век-

тор матрицы A не вычисляется и

'V'

— левый собственный вектор мат-

рицы A вычисляется;

–

jobvr

аналогичен

jobvl

, но задаёт настройки вычисления для правого

собственного вектора;

–

n

— размерность матрицы A;

–

lda

— главная размерность массива

A

; в общем случае

lda >= max(1,n)

,

но обычно

lda = n

;

–

wr

и

wr

— массивы, содержащие действительную и мнимую части вы-

численных собственных значений;

–

Vl

и

Vr

— левый и правый собственные векторы; в случае если j-е соб-

ственное значение действительно, то

u(j) = Vl(:,j)

, j-й столбец мас-

сива

Vl

; в случае же если j-е и j + 1-е собственные значения являются

комплексно сопряжённой парой, то

u(j) = Vl(:,j) + i*Vl(:,j+1)

и

u(j+1) = Vl(:,j) - i*Vl(:,j+1)

; то же справедливо и для

Vr

;

–

ldvl

и

ldvr

— главные размерности

Vl

и

Vr

; требуется, чтобы

ldvl,ldvr

>= n

;

–

work

и

lwork

— служебный массив, необходимый для вычислений, и

его размерность (для хорошей производительности рекомендуют брать

lwork > = 4*n

);

–

info

— индикатор успешности выполнения (может принимать следую-

щие значения: 0 — выполнение успешно, > 0 — ошибка и = i < 0, где

i означает, что неверно задан i-й аргумент).

Пример использования процедуры

sgeev

приведён ниже.

! Нахождение собственных значений

program lapack_ex3

implicit none

!---------------------------------------------------

real, allocatable :: A(:,:), Vl(:,:), Vr(:,:)

real, allocatable :: wr(:), wi(:), work(:)

character(len=1) :: jobvl, jobvr

integer :: n, i, j, lda, ldvl, ldvr, lwork, info

!---------------------------------------------------

write(*,*), "Введите размерность квадратной матрицы n"

read(*,*) n

lda = n

ldvl = n

ldvr = n

38

Глава 2. Библиотека LAPACK

jobvl = 'V'

jobvr = 'V'

lwork = 5*n

allocate(A(1:n,1:n))

allocate(Vl(1:n,1:n))

allocate(Vr(1:n,1:n))

allocate(work(1:lwork))

allocate(wr(1:n))

allocate(wi(1:n))

! Заполняем числами от 0 до 1

call random_number(A)

!!

call sgeev(jobvl, jobvr, n, A, lda, wr, wi, &

Vl, ldvl, Vr, ldvr, work, lwork, info)

!!

if (info .eq. 0) then

write(*,*) "Успешное выполнение"

else if(info .gt. 0) then

write(*,*) "Ошибка!"

else if(info .lt. 0) then

write(*,*) "Неправильно задан параметр №", abs(info)

end if

end program lapack_ex3

2.6.

Процедура

sgesvd

вычисления SVD-разложения

Ещё одним полезным разложением матриц, часто используемым на прак-

тике, является сингулярное разложение или, иначе, SVD-разложение (Singu-

lar Value Decomposition). Любая матрица M порядка m

× n, элементы ко-

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

M = UΣV

†

,

где U — унитарная матрица порядка m

× m, Σ — диагональная матри-

ца порядка m

× n с неотрицательными вещественными коэффициентами

на диагонали, V — унитарная матрица порядка n

× n, V

†

— эрмитово-

сопряжённая матрица к V.

Элементы σ

ii

диагонали матрицы Σ называются сингулярными числа-

ми матрицы M и определены с точностью до их перестановки. Обычно их

выстраивают по убыванию σ

1

≥ σ

2

≥ σ

3

≥ . . . ≥ σ

n

. Столбцы матриц U и

V называют соответственно левым и правым сингулярными векторами. В

случае действительной матрицы M матрицы U и V ортогональны.

SVD-разложение можно использовать для вычисления ранга матрицы,

нахождения обратной и псевдообратной матриц, нахождения приближён-

ного решения переопределённой системы уравнений и т.д.

Геворкян М. Н., Королькова А. В., Кулябов Д. С. Параллельное программирование

39

Значение ранга матрицы равно количеству ненулевых сингулярных чи-

сел. В случае квадратной матрицы, если хотя бы одно значение σ

ii

= 0, то

матрица вырожденная (det(M ) = 0).

Для вычисления обратной матрицы вначале следует удостовериться,

что матрица невырожденная. Если это так, то обратная матрица вычисля-

ется по формуле

A

−1

= VΣ

−1

U

T

,

Σ

−1

= diag

(

1

σ

11

, . . . ,

1

σ

nn

)

.

Рассмотрим теперь процедуру библиотеки

LAPACK

, осуществляющую

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

ми, то эрмитово сопряжение заменяется на простую транспозицию:

sgesvd (jobu, jobvt, m, n, A, lda, s, U, ldu, Vt,

ldvt, work, lwork, info)

character*1, intent(in) :: jobu

character*1, intent(in) :: jobvt

integer, intent(in) :: m,

integer, intent(in) :: n,

real, dimension(lda, *), intent(inout) :: A,

integer, intent(in) :: lda,

real, dimension(*), intent(out) :: s,

real, dimension(ldu, *), intent(out) :: U,

integer, intent(in) :: ldu,

real, dimension(ldvt, *), intent(out) :: VT,

integer, intent(in) :: ldvt,

real, dimension(*), intent(out) :: work,

integer, intent(in) :: lwork,

integer, intent(out) :: info

Поясним смысл передаваемых процедуре аргументов:

– Параметр

jobu

может принимать четыре значения:

'A'

— все m ко-

лонок матрицы U возвращаются в массиве

U

,

'S'

— первые min(m, n)

столбцов матрицы U (левые сингулярные векторы) возвращаются в мас-

сиве

U

,

'O'

— первые min(m, n) столбцов матрицы U (левые сингу-

лярные векторы) возвращаются в перезаписываемом массиве

A

,

'N'

—

столбцы U не вычисляются.

– Параметр

jobvt

задаёт точно такие же параметры, но для строк матрицы

V

T

(правые сингулярные векторы).

– Параметры

m, n

задают число строк и столбцов матрицы A, m, n

⩾ 0.

– Параметр

A

— изначальная матрица, которую следует подвергнуть раз-

ложению.

– Параметр

lda

— главная размерность

A

,

lda >= max(1,m)

.

– Параметр

S

— двухмерный массив, соответствующий матрице Σ.

– Параметр

U

и

Vt

— это матрицы U и V

T

.

– Параметры

ldu

и

ldvt

— главные размерности U и V

T

.

40

Глава 2. Библиотека LAPACK

– Параметры

work, lwork, info

— служебные.

!!

Использование процедуры для SVD разложения

!!

действительной матрицы одинарной точности

program lapack_ex4

implicit none

character(len=1) :: jobu, jobvt

real, allocatable :: A(:,:), U(:,:), Vt(:,:), s(:), work(:)

integer :: m, n, lda, ldu, ldvt, lwork, info, i, j

!---------------------------------------------------------

write(*,*), "Введите размерности матрицы m и n"

read(*,*) m, n

lda = m

ldu = m

ldvt = m

jobu = 'A'

jobvt = 'A'

! Служебная переменная, чем больше, тем лучше

lwork = 20*max(m,n)

allocate(A(1:m,1:n))

allocate(U(1:m,1:m))

allocate(Vt(1:n,1:n))

! Диагональные элементы сингулярной матрицы

allocate(s(1:min(n,m)))

! Служебный массив

allocate(work(1:lwork))

! Заполняем числами от 0 до 1

call random_number(A)

!!

call sgesvd(jobu, jobvt, m, n, A, lda, s, &

U, ldu, Vt, ldvt, work, lwork, info)

!!

end program lapack_ex4

2.7.

Задания для лабораторной работы

2.7.1.

Задание № 1

Проверьте работоспособность четырёх вышеприведённых примеров. До-

бавьте в программы код, организующий распечатку результатов. Проверь-

те правильность получаемых результатов на простых примерах матриц 3

×3

(т.е. требуется задать массив вручную и вычислить решение СЛАУ, соб-

ственные векторы, матрицы разложения и т.д.).

Геворкян М. Н., Королькова А. В., Кулябов Д. С. Параллельное программирование

41

2.7.2.

Задание № 2

Вычислить определитель матрицы, используя

LAPACK

тремя разными

способами: с помощью собственных векторов матрицы, с помощью SVD

и LUP разложений. Какой способ наиболее прост? Какой способ наиболее

быстр?

2.7.3.

Задание № 3

Вычислить обратную или псевдообратную матрицу с помощью SVD

разложения.

2.8.

Содержание отчёта

Отчёт должен включать:

1) титульный лист;

2) формулировку цели работы;

3) описание результатов выполнения задания:

– листинги программ;

– результаты выполнения программ (текст, снимок экрана, таблицы и

графики в в зависимости от задания);

4) выводы по каждому заданию лабораторной работы.

2.9.

Контрольные вопросы

Для ответов на некоторые вопросы следует воспользоваться докумен-

тацией с официального сайта

LAPACK http://www.netlib.org/

.

1. Как расшифровывается аббревиатура

LAPACK

?

2. Что такое

BLAS

и зачем он использован в

LAPACK

?

3. Процедуры

LAPACK

являются подпрограммами или функциями?

4. Какова схема названий процедур библиотеки

LAPACK

?

5. Какие языки программирования поддерживает

LAPACK

? На каком языке

написан сам

LAPACK

?

6. Как будет именоваться процедура по вычислению собственных векто-

ров для комплексной трёхдиагональной матрицы?

7. Есть ли в

LAPACK

процедура для вычисления определителя? Если есть,

то приведите её название и аргументы.

8. Есть ли в

LAPACK

процедура для решения переопределённой или неопре-

делённой системы алгебраических линейных уравнений? Если есть, то

приведите её название и аргументы.

9. В описании аргументов процедуры

sgesvd

об аргументе

jobvt

было

сказано кратко, так как он аналогичен

jobu

. Дайте его подробное опи-

сание по аналогии с

jobu

.

10. Что делает процедура

sgeqrf

? Дайте краткое описание.

42

3. OpenMP

3.1.

Введение

OpenMP (Open Multi-Processing) — открытый стандарт для создания

многопоточных программ. Поддерживаются языки

Fortaran

,

C

и

C++

.

OpenMP является реализацией многопоточности (multithreading). Мно-

гопоточность означает, что выполняемый процесс может состоять из несколь-

ких потоков (thread), выполняемых одновременно. Термин thread также

переводят буквально как «нить».

Все потоки выполняются в адресном пространстве одного процесса, из

чего следует, что многопоточная технология программирования в общем и

OpenMP в частности может применяться на многопроцессорных системах с

разделяемой памятью (например, на компьютерах с многоядерным процес-

сором). Это обстоятельство определяет как достоинства многопоточности,

так и недостатки.

К достоинствам можно отнести:

– время, требуемое для создания и уничтожения нити в рамках процесса,

мало по сравнению с временем, требуемым на создание и уничтожение

полноценного нового процесса;

– относительная простота программирования многопоточных программ;

– эффективное использование многоядерных процессоров.

К недостаткам можно отнести:

– многопоточная программа может работать лишь на системах с разде-

ляемой памятью (грубо говоря, в рамках одного компьютера, пусть и

оснащённого несколькими процессорами);

– сложность адаптации последовательных алгоритмов для многопоточно-

го выполнения.

Ввиду малых ресурсов, необходимых для создания и уничтожения по-

токов, в течение работы программы потоки могут создаваться и уничто-

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

один главный поток (master thread), который создаётся в самом начале и не

уничтожается в течение всего времени его работы. Главный поток создаёт

набор подчинённых (slave) потоков и задача распределяется между ними.

3.2.

OpenMP и Fortran

Для компиляции программы, использующей OpenMP директивы, необ-

ходимо указать опцию

–fopenmp

gfortran –fopenmp src-file.f90 -o prog-name

Геворкян М. Н., Королькова А. В., Кулябов Д. С. Параллельное программирование

43

где

src-file.f90

— имя исходного файла с кодом программы,

prog-name

— имя исполняемого файла. В случае если не указать опцию

-o

, испол-

няемый файл получит стандартное имя

a.out

. Для проверки поддержки

OpenMP можно использовать простейшую программу следующего вида:

program ex1

!$ print *, "OpenMP поддерживается!"

end program ex1

Области кода, которые следует выполнять параллельно, выделяются с

помощью специальных директив препроцессора соответствующего языка —

прагм. В случае языка Fortran все директивы в общем случае выглядят сле-

дующим образом:

!$omp <директива> [опции]

<код>

!$omp end <директива>

Следует также отметить, что программирование с помощью OpenMP

реализовано в рамках идеологии одна программа — множество данных

(Single Program Multiple Data, SPMD). В рамках этой идеологии для двух

нитей используется одинаковый код. В идеальном случае код параллель-

ной и последовательной версий программ должен отличаться лишь дирек-

тивами OpenMP, однако это случается нечасто.

Кроме директив предпроцессору в библиотеке OpenMP реализованы

некоторые вспомогательные функции, позволяющие получать отладочную

информацию. Для использования этих функций следует подключить заго-

ловочный файл

omp_lib.h

.

Рассмотрим применение функций для замера времени выполнения участ-

ков кода:

program ex2

implicit none

include "omp_lib.h"

real :: t1, t2, tick

t1 = omp_get_wtime()

call sleep(1)

t2 = omp_get_wtime()

write(*,*) "Спали ", t2-t1, " сек."

write(*,*) "Точность таймера: ", omp_get_wtick()

end program ex2

Функции

integer :: omp_get_wtime()

integer :: omp_get_wtick()

возвращают соответственно системное время в секундах и разрешение тай-

мера в секундах, которое можно рассматривать как меру точности таймера.

44

Глава 3. OpenMP

3.3.

Параллельные и последовательные области

При запуске программы порождается нить мастер, которая начинает

выполнение программы. Как только выполнение доходит до директивы

parallel

, порождаются дополнительные нити (число которых зависит от

настроек системы и заданных опций).

Директива, задающая начало и конец параллельной области, выглядит

следующим образом:

!$omp parallel [опции]

<код параллельной области>

!$omp end parallel

Рассмотрим следующий пример:

program ex3

implicit none

include "omp_lib.h"

integer :: i, n, tn

real :: t1, t2

write(*,*) "Последовательная область"

!$omp parallel num_threads(4)

write(*,*) "Параллельная область"

!$omp end parallel

write(*,*) "Последовательная область"

end program

ex3

В данном примере создаётся одна параллельная область. Опция

num_threads

позволяет явно указать количество создаваемых нитей. Ко-

личество нитей, создаваемых по умолчанию, задаётся переменной окруже-

ния

OMP_NUM_THREADS

. Оптимальное количество нитей зависит как от кон-

фигурации системы, так и от алгоритма, и в каждом случае подбирается

индивидуально.

Все нити, порождённые при входе в параллельную область, начнут вы-

полнять указанные в ней команды. Например, при выполнении программы

ex3

:

gfortran -fopenmp ex3.f90 && ./a.out

в консоли будут распечатаны следующие сообщения:

Последовательная область

Параллельная область

Параллельная область

Параллельная область

Параллельная область

Последовательная область

Кроме опции

num_threads

с директивой

parallel

можно использо-

вать следующие опции:

жүктеу/скачать 0,58 Mb.

Достарыңызбен бөлісу: