2

有没有一种聪明的方法来矢量化一个将元素分配给矩阵子矩阵的 for 循环？
最初，我有两个 for 循环： 

U=zeros(6*(M-2),M-2);
for k=2:M-3  
    i=(k-1)*6+1; 
    for j=2:M-3
        U(i:i+5,j)=A*temp(i:i+5,j)+B*temp(i:i+5,j-1)+C*temp(i:i+5,j+1)+D*temp(i-6:i-1,j)+E*temp(i+6:i+11,j);
    end
end

然后我对内部循环进行矢量化，这样代码现在读取

U=zeros(6*(M-2),M-2);
j=2:M-2;
for k=2:M-3
    i=(k-1)*6+1;
    U(i:i+5,j)=A*temp(i:i+5,j)+B*temp(i:i+5,j-1)+C*temp(i:i+5,j+1)+D*temp(i-6:i-1,j)+E*temp(i+6:i+11,j);
end

这将我的 CPU 时间减少了 90% 以上，所以我想知道是否可以对外部循环执行相同的操作，但这似乎有点棘手，因为我在 U 矩阵中分配了 (6x1) 矩阵。我试过

U=zeros(6*(M-2),M-2);
k=2:M-3;
i=(k-1)*6+1;
j=2:M-2;
U(i:i+5,j)=A*temp(i:i+5,j)+B*temp(i:i+5,j-1)+C*temp(i:i+5,j+1)+D*temp(i-6:i-1,j)+E*temp(i+6:i+11,j);

但这失败了，因为 i:i+5 只取出我想要的前 6 个索引。

我也尝试使用 reshape() 函数将矩阵转换为向量，但一次分配给多个元素块似乎仍然很困难。代码中总共有三个这样的 for 循环，所以我想另一种优化是以某种方式并行化它们。但是，如果无法访问并行工具箱，在我看来，如果可能的话，矢量化是一个很好的解决方案。

该代码是数值有限差分法中子程序的一部分，用于求解网格上的 6 个方程组，因此这个问题可能与任何在方程组（尤其是PDE ）上进行矩阵计算的人相关。对优化代码的建议将不胜感激！

4

2 回答 2

0

为了选择矩阵的非矩形部分，您需要使用线性索引：在 3x3 矩阵中 A、A(3,3)==A(9)和A([1 3 5 7 9])是无法通过行/列索引方法实现的向量。

该sub2ind函数将行/列索引转换为线性索引，因此您可以在表单中使用它sub2ind(size(U),i:i+5,j)来获取一个 U 块的线性索引。将循环更改为仅执行收集线性索引的工作，然后您可以在外面说循环：

U(ind_U) = A*temp(ind_A) + B*temp(ind_B) ...

此外，在处理 FDM 或 FEM 时，请考虑是否应该使用稀疏矩阵。

于 2011-03-27T11:51:42.303 回答
0

要了解如何在没有循环的情况下在一行中编写赋值，将数组绘制temp为矩形可能会有所帮助。然后，将组合成的不同加法只不过是（或子网格，如果您想跟踪单个元素，这将导致特定元素）的U子矩形，它们向左，向右，分别为顶部、底部。temptempU

%# define row, column shifts
rowShift = 6;
colShift = 1;

%# That's how we'd like to shift 
%# U(i:i+5,j)=A*temp(i:i+5,j)+B*temp(i:i+5,j-1)+C*temp(i:i+5,j+1)+
%# D*temp(i-6:i-1,j)+E*temp(i+6:i+11,j);

%# assign U
U = A * temp(rowShift+1 : end-rowShift, colShift+1 : end-colShift) +... 
    B * temp(rowShift+1 : end-rowShift, 1 : end-2*colShift) + ...
    C * temp(rowShift+1 : end-rowShift, 2*colShift+1 : end) + ...
    D * temp(1 : end-2*rowShift, colShift+1 : end-colShift) + ...
    E * temp(2*rowShift+1 : end, colShift+1 : end-colShift);
于 2011-03-27T18:26:39.307 回答