案例研究: 使用分布式优化框架在 Monte Carlo 欧式期权方面实现高级性能

1. 简介

Monte Carlo 使用统计计算方法解决复杂的科学计算问题。 它创新地使用随机数字模拟一个问题输入结果的不确定性,并通过处理重复的参数抽样获得一个确定的结果和解决一些以其他方式无法解决的问题。 该方法最早起源于上世纪 40 年代末,由参与“曼哈顿”计划的核物理学家们率先提出。 并采用摩纳哥最大的赌城 Monte Carlo 来命名。

来自加拿大的金融数学家 Phelim Boyle 被认为是在量化金融领域使用 Monte Carlo 方法的第一人。 他于 1977 年 5 月在金融经济周刊中发表了一篇题为期权 — Monte Carlo 方法的文章。在这篇文章中,他提出了采用模拟方法来解决欧式期权的价值问题,并且独立地验证了 Black-Scholes 公式。 自此,Monte Carlo 方法便在华尔街衍生工具定价及风险管理领域得到快速普及。 最近的一份报告显示,Monte Carlo 方法被用于解决高达 40% 的量化金融问题。 对任何期望实现高性能金融计算的公司来说,被称为 Monte Carlo 农场的大型集群计算机服务器安装几乎已经成为了一种必备配置。 1996 年,Mark Broadie 和 Paul Glasserman 将 Monte Carlo 应用到亚式期权领域。 2001 年,Francis Longstaff 和 Eduardo Schwartz 创建了第一个以 Monte Carlo 算法来定价美式期权。 2006 年,牛津大学的 Mike Giles 创建了多级别的 Monte Carlo 算法,旨在以经济高效的方法解决高维度问题,同时满足用户所指定的准确性要求。

尽管 Monte Carlo 利用精巧的数学表达式即可解决大量高维度问题,而且内存要求也低于直接 PDE 和 lattice方法,但是它同样需要强大的原始计算能力,尤其是在必须要确保确定性和准确性的情况下。 大多数情况下,它需要执行成百上千(有时候甚至是数百万)次抽样才能获得满意结果。 如果要获得多一位小数点的精确性,可能需要多执行 100 次的抽样。 由于许多问题在计算方面存在严峻的挑战,金融行业需要一款创新的架构来满足其原始计算能力需求。

1.1 英特尔® 集成众核架构 (英特尔® MIC)

英特尔® 集成众核架构 (英特尔® MIC) 为与 Monte Carlo 相关的应用提供了一种理想的执行平台。 由于采用与 IA32/英特尔 64 相同的指令集架构以及与英特尔® 至强™ 多核处理器相同的编程模式,英特尔® 至强融核™ 协处理器产品进一步扩展了并行执行基础设施,并支持在几乎不修改源代码的情况下,使高度并行化应用的性能远远超过经过改造的显卡(通常是指 GPGPU)。 用户获得的优势包括轻松的软件开发与部署,以及提升的软件应用性能。

英特尔至强融核协处理器(协处理器)是基于英特尔集成众核架构的首款产品。 它采用了一个全新的指令集,即英特尔® 初始众核指令(英特尔® IMCI)。 英特尔 IMCI 从 IA32/英特尔64 继承了包括 x87 指令在内的标量指令,同时采用了运行在全新设计的矢量处理单元 (VPU) 上的 512 位矢量指令。 然而,英特尔 IMCI 不支持英特尔 SSE,也不支持 SSE2 和英特尔 AVX。 对高度并行化应用而言,英特尔编译器可以生成在 VPU 上执行的矢量化代码。

1.1.1 处理器架构

在微处理器架构层面,协处理器就是一种 SMP 处理器,由 61 个以上的内核组成,这些内核被分成两个单向环。 处理器上有 8 个模上内存控制器,支持 16 个 GDDR5 通道,预计最高性能可达 5.5 GT/秒。 在内核层面,每个协处理器内核都具备完整的功能,并根据其自身权限进行排序,能够独立于其它内核之外运行 IA 指令。 此外,每个内核还支持硬件超线程技术及四个硬件环境或线程。 在任意时钟周期下,每个内核可从任意一个环境最多发出两条指令。

1.1.2 内核与矢量处理器单元

协处理器内核支持全部 16 个通用 64 位寄存器,以及与 64 位扩展相关的大多数新指令。 但是,只有初始英特尔® 众核指令集(英特尔® IMCI)是受支持的矢量指令。 尽管标量算术单元 (x87) 是集成式的并具备完整功能,但是协处理器内核中不支持英特尔® MMX™ 技术、英特尔® SSE 或英特尔® AVX。

其采用的是经过全新设计的矢量处理单元 (VPU)。 这是一种 512 位 SIMD 引擎,能够执行 16 宽单精度或 8 宽双精度浮点 SIMD 操作,并且完整支持全部四种舍入模式 (IEEE 754R)。 该 VPU 可以处理 32 位整数数据和 64 位长整数数据,与处理相同位的浮点数据的方式类似。

矢量处理单元内部有一个新型的可扩展数学单元,可加快单精度超越函数的实施速度:倒数、反平方根、以 2 为底的对数、以 2 为底的指数函数(使用极小极大二次多项式近似)。 这四个硬件执行的初等函数有着 4 到 8 个周期的延迟,可以实现 1 到 2 个周期的高吞吐率。 其他超越函数可以从这些初等函数中获取。

函数

数学

延迟周期

吞吐量. 周期

RECIP

1/x

 

1

RSQRT

1/√x

 

1

EXP2

2x

 

2

LOG2

log2x

 

1

1.1.3 高速缓存架构

英特尔至强融核协处理器的 1 级 (L1) 高速缓存设计旨在满足每个内核四种硬件环境中更高的工作集要求。 它具有一个 32 KB 的 L1 指令高速缓存和 32 KB 的 L1 数据高速缓存。 相关性为 8 路,带有 64 字节高速缓存行。 Bank(存储区)宽度为 8 字节。 数据返回可以是乱序的。 访问时间具有 3 个周期的延迟。

512 KB 统一二级 (L2) 高速缓存包括 64 个字节(每路,带 8 路相关性)、1024 个集、2 个存储区、32 GB(35 位)可缓存的地址范围。 预计闲置访问时间约为 80 个周期。

二级高速缓存具有一个流动的硬件预取器,它可以选择性地预取代码、读取高速缓存行,并将其 RFO 到二级高速缓存中。 它支持 16 个流,最多可提供 4-KB 的数据页面。 一旦检测到流方向,预取器便会发出最多 4 个多预取请求。 L2 高速缓存现在也获得了 ECC 支持。 一级和二级高速缓存的替换算法基于伪 LRU 实施。

下表概括了协处理器的高速缓存架构的主要参数:

参数

一级

二级

一致性

MESI

MESI

大小

32K (I) + 32K  (D)

512 K (合并)

相关性/行大小/库

8 路/64 字节/8

8 路/64 字节/8

访问时间(时钟数)

1

11

策略

伪 LRU

伪 LRU

工作周期

每时钟 1 个

每时钟 1 个

端口

读/写

读/写

1.2 英特尔® 至强™ 处理器和英特尔® 至强融核™ 协处理器的软件开发环境

1.2.1 英特尔® Parallel Studio XE 2013

英特尔于 2012 年 9 月发布了英特尔® Parallel Studio XE 2013。 这是一款面向英特尔多核与英特尔 MIC 架构产品线的集成软件开发工具套件。 在该软件包中,您可以找到英特尔® C++ 和 Fortran Composer XE 以及英特尔® VTune™ Amplifier XE 应用性能分析器。 同时包含在内的还有英特尔® 性能库,诸如英特尔® 线程构建模块(英特尔® TBB)和英特尔® 数学核心函数库(英特尔® MKL)。开发人员可以使用这些性能库构建应用和其他库。 集成式英特尔多核与英特尔® MIC 架构支持是最新版本所具备的主要特性之一。 自 2013 年 1 月至今,我们已经发布了 5 个不同的更新。 本文中提到的程序已经内置了第 5 个更新。 本文中提到的性能是通过运行构建有英特尔® C++ 英特尔® 64 Composer XE 的 13.1 更新 3 的可执行文件后测试所得。

