On Tue, Jun 29, 2010 at 10:23 AM, Marco Caliari <address@hidden> wrote:
Dear users,
I'm wondering if it is possible to vectorize the following code
Hx = rand(n); % n natural number
Hy = rand(n);
Hz = rand(n);
u = rand(n,n,n);
for k = 1:n
us(:,:,k) = zeros(n);
for j = 1:n
us(:,:,k) = us(:,:,k)+Hz(k,j)*(Hy*u(:,:,j)*Hx');
end
end
Thanks,
Marco
_______________________________________________
Help-octave mailing list
address@hidden
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave
This is basically a linear transform of u along each dimension. It can
be done with permute(), reshape(), and three matrix multiplications:
n = 100;
Hx = rand(n); % n natural number
Hy = rand(n);
Hz = rand(n);
u = rand(n,n,n);
tic;
for k = 1:n
us(:,:,k) = zeros(n);
for j = 1:n
us(:,:,k) = us(:,:,k)+Hz(k,j)*(Hy*u(:,:,j)*Hx');
end
end
toc
tic;
v = reshape (Hy*u(:,:), [n,n,n]);
v = reshape (permute (v, [1,3,2]), n^2, n) * Hx.';
v = permute (reshape (v, [n,n,n]), [1,3,2]);
v = reshape (reshape (v, n^2, n) * Hz.', [n,n,n]);
norm ((v-us)(:))
toc
run on my platform:
Elapsed time is 5.40943 seconds.
ans = 0
Elapsed time is 0.1 seconds.
I think I could make a general function for this, taking a
n-dimensional array and n matrices, doing the permute-transpose
gymnastics behind the veil. linear-algebra seems a suitable package.
function Y = nd_linear_transform (X, M1, M2, M3, ...)
returns Y such that
Y(i1,i2,i3,...) = sum (M1(i1,j1) * M2(i2,j2) * M3(i3,j3) * ... * X(j1,
j2, j3, ...)) over all j1, j2, j3...
Would you be interested?