英特尔多核平台编程优化大赛报告

代码优化前所需时间：4.765秒

代码优化后所需时间：0.25秒（保留小数点后7位精度）

前言

本次优化使用的CPU是Intel Xeon 5130 主频为2.0GHz 同Intel酷睿2一样是基于Core Microarchitecture 的双核处理器。本次优化在Intel的工具帮助下主要针对Core Microarchitecture 系列处理器进行优化。但是由于未知原因，Intel VTune Analyzers并不能在该系统下正常工作。所以，所有使用Intel VTune Analyzers的测试均使用另外一个奔腾D 820的系统测试。

第一章主要介绍了程序的串行优化。其中有关于Intel编译器使用，以及Intel Math Kernel Library使用，Intel VTune Analyzers使用的介绍。在借助Intel工具的帮助下，结合Intel Core Microarchitectured的特性。设计出了针对L1 Cache进行优化的，高效率的串行代码。程序的执行时间从优化前的4.765秒达到了优化后的0.765秒。

第二章主要介绍了程序的并行化。首先讨论了2种并行算法的优缺点。然后选择了适合本程序的并行算法进行优化。并且在最后分析了并行化时的性能瓶颈。通过并行化，程序达到了0.437秒。

第三章主要介绍了程序的汇编优化。首先介绍了计算的数学理论。然后介绍了汇编代码的编写。最后进行了性能分析。通过该步优化程序在保留小数点后7位精度的前提下达到了0.312秒的好成绩。并且在Intel酷睿2 E6600 上测试达到了0.25秒。

附录A 说明了本次报告的目录结构和优化方法。

附录B 列出了进行本次竞赛所参考的文献。

目录

一、串行优化

1.1 代码的基本修改和优化

1.2 基于Intel编译器的优化

1.3 使用Intel VTune Analyzers进行性能分析

1.3.1 Intel VTune Analyzers概述

1.3.2 基于SAMPLING方式的分析

1.3.3 对于本次程序的分析

1.4 优化computePot函数

1.5 使用Intel Math Kernel Library

1.6 根据Cache大小优化Intel Math Kernel Library调用

1.7 优化updatePositions函数

1.8 其他优化以及性能分析

二、并行优化

2.1 并行优化概述

2.2 优化方案一

2.3 优化方案二

2.4 并行实现

2.5 性能分析

三、汇编级优化

3.1 优化目标

3.2 数学理论

3.3 汇编码实现

3.4 性能分析

3.5 总结

附录A 目录结构和编译方法

附录B 参考文献

一、串行优化

1.1 代码的基本修改和优化

首先根据主办方的要求把代码的输出精度改为小数点后7位。

if (i%10 == 0) printf("%5d: Potential: %20.7f\n", i, pot);

在进行任何优化前代码的执行时间是4.765秒。

接着把项目转换成使用Intel C++ Compiler，代码的执行时间是4.531秒。

然后执行最基本的优化，把代码中的pow函数优化成乘法。代码如下：

distx = (r[0][j] - r[0][i])*(r[0][j] - r[0][i]);

disty = (r[1][j] - r[1][i])*(r[1][j] - r[1][i]);

distz = (r[2][j] - r[2][i])*(r[2][j] - r[2][i]);

执行时间依然为4.531秒。说明Intel编译器已经将pow函数优化掉了。

1.2 基于Intel编译器的优化

这里介绍本程序中基于Intel编译器优化技术。其中有些优化参数是可以确定的，有些优化参数需要在程序的不同阶段反复调试以确定最优方案，而有些优化技术是在后面的优化中使用的。

编译器优化级别

Intel的编译器共有如下一些主要的优化级别：

u /O1：实现最基本的优化

u /O2：基于代码速度实现常规优化，这个也是默认的优化级别

u /O3：在/O2的基础上实现进一步的优化，包括Cache预读，标量转换等等，但是在某些情况下反而会减慢代码的执行速度。

u /Ox：实现最大化的优化，包括自动内联函数的确定，全局优化，使用EBP作为通用寄存器等。

u /fast：等同于/O3, /Qipo, /Qprec-div-, and /QxP。

通过测试，目前选用/O3，但是随着代码的更改，需要重新测试，选择合适的优化级别。

针对特定处理器进行优化

Intel的编译器一共支持如下3种针对特定处理器的优化：

u /G：使用这个优化选项，Intel将针对特定的CPU进行优化，但是其代码依然可以在所有的CPU上执行。

u /Qx：使用这个优化选项，Intel将针对特定的CPU进行优化，并且产生的代码无法使用在不兼容的CPU上。

u /Qax：使用这个优化选项，Intel将针对特定的CPU进行优化，并且产生多份代码，在运行时根据CPU类型自动选择最优的代码。

由于本程序只需要运行在基于Core Microarchitecture 的处理器上，而无需考虑兼容性。所以本程序选择/Qx选项。并且针对运行时的酷睿2处理器，选择/QxT。但是在进行VTune测试时，由于测试平台为奔腾D 820，所以暂时使用/QxP的参数。

使用IPO

使用/Qipo可以启用Intel编译器的过程间优化(Interprocedural Optimizations)。通过过程间优化，编译器可以通过使用寄存器优化函数调用、内联函数展开、过程间常数传递、跨多文件优化等方式进一步优化程序。

此外，Intel编译器支持多文件的过程间优化，而由于本程序只有一个文件，所以并不需要使用。

但是IPO优化却会对本程序的调试带来极大的麻烦。所以本程序开发时不使用IPO优化，只有在最后的版本中才尝试使用IPO优化能否提高效率。

使用GPO

Intel编译器支持GPO(Profile-Guided Optimization)。GPO由一下三步组成。

第一步：使用/Qprof-gen编译程序，产生能记录运行细节的特殊程序。

第二步：运行第一步产生的程序，生成动态信息文件(.dyn)。

第三步，使用/Qprof-use，结合动态信息文件重新编译程序，产生更优化的程序。

通过使用GPO，Intel编译器可以更详细得了解程序的运行情况，从而根据实际情况产生更优化的代码。比如优化条件跳转，使得CPU分支预测的能力更准确，又如决定哪些函数需要内联，哪些不要内联等。

此外，基于GPO还有很多的工具方便用户开发程序。比如Code-Coverage Tool可以进行代码覆盖测试。

由于GPO收集的信息和特定的程序有关，而本程序一直在修改。所以本程序只在每个版本的最后部分使用GPO进行优化。

循环展开

循环展开(Loop Unrolling)通过在把循环语句中的内容展开从而使执行的代码速度更快。循环展开可以提高代码的并行程度，减少条件转移次数从而提高速度。另外，对于Pentium 4处理器，其分支预测功能可以精确得预测出16次迭代以内的循环，所以，如果能把循环展开到迭代次数在16次以内，对于特定的CPU可以提高分支预测准确度。