1.2.2 英特尔至强主机处理器系统

主机系统就是基于英特尔至强的双插槽平台,通过 PCI Express* (PCIe) 接口与一个英特尔至强融核协处理器连接在一起。 该主机的每个插座中均安装一路英特尔® 至强™ E5-2670 八核处理器 (2.6Ghz)。 英特尔至强 E5-2670 支持面向 SIMD 并行机制的英特尔高级矢量扩展指令集。

集成软件开发工具堆栈运行在基于英特尔至强的主机系统上。 应用开发人员可以为英特尔至强融核协处理器创建卸载应用或本地应用。

一个卸载应用启动后首先在主机系统上运行。 根据应用开发人员的指示或指令,程序的一部分可以在协处理器上运行。 编译器为主机和协处理器创建一个带有二进制代码的可执行文件。 集成工具的运行时库负责为协处理器提取和拷贝二进制文件,设置输入和输出数据,并调用协处理器程序。

本地应用启动后只在协处理器卡上运行。 程序设计员告知编译器整个程序计划将在协处理器上运行,采用的方法是在编译器调用行使用一个特殊的开关 (–mmic )。 程序设计员负责将可执行文件连同共享库从主机拷贝到协处理器上,准备输入数据和输出数据,同时调用协处理器上的程序。

本文仅讨论如何为协处理器创建本地应用。 卸载应用不在本文的讨论范围之内。

1.2.3 英特尔至强融核协处理器

协处理器是一个可以运行自己的操作系统的计算设备,该操作系统是 Linux* 操作系统的一个改进版本,被称为多核平台系统软件或 MPSS。 协处理器设备不能自行启动。 主机系统通过标准的 Linux 设备操作命令来启动或关闭协处理器。 可以与英特尔 Parallel Studio XE 2013 的更新 5 配合工作的 MPSS 版本是 2.1.6720-13。 在本文中,我们使用配有 61 个内核的 7120p 型协处理器。 每个内核配备 16 GB 的 GDDR 内存,运行频率达 1.238 GHz。 内核的内存速度达 5.5 GHz。

2. Monte Carlo 欧式期权定价算法

在本节中,我们会讲解一些基础的金融衍生工具,特别是股票期权定价。 我们将介绍使用 Monte Carlo 方法对输入的不确定性进行建模、如何抽取样本,以及如何根据期权规范来确定收益函数。 最后,我们使用标准 C++ 来实施 Monte Carle。 我们还将着重讲解这种实施方法的性能瓶颈问题。

2.1 期权定价背景

在金融领域,衍生品是一种金融工具,其价值取决于其它基础变量的价值。 通常,标的衍生品的变量是交易资产的价格。 例如,作为一种衍生品,股票期权的价值取决于股票价格本身。 并非所有衍生品变量都依赖于交易资产。 例如,某些变量可能是一个特定度假村的降雪,其它变量可能是一个特定时间间隔内的平均温度等。

期权是一种衍生品,它规定了双方以一个参考价格就某种资产所达成的未来交易合约(称作行使价)。 期权的买方获得行使该交易的权利(而不是义务),卖方则承担履行交易的相应义务。 期权分为两种类型: 看涨期权赋予合约持有者(即买方)按照约定的价格在约定日期前购买标的资产的权利。 看跌期权赋予合约持有者(即买方)按照约定的价格在约定日期前卖出标的资产的权利。 合约中的价格被称为行使价格。 合约中的日期被称为到期日。 欧式期权仅能在到期日当天行使。 美式期权可在到期日或之前的任一交易日行使。

2.2 使用 Monte Carlo 方法定价欧式股票期权

Monte Carlo 模拟使用风险中性估值方法来确定期权的价值。 它对一个路径进行取样以获得风险中性环境中的预计收益,然后使用无风险利率,根据现值对收益进行折算。 我们来考虑一个股票期权,该股票的当前价格为 S,并在时间 T 提供收益。假设利率是常量,我们可以对其估值如下:

  1. 在风险中性环境中对 S 的随机路径进行取样。
  2. 计算衍生品的收益。
  3. 重复第一步和第二步,在风险中性环境中获得大量衍生品收益价值的样本。
  4. 计算样本收益的平均值,获得风险中性环境中的预估收益。
  5. 使用无风险利率对预估收益进行折算,获得衍生品的预估价值。

假设风险中性环境中的标的市场变量后的过程为:

其中 dWt 代表 Wiener 过程、µ 代表风险中性环境中的预估回报,σ 代表波动性。 为了模拟 lnS 路径,我们可以将衍生品的寿命除以长度 t 的 N 个短间隔,近似方程 3.1 为:

或者为:

其中,ε 为来自 φ (0,1) 的随机样本。 这有助于根据 S 的初始值计算在 S 在时间 Δt 的值,以及根据时间 Δt 时的值计算 2 Δt 的值等等。

对于 1 和 M 之间的每个 k,此处每个 εi 均从标准正态分布提取。

我们知道欧式期权在到期日的价值。看涨期权和看跌期权的等式为:

借助 Monte Carlo,我们可以根据标的 Wiener 过程生成带有所需 (0, 1) 分布的 M 个数字样本,然后根据每个样本的值计算到期日股票收益的平均值:

此外还可以获得类似的欧式看跌期权的函数。

结果仍然是期权在时间 T 时的未来价值。使用因数对该值进行折现,即可获得欧式看涨期权的现值。

它遵循中心极限定理,因此可将标准偏差减半,所需的样本路径需要增加四倍。 换句话说,Monte Carlo 的标准误差为

Monte Carlo 模拟的优势在于,当收益取决于路径(后面为截至到期日时的基础变量 S)以及在期权的整个生命周期内多次发生收益时,可以使用该模拟方法。 当收益函数包含多个独立变量时,该方法的作用尤为明显。 当所有其他的分析方法都失效时,Monte Carlo 将成为您唯一的选择。

2.3 实施 Monte Carlo 欧式期权定价算法

只要我们具备上一章节中提到的收益函数,那么实施 Monte Carlo 欧式期权定价便会相对容易一些。 首先,我们要从φ (0, 1)找到一个随机数,然后使用收益函数计算标的期权的预估值和可信区间。 C/C++ 实施如下所示:


        for (int pos = 0; pos < RAND_N; pos++)
        {
            float callValue =  max(0.0, Sval *exp(MuByT + VBySqrtT *  gen()) - Xval);
            val  += callValue;
            val2 += callValue * callValue;
        }

我们来看一个假设情况 — 一家公司需要计算数百万种金融工具的欧式期权。 每种工具都有一套参数:当前价格、行使价格和期权到期时间。 该公司想要使用 Monte Carlo 模拟为每套数据计算欧式看涨期权。 在下面的例子中,我们展示了使用 265 K (或 266,144)路径长度针对 1500 万期权数据集对欧式看涨期权进行定价的代码片段。


#include "MonteCarlo.h"
#include <math.h>
#include <tr1/random>

#ifndef max
  #define max(a,b) (((a) > (b)) ? (a) ; (b))
#endif
  
void MonteCarlo(
    float *h_CallResult,
    float *h_CallConfidence,
    float *S,
    float *X,
    float *T
)
{
    typedef std::tr1::mt19937                     ENG;    // Mersenne Twister
    typedef std::tr1::normal_distribution<float> DIST;    // Normal Distribution
    typedef std::tr1::variate_generator<ENG,DIST> GEN;    // Variate generator

    ENG  eng;
    DIST dist(0,1);
    GEN  gen(eng,dist);


    for(int opt = 0; opt < OPT_N; opt++) 
    {
        float VBySqrtT = VOLATILITY * sqrt(T[opt]);
        float MuByT = (RISKFREE - 0.5 * VOLATILITY * VOLATILITY) * T[opt];
        float Sval = S[opt];
        float Xval = X[opt];
        float val = 0.0, val2 = 0.0;

        for(int pos = 0; pos < RAND_N; pos++)
        {
            float callValue =  max(0.0, Sval *exp(MuByT + VBySqrtT *  gen()) - Xval);
            val  += callValue;
            val2 += callValue * callValue;
        }

        float exprt = exp(-RISKFREE *T[opt]);
        h_CallResult[opt] = exprt * val / (float)RAND_N;
        float stdDev = sqrt(((float)RAND_N * val2 - val * val)/ ((float)RAND_N * (float)(RAND_N - 1)));
        h_CallConfidence[opt] = (float)(exprt * 1.96 * stdDev / sqrtf((float)RAND_N));
    }
}

