この記事は、The Parallel Universe Magazine 52 号に掲載されている「Solve Linear Systems Using oneMKL and OpenMP Target Offloading」の日本語参考訳です。原文は更新される可能性があります。原文と翻訳文の内容が異なる場合は原文を優先してください。
前号の記事、「Fortran、oneMKL、OpenMP*を使用して LU 因数分解を高速化」では、LU 因数分解と後続の反転をアクセラレーターにオフロードする方法を紹介しました。OpenMP* オフロード領域を統合してホストとデバイス間のデータ転送を最小限に抑えることでオーバーヘッドを削減するヒントをいくつか示しましたが、連立線型方程式を解くには至りませんでした。これは読者の演習として残しておいたものですが、『oneAPI GPU 最適化ガイド』の「データ転送とメモリー割り当てを最小化する」のいくつかのテクニックを説明するためにここで実行します。
OpenMP* を使用してアクセラレーターにオフロードされた oneMKL の LU 因数分解とソルバー関数を使用してバッチ処理された線形システムのグループを解く方法を示します。逆行列の因数は実際にはほとんど使用しないため、このデモではそのステップを省略し、因数分解から解法に進みます。前のサンプルコードは Fortran で記述されていました。ここでも引き続き Fortran を使用します。
まず、いくつかの小さな線形システムを読み込んで正しく解いていることを確認します。古い線形代数の教科書に適切な例がいくつか見つかりました (図 1)。行列は正方で、LU 因数分解を使用して解くことができます。
図 1. oneMKL と OpenMP* ターゲットオフロードを使用して、これらの 2 つの線形システムをバッチで解きます。これらの小さな問題はテストには適していますが、高速化の恩恵を受けるには小さすぎます。
このプログラムは、前号の記事のプログラムと似ています。まず OpenMP* オフロードをサポートする Fortran インターフェイスをインクルードし、次に 32 ビット整数型と 64 ビット整数型のどちらを使用するかを決定します (リスト 1) (詳細は、「ILP64 インターフェイスと LP64 インターフェイスの使用」 (英語) を参照)。この判断はコンパイル時に行われます (次のコンパイラー・コマンドを参照)。図 1 のテスト問題は、それぞれ 1 つの右辺 (RHS) がある 2 つの 3×3 線形システムで構成されています。batch_size、n、nrhs、および stride 変数がそれぞれ設定され、a、b、および ipiv 配列が割り当てられて、バッチ処理されたシステムが保持されます。最後の 4 つの文は、2 つの行列とその右辺をそれぞれ a と b にロードします。
include “mkl_omp_offload.f90”
program solve_batched_linear_systems
! 32 ビット整数型と 64 ビット整数型のどちらを使用するか決定
#if defined(MKL_ILP64)
use onemkl_lapack_omp_offload_ilp64 ! 64 ビット
#else
use onemkl_lapack_omp_offload_lp64 ! 32 ビット
#endif
integer, parameter :: n = 3, batch_size = 2, nrhs = 1
integer :: lda, stride_a, stride_ipiv
integer :: ldb, stride_b
real (kind=8), allocatable :: a(:,:), b(:,:)
integer, allocatable :: ipiv(:,:), info(:)
integer allocstat
character (len = 132) :: allocmsg
lda = n
stride_a = n * lda
stride_ipiv = n
ldb = n
stride_b = n * nrhs
! 必要なメモリーを割り当て
allocate (a(stride_a, batch_size), b(n, batch_size*nrhs), &
ipiv(stride_ipiv, batch_size), &
info(batch_size), stat = allocstat, errmsg = allocmsg)
if (allocstat > 0) stop trim(allocmsg)
! 行列を初期化。Fortran は列優先言語であることに注意。
a(:,1) = (/2.0, 1.0, -6.0, 4.0, -4.0, -9.0, -4.0, 3.0, 10.0/)
a(:,2) = (/0.0, 2.0, 6.0, 0.0, 4.0, 5.0, 3.0, -1.0, 5.0/)
b(:,1) = (/12.0, -21.0, -24.0/) ! x = (-4, 2, -3)
b(:,2) = (/ 6.0, -10.0, -7.0/) ! x = (-2, -1, 2)
リスト 1. サンプルプログラムのセットアップ。OpenMP* ターゲットオフロードをサポートする oneMKL ヘッダーとモジュールは、青色で表示しています。
計算を開始する準備が整いました。2 つの OpenMP* ターゲット領域を使用して LU 因数分解と解を線形システムにディスパッチする基本的な実装から始めましょう (リスト 2)。サンプルを次のようにコンパイルして実行し、正しい結果が得られることを確認します。
$ ifx -i8 -DMKL_ILP64 -qopenmp -fopenmp-targets=spir64 -fsycl -free \
> lu_solve_omp_offload_ex1_small.F90 -o lu_solve_ex1_small \
> -L${MKLROOT}/lib/intel64 -lmkl_sycl -lmkl_intel_ilp64 -lmkl_intel_thread \
> -lmkl_core -liomp5 -lpthread -ldl
$ ./lu_solve_ex1_small
=================================
Solutions to the linear systems
=================================
-4.0000 2.0000 -3.0000
-2.0000 -1.0000 2.0000
! OpenMP* オフロードを使用して LU 因数分解を計算
! 開始時に a は入力行列を格納
! 終了時に a は LU 因数分解された行列を格納
!$omp target data map(tofrom:a) map(from:ipiv) map(from:info)
!$omp dispatch
call dgetrf_batch_strided(n, n, a, lda, stride_a, &
ipiv, stride_ipiv, batch_size, info)
!$omp end target data
if (any(info .ne. 0)) then
print *, ‘Error: getrf_batch_strided returned with errors’
else
! 線形システムを解く。終了時に解を b に格納。
!$omp target data map(to:a) map(to:ipiv) map(tofrom: b) map(from:info)
!$omp dispatch
call dgetrs_batch_strided(‘N’, n, nrhs, a, lda, stride_a, &
ipiv, stride_ipiv, &
b, ldb, stride_b, batch_size, info)
!$omp end target data
if (any(info .ne. 0)) then
print *, ‘Error: getrs_batch_strided returned with errors’
else
print *, ‘Computation executed successfully’
endif
endif
リスト 2. oneMKL と OpenMP* ターゲットオフロードを使用してバッチ処理された線形システムの基本的なソリューション。OpenMP* ディレクティブは青で表示しています。LAPACK の呼び出しは緑色で表示しています。
target 構造は、アクセラレーター・デバイスに制御を転送します。最初の領域は、入力行列をデバイスメモリーに転送し［map(tofrom:a)］、インプレース LU 因数分解をデバイスにディスパッチし (dgetrf_batch_strided)、LU 因数分解された行列 (a に格納)、ピボット・インデックス［map(from:ipiv)］、ステータス［map(from:info)］をデバイスメモリーから取得します。因数分解が正常に実行されると、プログラムは次の OpenMP* 領域に進みます。LU 因数分解された行列とピボット・インデックスをデバイスメモリーに戻し［map(to:a) および map(to:ipiv)］、線形システムをデバイスで解きます (dgetrs_batch_strided)。RHS は解ベクトルで上書きされ、計算［map(from:info)］のステータスとともにデバイスメモリー［map(tofrom:b)］から取得されます。ホストとデバイス間のデータ転送を図式化したのが図 2 です。
