IPSJ SIG Technical Report Vol.2015-HPC-148 No /3/2 CUDA-BLAS GPU 1,3,a) 1 2,3 2,3, GPU CUDA-BLAS,., GPU Eigen-G, MAGMA CUDA- BLAS., CUDA-BLAS AS

CUDA-BLAS

等の選択による

最速

GPU

固有値ソルバーの性能評価

今村 俊幸

1,3,a)

_{椋木 大地}

1

_{山田 進}

2,3

_{町田 昌彦}

2,3 概要： 本研究では,既存の固有値計算で最速となるGPU固有値ソルバーに対してCUDA-BLASの選択による性 能の違いを検証し,最速ソルバーの組み合わせを発見することを目的とした議論と実験を行う.論点として, 既存のGPU固有値ソルバーEigen-G, MAGMAに対する分析とその内部で共通に使用されている CUDA-BLASの実装について最良の組み合わせを見つけ出すことを考える. 本報告では, CUDA-BLASの実装には 我々が開発を進めてきたASPEN.K2のdsymvが固有値ソルバーの性能に寄与し, MAGMA+ASPEN.K2 の実装による更なる高性能化の一例を示す.

1.

はじめに

我々はこれまでにGPU向けの高性能BLASを開発して きており([1], [2]など), CUDA環境[3]での高度な最適化 に加えて自動チューニング技術の適用によりCUBLAS[4] やMAGMABLAS[5]を凌ぐ性能を記録している. その中 でも,対称(エルミート)行列ベクトル積((SY|HE)MV)は HPC-138[6], HPC-146[7]と本研究会にて開発経過を報告 してきている. (SY|HE)MVは次式に示すものであり. 数値線形代数計 算の中でも重要な位置を占める計算パターンといえる.

y := αA‘U’ or ‘L’x + βy

where A(= A∗)∈ Kn×n, x∈ Kn, K = R or C. (1) 例えば, SYMVは対称密行列の固有値計算の前処理である ハウスホルダー三重対角化においてランク更新計算と並ぶ 高負荷計算部分であり,多くの場合SYMVの性能が対称行 列の固有値ソルバーの性能を決めるとも言われている. 高 性 能 計 算 を 志 向 し た GPUの 高 い 性 能 を 利 用 し た GPGPUにおいても同様のことが成り立ち, SYMVの最適 化が固有値ソルバーの性能を左右してきた. これまでに多 くの実装技術が報告されていることがそれを実証してい 1 _{理化学研究所 計算科学研究機構}

RIKEN Advanced Institute for Computational Science, Kobe, Hyogo

2 _{日本原子力研究開発機構}

Japan Atomic Energy Agency, Kashiwa, Chiba 3 _{科学技術振興機構CREST}

CREST JST, Kawaguchi, Saitama a) _{[email protected]}

る. GPUを使用する固有値ソルバーとしてMAMGA[5]の

magma dsyevdやmagma dsyevdx 2stage,更にEigen-G[8]

などがある. これらは, フルにGPU上でのコード開発が なされているわけでなく,ホストCPU側から適切な部分 をGPUにオフロードする形の実装形態をとる. 具体的に

はLevel2, Level3の主要なBLASをGPUにオフロードす

る形で開発されている. 従って,適切なメモリ管理や初期 化処理等がなされていれば, CUDA-BLASを自由に他の実 装に置き換えることが可能である. また, MAGMA並びに Eigen-Gはホスト側のCPUリソースに対するケアもされ ており, CPUとGPUの計算能力を最大限に活かすことが 可能である. つまり, GPU上の線形計算ライブラリは下位 のCUDA-BLASカーネルを選択したりCPU資源への計 算割付をする中位のドライバールーチンで決まることとな り, CPU+GPU環境に応じた最適な組み合わせが存在する ことになる. 本報告の目的は, CUDA-BLASならびにGPU固有値ソ ルバー先行研究に最新の技術を適用し最速となるライブラ リの組み合わせを調査しその性能測定結果の速報をまとめ ることである.

2.

CUDA-BLAS 実装

2.1 CUBLAS

CUBLAS[4]はNVIDIA社が開発するCUDA SDK[3]

