Index of /staff/pdc/E
Name Last modified Size Description
![[DIR]](/icons/back.gif)
Parent Directory 06-Mar-2000 11:06 -
![[TXT]](/icons/text.gif)
HouseCom.f 26-Oct-1999 22:01 6k
Houselib.f 26-Oct-1999 22:01 3k
MakeMat.f 26-Oct-1999 22:02 5k
Makefile 26-Oct-1999 18:30 2k
PBisec.f 26-Oct-1999 22:03 8k
![[ ]](/icons/unknown.gif)
PEigenSolver.f 26-Oct-1999 22:03 11k
PHouseInvTrs.f 26-Oct-1999 22:04 4k
PTred.f 26-Oct-1999 22:04 16k
PTriInvCGS.f 26-Oct-1999 22:05 15k
PTriInvMGS.f 26-Oct-1999 22:06 14k
Redist.f 26-Oct-1999 22:06 9k
TestSolver.f 26-Oct-1999 22:07 9k
TestSym.f 26-Oct-1999 22:07 9k
TriDiagLU.f 26-Oct-1999 22:08 7k
mppf.f 26-Oct-1999 22:09 8k
mppf.h 26-Oct-1999 18:30 6k
(e) 並列汎用固有値問題ルーチン
●[入力データ仕様]
入力データは、標準対称固有値問題 A x = λ x における行列 A の
要素である。この要素をサイクリックサイクリック分散方式で各PEに分散
させ、それを入力データとする。
サイクリックサイクリック分散方式について、以下に具体的な説明を示す。
まず p個のプロセッサの一次元的な番号を pe_num とする。このとき
pe_num : 0, 1, ...., p-1
の範囲をとる。また、p の ルートの値を ncelx とする。このとき
p = ncelx * ncelx
とする。ここで一次元的な番号 pe_num から二次元的なプロセッサ番号
(ncidx, ncudy) への変換は、以下のサブプログラムで行なっている。
subroutine cid2xy(pe_num, ncelx, ncidx, ncidy)
integer pe_num, ncelx, ncidx, ncidy
ncidx = modulo(pe_num, ncelx)
ncidy = (pe_num - modulo(pe_num, ncelx))/ncelx
return
end
いま、二次元的なプロセッサ番号 (ncidx, ncidy) が得られた。
さて 行列 n×n行列 A を p プロセッサに分配するには以下のように
書く (ただし実装はしやすいが、あまり実行効率は良くない方法)。
いま行列 A の 添字は 0 〜 n-1 まで変化するとする。
do i=0,n-1
li = LOC(i)
if (ncidx .eq. OWNER(i)) then
do j=0,n-1
if (ncidy .eq. OWNER(j)) then
A(li, LOC(j)) = n - i
endif
enddo
endif
enddo
ここで、OWNER(i)は第i列を所有するプロセッサ群の二次元的な番号を
もつ配列、LOC(i)は第i列を所有するプロセッサのローカルな添字番号をもつ
配列である。このリストの作成方法は
do i=0, n-1
LOC(i) = CYCLIC_LOC(i,ncelx)
OWNER(i) = CYCLIC_OWNER(i,ncelx)
enddo
とする。関数 CYCLIC_LOC(i,ncelx)、CYCLIC_OWNER(i,ncelx)は
integer function CYCLIC_OWNER(i, p)
integer i, p
CYCLIC_OWNER = modulo(i, p)
return
end
integer function CYCLIC_LOC(i, p)
integer i, p
real*8 ii, pp
ii = i
pp = p
CYCLIC_LOC = idfloor(ii/pp)
return
end
このような関数である。
第 i列を所有するプロセッサの二次元的な番号の求め方が分かれば、
第 j列を有するプロセッサの二次元的な番号も同様なので、列方向に関する
所有情報は求める必要はない。
具体的に書くと (p=4 ,n=3)
行列 A OWNER(i) 所有者情報 (ncidx, ncidy)
→j
↓ a_11 a_12 a_13 0 (0,0) (0,1) (0,0)
i a_21 a_22 a_23 1 (1,0) (1,1) (1,0)
a_31 a_32 a_33 0 (0,0) (0,1) (0,0)
LOCAL(i) 内部添字情報(local_i, local_j)
0 (0,0) (0,0) (0,1)
0 (0,0) (0,0) (0,1)
1 (1,0) (1,0) (1,1)
PE pe_num , (ncidx, ncidy)
0 , (0, 0)
1 , (1, 0)
2 , (0, 1)
3 , (1, 1)
となる。こここで 所有者情報は 行列 a_{ij} を所有している２次元的なPE
番号を意味する。また 内部添字情報は 行列 a_{ij} を代入するべきローカル
な添字番号を示す。
例えば 0番のプロセッサ (0, 0) は、行列Aの要素
a_{11}, a_{31}, a_{13}, a_{33}
を所有していることが所有者情報からわかり、ローカルに代入すべき場所も
それぞれ
(0,0), (1,0), (0,1), (1,1)
に代入すべきことが 内部添字情報から分かる。
●[出力データ仕様]
出力データは、実数対称固有値問題における全実数固有値とそれに対応した
固有ベクトルである。固有値については、全てのプロセッサで全固有値を
格納する。固有ベクトルは各プロセッサで分散して格納される。プロセッサ
台数を p、問題サイズを n とすると、プロセッサ番号の若い順に
n/p 個づつ (割り切れない場合は切り上げ) 所有することとする。
以下にその例を示す。
出力データの例：
プロセッサ台数 p=4、固有ベクトル個数 n=5 ,(v_i, i=1,...,5) とする。
n/p = 2 なので、
プロセッサ番号 0の所有： v_1, v_2、
プロセッサ番号 1の所有： v_3,v_4
プロセッサ番号 2の所有： v_5、
プロセッサ番号 3の所有： なし、
のように格納される。
● [本サブルーチンの入出力仕様]
本サブルーチンは、ユーザが通信に関する処理を初期化した上で、
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
以下に示すように Fortran90のサブーチンコールで呼ぶことを前提とする。
call PEigenSolver(A, eig, W, n, nx, ny, p, sw, cid, errcode)
ここで 各引数の説明を以下に示す。
________________________________________________________________
________________________________________________________________
| |
引数名 | 入出力方向 | 説明
----------------------------------------------------------------
A | 入力 | real*8型。入力仕様に伴う処理データを格納する。
| | 要素番号は (0:nx-1，0:ny-1)。
| | 記憶内容は壊される。
-----------------------------------------------------------------
eig | 出力 | real*8型。全固有値格納用。
| | 要素番号は (1:n)。
-----------------------------------------------------------------
W | 出力 | real*8型。分散された全固有ベクトル値格納用。
| | 出力形式は出力仕様に伴う。
| | 要素番号は (1:n，1:n/p+1)。
-----------------------------------------------------------------
n | 入力 | integer型。問題サイズ。
-----------------------------------------------------------------
nx | 入力 | integer型。(問題サイズ) / (x方向のPE数+1)。
| | (全PE数 = x方向のPE数 × y方向のPE数)
-----------------------------------------------------------------
ny | 入力 | integer型。(問題サイズ) / (y方向のPE数+1)。
| | (全PE数 = x方向のPE数 × y方向のPE数)
-----------------------------------------------------------------
p | 入力 | integer型。全PE数。
-----------------------------------------------------------------
sw | 入力 | integer型。逆反復法における再直交処理方式の
| | 切り替えスイッチ。1:MGS法 2:CGS法
-----------------------------------------------------------------
cid | 入力 | integer型。自PEのPE番号 (0〜 p-1)。
-----------------------------------------------------------------
errcode| 出力 | integer型。エラーコード格納用。
-----------------------------------------------------------------
●[テストプログラムについて]
本サブルーチンの利用例として、テストプログラム TestSolver.f
が提供されている。以降にその利用法を述べる。
○ コンパイル方法
コマンドライン上で
|> make TestSolver
と入力する。
○ 実行方法
4プロセッサを用いて、パーテイション名 ALL にて実行する場合：
コマンドライン上で
|> mpirun -np 4 -part ALL TestSolver
と入力する。
○ データ入力方法
テストプログラムの起動に成功すると、以下のような表示がなされる。
|=== Parallel Eigensolver Check Program ===
|Parallel Eigensolver on MPI Ver. 1.00 1997/10
|Composed by KATAGIRI Takahiro
|
|Memory 230.79 MByte OK.
|
|problem size >
|?
ここで、 1000次元 の問題を CGS法を用いた直交化方式で解く場合、
結果を画面上に出力しない時は以下の様に入力する。
|problem size >
|?
|1000
|Inverse Iteration
|1)MGS 2)CGS >
|?
|2
|Show Eigenvale&Eigenvector(y/n) >
|?
|n
計算が正常になされるならば、以下のような結果が出力される。
|正常終了
|=== Result
|================================================
|n = 1000 p = 4
|Test Matrix: Frank
|Inverse Iteration : CGS
|
|Eigenvalue Max Error = 0.228663089789461758E-008
|Orthogonality = 0.836387922105658582E-013
|Max Remainder = 0.253391179023594428E-009
|Total Calculation Time = 75.842 [sec]
======
1998年 1月31日