複素対称 (密) 行列対角化ルーチン並列化

2016 年 2 月 29 日
理化学研究所、仁科加速器研究センター
船木 靖郎 (funaki@riken.jp)• 依頼内容N × N 正定値実対称行列 B 及び N × N 複素対称行列 A(t A = A) に対し、Ax = λBx、あるいは Ax = λx の形の固有値問題を解く MPI 版のサブルーチンを開発すること。
• 支援、開発内容
1. 複素対称 (密) 行列に対し、Ax = λx の固有値問題を解き、全ての固有値、固有ベクトルを求める並列版サブルーチンの開発。
対角化のドライバルーチンは
PZDIAG(MAJOR, IP, N, A, NP, NB, A2, V2, INFO)
本ルーチンの引数仕様は Table 1 を参照。
Table 1: PZDIAG 引数仕様
引数MAJOR型
Character*1

種別INPUTIP
NA

Integer
IntegerComplex*16

INPUT
INPUT
INPUT

NP
NB

Integer
Integer

INPUT
INPUT

A2V2
INFO

Complex*16
Complex*16Integer

OUTPUTOUTPUTOUTPUT

引数の説明’R’ か’C’ を選択。
係数行列をブロック分割したプロセスマップの方向。
’R’(’C’) 行 (列) 方向に連続した MPI プロセス番号を割り当てる。
部分行列の分散方法は ScaLAPACK 準拠。
各プロセスのランク。
全体行列 A の次元数。N ≥ 2
A(N (N + 1)/2): 行列 A の要素を 1 次元配列で指定。
A(I, J) :⇒ A(I(I − 1)/2 + J) (I ≥ J),
A(J(J − 1)/2 + I) (J ≥ I)。
並列化数。 2 の偶数べき乗に設定すること (1, 4, 16, 64, 256, · · ·)。
係数行列をブロック分割した小行列の次元数。
8 以上の偶数に設定すること。
A2(N ): 1 次元配列。A の固有値。
V 2(N, N ): 2 次元配列。A の固有ベクトル。
リターンコード。

10-2-66-13
100

Table 2: リターンコード (INFO)
正常終了
N ≥ 2 が設定されていない
N B が 8 以上の整数でない
MAJOR=’R’ または MAJOR=’C’ が設定されていない
プロセスマップが正方かつ 1 辺が 2 のべき乗でない

PZDIAG における計算手順 (詳細はパッケージ内 Document/に置かれたマニュアル等参照)(1). 複素対称 t A = A(非エルミート) 行列 A を各プロセスに分配 (MATDIV)。
(2). A をハウスホルダ変換により 3 重対角化する (PZSYTRD)。
(3). ハウスホルダ逆変換により初期固有ベクトルを生成 (PZSYGHR)。
(4). (3). で計算した初期値をランク 0 のプロセスに集める (MARGEV)。
(5). ランク 0 のプロセスで QR 法により固有値、固有ベクトルを逐次計算 (QR)。
PZDIAG を使用したテストプラグラム例はパッケージ (complex diag.zip) 内の test.f。テス
トプラグラム中で呼び出している MATGEN は Input の複素対称行列 A を生成するサンプル∑プログラム。RESID は固有ベクトルの規格直交性、及び残差二乗和 i kAx(i) −λ(i) x(i) k/|λ(i) |をチェックするプログラム。
PZDIAG が使用するソースプログラムは SRC/に置かれている。PZDIAG を使用するためのライブラリ (libPZDIAG.a) を作成する makefile の例は make lib に書かれている。またパッケージ内の makefile には Fortran コンパイラに Intel コンパイラを利用する mpif90 を指定している。

2

2. 複素対称 (密) 行列に対し、Ax = λBx の一般化固有値問題を解き、全ての固有値、固有ベク
トルを求める並列版サブルーチンの開発。
対角化のドライバルーチンはPZDIAG GEN(MAJOR, IP, N, A, B, NP, NB, A2, V2, INCX, INFO)
本ルーチンの引数仕様は Table 3 を参照。本ルーチン内で、LAPACK、ScaLAPACK のルーチンを呼び出しているため、それらを使用できる環境が必須。

