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

6 views (last 30 days)
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?

Accepted Answer

James Tursa
James Tursa 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:

More Answers (1)

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

Categories

Find more on Linear Algebra in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!