здесь

Report
Лекция 11
Математические
библиотеки *
* www.intuit.ru
Оптимизация приложений с использованием библиотеки
Intel Math Kernel Library.
Обзор
Intel ® Math Kernel Library
Основная математическая библиотека
Производительность (для больших задач) !
Параллельное выполнение !
Оптимизирована под процессоры Intel® !
IMSL
International Mathematical and Statistical
Library
Международная математическая библиотека
подпрограмм.
Компоненты IMSL
Процедуры линейной алгебры;
Решение систем линейных уравнений;
Преобразования Фурье и Лапласа;
Решение нелинейных уравнений;
Задачи оптимизации;
Сортировка и поиск данных;
Интерполяция и аппроксимация;
Интегрирование и дифференцирование;
Решение ОДУ;
Решение ДУЧП.
Компоненты MKL
BLAS - базовые подпрограммы линейной алгебры;
Sparse BLAS - BLAS для разреженных матриц;
LAPACK - решатели задач линейной алгебры;
ScaLAPACK - параллельное расширение LAPACK;
DFT – дискретные преобразования Фурье;
Cluster DFT – кластерные DFT;
Компоненты MKL
Sparse Solvers - решатели систем линейных
уравнений с невырожденной разреженной матрицей;
VML (Vector Math Library) – математические
функции над векторами;
VSL (Vector Statistical Library) – библиотека
статистики;
PDEs (Partial Differential Equations) –
решатели уравнений в частных производных;
Optimization Solvers – подпрограммы нелинейной
оптимизации.
Настройка MS Visual Studio 2010
Расположение модулей (*.mod файлы)
C:\Program Files\Intel\ComposerXE-
2011\mkl\include\ia32
Библиотеки
mkl_intel_c.lib
libiomp5md.lib
mkl_intel_thread.lib
+
mkl_blas95.lib
mkl_lapack95.lib и др.
Расположение библиотек (*.lib файлы)
C:\Program Files\Intel\ComposerXE-
2011\mkl\lib\ia32
Настройки проекта
Используем MKL Parallel
Настройки проекта
Библиотеки (в зависимости от процедур)
Документация и примеры
Документация пользователя (User guide)
C:\Program Files\Intel\ComposerXE-2011
\Documentation\en_US\mkl\mkl_userguide\
index.htm
Справочные страницы (Manual pages)
C:\Program Files\Intel\ComposerXE-2011
\Documentation\en_US\mkl\mkl_manual\
index.htm
Примеры
C:\Program Files\Intel\ComposerXE2011\mkl\examples
Возможности на примере BLAS
Уровень 1 : Векторные операции
y = ax + y
a = x·y
Уровень 2 : Векторно-матричные операции
y = aAx + by
A = axT·y + A
Уровень 3 : Матричные операции
C = aAB + C
Хранение векторов
Элементы хранятся в одномерном массиве
n - количество элементов;
incx – шаг
incx > 0, хранение по возрастанию;
incx < 0, хранение по убыванию;
incx = 0, все элементы равны значению x(1)
x = (1.0, 3.0, 5.0, 7.0, 9.0, 11.0, 13.0)
Если
incx = 2, n = 3 --> (1.0, 5.0, 9.0)
incx = -2,n = 4 --> (13.0, 9.0, 5.0, 1.0)
incx = 0, n = 5 --> (1.0, 1.0, 1.0, 1.0, 1.0)
Хранение матриц
Full Storage – элементы матрицы хранятся в
двумерном массиве arr(i,j).
Для треугольных, симметричных или Эрмитовых матриц
задание элементов массива определяется параметром
uplo.
uplo = 'L' - нижняя часть, с диагональю, i >= j
uplo = 'U' - верхняя часть, с диагональю, i <= j
Остальные элементы необязательны.
character :: uplo = 'L'
real a(3,3)
a(1,:) = (/1,0,0/)
a(2,:) = (/2,5,0/)
a(3,:) = (/3,0,9/)
Хранение матриц
Packed Storage – хранение симметричных, Эрмитовых
или треугольных матриц : верхняя или нижняя
треугольные части упаковываются по столбцам в
одномерный массив.
uplo ='U', ---> arr(i+j*(j-1)/2), i <= j
uplo ='L', ---> arr(i+(2*n-j)*(j-1)/2), i >= j
2

5
7

3

5
8
9
3
7
9
9
8
3

3
8

3 
2
5 7 3 8 9 3 9 8 3
Хранение матриц
Band Storage - ненулевые диагонали, матрицы
размером (m x n), хранятся построчно в двумерном
массиве arrB размером (kl+ku+1 x n), где
ku – число диагоналей выше главной диагонали
kl – число диагоналей ниже главной диагонали
 a11

 a21
a
 31
 0
 0

a12
0
0
a22
a32
a23
a33
0
a34
a42
a43
a44
0
a53
a54
0 

0 
0 

a45 
a55 
 

 a11
a
 21
a
 31
a12
a23
a34
a22
a33
a44
a32
a43
a54
a42
a53

a45 

a55 
 

 
Хранение матриц
При выполнении LU-факторизации над ленточными
матрицами в массиве arrB следует добавить kl
дополнительных строк.
 a11

 a21
a
 31
 0
 0

a12
0
0
a22
a32
a23
a33
0
a34
a42
a43
a44
0
a53
a54
0 

0 
0 