值得注意的是,上述代码顺序使用的是 C++(使用 TR1(C++ 标准委员会技术报告 1)扩展)中的随机数生成工具。 C++ 随机数生成类别和函数在 <random> 标头中定义并包含在命名空间 std::tr1 中。

对于任意伪随机数生成软件来说,其核心是能够生成均匀分布的随机整数。 然后可以在 bootstrap 过程中使用这些软件来生成均匀分布的浮点数。 通过转换和舍选法等,这些均匀分布的浮点数列可用于从其它分布中生成。

C++ TR1 可使您选择多种核心生成器,即“引擎”。 GCC 4.3.x 和 Visual Studio* 2008 特性包中支持以下四种引擎。

  • linear_congruential
  • mersenne_twister
  • subtract_with_carry
  • subract_with_carry_01

C++ TR1 库通过分布类别支持非均匀的随机数生成。 这些类别通过 operator() 方法返回随机样本。

温度类别 variate_generator 描述了一个含有引擎和分布的对象,并通过将封装的引擎对象传输到分布对象的 operator() 来生成值。

最后使用 GCC 4.4.6 进行编译,在基于八核英特尔至强 E5-2670 处理器 (2.6 GHz) 的两路主机系统上首次实施。Monte Carlo 欧式期权定价的结果为 34 opt/sec 以上。


[sli@localhost step1]$ ./MonteCarlo
Monte Carlo European Option Pricing in Single Precision
Pricing 32768 options with path length of 262144.
Completed in 955.2525 seconds.
Computation rate -   34.303 options per second.

很明显,该实施具有很大的改进空间。 在下面的章节中,我们会通过一个结构化框架测试这个程序,并尝试改进其性能。

3. 分布式优化框架和性能优化

与其它所有科研过程一样,性能优化需要采用系统化和结构化方法。 在此章节中,我们重点介绍可帮助您针对应用性能优化采取整体方法的代码现代化优化框架。 该框架的目的是通过最佳的工具和现有的库实现最高的性能。 该框架包含五个优化步骤。 每个步骤都通过使用不同的技术,力争在一个方向上提高应用性能。

本方法的目标就是解决所有与应用性能有关的问题和障碍,并在英特尔架构上实现尽可能高的性能,同时利用所有可能的并行执行资源。

3.1 分布式代码现代化框架

分布式代码现代化由五个支持并行的步骤组成,可帮助性能工程师在尽可能短的时间内实现最佳的应用性能。 换句话说,它支持程序在执行环境中最大限度地使用所有的并行硬件资源。

  • 第一步是选择优化开发环境。 该环境应能够创建优化代码并支持您使用现有的优化库。
  • 第二步是优化您正在执行的操作。 它确保预定的操作以最佳的方式执行。
  • 第三步全部与矢量化有关,即我们如何在您的应用中利用 SIMD 指令以及探索同步数据并行化,并以编译器能够理解的方式表达出来。 我们的目标就是为使用一条线程发挥单个内核的最大性能。
  • 第四步即线程并行化,此时我们使用并行性解决同步操作的问题。 最终目标是充分利用所有的内核。
  • 第五步是把您的应用从多核扩展到英特尔 MIC 上。 这一步对高度并行化应用来说极为重要。 其目的是希望使用微架构的独有特性并根据这些特性执行其他优化,例如,内存带宽和处理能力间的差异、SIMD 对齐要求的变化、每条线程的高速缓存等。此外,还希望在执行目标从英特尔架构的一种风格(英特尔至强处理器)改变为另一种(英特尔至强融核协处理器)时能最大限度地减少变化和提高性能。

3.2 第 1 步: 利用优化工具和库

在开始优化项目时,您需要选择一个优化开发环境。 该选择对于后续步骤具有重要的影响。 它不仅会影响您得到的结果,还能大幅减少您的工作量。正确的优化开发环境可以为您提供出色的编译器工具、现成的优化库、调试工具和性能评测工具,帮助您准确地查看代码在运行时正在做什么。

有时候,单个环境可能很难满足您的所有需求。 在这种情况下,您必须根据不同的工具对解决方案进行整合。 例如,如果您想要优化 Java* 程序,而该程序的一些 C/C++ 代码对性能要求极高,那么您需要使用 Java SDK、C/C++ 编译器、JVM 环境,还可能需要一个系统级的分析工具。 其他时间,解决方案可能很明显。 英特尔 C++ Composer XE 集成了多核和众核开发环境,是目前唯一具备矢量化功能的应用编译器。 如果您的项目要求将高度并行化应用从英特尔至强处理器扩展到英特尔至强融核协处理器,英特尔 C++ Composer XE 2013 是目前唯一可选的开发工具。 集成式性能分析工具之所以成为最佳的工具,是因为它提供了性能监控事件,同时可准确告诉您如何分配您的优化资源。 它就像一个夜视镜,在您需要时为您照亮目标。

Linux* 上的英特尔 C++ Composer XE 2013 兼容随 RHEL 6.3 服务器版一起发布的 GCC 4.4.6。 英特尔 C++ Composer 提供了与 g++ 类似的工具 icpc,两者可以编译相同的内容。 例如,您可以发布 icpc -o MonteCarlo -O2 MonteCarloStage1.cpp,而不是使用 g++ -o MonteCarlo -O2 MonteCarloStage1.cpp。 当您运行英特尔 C++ Composer 创建的二进制代码时,您甚至无需使用源代码就能获得 1.37 倍的即时性能提升。 我们目标就是在确保计算准确性的同时提升性能。


[sli@localhost step1]$ make -B CXX=icpc
icpc -c -O2  -o Driver.o Driver.cpp
icpc -c -O2  -o MonteCarloStep1.o MonteCarloStep1.cpp
icpc Driver.o MonteCarloStep1.o -o MonteCarlo
[sli@cthor-knc14 step1]$ ./MonteCarlo
Monte Carlo European Option Pricing in Single Precision
Pricing 32768 options with path length of 262144.
Completed in 694.6562 seconds.
Computation rate -   47.172 options per second.

英特尔 MKL 经过预编译,包含基本的数学例程,例如 BLAS、LAPACK 和随机数生成函数。 MKL 为这些函数提供了 C 和 Fortran 接口。 您可以从 C/C++ 或 Fortran 调用这些函数。

随机数生成 (RNG) 效率反映了所有 Monte Carlo 模拟的性能。 如同 C++ TR1 一样,英特尔 MKL 也提供随机数生成机制。 像 C++ TR1 一样,如果适用,英特尔 MKL 还允许您独立地选择不同的基础 RNG 引擎和不同的分布(乘以两个浮点和两个 64 位和 32 位的整数数据)。 与 C++ TR1 不同,英特尔 MKL 的 RNG 接口允许您在 C、C++ 和 Fortran 中使用它们。 与其他编译库类似,您需要首先包含它的接口声明 .h 文件,然后在您的 C/C++ 实施文件中进行英特尔 MKL API 调用,最后再链接预编译的英特尔 MKL 库。 在我们的例子中,我们需要包含三个英特尔 MKL 库:mkl_intel_lp64、mkl_sequential 和 mkl_core。 您可以连同 –lpthread 一起列举这些库,诸如 -lmkl_intel_lp64 -lmkl_sequential -lmkl_core –lpthread,或您可以将 –mkl 传递到链接器,该链接器将为您进行处理。

