8. 連立一次方程式

8.1 NumPy,SciPy,CuPyのソースコード

リスト8-1にNumPy,SciPy,CuPyを用いて連立一次方程式を解くソースコードを示します。
NumPyとSciPyはCPU用、CuPyはGPU用です。 いずれもlinalg.solve関数を使用しています。
実数と複素数、単精度と倍精度を選択することができます。
行列は乱数で作成した非対称密行列です。
CuPyではデバイスメモリーを確実に解放するためにmemory poolを使用しています。

リスト8-1 NumP,SciPy,CuPyを用いて連立一次方程式を解くソースコード (leq.py)

"""
leq.py
solve linear equations
NumPy, SciPy, CuPy
"""

import numpy as np
import scipy as sp
import cupy as cp
import time

# parameters
N = 5000  # matrix size
L = 1  # repeat
dkind = 1      # 1=real, 2=complex
solver = 1     # 1=NumPy, 2=SciPy, 3=CuPy
precision = 1  # 1=single, 2=double

# dtype
f_dtype = 'f4' if precision == 1 else 'f8'
c_dtype = 'c8' if precision == 1 else 'c16'

# timer
t0 = time.time()

# real
if dkind == 1:
    a = np.random.rand(N, N).astype(f_dtype)
    b = np.random.rand(N).astype(f_dtype)
    x = np.zeros(N, f_dtype)
# complex
elif dkind == 2:
    a = np.random.rand(N, N).astype(f_dtype) \
      + np.random.rand(N, N).astype(f_dtype) * 1j
    b = np.random.rand(N).astype(f_dtype) \
      + np.random.rand(N).astype(f_dtype) * 1j
    x = np.zeros(N, c_dtype)

# well-condition
for i in range(N):
    a[i, i] *= 10

# alloc device memory
if solver == 3:
    mempool = cp.get_default_memory_pool()
    # copy to device
    d_a = cp.array(a)
    d_b = cp.array(b)
    d_x = cp.array(x)

# timer
t1 = time.time()

# solve
for _ in range(L):
    if   solver == 1:
        # NumPy (CPU)
        x = np.linalg.solve(a, b)
    elif solver == 2:
        # SciPy (CPU)
        x = sp.linalg.solve(a, b)
    elif solver == 3:
        # CuPy (GPU)
        d_x = cp.linalg.solve(d_a, d_b)

# timer
t2 = time.time()

# free device memory
if solver == 3:
    # copy to host
    x = cp.asnumpy(d_x)
    # free
    d_a = None
    d_b = None
    d_x = None
    mempool.free_all_blocks()

# check
res = np.dot(a, x) - b
ans = np.linalg.norm(res)

# output
device = 'CPU(NumPy)' if solver == 1 else 'CPU(SciPy)' if solver == 2 else 'GPU(CuPy)'
print('%s N=%d(%d), %.2f+%.2f[sec], err=%.3e'
    % (device, N, L, t1 - t0, (t2 - t1) / L, ans))

# free
a = None
b = None
x = None

8.2 実数の計算時間

表8-1に実数の連立一次方程式の計算時間を示します。
計算時間を正確に測定するために、 計算時間が短いときは多数回計算してその回数で割っています。

表8-1から以下のことがわかります。

• NumPyでは倍精度の計算時間は単精度より少し短い（理由不明）。
• SciPyでは倍精度の計算時間は単精度の約2倍である。
• CuPyでは単精度はNumPyとSciPyの約10倍速い。倍精度はNumPy,SciPyとほぼ同じである。 倍精度ではハードウェアの制限がかかっている影響が考えられる。
• 計算時間は行列の大きさNの3乗にほぼ比例する。

以上から、連立一次方程式を解くときの指針は以下のようになります。

• 単精度：計算時間の観点からは、CPUではSciPy、GPUではCuPyが適していますが、 残差を確認するとSciPyとCuPyはNumPyより2桁大きくなっており計算精度に注意が必要です。
  
• 倍精度：計算時間の観点からはCPUのNumPyが適しています。 倍精度の速いGPUではCuPyが適しています。
  

8.3 複素数の計算時間

表8-2に複素数の連立一次方程式の計算時間を示します。
表8-1と同じ傾向が見られます。 複素数の計算時間は実数の3～5倍になります。