同时分配给矩阵的多个子矩阵。通过矢量化索引可能进行的优化

Question

有没有一种聪明的方法来向量化 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() 函数将矩阵转换为向量，但一次分配给多个元素块似乎仍然很困难。代码中总共有三个这样的forops，所以我想另一种优化是以某种方式并行化它们。然而，如果可能的话，如果无法访问并行工具箱，矢量化在我看来是一个很好的解决方案。

该代码是用于求解网格上6个方程系统的数值有限差异方法的子例程的一部分，因此，此问题可能与对方程式系统进行矩阵计算的任何人有关，尤其是 偏微分方程. 。优化代码的建议将不胜感激！

Answer 1

要了解如何在一行中编写作业而无需循环，绘制数组可能会有所帮助 temp 作为一个矩形。然后，不同的加数将组合起来 U 只不过是子矩形 temp （或者子网格，如果你想跟踪单个元素 temp 这将导致特定元素 U）分别向左、右、上、下移动。

%# 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);

Answer 2

为了选择矩阵的非矩形部分，您需要使用线性索引：在 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 时，请考虑是否应该使用稀疏矩阵。