其中一个主要区别是,在 C++ TR1 中,借助 variate_generator 类的 gen() 方法调用,一次仅提供一个随机数。 而在英特尔 MKL 中,创建 RNG 引擎对象的流程由创建 RNG 流取代;创建分布对象的流程由选择正确的 RNG 接口 API 调用取代。 一个英特尔 MKL RNG 接口调用可提供任意数量的随机数。 在我们的例子中,256 K 随机数在一个 RNG 接口调用中提供。


#include "MonteCarlo.h"
#include "math.h"
#include "mkl_vsl.h"
#define RANDSEED 123

#ifndef max
  #define max(a,b) (((a) > (b)) ? (a) : (b))
#endif
  
void MonteCarlo(
    float *h_CallResult,
    float *h_CallConfidence,
    float *S,
    float *X,
    float *T
)
{
    float random [RAND_N];VSLStreamStatePtr Randomstream;vslNewStream(&Randomstream, VSL_BRNG_MT19937, RANDSEED);vsRngGaussian (VSL_METHOD_SGAUSSIAN_ICDF, Randomstream, RAND_N, random, 0.0, 1.0);

    for(int opt = 0; opt < OPT_N; opt++) 
    {
        float VBySqrtT = VOLATILITY * sqrt(T[opt]);
        float MuByT = (RISKFREE - 0.5*VOLATILITY*VOLATILITY)*T[opt];
        float Sval = S[opt];
        float Xval = X[opt];
        float val = 0.0, val2 = 0.0;

        for(int pos = 0; pos < RAND_N; pos++)
        {
            float callValue =  max(0.0, Sval*exp(MuByT+VBySqrtT*random[pos])-Xval);
            val  += callValue;
            val2 += callValue * callValue;
        }

        float exprt = exp(-RISKFREE *T[opt]);
        h_CallResult[opt] = exprt * val / (float)RAND_N;
        float stdDev = sqrt(((float)RAND_N*val2-val*val)/((float)RAND_N*(float)(RAND_N-1)));
        h_CallConfidence[opt] = (float)(exprt * 1.96 * stdDev / sqrtf((float)RAND_N));
    }
    vslDeleteStream(&Randomstream);
}

实际上,通过使用英特尔 MKL RNG,我们可以减少代码行的数量。 性能上的差异十分明显。 与原代码相比实现了 5.53 倍以上的性能提升。


[sli@cthor-knc14 .solution]$ ./MonteCarlo
Monte Carlo European Option Pricing in Single Precision
Pricing 32768 options with path length of 262144.
Completed in 125.5210 seconds.
Computation rate -  261.056 options per second.

与原始代码相比,将 RNG 从 C++ TR1 改变成英特尔 MKL 不仅能实现 5.53 倍的性能提升,还能支持编译器进行内联函数调用和对内层循环中的代码进行矢量化,我们将在第 3 步中讨论该内容。

目前为止,我们应能够将所提供的优化软件用作优化的软件开发环境。 这可以为我们节省大量的时间,同时有助于我们专注于开发自己的 IP 和解决新的问题。 更重要的是,如果我们希望实现更多的性能提升,我们就必须要改变一下源代码。 在我们应用任意并行机制之前,首先要精简和清理我们的串行代码。

3.3 第 2 步: 标量串行优化

使用标量串行优化检查是否存在冗余,确保您的代码运行时获得最高的效率。 在该步骤中,您需要确保根据现有的问题选择数据结构的精度和数值方法的准确性,同时避免进行过度选择,确保可以使用较低成本的运算(例如单精度)代替高成本的双精度运算。

C/C++ 依然保留了早期 C 的弱模式培养。在某些情况下,人们认为可以通过自动升级实现较高的准确性。 尽管该实践对提升准确性的程度还不确定,但是它肯定不利于性能提升。 下表列出了标量和串行优化阶段所要注意的其它事项。

在 Monte Carlo 应用中,我们有一个超越函数调用,内部循环中以自然数为底数的指数函数。 我们为避免该函数调用所做的任何改善都能最终促成较大的性能提升。


#include "MonteCarlo.h"
#include "math.h"
#include "mkl_vsl.h"
#define RANDSEED 123

#ifndef max
  #define max(a,b) (((a) > (b)) ? (a) : (b))
#endif
  
static const float  RVV = RISKFREE-0.5f*VOLATILITY*VOLATILITY;static const float  INV_RAND_N = 1.0f/RAND_N;static const float  F_RAND_N = static_cast<float>(RAND_N);static const float STDDEV_DENOM = 1 / (F_RAND_N * (F_RAND_N - 1.0f));static const float CONFIDENCE_DENOM = 1 / sqrtf(F_RAND_N);

void MonteCarlo(
    float *h_CallResult,
    float *h_CallConfidence,
    float *S,
    float *X,
    float *T
)
{
    float random [RAND_N];
    VSLStreamStatePtr Randomstream;
    vslNewStream(&Randomstream, VSL_BRNG_MT19937, RANDSEED);

    vsRngGaussian (VSL_METHOD_SGAUSSIAN_ICDF, Randomstream, RAND_N, random, 0.0f, 1.0f);

    for(int opt = 0; opt < OPT_N; opt++) 
    {
        float VBySqrtT = VOLATILITY * sqrtf(T[opt]);float MuByT = RVV * T[opt];
        float Sval = S[opt];
        float Xval = X[opt];
        float val = 0.0, val2 = 0.0;

        for(int pos = 0; pos < RAND_N; pos++)
        {
            float callValue =  max(0.0, Sval *expf(MuByT + VBySqrtT *  random[pos]) - Xval);
            val  += callValue;
            val2 += callValue * callValue;
        }

        float exprt = expf(-RISKFREE *T[opt]);h_CallResult[opt] = exprt * val * INV_RAND_N;float stdDev = sqrtf((F_RAND_N * val2 - val * val)* STDDEV_DENOM);h_CallConfidence[opt] = (exprt * 1.96f * stdDev * CONFIDENCE_DENOM);
    }
    vslDeleteStream(&Randomstream);
}

请注意,我们已经从循环中移除了“除以常数”运算。 这些运算在内层循环之外执行,因此编译器通常会错过它们。 性能得到了适度提升 (41.11%)。


[sli@localhost step3]$ make -B
icpc -c -O3 -ipo -fimf-precision=low -fimf-domain-exclusion=31 -fimf-accuracy-bits=11 -no-prec-div -no-prec-sqrt -o Driver.o Driver.cpp
icpc -c -O3 -ipo -fimf-precision=low -fimf-domain-exclusion=31 -fimf-accuracy-bits=11 -no-prec-div -no-prec-sqrt -o MonteCarloStep2.o MonteCarloStep2.cpp
icpc Driver.o MonteCarloStep2.o -o MonteCarlo -mkl
[sli@cthor-knc14 .solution]$ ./MonteCarlo
Monte Carlo European Option Pricing in Single Precision
Pricing 32768 options with path length of 262144.
Completed in  88.9513 seconds.
Computation rate -  368.381 options per second.

这表明,无论你在代码上进行多少努力,如果该代码在标量和串行模式中运行,其改进空间将会非常有限。 因此我们必须借助并行化。

3.4 第 3 步: 矢量化

矢量化在不同的环境中有着不同的含义。 在本文中,矢量化就是编译器生成的 SIMD 指令,这些指令可以使用 CPU 中的矢量寄存器。 这意味着可以利用一些多数据元素指令的同步执行。

您可以使用多种不同的方式将矢量化引入到您的程序中,从使用处理器固有函数到使用英特尔® Cilk™ Plus 数组符号。 由于编程人员对于生成代码的控制数量、语法表达以及串行程序所需的变化数量存在差异,因此这些基于编译器的矢量化技巧各有不同。

