NEC

MathKeisanユーザーズガイド

Solver

概要

MathKeisanのSolverは、行列 $ A $ が対称である疎行列の連立1次方程式 $ Ax=b $ を直接法で解きます。
行列$ A $ を $ LDL^T $ 分解($ LDL^T \leftarrow PAP^T $) して、 $ x $ を前進/後退三角解法により計算します。ここで、$ L $ は下三角疎行列で、$ D $ は対角行列、$ ^T $ は転置を意味します。$ P $ は、疎因数 $ L $ のfill-inを減らすための置換行列で、Nested Dissection法またはMultiple Minimum Degree法で計算されます。分解時の数値安定性のための軸選択は行いませんが、反復改良により解の精度を向上させるサブルーチンがあります。

ユーザインタフェース

ユーザインタフェース情報は、いくつかの箇所に記載されています。

Solverの処理段階

疎行列の連立1次方程式 $ Ax=b $ を解く手順は、大きく3つの処理に分かれます。

1. 前処理

ldlt_preprocess ( token, n, colptr, rowind, Lnnz, order )

この処理において、置換行列 $ P $ が作成されます。

元の疎行列および置換行列 $ P $ により、因数 $ L $ の非ゼロ構造はシンボリック分解により決定されます。

前処理段階では、以降の計算で行われる数値分解および三角解法を容易にするため、動的メモリ割り付けも行います。

2. 数値分解

ldlt_fact ( token, colptr, rowind, nzvals, flops )

下三角行列 $ L $ および対角行列 $ D $ が、$ PAP^T = LDL^T $ となるように作成されます。本処理が、最も計算が集中する処理です。

3. 三角解法

ldlt_solve ( token, x, rhs, flops )

行列が分解され、右側の $ b $ が与えられると、連立1次方程式は以下のように解法されます。

求解 $ Ly=Pb $ を $ y $ について解く
代入 $ y \leftarrow D^{-1}y $
求解 $ Lx=y $ を $ x $ について解く
代入 $ x \leftarrow P^{-T}x $

補助サブルーチン

他に3つのルーチンがあります。

ldlt_getdiag ( token, diag )
ldlt_refine (token, colptr, rowind, nzvals, rhs, x, newx, flops )
ldlt_free ( token )
これら3つのルーチンは、それぞれ"対角行列 $ D $ の抽出"、"解の改良"、"連立1次方程式の求解に関連する全ての作業領域の解放"を行うために提供しています。

行列オーダリング

MathKeisan Solverは、 Multiple Minimum Degree(MMD) [4]または Nested Dissection(ND) [3]行列オーダリングを用います。この選択は、ldlt_preprocessサブルーチンに適切な整数引数orderを渡すことにより行います。

order = 0 の時、修正MMDルーチンが使われます。本ルーチンは、オーダリング中にsupernode情報を抽出します。supernodeとは、同じ疎構造を共有する因数Lの列の集まりです。物理系の行列では、一般的にsupernodeの多い因数 $ L $ ができます。supernode情報により分解を効率的に行うことができます。

order = 1 の時、オリジナルのMMDアルゴリズムが使われます。

order = 2またはorder = 3の時、ノードベースまたはエッジベースのNDアルゴリズムがそれぞれ使われます。一般的に、ノードベースのNDアルゴリズム(order = 2)の方が、よりよいオーダリングが行えます。NDオーダリングの計算時間は一般的にMMDオーダリングの計算時間より長くなりますが、NDオーダリングの方が一般的にfill-inが小さくなり、かつorder + factorization + solveの合計時間がMMD使用時よりも短くなります。

order = 1、2または3の時、オーダリング中にsupernode情報は収集されません。supernodeを特定し、行列オーダリングを行った後の $ L $ に対する格納要件を推定するのにLiu、Ng および Peyton [6]が開発したルーチン群が使われます。

分解アルゴリズム

データのアクセスおよび変更パターンにより、$ LDL^T $分解はいくつかの異なる方法で行うことができます。最も広く使われているのは、left-looking [7]、right-looking [8]、multifrontal [5]です。これらの実装の比較は、Rothberg [8]で報告されています。MathKeisanでは、left-lookingアルゴリズムが使われています。

