2012-11-15 67 views
9

我無法想象如何向量化這組循環。任何指導將不勝感激。如何向量化嵌套循環

ind_1 = [1,2,3]; 
ind_2 = [1,2,4]; 
K = zeros(3,3,3,3,3,3,3,3,3); 
pp = rand(4,4,4); 

for s = 1:3 
for t = 1:3 
    for k = 1:3 
    for l = 1:3 
    for m = 1:3 
    for n = 1:3 
     for o = 1:3 
     for p = 1:3 
     for r = 1:3 
     % the following loops are singular valued except when 
     % y=3 for ind_x(y) in this case 
     for a_s = ind_1(s):ind_2(s) 
      for a_t = ind_1(t):ind_2(t) 
      for a_k = ind_1(k):ind_2(k) 
      for a_l = ind_1(l):ind_2(l) 
      for a_m = ind_1(m):ind_2(m) 
       for a_n = ind_1(n):ind_2(n) 
       for a_o = ind_1(o):ind_2(o) 
       for a_p = ind_1(p):ind_2(p) 
       for a_r = ind_1(r):ind_2(r) 
        K(s,t,k,l,m,n,o,p,r) = K(s,t,k,l,m,n,o,p,r) + ... 
        pp(a_s, a_t, a_r) * pp(a_k, a_l, a_r) * ... 
        pp(a_n, a_m, a_s) * pp(a_o, a_n, a_t) * ... 
        pp(a_p, a_o, a_k) * pp(a_m, a_p, a_l); 
       end 
       end 
       end 
       end 
      end 
      end 
      end 
      end 
     end 
     end 
     end 
     end 
    end 
    end 
    end 
    end 
end 
end 

編輯:

代碼由每個索引的pp■一個或兩次的產品的值求和產生秩9張量索引爲從1到3,根據值的ind_1ind_2

編輯:

這裏是一個三維例子,但請記住,一個事實,即pp指數只是排列在9D版本不會被保留:

ind_1 = [1,2,3]; 
ind_2 = [1,2,4]; 
K = zeros(3,3,3); 
pp = rand(4,4,4); 

for s = 1:3 
for t = 1:3 
    for k = 1:3 
    % the following loops are singular valued except when 
    % y=3 for ind_x(y) in this case 
    for a_s = ind_1(s):ind_2(s) 
    for a_t = ind_1(t):ind_2(t) 
    for a_k = ind_1(k):ind_2(k) 
     K(s,t,k) = K(s,t,k) + ... 
     pp(a_s, a_t, a_r) * pp(a_t, a_s, a_k) * ... 
     pp(a_k, a_t, a_s) * pp(a_k, a_s, a_t); 
    end 
    end 
    end 
    end 
end 
end 
+0

你能詳細介紹一下這段代碼在做什麼並添加註釋嗎? – igon

+0

@igon:我編輯過。 – erbridge

+0

你可以創建一個2D或3D例子來說明嗎?我們合作會更容易,創建這個例子可以幫助你自己找出一個策略。 – tmpearce

回答

5

窩!非常簡單的解決方案,但不容易找到。順便說一句,我想知道你的公式來自哪裏。

如果你不介意暫時失去一點記憶(以前是2^4^9陣列vs 3^9),那麼你可以推遲第3和第4超平面的積累。

測試用在Unix倍頻3.2.4,它23S(67MB)0.17s(98MB)下降。

function K = tensor9_opt(pp) 

    ppp = repmat(pp, [1 1 1 4 4 4 4 4 4]) ; 
    % The 3 first numbers are variable indices (eg 1 for a_s to 9 for a_r) 
    % Other numbers must complete 1:9 indices in any order 
    T = ipermute(ppp, [1 2 9 3 4 5 6 7 8]) .* ... 
     ipermute(ppp, [3 4 9 1 2 5 6 7 8]) .* ... 
     ipermute(ppp, [6 5 1 2 3 4 7 8 9]) .* ... 
     ipermute(ppp, [7 6 2 1 3 4 5 8 9]) .* ... 
     ipermute(ppp, [8 7 3 1 2 4 5 6 9]) .* ... 
     ipermute(ppp, [5 8 4 1 2 3 6 7 9]) ; 

    % I have not found how to manipulate 'multi-ranges' programmatically. 
    T1 = T (:,:,:,:,:,:,:,:,1:end-1) ; T1(:,:,:,:,:,:,:,:,end) += T (:,:,:,:,:,:,:,:,end) ; 
    T = T1(:,:,:,:,:,:,:,1:end-1,:) ; T (:,:,:,:,:,:,:,end,:) += T1(:,:,:,:,:,:,:,end,:) ; 
    T1 = T (:,:,:,:,:,:,1:end-1,:,:) ; T1(:,:,:,:,:,:,end,:,:) += T (:,:,:,:,:,:,end,:,:) ; 
    T = T1(:,:,:,:,:,1:end-1,:,:,:) ; T (:,:,:,:,:,end,:,:,:) += T1(:,:,:,:,:,end,:,:,:) ; 
    T1 = T (:,:,:,:,1:end-1,:,:,:,:) ; T1(:,:,:,:,end,:,:,:,:) += T (:,:,:,:,end,:,:,:,:) ; 
    T = T1(:,:,:,1:end-1,:,:,:,:,:) ; T (:,:,:,end,:,:,:,:,:) += T1(:,:,:,end,:,:,:,:,:) ; 
    T1 = T (:,:,1:end-1,:,:,:,:,:,:) ; T1(:,:,end,:,:,:,:,:,:) += T (:,:,end,:,:,:,:,:,:) ; 
    T = T1(:,1:end-1,:,:,:,:,:,:,:) ; T (:,end,:,:,:,:,:,:,:) += T1(:,end,:,:,:,:,:,:,:) ; 
    K = T (1:end-1,:,:,:,:,:,:,:,:) ; K (end,:,:,:,:,:,:,:,:) += T (end,:,:,:,:,:,:,:,:) ; 
endfunction 

pp = rand(4,4,4); 
K = tensor9_opt(pp) ; 
+0

非常好。謝謝。 – erbridge