在我们使用编译器对串行代码进行矢量化处理以及生成 SIMD 指令之前,我们必须确保内存对齐正确无误。 没有对齐的内存访问会导致高速缓存行隔离、对象代码重复,从而影响应用性能,严重的时候还会产生处理器故障。 确保内存对齐的一种方法是:始终请求并使用明确的对齐要求。 借助英特尔 C++ Composer XE 2011,您可以通过为内存定义添加 __attribute__(align(64)) 前缀的方法请求静态分配的内存。 对为英特尔至强融核协处理器矢量寄存器指派的内存来说,64 字节边界是最低的对齐要求。 此外,您还可以使用 _mm_malloc 和 _mm_free 来请求和提供动态分配的内存。 Linux 操作系统上存在英特尔和 GCC 支持的内存分配规则。

除内存对齐外,基于英特尔编译器的矢量化只有在热循环中的函数调用降至最低时才会发挥出最大作用,这主要是因为函数调用也会产生开销。 并非所有的函数调用都允许 CPU 的矢量引擎继续以 SIMD 模式执行。 您可以内联尽可能多的函数,这是因为函数内联能够避免函数调用开销,并支持基于编译器的矢量化器同时分析和矢量化被调用者和调用者代码。


#include "MonteCarlo.h"
#include "math.h"
#include "mkl_vsl.h"
#define RANDSEED 123

static const float  RVV = RISKFREE-0.5f*VOLATILITY*VOLATILITY;
static const float  INV_RAND_N = 1.0f/RAND_N;
static const float  F_RAND_N = static_cast<float>(RAND_N);
static const float STDDEV_DENOM = 1 / (F_RAND_N * (F_RAND_N - 1.0f));
static const float CONFIDENCE_DENOM = 1 / sqrtf(F_RAND_N);

void MonteCarlo(
    float *h_CallResult,
    float *h_CallConfidence,
    float *S,
    float *X,
    float *T
)
{
    __attribute__((align(4096))) float random [RAND_N];
    VSLStreamStatePtr Randomstream;
    vslNewStream(&Randomstream, VSL_BRNG_MT19937, RANDSEED);

    vsRngGaussian (VSL_METHOD_SGAUSSIAN_ICDF, Randomstream, RAND_N, random, 0.0f, 1.0f);
    for(int opt = 0; opt < OPT_N; opt++) 
    {
        float VBySqrtT = VOLATILITY * sqrtf(T[opt]);
        float MuByT = RVV * T[opt];
        float Sval = S[opt];
        float Xval = X[opt];
        float val = 0.0, val2 = 0.0;

#pragma vector aligned#pragma simd reduction(+:val) reduction(+:val2)#pragma unroll(4)
        for(int pos = 0; pos < RAND_N; pos++)
        {
            float callValue  = Sval * expf(MuByT + VBySqrtT * random[pos]) - Xval;
            callValue = (callValue > 0) ? callValue : 0;
            val  += callValue;
            val2 += callValue * callValue;
        }

        float exprt = expf(-RISKFREE *T[opt]);
        h_CallResult[opt] = exprt * val * INV_RAND_N;
        float stdDev = sqrtf((F_RAND_N * val2 - val * val)* STDDEV_DENOM);
        h_CallConfidence[opt] = (exprt * 1.96f * stdDev * CONFIDENCE_DENOM);
    }
    vslDeleteStream(&Randomstream);
}

#include <stdlib.h>
#include <stdio.h>
#include <sys/time.h>
#include <string.h>
#include <math.h>
#include "MonteCarlo.h"
#include <iostream>
using namespace std;

#define SIMDALIGN 64

double second()
{
    struct timeval tv;
    gettimeofday(&tv, NULL);
    return (double)tv.tv_sec + (double)tv.tv_usec / 1000000.0;
}

inline float RandFloat(float low, float high){
    float t = (float)rand() / (float)RAND_MAX;
    return (1.0f - t) * low + t * high;
}

///////////////////////////////////////////////////////////////////////////////
// Polynomial approximation of cumulative normal distribution function
///////////////////////////////////////////////////////////////////////////////

double CND(double d){
    const double       A1 = 0.31938153;
    const double       A2 = -0.356563782;
    const double       A3 = 1.781477937;
    const double       A4 = -1.821255978;
    const double       A5 = 1.330274429;
    const double RSQRT2PI = 0.39894228040143267793994605993438;

    double
        K = 1.0 / (1.0 + 0.2316419 * fabs(d));

    double
        cnd = RSQRT2PI * exp(- 0.5 * d * d) *
        (K * (A1 + K * (A2 + K * (A3 + K * (A4 + K * A5)))));

    if(d > 0)
        cnd = 1.0 - cnd;

    return cnd;
}

void BlackScholesFormula(
    double& callResult,
    double Sf, //Stock price
    double Xf, //Option strike
    double Tf, //Option years
    double Rf, //Riskless rate
    double Vf  //Volatility rate
){
    double S = Sf, X = Xf, T = Tf, R = Rf, V = Vf;

    double sqrtT = sqrt(T);
    double    d1 = (log(S / X) + (R + 0.5 * V * V) * T) / (V * sqrtT);
    double    d2 = d1 - V * sqrtT;
    double CNDD1 = CND(d1);
    double CNDD2 = CND(d2);

    double expRT = exp(- R * T);
    callResult   = (S * CNDD1 - X * expRT * CNDD2);
}

int main(int argc, char* argv[])
{
    float 
        *CallResultParallel,
        *CallConfidence,
        *StockPrice,
        *OptionStrike,
        *OptionYears;
    double
	sTime, eTime;
    int 
        mem_size, rand_size, verbose = 0; 
  
    const int RAND_N = 1 << 18;

    if (argc > 2)
    {
      	printf("usage: MonteCarlo <verbose>  where verbose = 1 for validtating result, the default is not to validate result. \n");
      	exit(1);
    }
    if (argc == 1)
	verbose = 0;
    if (argc == 2)
	verbose = atoi(argv[1]);

    printf("Monte Carlo European Option Pricing in Single Precision\n");	
    mem_size = sizeof(float)*OPT_N;
    rand_size = sizeof(float)*RAND_N;

    CallResultParallel = (float *)_mm_malloc(mem_size, SIMDALIGN);CallConfidence     = (float *)_mm_malloc(mem_size, SIMDALIGN);StockPrice         = (float *)_mm_malloc(mem_size, SIMDALIGN);OptionStrike       = (float *)_mm_malloc(mem_size, SIMDALIGN);OptionYears        = (float *)_mm_malloc(mem_size, SIMDALIGN);

    if (verbose)
    {
        printf("...generating the input data.\n");
    }
    for(int i = 0; i < OPT_N; i++)
    {
        CallResultParallel[i] = 0.0;
        CallConfidence[i]= -1.0;
        StockPrice[i]    = RandFloat(5.0f, 50.0f);
        OptionStrike[i]  = RandFloat(10.0f, 25.0f);
        OptionYears[i]   = RandFloat(1.0f, 5.0f);
    }

    printf("Pricing %d options with path length of %d.\n", OPT_N, RAND_N);
    sTime = second();
    MonteCarlo(
        CallResultParallel,
        CallConfidence,
        StockPrice,
        OptionStrike,
        OptionYears);
    eTime = second();
    printf("Completed in %8.4f seconds.\n",eTime-sTime );
    printf("Computation rate - %8.3f options per second.\n", OPT_N/(eTime-sTime));

    if (verbose)
    {
        double
           delta, sum_delta, sum_ref, L1norm, sumReserve;
        double CallMaster;

        sum_delta = 0;
        sum_ref   = 0;
        sumReserve = 0;

        for(int i = 0; i < OPT_N; i++)
        {
            BlackScholesFormula(CallMaster,
                (double) StockPrice[i], 
                (double) OptionStrike[i],
                (double) OptionYears[i],
                (double) RISKFREE,
                (double) VOLATILITY);
            delta = fabs(CallMaster - CallResultParallel[i]);
            sum_delta += delta;
            sum_ref   += fabs(CallMaster);
            if(delta > 1e-6) 
                sumReserve += CallConfidence[i] / delta;
        }
        sumReserve /= (double)OPT_N;
        L1norm = sum_delta / sum_ref;
        printf("L1 norm: %E\n", L1norm);
        printf("Average reserve: %f\n", sumReserve);

        printf("...freeing CPU memory.\n");
        printf((sumReserve > 1.0f) ? "PASSED\n" : "FAILED\n");
    }
    _mm_free(CallResultParallel);_mm_free(CallConfidence);_mm_free(StockPrice);_mm_free(OptionStrike);_mm_free(OptionYears);
    return 0;
}

