Лекция 11 Математические библиотеки * * Оптимизация приложений с использованием библиотеки 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 ;

Sparse Solvers - решатели систем линейных уравнений с невырожденной разреженной матрицей; VML (Vector Math Library) – математические функции над векторами; VSL (Vector Statistical Library) – библиотека статистики; PDEs (Partial Differential Equations) – решатели уравнений в частных производных; Optimization Solvers – подпрограммы нелинейной оптимизации. Компоненты MKL

Расположение модулей ( *.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 и др. Настройка MS Visual Studio 2010 Расположение библиотек ( *.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 : Векторные операции y = ax + y a = x·y Уровень 2 : Векторно-матричные операции y = aAx + by A = ax T ·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

Packed Storage – хранение симметричных, Эрмитовых или треугольных матриц : верхняя или нижняя треугольные части упаковываются по столбцам в одномерный массив. uplo ='U', ---> arr(i+j*(j-1)/2), i arr(i+(2*n-j)*(j-1)/2), i >= j Хранение матриц

Band Storage - ненулевые диагонали, матрицы размером (m x n), хранятся построчно в двумерном массиве arrB размером (kl+ku+1 x n), где ku – число диагоналей выше главной диагонали kl – число диагоналей ниже главной диагонали Хранение матриц

При выполнении LU -факторизации над ленточными матрицами в массиве arrB следует добавить kl дополнительных строк. * - элементы не используемые в вычислениях + - используются для хранения результатов Хранение матриц

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

Пример. Подготовить матрицу для LU-факторизации. ku = 1, kl = 1, m = 900, n=900. arrB(2*kl+ku+1 x n) = (4 х 900). Хранение матриц cba cba cba 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

Шаблон: ( ) s - real, single precision (вещественные данные одинарной точности) c - complex, single precision (комплексные данные одинарной точности) d - real, double precision (вещественные данные двойной точности) z - complex, double precision (комплексные данные двойной точности) Именование функций

Примеры ddot - скалярное произведение двух векторов двойной точности sgemv - матрично-векторное произведение, матрица обычная, одинарная точность ztrmm - матрично-матричное произведение, треугольная матрица, комплексные числа двойной точности Именование функций

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 = Для продолжения нажмите любую клавишу... Результат вычислений Пример Вычисление скалярного умножения

Решить одномерное уравнение теплопроводности методом конечных разностей. Самостоятельно разобрать процедуры gbtrf, gbtrs. Пример

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 ! --- верхняя Вариант программы (1)

!----- формируем 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 Вариант программы (2)