但是循环展开必须有一个度，并不是展开层数越多越好，展开层数多了，可能反而影响代码的执行速度。所以通常的做法是让编译器自己决定循环展开的层数。

Intel编译器对于循环展开有如下选项：

u /Qunrolln：执行循环展开n层。

u /Qunroll：让Intel编译器自己决定循环展开的层数。

此外Intel编译器还提供在了程序中使用编译制导语句规定某个特定循环的展开次数。如下例指示for循环展开n层。

#pragma unroll(n)

for(i=0;i<10000;i++){……}

所以本程序使用/Qunroll参数，让Intel编译器自己决定使用循环展开的层数。但是在程序的最终优化时，如果发现Intel编译器的循环展开并不是最优的，则通过在特定循环前加上编译制导语句，使用最佳的循环展开层数。

浮点计算优化

Intel编译器提供了很多基于浮点数的优化参数，有提供精度的，也有提高速度的。对于本程序，主要使用如下优化参数。

u /fp: fast或/fp: fast=1：这两个参数的等价的，同时也是默认的参数。他告诉编译器进行快速浮点计算优化。

u /fp: fast=2：这个参数比/fp: fast=1提供更高的优化级别，同时也可能带来更大的精度损失。

本程序使用/fp: fast=2优化，但是如果发生精度问题，可以考虑使用/fp: fast=1。

自动并行化

Intel的编译器支持自动并行化(Auto-parallelization)。通过/Qparallel可以打开编译器的自动并行化，编译器会在分析了用户的串行程序后，自动选择可以并行的部分进行并行化。自动并行化的有点是方便，不需要用户懂得专业知识，不需要更改原来的串行程序。但是缺点也是显而易见的，由于编译器并不知道用户的程序逻辑，所以无法很好得进行并行化。在对本程序试用/Qparallel后发现，效果并不好。所以本程序不只用/Qparallel进行自动并行化。

使用OpenMP并行化

OpenMP是一种通用的并行程序设计语言，其通过在源代码中添加编译制导语句，提示编译器如何进行程序的并行化。OpenMP具有书写方便，不需要改变源代码结构等多种优点。Intel的编译器支持OpenMP。本次程序并不打算使用OpenMP进行并行化，而打算使用Windows Thread。但是由于本程序需要使用到Intel Math Kernel Library，而Intel Math Kernel Library中的代码支持OpenMP并行化。所以有必要使用一些基本的OpenMP设置函数。

需要使用OpenMP，需要在编译时加上/Qopenmp选项。并且在源代码中包含” omp.h”文件。

OpenMP提供了函数omp_set_num_threads(nthreads)设置OpenMP使用的线程数，由于其设置会影响到Intel Math Kernel Library，所以将其设置成1，禁止Intel Math Kernel Library的自动并行化。

向量化

Intel的编译器支持向量化(Vectorization)。可以把循环计算部分使用MMX，SSE，SSE2，SSE3，SSSE3等指令进行向量化，从而大大提高计算速度。这也是本程序串行化时的主要优化点。前面提到的针对处理器的/QaxT优化选项已经打开了向量化。将代码向量化还有许多需要注意的地方，具体的注意点和方法将在后面具体的代码中说明。这里先给出一些对向量化有用的编译制导语句以及选项。

u /Qrestrict选项：当Intel编译器遇到循环中使用指针时，由于多个指针可能指向同一个地址，所以其无法保证指针指向内容的唯一性。故Intel编译器无法确定循环内数据是否存在依赖性。这是可以通过使用/Qrestrict选项与restrict关键字，指示某个指针指向内容的唯一性。从而能解决数据依赖性不确定的问题。

u #pragma vector编译制导语句：该编译制导语句一共包含3个。#pragma vector always用于指示编译器忽略其他因素，进行向量化。#pragma vector aligned用于指示编译器进行向量化时使用对齐的数据读写方式。#pragma vector unaligned用于指示编译器进行向量化时使用不对齐的数据读写方式。由于在使用SSE类指令进行向量化时，需要同时处理多个数据，所以每次读写的数据长度很长，可以达到128bit。所以将要处理的数据按照128bit(16byte)对齐，使用对齐的读写指令是可以提高程序运行速度的。但是需要注意的是对于实际没有对齐的数据使用#pragma vector aligned会造成程序运行错误。

使用变量对齐指示

Intel编译器提供了__declspec(align(n))用于在定义变量时指定其需要进行n字节对齐。变量对齐对于向量化计算的读取速度有很大关系。对于向量化计算一般使用__declspec(align(16))进行对齐。另外也可以使用__declspec(align(64))指定变量对齐到Cache的行首。关于Cache的行对齐的详细讨论请见后文的分析。

数据预读

通常数据是放在内存中，当要计算时才读入CPU进行计算。由于内存到CPU的传输需要很长时间，所以CPU中有多级Cache机制。Intel编译器支持数据预读优化选项。通过/Qprefetch打开数据预读优化，编译器会在使用数据前先插入预读指令，让CPU先把数据预读到Cache中，从而加快数据的访问速度。该选项默认情况下是打开的。此外Intel还提供了数据预读的编译制导语句，通过使用#pragma prefetch语句，用户可以人为得在程序中增加数据预读指令。但是需要注意的是，数据预读指令并不是越多越好的。不恰当的数据预读指令会占用内存带宽，把有用的数据从Cache中挤出去，反而影响速度。并且Core Microarchitecture体系结构已经支持给予硬件的数据预读指令。所以本程序倾向于使用给予硬件的数据预读机制。而由于/Qprefetch默认的打开的，也没有必要特意关闭该选项，Intel编译器有能力判断哪些地方可以通过合适的数据访问模式激活硬件数据预读机制，哪些地方需要额外添加数据预读指令。

产生调试信息

通过使用/Zi选项产生调试信息以帮助调试。默认为关闭。在本程序的开发阶段，打开此选项。在开发完成后关闭此选项。

使用全局优化

通过使用/Og选项打开编译器的全局优化功能。改选项需要在本程序不同的开发阶段分别尝试是否打开以确定最优优化选项。

针对Windows程序优化

通过使用/GA选项可以打开Intel编译器的针对Windows程序优化的功能。其实通过打开/GA选项，Intel可以提高访问Windows下thread-local storage(TLS)变量的速度。TLS变量通过__declspec(thread)来定义。在本程序中，并不打算使用TLS变量。但还是打开/GA选项。

内联函数扩展

Intel编译器可以通过/Obn来定义内联函数的扩展级别。当n为0禁止用户定义的内核函数的扩展。当n为1时，根据用户定义的inline关键字进行扩展。当n为2时，根据Intel编译器的自动判断进行扩展。本次程序使用/Ob2选项。

