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


program lapack_ex2 implicit none



Pdf көрінісі
бет5/7
Дата21.11.2019
өлшемі0,58 Mb.
#52198
1   2   3   4   5   6   7
Байланысты:
параллель


program lapack_ex2

implicit none

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

realallocatable :: A(:,:)

integerallocatable :: pivot(:)

integer :: m, n, i, j, lda, info

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

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

read(*,*) m, n

lda = n


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

allocate(pivot(1:m))

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

call random_number(A)

!!

call sgetrf(m, n, A, m, pivot, info)

!!

end program lapack_ex2

2.5.

Процедура

sgeev

вычисления собственных

векторов и собственных значений

Собственным вектором линейного преобразования (матрицы) A назы-

вается такой ненулевой вектор x, что для некоторого λ



∈ RC справедливо

Aλx.



Собственным значением линейного преобразования A называется та-

кое число λ



∈ RC, для которого существует собственный вектор, т.е. урав-

нение Aλимеет ненулевое решение x. Различают также правый Aи

левый yA собственные векторы. Для нахождения собственных векторов и

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

sggev

:

subroutine sggev(jobvl, jobvr, n, A, lda, wr, wi, Vl,



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

character(len=1), intent(in) :: jobvl,

character(len=1), intent(in) :: jobvr,

integerintent(in) :: n,

realdimension(lda,*), intent(inout) :: A,

integerintent(in) :: lda,

realdimension(*), intent(out) :: wr,

realdimension(*), intent(out) :: wi,

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

37

realdimension(ldvl,*), intent(out) :: Vl,



integerintent(in) :: ldvl,

realdimension(ldvr,*), intent(out) :: Vr,

integerintent(in) :: ldvr,

realdimension(*), intent(out) :: work,

integerintent(in) :: lwork,

integerintent(out) :: info

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



jobvl


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

'N'


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

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

'V'

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



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

jobvr


аналогичен

jobvl


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

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



n

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



lda


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

A

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



lda >= max(1,n)

,

но обычно



lda = n

;

wr

и

wr



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

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



Vl

и



Vr

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

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

u(j) = Vl(:,j)

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

сива


Vl

; в случае же если 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-й аргумент).

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

sgeev

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



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

program lapack_ex3

implicit none

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

realallocatable :: A(:,:), Vl(:,:), Vr(:,:)

realallocatable :: 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() = 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

integerintent(in) :: m,

integerintent(in) :: n,

realdimension(lda, *), intent(inout) :: A,

integerintent(in) :: lda,

realdimension(*), intent(out) :: s,

realdimension(ldu, *), intent(out) :: U,

integerintent(in) :: ldu,

realdimension(ldvt, *), intent(out) :: VT,

integerintent(in) :: ldvt,

realdimension(*), intent(out) :: work,

integerintent(in) :: lwork,

integerintent(out) :: info

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



– Параметр

jobu


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

'A'


— все ко-

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

U

,

'S'



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

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

сиве

U

,



'O'

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

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

A

,



'N'

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



– Параметр

jobvt


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

V

T

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



– Параметры

m, n


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

⩾ 0.


– Параметр

A

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



ложению.

– Параметр

lda


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

A

,



lda >= max(1,m)

.

– Параметр

S

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



– Параметр

U

и



Vt

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



T

.

– Параметры

ldu

и

ldvt



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

T

.


40

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



– Параметры

work, lwork, info

— служебные.

!!

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

!!

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

program lapack_ex4

implicit none

character(len=1) :: jobu, jobvt

realallocatable :: 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

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

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



Достарыңызбен бөлісу:
1   2   3   4   5   6   7




©engime.org 2024
әкімшілігінің қараңыз

    Басты бет