SODECL：一个并行计算随机微分方程组多轨道的开源库 SODECL公司

2.2.1 集成方法。SODECL使用Euler-Maruyama方法[Higham2001; 斋藤和三井1996]在方程式中整合方程式(1). Euler-Maruyama近似$\lbrace\mathbf{X}^k=（X_1^k，\ldots，X_n^k）：k\ge 1\rbrace$找到真正的解决方案$\lbrace\mathbf{X}（t）：\mathbf}（t_0）=\mathbf{十} _0（0），t\ge t_0\r痕迹$由递归定义

2.2.2 噪声发生器。生成噪波变量$N_j（0,1）$用于方程式的每个步骤(2)，SODECL使用Random123库[Salmon等人。2011]. Random123生成均匀分布的随机数，然后使用Box-Muller算法将其转换为正态分布的随机数字[Box和Muller1958]:

2.2.3 OpenCL代码生成。SODECL在运行时创建并构建OpenCL源代码字符串。这允许用户更改f_i美元$和$g_{ij}$方程式中的函数(1)无需编译源代码。此外，SODECL允许实现包装器，例如，对于MATLAB或Python，使用运行SODECL可执行文件的脚本/函数，或使用MATLAB.mex文件而不重新编译。形成OpenCL源代码字符串的过程通过下面概述的四个步骤执行。

第1步。最初，SODECL将五个命名常量的定义附加到空OpenCL源代码字符串的开头（参见表1). 这些常数的值取决于SDE系统的一般属性（方程数、噪声项和参数）、积分器的时间步长以及积分器在每个OpenCL内核调用中要执行的步骤数。这些值由用户传递给SODECL库。OpenCL源代码字符串的其余部分是使用位于三个单独文件中的代码生成的，如下面的步骤2-4所述。

第2步第一个文件，积分器调用器.cl，包含主机调用以在OpenCL设备上运行的内核函数。此功能如图所示1。内核函数的参数是指针，显示数组在设备全局内存中的位置。这些数组包含SDE系统的时间点，用于近似求解，SDE系统因变量的值以及要计算不同轨道的参数组合。另一个计数器数组被传递给内核，供Random123库使用。内核生成SDE系统集成所需的独立、正态分布的数字。随后，相同的内核调用集成SDE系统的函数一个时间步长（参见步骤3）。调用总数由常量定义_numsteps（步数）_.已集成SDE系统_numsteps（步数）_次，内核将最后的集成状态值复制到OpenCL设备的全局内存中，以便在下一个内核步骤中使用并由主机访问。

图1。 SODECL内核函数的定义。函数参数是指针，显示OpenCL设备全局内存中数组的位置。这些数组是SDE系统的解决方案的近似时间点（通过t吨)，因变量的相应值（通过访问年)，以及系统参数的组合（由访问参数_g). 一个额外的计数器数组（由访问计数器（_g）)传递给内核供Random123库使用。

步骤3。创建OpenCL源时使用的第二个文件名为随机euler.cl。它包含函数系统集成器（如图所示2)，它实现了Euler-Muryama集成方法。系统集成器它调用计算SDE系统确定性和随机性组件的函数（参见步骤4）。一旦对这些进行了评估，随机euler.cl在下一个时间步计算所有因变量的值。

图2。 SODECL的定义系统集成器功能。功能参数是SDE系统的当前时间值(t吨)，因变量的对应值(年)，参数值(第页)，当前噪声值(噪音)，以及通过Euler-Maruyama方法计算的因变量的新值(你).

第4步。第三个文件包含两个名为草皮系统和sode_system_stoch公司计算漂移和扩散项${f_i}$和$g_{ij}$在数值格式中(2)分别是。它们的定义如图所示三，并且必须由用户指定（图中显示了一个示例5).

图3。 指定SDE系统的SODECL功能的定义。功能参数为当前时间值(t吨)，SDE系统因变量的对应值(年)SDE确定性分量因变量的新值(你)、模型参数(第页)、噪声值(噪音)和SDE随机分量的新值(斯托克).

2.2.4 SDE积分器执行算法。SODECL集成执行算法的逻辑流主要是一个循环。在每次循环迭代期间，所有轨道的因变量的相应值都会复制到主机设备，并且复制的值会并行保存到数组中，为下一次内核调用做准备。调用OpenCL内核来集成模型_numsteps（步数）_长度步数_m_dt（m_dt）_，这意味着集成的输出每(_numsteps（步数）_ $\次$ _m_dt（m_dt）_)积分时间的秒数。这允许以以下频率存储数据$1/$(_numsteps（步数）_ $\次$ _m_dt（最小值）_)赫兹。一旦完成所有参数组合在所需时间跨度内的积分，算法结束。循环迭代次数，千美元$，使用以下公式计算

