疎行列固有値解法における4倍精度演算とその性能評価

全文

(1)情報処理学会研究報告 IPSJ SIG Technical Report. Vol.2010-ARC-192 No.18 Vol.2010-HPC-128 No.18 2010/12/17. 規 模並列環境 へ の 対応 を 中心に 研 究 を 行う と と も に ,同年 11 月よ り 疎行列固有値解法に 対 応 し た 新版を 公開し て いる .. 疎行列固有値解法に お け る 4 倍精度演 算と その 性能評価. Lis は , Krylov 部分空 間 法を 中心と す る 多様な 反復解法を 実装し た ライ ブ ラリで あ る .同 様な 目的の ライ ブ ラリと し て ,線型方程式に つ いて は Argonne 米国 立研 究 所の 並列反復解. 西. 田. 晃†1. 法ライ ブ ラリ PETSc ラリ Hypre. 計算科 リ性能と 研 究 で は 度演 算を を 得た が す る .. よ る SLEPc. 学に お いて ,計算精度の 向上は 本質的に 重要な 目標で あ る .演 算性能は メ モ 比較し て 速い速度で 向上し て お り ,こ の ギ ャ ップ が 問題と な る こ と が あ る .本 並列反復解法ライ ブ ラリ Lis を 対象に ,double-double 精度に よ る 4 倍精 実装し ,こ の 問題を 計算精度を 向上させ る こ と に よ っ て 緩和で き る と の 知見 ,こ こ で は こ れ を 固有値解法に 適用し ,4 倍精度演 算の 有効性に つ いて 検 討. ?3 ?4. ?2. や Lawrence Berkeley 米国 立研 究 所に よ る 並列反復解法ライ ブ. な ど を 挙 げる こ と が で き る .固有値解法に つ いて は ,Valencia 工科 大学に. (PETSc を 用いて 開発され て いる ) や , Colorado 大学に よ る BLOPEX. (Hypre を 用いて いる ) な ど に 疎行列を 対象と し た 固有値解法が 実装され て いる. ?6. ?5. .. Lis で は こ れ ら の 機能を 単一 の ライ ブ ラリに お いて 実現す る と と も に ,多様な 前処理ア ル ゴ リズ ム を 実装し て いる .ま た ,近 年一 般的と な っ た マ ルチ コ ア 環境 に 適し た ハ イ ブ リッド 並列処理に 対応 し て いる 点も 特徴と し て 挙 げる こ と が で き る .表 1-3 に 現時点で 対応 し て いる 固有値解法, 線型方程式解法, 行列格納形式の 一 覧を 示す .. Quadruple Precision Operation in Eigensolvers. 表 1 Lis で 利用可 能な 固有値解法 Power Iteration Inverse Iteration Approximate Inverse Iteration Subspace Iteration Lanczos Iteration Conjugate Gradient Conjugate Residual. for Sparse Matrices and its Performance Evaluation Akira Nishida. †1. The improvement of the accuracy is a critically important target of computational science, FLOPS performance is increasing with faster speed than memory performance. This gap sometimes raises serious problem, but our previous study have shown that the problem can be reduced by improving floating point precision. In this paper, we discuss the validity of the quadruple precision arithmetics for eigensolvers, based on the double-double precision floating point operations implemented on the parallel iterative solver library Lis.. 線型方程式の 求 解に 用いら れ る 共 役勾配法等の 反復解法は ,理論的に は た か だ か n 回の 反復で 収束す る こ と が 知ら れ て いる .し か し な が ら ,実際に は 丸 め 誤差の 影響 に よ り ,有限 精度計算で は 一 般に 収束ま で に は よ り 多く の 反復回数を 要し ,ま た 収束が 停滞す る 場合も あ る .こ の よ う な 収束の 特性を 改善す る 上で ,多倍長演 算は 有効な 手法で あ る と 考えら れ る が ,. 1. 背. ハ ー ド ウ ェ ア で 実装され て いる 倍精度演 算等に 比べ ,計算コ ス ト の 大き さが 問題と な る .. 景. こ の 問題を 解決す る た め ,Lis で は ,近 年の プ ロセ ッサ に 搭載され て いる SIMD 命令を. 本研 究 で は , 平成 14-19 年度科 学技 術振興機構 CREST 事業 の 一 環と し て , 反復解法ライ ブ ラリ Lis. ?1. http://www-unix.mcs.anl.gov/petsc/ http://computation.llnl.gov/casc/hypre/ http://www.grycap.upv.es/slepc/ http://www-math.cudenver.edu/~aknyazev /software/BLOPEX/ ?6 http://www.netlib.org/utk/people /JackDongarra/la-sw.html. を 開発, 配布し , 様々 な 並列計算機上で 大規 模な 線型方程式を 解く た め の 環境. ?2 ?3 ?4 ?5. を 提供 し て き た . ま た 平成 20 年度か ら は 九 州大学情報基盤研 究 開発セ ンタ ー に お いて ,大 †1 九 州大学情報基盤研 究 開発セ ンタ ー Research Institute for Information Technology, Kyushu University ?1 http://www.ssisc.org/lis/. 1. ⓒ2010 Information Processing Society of Japan.

