MATLAB中依赖型嵌套for循环的并行计算优化方案咨询

你当前的嵌套循环计算有限差分，数组依赖前一时间步的全局值和当前空间步的相邻值，parfor确实无法直接处理这种有依赖关系的循环，以下是几个可行的替代优化方案：

原代码 #

for l=2:timesteps         %start time loop
    
    %now calculate new concentration values for the current timestep
    for k=2:nodes-1         %start spatial loop over inner nodes
        conc(k,l)=conc(k,l-1)+(r*(conc(k-1,l-1)-(2*conc(k,l-1))+conc(k+1,l-1)));
    end                   %end spatial loop over inner nodes
    
    %calculate boundary nodes
    conc(1,l)= c_surface;
    conc(nodes,l)=conc(nodes)+(r*((2*conc(nodes-1))-(2*conc(nodes))));
      
end                      %end time loop

注：原代码边界计算中conc(nodes)和conc(nodes-1)应该是上一时间步的值（conc(nodes,l-1)、conc(nodes-1,l-1)），否则会用当前时间步未更新的值，已在后续方案中修正。

1. 向量化优化（优先尝试，无需并行工具箱） #

Matlab对向量/矩阵运算的底层优化远好于显式循环，把内层空间循环改成向量操作，直接利用Matlab的BLAS加速：

for l=2:timesteps
    % 内层空间循环向量化，一次性计算所有内部节点
    conc(2:nodes-1,l) = conc(2:nodes-1,l-1) + r*(conc(1:nodes-2,l-1) - 2*conc(2:nodes-1,l-1) + conc(3:nodes,l-1));
    
    % 修正后的边界计算
    conc(1,l) = c_surface;
    conc(nodes,l) = conc(nodes,l-1) + r*(2*conc(nodes-1,l-1) - 2*conc(nodes,l-1));
end

这种方式不用改动并行逻辑，就能获得数倍甚至数十倍的速度提升，是最划算的优化手段。

2. GPU加速（需GPU硬件支持） #

如果有NVIDIA GPU，Matlab的gpuArray可以把数组放到GPU上计算，同一时间步的所有空间点计算可以并行执行，完美匹配你的计算模式：

% 把数组转移到GPU
conc = gpuArray(conc);

for l=2:timesteps
    conc(2:nodes-1,l) = conc(2:nodes-1,l-1) + r*(conc(1:nodes-2,l-1) - 2*conc(2:nodes-1,l-1) + conc(3:nodes,l-1));
    conc(1,l) = c_surface;
    conc(nodes,l) = conc(nodes,l-1) + r*(2*conc(nodes-1,l-1) - 2*conc(nodes,l-1));
end

% 把结果从GPU转回CPU
conc = gather(conc);

GPU擅长这种大规模并行的数组运算，当节点数nodes很大时，提速效果非常明显。

3. 基于spmd的CPU空间域并行（需并行计算工具箱） #

如果只能用CPU并行，可以用spmd将空间域拆分给多个工作进程，每个进程负责一部分节点的计算，时间步之间通过进程通信传递相邻边界数据：

spmd
    % 每个worker分配专属的空间节点区间
    local_nodes = partition(1:nodes, numlabs, labindex);
    local_conc = conc(local_nodes, :);
    
    for l=2:timesteps
        left_bdy = [];
        right_bdy = [];
        
        % 和相邻worker交换上一时间步的边界数据
        if labindex > 1
            send(local_conc(1,l-1), labindex-1);
            left_bdy = receive(labindex-1);
        end
        if labindex < numlabs
            send(local_conc(end,l-1), labindex+1);
            right_bdy = receive(labindex+1);
        end
        
        % 计算本地内部节点
        for k=2:length(local_nodes)-1
            local_conc(k,l) = local_conc(k,l-1) + r*(local_conc(k-1,l-1) - 2*local_conc(k,l-1) + local_conc(k+1,l-1));
        end
        
        % 处理本地边界（含全局边界和进程间边界）
        if local_nodes(1) == 1
            local_conc(1,l) = c_surface;
        else
            local_conc(1,l) = local_conc(1,l-1) + r*(left_bdy - 2*local_conc(1,l-1) + local_conc(2,l-1));
        end
        
        if local_nodes(end) == nodes
            local_conc(end,l) = local_conc(end,l-1) + r*(2*local_conc(end-1,l-1) - 2*local_conc(end,l-1));
        else
            local_conc(end,l) = local_conc(end,l-1) + r*(local_conc(end-1,l-1) - 2*local_conc(end,l-1) + right_bdy);
        end
    end
end

% 合并所有worker的计算结果
conc = gather(local_conc);

注意：这种方式需要额外的进程通信开销，只有当nodes非常大时，并行收益才会超过通信成本。

内容的提问来源于stack exchange，提问作者Thilo