与CPU相比，结果的存储频率在GPU提供的加速中起着重要作用。这是由于与GPU内存之间的传输造成的瓶颈，以及将数据存储在主机内存中以供进一步分析所需的时间。低频（250–1000 Hz）可避免主机内存中的大量数据存储。此外，模型大小也是一个重要因素，因为它设置了要保存在计算机RAM上的模型集成数量。此外，返回的数据点数量越多，分析时间就会越长；然而，数量少可能不足以进行准确的分析[阿夫拉米迪斯和阿克曼2017].

2.2.5 Python接口。SODECL的Python包接口函数如图所示4，而表2描述了函数传递给SODECL可执行文件的参数。Python库pybind11(https://github.com/pybind/pybind11)用于实现接口。SODECL还提供了一个MATLAB函数接口，但需要注意的是，这目前是实验性的。

图4。 SODECL Python包装函数的定义。该函数用于从Python中执行SODECL可执行文件。表中定义了每个参数2.

2.3.1 协议。为了评估SODECL的执行速度，我们使用了成熟的Kuramoto模型的随机变量来描述弱耦合相位振荡器种群的动力学[Acebrón等人。2005; Bick等人。2011; 大岛1992; Nakao公司2016; 斯特罗加茨2000]. Kuramoto模型通常用于模拟生物振荡器的同步[Acebrón等人。2005; Bick等人。2011]并且实现和定制相对简单（例如，振荡器种群大小可以很容易地改变）。我们在这里使用的SDE Kuramoto系统具有以下形式

图5。 随机Kuramoto模型的SODECL实现（方程(5)). 第一个函数草皮系统用于计算模型的确定性组件，而第二个函数sode_system_stoch公司用于计算随机部分。

我们比较了使用不同多核和多核计算设备时SODECL的执行速度与仅使用HPC CPU的一个核时的执行速度。表中列出了用于评估SODECL性能的计算设备三。ID为I4790K的CPU是高端台式机CPU，而2X6142是具有两个多核CPU的HPC节点。W8100和P100是工作站/服务器GPU，具有科学计算所需的高双精度计算能力。为了更好地理解SODECL运行时如何受到系统大小的影响，我们改变了Kuramoto模型中方程的数量（即，我们改变振荡器的数量亿美元$网络中）。此外，我们还研究了轨道数百万美元$集成会影响SODECL在不同计算设备上的性能。针对每个选项执行八次独立的求解器运行N美元$和百万美元$.

我们还将SODECL在桌面CPU上的执行速度与在同一CPU上运行的等效并行MATLAB实现进行了比较。为了进一步加速此实现，如图所示6，我们使用MATLAB编码器将MATLAB函数转换为C++代码。MATLAB的并行化是通过OpenMP（Open Multi-Processing）自动执行的，OpenMP是一个应用程序编程接口，用于使用CPU和协处理器内核的并行应用程序[Dagum和Menon1998].

图6。 MATLAB函数用于并行计算随机Kuramoto模型的轨道。第一个功能kuramoto平行调用第二个函数库拉莫托系统-它计算模型的一个轨道-多次。

2.3.2 结果。图中显示了使用SODECL与不同计算设备、轨道数和系统大小获得的比较执行速度和加速比7（见表S1（第一阶段）–第6页对应数值）。对于每个计算设备，加速比是相对于2X6142节点的一个内核计算的（即，2X6142-单核运行时与计算设备运行时的比率）。用于每个设备和系统大小的本地组大小(N美元$)如表所示4。本地组大小是将并行运行并可以通信的工作项的数量。在令人尴尬的并行应用程序的情况下，尽管没有通信，但组大小可能会对运行时产生影响。显示的值是通过探索性测试发现的，因为据我们所知，没有计算最佳局部组大小的分析公式。

图7。 随机Kuramoto模型的运行时间（秒）和加速比，作为不同系统尺寸轨道数的函数(N美元$)，当使用不同的并行计算机硬件实现时。图A显示了$N=5$，图C$N=10$，并为绘制EN=15美元$。这些值是八次运行的平均值。图B、D和F显示了相应的加速比$N=5$,$N=10$、和N=15美元$分别相对于2X6142的一个逻辑处理器。表中列出了每种情况下使用的硬件三.

对于测试的CPU，结果表明，当使用2X6142节点的所有32个内核时，SODECL介于$\大约$1.12和大约$\$在不同的系统尺寸和轨道数下，比同一节点的一个内核快6.65倍。相比之下，当使用I4790K CPU的所有内核并将其速度与2X6142的一个内核进行比较时，我们观察到从大约$\$0.69至大约$\$2.66随着系统规模从5增加到15，轨道数从512增加到163840。2X6142总体上获得了更高的加速比，这是因为其内核数量更多，体系结构更新。

