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

5 views (last 30 days)

Show older comments

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?

##### 0 Comments

### Accepted Answer

James Tursa
on 7 Aug 2013

### More Answers (1)

Richard Brown
on 7 Aug 2013

Edited: Richard Brown
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.

##### 0 Comments

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!