单指令多数据（SIMD）基础
4.88/5 (9投票s)
使用 C# 通过 Intel MKL 指令集。
通过托管代码访问 SIMD 指令：Intel MKL
前言：试着耐心理解一些基础知识和背景
实数 3E-5 以科学计数法的形式书写。它的意思是 3*10-5（或 10 的负 5 次方乘以 3）。计算机是整数机器，只能通过复杂的编码来表示实数。最流行的实数表示编码是 IEEE 浮点标准。浮点一词来源于小数点前后的位数没有固定；也就是说，小数点可以浮动。也有固定小数点位数表示的，称为定点表示。总的来说，浮点表示比定点表示速度慢且精度低，但它们可以处理更广泛的数字。请注意，计算机能表示的大多数浮点数只是近似值。使用浮点值编程的一个难点在于确保这些近似值能产生合理的结果。由于浮点数运算需要大量的计算能力，许多微处理器都配有专门用于执行浮点运算的芯片，称为浮点单元（FPU）。FPU 也称为数学协处理器或数字协处理器。
由于早期的微处理器实际上没有任何浮点能力，它们只处理整数。浮点计算是在独立的专用硬件上完成的，通常是数学协处理器的形式。但是，随着集成电路芯片技术的发展，计算行业发现晶体管的尺寸可以减小——因此可以在物理芯片内的半导体材料上蚀刻更多的晶体管。晶体管尺寸的减小使得浮点单元可以直接集成在主 CPU 芯片上。将这些单元添加到主 CPU 芯片（作为中央处理单元运行的物理核心）中，就增加了硬件和浮点指令。最终引入了一组扩展的（操作码）指令，称为 SIMD。
如今，现代微处理器可以对多个数据执行相同的指令。这称为单指令多数据（SIMD）。SIMD 指令处理浮点实数，也能为算法提供重要的速度提升。由于 SIMD 指令的执行单元通常属于物理核心，因此可以并行执行的 SIMD 指令数量与可用的物理核心数量相同。如前所述，并行使用这些矢量处理能力可以在某些算法中提供重要的速度提升。将 SIMD 指令及硬件添加到现代多核 CPU 中，比添加浮点能力要复杂得多。自从诞生以来，微处理器一直是 SISD 设备（单指令流，单数据流）。SIMD 因其基本组织单位是矢量而被称为矢量处理。
常规 CPU 一次处理一个标量。（超标量 CPU 一次处理多个标量，但对每条指令执行不同的操作。）另一方面，矢量处理器会将一整行这些同类型的标量排列起来，并将它们作为一个单元进行处理。因此，让我们来看一个例子，它将帮助我们理解 SIMD 指令的强大之处。下图代表了 PABSD 指令。这条指令是 Intel Core 2 架构引入的增强型流式 SIMD 扩展 3 (SSSE 3) 的一部分。
PABSD 的助记符表示“打包双字绝对值”。这条汇编指令接收一个 128 位输入参数，其中包含四个 32 位有符号整数。该指令返回一个 128 位输出，其中包含这四个 32 位有符号整数的绝对值，打包在 128 位输出中。您只需一次调用 PABSD 指令即可计算四个 32 位有符号整数的绝对值。如果您需要计算 1,000 个 32 位有符号整数的绝对值，可以使用 250 次此指令调用来完成，而不是为每个 32 位有符号整数使用一条指令。您可以实现非常重要的速度提升。但是，由于在调用 SIMD 指令之前需要打包数据，并在之后解包输出，因此测量这些会增加一些代码的开销也很重要。如果您只需要计算四个 32 位有符号整数的绝对值，额外的开销会降低整体速度提升。然而，如果您需要计算 100 个 32 位有符号整数的绝对值，通常会从使用此类 SIMD 指令中获益。
大多数现代微处理器都可以执行 SIMD 指令，但这些指令属于不同的扩展指令集。在撰写本文时，最先进的 Intel CPU 支持以下 SIMD 指令集：
- MMX - 多媒体扩展
- SSE - 流式 SIMD 扩展
- SSE2 - 流式 SIMD 扩展 2
- SSE3 - 流式 SIMD 扩展 3
- SSSE3 - 增强型流式 SIMD 扩展 3
- SSE4.1 - 流式 SIMD 扩展 4.1
- SSE4.2 - 流式 SIMD 扩展 4.2
- AES-NI - 高级加密标准新指令
- AVX - 高级矢量扩展
好的。那么，我如何确定我的操作系统支持哪些 SIMD 指令？我正在一台运行 Intel Core 2 Duo 处理器的 SONY VAIO 笔记本电脑上撰写这篇文章。我只需访问 www.cpuid.com 下载 CPU-Z，这是 CPUID 开发的一款免费实用工具。这是运行 CPU-Z 获得的输出。
在此示例中，指令字段包含以下值：MMX、SSE（1、2、3、3S）、EM64T、VT-x。输出告诉我，我可以运行托管代码，同时使用性能库调用将使用 SIMD 指令的函数。现在，让我们看看 SIMD 如何并行化数据。请看下图。
SIMD 机器利用了数据流的一个称为数据并行性的特性。在这种情况下，当您有一个需要执行相同指令的大量同类型数据时，就会产生数据并行性。数据并行性的一个典型例子是反转 RGB 图片以产生其负片。您需要遍历一个同类型整数值（像素）的数组，并对每个像素执行相同的操作（反转）——多个数据点，一个操作。现代的超标量 SISD 机器利用了指令流的一个称为指令级并行性（ILP）的特性。这意味着您在同一数据流上同时执行多个指令。因此，SIMD 机器与普通微处理器是不同类的机器。SIMD 关注利用数据流中的并行性，而超标量 SISD 关注利用指令流中的并行性。.NET Framework 4.0 idlertid不提供调用 SIMD 指令的机制。这意味着 TPL 中不包含任何调用 SIMD 指令的构造。
但是，我们仍然可以通过集成一个独立的库在 C# 代码中调用 SIMD 指令。Intel 提供了集成性能基元库 (IPP) 和 Intel 数学核心库。Intel MKL 是一个高度优化的数学例程库，用于科学、工程和金融应用。这些数学例程利用多线程和 SIMD 指令，根据底层硬件实现最佳性能。因此，让我们配置我们的开发环境以与此库集成。
使用 C# 源代码与 Intel 数学库进行互操作
读者须知：本文的这一部分不解释 Intel MKL 的构造——它仅演示了如何将 Intel 数学核心库 (Intel MKL) 函数调用并链接到 C# 程序。这提供了调用 SIMD 指令的函数。有兴趣的用户应下载 Intel C++ ComposerXE-2011 并将其安装到已安装的 Microsoft Visual Studio 2010 中。Intel C++ Composer 2011 附带 Intel C++ 编译器，以及一些性能库，如 Intel TBB、集成性能基元 (IPP) 和 Intel MKL。提供的源代码文件来自 Intel 软件网站。
Visual Studio 提供了一个命令提示符，其中已设置好运行 Microsoft Visual C++ 编译器的环境变量。如果您转到“开始”菜单并单击“Visual Studio”文件夹，您会看到一个名为“Visual Studio 工具”的子文件夹。启动命令提示符并输入“vcvars32.bat”。该命令的输出将验证您可以运行 cl.exe 编译器，更重要的是，还可以运行 nmake 工具。在提示符下，通过键入“md c:\Examples”创建一个新目录。将可下载的 C# 文件提取到此目录中。提供了调用 BLAS、LAPACK、DFTI（FFT 接口）、PARDISO 直接稀疏求解器和矢量数学库（VML）的示例。现在，因为这些是 C# 源代码文件，我们需要为 C# 编译器设置路径。
set PATH=%PATH%;.;C:\Windows\Microsoft.NET\Framework\v4.0.30319
此时如果启动 nmake 工具，会在构建过程中间遇到错误。因此，即使我们的文件位于 VC++ 命令行环境中，并且我们的环境变量已设置为使用 csc.exe C# 编译器驱动程序，我们也必须设置路径才能使用 Intel C++ 命令行编译器 icl.exe。Intel ComposerXE-2011 套件的典型安装会创建一个类似这样的目录路径：
C:\Program Files\Intel\ComposerXE-2011\bin\ia32
这里是 icl.exe 编译器的位置。但是，我们必须执行 iclvars.bat 命令，该命令位于 bin 目录中。我们还必须设置该环境变量并运行命令“iclvars.bat ia32”。请注意下图。
set PATH=%PATH%;.;C:\Program Files\Intel\ComposerXE-2011\bin
我们运行“iclvars.bat ia32”命令。之后，我们使用 nmake 运行构建脚本。
nmake ia32 MKLROOT="C:\Program Files\Intel\ComposerXE-2011\mkl"
这是执行构建脚本的输出。
文本已换行以避免页面滚动。
Add path of the MKL libraries to the lib environment variable
set lib=%MKLROOT%\ia32\lib;%lib%
MKL entries for custom dll
Workaround for pardiso
_fseeki64.c
Build MKL custom dll
nmake mkl.dll MACHINE=IX86
MKL_LIB="mkl_intel_c_dll.lib mkl_intel_thread_dll.lib
mkl_core_dll.lib" MSLIB=user32.lib
'mkl.dll' is up-to-date
Add path of the mkl.dll and path of the MKL binaries
to the path environment variable
set path=c:\codeExamples;%MKLROOT%\ia32\bin;%path%
Build and run examples
nmake /a dgemm.exe dgeev.exe dfti_d1.exe pardiso.exe vddiv.exe
Compile dgemm.cs
csc .\dgemm.cs
Microsoft (R) Visual C# 2010 Compiler version 4.0.30319.1
Copyright (C) Microsoft Corporation. All rights reserved.
warning CS1668: Invalid search path
'C:\Program Files\Intel\ComposerXE-2011\mkl\ia32\lib'
specified in 'LIB environment variable' --
'The system cannot find the path specified. '
Run dgemm example
dgemm.exe
MKL cblas_dgemm example
alpha=1
beta=-1
Matrix A
1 2 3
4 5 6
Matrix B
0 1 0 1
1 0 0 1
1 0 1 0
Initial C
5 1 3 3
11 4 6 9
Resulting C
0 0 0 0
0 0 0 0
TEST PASSED
Compile dgeev.cs
csc .\dgeev.cs
Microsoft (R) Visual C# 2010 Compiler version 4.0.30319.1
Copyright (C) Microsoft Corporation. All rights reserved.
warning CS1668: Invalid search path
'C:\Program Files\Intel\ComposerXE-2011\mkl\ia32\lib'
specified in 'LIB environment variable' --
'The system cannot find the path specified. '
Run dgeev example
dgeev.exe
MKL LAPACK dgeev example
Matrix A:
-2 -5 1 -1 -3
5 0 -4 -1 9
8 -1 1 6 0
2 0 3 1 -2
-4 -3 4 -4 -2
Matrix A on exit:
-3.74979687498421 3.30457083287809 0 0 0.272727272727273
-10.4969456995383 -3.74979687498421 0 0 0
6.28146648851724 6.68544105675216 -2.34478199802706 0 0
0.961154157305769 -2.58938020040906 -2.97051040398257
2.05286581472406 0
3.08511288713447 -0.862950290140918 -1.09233902416012
6.28154038161816 5.79150993327142
Wr on exit:
-3.74979687498421 -3.74979687498421 -2.34478199802706
2.05286581472406 5.79150993327142
Wi on exit:
5.88964350304832 -5.88964350304832 0 0 0
info on exit: 0
TEST PASSED
Compile dfti_d1.cs
csc .\dfti_d1.cs
Microsoft (R) Visual C# 2010 Compiler version 4.0.30319.1
Copyright (C) Microsoft Corporation. All rights reserved.
warning CS1668: Invalid search path
'C:\Program Files\Intel\ComposerXE-2011\mkl\ia32\lib' specified
in 'LIB environment variable' --
'The system cannot find the path specified. '
Run dfti_d1 example
dfti_d1.exe
Backward transform scale: 0.166666666666667
Initial data:
0 1 2 3 4 5
Resulting data:
0 1 2 3 4 5
TEST PASSED
Compile pardiso.cs
csc .\pardiso.cs
Microsoft (R) Visual C# 2010 Compiler version 4.0.30319.1
Copyright (C) Microsoft Corporation. All rights reserved.
warning CS1668: Invalid search path
'C:\Program Files\Intel\ComposerXE-2011\mkl\ia32\lib'
specified in 'LIB environment variable'
-- 'The system cannot find the path specified. '
Run pardiso example
pardiso.exe
=== PARDISO is running in In-Core mode, because iparam(60)=0 ===
================ PARDISO: solving a symmetric indef. system ================
The local (internal) PARDISO version is : 103000115
1-based array indexing is turned ON
PARDISO double precision computation is turned ON
METIS algorithm at reorder step is turned ON
Single-level factorization algorithm is turned ON
Scaling is turned ON
Summary PARDISO: ( reorder to reorder )
================
Times:
======
Time spent in calculations of symmetric matrix portrait(fulladj): 0.000014 s
Time spent in reordering of the initial matrix(reorder) : 0.000066 s
Time spent in symbolic factorization(symbfct) : 0.000239 s
Time spent in in allocation of internal data structures(malloc) : 0.000108 s
Time spent in additional calculations : 0.002201 s
Total time spent : 0.002628 s
Statistics:
===========
< Parallel Direct Factorization with #processors: > 2
< Numerical Factorization with Level-3 BLAS performance >
< Linear system Ax = b>
#equations: 8
#non-zeros in A: 18
non-zeros in A (): 28.125000
#right-hand sides: 1
< Factors L and U >
#columns for each panel: 128
#independent subgraphs: 0
< Preprocessing with state of the art partitioning metis>
#supernodes: 4
size of largest supernode: 4
number of nonzeros in L 31
number of nonzeros in U 1
number of nonzeros in L+U 32
Reordering completed ...
Number of nonzeros in factors = 32
Number of factorization MFLOPS = 0
Percentage of computed non-zeros for LL^T factorization
9 %
19 %
45 %
96 %
100 %
================ PARDISO: solving a symmetric indef. system ================
Summary PARDISO: ( factorize to factorize )
================
Times:
======
Time spent in copying matrix to internal data structure(A to LU): 0.000000 s
Time spent in factorization step(numfct) : 0.000126 s
Time spent in in allocation of internal data structures(malloc) : 0.000053 s
Time spent in additional calculations : 0.000001 s
Total time spent : 0.000179 s
Statistics:
===========
< Parallel Direct Factorization with #processors: > 2
< Numerical Factorization with Level-3 BLAS performance >
< Linear system Ax = b>
#equations: 8
#non-zeros in A: 18
non-zeros in A (): 28.125000
#right-hand sides: 1
< Factors L and U >
#columns for each panel: 128
#independent subgraphs: 0
< Preprocessing with state of the art partitioning metis>
#supernodes: 4
size of largest supernode: 4
number of nonzeros in L 31
number of nonzeros in U 1
number of nonzeros in L+U 32
gflop for the numerical factorization: 0.000000
gflop/s for the numerical factorization: 0.000541
Factorization completed ...
========= PARDISO: solving a symmetric indef. system =========
Summary PARDISO: ( solve to solve )
================
Times:
======
Time spent in direct solver at solve step (solve) : 0.000032 s
Time spent in additional calculations : 0.000068 s
Total time spent : 0.000100 s
Statistics:
===========
< Parallel Direct Factorization with #processors: > 2
< Numerical Factorization with Level-3 BLAS performance >
< Linear system Ax = b>
#equations: 8
#non-zeros in A: 18
non-zeros in A (): 28.125000
#right-hand sides: 1
< Factors L and U >
#columns for each panel: 128
#independent subgraphs: 0
< Preprocessing with state of the art partitioning metis>
#supernodes: 4
size of largest supernode: 4
number of nonzeros in L 31
number of nonzeros in U 1
number of nonzeros in L+U 32
gflop for the numerical factorization: 0.000000
gflop/s for the numerical factorization: 0.000541
Solve completed ...
The solution of the system is:
x [0] = -0.0418602012868094
x [1] = -0.0034131241592791
x [2] = 0.117250376805018
x [3] = -0.11263957992318
x [4] = 0.0241722444613714
x [5] = -0.107633340356223
x [6] = 0.198719673273585
x [7] = 0.190382963551205
TEST PASSED
Compile vddiv.cs
csc .\vddiv.cs
Microsoft (R) Visual C# 2010 Compiler version 4.0.30319.1
Copyright (C) Microsoft Corporation. All rights reserved.
warning CS1668: Invalid search path
'C:\Program Files\Intel\ComposerXE-2011\mkl\ia32\lib' specified
in 'LIB environment variable' -- 'The system cannot find the path specified. '
Run vddiv example
vddiv.exe
Vector A
1 2 3
Vector B
2 0 2
CallBack, in function: vdDiv errcode=2 index=1 arg1=2 arg2=0 res=Infinity
CallBack, result correction: 777
Vector Y
0.5 777 1.5
TEST PASSED
如果您执行“dir”DOS 命令，您将获得运行构建脚本处理 C# 文件后生成的文件的列表。
参考文献
- 《专业并行编程》作者：Gaston C. Hiller
- Intel 网站知识库：Intel 示例代码的来源
- Microsoft MSDN 库 - 编译器内在函数部分