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

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

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

コードが遅い具体的な原因
高速化の主な方針
コードのbefore & after
- 改善前
- 改善後
計算時間の計測結果
改善後のソースコード
次回

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

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

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

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

高速化の主な方針

forは行列で一括に

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

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

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

ifはLogical行列との積で

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

予めこの要素は-1，他は-2というようなフラグ行列 (flag) を用意する
場合分けした全ての計算方法を行列の全ての要素に適用し，各場合での計算結果の行列（X_1, X_2とする）を用意する
flag==-1とすると-1とフラグを立てた要素のみが1，他の要素が0となる行列（Logical行列）を作れる
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;