(2) 情報処理学会研究報告 IPSJ SIG Technical Report. Vol.2010-ARC-192 No.18 Vol.2010-HPC-128 No.18 2010/12/17. 表2. in the last place” を 意 味す る . 12 ulp は 丸 め 誤差の 上限で あ る .. Lis で 利用可 能な 線型方程式解法 CG CR BiCG BiCR CGS CRS BiCGSTAB BiCRSTAB GPBiCG GPBiCR BiCGSafe BiCRSafe BiCGSTAB(l) BiCRSTAB(l) Jacobi Gauss-Seidel SOR Orthomin(m) TFQMR MINRES GMRES(m) FGMRES(m) IDR(s). す べ て の 演 算は ,IEEE 倍精度演 算で round-to-even 丸 め と 仮 定す る .x と y を 倍精度と し ,x+y の 倍精度加 算の 結果 を fl(x+y)  と 表す .err(x+y) は x+y = fl(x+y)+err(x+y) を 満た す も の と す る .こ の 時,以下 の よ う に し て 4 倍精度加 算 a = b + c が 計算で き る .た だ し a = (a.hi, a.lo), b = (b.hi, b.lo), c = (c.hi, c.lo) と す る . ま ず b と c の 上位 b.hi と c.hi に ,丸 め 誤差の な い加 算. b.hi + c.hi = fl(b.hi + c.hi) + err(b.hi + c.hi). (1). を 行い,. fl(b.hi + c.hi) + fl(err(b.hi + c.hi) + b.lo + c.lo) 表 3 Lis で 利用可 能な 行列格納形式 Compressed Row Storage Compressed Column Storage Modified Compressed Sparse Row Diagonal Ellpack-Itpack generalized diagonal Jagged Diagonal Block Sparse Row Block Sparse Column Variable Block Row Dense Coordinate. (2). を a + b の 近 似値と す る .今回の 実装で は 高速性を 重視し て 下 位の 誤差 fl(err(b.hi + c.hi) +. (CRS) (CCS) (MSR) (DIA) (ELL) (JDS) (BSR) (BSC) (VBR) (DNS) (COO). b.lo + c.lo) は 無視し て いる .な お ,下 位の 足し 合わ せ に よ っ て 繰り 上が り の 可 能性が あ る た め ,こ の 計算は 丸 め 誤差の な い加 算に よ っ て 行う .. 4 倍精度乗算に つ いて も x × y = fl(x × y) + err(x × y) と す る .こ こ で. 1 ulp(fl(x × y)) 2. ≥. |err(x × y)| で あ る .4 倍精度乗算は ,ま ず b と c の 上位 b.hi と c.hi に 丸 め 誤差の な い乗 算を 行い,. b.hi + c.hi = fl(b.hi × c.hi) + err(b.hi × c.hi). (3). と す る .次に ,b の 上位と c の 下 位,b の 下 位と c の 上位の 乗算結果 と err(b.hi × c.hi) と の 間 で 丸 め 誤差の な い加 算を 行い,b × c の 近 似値を 得る .. 活用す る こ と に よ り ,4 倍精度演 算を 効率よ く 実行で き る よ う 実装を 行っ て いる . 反復解法 の 丸 め 誤差を 減少させ る た め ,線型方程式に お け る 係数行列と 右辺ベ ク ト ル,初期ベ ク ト. 4 倍精度ベ ク ト ルの デ ー タ 構造と し て は ,上位と 下 位を それ ぞ れ 別な 配列に 格納す る 方法,. ルに 倍精度を 使用し ,内部処理の み を 4 倍精度化 す る こ と を 考える .Lis で は ,こ の 目的の. 交互に 格納す る 方法が 考えら れ る が ,前者は 上位を 格納す る 配列を 用いて 倍精度演 算を 行う. た め に Bailey ら に よ っ て 提案され た 倍精度浮動小数点数を 2 個使用す る “double-double”. こ と が で き る た め ,こ こ で は 前者の 方法を 採用す る .. 精度ア ルゴ リズ ム. 1). 反復解法の 計算の 主要部は ,疎行列-ベ ク ト ル間 演 算,内積,及 び ベ ク ト ル間 演 算か ら な. を 用い,Intel プ ロセ ッサ に 実装され て いる SIMD 命令の 一 種で あ る. SSE2 に よ り ,こ れ を 実装し て いる .“double-double” 精度の 実装に つ いて は Hida ら. 2). る .こ れ ら を 4 倍精度演 算に 置き 換える .た だ し 倍精度演 算と 同一 の イ ンタ フ ェ ー ス を 保つ. に. よ る QD ライ ブ ラリな ど が あ り ,こ れ を 用いた 反復解法ライ ブ ラリ GMM++3) が 存在す. た め ,. る .本研 究 は ,倍精度演 算に よ っ て 実装され る 4 倍精度演 算を , SIMD 命令に よ っ て 高速. • 係数行列,右辺ベ ク ト ルは 倍精度. 化 し た 点で こ れ ら の 研 究 と 異な っ て いる .以下 で は ,具 体的な 実装手法に つ いて 説明す る .. • 解ベ ク ト ルの 入出力は 倍精度,内部演 算は 4 倍精度 • 反復解法中の ベ ク ト ルと ス カ ラー は 4 倍精度. 2. 4 倍精度演 算の 実装. と す る .こ の た め ,疎行列-ベ ク ト ル間 演 算で は 倍精度-4 倍精度間 演 算,内積,ベ ク ト ル間 演 算,ス カ ラ間 演 算で は 4 倍精度-4 倍精度間 演 算を 実装す る 必要が あ る .. Bailey の ア ルゴ リズ ム で は ,double-double 精度浮動小数 a を a = a.hi + a.lo, 1 ulp(a.hi) 2. こ れ ら の 演 算で は 実行時間 が 問題と な る た め ,SIMD 命令の 利用に よ る 高速化 を 検 討し た .. ≥ |a.lo| (上位の a.hi, 下 位の a.lo と も 倍精度) と し て ,4 倍精度演 算を 倍精. 度の 四則演 算の 組み 合わ せ で 実現す る .な お ulp(x) は x の 仮 数部の 誤差範囲 を 示す “unit. Intel プ ロセ ッサ の 場合を 考える と ,SIMD 命令と し て ,SSE (Streaming SIMD Extensions). 2. ⓒ2010 Information Processing Society of Japan.