对于测试的GPU，SODECL与使用一个2X6142内核相比，其加速比在0.6-4.5范围内，这取决于Kuramoto模型的大小和轨道数。更多的模型方程和轨道允许更多地使用GPU资源，这意味着提高了速度。当将GPU运行时间与使用2X6142的所有内核时获得的运行时间进行比较时，可以看出GPU通常较慢。例如，这可能是由于较高的初始化时间开销造成的。

最后，图8比较了I4790K CPU上Kuramoto模型求解器的SODECL和MATLAB实现所获得的运行时间（相应的数值在表中给出第7部分–第9部分). 对于最小的系统尺寸和最少的轨道数，MATLAB实现只比SODECL快（见图8（B） （E）和（H））。在所有其他情况下，SODECL$\大约$速度提高5.8倍。

图8。 不同系统尺寸的随机Kuramoto模型的运行时间作为轨道数的函数(N美元$)，在I4790K CPU的所有逻辑处理器上使用MATLAB和SODECL实现。运行时间以秒为单位。图A和图B显示了$N=5$、D和E代表$N=10$，G和H代表$N=15美元$图B、E和H显示了前5000个轨道的运行时间。图C、图F和图I显示了MATLAB和SODECL实现的相对速度$N=5，\，10$和15个。

2.4.1 协议。为了评估SODECL的准确性和数值稳定性，我们使用N=100美元$,$\omega_i\sim\mathcal{U}（0.2,0.4）$,$\phi_i（0）\sim\mathcal{U}（-\pi，\pi）$、和$p_i\sim\mathcal{U}（0.01,0.03）$对于两个值$K_1=0.02$和$K_2=0.2$整体耦合强度的千美元$这些值被选在连续极限下确定性系统中发生的Kuramoto跃迁的任一侧($N\rightarrow\infty（右箭头）$)在$K=K_c=\压裂{2}{5\pi}\近似值$0.1273. 作为千美元$通过这个临界值增加，宏观相互夹带（MME）-一个共同频率的集体同步节律-在振荡器群中观察到[Daido1992; 斯特罗加茨2000]. 对于的两个值千美元$，使用以下积分时间步长，在400 s内进行了64次独立积分（平均周期20.94 s的约19个周期）：$\增量t=\压裂{2^{l-5}}{5}，1\le l\le 5$对于每个实现，使用复数阶参数[Nakao量化模拟振荡器群的同步程度2016; 斯特罗加茨2000]定义如下：

2.4.2 结果。图9绘制平均相位相干性$\langle r（t）\rangle美元$相位相干性的标准差$\西格玛\左（r（t）\右）$随着时间的变化千美元$值，当$\增量t=0.05$。如预期$K=K_1 \lt K_c$,$\langle r（t）\rangle美元$显示接近0的小振幅振荡，而对于$K=K_2\gt K_c$,$\langle r（t）\rangle美元$渐近线到接近1的值。Kymograph图显示了每一个生成的典型实现的种群阶段的时间演变千美元$值绘制在图中10这些与相位相干图一致-振荡器在整个积分间隔内失同步$K=K_1$，而对于$K=k2$，振荡器在大约50秒后同步，然后保持同步。最后，图11绘制的值$\langle r（t）\rangle美元$和$\西格玛\左（r（t）\右）$积分间隔结束时$t_{MAX}=400$随集成时间步长变化$\增量t$可以看出，这两种情况千美元$值，$\langle r（t_{MAX}）范围$和$\西格玛\左（r（t_{MAX}）\右）$表现出对$\增量t$表明积分算法在数值上是稳定的[Kloeden和Platen2011].

图9。 相位相干的时间演化美元（吨）$在随机Kuramoto模型中（方程(5))何时$N=100美元$对于K美元=0.02$（黑线）和K美元=0.2$（红线）。$\langle r（t）\rangle美元$和$\西格玛\左（r（t）\右）$表示的平均值和标准偏差美元（吨）$分别根据64个独立实现的方程式计算(6)具有集成时间步长$\增量t=0.05$.振荡器频率分布$\欧米茄_i$，初始阶段$\θ_i（0）$和噪声强度$p_i$如文中所述。

图10。 振荡器相位的时间演变$θ_i（t）$在随机Kuramoto模型中。每个曲线图对应于方程式的实现之一(6)用于计算图中所示的相位相干统计9图A和图B显示了以下方面的模拟K美元=0.02$和K美元=0.2$分别是。

图11。 对积分时间步长的依赖性$\增量t$相位相干统计$\langle r（t）\rangle美元$和$\西格玛\左（r（t）\右）$积分间隔结束时$t_{MAX}=400$。黑线和轴显示的结果K美元=0.02$; 红线和轴显示的结果K美元=0.2$在每种情况下，$\langle r（t_{MAX}）范围$和$\西格玛\左（r（t_{MAX}）\右）$由64个独立实现的方程式计算得出(6)带有N=100美元$.振荡器频率分布$\欧米茄_i$，初始阶段$\θ_i（0）$和噪声强度$p_i$如文中所述。注意上的对数刻度x-轴。