Transcript 第八章：MKL概览

性能库:Intel 数学核心库(MKL)
杨全胜
http://www.njyangqs.com/
东南大学成贤学院计算机系
多核结构与程序设计
内容




介绍
性能特性
使用库
库内容
Southeast University
东 南 大 学
2
http://www.njyangqs.com/
多核结构与程序设计
库内容

性能库: Intel® MKL
 Intel®
数学核心库是一个广泛的科学/工程数学
库
 为Intel® processors优化
 他的多线程可以在SMP的机器上有效的使用
Southeast University
东 南 大 学
3
http://www.njyangqs.com/
多核结构与程序设计
介绍

Intel®数学核心库目的
性能，性能，还是性能!
 Intel公司的工程、科学和金融数学库
 访问:

(BLAS, LAPACK)
 特征向量/特征值求解(BLAS, LAPACK)
 一些量子化学的需要(dgemm)
 偏微分方程, 信号处理, 地震, 固态物理(FFTs)
 常规科学，金融科学[向量超越函数 (VML)和向量随机
数发生器(VSL)]
 解算程序

为Intel的当前和未来处理器做调整
Southeast University
东 南 大 学
4
http://www.njyangqs.com/
多核结构与程序设计
介绍

Intel® 数学核心库目的 – 不能做的
 但是不能用Intel®
做 …
X’
Y’
Z’
W’
数学核心库(Intel® MKL) 来
=
4x4
变换矩阵
X
Y
Z
W
几何变换
MKL
 不要在小n上调用向量数学函数.
 不要在“小”计算上用Intel®
Southeast University
东 南 大 学
5
http://www.njyangqs.com/
多核结构与程序设计
介绍

Intel® 数学核心库内容

BLAS (基本线性代数子例程)

Level 1 BLAS – 向量-向量操作
• 15个函数类型
• 48个函数

Level 2 BLAS – 矩阵-向量操作
• 26个函数类型
• 66个函数

Level 3 BLAS – 矩阵-矩阵操作
• 9个函数类型
• 30个函数

Southeast University
东 南 大 学
Extended BLAS – 对稀疏向量的level 1 BLAS
• 8个函数类型
• 24个函数
6
http://www.njyangqs.com/
多核结构与程序设计
介绍

Intel® 数学核心库内容

LAPACK (线性代数包)

解算程序进而特征值求解. 全部数百个子例程
 有总数超过1000的用户可调用和支持例程

DFTs (离散福利叶变换)
混合基数, 多维变换
 多线程


VML (矢量数学库)
矢量超越函数集
 大多数的libm函数，但是更快


VSL (矢量统计学库)

Southeast University
东 南 大 学
矢量随机数发生器集
7
http://www.njyangqs.com/
多核结构与程序设计
介绍

Intel® 数学核心库内容
 BLAS和LAPACK*
都是Fortran.
 传统的高性能计算
 VSL和VML有Fortran和C接口
 DFTs有Fortran
95和C接口
 cblas接口：更便于C/C++程序员调用BLAS
Southeast University
东 南 大 学
8
http://www.njyangqs.com/
多核结构与程序设计
介绍

Intel® 数学核心库(Intel® MKL) 环境
编译器
库



Linux*
Intel, Gnu
.a, .so
支持32位和64位Intel处理器
一个例子和测试的大集合
大量的文档
Southeast University
东 南 大 学
Windows*
Intel, CVF,
Microsoft
.dll, .lib
9
http://www.njyangqs.com/
多核结构与程序设计
内容




介绍
性能特性
使用库
库内容
Southeast University
东 南 大 学
10
http://www.njyangqs.com/
多核结构与程序设计
性能特性

资源的有限优化


所有优化的目标是速度最大化
有限的资源优化 – 消耗一个或者更多的系统资源:
CPU: 寄存器使用, 浮点单元
 Cache: 尽可能久地保持cache中的数据; 处理cache更迭
 TLBs: 最大限度地使用每页上的数据
 存储带宽: 尽可能少地访问存储器.
 计算机: 用所有可用的处理器执行线程
 系统: 使用所有的可用节点(机群软件).

Southeast University
东 南 大 学
11
http://www.njyangqs.com/
多核结构与程序设计
性能特性

线程化

大部分Intel® 数学核心库(Intel® MKL) 能够线程化，但
是：
内存带宽是有限的资源
 线程化level 1和level 2 BLAS 基本上没效果( O(n) )