に標準で提供され ているBLASのCUDA 実装である.

NVIDIA社が長い年月をかけLevel1からLevel3まで（ほ

ぼ)全てのカーネルをサポートし, 最適化されたオープン ソースを取り込む形で成長している. 現在では, CUDAで

のBLAS実装ではほぼ最速であることが知られている. 特 に, DGEMMについてはGPUの実効最高性能を測る目 的もあろうが, 後述する他のCUDA実装よりも高速であ る. 例えば,本実験で使用するTesla K20cでは倍精度計算 (DGEMM)において1TFLOPSを超える. 2.2 MAGMABLAS MAGMABLASはテネシー大学で開発されているGPU やメニイコアアクセラレータを混載する環境下での数値線 形計算をサポートするMAGMAライブラリ[5]に含まれ ているBLASの実装である. 原稿執筆現在でバージョン 1.6.1がリリースされている. MAGMABLASはCUBLAS の様にすべてのカーネルで最新の最適化がなされているわ けではないが, 最適化技術の会議報告は多数あり([9]ほか 多数), CUDA研究の初期段階では唯一の数値線形代数向 けオープンソースということもあり基本的なCUDA上で の並列アルゴリズムのリファレンスといえる. 本研究でも 幾つかの実装技術を参考にしている. 2.3 KBLAS

KBLAS[10], [11]はKAUSTで開発されているCUDA

BLASの一実装である. Level2であるGEMVとSYMV系 統をサポートし, 我々と同様にアトミック演算を使用した 実装によりSYMVで高い性能を達成している. 現在,バー ジョン1.2(1.3-betaも存在)がリリースされている. 現在, MAGMAが上位のライブラリ関数群の開発を精力的に進 めるのに対して, KBLASは下位の実装を進めるという棲 み分けになっている. 2.4 ASPEN.K2 ASPEN.K2[1]は 今 村 ら が 開 発 を 進 め て き た CUDA-BLASである. もともとは, GEMVの周期的な性能の 揺れを解決する高性能なライブラリの開発に主眼を置いて きたが,近年はSYMVの高性能化を中心に研究開発を進め ている([6], [7]など). SYMVの実装の特徴については次節 以降で詳しく説明を行っていく. 2.5 その他 上記に挙げた他にもいくつかの実装例は存在しており, 商用で最も有名なものとしてEM Photonics社が開発する CULA[16]に梱包されているCULABLASがある. 商用と して古くからの老舗的な位置づけであるが, BLAS内部に も一部手を入れているようで,過去の報告ではCUBLAS の未チューニング部分を独自実装により補っているものと 考えられる（ただしCUBLASとはライブラリ名レベルで 衝突するため同時に利用することはできない).

他にも, GLAS[14],が存在している. GLASはSørenssen

らGPUlab DTUで開発されたLevel1, Level2に対して自

動チューニングによって高い性能を出す最適化がなされ ている. 椋木[2]によるGEMVに対する自動最適パラメ タ選定機能(自動チューニングの一種)の研究も CUDA-BLASの高度な実装技術の一例といえる. 本研究報告では

CUBLAS, MAGMA, KBLAS, ASPEN.K2を使用する.

3.

CUDA 上の固有値ソルバー実装

CUDA上でGPUを使用する固有値ソルバーとして以下

3種類について説明する. 本報告ではこれらのうち商用で

はないMAGMAとEigen-Gを性能評価に用いる.

3.1 CULA

CULA[16]はLAPACKのCUDA実装の商用版である.

固有値ソルバーとしてサポートするのはQR法,二分法+逆 反復法に基づく計算ルーチンsyev, syevxが主である. バ ンド行列フォーマットにも対応するなど商用として幅広い 関数をサポートしている. しかしながら,後述のMAGMA やASPEN.K2で採用されている分割統治法をサポートし ていない. QR法や二分法は反復法であり, 分割統治法に 対して計算時間の面で大きく不利である. CULAのベンチ マーク情報[17]ではCPU単体に対して高い性能加速率を 示してはいるが,本研究での目的とは異なる視点での評価 といえる. また, 商用であるためソースコードを組み合わ せたり内部パラメータを細かく変更するなどのことはでき ない. 3.2 MAGMA