FTZ与DAZ

在计算机内浮点数是由尾数和指数组成的。尾数通常被规范化成[1,2)之间。但是当数字接近0时，由于其指数已经无法将尾数规范成[1,2)之间，所以需要在尾数表示成0.0000xx的形式。这种表示形式称为不规范的形式。其会影响CPU的浮点计算速度。并且由于这种数非常接近0，所有有时将其表示成0并不会影响计算的结果。所以CPU的浮点控制器有2个用于控制对于不规范数处理的选项。FTZ用于将计算结果中的不规范数表示成0，DAZ用于在读入不规范数时将其表示成0。Intel编译器提供了内置的宏来方便用户设置这两个模式。这两个宏分别是_MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON)和_MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON)。用户在程序中设置了这两个模式将有助于提高浮点计算速度。但是实际上对于本程序，由于已经使用了/O3以及SSE指令集优化。所以Intel编译器已经设置好了FTZ模式，用户不必另外设置FTZ。并且由于本程序中所有的数都是计算得来的，所以只要计算时使用了FTZ，那读取数据时就不会碰到不规范的数据，所以用户也没必要设置DAZ。

编译器报告

编译器报告虽然不能直接提供优化，但是却可以让用户了解编译器处理程序的信息，给用户更改源代码提供了很多有用的信息。对于本程序，向量化是非常重要的一步，而编译器报告可以指出某个地方是由于什么原因造成没有向量化。所以本使用使用/Qvec-report3参数对向量优化进行报告。

使用Intel编译器函数进行精确时间测量

Intel编译器提供了许多特殊的函数。这类函数一般都对应一条或者几条汇编语言。其可以让用户以比汇编语言方便的方式写出性能接近汇编语言的代码。其中最主要的是对SIMD类指令的支持。当然其中还有很多其他功能的函数。比如_rdtsc()函数。

需要注意的是要使用这些函数必需打开/Oi选项。这个选项默认是打开的。

当程序需要进行精确时间测量，比如优化后需要知道某段特定的代码到底快了多少毫米时，使用Windows的时间函数已经无法满足精度要求。这是用户可以使用Intel VTune Analyzers进行测量（具体使用方法将在后面介绍）。其实CPU已经提供了一个特殊的机器指令rdtsc，使用这条指令可以读出CPU自从启动以来的时钟周期数。由于现在的CPU主频已经是上GHz了。所以，其计时精度可以达到纳秒级。Intel提供的_rdtsc()函数使得用户不必再使用汇编语言，可以像调用函数一样得到CPU的时钟周期数。例子代码如下：

注：以下代码摘自“Intel C++ Compiler Documentation”

#include <stdio.h>

int main()

{

__int64 start, stop, elaspe;

int i;

int arr[10000];

start= _rdtsc(); 

for(i=0; i<10000; i++)

{

    arr[i]=i;

}

stop= _rdtsc();

elaspe = stop -start;

printf("Processor cycles\n %I64u\n", elaspe);

return 0;

}

优化结果

经过以上编译器选项的调整，程序的运行速度已经达到了2.25秒。

1.3 使用Intel VTune Analyzers进行性能分析

1.3.1 Intel VTune Analyzers概述

Intel VTune Analyzers用于监视程序或者系统的各种性能，从而为用户优化程序提供有价值的数据。同时Intel VTune Analyzers也能分析其收集的信息，给出用户优化程序的建议。Intel VTune Analyzers即支持本地的数据收集，也支持远程的数据收集。在本程序中，我们只需使用其本地数据收集功能。Intel VTune Analyzers共支持3种数据收集机制。每种机制都有其自己的适用范围，详细介绍如下：

u SAMPLING：其通过使用CPU内部的监视功能来检测系统底层的各种性能事件。使用这个功能无需在执行代码中插入特定的指令，因此其几乎没有探针效应。其无法给出函数间的调用关系。但是可以把相应的事件关联到程序中某行源代码或者汇编代码上。该方法通常适用于对某段程序的微调或者针对特定性能事件的调整上。

u CALL GRAPH：其通过在程序中插入特殊的指令，来记录每个函数执行的时间。函数间的调用关系等。其有一定的探针效应。该方法通常用于对于整个比较庞大的程序，进行分析，找出其中具有性能瓶颈的函数。

u COUNTER MONITOR：其无需在程序内部插入特殊的指令，因此其几乎没有探针效应。该方法即无法显示函数间的调用关系，也没法把事件定位到具体的某行代码中。该方式是用于测试整个系统的某些性能，比如CPU占用率，内存带宽等。通常用于系统级的调试。

对于本程序。由于程序结构简单。无需进行函数间调用的分析。而主要需要进行基于特定代码的分析。特别是后期需要针对CPU内部的事件特性进行源代码级甚至是汇编级的调试。所以本次优化主要采用SAMPLING方式。

1.3.2基于SAMPLING方式的分析

原理：Intel的CPU有一组性能检测寄存器，由于记录各种影响性能的事件。程序首先通过编程设定需要检测的事件，并且设定触发中断的计数值。当CPU中被检测的事件达到预设的值后触发相应的中断。Intel VTune Analyzers中的SAMPLING就是使用CPU的性能检测功能帮助用户分析程序的性能。其中有关于内存访问的事件，分支预测的事件，指令执行的事件等等。由于不同的CPU支持不同的性能事件，所以在不同的CPU上使用VTune时，所能监视的事件并不相同。

使用注意事项：SAMPLING一共支持2种统计。一种是Event，其是直接测量得到的值。另外一种是Event Ratio，其是基于多个Event计算得到的，有时更有实际意义，更直观。需要注意的是，每个Event都有一个预设的值，当这个预设的值到了以后，CPU引起中断，VTune进行统计。而这个值的设置不能太大，否则统计到的事件不够多，无法分析。也不能太小，否则频繁引起中断，会加大探针效应。用户可以在每个Event上手工设置合适的Sample After值，也可以通过选项卡上的选项，让VTune先运行一遍程序，然后根据实际的事件数量来校准触发值。对于本程序，这点尤其需要引起注意。因为本程序优化到后面时间非常短，如果不校准触发值，分析的效果会不理想。需要注意的是Clockticks和Instructions Retired这两个最基本的事件，默认是不校准触发值的，我们需要把他们调整成自动校准。此外对于某个Event的发生，大部分的中断点并不是精确的。即真正发生该事件的指令在所记录事件指令的前几条。但是有一部分属于精确事件，引起这类事件的指令正好是发生中断的前一条。

1.3.3对于本次程序的分析

本程序首先使用VTune最基本的3个事件(Clockticks、Instructions Retired和CPI)进行程序耗时分析。其结果如图：