a45 
a55 
 *

 *
 *

 a11
a
 21
a
 31
*
*

*


a12
a22
a23
a33
a34
a44
a32
a43
a54
a42
a53
*
 

 
a45 

a55 
* 
* 
* - элементы не используемые в вычислениях
+ - используются для хранения результатов
Хранение матриц
После LU-факторизации входной массив arrB будет
изменён следующим образом:
 *

 *
 *

 u11
m
 21
m
 31
*
*
u14
*
u13
u24
u12
u22
u23
u33
u34
u44
m32
m43
m54
m42
m53
*
arrB
u25 

u35 
u45 

u55 
* 
* 
 u11 u12

 0 u 22
 0
0

0
 0
 0
0

u13
u14
u 23
u33
u 24
u34
0
u 44
0
0
0 

u 25 
u35 

u 45 
u55 
верхняя треугольная
матрица
ui , j - элементы верхней треугольной матрицы;
mi , j
- множители факторизации.
Хранение матриц
Пример. Подготовить матрицу для LU-факторизации.
1

a
0

 ...
0

0
0 0 0 0 0 0

b c 0 0 0 0
a b c 0 0 0 

0 0 0 a b c

0 0 0 0 0 1
0

0
1

a

0
0
b
a
0
c
b
a
...
...
...
...
0
c
b
a
0
c
b
0
ku = 1, kl = 1, m = 900, n=900.
arrB(2*kl+ku+1 x n) = (4 х 900).
integer, parameter :: m = 900
real :: arrB(4,m) = 0, a,b,c
arrB(2,3:m)
= c ! вторая, первую строку пропускаем
arrB(3,2:m-1) = b ! третья
arrB(4,1:m-2) = a ! четвертая
arrB(3,1)=1; arrB(3,m)=1
0

c
1

0 
Именование функций
Шаблон: <character> <name> <mod> ( )
<character>
s
- real, single precision
(вещественные данные одинарной точности)
c
- complex, single precision
(комплексные данные одинарной точности)
d
- real, double precision
(вещественные данные двойной точности)
z
- complex, double precision
(комплексные данные двойной точности)
Именование функций
Примеры
ddot <d> <dot> - скалярное произведение двух
векторов двойной точности
sgemv <s> <ge> <mv> - матрично-векторное
произведение, матрица обычная, одинарная точность
ztrmm <z> <tr> <mm> - матрично-матричное
произведение, треугольная матрица, комплексные
числа двойной точности
Пример
Вычисление скалярного умножения
program prog
use mkl95_blas
integer, parameter :: n = 5
! количество элементов
real :: a(n) = (/1,2,3,4,5/),& ! векторы
b(n) = (/2,4,6,8,9/)
real res
res = dot(a,b) ! скалярное произведение
write(*,*) "Scalar product = ", res
end
Результат вычислений
Scalar product =
105.0000
Для продолжения нажмите любую клавишу . . .
Пример
Решить одномерное уравнение теплопроводности
методом конечных разностей.
Самостоятельно разобрать процедуры gbtrf, gbtrs.
T  2T
 2 , T 0  0, T 1  5t.
t
x
T1  0, TN  5t , i  1, N
2
1
t
x
t
3
N 1
N
x
Пример
n 1
n 1
i 1
n 1
n 1
i 1
 Ti
T  2  Ti  T

t
x 2
1
1  n 1
1
1 n
 2
n 1
n 1
 2  Ti 1   2    Ti  2  Ti 1   Ti
x
t 
x
t
 x
Ti
n
a
b
1

a
0

0
 ...

0

0
c
0 0 0 0 0 0

b c 0 0 0 0
a b c 0 0 0

0 a b c 0 0


0 0 0 a b c

0 0 0 0 0 1
Вариант программы (1)
program teplo_1D
use lapack95
use f95_precision
integer, parameter :: N = 64
integer k
real DL(N-1), D(N), DU(N-1) ! DL - нижняя, D - основная, DU - верхняя
real AB(4,N)
real :: UN(N) = 0
real :: dt = 0.01
real, parameter :: H = 1.0
real :: dh = H/(N-1)
! --- шаг по времени
! --- длина области
! шаг по области
integer IPIV(N), INFO
!======================================================================
open(1,file = "res1.txt")
open(3,file = "res2.txt")
open(5,file = "res3.txt")
!------ заполняем диагонали
!------ так как на границах Г.У. первого рода
DL = -dt/dh**2;
DL(N-1) = 0
! --- нижняя
D = 2*dt/dh**2+1;
D(1) = 1; D(N) = 1
! --- основная
DU = -dt/dh**2;
DU(1) = 0
! --- верхняя
Вариант программы (2)
!----- формируем AB-матрицу
AB(1,:) = 0
! не заполняем, нужна для процедур
AB(2,2:N) = DU
AB(3,1:N) = D
AB(4,1:N-1) = DL
call gbtrf(AB, ipiv = IPIV) ! ---- делаем факторизацию
do k = 1,500
UN(N) = 5*k*dt ! ----- правая часть
call gbtrs(AB, UN, IPIV, info = INFO) ! --- вызываем решатель
if ((k == 100).OR.(k == 300).OR.(k == 500)) then
do i = 1,N
x = (i-1)*dh
write(k/100,*) x, UN(i)
end do
close(k/100)
end if
end do
close(1)
close(3)
close(5)
end

similar documents