MAGMA[5]に収録される2サブルーチンmagma dsyevd

とmagma dsyevdx 2stageは実対称行列の固有値計算を行

うものである.

magma dsyevdはLAPACKに収録されるdsyevdに対応

しており,計算アルゴリズムは

( 1 )ハウスホルダー三重対角変換(magma dsytrd)

( 2 )分割統治法(magma dstedx)

( 3 )逆変換(magma dormtr)

の三ステップで構成される. ソースコードの分析から, 1)

ハウスホルダー変換(magma dsytrd)でdsymvとdsyr2k

をGPUにオフロードする構成になっている. 2) 分割統治 法(magma dstedx)では, 各統治フェーズにおける固有ベ クトルの相似変換で行列行列積カーネルdgemmをGPU にオフロードする. 3)逆変換(magma dormtr)では,コン パクトWY表現の構成まではCPUで実施し,固有ベクト ルへの作用が複数登場するため, dgemmをバッチモードで GPUにオフロードしている.

一方, magma dsyevdx 2stageはハウスホルダー三重対

角化による変換を行わず,同じハウスホルダー変換ではあ るがまず帯行列化を行ってから,三重対角行列に変換する という２ステップアルゴリズムを採用する. 固有ベクトル

計算までには以下の5ステップを必要とする.

( 1 )密_→帯変換(magma dsytrd sy2sb)

( 2 )帯の三重対角化(magma dsytrd sb2st)

( 3 )分割統治法(magma dstedx)

( 4 )三重_→帯の逆変換(magma dbuldge back)

( 5 )帯_→密の逆変換(magma dormqr 2stages)

magma dsytrd sy2sbは基本的な演算はハウスホルダー三

重対角化と同様であるが,計算の単位が1本のベクトルか ら複数本のベクトルの同時処理になる. 従って主要な処理

がdgemm, dsymm, dsyr2kとなりこれらをGPUにオフ

ロードしている. 三重対角化に比べて,オフロードの範囲

とLevel3 BLASのみになっている点が異なる. GPUにお

いてもLevel2よりもLevel3の方が基本性能は高い. 故に

sy2sbは性能向上が期待できる. magma dsytrd sb2stは複

数のバルジ追跡をホスト側のマルチコアで並列処理する.

magma dbuldge backも同様にホスト側でのマルチコアによ

るマルチスレッド並列処理による. magma dormqr 2stages

はAPIが異なるのみで基本的にはハウスホルダー逆変換 と同様である.

以上の様にMAGMAの固有値ソルバーはGPU, CPUを うまく使い分けた実装になっている.

3.3 Eigen-G

Eigen-G[8]は今村らが開発するEigenK, EigenExa[18]

のスタンドアローン版をGPU化したものである. 実装上 の技術要件は文献[8]に詳しく書かれている.

Eigen-Gが採用するアルゴリズムは, magma dsyevdと大

きな違いがなく,

( 1 )ハウスホルダー三重対角化

( 2 )分割統治法

( 3 )逆変換

の3ステップである. CPU・GPU間での非同期データ通

信の実現方法や, Eigen-GではCPUとGPUのdgemmを 同時に起動して計算システム全体の演算能力を使用するな どの違いがある. 2013年9月時点でのMAGMA1.4.0の実

装では, Eigen-Gはmagma dsyevdの計算時間の2/3程度

で固有値計算が可能であった. 当時の実装を分析すると, ( 1 ) DSYMVの性能差が際立っていた. MAGMAは内部 でCUBLASを呼んでおり, ASPEN.K2の3割程度の 性能しか出ていない. ( 2 ) async系の非同期通信よりも通信専用スレッドを立て た同期通信の方が有利であった点 などの理由が考えられた.

4.

固有値ソルバーの性能向上策