icpc -c -O3 -ipo -xAVX -fimf-precision=low -fimf-domain-exclusion=31 -fimf-accuracy-bits=11 -no-prec-div -no-prec-sqrt -fno-alias -vec-report2 -o Driver.o Driver.cpp
icpc: remark #10346: optimization reporting will be enabled at link time when performing interprocedural optimizations
icpc -c -O3 -ipo -xAVX -fimf-precision=low -fimf-domain-exclusion=31 -fimf-accuracy-bits=11 -no-prec-div -no-prec-sqrt -fno-alias -vec-report2 -o MonteCarloStep4.o MonteCarloStep4.cpp
icpc: remark #10346: optimization reporting will be enabled at link time when performing interprocedural optimizations
icpc Driver.o MonteCarloStep4.o -o MonteCarlo -mkl
Driver.cpp(130): (col. 5) remark: loop was not vectorized: existence of vector dependence.
Driver.cpp(141): (col. 5) remark: SIMD LOOP WAS VECTORIZED.
Driver.cpp(141): (col. 5) remark: loop was not vectorized: not inner loop.
Driver.cpp(161): (col. 9) remark: LOOP WAS VECTORIZED.
[sli@localhost step4]$ ./MonteCarlo
Monte Carlo European Option Pricing in Single Precision
Pricing 32768 options with path length of 262144.
Completed in  10.9624 seconds.
Computation rate - 2989.117 options per second.

矢量化可实现 8.11 倍的性能提升 — 比预期的 8 倍还要高一些。 编译器变化、英特尔 MKL、串行化和矢量化共帮助我们实现将近 87.14 倍的性能提升。 而且,我们仅使用一条基于一个内核运行的线程就实现了这些优势。

3.5 第 4 步: 并行化

现在的微处理器全部都是多核处理器。 它们从超线程技术演进而来,包括双核、四核、八核,甚至十核、十二核。 英特尔至强融核协处理器家族的第一款产品使用的核心数多达 61 个。 设计一款可随着内核数量进行扩展的程序已经成为高性能软件开发的基本要求。 通常的做法是将整个任务划分为多个小型任务,隐含或明确地创建较小的轻型进程或线程来执行这些任务,并将每条线程与硬件线程或操作系统调度程序实体绑定。 该进程也被称作线程化或多线程化,而执行它的程序也被称为多线程程序。

借助英特尔 C++ Composer XE 2013,你可以根据您需要的控制数量和现有问题的要求选择通过不同方法创建多线程程序。

对 Monte Carlo 实施进行线程化处理时,我们正在寻找一种理想的方式,以避免我们在第 3 步中获得的矢量化串行代码发生改变。 我们的 Monte Carlo 欧式期权也非常容易划分成子任务。 由于存在将近 32K 数据集,在双路至强处理器上,我们可以将代码分为 16 个子任务,每个子任务负责 2K 数据集。 OpenMP* 非常符合我们的标准。


#include "MonteCarlo.h"

#include "math.h"
#include "mkl_vsl.h"
#define RANDSEED 123

static const float  RVV = RISKFREE-0.5f*VOLATILITY*VOLATILITY;
static const float  INV_RAND_N = 1.0f/RAND_N;
static const float  F_RAND_N = static_cast<float>(RAND_N);
static const float STDDEV_DENOM = 1 / (F_RAND_N * (F_RAND_N - 1.0f));
static const float CONFIDENCE_DENOM = 1 / sqrtf(F_RAND_N);

void MonteCarlo(
    float *h_CallResult,
    float *h_CallConfidence,
    float *S,
    float *X,
    float *T
)
{
    __attribute__((align(4096))) float random [RAND_N];
    VSLStreamStatePtr Randomstream;
    vslNewStream(&Randomstream, VSL_BRNG_MT19937, RANDSEED);

    vsRngGaussian (VSL_METHOD_SGAUSSIAN_ICDF, Randomstream, RAND_N, random, 0.0f, 1.0f);
#pragma omp parallel for
    for(int opt = 0; opt < OPT_N; opt++) 
    {
        float VBySqrtT = VOLATILITY * sqrtf(T[opt]);
        float MuByT = RVV * T[opt];
        float Sval = S[opt];
        float Xval = X[opt];
        float val = 0.0, val2 = 0.0;

#pragma vector aligned
#pragma simd reduction(+:val) reduction(+:val2)
#pragma unroll(4) 
        for(int pos = 0; pos < RAND_N; pos++)
        {
            float callValue  = Sval * expf(MuByT + VBySqrtT * random[pos]) - Xval;
            callValue = (callValue > 0) ? callValue : 0;
            val  += callValue;
            val2 += callValue * callValue;
        }

        float exprt = expf(-RISKFREE *T[opt]);
        h_CallResult[opt] = exprt * val * INV_RAND_N;
        float stdDev = sqrtf((F_RAND_N * val2 - val * val)* STDDEV_DENOM);
        h_CallConfidence[opt] = (exprt * 1.96f * stdDev * CONFIDENCE_DENOM);
    }
    vslDeleteStream(&Randomstream);
}

[sli@localhost step5]$ make -B
icpc -c -g -O3 -ipo -openmp -xAVX -fimf-precision=low -fimf-domain-exclusion=31 -fimf-accuracy-bits=11 -no-prec-div -no-prec-sqrt -fno-alias -vec-report2 -o Driver.o Driver.cpp
icpc: remark #10346: optimization reporting will be enabled at link time when performing interprocedural optimizations
icpc -c -g -O3 -ipo -openmp -xAVX -fimf-precision=low -fimf-domain-exclusion=31 -fimf-accuracy-bits=11 -no-prec-div -no-prec-sqrt -fno-alias -vec-report2 -o MonteCarloStep5.o MonteCarloStep5.cpp
icpc: remark #10346: optimization reporting will be enabled at link time when performing interprocedural optimizations
icpc Driver.o MonteCarloStep5.o -o MonteCarlo -L /opt/intel/composerxe/mkl/lib/intel64 -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -openmp
Driver.cpp(130): (col. 5) remark: loop was not vectorized: existence of vector dependence.
Driver.cpp(161): (col. 9) remark: loop was not vectorized: existence of vector dependence.
MonteCarloStep5.cpp(69): (col. 9) remark: SIMD LOOP WAS VECTORIZED.
MonteCarloStep5.cpp(58): (col. 5) remark: loop was not vectorized: not inner loop.
[sli@localhost step5]$ setenv KMP_AFFINITY "compact,granularity=fine" 
[sli@localhost step5]$ ./MonteCarlo
Monte Carlo European Option Pricing in Single Precision
Pricing 1998848 options with path length of 262144.
Completed in  43.9896 seconds.
Computation rate - 45439.090 options per second.

在双插槽 Sandy Bridge 系统 (2.6 GHz)上,我们实现了 15.20 倍的性能提升。 在配备 16 个内核、32 条线程的执行环境中,15.20 倍的性能提升意味着每条线程运行的大部分代码都采用了并行化的形式。 请注意,我们必须将线程的相似性设置成紧凑模式,以最大限度地利用共享高速缓存,同时将粒度设置为"精细"。 相对于使用 GCC 4.4.6 建立的基准,系统性能提升 1324.6 倍。

以下总结了每种技术所实现的性能改进。 测量应用这些技术前后的性能,并计算出比率。 例如,以使用英特尔编译器作为基准,添加英特尔 MKL 将使性能提升 5.53 倍。

3.6 第 5 步: 从英特尔® 多核扩展到英特尔® 集成众核架构

