3

この質問は、 MATLAB が Numpy の 2 倍の速さであるという最近の質問へのフォローアップです。

私は現在、MATLAB と Numpy の両方に実装された Gauss-Seidel ソルバーを持っています。これは、2D 軸対称ドメイン (円柱座標) に作用します。コードはもともと MATLAB で作成され、その後 Python に移行されました。MATLAB コードは ~20 秒で実行されますが、Numpy コードは ~30 秒かかります。Numpy を使用したいのですが、このコードはより大きなプログラムの一部であるため、シミュレーション時間がほぼ 2 倍長くなることは重大な欠点です。

このアルゴリズムは、離散化されたラプラス方程式を長方形のメッシュ (円筒座標) で単純に解きます。メッシュの更新間の最大差が指定された許容値より小さくなると終了します。

Numpy のコードは次のとおりです。

import numpy as np
import time

T = np.transpose

# geometry
length = 0.008
width = 0.002

# mesh
nz = 256
nr = 64

# step sizes
dz = length/nz
dr = width/nr

# node position matrices
r = np.tile(np.linspace(0,width,nr+1), (nz+1, 1)).T
ri = r/dr

# equation coefficients
cr = dz**2 / (2*(dr**2 + dz**2))
cz = dr**2 / (2*(dr**2 + dz**2))

# initial/boundary conditions
v = np.zeros((nr+1,nz+1))
v[:,0] = 1100
v[:,-1] = 0
v[31:,29:40] = 1000
v[19:,54:65] = -200

# convergence parameters
tol = 1e-4

# Gauss-Seidel solver
tic = time.time()
max_v_diff = 1;
while (max_v_diff > tol):

    v_old = v.copy()

    # left boundary updates
    v[0,1:nz] = cr*2*v[1,1:nz] + cz*(v[0,0:nz-1] + v[0,2:nz+2])
    # internal updates
    v[1:nr,1:nz] = cr*((1 - 1/(2*ri[1:nr,1:nz]))*v[0:nr-1,1:nz] + (1 + 1/(2*ri[1:nr,1:nz]))*v[2:nr+1,1:nz]) + cz*(v[1:nr,0:nz-1] + v[1:nr,2:nz+1])
    # right boundary updates
    v[nr,1:nz] = cr*2*v[nr-1,1:nz] + cz*(v[nr,0:nz-1] + v[nr,2:nz+1])
    # reapply grid potentials
    v[31:,29:40] = 1000
    v[19:,54:65] = -200

    # check for convergence
    v_diff = v - v_old
    max_v_diff = np.absolute(v_diff).max()

toc = time.time() - tic
print(toc)

これは実際に私が使用する完全なアルゴリズムではありません。完全なアルゴリズムでは、逐次過緩和とチェッカー盤反復スキームを使用して速度を向上させ、ソルバーの方向性を取り除きますが、単純にするために、このわかりやすいバージョンを提供しました。Numpy の速度の欠点は、フル バージョンではより顕著です (Numpy と MATLAB でそれぞれ 17 秒対 9 秒のシミュレーション時間)。

v を列優先の配列に変更して、前の質問の解決策を試しましたが、パフォーマンスは向上しませんでした。

助言がありますか？

編集: 参照用の Matlab コードは次のとおりです。

% geometry
length = 0.008;
width = 0.002;

% mesh
nz = 256;
nr = 64;

% step sizes
dz = length/nz;
dr = width/nr;