前節Eigen-Gの性能に関する分析の部分でも述べたよう に, GPU上の線形計算ソフトウェアの性能向上の最短手と しては, 既存の最良実装の組み合わせを考慮して置く必要 がある. 4.1 最速のCUDA-BLASの選択 MAGMAはバージョン1.6.1から上三角フォーマットの SYMVをサポートすることに伴い,内部のSYMV呼び出 しをCUBLASからMAGMABLASに変更している. これ により,従来のCUBLAS呼び出しよりも性能が倍以上よく なっている. この様に,呼び出し可能なCUDA-BLASの中 から最速な実装を呼び出すことがGPUの基本性能を出す ためには必要不可欠である. 表1ならびに表2はdgemm, dsymvの各CUDA-BLAS実装での性能をベンチマークし た結果である. 明らかにCUBLASがdgemmでの最速を 示している. GTX980では倍精度浮動小数点の性能が抑え られているため, MAGMABLASとの性能差は殆どないが 単精度浮動小数点(sgemm)では圧倒的な性能差があるこ とはよく知られている. 一方, dsymvはASPEN.K2がい ずれも最高性能を示しているがCUBLASのAtomic演算 版やKBLASもかなり近い性能である. これは, CUBLAS のAtomic演算版がKBLASを取り込んだものに対してパ ラメータチューニングをしたものであるからである. 一 方,幾つかの指摘もあるように, KBLASのAtomic演算の 実装方式ではAtomicAddの結果がスレッドの呼び出し順 序に依存してしまうためビットレベルでの丸め誤差相当 で結果が異なることが知られている. この問題について は, HPC-146においてmutexによる排他制御の工夫によ りASPEN.K2では起こらないことが報告されている. そ の意味でも, ASPEN.K2をdsymvに利用することが現時 点での最良の判断であるといえる. 更に詳細なdsymvの性能プロファイルを図1, 図2,図 3に示す. 少々の性能の揺れやスパイクは見受けられる が, 明らかにASPEN.K2は最高性能を示している. また,

ASPEN.K2以外は,行列のLower(下三角), Upper(上三角)

のフォーマットによって性能の違いが見受けられるが AS-PEN.K2にはほとんど違いはない. これらから,性能面か らのASPEN.K2の優位性は判断することはできる. また, ハウスホルダー三重対角化(dsytrdに相当する関数)では 行列のサイズを1ずつ減らしながら, dsymvを呼び出して いる. 問題サイズが小さいときには, dsymvの効率面での 差は大きく見えることはないが,問題が大きくなるとこれ は無視することはできない. ここで,固有値計算ライブラリにおけるCUDA-BLASの 選択方針は • dgemmはCUBLASを使用 • dsymvはASPEN.K2を使用する. なお, dsyr2kは内部的にはdgemm相当の実装になってい るため上記選択方針をとることにする.

表1 DGEMMの性能測定ベンチマーク結果(表中数字の単位はGFLOPS) 1088 2112 3136 4160 5184 6208 7232 8256 K20c CUBLAS 6.5 878.50 1014.45 1013.47 1024.86 1032.09 1040.63 1041.29 1042.92 MAGMABLAS 1.6.1 538.65 575.48 577.06 575.13 571.95 565.84 563.17 563.88 MKL 121.32 23.87 142.00 146.54 147.87 145.89 145.57 145.44 GTX980 CUBLAS 7.0 125.16 137.61 142.58 143.62 143.69 143.80 143.91 143.87 MAGMABLAS 1.6.1 121.34 126.19 141.74 142.75 142.67 142.89 143.02 142.96 MKL 47.18 21.77 66.08 65.72 60.91 64.71 65.72 69.71 表2 DSYMVの性能測定ベンチマーク結果(表中数字の単位はGFLOPS) 1088 2112 3136 4160 5184 6208 7232 8256 K20c CUBLAS 6.5 9.67 13.63 15.59 16.72 17.38 17.90 18.26 18.48 (Atomic) 17.56 34.22 46.62 55.04 56.23 56.94 59.21 59.54 MAGMABLAS 1.6.1 13.32 26.57 37.56 40.91 43.63 46.50 48.91 50.97 MKL 10.57 15.31 7.77 8.04 7.91 7.97 7.76 7.24 KBLAS 1.2 23.25 43.15 49.81 54.39 56.44 57.66 58.38 59.26 ASPEN.K2 1.5p2 26.13 46.71 54.88 58.82 59.21 61.73 62.39 62.64 GTX980 CUBLAS 7.0 16.35 26.47 30.13 30.91 33.29 34.01 34.83 35.47 (Atomic) 20.08 40.74 54.51 64.97 73.04 76.64 78.26 78.09 MAGMABLAS 1.6.1 17.81 30.46 39.04 44.66 47.53 49.64 51.29 53.22 MKL 4.46 4.30 3.98 4.00 4.03 4.15 3.85 4.08 KBLAS 1.2 23.75 48.53 63.16 73.13 77.62 79.27 79.75 79.93 ASPEN.K2 1.5p3 31.72 56.18 65.94 73.13 77.33 79.47 82.19 81.53