基于集成众核架构的英特尔至强融核协处理器可为 Monte Carlo 等高度并行化应用带来立竿见影的优势。 基于英特尔 MIC 架构的第一款协处理器拥有多达 61 个内核和 244 条线程,有助于进一步提升可扩展软件的性能。 而且,在这 61 个内核中,每个都拥有 512 位的 SIMD 引擎,该引擎可以处理 8 种 64 位数据、或 16 种 32 位数据整数或浮点数。

3.6.1 重建英特尔 MIC 架构

英特尔 C++ Composer XE 2013 的一个最重要的特性是集成了英特尔多核与英特尔 MIC 开发环境。 您需要安装 SKU。 您可以编译和运行原生多核应用;或者编译相同的源代码以构建面向英特尔至强融核协处理器的可执行程序。

按照分布式优化框架获得软件后,您可以对该软件进行重新编译,创建一个可以在本地协处理器上运行的二进制文件。


icpc -c -g -O3 -ipo -openmp -mmic -fimf-precision=low -fimf-domain-exclusion=31 -fimf-accuracy-bits=11 -no-prec-div -no-prec-sqrt -fno-alias -vec-report2 -opt-threads-per-core=4 -o Driver.o Driver.cpp
icpc: remark #10346: optimization reporting will be enabled at link time when performing interprocedural optimizations
icpc -c -g -O3 -ipo -openmp -mmic -fimf-precision=low -fimf-domain-exclusion=31 -fimf-accuracy-bits=11 -no-prec-div -no-prec-sqrt -fno-alias -vec-report2 -opt-threads-per-core=4 -o MonteCarloStep5.o MonteCarloStep5.cpp
icpc: remark #10346: optimization reporting will be enabled at link time when performing interprocedural optimizations
icpc -mmic Driver.o MonteCarloStep5.o -o MonteCarlo -L /opt/intel/composerxe/mkl/lib/mic -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -openmp
Driver.cpp(130): (col. 5) remark: loop was not vectorized: existence of vector dependence.
Driver.cpp(161): (col. 9) remark: loop was not vectorized: existence of vector dependence.
MonteCarloStep5.cpp(69): (col. 9) remark: SIMD LOOP WAS VECTORIZED.
MonteCarloStep5.cpp(58): (col. 5) remark: loop was not vectorized: not inner loop.

为了运行由之前的生成文件创建的二进制 Monte Carlo,我们必须使用 scp 将它从主机系统拷贝到卡上。 同时还要拷贝三个英特尔 MKL 共享库。 完成之后,我们必须使用 ssh 打开设备操作系统命令窗口。 必须对部分环境变量进行设置,以确保正确的 OpenMP 运行时行为。


[sli@cthor-knc14 step5]$ scp MonteCarlo cthor-knc14-mic1:
MonteCarlo                                                  100%  493KB 492.8KB/s   00:00
[sli@cthor-knc14 step5]$ ssh cthor-knc14-mic1
~ $ export LD_LIBRARY_PATH=.
~ $ export KMP_AFFINITY="compact,granularity=fine"
~ $ ./MonteCarlo
Monte Carlo European Option Pricing in Single Precision
Pricing 1998848 options with path length of 262144.
Completed in   4.9851 seconds.
Computation rate - 400964.798 options per second.

以前在主机系统上需要花 44 秒才能完成的程序,现在只需不到 5 秒的时间,实现了 8.82 倍的速度提升。 在本节中,我们研究了各种可以进一步提升应用性能的方法。

3.6.2 目标架构优化

不同的英特尔架构实施方式可以创建包含独有性能配置文件的 ISA,跨 ISA 移动软件时我们可以使用该配置文件。 在协处理器上,英特尔实施了一个名为“可扩展数学单元”的硬件函数近似单元。 四个初等函数(log2()、exp2()、rsqrt() 和 rcp())直接在具备 1-2 周期吞吐率和 4-8 周期延迟的硬件中实施。 快捷的初等函数可以帮助使用初等函数和简单转换构建的四种其他函数加快速度。

例如,底数 2 的指数和底数 e 的指数在以下称为换底公式的数学等式中关联起来。


 ex = 2x*log2E
ln(x) = log2(x)/log2E  = log2(x)*ln2 = log2(x)*M_LN2
exp(x)= exp2(x * log2E)= exp2(x * M_LOG2E) 


请注意 LOG2E,e 的底数 2 对数是在数学 .h 中定义为 M_LOG2E 的常数。

初等函数

延迟

吞吐量

exp2()

8

2

log2()

4

1

rcp()

4

1

rsqrt()

4

1

导函数

延迟

吞吐量

pow()

16

4

sqrt()

8

2

div()

8

2

ln()

8

2

同样,pow() 可以由 exp2()、乘法和 log2() 组成;sqrt() 由 rsqrt() 和 rcp() 组成;div() 由 rcp() 和乘法组成。

这有助于我们对 MonteCarlo 代码进行调整,以便调用以 2 为底数的函数,而不是以自然数为底数的函数。


float VBySqrtT = VLOG2E * sqrtf(T[opt]);float MuByT = RVVLOG2E * T[opt];
        float Sval = S[opt];
        float Xval = X[opt];
        float val = 0.0, val2 = 0.0;

#pragma vector aligned
#pragma simd reduction(+:val) reduction(+:val2)
#pragma unroll(4) 
        for(int pos = 0; pos < BLOCKSIZE; pos++)
        {
            float callValue  = Sval * exp2f(MuByT + VBySqrtT * random[pos]) - Xval;
            callValue = (callValue > 0) ? callValue : 0;
            val  += callValue;
            val2 += callValue * callValue;
        }

在本案例中,在 exp2f() 中有两种表达式,我们必须同时对它们进行调整,同时提高内层循环之外的调整。

您可能认为我们需要做的就是使用 exp2(x*M_LOG2E) 和 LOG2(x)*M_LN2 (或 LOG2(x)/M_LOG2E) 对 exp(x) 和 log(x) 进行替换。 但是,由于存在其他乘法,仅做这些还不足以获得 EMU 函数的所有优势。 在大多数情况下,当可以提前调整内层循环之外的多重时,乘法的成本将消失,而使用 EMU 函数的性能优势就会出现。

3.6.3 内存访问模式优化

在我们当前的 Monte Carlo 实施中,所有线程共享了提前生成的 256K 随机数,这要求单精度浮点占用 1 MB 的空间。 这些随机数与工作线程的二级高级缓存不匹配。 每个线程只能将一部分数据带入到高速缓存中。 当它完成了当前的部分时,它会生成二级高速缓存缺失事件,该事件将触发环形互连进而将更多的数据带入到二级高级缓存中。 如果在访问这些数据时线程不一致,每个线程可能会访问这 1 MB 数据的任意部分,此时环形互连会因为太过繁忙而不能满足来自所有 244 条线程的 二级高速缓存缺失。 从本质上来说,该环形互连的带宽成为了瓶颈,即使此时只包含了 1 MB 的数据。 不一致的内存访问使该 1 MB 数据被访问了 244 次,创建了等同于 244 MB 的带宽需求。

而一致的访问将确保在第一条线程请求数据时数据一次就能从内存中抓取出来。 然后,数据将越过环形互连并寄存在第一条线程的 二级高速缓存中。 当其他线程需要数据时,它们就能从第一条线程的 二级高级缓存中获取数据。 每条线程只需停留在二级高级缓存中就能访问所有的数据。 由于内存访问是一致的,1 MB 的数据只需跨越环形互连一次,因此减轻了内存带宽的压力。

为了最大限度地利用二级高级缓存的容量,我们需要计算有多少数据可以匹配每条工作线程的二级高级缓存。 每个内核拥有 512 K 的高级缓存和 4 条线程。 每条线程平均拥有 128 K。 减去运行时需求,只有 64 K 可供应用程序控制。 64 KB 是 16 K SP 和 8 K DP 浮点数。

我们可以对 Monte Carlo 的结构进行调整,使其可以一次处理 16 K 的随机数据块,然后当所有线程完成时可以在下一个数据块中处理更多任务。 我们必须在内存中保存部分中间结果。 当将所有数据块处理完成后,我们将处理中间结果,并将其转变为最终结果。


