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

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

代码是数值有限差分法中子程序的一部分，用于求解网格上的 6 个方程组，因此这个问题可能与任何从事矩阵计算的人相关 方程组，尤其是PDEs。非常感谢您提出优化代码的建议！

Answer 1

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

%# 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 时，请考虑是否应使用稀疏矩阵。