Asked by AP
on 7 Aug 2013

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?

Answer by James Tursa
on 7 Aug 2013

Accepted Answer

AP
on 7 Aug 2013

Thanks @James. It is of size 195x156x150x30.

Richard Brown
on 7 Aug 2013

That is a neat trick!

Sign in to comment.

Answer by Richard Brown
on 7 Aug 2013

Edited by 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.

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.