并行蒙特卡洛算法 #1
如何充分利用当前技术处理计算密集型应用程序？
感谢您提名我为月度文章。
引言
这是我希望撰写的关于计算密集型应用程序的并行化和分布式处理系列文章的第一篇。
在本文中，我将介绍 .NET 开发人员在单台机器上并行化蒙特卡洛例程可用的几种技术。我们发现，虽然 TPL 有其优势，但编写能超越线程池和混合 C#/C++ 实现的高效例程也很困难。
背景
简而言之，蒙特卡洛例程是任何迭代例程，它使用随机数生成来计算一个结果，该结果的估计精度随着迭代次数的增加而提高。此外，由于每次迭代结果都是独立的，我们可以将计算分散到任意数量的位置，然后再将结果聚合。
不明白？那请点击这里 获取一个合适的描述。
在本文中，我使用了 Press、Flannery、Teukolsky 和 Vetterling 在《Numerical Recipes in C》第 239-240 页中提出的算法版本。该例程旨在找到具有不同密度的圆环截面的重量和质心。通俗地说：我们正在切割一个有许多凹凸的甜甜圈，并找到单个块的质心。
Using the Code
附加的解决方案应该可以在 Visual Studio Express for Windows Desktop 2012 中编译和运行，您应该在 Release 模式下构建。它包含 4 个项目
Montecarlo- 仅包含 .NET 代码ComponentWin32- 一个具有静态导出的本机 C++ DLLComponentCLR- 一个 C++/CLR 程序集，可以被任何其他 .NET 程序集引用CppConsole- 一个简单的本机蒙特卡洛实现运行器
关注点
演示了几种技术，这些技术是我在我的 4 核台式机上开发和测试的
- 纯 C# 单线程，我们将用它作为准确性的基准
- 本机 C++ 单线程作为速度基准
- 混合 C#/C++ 版本，用于演示一种集成方法
- 单线程 PInvoke 版本，我们创建了一个非托管对象并调用其方法
- 一个
Threadpool版本，其中工作被分割，由池化线程执行，然后进行聚合 - TPL
Parallel.For方法，通过每次演进来提高性能和结果的准确性 - 线程池化和任务化的 PInvoke，我们尝试使用最快的执行方法配合 .NET 并行化。
单线程 C#
从 C# 的单线程方法开始，我们可以快速实现我们的基本例程，并用它来验证任何后续演进的结果。
public void Calculate(ref ResultsdotNet results)
{
double w,x,y,s,z,dw,dx,dy,dz,sw,swx,swy,swz, varw,varx,vary,varz,ss,vol;
w=x=y=s=z=dw=dx=dy=dz=sw=swx=swy=swz=varw=varx=vary=varz= 0.0;
ss=0.2*(Math.Exp(5.0)-Math.Exp(-5.0));
vol=3.0*7.0*ss;
long count = 0;
long n = results.Iterations;
for(int j=1; j<= n; j++)
{
x=1.0+3.0* r.NextDouble();
y=(-3.0)+7.0* r.NextDouble();
s = 0.00135 + ss * r.NextDouble();
z=0.2*Math.Log(5.0*s);
double b = Math.Sqrt((x * x ) + (y * y) ) - 3.0;
double a = (z * z) + ( b * b);
if (a < 1.0)
{
sw += 1.0;
swx += x;
swy += y;
swz += z;
varw = 1.0;
varx += x * x;
vary += y * y;
varz += z * z;
}
count++;
}
results.W = vol * sw / n;
results.X = vol * swx/n;
results.Y = vol * swy/n;
results.Z = vol * swz/n;
results.dW = vol * Math.Sqrt((varw / n - Math.Pow((sw /n), 2.0)) / n);
results.dX = vol * Math.Sqrt((varx / n - Math.Pow((swx / n), 2.0)) / n);
results.dY = vol * Math.Sqrt((vary / n - Math.Pow((swy / n), 2.0)) / n);
results.dZ = vol * Math.Sqrt((varz / n - Math.Pow((swz / n), 2.0)) / n);
}
在 Release 模式下运行 1000 万次迭代得到
我们知道可以做得更快，但只要其他方法给我们类似的数字，我们就知道我们走在正确的轨道上。
单线程本机 C++
让我们在 C++ 中尝试同样的事情，看看结果如何。
vector<double> TorusMontecarlo::Calculate(int n)
{
double w,x,y,s,z,dw,dx,dy,dz,sw,swx,swy,swz, varw,varx,vary,varz,ss,vol;
w=x=y=s=z=dw=dx=dy=dz=sw=swx=swy=swz=varw=varx=vary=varz= 0.0;
ss=0.2*(exp(5.0)-exp(-5.0));
vol=3.0*7.0*ss;
for(int j=1; j<= n; j++)
{
x=1.0+3.0* RandomNumbers::ran3(&idum);
y=(-3.0)+7.0* RandomNumbers::ran3(&idum);
s=0.00135+ss* RandomNumbers::ran3(&idum);
z=0.2*log(5.0*s);
if (SQR(z) + SQR( sqrt(SQR(x)+SQR(y)) -3.0) < 1.0)
{
sw += 1.0;
swx+= x;
swy+= y;
swz += z;
varw = 1.0;
varx += SQR(x);
vary += SQR(y);
varz += SQR(z);
}
}
w = vol * sw/n;
x = vol * swx/n;
y = vol * swy/n;
z = vol * swz/n;
dw=vol*sqrt((varw/n-SQR(sw/n))/n);
dx=vol*sqrt((varx/n-SQR(swx/n))/n);
dy=vol*sqrt((vary/n-SQR(swy/n))/n);
dz=vol*sqrt((varz/n-SQR(swz/n))/n);
double res [] = {w,dw,x,dx,y,dy,z,dz};
vector<double> result (res, res + sizeof(res)/ sizeof(double) );
return result;
}
在没有太多关于如何优化 MS C++ 编译器知识的情况下，并使用一个相当原始的秒表实现，我们获得了相似的准确性，但耗时约为 2/3。
显然，C++ 不是每个人都喜欢，并且开发完整的应用程序而不是仅仅库可能不适合业务应用程序。所以让我们看看在混合 C++ 和托管世界时我们表现如何。
混合 C# 和 C++
让我们在编译器中打开 /clr 开关，添加一些支持类，并从 C# 中调用此例程（实现在 ComponentCLR.dll 中）。
效果不是特别惊人，但仍然优于纯 C# 实现。
如果我们通过 PInvoke 公开我们的本机 Win32 C++ 实现并从中静态调用它，会怎么样？
我们的非托管对象表现与我们的本机实现一样好，那么何必费心使用托管 C++ 呢？
C# 线程池
从 .NET 1.1 开始，C# 开发人员就可以利用线程池并行处理工作，而无需手动实例化线程。但要做到最优，我们需要弄清楚特定机器上的核心数量，并将要完成的工作平均分配给它们。一旦工作被发送到线程池并启动，父线程就需要等待直到每个批次发出完成信号。
下面这样的代码应该可以工作
public ResultsdotNet CalculateThreadpool(ResultsdotNet total)
{
ManualResetEvent[] events = new ManualResetEvent[Environment.ProcessorCount];
List<ResultsdotNet> results = new List<ResultsdotNet>();
for (int i = 1; i <= Environment.ProcessorCount; i++)
{
ResultsdotNet batchResult = new ResultsdotNet()
{ Iterations = total.Iterations / Environment.ProcessorCount,
ManualResetEvent = new ManualResetEvent(false) };
results.Add(batchResult);
events[i - 1] = batchResult.ManualResetEvent;
ThreadPool.QueueUserWorkItem(new WaitCallback(CalculateCallback), batchResult);
}
WaitHandle.WaitAll(events);
foreach (ResultsdotNet batchResult in results)
{
total.dW += batchResult.dW;
total.dX += batchResult.dX;
total.dY += batchResult.dY;
total.dZ += batchResult.dZ;
total.W += batchResult.W;
total.X += batchResult.X;
total.Y += batchResult.Y;
total.Z += batchResult.Z;
total.Iterations += batchResult.Iterations;
}
total.dW = total.dW / results.Count;
total.dX = total.dX / results.Count;
total.dY = total.dY / results.Count;
total.dZ = total.dZ / results.Count;
total.W = total.W / results.Count;
total.X = total.X / results.Count;
total.Y = total.Y / results.Count;
total.Z = total.Z / results.Count;
return total;
}
代码虽然不复杂，但使用 TPL 可以使其更加优雅。当我们纯粹在 C# 中运行 1000 万次迭代，平均分配给 CPU 核心时，我们得到
这比单线程快了将近 4 倍。我们可以假设，除去任务切换的少量开销，蒙特卡洛例程的扩展性会非常好！
C# 任务并行库
让我们尝试 TPL，特别是 Parallel.For 方法，看看是否能获得良好的性能和更优雅的代码。例如，您可能期望以下代码的效果与使用 Threadpool 一样好，甚至更好，毕竟它更新一些
public void CalculateParallel(ref ResultsdotNet results)
{
double w, x, y, s, z, dw, dx, dy, dz, swx, swy, swz,
varw, varx, vary, varz, ss, vol, sw;
w = x = y = s = z = dw = dx = dy = sw = dz = swx =
swy = swz = varw = varx = vary = varz = 0.0;
ss = 0.2 * (Math.Exp(5.0) - Math.Exp(-5.0));
vol = 3.0 * 7.0 * ss;
long count = 0;
long n = results.Iterations;
Parallel.For( 0, n, f =>
{
Random r = new Random(Environment.TickCount);
x = 1.0 + 3.0 * r.NextDouble();
y = (-3.0) + 7.0 * r.NextDouble();
s = 0.00135 + ss * r.NextDouble();
z = 0.2 * Math.Log(5.0 * s);
double b = Math.Sqrt((x * x) + (y * y)) - 3.0;
double a = Math.Pow(z, 2.0) + (b * b);
if (a < 1.0)
{
sw += 1;
swx += x;
swy += y;
swz += z;
varw = 1.0;
varx += Math.Pow(x, 2.0);
vary += Math.Pow(y, 2.0);
varz += Math.Pow(z, 2.0);
}
Interlocked.Increment(ref count);
});
results.W = vol * sw / n;
results.X = vol * swx / n;
results.Y = vol * swy / n;
results.Z = vol * swz / n;
results.dW = vol * Math.Sqrt((varw / n - Math.Pow((sw / n), 2.0)) / n);
results.dX = vol * Math.Sqrt((varx / n - Math.Pow((swx / n), 2.0)) / n);
results.dY = vol * Math.Sqrt((vary / n - Math.Pow((swy / n), 2.0)) / n);
results.dZ = vol * Math.Sqrt((varz / n - Math.Pow((swz / n), 2.0)) / n);
}
结果
效果并不如 threadpool 那么惊人，因为它表现更差。看起来我们没有在执行工作的线程之间正确地分割工作。
为了解决劳动分工问题，让我们使用 Parallel.ForEach ，并利用分区来更有效地批处理工作
long n = results.Iterations;
OrderablePartitioner<Tuple<long, long>> chunkPart =
Partitioner.Create(0, n, n/100);
long count = 0;
Parallel.ForEach(chunkPart, chunkRange =>
{});
请注意第 3 个参数，我将总迭代次数分区了 100 次。当然，您可能会认为，根据核心数量对分区进行大小调整会与上面的 Threadpool 示例类似。不幸的是，似乎并非如此。随意调整代码并证明这一点。
接下来，让我们使用 interlocking 来确保对状态变量的操作是原子的。由于 .NET 中没有 Interlocked 等价物来添加 double，我们必须自己实现
public static double Add(ref double location1, double value)
{
double newCurrentValue = 0;
while (true)
{
double currentValue = newCurrentValue;
double newValue = currentValue + value;
newCurrentValue = Interlocked.CompareExchange(ref location1, newValue, currentValue);
if (newCurrentValue == currentValue)
return newValue;
}
}
在添加了分区和 interlocking 后，我们得到了什么？
这次，结果看起来是正确的，并且我们在速度上取得了显著的提升，但不足以超越 Threadpool 或我们的混合 C#/C++ 解决方案。
显然，我们如何划分要完成的工作，然后如何将结果重新组合是一个重要的因素，所以让我们比较传统的 Threadpool 和基于任务的方法，使用类似的逻辑来划分工作
public void Run(ResultsdotNet results, int partitions)
{
long iterationsPerTask = results.Iterations / partitions;
Task<ResultsdotNetValue>[] tasks = new Task<ResultsdotNetValue>[partitions];
for (int n = 1; n <= partitions; n++)
tasks[n - 1] = Task.Factory.StartNew < ResultsdotNetValue >(
() => Calculate(iterationsPerTask, new Random(Environment.TickCount))
);
Task.WaitAll(tasks);
GetResults(tasks, ref results, partitions);
}
请注意，（由于缺乏更好的技术）我们再次使用处理器核心数量来划分工作，此外，我们还可以为每个批次使用值类型。
这要好得多，结果可与 Threadpool 示例相媲美，甚至略好。然而，当我再次在我的笔记本电脑（另一个 4 核机器）上运行时，结果却颠倒了。Threadpool 运行速度比 Tasks 稍快，所以我预计您的结果也会有所不同。
混合多线程和互操作
基于我们的单线程测试，我们可能会认为混合使用多线程和 C#/PInvoke 可以为我们带来最快的性能。所以让我们对此进行测试。
不幸的是，Numerical Receipes 中的 ran3() 的原始 C++ 实现不是线程安全的。我们必须对其进行一些重构
#define MBIG 1000000000
#define MSEED 161803398
#define MZ 0
#define FAC (1.0/MBIG)
RandomNumbers::RandomNumbers()
{
iff=0;
}
double RandomNumbers::ran3(long *idum)
{
long mj,mk;
int i,ii,k;
if (*idum < 0 || iff == 0) {
iff=1;
mj=labs(MSEED-labs(*idum));
mj %= MBIG;
ma[55]=mj;
mk=1;
for (i=1;i<=54;i++) {
ii=(21*i) % 55;
ma[ii]=mk;
mk=mj-mk;
if (mk < MZ) mk += MBIG;
mj=ma[ii];
}
for (k=1;k<=4;k++)
for (i=1;i<=55;i++) {
ma[i] -= ma[1+(i+30) % 55];
if (ma[i] < MZ) ma[i] += MBIG;
}
inext=0;
inextp=31;
*idum=1;
}
if (++inext == 56) inext=1;
if (++inextp == 56) inextp=1;
mj=ma[inext]-ma[inextp];
if (mj < MZ) mj += MBIG;
ma[inext]=mj;
return mj*FAC;
}
#undef MBIG
#undef MSEED
#undef MZ
#undef FAC
当我们通过 Threadpool 运行时，我们得到
同样通过 Tasks
令人惊讶的是，我们得到的结果比预期的要差。为什么？最明显的原因是每次线程都创建一个新的 C++ 计算器实例。我们可能可以通过创建计算对象池来做得更好，但这将留待下次更新。
结论
显然，要为“简单”的迭代例程（如蒙特卡洛积分）使用多线程技术并超越单线程方法并不容易。正如我们在 Parallel.For 方法中所见，关键在于有效地划分每个执行单元的工作，然后在最后正确地聚合。
到目前为止，我们发现 Tasks 和 Threadpooling 是我们扩展时的最佳选择，但您选择使用哪一种将很大程度上取决于个人情况。几乎可以肯定的是，通过将此与基于 PInvoke 的计算池以及此处解释的技术结合起来 ，我们将实现最有效的资源利用。
到目前为止，我们只研究了在单台机器上利用 CPU。但是，显卡或者您组织中所有闲置的高性能台式机怎么办？在下一篇文章中，我们将研究构建计算池，并将其与其他技术（如 CUDA 和分布式网格处理）进行比较和组合，以了解如何进一步扩展。
所有改进建议均被愉快接受。
历史
- 创建
- 添加了 Andy Harman 和 Alexey 提出的一些性能改进。
- 标记为 C++ 以便该受众能够看到，以防有一些有用的建议
- 修复了 zip 文件，更新了单线程和
Parallel.For的图像和评论 - 添加了
Task示例，更新了代码 - 添加了
Task和Threadpool与 PInvoke 示例，更新了代码 - 编辑了文章