(3) 情報処理学会研究報告 IPSJ SIG Technical Report. Vol.2010-ARC-192 No.18 Vol.2010-HPC-128 No.18 2010/12/17. が 搭載され て いる .こ こ で は ,128 ビ ット デ ー タ に 対す る 処理が 可 能な SSE2 を 用いて ,倍. 0=4.46e-2. 1=1.11e-1. 2=1.11e-1. 1=2.20e-1. 5=2.20e-1. 精度浮動小数に 対し て 同時に 2 つ の 演 算を 行う こ と を 考える .計算の 依 存関 係に よ る 性能 の 低下 を 考慮し ,2 段の ルー プ ア ンロー リング を 併用し た .ま た 各演 算の ルー プ 内で 使用す る ス カ ラ変数に つ いて は ,必要に 応 じ て 128 ビ ット XMM レジス タ に 格納す る と と も に ,. 16 バ イ ト の ア ライ ンメ ント 境 界を 考慮し た デ ー タ 処理を 行っ て いる . 3=1.78e-1. 3. 固有値解法へ の 適用 疎行列固有値解法の 適用例と し て ,2 次元 Helmholtz 問題を 考える .矩 形領域 [0, l]×[0, m] の 膜の 振動を 表す 2 次元 Laplace 作用素を 5 点中心差分に よ り 離散化 し た 場合,対応 す る 行列の 固有値は. π2. . σ2 τ2 + 2 2 l m. . 図1. , σ, τ ∈ N. % cumulative time seconds 66.67 0.04 16.67 0.05 16.67 0.06. で 与えら れ る .l = m = 20,す な わ ち 行列サ イ ズ 202 × 202 ,と し て ,Lis に 実装し た 部分 空 間 反復法に よ り 絶対値最小の も の か ら 順に 6 個の 固有対を 求 め た 結果 を 図 1 に 示す .な お 内部の 解法に は 逆 反復法,前処理な し 共 役勾配法を 使用し て いる . Intel Xeon 5570 サ ー バ (2.93GHz ク ア ッド コ ア プ ロセ ッサ ×2) の 1 コ ア 上で 倍精度,4 倍精度で の 逐次計算を. 図2. 行っ た 場合の 計算時間 の 内訳の う ち ,1%を 越 える も の を 図 2-3,同様の 計算を 行列サ イ ズ を. 2002 × 2002 と し て 行っ た 場合の 内訳を 図 4-5 に 示す .行列サ イ ズ を 2002 × 2002 と し た 場 の 数を 1 と し た 場合に は ,倍精度,4 倍精度で の 計算時間 は それ ぞ れ 1.75 秒,1.74 秒と な っ た .ま た ,内訳か ら CRS 形式で の 疎行列ベ ク ト ル積を 計算す る ルー チ ン lis matvec crs() の 計算時間 の 減少分が 大き いこ と が 分か る が ,こ の こ と か ら ,4 倍精度演 算は 直交化 に お け る 反復回数の 減少に 関 連が あ る こ と が 見て 取れ る . 各固有値の 計算に 要す る 反復回数を ,表 5 に 示す .モ ー ド 3 に つ いて は ,残差の 閾値 10−12 −12. self seconds 0.04 0.01 0.01. calls 10374 12154 1768. self s/call 0.00 0.00 0.01. total s/call 0.00 0.00 0.01. name lis_matvec_crs lis_vector_copy lis_vector_set_all. 部分空 間 反復法に お け る 各処理の 割合( Intel Xeon 5570 サ ー バ 上,問題サ イ ズ 202 × 202 ,倍精度. Lis の プ ロフ ァ イ リング オ プ シ ョ ンに よ る . ). % cumulative time seconds 55.56 0.05 11.11 0.06 11.11 0.07 11.11 0.08 11.11 0.09. 合,倍精度,4 倍精度で の 計算時間 は それ ぞ れ 63.32 秒,44.64 秒で あ っ た .求 め る 固有値. に 達し て いな いた め ,反復回数の 上限 1000 回 (残差 1.04 × 10. 2 次元 Helmholtz 問題( 膜の 振動) に お け る モ ー ド. (4). self seconds 0.05 0.01 0.01 0.01 0.01. calls 9398 17252 4650 2849 772. self ms/call 0.01 0.00 0.00 0.00 0.01. total ms/call 0.01 0.00 0.00 0.00 0.01. name lis_matvec_crs lis_vector_dot lis_vector_duplicateex lis_vector_axpyex_mmm lis_cg_check_params. 図 3 部分空 間 反復法に お け る 各処理の 割合( Intel Xeon 5570 サ ー バ 上,問題サ イ ズ 202 × 202 ,4 倍精度. ). ) で 打ち 切っ て いる .よ っ. て ,こ の 固有値の 計算に お け る 残差の 停滞が ,全体の 計算時間 に 影響 を 与えて いる こ と が 分. ら に ,MPI 版を 使用し て 使用コ ア 数を 8 に 増や し ,30 個ま で の 固有値計算を 行っ た 結果 を. か る .. ??に 示す .倍精度で の 計算時間 は 123 秒,4 倍精度で の 計算時間 は 489 秒で あ っ た .モ ー. 一 方,計算す る 固有値の 数を 20 個に 増や し た 場合の 結果 を 表?? に 示す .モ ー ド の 小さい. ド が 大き く な る と 収束が 停滞す る 傾向は よ り 顕著に な っ て いる が ,こ の 場合に は 4 倍精度. 固有値に つ いて は 倍精度で の 計算と 比較し て 少な い反復回数で 収束し て いる が ,モ ー ド が 大. 演 算で も 十分な 効果 が 得ら れ て いな いこ と が 分か る .. き く な る と 4 倍精度で あ っ て も 収束が 停滞す る も の が 出て き て いる .な お ,こ の 場合の 倍精. 以上の 結果 か ら ,計算す る 固有値の 量が 増える に つ れ ,要求 され る 精度が 高く な り ,4 倍. 度で の 計算時間 は 298 秒で あ っ た の に 対し ,4 倍精度で の 計算時間 は 396 秒で あ っ た .さ. 精度演 算で も 十分な 精度が 得ら れ な く な っ て いる と 考えら れ る .4 倍精度演 算が 有効な の は ,. 3. ⓒ2010 Information Processing Society of Japan.