说明程序中耗时最长的是computePot函数。

1.4 优化computePot函数

在对computePot函数向量化前，我们可以注意到distx，disty，distz三个变量都是临时变量。先将这3个变量去掉，从而可以使得Intel编译器能够更灵活得进行中间结果优化。另外最完成循环的i虽然是从0开始的，但是实际0和1并不进行计算，所以把外层循环的i设置层从2开始。代码如下：

   for( i=2; i<NPARTS; i++ ) {

      for( j=0; j<i-1; j++ ) {

        dist = sqrt( (r[0][j] - r[0][i])*(r[0][j] - r[0][i]) + (r[1][j] - r[1][i])*(r[1][j] - r[1][i]) + (r[2][j] - r[2][i])*(r[2][j] - r[2][i]) );

        pot += 1.0 / dist;

      }

   }

此时编译器显示内层循环已经向量化了。但是这个绝非我们的目标。为了提高计算开根号倒数的速度，为了使用Intel Math Kernel Library，我们需要把开根号倒数的计算先存在一组向量中，再一同计算。既将dist变量变成，dist数组，然后再对dist数组统一计算，再求和。代码如下：

   for( i=2; i<NPARTS; i++ ) {

      for( j=0; j<i-1; j++ ) {

           dist[j] = (r[0][j] - r[0][i])*(r[0][j] - r[0][i]) + (r[1][j] - r[1][i])*(r[1][j] - r[1][i]) + (r[2][j] - r[2][i])*(r[2][j] - r[2][i]);

      }

      for( j=0; j<i-1; j++ ) {

              dist[j] = 1.0 / sqrt(dist[j]);

      }

      for( j=0; j<i-1; j++ ) {

              pot += dist[j];

      }

   }

Intel编译器提示，内部的3个循环都进行了向量化。此时出现了令人惊喜的成绩。程序的执行时间突然降到了1.453秒。使用VTune进行分析，发现Intel编译器对于开根号倒数的计算自动调用了内部的向量化代码库。注意此时，还没有使用Intel Math Kernel Library，所以这个向量代码库是Intel编译器内置的，虽然效率没有使用Intel Math Kernel Library高，但是速度已经提高了很多。调用Intel编译器内置的向量库的结果如图：

1.5 使用Intel Math Kernel Library

Intel Math Kernel Library中提供了一部分的向量函数(Vector Mathematical Functions)。这类函数提供了对于普通数学计算函数的快速的向量化计算。VML中有一个向量函数就是计算开根号倒数的。

Intel的VML库中提供了如下函数来计算整个向量中各个数的开根号倒数：

vdInvSqrt( n, a, y )

其中n表示计算的元素个数。a是指向输入计算数据数组的头指针。y是指向输出计算数据数组的头指针。其中a和y可以相同。

要使用该函数，首先需要在头文件中包含”mkl.h”，并且链接mkl_c.lib文件和libguide40.lib文件。

除了基本计算功能外，VML还提供了一个设置模式的函数，用于设置特定的计算模式：

vmlSetMode ( mode )

其中的mode是一个预定义宏。在我们的程序中，需要设置如下模式：

VML_LA：VML的所有向量函数都提供了2个精度的版本。精度低的版本计算速度也相对比较快。本程序只需要保留小数点后7位精度。低精度的版本符合要求，所以设定VML使用低精度的版本。

VML_DOUBLE_CONSISTENT：该选项用于控制FPU的计算精度为double，其实由于我们这次使用的函数基本上是使用SSE2指令集进行计算的，和FPU没什么关系。但是也可能存在使用FPU的可能，所以设定VML使FPU的精度为double。

VML_ERRMODE_IGNORE：该选项用于关闭VML的错误处理功能，本程序不需要进行错误处理。

VML_NUM_THREADS_OMP_FIXED：VML函数都能使用OpenMP，根据特定的硬件环境进行并行化。而我们并不需要其进行并行化。所以使用该选项和前面提到的omp_set_num_threads(1)结合。关闭VML的自动并行化功能。

具体的代码如下：

   for( i=2; i<NPARTS; i++ ) {

      for( j=0; j<i-1; j++ ) {

           dist[j] = (r[0][j] - r[0][i])*(r[0][j] - r[0][i]) + (r[1][j] - r[1][i])*(r[1][j] - r[1][i]) + (r[2][j] - r[2][i])*(r[2][j] - r[2][i]);

      }

        vdInvSqrt(i-1,dist,dist);

      for( j=0; j<i-1; j++ ) {

              pot += dist[j];

      }

   }

优化后出现了令人可惜可贺的成绩：0.796秒。

1.6 根据Cache大小优化Intel Math Kernel Library调用

在上面的程序中对于MKL函数的调用是每次内部循环都执行一次调用，我们知道每次执行函数的调用都是需要开销的，那是否有更优化的调用MKL方法那？下面这段话摘自Intel Math Kernel Library的说明文档上：

There are two extreme cases: so-called "short" and "long" vectors (logarithmic scale is used to show both cases). For short vectors there are cycle organization and initialization overheads. The cost of such overheads is amortized with increasing vector length, and for vectors longer than a few dozens of elements the performance remains quite flat until the L2 cache size is exceeded with the length of the vector.

下面这副性能分析图片摘自Intel Math Kernel Library的网站上：

从这段文字和这副图片中，我们了解到对于MKL函数的调用时，所处理的向量不能太短，否则函数的建立时间开销将是非常大的，也不能太长，操作了L2 Cache，否则函数执行时访问内存的开销是很大的。并且通过图片了解到不合适的长度对于函数的性能将产生指数级影响。

根据理论计算：每次执行computePot函数，总共需要执行的计算量为(1+998)*998/2=498501个。每个double类型占用8个字节，所有总共需要占用的空间为498501*8=3988008byte=3894KB。而这次进行竞赛的测试平台的CPU的L2 Cache大小为2M，由于有2个线程同时计算，平均每个线程分到的L2 Cache为1M。由于L2 Cache可能还被其他数据占据。所以为了保证所计算的数据在L2 Cache中，最好每次计算的向量长度在512KB左右。故把整个computePot函数的计算量分成8份。每份计算量的中间结果向量长度为3894KB/8=486KB。

但是实际情况并非如此，进行这种优化后，程序的执行速度反而降低了。通过分析发现原来CPU中的L1 Cache大小为32KB。数组r有3000个元素，如果每次迭代都进行vdInvSqrt调用。那dist的长度为1000个元素左右。加起来正好可以全部在L1 Cache中。而如果合并起来调用vdInvSqrt，则由于vdInvSqrt过长。其L1 Cache中存放不下，需要存放在L2 Cache中，从而反而影响了速度。看来，对于本程序，不应该根据L2 Cache进行优化，而应该根据L1 Cache进行优化。但是对于只有几个或者几十个数据就调用MKL函数，其开销还是很大的。因此本程序使用了折中的方法，对于前面非常小的几十个数据，凑足1000个放在一起进行计算，而后面的数据还是按照原来的方式计算。具体实现的代码如下：