表3 Tesla K20cでの個々のルーチンの実測結果(秒) 1088 2122 3136 4160 MAGMA (1) trd 0.15 0.49 1.05 1.94 (2) ed 0.04 0.14 0.31 0.52 (3) tbk 0.01 0.07 0.14 0.28 Eigen-G (1) trd 0.11 0.33 0.74 1.42 (2) ed 0.04 0.15 0.24 0.65 (3) tbk 0.02 0.06 0.15 0.29 表4 GTX980での個々のルーチンの実測結果(秒) 1088 2122 3136 4160 MAGMA (1) trd 0.13 0.49 1.15 2.22 (2) ed 0.08 0.24 0.58 1.10 (3) tbk 0.03 0.18 0.56 1.22 Eigen-G (1) trd 0.11 0.38 0.86 1.69 (2) ed 0.06 0.28 0.73 1.39 (3) tbk 0.04 0.18 0.52 1.12 4.2 最速のサブルーチンの組合せ MAGMAとEigen-Gは基本的に同様のアルゴリズムを 採用しており, 中間データ保持などが共通化できれば相互 に最速なサブルーチンを組み合わせたより高速な固有値 ソルバーを構成することができる. ここでの議論には2ス テップアルゴリズム(magma dsyevdx 2stage)は相互に利 用可能な部分が少ない為, magma dsyevdでの議論を行う

こととする.

表3,表4に異なるプラットフォームでの固有値ソルバー の３ステップの実行時間の測定結果を示す. それぞれの環

境で,行列の次元を1088から4160次元と変化させた結果

である. Tesla K20c側はCPUもGPUも強力でGTX980

側はどちらも相対的に非力なものに相当する(あくまでも Tesla K20cとGTX980との倍精度浮動小数点の計算能力 での比較のもとでの話である). 固有値ソルバーの3ステッ プに共通するのはdgemmであり,その能力差が全体の実 測結果に反映されている形になるが, 共通していえること は(1)ではEigen-G, (2)はMAGMAがよりよい結果を示 している. (3)については環境に依存している. (1)のハウスホルダー三重対角化の主な違いは前節でも 示した, dsymvの実装によるところである. dgemmの性能 差はそれをひっくり返すことはできない規模であることに なる. (2)はMAGMAの実装が優れているということに尽 きる. Eigen-Gは単にdgemmをオフロードするのみであ り特別な実装を行ってはいない. (3)の違いを説明するに は更なる詳細な分析が必要であるが, MAGMAではCPU 側の能力をdgemmに生かしていないことによるものと 考えられる. Eigen-GではCPUも合わせた性能を出すた

め, GTX980でdgemmのGPU性能が140GFLOPS程度

であっても, CPUも合わせて210GFLOPS相当になるこ とから3:2程度の性能差が生まれていることが予想される. 0 10 20 30 40 50 60 70 80 90 0 2000 4000 6000 8000 10000 12000 14000 GFlops N DSYMV (GeForce GTX 980)

ASPEN.K2 1.5p3x DSYMV, Upper ASPEN.K2 1.5p3x DSYMV, Lower

0 10 20 30 40 50 60 70 80 90 0 2000 4000 6000 8000 10000 12000 14000 GFlops N DSYMV (GeForce GTX 980)

CUBLAS 7.0RC, DSYMV, Upper (atomics mode) CUBLAS 7.0RC, DSYMV, Lower (atomics mode)

