2

要素を行列の部分行列に割り当てる for ループをベクトル化する賢い方法はありますか?
最初は、2 つの for ループがありました。 

U=zeros(6*(M-2),M-2);
for k=2:M-3  
    i=(k-1)*6+1; 
    for j=2:M-3
        U(i:i+5,j)=A*temp(i:i+5,j)+B*temp(i:i+5,j-1)+C*temp(i:i+5,j+1)+D*temp(i-6:i-1,j)+E*temp(i+6:i+11,j);
    end
end

次に、内側のループをベクトル化して、コードが読み取れるようにしました。

U=zeros(6*(M-2),M-2);
j=2:M-2;
for k=2:M-3
    i=(k-1)*6+1;
    U(i:i+5,j)=A*temp(i:i+5,j)+B*temp(i:i+5,j-1)+C*temp(i:i+5,j+1)+D*temp(i-6:i-1,j)+E*temp(i+6:i+11,j);
end

これにより CPU 時間が 90% 以上短縮されたので、外側のループでも同じことができるかどうか疑問に思いましたが、U マトリックス内の (6x1) マトリックスに割り当てるため、少し難しいようです。私は試した

U=zeros(6*(M-2),M-2);
k=2:M-3;
i=(k-1)*6+1;
j=2:M-2;
U(i:i+5,j)=A*temp(i:i+5,j)+B*temp(i:i+5,j-1)+C*temp(i:i+5,j+1)+D*temp(i-6:i-1,j)+E*temp(i+6:i+11,j);

i:i+5 は必要な最初の 6 つのインデックスしか取得しないためです。

また、 reshape() 関数を使用して行列をベクトルに変換しようとしましたが、要素の複数のブロックを一度に割り当てるのはまだ難しいようです。コードにはそのような for ループが全部で 3 つあるので、代替の最適化は何らかの形でそれらを並列化することだと思います。ただし、並列ツールボックスにアクセスできない場合は、可能であればベクトル化が適切なソリューションのように思えます。

このコードは、グリッド上の 6 つの連立方程式を解くための数値有限差分法のサブルーチンの一部であるため、この質問は、連立方程式、特に偏微分方程式で行列計算を行うすべての人に関連する可能性があります。コードを最適化するための提案をいただければ幸いです。

4

2 に答える 2

0

行列の四角形以外の部分を選択するには、線形インデックスを使用する必要がありA(3,3)==A(9)ますA([1 3 5 7 9])。

このsub2ind関数は行/列インデックスを線形インデックスに変換するため、これを形式で使用して、sub2ind(size(U),i:i+5,j)U の 1 つのブロックの線形インデックスを取得できます。線形インデックスを収集する作業のみを行うようにループを変更すると、outside と言うことができます。ループ:

U(ind_U) = A*temp(ind_A) + B*temp(ind_B) ...

また、FDM または FEM を扱っているときはいつでも、疎行列を使用する必要があるかどうかを検討してください。

于 2011-03-27T11:51:42.303 に答える
0

割り当てをループなしで 1 行に記述する方法を理解するには、配列tempを四角形として描画すると役立つ場合があります。次に、結合されるさまざまな被加数は、左U、右、上、下、それぞれ。temptempU

%# define row, column shifts
rowShift = 6;
colShift = 1;

%# That's how we'd like to shift 
%# U(i:i+5,j)=A*temp(i:i+5,j)+B*temp(i:i+5,j-1)+C*temp(i:i+5,j+1)+
%# D*temp(i-6:i-1,j)+E*temp(i+6:i+11,j);

%# assign U
U = A * temp(rowShift+1 : end-rowShift, colShift+1 : end-colShift) +... 
    B * temp(rowShift+1 : end-rowShift, 1 : end-2*colShift) + ...
    C * temp(rowShift+1 : end-rowShift, 2*colShift+1 : end) + ...
    D * temp(1 : end-2*rowShift, colShift+1 : end-colShift) + ...
    E * temp(2*rowShift+1 : end, colShift+1 : end-colShift);
于 2011-03-27T18:26:39.307 に答える