Matlab для векторизации циклов и памяти

X, Y и z — координаты поверхности. Чтобы рассчитать некоторую величину, назовем ее потоком, в точке i,j поверхности, мне нужно вычислить вклад от всех других точек (i0,j0). Для этого мне нужно, например, знать угол между точкой i0, j0 и всеми другими точками (альфа). Тогда все вклады от i0,j0 нужно перемножить на некоторые константы и сложить. zv0 в каждой точке i,j является окончательным требуемым результатом.

Я придумал код, написанный ниже, и он кажется крайне неуместным. Во-первых, это замедляет остальную часть программы и, кажется, использует всю доступную память. Моя система имеет 4 ГБ физической памяти и 12 ГБ файла подкачки, и ей всегда не хватает памяти, хотя размеры всех переменных не превышают 10 КБ. Пожалуйста, помогите с проблемами ускорения/векторизации и памяти.

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;

matlab memory vectorization
person user1364012    schedule 06.12.2012    source источник
comment
Насколько велики x00 и y00?   -  person HerrKaputt    schedule 06.12.2012
comment
так что X0=сетка(x00,y00). Если вам нужны цифры, то они не менее 50.   -  person user1364012    schedule 06.12.2012

Ответы (1)


arrow_upward
3
arrow_downward

Очевидно, что для этого нужно использовать продукт Kroneker. Функция Matlab — kron(A,B) для матриц размерностей nAxmA и nBxmB. Эта функция вернет матрицу размерности (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);

Затем вы найдете VEC для любых i0,j0 как (если вы определяете n и m как компоненты размера X)

VEC((1+n*(i0-1)):(n*i0),(1+m*(j0-1)):(m*j0))
person Krivoi    schedule 06.12.2012
comment
солидный ответ. я обязательно попробую - person user1364012; 06.12.2012
comment
я пытался использовать этот подход, жаль, что он работает еще медленнее, чем с циклом for - person user1364012; 07.12.2012
comment
Это как-то грустно. Единственное, что я могу вам порекомендовать, это использовать repmat(X,size(X)) вместо kron(I,X). Возможно, это немного улучшит скорость. - person Krivoi; 08.12.2012