如何对Scilab中的所有矩阵元素进行运算？

Question

我正在尝试模拟无限板上随时间的热分布。为此，我编写了一个 Scilab 脚本。现在，它的关键点是计算所有板点的温度，并且必须为我想要观察的每个时间实例完成：

for j=2:S-1
    for i=2:S-1
        heat(i, j) = tcoeff*10000*(plate(i-1,j) + plate(i+1,j) - 4*plate(i,j) + plate(i, j-1) + plate(i, j+1)) + plate(i,j);          
    end;
end

问题是，如果我想为 100x100 点的板做这件事，这意味着，在这里（它只针对内部，没有边界条件），我将不得不循环 98x98 = 9604次，每次都计算给定 i,j 点的热量。如果我想观察 100 秒，每步 1 秒，我必须重复 100 次，总共进行 960,400 次迭代。这需要相当长的时间，我想避免它。高达 50x50 的印版，这一切都发生在合理的 4-5 秒时间范围内。

现在我的问题是 - 是否有必要使用 for 循环来完成所有这些工作？ Scilab 中是否有任何内置的聚合函数，可以让我对矩阵的所有元素执行此操作？我还没有找到方法的原因是每个点的结果都取决于其他矩阵点的值，这让我用嵌套循环来做。关于如何使其更快的任何想法表示赞赏。

Answer 1

在我看来，您想计算热场和特定扩散模式的二维相互关系。这种模式可以被认为是一个 "filter" 核，这是一种用线性滤波器矩阵修改图像的常用方法。您的 "filter" 是：

F=[0,1,0;1,-4,1;0,1,0];

如果您安装图像处理工具箱 (IPD)，您将有一个 MaskFilter 函数来执行此 2D 互相关。

S=500;
plate=rand(S,S);
tcoeff=1;

//your solution with nested for loops
t0=getdate();
for j=2:S-1
  for i=2:S-1
    heat(i, j) = tcoeff*10000*(plate(i-1,j)+plate(i+1,j)-..
    4*plate(i,j)+plate(i,j-1)+plate(i, j+1))+plate(i,j);          
  end
end
t1=getdate();
T0=etime(t1,t0);
mprintf("\nNested for loops: %f s (100 %%)",T0);

//optimised nested for loop
F=[0,1,0;1,-4,1;0,1,0];   //"filter" matrix
F=tcoeff*10000*F;
heat2=zeros(plate);
t0=getdate();
for j=2:S-1
  for i=2:S-1
    heat2(i,j)=sum(F.*plate(i-1:i+1,j-1:j+1));
  end
end
heat2=heat2+plate;
t1=getdate();
T2=etime(t1,t0);
mprintf("\nNested for loops optimised: %f s (%.2f %%)",T2,T2/T0*100);

//MaskFilter from IPD toolbox
t0=getdate();
heat3=MaskFilter(plate,F);
heat3=heat3+plate;
t1=getdate();
T3=etime(t1,t0);
mprintf("\nWith MaskFilter: %f s (%.2f %%)",T3,T3/T0*100);

disp(heat3(1:10,1:10)-heat(1:10,1:10),"Difference of the results (heat3-heat):");

请注意，MaskFilter 在应用过滤器之前填充图像（原始矩阵），据我所知，它在边界上使用了 "mirror" 数组。您应该检查此行为是否适合您。

速度提升了大约*320（执行时间是你原来代码的0.32%）。够快吗？

理论上它可以用两个 2D 傅里叶变换完成（也许内置 Scilab mfft）但它可能不会比这更快。看这里：http://mailinglists.scilab.org/Image-processing-filter-td2618144.html#a2618168

Answer 2

请注意，向量化运算与并行计算之间存在很大差异，正如我所解释的 here。尽管矢量化可能会稍微提高性能，但这无法与通过 GPU 计算（例如 OpenCL）等实现的性能相提并论。我将尝试解释您代码的矢量化形式，而不会过多介绍细节。考虑这些：

S = ...;
tcoeff = ...;

function Plate = plate(i, j)
    ...;
endfunction

function Heat = heat(i, j)
    ...;
endfunction

现在你可以定义一个 meshgrid:

x = 2 : S - 1;
y = 2 : S - 1;
[M, N] = meshgrid(x,y);
Result = feval(M, N, heat);

feval 是这里的关键，它将在 M 和 N 矩阵上广播 feval 函数。

Answer 3

您的方案是拉普拉斯算子应用于矩形网格的有限差分方案。如果您选择自由度的按行或按列编号（此处为板（i，j））以便将它们视为向量，则应用 "discrete" 拉普拉斯算子可以通过乘以左侧的稀疏矩阵（非常快）这在以下文档中有特别好的解释：

https://www.math.uci.edu/~chenlong/226/FDMcode.pdf。

实现在 Matlab 中进行了描述，但在 Scilab 中很容易翻译。

如何对Scilab中的所有矩阵元素进行运算？

How to perform operation for all matrix elements in Scilab?

performance

loops

matrix

scilab

heat