0 10 20 30 40 50 60 70 80 90 0 2000 4000 6000 8000 10000 12000 14000 GFlops N DSYMV (GeForce GTX 980)

MAGMA 1.6.1, DSYMV, Upper MAGMA 1.6.1, DSYMV, Lower MAGMA 1.6.1, DSYMV, Upper (work) MAGMA 1.6.1, DSYMV, Lower (work)

0 10 20 30 40 50 60 70 80 90 0 2000 4000 6000 8000 10000 12000 14000 GFlops N DSYMV (GeForce GTX 980)

KBLAS 1.2, DSYMV, Upper KBLAS 1.2, DSYMV, Lower

図2 GeForce GTX980でのSYMVの性能(それぞれ行列は32次 元毎に測定)

表5 (補足)実験に使用したGPUの諸性能ならびにホストCPUのハードウェア/ソフトウェ ア環境

Tesla K20c GTX980 GPU Name GK110 GM204 Compute Capability 3.5 5.2

GPU Clock (MHz) 706(boost NA) 1126(boost 1216) Multiprocessors 13 16

CUDA Cores 2496 (=13*192) 2048 (=16*128) Memory Capacity (MByte) 5120 (GDDR5) 4096 (GDDR5) Memory Clock (MHz) 5200(320bit) 7012(256bit) Memory Bandwidth (GB/s) 208 224

ECC Support Enabled NA (ECC onで実行)

PCI bus PCIe 2.0×16 PCIe 3.0×16 (host側はPCIe2で接続)

Host (a) (b)

Host (a) Host (b) CPU AMD FX-8120 Intel Core i7-3930K

CPU Core数 8 6

(4FPUs) (AVX available) CPU Clock (GHz) 3.1 3.2 Memory Capacity (GB) 16 16

Linux Kernel version 3.6.11-4 3.11.10-100 CUDA Version 7.0RC 6.5 Driver Version 346.29 343.19 GNU gcc Version 4.6.3 4.7.2 Intel MKL Version 13.0.1 13.0.1

0 10 20 30 40 50 60 70 80 90 0 2000 4000 6000 8000 10000 12000 14000 GFlops N DSYMV (GeForce GTX 980)

ASPEN.K2 1.5p3x DSYMV, Upper ASPEN.K2 1.5p3x DSYMV, Lower KBLAS 1.2, DSYMV, Upper KBLAS 1.2, DSYMV, Lower MAGMA 1.6.1, DSYMV, Upper MAGMA 1.6.1, DSYMV, Lower MAGMA 1.6.1, DSYMV, Upper (work) MAGMA 1.6.1, DSYMV, Lower (work) CUBLAS 7.0RC, DSYMV, Upper (atomics mode) CUBLAS 7.0RC, DSYMV, Lower (atomics mode)

0 10 20 30 40 50 60 70 0 2000 4000 6000 8000 10000 12000 14000 GFlops N DSYMV (Tesla K20c)

ASPEN.K2 1.5p2 DSYMV, Upper ASPEN.K2 1.5p2 DSYMV, Lower KBLAS 1.2, DSYMV, Upper KBLAS 1.2, DSYMV, Lower MAGMA 1.6.1, DSYMV, Upper MAGMA 1.6.1, DSYMV, Lower MAGMA 1.6.1, DSYMV, Upper (work) MAGMA 1.6.1, DSYMV, Lower (work) CUBLAS 6.5, DSYMV, Upper (atomics mode) CUBLAS 6.5, DSYMV, Lower (atomics mode)

図1 DSYMVの各CUDA-BLAS実装における性能

0 10 20 30 40 50 60 70 0 2000 4000 6000 8000 10000 12000 14000 GFlops N DSYMV (Tesla K20c)

ASPEN.K2 1.5p2 DSYMV, Upper ASPEN.K2 1.5p2 DSYMV, Lower

0 10 20 30 40 50 60 70 0 2000 4000 6000 8000 10000 12000 14000 GFlops N DSYMV (Tesla K20c)

CUBLAS 6.5, DSYMV, Upper (atomics mode) CUBLAS 6.5, DSYMV, Lower (atomics mode)

0 10 20 30 40 50 60 70 0 2000 4000 6000 8000 10000 12000 14000 GFlops N DSYMV (Tesla K20c)