for( i=2,k=0; i<47; i++ ) {

for( j=0; j<i-1; j++,k++ ) {

         dist[k] = (r[0][j] - r[0][i])*(r[0][j] - r[0][i]) + (r[1][j] - r[1][i])*(r[1][j] - r[1][i]) + (r[2][j] - r[2][i])*(r[2][j] - r[2][i]);

      }

   }

   vdInvSqrt(k,dist,dist);

for( j=0; j<k; j++ ) {

      pot += dist[j];

   }

for( i=47; i<NPARTS; i++ ) {

for( j=0; j<i-1; j++ ) {

         dist[j] = (r[0][j] - r[0][i])*(r[0][j] - r[0][i]) + (r[1][j] - r[1][i])*(r[1][j] - r[1][i]) + (r[2][j] - r[2][i])*(r[2][j] - r[2][i]);

      }

      vdInvSqrt(i-1,dist,dist);

for( j=0; j<i-1; j++ ) {

         pot += dist[j];

      }

   }

通过该不优化，程序的性能略微有所提高，达到了0.781秒。

1.7 优化updatePositions函数

虽然updatePositions函数执行的时间非常短。但还是值得优化的。

首先进行的是基于数学的优化。我们发现在updatePositions和initPositions中，都有加0.5的计算。但是从后面的computePot的相减计算中发现，这个0.5是被抵消的，既不加0.5对结果没有影响。故去掉该加0.5的计算。另外updatePositions和initPositions中都有除以RAND_MAX的计算。而通过提取公因子的变换发现，如果此处不除以RAND_MAX而将最后的pot乘以RAND_MAX，则最后结果相同。故去掉该处的除以RAND_MAX的计算，而以在pot上一次乘以RAND_MAX为替换。具体代码如下：

void initPositions() {

int i, j;

for( i=0; i<DIMS; i++ )

for( j=0; j<NPARTS; j++ )

         r[i][j] = (double) rand();

}

void updatePositions() {

int i, j;

for( i=0; i<DIMS; i++ )

for( j=0; j<NPARTS; j++ )

         r[i][j] -= (double) rand();

}

在main函数中：

      pot = 0.0;

      computePot();

      pot*=(double)RAND_MAX;

if (i%10 == 0) printf("%5d: Potential: %20.7f\n", i, pot);

其次需要进行updatePositions内rand函数的优化。虽然rand函数本身的执行时间非常短，但是其频繁得进行调用却影响了性能。通过查找Microsoft Visual Studio .NET 2005中提供的源代码。将其中的rand函数提取出来，进行必要的修改，并且加上inline属性。从而加快程序的调用速度。具体代码如下：

int holdrand=1;

inline int myrand (){

return( ((holdrand = holdrand * 214013L+ 2531011L) >> 16) & 0x7fff );

}

经过上述优化，代码的执行速度已经达到了0.765秒。

1.8 其他优化以及性能分析

至此，该程序串行优化部分已经一本完成。但是还有一点细小的地方需要优化。

变量对齐对于数据读取速度是非常重要的。尤其是使用SIMD指令集进行优化后，对于对齐的变量，可以使用对齐的读写指令提高速度。一般对于SIMD指令需要进行16字节对齐。但是对于本程序，由于后面要进行多线程优化，而多线程执行时基于Cache Line的共享冲突会对读写造成很大的损失。故本程序使用64字节对齐。代码如下：

__declspec(align(64)) int holdrand=1;

__declspec(align(64)) double r[DIMS][NPARTS];

__declspec(align(64)) double pot;

__declspec(align(64)) double dist[1048];

在computePot函数的第一次迭代中。有一处进行pot累加的地方，使用了k变量作为循环条件。但是其实该变量的确切值是可以计算出来的。通过计算出该变量的确切值，可以让Intel编译器在编译时就知道循环的次数，从而有助于优化。具体代码如下（注意1035这个值）：

for( i=2,k=0; i<47; i++ ) {

for( j=0; j<i-1; j++,k++ ) {

         dist[k] = (r[0][j] - r[0][i])*(r[0][j] - r[0][i]) + (r[1][j] - r[1][i])*(r[1][j] - r[1][i]) + (r[2][j] - r[2][i])*(r[2][j] - r[2][i]);

      }

   }

   vdInvSqrt(k,dist,dist);

for( j=0; j<1035; j++ ) {

      pot += dist[j];

   }

此外再调整以下编译器的某些优化参数，选择合适的使用。比如使用哪个编译级别，是否打开全局优化，使用IPO，使用GPO等。

至此本程序的串行优化全部完成。使用Intel VTune Analyzers的分析结果为：

Full Name	CPI	Clockticks events	Clockticks %
void updatePositions(void)	3.214080375	8274621	0.287907869
int computePot(void)	1.294881302	926757552	32.24568138
mkl_vml_core_t7_vml_dInvSqrt_50	0.91981472	1925228486	66.9865643

（注：此分析数据是在奔腾D 820上测得）

从以上数据上表明updatePositions函数说执行的事件非常短，低于1%，computePot函数的执行时间在三分之一左右。mkl_vml_core_t7_vml_dInvSqrt_50的执行时间在三分之二左右。这些数据对下面一步并行化采用的策略是非常重要的。

二、并行优化

2.1 并行优化概述

在进行本程序的并行优化前先谈谈并行优化需要注意的问题。在并行优化中经常用到数据重复和计算重复的方法。所谓数据重复，就是为了保证多个线程能同时进行计算，就把数据复制多份来提高并行度。所谓计算重复，就是有时使用计算换通信的方法，提高并行度。

在对本程序进行优化前需要注意的是。测试平台使用的是基于Core Microarchitecture结构的。这个结构的双核CPU是共享L2 Cache的。但是当数据在一个核中进行修改，另外一个核去读他时，需要消耗几十个时钟周期的延迟。其代价的非常高的。这里需要注意的是，数据在Cache中是按行进行存放的，也就是说，CPU看待数据有没有被修改过是根据Cache Line的。所以2个分别被不同的核修改的数据如果存在于同一行Cache中，访问时的效率就会非常低。也就是发生了共享冲突。所以在分配变量时要尽量把不同性质的变量分配到不同的Cache Line中。我们的测试平台的L1 Cache和L2 Cache都是每行64byte的。所以前一章中的变量对齐都使用了64byte对齐。同样，在程序并行化时也需要考虑这种情况。