(4) 情報処理学会研究報告 IPSJ SIG Technical Report. % cumulative time seconds 64.62 40.92 12.00 48.52 9.16 54.32 5.50 57.80 4.34 60.55 2.91 62.39 1.06 63.06. self seconds 40.92 7.60 5.80 3.48 2.75 1.84 0.67. Vol.2010-ARC-192 No.18 Vol.2010-HPC-128 No.18 2010/12/17. calls 124700 251966 253778 124700. self s/call 0.00 0.00 0.00 0.00. total s/call 0.00 0.00 0.00 0.00. 130137. 0.00. 0.00. 表5. name lis_matvec_crs lis_vector_axpy lis_vector_dot lis_vector_xpay __intel_new_memcpy lis_vector_nrm2 __intel_new_memset. 問題サ イ ズ 2002 × 2002 の 場合に 各固有値の 計算に 要す る 反復回数.計算す る 固有値の 数を 20 個と し た 場合.. Mode 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19. 図 4 部分空 間 反復法に お け る 各処理の 割合( Intel Xeon 5570 サ ー バ 上,問題サ イ ズ 2002 × 2002 ,倍精度. ). % cumulative time seconds 62.19 27.76 11.20 32.76 8.92 36.74 4.86 38.91 3.67 40.55 3.16 41.96 2.51 43.08 1.90 43.93. self seconds 27.76 5.00 3.98 2.17 1.64 1.41 1.12 0.85. calls 84732 167724 167724 84732. self s/call 0.00 0.00 0.00 0.00. total s/call 0.00 0.00 0.00 0.00. 87343 4203 3333. 0.00 0.00 0.00. 0.00 0.00 0.00. name lis_matvec_crs lis_vector_axpy lis_vector_dot lis_vector_xpay __intel_new_memcpy lis_vector_nrm2 lis_vector_dotex_mmm lis_vector_axpyex_mmm. 図 5 部分空 間 反復法に お け る 各処理の 割合( Intel Xeon 5570 サ ー バ 上,問題サ イ ズ 2002 × 2002 ,4 倍精度. ). 少数の 固有値を 求 め る 場合に 限ら れ る 可 能性が 高いが ,適用範囲 に つ いて は 検 討が 必要で. 4. ま. あ る .. Iteration (double precision) 18 78 122 1000 102 492. Iteration (quad precision) 18 60 117 120 96 459 232 118 604 234 253 401 312 345 1000 1000 253 1000 1000 1000. め. 本稿で は ,並列反復解法ライ ブ ラリ Lis を 対象に ,double-double 精度に よ る 4 倍精度演. 表 4 問題サ イ ズ 2002 × 2002 の 場合に 各固有値の 計算に 要す る 反復回数.. Mode 0 1 2 3 4 5. と. Iteration (double precision) 18 78 122 1000 102 492 252 114 667 1000 249 439 272 299 995 300 263 1000 1000 1000. 算を 固有値解法に 適用し ,4 倍精度演 算の 有効性に つ いて 検 討し た .直交化 を 行いな が ら 複. Iteration (quad precision) 18 60 117 120 96 459. 数の 固有値を 求 め る 必要の あ る 問題の 場合,数値実験で 示し た よ う に ,計算を 進め る 過程 で 特定の 固有値に 関 し て 残差の 収束が 停滞す る こ と が あ り ,その よ う な 場合の 対策と し て , 本手法は 有効で あ る と 考えら れ る .た だ ,計算す る 固有値が 多い場合に は ,内部で の 4 倍精 度演 算に は 限界が あ る 可 能性が あ り ,適用範囲 に つ いて は よ り 詳細に 検 討す る 必要が あ る .. 参 考. 文. 献. 1) Bailey, D. H.: QD: C++/Fortran-90 double-double and quad-double package, http://crd.lbl.gov/dhbailey/mpdist/ (2010).. 4. ⓒ2010 Information Processing Society of Japan.