MAGMA 1.6.1, DSYMV, Upper MAGMA 1.6.1, DSYMV, Lower MAGMA 1.6.1, DSYMV, Upper (work) MAGMA 1.6.1, DSYMV, Lower (work)

0 10 20 30 40 50 60 70 0 2000 4000 6000 8000 10000 12000 14000 GFlops N DSYMV (Tesla K20c)

KBLAS 1.2, DSYMV, Upper KBLAS 1.2, DSYMV, Lower

図3 Tesla K20cでのSYMVの性能(それぞれ行列は32次元毎に 測定) 表6 Tesla K20cにおけるMAGMA+ASPEN.K2の(1)trd部分 の実測結果(秒) 1088 2122 3136 4160 MAGMA+ASPEN.K2 0.11 0.33 0.77 1.46 MAGMA only 0.15 0.49 1.05 1.94 Eigen-G 0.11 0.33 0.74 1.42 4.3 MAGMA+ASPEN.K2 最後に, MAGMA+ASPEN.K2での異なるライブラリの 組み合わせによりどの程度まで性能が向上するか前節ま での議論の検証の意味で紹介する. magma dsyevdの(1) のサブルーチン（magma dsytrd)中から呼び出されてい

るdsymvに対して, ASPEN.K2のdsymv実装を呼び出

す改良を加えて前節の条件と同様に実行した結果が以下 のようになる((2)(3)は変更がない為(1)のみを示した). MAGMA+ASPEN.K2の計算時間はEigen-G単体のもの と殆ど同じとなっている. このことからも, CUDA-BLAS の選択が性能に大きく影響していることは納得できる. MAGMAならびにEigen-G両者は同様のアルゴリズムに 基づくためこれ以上の差は実装の僅かな差によるところが 多いと推測される. 後は, MAGMAではBLAS以外の部分 もGPUにオフロードするなどの積極的な実装があること も評価しなくてはいけない. なお, 現時点でMAGMAと Eigen-Gを同時リンクし相互に利用するまでには至ってい ない. 今後の開発でこれらの問題を解決し,真の意味での 最速固有値ソルバーの実現を目指したい. 4.4 その他の議論 最 後 に, 2 ス テ ッ プ ア ル ゴ リ ズ ム を 採 用 す る

magma dsyevdx 2stageにも言及しておかなくてはなら

ない. 2ステップアルゴリズムはメモリバンド幅にシビア なハウスホルダー三重対角化の性能を改善する特効薬とし て注目を集めている. しかしながら,全ての固有値・固有 ベクトルを計算するシナリオの下ではハウスホルダー三重 対角化を経由する1ステップアルゴリズムと計算時間の 面で大きな優位性を持てない状況にある. 具体的な実験例 は示さないが, 現状のMAGMAライブラリを実行しても MAGMA+ASPEN.K2のものと大きな差はない. 環境に よっては有利な状況もありうるかもしれない. 性能向上の ボトルネックは帯_→三重対角の逆変換部分の実装にある, 今後の新しい実装方法の出現による性能改善に期待したい.

5.

まとめ

本研究では,既存の固有値計算で最速となる固有値ソル バーに対してCUDA-BLASの選択による性能の違いを検 証し,最速ソルバーの組み合わせを発見することを目的とし た議論と実験を行った. Eigen-G, MAGMA共に共通のアル ゴリズムを採用しており,本質的な性能差はCUDA-BLAS の実装部分に由来することは確認できた. 本研究で確認

される範囲では, MAGMA+ASPEN.K2の実装組み合わせ が, CPU+single GPUの実行環境においては最速の組み合 わせとなるはずである. 本報告ではMAGMAに収録され た2stageアルゴリズムの詳細を殆ど紹介するできなかった が, MAGMAチームの努力により性能向上が目覚ましい. 次回の研究会発表の際にご紹介できればと考えている. 最 後 に, 本 研 究 は 科 研 費 新 学 術 領 域 研 究(課 題 番 号: 22104003) ならびに計算科学振興財団(研究助成研究教 育拠点(COE)形成事業)の支援を受けている. 参考文献