2.2 优化方案一

此方案使用数据重复的方法。程序可以定义2个r数组。以及2个pot数组。通过定义2个r数组，使得主线程可以在从线程使用一个r数组计算时同时更新第二个r数组。即主线程先更新r数组，然后主线程和从线程同时开始计算。但是从线程的计算量比主线程大一点。这样当主线程计算完后，可以继续更新第二个r数组，而此时从线程还在计算原来r数组的内容。当主线程更新完第二个r数组时，从线程正好完成前面的计算，并和主线程一同计算第二个r数组，依次类推。同时2个pot数组，一个给主线程计算每步的中间结果，另一个给从线程计算每步的中间结果。等计算结束后，再将其结果相加，打印。

优点：使用该方法的优点是显而易见的，理论上线程可以做到完全同步。

缺点：使用该方法的缺点是，从线程每次计算需要从主线程计算好的r数组中读取内容，由于是2个核，所以其访问延迟非常大。此外使用2个数值，每次迭代都需要将指针指向使用的数组，增加了程序的设计难度。同时计算任务分配的调优也是非常繁琐的。

由于在前一章中，我们发现updatePositions函数所花费的时间非常短。所以做到线程间的完全平衡意义并不大。

2.3 优化方案二

在前一个方案中，我们提到了线程的完全平衡的算法。同时我们发现完全平衡的意义不大。因此我们设计适合本程序的更优的方案。既然updatePositions函数所花费的时间非常短。那2个线程同时执行updatePositions造成的额外开销也是可以忽略的。本方案使用了数据重复和计算重复的方法。同样使用2个r数组，但是2个线程同时进行重复计算，并且2个线程分区完成不同的迭代步骤的computePot计算。即主线程完成整个r数组的更新，但是只计算其中的奇数次迭代。从线程同样完成整个r数组的更新，但是只进行偶数次迭代。并且同样使用了一个pot数组，2个线程分别将自己的计算结果先存储到pot数组中。等最后同步的时候再打印。

优点：使用该方案，程序的设计相对来说比较简单，负载均衡的调整也很容易。程序只需要很少的同步操作（在本程序中，只使用了2次同步）。并且重要的是。由于2个线程都在各自的CPU上使用各自的数据进行计算，所以最大化得避免了共享冲突的发生。同时也保留了前一章优化中针对L1 Cache的命中率。

缺点：该方案的缺点是存在重复计算。但是通过前面VTune的测试，已经发现其重复计算量非常小，可以忽略。

2.4 并行实现

本程序使用方案二进行并行化。首先将所有需要计算的数据和函数都复制2份，代码如下：

int computePot1(void);

void initPositions1(void);

void updatePositions1(void);

int computePot2(void);

void initPositions2(void);

void updatePositions2(void);

__declspec(align(64)) int holdrand1=1;

__declspec(align(64)) double r1[DIMS][NPARTS];

__declspec(align(64)) double pot1;

__declspec(align(64)) double dist1[1048];

__declspec(align(64)) int holdrand2=1;

__declspec(align(64)) double r2[DIMS][NPARTS];

__declspec(align(64)) double pot2;

__declspec(align(64)) double dist2[1048];

__declspec(align(64)) double potfinal[264];

其中的potfinal数组记录每次迭代的计算结果，用于最后的数组。

在主函数的并行中。我们发现由于偶数次迭代比奇数次迭代需要多算一次。故本程序的偶数次迭代在进行到快完成前先释放一个同步锁。使得主线程可以先输出一部分数据。而从线程在执行完所有的偶数次迭代后再释放一个同步锁，使主线程输出剩余的数据。由于输出数据也有一点的耗时。所以使用这种方法可以提高一点并行度。另外在本代码中使用了SetThreadAffinityMask分别设置不同的线程对应各自的CPU，以防止线程在不同的CPU中切换从而影响L1 Cache命中率。具体代码如下：

DWORD WINAPI mythread( void *myarg ){

int i;

     SetThreadAffinityMask(GetCurrentThread(), 2);

     initPositions2();

     updatePositions2();

for(i=0;i<=190;i+=2){

        pot2 = 0.0;

        computePot2();

        pot2*=(double)RAND_MAX;

        potfinal[i]=pot2;

        updatePositions2();

        updatePositions2();

     }

     ReleaseSemaphore(semmiddle, 1, NULL);

for(i=192;i<=NITER;i+=2){

        pot2 = 0.0;

        computePot2();

        pot2*=(double)RAND_MAX;

        potfinal[i]=pot2;

        updatePositions2();

       updatePositions2();

     }

     ReleaseSemaphore(semafter, 1, NULL);

return 0;

}//从线程

int main() {

int i;

int myarg=0;

   clock_t start, stop;

   omp_set_num_threads(1);

   vmlSetMode(VML_LA);

   vmlSetMode(VML_DOUBLE_CONSISTENT);

   vmlSetMode(VML_ERRMODE_IGNORE);

   vmlSetMode(VML_NUM_THREADS_OMP_FIXED);

   semmiddle = CreateSemaphore(NULL, 0, 1, NULL);

   semafter = CreateSemaphore(NULL, 0, 1, NULL);

   CreateThread(0, 8*1024, mythread, (void *)&myarg, 0, NULL);

   SetThreadAffinityMask(GetCurrentThread(), 1);

   initPositions1();

   start=clock();

for(i=1;i<NITER;i+=2){

        pot1 = 0.0;

        updatePositions1();

        updatePositions1();

        computePot1();

        pot1*=(double)RAND_MAX;

        potfinal[i]=pot1;

   }

   WaitForSingleObject(semmiddle, INFINITE);

for(i=0;i<=190;i+=10)

        printf("%5d: Potential: %20.7f\n", i, potfinal[i]);

   WaitForSingleObject(semafter , INFINITE);

   i=200;

   printf("%5d: Potential: %20.7f\n", i, potfinal[i]);

   stop=clock();

   printf ("Seconds = %10.9f\n",(double)(stop-start)/ CLOCKS_PER_SEC);

}//主线程

2.5 性能分析

并行化后的性能并不没有像理论中这么高只有0.437秒。于是我们开始查找原因。通过使用Intel Threading Checker我们发现，VML库中存在着访问冲突。图片如下：

当然这个错误有可能是Intel Threading Checker的误报。因为程序每次执行都没有发现不正确的结果，并且VML函数的文档上说明是线程安全性的。

由于兼容性原因，本系统无法使用Intel VTune Analyzers进行每个函数的耗时分析。于是使用Intel编译器提供的内置函数_rdtsc()记录不同部分所花费的CPU时钟周期。结果发现VML函数的总执行时间大概增加了0.088秒左右。说明VML函数在用户使用Windows Thread函数并行化访问时，其同步开销可能有一定的影响。