(5) 情報処理学会研究報告 IPSJ SIG Technical Report. Vol.2010-ARC-192 No.18 Vol.2010-HPC-128 No.18 2010/12/17. 表 6 問題サ イ ズ 2002 × 2002 の 場合に 各固有値の 計算に 要す る 反復回数.計算す る 固有値の 数を 30. 2) Hida, Y., Li, X.S. and Bailey, D.H.: Algorithms for quad-double precision floating point arithmetic, Proceedings of 15th Symposium on Computer Arithmetic, pp. 155–162 (2001). 3) Renard, Y. and Pommier, J.: GMM++ User Guide, http://home.gna.org/getfem/gmm int (2010).. 個と し ,MPI 版で 8 並列で の 計算を 行っ た 場合.. Mode 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29. Iteration (double precision) 18 82 112 1000 103 449 242 120 598 241 243 409 316 330 789 285 257 1000 1000 1000 1000 440 451 1000 1000 449 1000 1000 1000 917. Iteration (quad precision) 18 64 122 99 97 451 245 118 660 225 243 471 331 322 917 293 247 1000 1000 1000 1000 471 441 1000 1000 445 1000 1000 1000 1000. 5. ⓒ2010 Information Processing Society of Japan.

(6)

表 2 Lis で 利用可 能な 線型方程式解法 CG CR BiCG BiCR CGS CRS BiCGSTAB BiCRSTAB GPBiCG GPBiCR BiCGSafe BiCRSafe BiCGSTAB(l) BiCRSTAB(l) Jacobi Gauss-Seidel SOR Orthomin(m) TFQMR MINRES GMRES(m) FGMRES(m) IDR(s) 表 3 Lis で 利用可 能な 行列格納形式 Compressed Row Storage (CRS) Compr

