## vectorizing calculation of eigen values of a large multi-dimensional array

on 7 Aug 2013

### James Tursa (view profile)

Dear All,
I have a 3D rectangular domain. Each point is associated with a vector(u,v,w) through time. u, v, and w are of size nx×ny×nz×nt, which represents the resolution of the domain in x, y, z, and time. I want to calculate eigen values of a tensor which is obtained from the vector field at each point and for all the times in a vectorized fashion. It is very easy to use for-loop over time and position but u, v, w are large. The following is just a working example.
nx = 10; ny = 10; nz = 10; nt = 5;
u = ones(nx, ny, nz, nt);
v = ones(nx, ny, nz, nt);
w = ones(nx, ny, nz, nt);
all_eigen_vals = zeros(nx,ny,nz,nt,3);
for t=1:nt
for k=1:nz
for j=1:ny
for i=1:nx
tensor=[u(i,j,k,t)^2, u(i,j,k,t)*v(i,j,k,t), u(i,j,k,t)*w(i,j,k,t);
u(i,j,k,t)*v(i,j,k,t), v(i,j,k,t)^2, v(i,j,k,t)*w(i,j,k,t);
u(i,j,k,t)*w(i,j,k,t), v(i,j,k,t)*w(i,j,k,t), w(i,j,k,t)^2]
all_eigen_vals(i,j,k,t,:) = eig(tensor);
end
end
end
end
Could someone help me?

### James Tursa (view profile)

on 7 Aug 2013

How large is large? You could form the tensor as a 3x3xN matrix and then use this FEX submission by Bruno Luong:

AP

### AP (view profile)

on 7 Aug 2013
Thanks @James. It is of size 195x156x150x30.
Richard Brown

### Richard Brown (view profile)

on 7 Aug 2013
That is a neat trick!

### Richard Brown (view profile)

on 7 Aug 2013
Edited by Richard Brown

### Richard Brown (view profile)

on 7 Aug 2013

You'll get a bit of a performance boost simply by looping over linear indices rather than nesting (this is more than twice as fast on my system):
evals = zeros(3, numel(u));
for i = 1:numel(u)
tensor=[u(i)^2, u(i)*v(i), u(i)*w(i);
u(i)*v(i), v(i)^2, v(i)*w(i);
u(i)*w(i), v(i)*w(i), w(i)^2];
evals(:,i) = eig(tensor);
end
evals = reshape(evals,3,nx,ny,nz,nt);
note that I've put changed the order of indexing in the evals matrix so that each set of 3 eigenvalues ought to be contiguous in memory. You'll probably want to double check the right bits are going in the right places.