三、汇编级优化

3.1 优化目标

本程序主要的执行时间在computePot函数与VML库中。对于computePot函数，通过查看Intel编译器产生的汇编码发现其已经很优了。而对于VML函数由于其需要满足通用性，所以本程序应该可以设计出最适合本程序的计算函数来。

3.2 数学理论

Intel的CPU支持的SSE2指令中，有2条是用于计算双精度浮点的开根号倒数的。sqrtpd指令可以同时计算2个double型的开根号，其吞吐率为28个时钟周期。divpd指令用于计算2个数的除法，即用于计算倒数，其吞出率为17个时钟周期。由此可以计算出，如果当当使用这2条指令计算双精度数的开根号倒数，那即使使用汇编语言，忽略其他开销。计算每个元素的时钟周期也有(17+28)/2=22.5。而Intel的VML库计算每个元素的只需要10多个时钟周期，说明其肯定是通过其他快速的数学计算方法得到的。所以要优化vdInvSqrt函数，关键是找到更快速的数学计算方法。在Quake 3在源代码中有如下一段具有传奇色彩的代码：

float InvSqrt(float x){

       float xhalf = 0.5f*x;

       int i = *(int*)&x; // get bits for floating value

       i = 0x5f3759df - (i>>1); // gives initial guess y0

       x = *(float*)&i; // convert bits back to float

       x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy

       return x;

}

（注：以上代码的注释摘自CHRIS LOMONT的《FAST INVERSE SQUARE ROOT》文章中）

在上面的代码中最后一条是典型的牛顿迭代，可以根据精度要求进行多次迭代。这段代码神奇的地方在于初始值的估算上，只用了减法和移位2个简单的操作，达到了非常接近的估算值。我们称0x5f3759df为幻数(magic number)。CHRIS LOMONT在他的《FAST INVERSE SQUARE ROOT》文章中给出了对于这个幻数的解释和计算方法。并且计算出了理论上最优的适用于double类型的幻数为0x5fe6ec85e7de30da。说们我们的代码中可以使用该方法进行计算，示例代码如下：

double myinvsqrt (double x)

{

double xhalf = 0.5*x;

__int64 i = *(__int64*)&x;

     i = 0x5fe6ec85e7de30da - (i>>1);

     x = *(double*)&i;

     x = x*(1.5-xhalf*x*x);

     x = x*(1.5-xhalf*x*x);

     x = x*(1.5-xhalf*x*x);

     x = x*(1.5-xhalf*x*x);

return x;

}

但是不幸的是，根据调试，需要达到比赛要求的小数点后7位精度，必需进行4此牛顿迭代也行。而4此牛顿迭代的计算量使得这个方法对于Intel的VML函数来说毫无优势可言。那能否降低牛顿迭代的次数那？

我们发现如果以上代码只进行3次牛顿迭代，那误差只有小数点最后的1，2位。CHRIS LOMONT在他的文中提到他说计算出来的理论最优值，而这个幻数只是在线性估计时是最优的。在多次牛顿迭代中，这个值并不是最优的。CHRIS LOMONT并没有给出对于多次牛顿迭代最优幻数的计算方法，他在文章中对于float类型的实际最优值也是穷举得到的。我们同样在理论最优值0x5fe6ec85e7de30da的基础上进行了一定的穷举操作，发现的确有更优的幻数。但是即使使用了更优的幻数，还是无法在3次牛顿迭代基础上达到精度要求。但是我们发现所有的数值都偏小。于是我们可以在三次牛顿迭代后再乘一个比1大一点点的偏移量。从而能做到3次牛顿迭代就能达到精度要求。示例代码如下：

double myinvsqrt (double x)

{

double xhalf = 0.5*x;

__int64 i = *(__int64*)&x;

     i = newmagicnum - (i>>1);

     x = *(double*)&i;

     x = x*(1.5-xhalf*x*x);

     x = x*(1.5-xhalf*x*x);

     x = x*(1.5-xhalf*x*x);

     x = x*offset

return x;

}

由于时间原因，这里并没有对newmagicnum和offset进行详细的计算与统计。只给出一个对于本程序相对较优的newmagicnum值0x5fe6d250b0000000。

在上面的代码中只进行了3次牛顿迭代。对于Intel的VML来说也没有什么优势可言。那能不能再减少一次牛顿迭代，只进行2次迭代就达到精度要求那？

我们知道要进行2次牛顿迭代就达到精度要求就必须对其初始值的估计更加准确。而使用上面的方法估计的初始值已经无法满足该准确性。这是通过查找《Intel 64 and IA-32 Architectures Optimization Reference Manual》，我们发现SSE指令集中有一条RSQRTPS的指令用于同时计算四个单精度浮点数的开根号倒数，而其在Core Microarchitecture上的延迟为3个周期，吞吐率为2个周期。也就是说我们可以在极短的时间内就算出单精度类型的开根号倒数值（看来在现在的CPU上，当初Quake 3那段具有传奇色彩的代码已经没有用了）。于是我们想到了先使用单精度类型精度初值估算，然后再使用牛顿迭代。实验结果表明该方法只需要进行2次牛顿迭代就能满足小数点后7位的精度要求。示例代码如下：

double myinvsqrt (double x)

{

double xhalf = 0.5*x;

float xf=(float)x;

__asm{

         movss xmm1,xf;

         rsqrtss xmm1,xmm1;

         movss xf,xmm1;

     }

     x=(double)xf;

     x = x*(1.5-xhalf*x*x);

     x = x*(1.5-xhalf*x*x);

return x;

}

不幸的是由于该代码涉及到了复杂的算法以及类型转换，Intel的编译器并无法将其很好的并行化。所以只有依靠手工使用汇编语言将其优化。

3.3 汇编码实现

在实现汇编码前先要将原来的代码进行优化，将牛顿迭代中的减法变成加法，代码如下：

double myinvsqrt (double x)

{

double xhalf = -0.5*x;

float xf=(float)x;

__asm{

         movss xmm1,xf;

         rsqrtss xmm1,xmm1;

         movss xf,xmm1;

     }

     x=(double)xf;

     x = x*(1.5+xhalf*x*x);

     x = x*(1.5+xhalf*x*x);

return x;

}

进行这种转变是一点都不影响计算结果的。但是确可以提高计算速度。这是因为，如果执行的是减法，汇编语言的减法指令会将结果存在原来存放被减数（即1.5）的寄存器中。从而覆盖掉了原来的常数1.5，使得每次计算必须重新读入该参数。而优化成加法后则没有这个问题。

下面列出了本次汇编语言优化时使用的主要的汇编指令及其延迟，吞吐率和使用的计算部件。这些数据对优化汇编代码有帮助。

