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

4 次查看（过去 30 天）
AP
AP 2013-8-7
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
James Tursa 2013-8-7
How large is large? You could form the tensor as a 3x3xN matrix and then use this FEX submission by Bruno Luong:
http://www.mathworks.com/matlabcentral/fileexchange/27680-multiple-eigen-values-for-2x2-and-3x3-matrices

更多回答（1 个）

Richard Brown
Richard Brown 2013-8-7
编辑：Richard Brown 2013-8-7
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.

类别

在 Help Center 和 File Exchange 中查找有关 Logical 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by