并行计算基于计算数字pi的例子
在当今世界，并行化已经不仅仅是一种优化，而是一种必要性。 随着数据量和计算复杂度的增加，顺序算法的效率往往不够高。 专为高性能计算而设计的Julia语言为并行工作提供了强大的工具，可以显着加快任务的执行速度。 本文通过使用Monte Carlo方法计算数pi的实际示例，演示了Julia中并行计算的各种方法。 我们会想办法的:
*用于性能比较的基本顺序实现
*使用Julia内置工具的几种并行化方法
*并行计算组织不当时出现的问题
*比较不同方法的有效性
*选择并行化策略的实用建议
本文介绍的代码演示了如何正确组织并行计算，避免常见错误，并演示了在Julia中正确使用多线程可以实现的真正加速。
蒙特卡罗计算pi的方法被选为一个清晰的例子，它易于理解，并且可以很好地与并行化进行扩展，使您能够清楚地展示组织并行计算的不同方法的优势。
代码中使用的库:
using Base.Threads-用于多线程计算（创建和管理线程）using Random-用于生成随机数（蒙特卡洛方法的基础）using Statistics-用于统计操作（求和，结果比较）
using Base.Threads
using Random
using Statistics
计算pi的蒙特卡洛方法：朴素顺序版本
此函数返回数字π的近似值，其精度随着随机点数量的增加而提高。 当使用1亿个点时，我们可以预期结果约为3.14159，典型误差约为0.0001。 生成的点越多，结果越接近π的真值。
该算法基于几何概率，并使用简单但强大的面积比思想。 我们正在考虑一个单位正方形，两个轴上的坐标从0到1。 以原点为中心的单位圆的四分之一适合这个正方形。
该函数根据以下原理工作：对于n个点中的每一个，在从0到1的范围内生成两个随机的x和y坐标。 然后通过检查坐标平方和不超过1的条件来检查点是否落在四分之一圆内。 如果点落在圆内，则计数器递增。
Π的最终值被计算为圆中的点的数量与产生的点的总数的四倍比率。 这是从面积的数学比率得出的：正方形的面积等于1，圆的四分之一的面积等于π/4，因此圆中的点数与总点数的比率接近π/4。
该方法的基础是几何形状的区域之间的比率。 单位正方形的面积为1，单位圆的四分之一的面积为π/4。 当点在正方形上均匀随机分布时，点落入圆的概率等于面积之比，即π/4。 因此，这个概率乘以4的统计估计给出π的近似值。
这种方法非常容易实现和理解。 它不需要复杂的数学计算或特殊知识。 然而，该方法的准确性直接取决于随机数生成器的质量和点数。 该算法具有线性时间复杂度和恒定内存复杂度，这使得它对于演示目的是有效的。
该方法的收敛率大约为1/√n，这意味着要将误差减少10倍，就必须将点的数量增加100倍。 尽管收敛速度相对较慢，但该方法非常适合演示统计建模和并行计算的原理，因为它易于并行化和修改。
function compute_pi_serial_naive(n)
count = 0
for i in 1:n
x = rand()
y = rand()
if x^2 + y^2 <= 1.0
count += 1
end
end
return 4.0 * count / n
end
向量实现
它的工作原理与以前的方法一样，但在实现方面存在关键差异。：此函数使用相同的蒙特卡洛方法计算数字pi，但使用矢量化方法而不是逐个元素循环。
该函数不是在循环中一次处理一个点，而是首先生成两个具有所有随机x和y坐标的大型数组。 然后，它使用元素向量运算同时计算所有点的距离平方。 之后，在一次运算中计算满足落入圆的条件的点数。
这种方法更好地利用了现代处理器的并行数据处理功能，并针对处理大型阵列进行了优化，但需要更多内存来存储中间阵列。
function compute_pi_serial_vectorized(n)
x_arr = rand(n)
y_arr = rand(n)
distances_sq = x_arr.^2 .+ y_arr.^2
count = sum(distances_sq .<= 1.0)
return 4.0 * count / n
end
它的工作原理与蒙特卡洛原理相同，但与并行计算：此版本使用多线程通过在多个处理器内核之间分配工作来加快计算速度。
该算法基于可用流的数量将点的总数分成近似相等的部分（块）。 每个线程都使用自己的数据部分，使用自己的带有唯一种子的随机数生成器来确保序列独立性。 线程在它们的部分中独立地计数圆中的点的数量，然后将结果汇总到一个公共计数器中。 因为每个线程都与局部变量一起工作，并且只有最终的减少需要组合结果。 对每个流使用单独的随机数生成器保证了并行执行时统计建模的正确性。
function compute_pi_parallel_reduction(n)
total_count = 0
@threads for chunk in 1:Threads.nthreads()
local_count = 0
chunk_size = n ÷ Threads.nthreads()
start_idx = (chunk - 1) * chunk_size + 1
end_idx = chunk == Threads.nthreads() ? n : chunk * chunk_size
rng = MersenneTwister(chunk + 42)
for i in start_idx:end_idx
x = rand(rng)
y = rand(rng)
if x^2 + y^2 <= 1.0
local_count += 1
end
end
total_count += local_count
end
return 4.0 * total_count / n
end
它使用相同的蒙特卡洛方法工作，但使用批量分离和结果的显式存储。：此版本也使用多线程，但在线程之间聚合结果的方式有所不同。
每个线程将其部分计数器写入预先创建的结果数组，而不是直接将结果求和到公共变量中。 数组的每个元素对应于它自己的流，这消除了在记录期间发生冲突的可能性。 所有线程完成后，数组元素的最终求和发生。 这种方法确保了线程之间数据的清晰分离，并完全消除了计算过程中同步的需要，因为每个线程只使用自己的数组元素。 该方法在线程数量预先已知并且在计算期间保持恒定时特别有效。
function compute_pi_parallel_batch(n)
nthreads = Threads.nthreads()
chunk_size = n ÷ nthreads
results = zeros(Int64, nthreads)
@threads for thread in 1:nthreads
local_count = 0
start_idx = (thread - 1) * chunk_size + 1
end_idx = thread == nthreads ? n : thread * chunk_size
rng = MersenneTwister(thread + 100)
for i in start_idx:end_idx
x = rand(rng)
y = rand(rng)
local_count += (x^2 + y^2 <= 1.0)
end
results[thread] = local_count
end
return 4.0 * sum(results) / n
end
它使用相同的蒙特卡洛方法，但使用异步任务而不是直接流控制。：该版本通过tasks机制实现并行计算，为并行化提供了更灵活的方法。
而不是使用指令跨线程显式分发工作 @threads，该函数通过宏创建独立的异步任务 @spawn. 每个任务都计算自己的点部分并返回一个本地计数器。 所有任务都是竞争性启动的，其结果是使用操作收集的 fetch，等待任务完成并检索结果。 这种方法允许Julia调度程序更灵活地跨可用线程分配任务，并更好地利用系统资源，特别是当任务可能执行不均匀时。 该方法还可以更轻松地扩展和组织更复杂的并行工作流。
function compute_pi_parallel_spawn(n)
nchunks = Threads.nthreads()
chunk_size = n ÷ nchunks
tasks = []
for chunk in 1:nchunks
start_idx = (chunk - 1) * chunk_size + 1
end_idx = chunk == nchunks ? n : chunk * chunk_size
task = Threads.@spawn begin
local_count = 0
rng = MersenneTwister(chunk + 200)
for i in start_idx:end_idx
x = rand(rng)
y = rand(rng)
local_count += (x^2 + y^2 <= 1.0)
end
local_count
end
push!(tasks, task)
end
total_count = sum(fetch.(tasks))
return 4.0 * total_count / n
end
接下来，实现测试和比较所有版本算法的功能。：此脚本对计算pi的蒙特卡洛方法的所有实现进行全面比较，评估每种方法的性能和准确性。
该脚本依次运行具有相同点数的所有五个版本的算法，测量每个函数的执行时间。 对于每个实现，存储π的近似值和计算时间。 在所有测试完成后，进行比较性能分析，其中相对于天真的顺序版本计算加速度，将其作为基线。
该函数为每个版本提供两个关键指标：相对于基法的加速度系数和相对于π精确值的绝对计算误差。 这使我们能够同时评估并行化的有效性和每种方法的统计可靠性，确定不同用例的最佳方法。
function benchmark_all_versions_fixed(n = 100_000_000)
println("线程数: ", nthreads())
println("迭代次数: ", n)
println()
versions = [
("天真和一致", compute_pi_serial_naive),
("矢量化顺序", compute_pi_serial_vectorized),
("平行(减少)", compute_pi_parallel_reduction),
("并行（批量）", compute_pi_parallel_batch),
("平行（产卵）", compute_pi_parallel_spawn)
]
results = []
times = []
for (name, func) in versions
print("测试$名称。.. ")
time_taken = @elapsed begin
pi_approx = func(n)
end
push!(results, pi_approx)
push!(times, time_taken)
println("π≈ $pi_approx,时间:$(round(time_taken,digits=3))秒")
end
println("\n" * "="^50)
println("性能比较:")
base_time = times[1]
for (i, (name, _)) in enumerate(versions)
speedup = base_time / times[i]
println("名称:加速度(round(speedup,digits=2))x")
end
println("\n" * "="^50)
println("检查准确性:")
actual_pi = π
for (i, (name, _)) in enumerate(versions)
error = abs(results[i] - actual_pi)
println("名称：错误（圆（错误，数字=6））")
end
end
benchmark_all_versions_fixed(100_000_000)
输出
性能优胜者：并行版本与 @spawn 它显示出1.82x的最佳加速度，接近7个线程上的理想加速度。
方法的准确性：所有实现都提供了可接受的精度，误差为10⁻⁴，对应于蒙特卡洛方法1亿点的理论误差。
并行化效率：与顺序实现相比，所有三个并行版本都显示出显着的加速，展示了Julia中多线程的实际好处。
为什么矢量化版本最慢？
内存拷贝过多的问题:
-正在创建临时数组 x_arr.^2 和 y_arr.^2
-正在形成中间阵列 distances_sq
-为条件创建布尔数组 .<= 1.0
低效的缓存使用:
-大型阵列不适合处理器缓存
-不断访问RAM，而不是快速寄存器
-缺乏数据本地化
内存管理的开销:
-临时数组的分配和解除分配
-冗余数据移动操作
矢量化操作隐藏条件转换，这些条件转换由现代处理器在逐元素处理中有效预测。 矢量化对于复杂的数学运算是有效的，但是对于具有条件的简单算术运算，内存管理的开销超过了并行数据处理的好处。