表 2

Lis で 利用可 能な 線型方程式解法 CG CR BiCG BiCR CGS CRS BiCGSTAB BiCRSTAB GPBiCG GPBiCR BiCGSafe BiCRSafe BiCGSTAB(l) BiCRSTAB(l) Jacobi Gauss-Seidel SOR Orthomin(m) TFQMR MINRES GMRES(m) FGMRES(m) IDR(s) 表 3 Lis で 利用可 能な 行列格納形式 Compressed Row Storage (CRS) Compr p.2
図 2 部分空 間 反復法に お け る 各処理の 割合( Intel Xeon 5570 サ ー バ 上,問題サ イ ズ 20 2 × 20 2 ,倍精度.

図 2

部分空 間 反復法に お け る 各処理の 割合( Intel Xeon 5570 サ ー バ 上,問題サ イ ズ 20 2 × 20 2 ,倍精度. p.3
表 6 問題サ イ ズ 200 2 × 200 2 の 場合に 各固有値の 計算に 要す る 反復回数.計算す る 固有値の 数を 30 個と し , MPI 版で 8 並列で の 計算を 行っ た 場合.

表 6

問題サ イ ズ 200 2 × 200 2 の 場合に 各固有値の 計算に 要す る 反復回数.計算す る 固有値の 数を 30 個と し , MPI 版で 8 並列で の 計算を 行っ た 場合. p.5

参照

Updating...

関連した話題 :

Scan and read on 1LIB APP