% node position matrices
r = repmat(linspace(0,width,nr+1)', 1, nz+1);
ri = r./dr;

% equation coefficients
cr = dz^2/(2*(dr^2+dz^2));
cz = dr^2/(2*(dr^2+dz^2));

% initial/boundary conditions
v = zeros(nr+1,nz+1);
v(1:nr+1,1) = 1100;
v(1:nr+1,nz+1) = 0;
v(32:nr+1,30:40) = 1000;
v(20:nr+1,55:65) = -200;

% convergence parameters
tol = 1e-4;
max_v_diff = 1;

% Gauss-Seidel Solver
tic
while (max_v_diff > tol)
    v_old = v;

    % left boundary updates
    v(1,2:nz) = cr.*2.*v(2,2:nz) + cz.*( v(1,1:nz-1) + v(1,3:nz+1) );
    % internal updates
    v(2:nr,2:nz) = cr.*( (1 - 1./(2.*ri(2:nr,2:nz))).*v(1:nr-1,2:nz) + (1 + 1./(2.*ri(2:nr,2:nz))).*v(3:nr+1,2:nz) ) + cz.*( v(2:nr,1:nz-1) + v(2:nr,3:nz+1) );    
    % right boundary updates
    v(nr+1,2:nz) = cr.*2.*v(nr,2:nz) + cz.*( v(nr+1,1:nz-1) + v(nr+1,3:nz+1) );
    % reapply grid potentials
    v(32:nr+1,30:40) = 1000;
    v(20:nr+1,55:65) = -200;

    % check for convergence
    max_v_diff = max(max(abs(v - v_old)));

end
toc
4

2 に答える 2

2

私のラップトップでは、コードは約 45 秒で実行されます。事前に割り当てられた作業配列の再利用を含め、中間配列の作成を最小限に抑えようとすることで、その時間を 27 秒に短縮することができました。これにより、MATLAB のレベルに戻るはずですが、コードは読みにくくなります。とにかく、以下のコードを見つけて、# Gauss-Seidel solverコメントの下のすべてを置き換えます。

# work arrays
v_old = np.empty_like(v)
w1 = np.empty_like(v[0, 1:nz])
w2 = np.empty_like(v[1:nr,1:nz])
w3 = np.empty_like(v[nr, 1:nz])

# constants
A = cr * (1 - 1/(2*ri[1:nr,1:nz]))
B = cr * (1 + 1/(2*ri[1:nr,1:nz]))

# Gauss-Seidel solver
tic = time.time()
max_v_diff = 1;
while (max_v_diff > tol):

    v_old[:] = v

    # left boundary updates
    np.add(v_old[0, 0:nz-1], v_old[0, 2:nz+2], out=v[0, 1:nz])
    v[0, 1:nz] *= cz
    np.multiply(2*cr, v_old[1, 1:nz], out=w1)
    v[0, 1:nz] += w1
    # internal updates
    np.add(v_old[1:nr, 0:nz-1], v_old[1:nr, 2:nz+1], out=v[1:nr, 1:nz])
    v[1:nr,1:nz] *= cz
    np.multiply(A, v_old[0:nr-1, 1:nz], out=w2)
    v[1:nr,1:nz] += w2
    np.multiply(B, v_old[2:nr+1, 1:nz], out=w2)
    v[1:nr,1:nz] += w2
    # right boundary updates
    np.add(v_old[nr, 0:nz-1], v_old[nr, 2:nz+1], out=v[nr, 1:nz])
    v[nr, 1:nz] *= cz
    np.multiply(2*cr, v_old[nr-1, 1:nz], out=w3)
    v[nr,1:nz] += w3
    # reapply grid potentials
    v[31:,29:40] = 1000
    v[19:,54:65] = -200

    # check for convergence
    v_old -= v
    max_v_diff = np.absolute(v_old).max()

toc = time.time() - tic
于 2013-07-10T21:53:24.253 に答える
2

次のプロセスに従って、ラップトップでの実行時間を 66 秒から 21 秒に短縮できました。

  1. ボトルネックを見つけます。IPythonコンソールからline_profilerを使用してコードをプロファイリングし、最も時間がかかった行を見つけました。時間の 80% 以上が「内部更新」を行う行に費やされていることが判明しました。

  2. 最適化する方法を選択します。numpy でコードを高速化するためのツールがいくつかあります (Cython、numexpr、weave...)。特に、 scipy.weave.blitz問題のある行などの numpy 式を高速なコードにコンパイルするのに適しています。理論的には、その行をラップし"..."て実行することができますweave.blitz("...")が、更新されている配列が計算で使用されるため、ドキュメントのポイント 4 で述べられているように、一時配列を使用して同じ結果を保持する必要があります。

    expr = "temp = cr*((1 - 1/(2*ri[1:nr,1:nz]))*v[0:nr-1,1:nz] + (1 + 1/(2*ri[1:nr,1:nz]))*v[2:nr+1,1:nz]) + cz*(v[1:nr,0:nz-1] + v[1:nr,2:nz+1]); v[1:nr,1:nz] = temp"
temp = np.empty((nr-1, nz-1))
...
while ...
    # internal updates
    weave.blitz(expr)

  3. 結果が正しいことを確認した後、 を使用してランタイム チェックを無効にしweave.blitz(expr, check_size=0)ます。コードは 34 秒で実行されます。

  4. Jaime の作業に基づいて、定数係数AとB式を事前に計算します。コードは 21 秒で実行されます (変更は最小限ですが、コンパイラが必要になります)。

これはコードの核心です：

from scipy import weave

# [...] Set up code till "# Gauss-Seidel solver"

tic = time.time()
max_v_diff = 1;
A = cr * (1 - 1/(2*ri[1:nr,1:nz]))
B = cr * (1 + 1/(2*ri[1:nr,1:nz]))
expr = "temp = A*v[0:nr-1,1:nz] + B*v[2:nr+1,1:nz] + cz*(v[1:nr,0:nz-1] + v[1:nr,2:nz+1]); v[1:nr,1:nz] = temp"
temp = np.empty((nr-1, nz-1))
while (max_v_diff > tol):
    v_old = v.copy()
    # left boundary updates
    v[0,1:nz] = cr*2*v[1,1:nz] + cz*(v[0,0:nz-1] + v[0,2:nz+2])
    # internal updates
    weave.blitz(expr, check_size=0)
    # right boundary updates
    v[nr,1:nz] = cr*2*v[nr-1,1:nz] + cz*(v[nr,0:nz-1] + v[nr,2:nz+1])
    # reapply grid potentials
    v[31:,29:40] = 1000
    v[19:,54:65] = -200
    # check for convergence
    v_diff = v - v_old
    max_v_diff = np.absolute(v_diff).max()
toc = time.time() - tic
于 2013-07-10T23:35:11.677 に答える