有大量的线程化的机会:
Level 3 BLAS ( O(n3) )
 LAPACK* ( O(n3) )
 FFTs ( O(n log(n) )
 VML, VSL ? 依赖处理器和函数



所有的线程化都是使用OpenMP*.
所有Intel MKL都被设计和编译成对线程化是安全的
Southeast University
东 南 大 学
12
http://www.njyangqs.com/
多核结构与程序设计
内容




介绍
性能特性
使用库
库内容
Southeast University
东 南 大 学
13
http://www.njyangqs.com/
多核结构与程序设计
使用库

连接Intel® 数学核心库(Intel® MKL)
ifort, BLAS, IA-32处理器:
 方案1:
 ifort
 方案2:
 f77
myprog.f mkl_c.lib
CVF, LAPACK, IA-32处理器:
myprog.f mkl_s.lib
 方案3:
使用一个在运行时连接的DLL来静态连接
一个C程序:
 link
myprog.obj mkl_c_dll.lib
 注意:基于处理器的优化二进制代码将在运行的时
候执行
Southeast University
东 南 大 学
14
http://www.njyangqs.com/
多核结构与程序设计
使用库

矩阵乘法
 自己做卷积或用点积
Roll Your Own
for( i = 0; i < n; i++ ){
for( j = 0; j < m; j++ ){
for( k = 0; k < kk; k++ )
c[i][j] += a[i][k] * b[k][j]; }}
ddot
for( i = 0; i < n; i++ ){
for( j = 0; j < m; j++ )
c[i][j] = cblas_ddot( n, &a[i], incx,&b[0][j], incy); }
Southeast University
东 南 大 学
15
http://www.njyangqs.com/
多核结构与程序设计
使用库

矩阵乘法
 DGEMV/DGEMM
dgemv
for( i = 0; i < n; i++ )
cblas_dgemv( CBLAS_RowMajor, CBLAS_NoTrans, m, n,
alpha, a, lda, &b[0][i], ldb, beta, &c[0][i], ldc );
dgemm
Cblas_dgemm( CblasColMajor, CblasNoTrans,
CblasNoTrans, m, n, kk, alpha, b, ldb, a, lda, beta, c, ldc );
Southeast University
东 南 大 学
16
http://www.njyangqs.com/
多核结构与程序设计
使用库

实验 1: DGEMM
 比较分别用C源码，DDOT,DGEMG和DGEMM实
现矩阵乘法的性能
 练习在MKL/BLAS中控制线程的能力
Southeast University
东 南 大 学
17
http://www.njyangqs.com/
多核结构与程序设计
内容




介绍
性能特性
使用库
库内容
Southeast University
东 南 大 学
18
http://www.njyangqs.com/
多核结构与程序设计
库内容

LAPACK*中优化Intel® 数学核心库
 很多重要的LAPACK优化:
– 有效地使用多处理器和多核
 递归分解
 线程化
• 减少串行时间(Amdahl’s law: t = t串行 + t并行/p)
• 代码中进一步扩展阻塞到整个列的块，以提高局部性
 无需运行时库的支持
 NETLIB
LAPACK有很多调用函数提高了对运行时
库函数的支持，这些函数都在Intel MKL中，所以
不再需要运行时库，可以和任何编译器方便的链接、
Southeast University
东 南 大 学
19
http://www.njyangqs.com/
多核结构与程序设计
库内容

离散福利叶变换(DFTs)
 1维，2维，3维…（受编程语言限制）
 支持多线程
 为内存的有效利用而进行混合基数变换，并满足
许多物理问题的需要。可以指定任何规模的变换，
但并不是所有规模的变换运行的同样快
 用户指定比例，变换的符号
 嵌入矩阵的变换，一次允许调用多个一维变换
 C和F90接口
Southeast University
东 南 大 学
20
http://www.njyangqs.com/
多核结构与程序设计
库内容

用Intel® 数学核心库做离散福利叶变换
 基本上是3个步骤
 创建一个描述器.
•
Status = DftiCreateDescriptor(MDH, …)
 提交描述器(实例化它).
•
Status = DftiCommitDescriptor(MDH)
 执行转换.
•
Status = DftiComputeForward(MDH, X)
 释放描述器（可选地）
Southeast University
东 南 大 学
21
http://www.njyangqs.com/
多核结构与程序设计
库内容

矢量数学库(VML)的功能/问题
 矢量数学库: 矢量超越函数 – 像libm，但是更好
（更快）
 接口: 有Fortran和C两个接口
 多种精确度
< 1 ulp )
 较低精度, 更快( < 4 ulps )
 高精度(
sin(0),等
 错误处理– 这里不能复制libm
 特殊值的处理√(-a),
Southeast University
东 南 大 学
22
http://www.njyangqs.com/
多核结构与程序设计
库内容

VML: 是怎样的?
 它对金融代码很重要(蒙托卡罗模拟法).
 指数,
对数
 依赖于超越函数的其他科学代码
 误差函数能够在一些代码中起到很重要的作用.
Southeast University
东 南 大 学
23
http://www.njyangqs.com/
多核结构与程序设计
库内容

矢量统计学库(VSL)
 随机数发生器组(RNGs)
 大量的非均匀分布
 VSL广泛地使用在各种变换上
– 一些函数
 用户能够提供自己的BRNG或变换
 5个基本的RNGs (BRNGs) – 位，整数，浮点数
 并行计算的支持
 MCG31,
Southeast University
东 南 大 学
R250, MRG32, MCG59, WH
24
http://www.njyangqs.com/
多核结构与程序设计
库内容

非均匀RNGs
 Gaussian
(two methods) 高斯
 Exponential 指数
 Laplace 拉普拉斯变换
 Weibull 韦伯分布
 Cauchy 柯西分布
 Rayleigh 瑞利
 Lognormal 对数正态
 Gumbel Gumbel分布
Southeast University
东 南 大 学
25
http://www.njyangqs.com/
多核结构与程序设计
库内容

使用VSL
 基本是三个步骤
 创建一个流指针
• VSLStreamStatePtr stream;
 创建一个流.
• vslNewStream(&stream,VSL_BRNG_MC_G31,
seed );
 生成一组RNGs.
• vsRngUniform( 0, &stream, size, out, start,
end );
 删除一个流(可选).
• vslDeleteStream(&stream);
Southeast University
东 南 大 学
26
http://www.njyangqs.com/
多核结构与程序设计
库内容

实验: 用蒙特卡罗法计算π
 根据圆面积公式，S=πr2
 看一个单位圆，其中，1/4个
1
x
单位圆的面积是单位矩形面积
y
的一部分，单位矩形面积为1，
-1
1
现在在单位矩形内产生大量
0
随机的点，则落在1/4圆内的
点所占的百分比就是1/4的单位
圆面积。
-1
 一个点是否在1/4单位圆内的判断方法就是该点的
坐标是否满足 x2+y2≤1
Southeast University
东 南 大 学
27
http://www.njyangqs.com/
多核结构与程序设计
库内容
int main(){
 实验: 用蒙特卡罗法计算π
unsigned int iter=200000000;
int i,j;
double x, y;
double dUnderCurve=0.0;
double pi=0.0;
srand( 0);
for (i=0;i<iter;i++) {
x=(double)rand()/(double)RAND_MAX;
y=(double)rand()/(double)RAND_MAX;
if (x*x + y*y <= 1.0) {
dUnderCurve++;
}
}
pi = dUnderCurve / (double) iter * 4 ;
return 0;
Southeast University
}
东 南 大 学
1
x
-1
y
0
1
-1
28
http://www.njyangqs.com/
多核结构与程序设计
库内容
unsigned int iter=200000000;
 int i,j;
double x, y;
double dUnderCurve=0.0;
double pi=0.0;
1
double r[BLOCK_SIZE*2];
x
VSLStreamStatePtr stream;
vslNewStream( &stream, BRNG, (int) clock() );
y
for(j=0;j<iter/BLOCK_SIZE; j++) {
-1
1
vdRngUniform( METHOD, stream, BLOCK_SIZE*2, r, 0.0, 1.0 );
0
for (i=0; i<BLOCK_SIZE; i++) {
x=r[i];
y=r[i+BLOCK_SIZE];
if (x*x + y*y <= 1.0) {
dUnderCurve++;
-1
}
}
}
vslDeleteStream( &stream );
Southeast University
29
pi = dUnderCurve / (double) iter * 4 ;
http://www.njyangqs.com/
东 南 大 学
实验: 用蒙特卡罗法计算π