[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: loops over a 3d array
From: |
Jaroslav Hajek |
Subject: |
Re: loops over a 3d array |
Date: |
Tue, 29 Jun 2010 13:13:01 +0200 |
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?
--
RNDr. Jaroslav Hajek, PhD
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz