Математические библиотеки*
Презентацию к лекции Вы можете скачать здесь.
Обзор
Intel ® Math Kernel Library Основная математическая библиотека
Производительность (для больших задач)! |
Параллельное выполнение! |
Оптимизирована под процессоры Intel®! |
IMSL International Mathematical and Statistical Library Международная математическая библиотека подпрограмм.
Компоненты IMSL
- Процедуры линейной алгебры;
- Решение систем линейных уравнений;
- Преобразования Фурье и Лапласа;
- Решение нелинейных уравнений;
- Задачи оптимизации;
- Сортировка и поиск данных;
- Интерполяция и аппроксимация;
- Интегрирование и дифференцирование;
- Решение ОДУ;
- Решение ДУЧП.
Компоненты MKL
BLAS | - базовые подпрограммы линейной алгебры; |
Sparse BLAS | - BLAS для разреженных матриц; |
LAPACK | - решатели задач линейной алгебры; |
ScaLAPACK | - параллельное расширение LAPACK; |
DFT | – дискретные преобразования Фурье; |
Cluster DFT | – кластерные DFT; |
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 mkl_intel_thread.lib libiomp5md.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\ComposerXE- 2011\mkl\examples |
Возможности на примере BLAS
Уровень 1 : Векторные операции

Уровень 2 : Векторно-матричные операции

Уровень 3 : Матричные операции

Хранение векторов
Элементы хранятся в одномерном массиве
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.
Остальные элементы необязательны.
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
Band Storage - ненулевые диагонали, матрицы размером (m x n), хранятся построчно в двумерном массиве arrB размером (kl+ku+1 x n), где
При выполнении LU-факторизации над ленточными матрицами в массиве arrB следует добавить kl дополнительных строк.

После LU-факторизации входной массив arrB будет изменён следующим образом:
Пример. Подготовить матрицу для LU-факторизации.

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
Именование функций
Шаблон:
<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.

Пример



Вариант программы (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