#include "MonteCarlo.h"
#include "math.h"
#include "mkl_vsl.h"
#include "omp.h"
#define RANDSEED 123

static const float  RVVLOG2E = (RISKFREE-0.5f*VOLATILITY*VOLATILITY)*M_LOG2E;
static const float  INV_RAND_N = 1.0f/RAND_N;
static const float  F_RAND_N = static_cast<float>(RAND_N);
static const float STDDEV_DENOM = 1 / (F_RAND_N * (F_RAND_N - 1.0f));
static const float CONFIDENCE_DENOM = 1 / sqrtf(F_RAND_N);
static const int BLOCKSIZE = 32*1024;
static const float  RLOG2E = RISKFREE*M_LOG2E;
static const float  VLOG2E = VOLATILITY*M_LOG2E;

void MonteCarlo(
    float *h_CallResult,
    float *h_CallConfidence,
    float *S,
    float *X,
    float *T
)
{
    __attribute__((align(4096))) float random [BLOCKSIZE];
    VSLStreamStatePtr Randomstream;
    vslNewStream(&Randomstream, VSL_BRNG_MT19937, RANDSEED);
#ifdef _OPENMP
    kmp_set_defaults("KMP_AFFINITY=compact,granularity=fine");
#endif
#pragma omp parallel for
    for(int opt = 0; opt < OPT_N; opt++)
    {
        h_CallResult[opt]     = 0.0f;
        h_CallConfidence[opt] = 0.0f;
    }

    const int nblocks = RAND_N/BLOCKSIZE;for(int block = 0; block < nblocks; ++block){
        vsRngGaussian (VSL_METHOD_SGAUSSIAN_ICDF, Randomstream, BLOCKSIZE, random, 0.0f, 1.0f);
#pragma omp parallel for
        for(int opt = 0; opt < OPT_N; opt++) 
        {
           float VBySqrtT = VLOG2E * sqrtf(T[opt]);
           float MuByT = RVVLOG2E * T[opt];
           float Sval = S[opt];
           float Xval = X[opt];
           float val = 0.0, val2 = 0.0;

#pragma vector aligned
#pragma simd reduction(+:val) reduction(+:val2)
#pragma unroll(4) 
           for(int pos = 0; pos < BLOCKSIZE; pos++)
           {
              float callValue  = Sval * exp2f(MuByT + VBySqrtT * random[pos]) - Xval;
              callValue = (callValue > 0) ? callValue : 0;
              val  += callValue;
              val2 += callValue * callValue;
           }

           h_CallResult[opt]     +=  val;h_CallConfidence[opt] +=  val2;
        }
    }
#pragma omp parallel for
    for(int opt = 0; opt < OPT_N; opt++)
    {
        const float val      = h_CallResult[opt];
        const float val2     = h_CallConfidence[opt];
        const float  exprt    = exp2f(-RLOG2E*T[opt]);
        h_CallResult[opt]     = exprt * val * INV_RAND_N;
        const float  stdDev   = sqrtf((F_RAND_N * val2 - val * val) * STDDEV_DENOM);
        h_CallConfidence[opt] = (float)(exprt * stdDev * CONFIDENCE_DENOM);
    }
    vslDeleteStream(&Randomstream);
}

减缓带宽压力将使 Monte Carlo 成为真正与计算绑定的工作负载。 它可将总运行时节省一秒,同时将性能提升 1.24 倍。


~ $ ./MonteCarlo
Monte Carlo European Option Pricing in Single Precision
Pricing 1998848 options with path length of 262144.
Completed in   4.0256 seconds.
Computation rate - 496530.481 options per second.

内存访问模式优化在英特尔至强处理器和英特尔至强融核协处理器上都可以实施。 就 IA 的两种风格而言,我们针对每线程可用 L2 容量的计算恰好是一样的。 因此,我们希望通过编译英特尔 MIC 优化程序使其可以运行在英特尔多核处理器上。 如您所见,它的运行速度相比以往提升了 1.11 倍。


[sli@localhost .solution]$ make -B
icpc -c -O3 -openmp -xAVX -fimf-precision=low -fimf-domain-exclusion=31 -fimf-accuracy-bits=11 -no-prec-div -no-prec-sqrt -fno-alias -vec-report2  -o Driver.o Driver.cpp
Driver.cpp(130): (col. 5) remark: loop was not vectorized: existence of vector dependence.
Driver.cpp(161): (col. 9) remark: LOOP WAS VECTORIZED.
icpc -c -O3 -openmp -xAVX -fimf-precision=low -fimf-domain-exclusion=31 -fimf-accuracy-bits=11 -no-prec-div -no-prec-sqrt -fno-alias -vec-report2  -o MonteCarloStep5.o MonteCarloStep5.cpp
MonteCarloStep5.cpp(60): (col. 5) remark: loop was not vectorized: existence of vector dependence.
MonteCarloStep5.cpp(53): (col. 5) remark: LOOP WAS VECTORIZED.
MonteCarloStep5.cpp(75): (col. 12) remark: SIMD LOOP WAS VECTORIZED.
MonteCarloStep5.cpp(64): (col. 5) remark: loop was not vectorized: not inner loop.
MonteCarloStep5.cpp(88): (col. 5) remark: LOOP WAS VECTORIZED.
icpc -xAVX -openmp Driver.o MonteCarloStep5.o -o MonteCarlo -L /opt/intel/composerxe/mkl/lib/intel64 -lmkl_intel_lp64 -lmkl_sequential -lmkl_core
[sli@localhost .solution]$ ./MonteCarlo
Monte Carlo European Option Pricing in Single Precision
Pricing 1998848 options with path length of 262144.
Completed in  39.6060 seconds.
Computation rate - 50468.304 options per second.

4. 结论

本文简单介绍了 Monte Carlo 方法及其在衍生品定价应用中的应用。 此外,我们还介绍了最新版的英特尔® 至强融核™ 协处理器和集成软件开发环境英特尔 C++ Composer XE 2013。 我们衍生出 Monte Carlo 欧式看涨期权定价的 C/C++ 实施,并将其用作我们的案例研究,同时还展示了分布式优化框架。

按照使用英特尔® C++ Composer XE 2013 的性能优化框架完成五个步骤之后,Monte Carlo 欧式期权定价在英特尔至强融核协处理器上实现的性能提升超过了 1,471 倍。 我们还展示了您可以重新编译和重构经过英特尔至强处理器优化的应用程序,使其在本地英特尔至强融核协处理器上运行。 您可以在英特尔至强融核协处理器上进一步优化该应用,并解决可扩展性问题。 生成的代码仍然可以重新编译,最终可以在英特尔至强处理器上运行并获得更高的性能。

事实证明,分布式优化框架不仅适用于金融数字应用,而且对于一般科学计算来说也同样有效。 借助英特尔 C++ Composer XE 2013,科学领域的编程人员能够更容易地以结构化活动形式进行性能优化,同时快速了解程序的局限性,并尽可能地实现各个层面上的并行化。

其他资源

  • 面向 Linux* 的英特尔® Composer XE 2013,包括英特尔® MIC 架构**
  • 英特尔® C++ 编译器 XE 12.0 用户和参考指南**
  • 英特尔® MIC 架构优化指南***
  • C++ 技术报告 1

** 本文假设安装了英特尔® Composer。

*** 本文在 MIC 开发人员网站上已有提供:https://mic-dev.intel.com。

关于作者

Shuo Li 任职于英特尔公司的软件和服务事业部。 他在软件开发领域拥有 24 年的丰富经验。 其关注的重点包括并行编程、计算金融和应用性能优化。 作为软件性能工程师,Shuo 最近与金融服务业的软件开发人员与建模人员密切合作,帮助其在英特尔平台上实现了最佳的性能。 Shuo 拥有俄勒冈大学计算科学专业硕士学位以及杜克大学 MBA 学位。

有关编译器优化的更完整信息,请参阅优化通知