1

X、Y 和 z 是表示表面的坐标。为了计算一些量，我们称之为流，在表面的 i,j 点，我需要计算所有其他点 (i0,j0) 的贡献。为此，我需要例如知道点 i0、j0 和所有其他点 (alpha) 之间的角度 cos。然后来自 i0,j0 的所有贡献必须乘以一些常数并相加。zv0 在每个点 i,j 是最终需要的结果。

我想出了下面写的一些代码，它似乎非常不合适。首先，它减慢了程序的其余部分，并且似乎使用了所有可用内存。我的系统有 4gb 物理内存和 12gb 交换文件，它总是内存不足，尽管所有变量大小都不大于 10kb。请帮助解决加速/矢量化和内存问题。

parfor i0=2:1:length(x00);
  for j0=2:1:length(y00);
    zv=red3dfunc(X0,Y0,f,z0,i0,j0,st,ang,nx,ny,nz);
    zv0=zv0+zv;
  end
end


function[X,Y,z,zv]=red3dfunc(X,Y,f,z,i0,j0,st,ang,Nx,Ny,Nz)
x1=X(i0,j0);
y1=Y(i0,j0);
z1=z(i0,j0);
alpha=zeros(size(X));
betha=zeros(size(X));
r=zeros(size(X));
XXa=X-x1;
YYa=Y-y1;
ZZa=z-z1;
VEC=((XXa).^2+(YYa).^2+(ZZa).^2).^(1/2);
VEC(i0,j0)=VEC(i0-1,j0-1);
XXa=XXa./VEC;
YYa=YYa./VEC;
ZZa=ZZa./VEC;
alpha=-(Nx(i0,j0).*XXa+Ny(i0,j0).*YYa+Nz(i0,j0).*ZZa);
betha=Nx.*XXa+Ny.*YYa+Nz.*ZZb;
r=VEC;
zv=(1/pi)*st^2*ang.*f.*(alpha).*betha./r.^2;
4

1 回答 1

3

要做到这一点，显而易见的事情是使用 Kroneker 产品。对于维度为 nAxmA 和 nBxmB 的矩阵，matlab 函数是 kron(A,B)。这个函数将返回维度矩阵 (nA*nB)x(mA*mB)，看起来像

[a11*B a12*B ... a1mA*B;
.......................;
anA1*B ........ anAmA*B]

所以你的问题可以通过引入矩阵 I=ones(size(X)) 来解决。然后，您将定义您的 XXa、YYa、ZZa 和 VEC 矩阵，没有任何循环为

XXa = kron(I,X)-kron(X,I);
YYa = kron(I,Y)-kron(Y,I);
ZZa = kron(I,Z)-kron(Z,I);
VEC=((XXa).^2+(YYa).^2+(ZZa).^2).^(1/2);

然后，您将找到任何 i0,j0 的 VEC（如果您将 n 和 m 定义为 X 的大小分量）

VEC((1+n*(i0-1)):(n*i0),(1+m*(j0-1)):(m*j0))
于 2012-12-06T17:47:15.970 回答