Fortran、oneMKL、OpenMP*を使用して LU 因数分解を高速化

2023.03.29
インテル® oneMKL
FortranOpenMPParallel Universe

この記事は、The Parallel Universe Magazine 51 号に掲載されている「Accelerate Lower Upper (LU) Factorization Using Fortran, oneMKL, and OpenMP*」の日本語参考訳です。原文は更新される可能性があります。原文と翻訳文の内容が異なる場合は原文を優先してください。

parallel_v51_04

BASIC、Pascal、C は私が最初に学んだプログラミング言語ですが、初めて本格的な (教室外での) コーディングをしたのは Fortran でした。インテルに入社して C と C++ に乗り換えるまでの 12年間は、ほとんど Fortran しか使いませんでした。この 20年間、Fortran を 1 行も書いていないことを、私はいくつかの著作で嘆いています。旧友の Jim Cownie が「Fortran は過去の言語なのか?」 (英語) において、そうではないという説得力のある答えを出すまで、Fortran がまだ使われているかどうか疑問に思っていたくらいです。そのため、同僚から最新の Fortran コンパイラーの OpenMP* アクセラレーター・オフロード・サポートを試してみないかと誘われたとき、私はそのチャンスに飛びつきました。

私は、Fortran と同じくらい OpenMP* が好きです。OpenMP* は常にコードを並列化する簡単な方法であり、1997年に並列コンピューティング分野に登場して爆発的に広がって以来、OpenMP* 仕様は長い道のりを歩んできました。現在では、アクセラレーターへのオフロードもサポートしています。そこで、いつも使っている SYCL* ではなく、並列プログラミングの原点に戻って、Fortran、OpenMP*、インテル® oneAPI マス・カーネル・ライブラリー (インテル® oneMKL) を使ってアクセラレーターに計算をオフロードする方法を説明したいと思います。

以前の記事 (記事 (英語)、ウェビナー(英語)、サンプルコード (英語)) を読んでいただいた方は、私がフーリエ相関アルゴリズムを気に入っていることをご存じだと思います。私はこれを生体分子のドッキングの研究に使い、最近では 2D 画像と 3D 画像の調整に使用しています。これまで、このアルゴリズムにはうんざりするくらい取り組んできたので、ここでは線形代数計算 (特に LU 分解) をアクセラレーターにオフロードすることがいかに簡単か見てみましょう。LU 分解について詳しく説明することがこの記事の目的ではないので、線形方程式を解くために使われるものだと理解いただければ十分です。詳しく知りたい場合は、インテル® MKL クックブック (英語) の「LU 因数分解されたブロック三重対角係数行列を含む連立線形方程式を解」 と「一般的なブロック三重対角行列の因数分解」を参照してください。ここでは、Fortran、インテル® oneMKL、および OpenMP* の実装を通して、オフロード構文を説明し、ヘテロジニアス並列処理におけるホストとデバイス間のデータ転送に注目します。まず、インテル® oneMKL に含まれている dgetri_oop_batch_strided.f90 の例から開始し (「サンプルコードの使用」)、効率を改善するためこれを変更します。

このプログラムは、まず OpenMP* オフロードをサポートする Fortran インターフェイスをインクルードし、次に 32 ビット整数型と 64 ビット整数型のどちらを使用するかを決定します (詳細は、「ILP64 インターフェイスと LP64 インターフェイスの使用」 (英語) を参照) (図 1)。この判断は、この記事の最後にあるコンパイラー・コマンドで示すように、コンパイル時に行います。セットアップを完了するため、因数分解する行列の数とそのサイズを設定し、インテル® oneMKL サブルーチンが必要とする次元とストライドのパラメーターを設定し、ワーク配列を割り当てます。これで計算を開始する準備が整いました。

include "mkl_omp_offload.f90"

program sgetri_oop_batch_strided_example

#if defined(MKL_ILP64)
   use onemkl_lapack_omp_offload_ilp64   ! 64-bit integer type
#else
   use onemkl_lapack_omp_offload_lp64    ! 32-bit integer type
#endif

   integer,   parameter :: n = 10, batch_size = 4
   integer              :: lda, ldainv, stride_a, stride_ainv, stride_ipiv
   real,    allocatable :: a(:,:), ainv(:,:)
   integer, allocatable :: ipiv(:,:), info(:)

   lda         = n
   ldainv      = n
   stride_a    = lda    * n
   stride_ainv = ldainv * n
   stride_ipiv = n

   allocate(a(stride_a, batch_size), ainv(stride_ainv, batch_size), &
            ipiv(stride_ipiv, batch_size), info(batch_size))

図 1. サンプルプログラムのセットアップ。OpenMP* ターゲットオフロードをサポートするインテル® oneMKL ヘッダーとモジュールは、青色で表示しています。ハードコードされた行列の次元とバッチサイズのパラメーターは、テストには適していますが、高速化には小さすぎるということに注意してください。

手元にある GPU は単精度の計算しかサポートしていないので、Fortran REAL 型と適切な LAPACK (英語) サブルーチンを使用して、LU 因子分解 (sgetrf_batch_strided と LU 因子分解した行列のアウトオブプレース反転 (sgetri_oop_batch_strided) を計算します (倍精度に変換するには、Fortran DOUBLE PRECISION 型に切り替え、LAPACK の dgetrf_batch_strided サブルーチンと dgetri_oop_batch_strided サブルーチンを使用するだけです)。オリジナルのコードでは、2 つの OpenMP* target data 領域を使用しています (図 2)。ターゲット構造は、デバイスに制御を転送します。最初の領域は、入力行列をデバイスメモリーに転送し [map(tofrom:a)] 、インプレース LU 分解 (sgetrf_batch_strided) をデバイスにディスパッチし、LU 分解した行列 (a に格納)、ピボット・インデックス [map(from:ipiv)] 、ステータス [map(from:info)] をデバイスメモリーから取得します。因数分解が正常に実行されると、プログラムは次の OpenMP* 領域に進みます。LU 因数分解された行列とピボット・インデックスはデバイスメモリーに戻され [map(to:a) と map(to:ipiv)] 、アウトオブプレース反転はデバイス上で計算され (sgetri_oop_batch_strided)、逆行列 [map(from:ainv)] とステータス [map(from:info)] はデバイスメモリーから取得されます。

! OpenMP* オフロードで LU 因数分解を計算
! 開始時に A は入力行列を格納
! 終了時に A は LU 因数分解された行列を格納
!$omp target data map(tofrom:a) map(from:ipiv) map(from:info)

    !$omp dispatch
    call sgetrf_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: sgetrf_batch_strided returned with errors.'
else
    ! OpenMP* オフロードで行列の反転を計算し、終了時に逆行列は Ainv に格納
    !$omp target data map(to:a) map(to:ipiv) map(from:ainv) map(from:info)

        !$omp dispatch
        call sgetri_oop_batch_strided(n, a, lda, stride_a, ipiv, stride_ipiv, &
                                      ainv, ldainv, stride_ainv, batch_size, info)
    !$omp end target data
endif

if (any(info .ne. 0)) then
    print *, 'Error: sgetri_oop_batch_strided returned with errors.'
else
    print '("The calculations executed successfully.")'
endif

図 2. 2 つの LAPACK サブルーチンそれぞれに OpenMP* ターゲット領域を使用して、LAPACK 計算をアクセラレーターにオフロード。OpenMP* ディレクティブは青で表示しています。LAPACK の呼び出しは緑色で表示しています。

続きはこちら
タイトルとURLをコピーしました