指令名	延迟	吞吐率	计算部件
movapd	1	0.33	FP_MOVE
cvtpd2ps	4	1	FP_ADD,MMX_SHFT
cvtps2pd	2	2	FP_ADD,MMX_SHFT,MMX_ALU
shufps	2	1	MMX_SHFT
rsqrtps	3	2	MMX_MISC
mulpd	5	1	FP_MUL
addpd	3	1	FP_ADD

（注：以上数据摘自《Intel 64 and IA-32 Architectures Optimization Reference Manual》）

在进行优化前，还有一点需要注意的是。rsqrtps函数是4个元素一算的，所以本程序使用4个元素作为一次计算单元来向量化。而用户输入的数据并不可能是正好4个元素。对于Intel编译器以及VML函数库来所，其使用的解决方法称为” Strip-mining and Cleanup”。即先按照4个数据一组进行计算。对于剩下的个别数据再进行单独计算。这对于通用化的程序来说是必须的。但是在我们的程序中，多计算几个并不会影响结果。而对于单独几个的数据如果另外处理不但会增加程序设计的复杂性，而且性能也可能会降低。所以本程序使用过渡计算的方法。即对于需要计算的数据中不足4个的，补满4个将其后面的数据计算掉。但是此时需要注意，由于dist变量是全局变量，默认的值为全0。如果过渡计算遇到0的值，速度可能会受到影响。所以本程序需要在一开始，将会被过渡计算使用到，但是从来不会被初始化的存储单元，初始化成1。具体代码如下：

void myinvsqrt (double *start,double *end)

{

__asm{

         mov esi,start;

         mov edi,end;

         test edi,0x0000001f;

         jz myalign;

         and edi,0xffffffe0;

         add edi,32;

myalign:

myagain:

         movapd xmm0,[esi];

         movapd xmm3,[esi+16];

         cvtpd2ps xmm6,xmm0;

         cvtpd2ps xmm7,xmm3;

         shufps xmm6,xmm7,01000100b;

         rsqrtps xmm6,xmm6;

         cvtps2pd xmm1,xmm6;

         shufps xmm6,xmm6,01001110b;

         cvtps2pd xmm4,xmm6;

         mulpd xmm0,mulcc;

         mulpd xmm3,mulcc;

         movapd xmm2,xmm1;

         movapd xmm5,xmm4;

         mulpd xmm1,xmm1;

         mulpd xmm4,xmm4;

         mulpd xmm1,xmm0;

         mulpd xmm4,xmm3;

         addpd xmm1,addcc;

         addpd xmm4,addcc;

         mulpd xmm1,xmm2;

         mulpd xmm4,xmm5;//前半段

         movapd xmm2,xmm1;

         movapd xmm5,xmm4;

         mulpd xmm1,xmm1;

         mulpd xmm4,xmm4;

         mulpd xmm1,xmm0;

         mulpd xmm4,xmm3;

         addpd xmm1,addcc;

         addpd xmm4,addcc;

         mulpd xmm1,xmm2;

         mulpd xmm4,xmm5;

         movapd [esi],xmm1;

         movapd [esi+16],xmm4;

         add esi,32;

         cmp esi,edi;

         jne myagain;

     }

}

//后半段

myinvsqrt(dist1,dist1+k); //调用方法

对于本函数的调用方法为分别传入其需要计算数据的头指针和尾指针。

3.4 性能分析

使用汇编语言优化后，程序跑出了惊人的0.312秒的好成绩。并且所有的输出数据全部都满足小数点后7位的精度要求。在使用Intel Threading Checker和Intel Threading Profiler分析程序时也得到了相对比较好的结果。如下图：

在Intel Threading Checker的检测中，没有发现程序有任何冲突。在使用Intel Threading Profiler的分析中，表现出了程序良好的并行性。

最后，在另外一台Intel酷睿2 E6600的机器上测试时，程序达到了0.25秒的好成绩，并且所有数据输出精度都达到了小数点后7位。

3.5 总结

在本次优化比赛中。我花了几个星期仔细钻研Intel的工具使用方法，并且结合Intel的CPU特性对源代码进行优化。在经过了漫长的调优后，终于在保留小数点后7位精度的七位精度的情况下达到了0.25秒的成绩。这里需要说明的是，在本程序的最后优化时虽然没有使用Intel的VML库，但并不是意味着VML库不好。VML库的通用化和高效率是有目共睹的。而是由于VML库是通用库，其需要考虑很多情况，而针对本程序自己设计的计算函数却不用考虑各自情况。所以设计有针对性的函数才能提高速度。当然要设计这种函数对用户的要求太高，需要了解数学理论，汇编语言，以及优化的方法。所以对于一般的用户还是使用VML库比较好。

最后需要说明的是，由于本次竞赛的时间有限。很多地方都没有得到更好的优化。比如那段汇编语言就能针对Intel CPU的指令延迟特性进一步优化。如果大赛能给出更多的时间，那有可能可以优化得更好。

附录A 目录结构和编译方法

本报告压缩文件的根目录下分别有version1，version2，version3三个目录，分别对应第一，第二，第三章中的，串行优化，并行优化，汇编优化三个不同阶段的版本。其中version3是最终版。

在每个version目录下，有Microsoft Visual Studio的项目文件，可以使用Microsoft Visual Studio直接打开。在对应的Release目录下有已经编译好的可执行文件。程序的优化选项可以在Microsoft Visual Studio中看到，具体的解释可以在第一章中查找。

此外，在压缩文件的根目录下有最终板的可执行文件potential_serial(final).exe方便进行测试。

附录B 参考文献

Intel C++ Compiler Documentation

Intel MKL Reference Manual

Intel MKL Technical User Notes

Getting Started Guide for Intel MKL

Getting Started with the VTune Performance Analyzer

VTune Performance Environment Help

Intel Thread Profiler for Windows

Intel Thread Checker for Windows

Intel 64 and IA-32 Architectures Software Developer’s Manual Volume 1 Basic Architecture

Intel 64 and IA-32 Architectures Software Developer’s Manual Volume 2A Instruction Set Reference, A-M

Intel 64 and IA-32 Architectures Software Developer’s Manual Volume 2B Instruction Set Reference, N-Z

Intel 64 and IA-32 Architectures Software Developer’s Manual Volume 3A System Programming Guide

Intel 64 and IA-32 Architectures Software Developer’s Manual Volume 3B System Programming Guide

Intel 64 and IA-32 Architectures Optimization Reference Manual

Using Spin-Loops on Intel Pentium 4 Processor and Intel Xeon Processor

FAST INVERSE SQUARE ROOT    CHRIS LOMONT

 

Source:http://blog.csdn.net/flyingdog/archive/2007/01/21/1489183.aspx