[1] Imamura, T., ASPEN-K2: Automatic-tuning and Stabilization for the Performance of CUDA BLAS Level 2 Kernels, 15th SIAM Conference on Par-allel Processing for Scientific Computing (PP2012), http://www.siam.org/meetings/pp12/

[2] 椋木大地, KeplerアーキテクチャGPUにおける高速な SGEMVの実装, GTC Japan 2014.

[3] NVIDIA Corporation, CUDA C Programming guide, http://docs.nvidia.com/cuda/pdf/CUDA C Programm ing Guide.pdf (2014).

[4] NVIDIA Corporation, The NVIDIA CUDA Basic Linear Algebra Subroutines,

http://developer.nvidia.com/cublas

[5] Innovative Computing Laboratory, University of Ten-nessee, Matrix Algebra on GPU and Multicore Archi-tectures, http://icl.cs.utk.edu/magma

[6] 今村俊幸, 内海貴弘, 林熙龍, 山田進, 町田昌彦, Fermi, Kepler複数世代GPUに対するSYMVカーネルの性能 チューニング,情報処理学会研究報告,「ハイパフォーマ ンスコンピューティング（HPC）」, Vol. 2012-HPC-138, No. 8 (2012) 1–7. [7] 今村俊幸,内海貴弘,山田進,町田昌彦, CUDA-xSYMV の実装と評価,情報処理学会研究報告,「ハイパフォーマ ンスコンピューティング（HPC）」, Vol. 2014-HPC-146, No. 14 (2014) 1–12.

[8] Imamura, T., Yamada, S., Machida, M., Eigen-G: GPU-based eigenvalue solver for real-symmetric dense matri-ces, 10th International Conference on Parallel Processing and Applied Mathematics (PPAM2014), LNCS8384, pp. 673-682, 2014

[9] Nath, R., Tomov, S., Dong, T. T., and Dongarra, J., Optimizing Symmetric Dense Matrix-vector Multiplica-tion on GPUs, in Proceedings of 2011 InternaMultiplica-tional Con-ference for High Performance Computing, Networking, Storage and Analysis, SC’11 (2011) 6:1–6:10.

[10] Abdelfattah, A., Keyes, D., and Ltaief, H., KAUST BLAS (KBLAS),

http://cec.kaust.edu.sa/Pages/kblas.aspx [11] Abdelfattah, A., Keyes, D., and Ltaief, H., KBLAS: High

Performance Level-2 BLAS on Multi-GPU Systems, http://ondemand.gputechconf.com/gtc/2014/poster /pdf/P4168 KBLAS GPU computing optimization.pdf, GTC2014 (2014).

[12] Sørensen, H. H. B., Auto-tuning Dense Vector and Matrix-Vector Operations for Fermi GPUs, Parallel Pro-cessing and Applied Mathematics, LNCS 7203 (2012) 619–629.

[13] Sørensen, H. H. B.. Auto-Tuning of Level 1 and Level 2 BLAS for GPUs, Concurrency Computat.: Pract. Ex-per., Wiley (2012) 1183–1198.

[14] GPUlab: GLAS library version 0.0.2, http://gpulab.imm.dtu.dk/docs/

glas v0.0.2 C2050 cuda 4.0 linux.tar.gz [15] Imamura, T., Yamada, S., and Machida, M., A High

Performance SYMV Kernel on a Fermi-core GPU, High Performance Computing for Computational Science — VECPAR 2012, LNCS 7851 (2013) 59–7.

[16] Humphrey, J.R., Price, D. K., Spagnoli, D. K., Paolini, A. L., Kelmelis, E. J., CULA: Hybrid GPU Acceler-ated Linear Algebra Routines, SPIE Defense and Secu-rity Symposium (DSS), April, 2010.

[17] CULA の ベ ン チ マ ー ク に 関 し て は 以 下 を 見 よ. http://www.culatools.com/dense/ [18] EigenExa ホ ー ム ペ ー ジ: http://www.aics.riken.jp/labs/lpnctrt/EigenExa.html EigenK ホ ー ム ペ ー ジ http://ccse.jaea.go.jp/ja/download/eigenk.html