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

6 ビュー (過去 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:
  2 件のコメント
AP
AP 2013 年 8 月 7 日
Thanks @James. It is of size 195x156x150x30.
Richard Brown
Richard Brown 2013 年 8 月 7 日
That is a neat trick!

サインインしてコメントする。

その他の回答 (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 ExchangeLinear Algebra についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by