LAPACKなどのDense Linear Algebraライブラリでは、性能最適化のためアルゴリズムはブロック化されます。このため、レベル2および3のBLASへの呼び出し時に効率的にデータの流用ができます。MathKeisanのSolverはsupernodeを使用しています。これにより、疎行列に対しブロック化されたアルゴリズムを密行列の場合と同じデータ流用率で使用でき、レベル2および3のBLASへの呼び出しが行えます。supernodeの使用により、間接アドレス指定も減少します。

並列分解

通常、疎行列の連立1次方程式の直接解法において、数値分解は最も時間のかかる処理です。MathKeisanのSolverでは、数値分解をOpenMPによって並列化することができます。

環境変数OMP_NUM_THREADSがnpに設定されている場合、分解は np 個のスレッド上で実行されます。OMP_NUM_THREADS が設定されていない場合には、分解は mp 個のスレッド上で実行されます。この時、mp はリソースグループまたはノード内の最大プロセス数です。

反復改良

分解時、数値安定化のための軸選択は行われません。行列が適切に条件づけられており、かつ対称正定値(SPD)行列である場合、解は安定します。行列が非正定値または適切に条件づけられていない場合、反復改良により解の精度を向上させることができます。

解$ x $が、$ \hat{x} $ が真の解、$ \delta x $ がエラーである$ x=\hat{x}+\delta x $により与えられると、$ A(\hat{x}+\delta x)=b+\delta b $となり、この時$ \delta b $ が$ \delta b=Ax-b $ により与えられる残差です。エラー$ \delta x $ は、$ A \delta x = \delta b $を解くことにより計算できます。

反復改良サブルーチンは、3つのステップで実行されます。

  1. 残差 $ \delta b=Ax-b $ を計算する。他のいくつかの実装では、この計算は拡張精度で行われています。MathKeisanでは、拡張精度を使いません。
  2. $ A \delta x = \delta b $ を解くことによりエラー $ \delta x $ を計算する。
  3. 改良した解 $ x \leftarrow x + \delta x $ を計算する。

Solverルーチン一覧

名前 説明
ldlt_preprocess 行列を並べ替え、ldlt_fact にメモリを割り当てる
ldlt_fact ldlt_preprocess からの並び替え行列を使って $ LDL^T $ 分解を計算
ldlt_solve ldlt_fact からの分解行列を使って連立1次方程式を解く
ldlt_refine 反復改良
ldlt_getdiag ldlt_fact からの分解行列の対角を取得
ldlt_free ldlt_preprocess に割り当てられたメモリを解放

参考文献

  1. J. W. Demmel. Applied Numerical Linear Algebra. SIAM, Philadelphia, 1997.
  2. N. J. Higham. Accuracy and Stability of Numerical Algorithms. SIAM, Philadelphia, 1996.
  3. G. Karypis and V. Kumar. "A fast and high quality multilevel scheme for partitioning irregular graphs." SIAM Journal on Scientific Computing, 20(1):359-392, 1999.
  4. J. W. Liu. "Modification of the minimum-degree algorithm by multiple elimination." ACM Transactions on Mathematical Software, 11(2):141-153, 1985.
  5. J. W. Liu. "The multifrontal method for sparse matrix solution: Theory and practice." SIAM Review, 34(1):82-109, 1992.
  6. J. W. Liu, E. Ng, and B. Peyton. "On finding supernodes for sparse matrix computations." SIAM J. Matrix Ahal. Apply., 14:242-252, 1993.
  7. E. G. Ng and B. W. Peyton. "Block sparse Cholesky algorithms on advanced uniprocessor computers." SIAM J. SCI. COMPUT., 14:1034-1056, 1993.
  8. E. Rothberg and A. Gupta. "An evaluation of left-looking, right-looking and multifrontal approaches to sparse Cholesky factorization on hierarchical-memory machines." Int. J. High Speed Computing, 5:537--593, 1993.