【MATLAB】コードを改善(ベクトル化)して高速化

前回,MATLABを用いて差分法による電位分布計算をしたが,格子数が多いと計算時間が1分以上かかってしまった

今回はforやifを無くしてベクトル化するなどソースコードを改善し,高速化を試みた

コードが遅い具体的な原因

行列の要素ごとに計算しているため非常に時間がかかっていた

MATLABでは行列やベクトルに関する操作に最適化されているため,要素ごとに計算するforやifを使うとMATLABの本領を発揮できない*1*2

そこで,forとifを行列で表現することで並列処理しやすくした

高速化の主な方針

forは行列で一括に

行列の中で隣り合った要素の値から計算するとき,行列の要素ごとに計算するのではなく,4つの方向にずらした行列を予め用意すればよい

例えば,行列の要素を1つ右にずらして元の行列に足すことは,元の行列のすべての要素(1列目を除く)に対して,1つ左の要素の値を足していることに等しい

行列の要素をずらすときは関数(circshift()*3)を用いることで実現できる

ifはLogical行列との積で

行列の境界部分など,行列のある要素だけ違う計算をして欲しいときは,次のようにベクトル化する

  1. 予めこの要素は-1,他は-2というようなフラグ行列 (flag) を用意する

  2. 場合分けした全ての計算方法を行列の全ての要素に適用し,各場合での計算結果の行列(X_1, X_2とする)を用意する

  3. flag==-1とすると-1とフラグを立てた要素のみが1,他の要素が0となる行列(Logical行列)を作れる

  4. 3より,X=(flag==-1).*X_1+(flag==-2).*X_2とすれば求めたい行列が求まる.".*"とは,行列の要素ごとの積算を意味している

計算に使わない要素も計算しているため,十分に並列処理ができていないとかえって遅くなるので注意する

ループの中は最低限

これはどのプログラミング言語にも言える

今回の場合は電位分布の行列を更新していく中で,各点でのz軸からの距離(r座標)は変わらないためwhile ループの外に出すという改善をした

MATLABで用意されいる関数を使う

MATLABは非常に多くの関数を提供している

何かしたいことがあれば「MATLAB 行列 要素 ずらす」というように毎度ネットで検索すれば便利な関数が見つかる

頑張ってphi_l=[zeros(1,M);phi(end-1,:))]; なんて手作りする必要はない

用意されいる関数は最適化されているのでMATLABの本領を発揮できる

コードのbefore & after

高速化したコードでしていることを簡単に説明すると,非境界部分は隣り合った4つの点から,境界部分は隣り合った3つの点から電位を計算している

それを何度も繰り返していくうちに値が収束し(前回の電位分布とあまり変わらなくなり),電位分布が求まる

本記事ではコードの高速化を目的としているため,電位の計算の原理については前回の記事を参考されたい

改善前

$ N\times M $行列の各要素について,場合分けして計算していた

原理を忠実に表しているため理解しやすいが遅い

    for i = 1:N
        for j = 1:M 
            r = i*dr; %格子点(r,z)のrの値
            %非境界のphi
            if (flag(i,j)==-10)
                phi(i,j)=0.25*(phi(i+1,j)+phi(i-1,j)+phi(i,j+1)+phi(i,j-1)+0.5*dr*(phi(i+1,j)-phi(i-1,j))/r);
            %z方向ノイマン境界
            elseif (flag(i,j)==-1) 
                phi(i,j)=0.25*(phi(i+1,j)+phi(i-1,j)+2*phi(i,j-1)+0.5*dr*(phi(i+1,j)-phi(i-1,j))/r);
            %r方向ノイマン条件
            elseif (flag(i,j)==-2)
                phi(i,j)=0.25*(2*phi(i+1,j)+phi(i,j+1)+phi(i,j-1));
            end
        end
    end

改善後

電位分布の計算を繰り返し行う前に,点 (r, z) でのr の値を先に行列にして保存しておく

r=repmat((dr:dr:N*dr)',1,M);%点(r,z)のrの値の行列

電位分布の行列の値を更新していくループの中のコード

見た目はきれい

行列演算に慣れていないと理解に時間がかかりそう

実行される行はたったの8行で,改善前の12行よりさらに簡潔になっている

    %ずらした行列を生成
    %行列のr行z列目の要素が点 (r,z) に対応していることに注意
    phi_l=circshift(phi,1,1); %r方向に1ずらした行列
    phi_h=circshift(phi,-1,1); %-r方向に1ずらした行列
    phi_j=circshift(phi,-1,2); %-z方向に1ずらした行列
    phi_k=circshift(phi,1,2); %z方向に1ずらした行列

    %各条件での計算結果を予め作成
    %非境界点
    phi_1=0.25*((phi_l+phi_h+phi_j+phi_k)+0.5*dr*(phi_h-phi_l)./r);    
    %z方向ノイマン条件
    phi_2=0.25 *((phi_l+phi_h+2*phi_k)+0.5*dr*(phi_h-phi_l)./r);     
    %r方向ノイマン条件
    phi_3=0.25*(2*phi_h+phi_j+phi_k);
    
    %flagは非境界部分の点:-10, z方向ノイマン条件:-1, r方向ノイマン条件:-2
    %ディリクレ条件 (flag: -3) の点の電位は固定
    phi=(flag==-10).*phi_1+(flag==-1).*phi_2+(flag==-2).*phi_3+(flag==-3).*Prev_phi;

計算時間の計測結果

400行300列

r方向の格子数400,z方向の格子数300で改善前と改善後のプログラムをプロファイラーで計測しながら実行した結果,下のような結果になった

改善前は2136.810秒,改善後は110.447秒

実行時間は20分の1にまで減らすことが出来ている

f:id:monhime:20191219210542p:plainf:id:monhime:20191219210628p:plain

計測しない場合は改善前は3分56秒,改善後は1分50秒で,プロファイラーで計測するための割り込みっぽいのがない分実行時間時間は速くなった

やはり改善後の方が速くなっていた

出力

どちらのコードでも大きな違いは見られなかった

改善したコードが計算を間違えている訳ではないことを確認できた

差分法による電界分布計算 差分法による電界分布計算

800行800列

さらに行列のサイズを増やしてプロファイラーではなくストップウォッチで計測した結果,改善前は42分42秒であるのに対し,改善後は22分27秒と約20分速くなった

計算時間
改善前 42分42秒
改善後 22分27秒

余談

改善後のコードのうち行列の要素をずらす部分で,初めは関数circshift()を使わずphi_l=[zeros(1,M);phi(end-1,:))]というように行列を組み合わせて定義していた

すると,全体の計算時間が改善前の場合と同じくらいかかった

最適化されているMATLABのコードをわざわざ他の方法で実装するのはナンセンス

改善後のソースコード

差分法による電位分布計算のコードの改良版