Table 3: PZDIAG GEN 引数仕様
引数MAJOR型
Character*1

種別
INPUT

IP
N
A

Integer
IntegerComplex*16

INPUT
INPUT
INPUT

B

Real*8

INPUT

NP
NB

Integer
Integer

INPUT
INPUT

A2
V2
INCX
INFO

Complex*16
Complex*16
Integer
Integer

OUTPUTOUTPUTINPUTOUTPUT

引数の説明’R’ か’C’ を選択。
係数行列をブロック分割したプロセスマップの方向。
’R’(’C’) 行 (列) 方向に連続した MPI プロセス番号を割り当てる。
部分行列の分散方法は ScaLAPACK 準拠。
各プロセスのランク。
全体行列 A の次元数。N ≥ 2
A(N (N + 1)/2): 行列 A の要素を 1 次元配列で指定。
A(I, J) :⇒ A(I(I − 1)/2 + J) (I ≥ J),
A(J(J − 1)/2 + I) (J ≥ I)。
B(N (N + 1)/2): 行列 B の要素を 1 次元配列で指定。
B(I, J) :⇒ A(I(I − 1)/2 + J) (I ≥ J),
B(J(J − 1)/2 + I) (J ≥ I)。
並列化数。 2 の偶数べき乗に設定すること (1, 4, 16, 64, 256, · · ·)。
係数行列をブロック分割した小行列の次元数。
8 以上の偶数に設定すること。
A2(N ): 1 次元配列。A の固有値。
V 2(N, IN CX): 2 次元配列。A の固有ベクトル。
求める固有ベクトルの個数
リターンコード。

30-2-66-13100

Table 4: リターンコード (INFO)
正常終了
N ≥ 2 が設定されていない
N B が 8 以上の整数でない
MAJOR=’R’ または MAJOR=’C’ が設定されていない
プロセスマップが正方かつ 1 辺が 2 のべき乗でない

PZDIAG GEN における計算手順 (詳細はパッケージ内 Document/に置かれたマニュアル等参照)(1). 正定値実対称行列 B 、複素対称 t A = A(非エルミート) 行列 A を各プロセスに分配 (それぞれ、MATDIV R、MATDIV)。
(2). B をコレスキー分解 B = LLT (ScaLAPACK ルーチン: PDPOTRF)。
(3). 下三角行列 L の逆行列 L−1 を求める (ScaLAPACK ルーチン: PDTRTRI)。
(4). A を標準化変換 C = L−1 AL−T (ScaLAPACK ルーチン: PZTRMM)。
(注: Ax = λBx ⇔ CLT x = λLT x)(5). C をハウスホルダ変換により 3 重対角化する (PZSYTRD)。
(6). ハウスホルダ逆変換により初期固有ベクトルを生成 (PZSYGHR)。
(7). (6). で計算した固有ベクトルの初期値をランク 0 のプロセスに集める (MARGEV)。
(8). ランク 0 のプロセスで QR 法により固有値及び (C の) 固有ベクトル V を逐次計算 (QR)。
(9). (2). で計算した下三角行列 L をランク 0 のプロセスに集める (MARGEL)。
(10). x = L−T V により A の固有ベクトルを求める (絶対値の小さな固有値に属するものからINCX 個)(LAPACK ルーチン: ZTRMV)。
PZDIAG GEN を使用したテストプラグラム例はパッケージ (complex diag gen.zip) 内のtest gen.f。テストプラグラム中で呼び出している MATGEN NORM、MATGEN HAMILはそれぞれ Input の正定値実対称行列 B 、複素対称行列 A を生成するサンプルプログラム。
∑RESID G は固有ベクトルの規格直交性、及び残差二乗和 i kAx(i) −λ(i) Bx(i) k/|λ(i) | をチェックするプログラム。
PZDIAG GEN が使用するソースプログラムは SRC/に置かれている。PZDIAG GEN を使用するためのライブラリ (libPZDIAG.a) を作成する makefile の例は make lib に書かれている。またパッケージ内の makefile には Fortran コンパイラに Intel コンパイラを利用するmpif90 を指定している。

4