当講習会の目標

第12回お試しアカウント付き
並列プログラミング講習会（試行）「ライブラリ利用：
高性能プログラミング初級入門」
PartⅡ

１．６０分後には，線形方程式に対して，１００種類
近くのアルゴリズムで求解できるようになる．
２．９０分後には，線形方程式に対して，並列計算による求解ができるようになる．
３．１８０分後には，反復解法に基づく求解アルゴ
リズムの特徴と注意点を理解している．

東京大学 情報基盤センター

こんな夢のようなことを実現するのが当講習会！

伊藤祥司
itosho@cc.u‐tokyo.ac.jp

目
前半の内容

その鍵となるのが・・・

後半の内容◆ Lis [ http://www.sscsc.org/ ]9. Lisの概観【41】1. HA8000での実習環境設定【5-7】10. 線形方程式用Lisの概要【42】2. Lisインストールの準備作業【8,9】11. 当実習で用いた評価条件【43】3. Lis逐次版のインストール【10-13】

12. 非定常反復法に関する盲点【44-53】

5. ジョブ実行（逐次版） 【15-17】→ Lis の「・詳細情報」をクリック→ Downloadの「Version 1.2.xx」からソースコード，ユーザマニュアル
をダウンロード※ ソースコードは次画面の「ダウンロードページ」→ 「Software License Agreement」へと進む．
ユーザ登録する／しないは各自の自由※ ただし，当資料，当講習会では，lis‐1.2.49 版を前提として説明，実習します．

13. PETScのインストールと使用要領【54-56】

6. Lis使用に関する解説【18-22】◆ テスト問題 Matrix Market [ http://math.nist.gov/MatrixMarket ]→ BrowseやSearchの項目のリンクをクリック

7. 入力データ形式について【23】8. Lis並列版の準備・インストール【24-39】2

0. 準備：Lisとテスト問題の入手

次

0. 準備：Lisとテスト問題の入手【4】4. 補足事項【14】LisとPETScだ！

【 】：当資料中のスライド番号34

1. HA8000での実習環境設定(1/3)

1. HA8000での実習環境設定(2/3)

各自の$HOME下に実習用sampleファイル一式をコピー：$ cd[Enter]$ cp –rp /home/t00003/sample$ ls -FR sample
sample:library/mat/qsub_script/

sample/library:lis-1.2.49.tar.gz
lis-manual-1.2.49_ja.pdf
sample/mat:494_bus.mtx
494_bus.rsa
add32.mtx.

setup.env$ cd sample
$ ls

sh/

test/

library

以後 [Enter]キー
押下の記述は省略

matqsub_script

setup.env

sh

test

# 実習環境のセットアップ$ source setup.env# 今回の実習で使用するライブラリ等

ログインの度に実行する（☆１）

petsc-lite-3.1-p3.tar.gz
petsc-manual-3.1.pdf# テスト用入力データ（係数行列，右辺項データ）

add32_rhs1.dat
add32_rhs1.mtx
add32_rhs2.mtx

add32.rua
fs_183_1.mtx
fs_183_1.rua

sample/qsub_script:
qsub_mpi.nqs qsub_omp.nqs

qsub_seq.nqs

sample/sh:
cdlis_kernel.sh*sample/test:
test1_mod2.c

各自の$HOME/sample下で実習用環境をセットアップ：# ホームディレクトリ($HOME)へ移動# 実習用sampleファイル一式をコピー

nos3.mtx
nos3.rsa# JOB実行用スクリプト# 今回の実習で役立つユーティリティ

morelis*# Lis実行で使用するmainルーチン

test1_mod.c# このような環境を設定$ cat setup.env

alias ls='ls -F'
alias lrt='ls -lrt'
alias cp='cp -i'
alias mv='mv -i'
alias rm='rm -i'
alias cdlis='source $HOME/sample/sh/cdlis_kernel.sh'
PATH=$PATH:$HOME/sample/sh##
20100912 created by itosho (ITOH, Shoji)#

test1_org.c● コマンド実行する部分は太字で表示してます（左端の$はUNIXのコマンドプロンプト）．
● 画面表示は，多少，文字サイズを小さくしてます（当テキスト編集の都合上）．
5● 特に注意の無い限り，# 記号から右部分の青文字はコメントです．

1. HA8000での実習環境設定(3/3)

2. Lisインストールの準備作業(1/2)

各自のNFS下へ実習環境を構築（T2K TODAI固有の環境セットアップ）：
$ cd
$ pwd

Lis実習環境の準備（後半のtest1.cに関する作業は，当実習ならではの事柄）：# ホームディレクトリ($HOME)へ移動$ tar zxf $HOME/sample/library/lis-1.2.49.tar.gz# その他の方法の一例として：（ $ gunzip -c $HOME/sample/library/lis-1.2.49.tar.gz | tar xf - ）# 各自のnfsのディレクトリをシンボリックリンク# /nfs/all/下への毎回の移動は面倒$ ls
lis-1.2.49//home/{各自の利用者番号}$ ln –s/nfs/all/$USER

nfs＃ 以後，{各自の利用者番号}を t00003 として説明

$ cd
$ ls

$ ls –l nfs
lrwxrwxrwx 1 t00003 t00 15 Sep 14 14:17 nfs -> /nfs/all/t00003/$ cd nfs$ pwd/home/t00003/nfs

6

lis-1.2.49/test/

defs*
etest1.cetest1f.F
etest2.cetest3.c# /home配下であることを確認$ lsetest4.cetest4f.F
etest5.c
Makefile.am
Makefile.in

spmvtest5.c
test1.ctest1f.F
test2.ctest3.ctest4.ctest4f.F
test5.ctest6.ctestmat.mtx

test.sh*$HOME/sample/test/

test1_mod2.c$ cp

matfiles.txt
spmvtest1.c
spmvtest2.c
spmvtest3.c
spmvtest4.c–p

test1_mod.c# test1_org.c と test1.c とは同一

test1_org.c$HOME/sample/test/test1_mod.c.

# ここでは test1_mod.c のみコピー$ mv
test1.ctest1_org.c$ ln –s test1_org.c test1.c# サンプルの test1.c を準備$ ls –l test1.c
lrwxrwxrwx 1 t00003 t00 11 Sep 14 14:55 test1.c -> test1_org.c
7$ cd$HOME/nfs

8

2. Lisインストールの準備作業(2/2)

3. Lis逐次版のインストール(1/3)
Lis逐次(seq)版のconfigure：

Lis逐次(seq)版の準備：
$ ls

$ pwd

lis-1.2.49//home/t00003/nfs/lis-1.2.49_seq$ cp –rp lis-1.2.49 lis-1.2.49_seq# Lis逐次(seq)版の環境準備$ ./configure
・・・ 途中省略 ・・・$ cd lis-1.2.49_seq$ lsaclocal.m4
AUTHORS
ChangeLog

config/configure*configure.in--enable-saamg# SAAMG前処理を有効にする● Lisを正式にインストールする場合，--prefix オプションでディレクトリも指定．

・・・ 出力の最後部 ・・・
COPYINGdoc/
include/

INSTALL
Makefile.am
Makefile.in

NEWS
README
src/

test/
win32/● 以後，コンパイラ等の処理系のバージョンアップに起因して，出力情報が異なる場合もありますが，それらは適宜読み替えて下さい．
● ただし，コンパイラ等のバージョンが異なることによる非互換などはあり得ます．
9

Using
C
Compiler
= icc
C
Flags
= -DHAVE_CONFIG_H -O3 -ansi_alias
赤文字の５項目
C
Libs
=
Libs
=
-lm
をチェック
F90 Compiler= ifortF90 Flags= -O3 -fpp -g -DZERO_ORIGIN=1
F90 Libs=
enable Fortran Wrapper= noenable SA-AMG Preconditioner= yesenable MPI= noenable OpenMP= noenable Quad Precision
= no
MPICC
=
MPIF77=● configure実行後は，作業ディレクトリやパスを
MPIFC=
変更しないように（生成されるMakefile中では絶対
MPILIBS=MPIRUN=
パスが用いられている）．
MPFLAG=
10

3. Lis逐次版のインストール(2/3)

3. Lis逐次版のインストール(3/3)＃ 今回の実習では，Intelコンパイラを使用する$ source /opt/itc/mpi/mpiswitch.sh$ which icc ifort/opt/intel/cce/10.1.018/bin/icc/opt/intel/fce/10.1.018/bin/ifort

mpich-mx-intel

ログインの度に実行する（☆２）# Intelコンパイラの環境設定# Intelコンパイラのiccのパス# Intelコンパイラのifortのパス

Lis逐次(seq)版のmake：Lis逐次(seq)版のインストールチェック(make check)：$ make check$ make・・・

途中省略

・・・ 出力の途中部分 ・・・・・・・・・ 出力の最後部 ・・・/bin/sh ../libtool --mode=link ifort -Vaxlib -nofor_main -o spmvtest5 -lm spmvtest5.o L../src -llis -lm
libtool: link: ifort -Vaxlib -nofor_main -o spmvtest5 spmvtest5.o -L/home/t00003/nfs/lis1.2.49_seq/src /home/t00003/nfs/lis-1.2.49_seq/src/.libs/liblis.a -lm
make[1]: Leaving directory `/nfs/all/t00003/lis-1.2.49_seq/test'
Making all in win32
make[1]: Entering directory `/nfs/all/t00003/lis-1.2.49_seq/win32'
make[1]: Nothing to be done for `all'.
make[1]: Leaving directory `/nfs/all/t00003/lis-1.2.49_seq/win32'
make[1]: Entering directory `/nfs/all/t00003/lis-1.2.49_seq'
make[1]: Nothing to be done for `all-am'.
make[1]: Leaving directory `/nfs/all/t00003/lis-1.2.49_seq'

make check-TESTS
make[2]: Entering directory `/nfs/all/t00003/lis-1.2.49_seq/test‘=== Running test test.shchecking solvers...
number of processes = 1
100 x 100 matrix 460 entries
Initial vector x = 0
PRECISION : DOUBLESOLVER: BiCG 2
PRECON: NoneSTORAGE: CRS
lis_solve : normal end

赤文字の項目をチェック
下線の項目は計測時間なので，一致して無くても良い

BiCG: iter
= 15 iter_double = 15 iter_quad = 0
BiCG: times= 2.660751e-04
BiCG: p_times = 1.788139e-05 (p_c = 2.145767e-06 p_i = 1.573563e-05 )BiCG: i_times = 2.481937e-04
BiCG: residual = 6.399933e-15

最後部でErrorが出ていなければOK

checking eigensolvers...
11

・・・

12

4. 補足事項：main関数の入替え

3. Lis逐次版のインストール(補足)
当実習の中では，make installは実行しなくても利用可能：

当実習ならではの事柄だが，反復法の収束判定に関する重要事項：$ cd test$ ls -l test1*.c$ make install
マニュアルには，最後に上記を実行するように記述されているが，これは，正式にシステムにインストールする（例えば，/usr/local/下など）場合に実行すれば良い．
この実習では，既に作成されている “test1 ” というプログラムを用いて説明するため，特にインストール先を指定しないでも良い．
現状で，生成済みの関連ファイル群は下記のディレクトリ構成となっている：

lrwxrwxrwx 1 t00003 t00
11 Sep 14 16:43 test1.c -> test1_org.c-rw-r--r-- 1 t00003 t00 5309 Sep 12 16:01 test1_mod.c-rw-r--r-- 1 t00003 t00 5157 Jun 18 10:48 test1_org.c# test1.cのリンク先$ mv test1 test1_org ; ls -l test1_org# 必須ではないが-rwxr-xr-x 1 t00003 t00 1669881 Sep 14 16:50 test1_org*$ rm test1.c test1.o ; ln –s test1_mod.c test1.c# main関数の入替え

rm: remove symbolic link `test1.c'? y
rm: remove regular file `test1.o'? y$HOME/nfs/lis-1.2.49_seq/ 下：
include/config.h，lis.h，lisf.h
src/.libs/
liblis.a$ ls -l test1*.c
lrwxrwxrwx 1 t00003 t00
11 Sep 14 16:53 test1.c -> test1_mod.c-rw-r--r-- 1 t00003 t00 5309 Sep 12 16:01 test1_mod.c-rw-r--r-- 1 t00003 t00 5157 Jun 18 10:48 test1_org.c# test1.cのリンク先$ make

インストール先のデフォルトは /usr/local/ 下であるが，それを変更する場合には，configure のステップまで戻って，icc -DHAVE_CONFIG_H -I. -I../include -I.. -I../include-O3 -ansi_alias -c test1.c/bin/sh ../libtool --mode=link ifort -Vaxlib -nofor_main -o test1 -lm test1.o -L../src-llis -lm
libtool: link: ifort -Vaxlib -nofor_main -o test1 test1.o -L/home/t00003/nfs/lis1.2.49_seq/src /home/t00003/nfs/lis-1.2.49_seq/src/.libs/liblis.a –lm$ ./configure --prefix={インストールディレクトリ}とする．

$ ls -l test1 test1_org
13

5. ジョブ実行スクリプトの用意(逐次版)$HOME/sample/qsub_script/

qsub_mpi.nqs

qsub_omp.nqs

qsub_seq.nqs$ cp –p $HOME/sample/qsub_script/qsub_seq.nqs .
$ cat qsub_seq.nqs#@$-q lecture
講習会後のキュー名．本日は#@$-N 1tutorial (sampleは変更済み)#@$-lm 2GB#@$-lT 0:01:00#@$-o OUT#@$-e ERR#@$-s /bin/bash# 標準出力のファイル名# 標準エラー出力のファイル名# スクリプト内のシェルを指定

cd $PBS_O_WORKDIR
pwd# スクリプト内もカレントディレクトリへ移動する# pwdは本来は不要だが，確認のため

5. HA8000でのジョブ実行(逐次版)$ qsub qsub_seq.nqs● NQSのスクリプトでは左端の“ #@$ ”記号は，これに続くハイフンで始まる文字列がNQSオプ
ションであることを宣言している．
● 赤文字の項目をチェック（NQSオプション“-N”でノード数を指定している）./test1 testmat.mtx 0 sol1.txt res1.txt# 求解1（デフォルト： -i bicg，-p none）./test1 testmat.mtx 0 sol2.txt res2.txt -i 3 -p 1
# 求解2#./test1 ./mat/add32.mtx 2 sol.txt res.txt -i 5 -p 2#./test1 ./mat/add32.mtx ./mat/add32_rhs1.dat sol.txt res.txt -i 5 -p 2#./test1 ./mat/add32.rua 2 sol.txt res.txt -i 5 -p 2#./test1 ./mat/nos3.mtx 2 sol.txt res.txt -i 5 -p 2
# see manual subsection 3.4##
20100914 created by itosho (ITOH, Shoji)#

14

ジョブ実行（qsub）と実行中・実行後の状態確認方法：

HA8000では，PBS（NQSの一種）でジョブ実行するスクリプトを使用する：$ ls-rwxr-xr-x 1 t00003 t00 1669905 Sep 14 16:54 test1*-rwxr-xr-x 1 t00003 t00 1669881 Sep 14 16:50 test1_org*

15

Request 606870.batch1 submitted to queue: lecture.

リクエストID（ジョブ固有の番号）$ qstat# ジョブ実行中の表示例

2010/09/14 (Tue) 19:30:26:
REQUESTS on HA8000 cluster
NQS schedule stop time : 2010/09/25 (Sat) 8:55:00 (Remain: 253h 24m 34s)
REQUEST
NAME
OWNER
QUEUE
PRI NODE E-TIME
MEM
STATE
606870.batch1
qsub_set00003lecture63160s
28GB RUNNING$ qstat# ジョブ終了後の表示例

2010/09/14 (Tue) 19:30:37:
REQUESTS on HA8000 cluster
NQS schedule stop time : 2010/09/25 (Sat) 8:55:00 (Remain: 253h 24m 23s)
REQUEST
NAME
OWNER
QUEUE
PRI NODE E-TIME
MEM
STATE
No requests.

$ type lrt
lrt is aliased to `ls -lrt‘$ lrt
・・・ 省略 ・・・-rw------- 1 t00003 t00-rw-r--r-- 1 t00003 t00-rw-r--r-- 1 t00003 t00-rw-r--r-- 1 t00003 t00-rw-r--r-- 1 t00003 t00-rw------- 1 t00003 t00

0 Sep 14 19:30 ERR# 標準エラー出力のファイル
3242 Sep 14 19:30 sol2.txt # 解ベクトル2(solution)
3242 Sep 14 19:30 sol1.txt # 解ベクトル1(solution)
208 Sep 14 19:30 res2.txt # 相対残差2(residual)のノルム収束履歴
208 Sep 14 19:30 res1.txt # 相対残差1(residual)のノルム収束履歴
864 Sep 14 19:30 OUT# 標準出力のファイル
165. 実行したジョブの結果確認(逐次版)$ cat OUT

※ 編集の都合上，出力の一部を折り返し表示している．

6. いま，Lisで何をやったのだろうか$ ./test1/home/t00003/nfs/lis-1.2.49_seq/test# 確認のpwdの出力
number of processes = 1
100 x 100 matrix 460 entries
Initial vector x = 0
PRECISION : DOUBLESOLVER: BiCG 2
PRECON: NoneSTORAGE: CRS
lis_solve : normal end

number of processes = 1
100 x 100 matrix 460 entries
Initial vector x = 0
PRECISION : DOUBLESOLVER: CGS 3
PRECON: Jacobi
STORAGE: CRS
lis_solve : normal end

係数行列のファイルを指定
右辺項ベクトルファイル，あるいは，自動生成機能（主にアルゴリズムの性能評価用途）のIDを指定
0 : 拡張Matrix Market形式（右辺ベクトルが行列データに含まれている）を用いる
1 : b = (1,1,...,1)T を用いる
2 : b = A ×(1,1,...,1)T を用いる
solution_filename: 解ベクトルのファイル名
residual_filename: （アルゴリズム中のベクトルを用いた）相対残差ベクトルの２ノルムの履歴[options]:
解法と前処理，その他のオプション指定

BiCG: iter
= 15 iter_double = 15
iter_quad = 0
BiCG: times= 1.130104e-04
BiCG: p_times = 1.811981e-05 (p_c =9.536743e-07 p_i = 1.716614e-05 )BiCG: i_times = 9.489059e-05
BiCG: residual = 6.399933e-15

CGS: iter= 15 iter_double = 15
iter_quad = 0
CGS: times= 1.420975e-04
CGS: p_times = 3.552437e-05 (p_c =9.536743e-07 p_i = 3.457069e-05 )CGS: i_times = 1.065731e-04
CGS: residual = 1.091864e-14# 参考： マニュアル 「テストプログラム」の「test1」小節

BiCG: rel_resid = 9.237507e-16

CGS: rel_resid = 1.575969e-15./test1 testmat.mtx 0 sol1.txt res1.txt# 求解1の結果（デフォルト： -i bicg，-p none）

Usage: test1 matrix_filename rhs_setting solution_filename residual_filename [options]matrix_filename:rhs_setting:$ more testmat.mtx./test1 testmat.mtx 0 sol2.txt res2.txt -i 3 -p 1
# 求解2の結果
17

6. 色々な解法で解いてみよう%%MatrixMarket matrix coordinate real general
100 100 460 1 0
1 2 -1.00000000000000000000e+000# 参考1： マニュアル 付録「ファイル形式」の1 11 -1.00000000000000000000e+000「拡張Matrix Market形式」小節
1 1 4.00000000000000000000e+000
testmat.mtx も係数行列と右辺項の両方の
2 3 -1.00000000000000000000e+000
2 12 -1.00000000000000000000e+000
データを格納した拡張形式．
2 1 -1.00000000000000000000e+000# 参考2：
2 2 4.00000000000000000000e+000
http://math.nist.gov/MatrixMarket/formats.html
・・・

18

6. 色々な前処理を組合せてみよう

解法

選択オプション

補助オプション

前処理

選択オプション

補助オプション

デフォルトはBiCG{ a | n} : a，nどちらでも良い[ ]内はデフォルト値

デフォルトは“なし”{ a | n} : a，nどちらでも良い[ ]内はデフォルト値‐i {bicg | 2}なし
Jacobi‐p {jacobi | 1}CGS

‐i {cgs | 3}

ILU( k )‐p {ilu | 2}

‐ilu_fill [0]フィルインレベル k

BiCGSTAB‐i {bicgstab | 4}SSOR

‐p {ssor | 3}

‐ssor_w [1.0]緩和係数ωの値(0<ω<2)

BiCGSTAB( l )‐i {bicgstabl | 5}

Hybrid‐p {hybrid | 4}

‐hybrid_i [sor] 線型方程式系解法，他

CG‐i {cg | 1}

BiCG

‐ell [2]解法の l の値

‐p {none | 0}

GPBiCG‐i {gpbicg | 6}

I+S

‐p {is | 5}

‐is_alpha [1.0] 前処理中の変数α，他

TFQMR‐i {tfqmr | 7}

SAINV

‐p {sainv | 6}

‐sainv_drop [0.05]

Orthomin(m)‐i {orthomin | 8}

‐restart [40]リスタート m の値

SA‐AMG‐p {saamg | 7}

‐saamg_unsym [false] 非対称版選択，他

GMRES(m)‐i {gmres | 9}

‐restart [40]リスタート m の値

Crout ILU

‐p {iluc | 8}

‐iluc_drop [0.05]

ドロップ基準，他

Jacobi‐i {jacobi | 10}ILUT‐p {ilut | 9}

‐ilut_drop [0.05]

ドロップ基準，他

Gauss‐Seidel‐i {gs | 11}

Additive Schwarz‐adds true

‐adds_iter [1]

繰り返し回数SOR‐i {sor | 12}‐omega [1.9]● 上記も含め，全22種類の解法が提供されている．

ドロップ基準緩和係数ωの値(0<ω<2)19● 上記の補助オプションは一部のみ掲載している．詳細はマニュアル参照．

20

6. 色々なオプションを試してみよう
オプション

内容

$ pwd

[ ]内はデフォルト値/home/t00003/nfs/lis-1.2.49_seq/test‐maxiter [1000]

最大反復回数‐tol [1.0e‐12]

収束判定基準（相対残差）‐print [0]

残差の画面表示‐scale [0]スケーリングの選択．
※ スケーリングされた係数行列，右辺項ベクトルは，元の‐initx_zeros [true]$ ls‐storage_block [2]

$ ln
$ ls

mat/

-s
-l

qsub_script/

setup.env

sh/

test/

$HOME/sample/mat .
matlrwxrwxrwx 1 t00003 t00 23 Sep 14 19:39 mat -> /home/t00003/sample/mat/$ ls mat/494_bus.mtx
494_bus.rsa
add32.mtx

データに上書きされて，プログラムで実行される．
初期ベクトル x0 の設定（デフォルトは全要素の値が0）

add32_rhs1.dat
add32_rhs1.mtx
add32_rhs2.mtx

add32.rua
fs_183_1.mtx
fs_183_1.ruanos3.mtx
nos3.rsa

参考までに，matの後にディレクトリを示す「/」を付け忘れると：$ ls mat# これは，☆１のsetup.env 影響● 拡張子
mtx : Matrix Market形式
rua, rsa : Harwell-Boeing形式
dat : ベクトルのPLAIN形式(拡張子自体は任意)

mat@

行列格納形式
BSR，BSCのブロックサイズ● 上記のオプションの内容の詳細はマニュアル参照．

$HOME/sample/

library/‐omp_num_threads [t] 実行スレッド数．t は最大スレッド数‐storage [0]6. 色々なデータで実行してみよう$ cat qsub_seq.nqs#./test1 ./mat/add32.mtx 2 sol.txt res.txt -i 5 -p 2#./test1 ./mat/add32.mtx ./mat/add32_rhs1.dat sol.txt res.txt -i 5 -p 2#./test1 ./mat/add32.rua 2 sol.txt res.txt -i 5 -p 2#./test1 ./mat/nos3.mtx 2 sol.txt res.txt -i 5 -p 2
21

7. 入力データ形式について● 各々を実行するときは，左端冒頭の「#」（コメントアウト）記号を削除する．

22

8. Lis並列(OMP, MPI)版の準備作業
逐次版と同様，main関数のtest1.c を入替える（当実習ならではの事柄）：$ cd $HOME/nfs/lis-1.2.49/test/$ ls -l test1*.c● Matrix Market(MM)形式は，データの圧縮度は低いものの，行列のインデックスと値との対応を確認し易く，また，多くのアプリケーションでもMM形式をサポートしたりツール提供されているので，初級者や中級者にとっては都合が良いと思われる．
● Lisのtest1.cでは，MM形式のデータを読込み，CRS形式に変換する．MPIのデータ分散にも対応している．
● 対称行列の場合には，対角要素から下三角の要素のみが格納されている．この場合，ヘッダには，半角スペースに続いて「symmetric」の文字列が記載されている．
→ mat/下の「494_bus.mtx」「nos3.mtx」が該当

lrwxrwxrwx 1 t00003 t00
11 Sep 14 14:55 test1.c -> test1_org.c-rw-r--r-- 1 t00003 t00 5309 Sep 12 16:01 test1_mod.c-rw-r--r-- 1 t00003 t00 5157 Jun 18 10:48 test1_org.c$ rm test1.c ; ln -s test1_mod.c test1.c# 逐次版同様に test1.c を入替える

rm: remove symbolic link `test1.c'? y$ ls -l test1*.c
lrwxrwxrwx 1 t00003 t00
11 Sep 15 15:06 test1.c -> test1_mod.c-rw-r--r-- 1 t00003 t00 5309 Sep 12 16:01 test1_mod.c-rw-r--r-- 1 t00003 t00 5157 Jun 18 10:48 test1_org.c● 非対称行列の場合には，「general」の文字列が記載されている．
→ mat/下の「add32.mtx」「fs_183_1.mtx」が該当$ ln -s $HOME/sample/mat .
$ cd $HOME/nfs$ cp –rp lis-1.2.49 lis-1.2.49_omp$ cp –rp lis-1.2.49 lis-1.2.49_mpi

add32.mtx: MMサイトで提供されているMM形式の係数行列データ
add32_rhs1.mtx :〃MM形式の右辺項ベクトルデータ
add32.rua:〃
Harwell-Boeing形式の係数行列データ# ついでにmatもシンボリックリンク# LisのOpenMP(omp)版の環境準備# LisのMPI(mpi)版の環境準備$ which icc ifort

実は，MMサイトで提供されている add32_rhs1.mtx のヘッダとフォーマットがおかしいようなので，当実習では以下の２通りの修正データを用意してます．
問題点と修正内容は，各自で確認して下さい．
add32_rhs1.dat : add32_rhs1.mtxを元にPLAIN形式にしたファイル
add32_rhs2.mtx : 本来のMM形式のベクトルデータ# lis-1.2.49_seq ではないことに注意/opt/intel/cce/10.1.018/bin/icc/opt/intel/fce/10.1.018/bin/ifort# Intelコンパイラのiccのパス# Intelコンパイラのifortのパス＃ Intelコンパイラの環境設定されていない場合には：$ source /opt/itc/mpi/mpiswitch.sh mpich-mx-intel
23

$ cd

lis-1.2.49_omp# Intelコンパイラの環境設定
24

8. Lis並列(OMP)版のインストール(1/3)
Lis並列(omp)版のconfigure：

8. Lis並列(OMP)版のインストール(2/3)
Lis並列(omp)版のmake：$ pwd$ make/home/t00003/nfs/lis-1.2.49_omp$ ./configure
・・・ 途中省略 ・・・--enable-omp --enable-saamg# ここでompを指定する● Lisを正式にインストールする場合，--prefix オプションでディレクトリも指定．

・・・ 出力の最後部 ・・・
Using
C
Compiler
= icc
C
Flags
= -DHAVE_CONFIG_H -O3 -ansi_alias -openmp
赤文字の６項目
C
Libs
=
Libs
=
-lm
をチェック
F90 Compiler= ifortF90 Flags= -O3 -fpp -g -DZERO_ORIGIN=1 -openmp
F90 Libs=
enable Fortran Wrapper= noenable SA-AMG Preconditioner= yesenable MPI= noenable OpenMP= yesenable Quad Precision
= no
MPICC
=
MPIF77=● configure実行後は，作業ディレクトリやパスを
MPIFC=
変更しないように（生成されるMakefile中では絶対
MPILIBS=MPIRUN=
パスが用いられている）．
MPFLAG=25

8. Lis並列(OMP)版のインストール(3/3)$ make check・・・/bin/sh ../libtool --mode=link ifort -Vaxlib -openmp -nofor_main -o spmvtest5 -lm
spmvtest5.o -L../src -llis -lm
libtool: link: ifort -Vaxlib -openmp -nofor_main -o spmvtest5 spmvtest5.o L/home/t00003/nfs/lis-1.2.49_omp/src /home/t00003/nfs/lis-1.2.49_omp/src/.libs/liblis.a lm
make[1]: Leaving directory `/nfs/all/t00003/lis-1.2.49_omp/test'
Making all in win32
make[1]: Entering directory `/nfs/all/t00003/lis-1.2.49_omp/win32'
make[1]: Nothing to be done for `all'.
make[1]: Leaving directory `/nfs/all/t00003/lis-1.2.49_omp/win32'
make[1]: Entering directory `/nfs/all/t00003/lis-1.2.49_omp'
make[1]: Nothing to be done for `all-am'.
make[1]: Leaving directory `/nfs/all/t00003/lis-1.2.49_omp'

最後部でErrorが出ていなければOK26

8. ジョブ実行スクリプトの用意(OMP版)$ cd test$ ls $HOME/sample/qsub_script/● NQSのスクリプトでは左端の“ #@$ ”記号は，これに続くハイ$ cp –p $HOME/sample/qsub_script/qsub_omp.nqs . フンで始まる文字列がNQSオプ$ cat qsub_omp.nqs
ションであることを宣言している．
#@$-q lecture・・・ 出力の途中部分 ・・・
make check-TESTS
make[2]: Entering directory `/nfs/all/t00003/lis-1.2.49_omp/test'=== Running test test.shqsub_mpi.nqs

checking solvers...

赤文字の項目をチェック
下線の項目は計測時間なので，一致して無くても良いqsub_omp.nqs#@$-N 1#@$-lm 2GB#@$-lT 0:01:00#@$-o OUT#@$-e ERR#@$-s /bin/bash

qsub_seq.nqs# 標準出力のファイル名# 標準エラー出力のファイル名# スクリプト内のシェルを指定● 赤文字の項目をチェック（NQSオプション“-N”でノード数を指定している）

export OMP_NUM_THREADS=4 # 環境変数によるスレッド数の指定
cd $PBS_O_WORKDIR# スクリプト内もカレントディレクトリへ移動
pwd# pwdは本来は不要だが，確認のため./test1 testmat.mtx 0 sol1.txt res1.txt# 求解1（デフォルト： -i bicg，-p none）./test1 testmat.mtx 0 sol2.txt res2.txt -i 3 -p 1
# 求解2#./test1 ./mat/add32.mtx 2 sol.txt res.txt -i 5 -p 2#./test1 ./mat/add32.mtx ./mat/add32_rhs1.dat sol.txt res.txt -i 5 -p 2#./test1 ./mat/add32.rua 2 sol.txt res.txt -i 5 -p 2#./test1 ./mat/nos3.mtx 2 sol.txt res.txt -i 5 -p 2

BiCG: iter
= 15 iter_double = 15 iter_quad = 0
BiCG: times= 1.162767e-03
BiCG: p_times = 1.111031e-04 (p_c = 1.907349e-06 p_i = 1.091957e-04 )BiCG: i_times = 1.051664e-03
BiCG: residual = 8.589712e-15
BiCG: rel_resid = 1.239818e-15
・・・

途中省略

・・・ 出力の最後部 ・・・

HA8000では，PBS（NQSの一種）でジョブ実行するスクリプトを使用する：

Lis並列(omp)版のインストールチェック(make check)：number of processes = 1
max number of threads = 16
number of threads = 2
100 x 100 matrix 460 entries
Initial vector x = 0
PRECISION : DOUBLESOLVER: BiCG 2
PRECON: NoneSTORAGE: CRS
lis_solve : normal end

・・・

27

# see manual subsection 3.4#
20100914 created by itosho (ITOH, Shoji)28

8. HA8000でのジョブ実行(OMP版)

8. 実行したジョブの結果確認(OMP版)

ジョブ実行（qsub）と実行中・実行後の状態確認方法：※ 編集の都合上，出力の一部を折り返し表示している．

$ qsub qsub_omp.nqs$ cat OUT

Request 607327.batch1 submitted to queue: lecture.

/home/t00003/nfs/lis-1.2.49_omp/test# 確認のpwdの出力
number of processes = 1
max number of threads = 16
number of threads = 4
100 x 100 matrix 460 entries
Initial vector x = 0
PRECISION : DOUBLESOLVER: BiCG 2
PRECON: NoneSTORAGE: CRS
lis_solve : normal end

number of processes = 1
max number of threads = 16
number of threads = 4
100 x 100 matrix 460 entries
Initial vector x = 0
PRECISION : DOUBLESOLVER: CGS 3
PRECON: Jacobi
STORAGE: CRS
lis_solve : normal end

BiCG: iter
= 15 iter_double = 15
iter_quad = 0
BiCG: times= 6.790161e-04
BiCG: p_times = 8.130074e-05 (p_c =9.536743e-07 p_i = 8.034706e-05 )BiCG: i_times = 5.977154e-04
BiCG: residual = 7.177916e-15

CGS: iter= 15 iter_double = 15
iter_quad = 0
CGS: times= 7.081032e-04
CGS: p_times = 1.065731e-04 (p_c =9.536743e-07 p_i = 1.056194e-04 )CGS: i_times = 6.015301e-04
CGS: residual = 1.359740e-14

BiCG: rel_resid = 1.036043e-15

CGS: rel_resid = 1.962616e-15# ジョブ実行中の表示例$ qstat2010/09/15 (Wed) 16:04:56:
REQUESTS on HA8000 cluster
NQS schedule stop time : 2010/09/25 (Sat) 8:55:00 (Remain: 232h 50m 4s)
REQUEST
NAME
OWNER
QUEUE
PRI NODE E-TIME
MEM
STATE
607327.batch1
qsub_omt00003lecture63160s
28GB RUNNING# ジョブ終了後の表示例$ qstat2010/09/15 (Wed) 16:05:06:
REQUESTS on HA8000 cluster
NQS schedule stop time : 2010/09/25 (Sat) 8:55:00 (Remain: 232h 49m 54s)
REQUEST
NAME
OWNER
QUEUE
PRI NODE E-TIME
MEM
STATE
No requests.

$ type lrt
lrt is aliased to `ls -lrt‘$ lrt
・・・ 省略 ・・・-rw------- 1 t00003 t00-rw-r--r-- 1 t00003 t00-rw-r--r-- 1 t00003 t00-rw-r--r-- 1 t00003 t00-rw-r--r-- 1 t00003 t00-rw------- 1 t00003 t00

0 Sep 15 16:04 ERR# 標準エラー出力のファイル
3242 Sep 15 16:04 sol2.txt # 解ベクトル2(solution)
3242 Sep 15 16:04 sol1.txt # 解ベクトル1(solution)
208 Sep 15 16:04 res2.txt # 相対残差2(residual)のノルム収束履歴
208 Sep 15 16:04 res1.txt # 相対残差1(residual)のノルム収束履歴
962 Sep 15 16:04 OUT# 標準出力のファイル
29

8. Lis並列(MPI)版のインストール(1/3)
Lis並列(mpi)版のconfigure：$ pwd./test1 testmat.mtx 0 sol1.txt res1.txt# 求解1の結果（デフォルト： -i bicg，-p none）./test1 testmat.mtx 0 sol2.txt res2.txt -i 3 -p 1
# 求解2の結果
30

8. Lis並列(MPI)版のインストール(2/3)
Lis並列(mpi)版のmake：$ make/home/t00003/nfs/lis-1.2.49_mpi$ ./configure
・・・ 途中省略 ・・・--enable-mpi --enable-saamg# ここでmpiを指定する● Lisを正式にインストールする場合，--prefix オプションでディレクトリも指定．

・・・ 出力の最後部 ・・・
Using
C
Compiler= mpiccC
Flags= -DHAVE_CONFIG_H -O3 -ansi_alias -DUSE_MPI
赤文字の１０項目
C
Libs
=
Libs
=
-lm
をチェック
F90 Compiler= mpif90
F90 Flags= -O3 -fpp -g -DZERO_ORIGIN=1 -DUSE_MPI
F90 Libs=
enable Fortran Wrapper= noenable SA-AMG Preconditioner= yesenable MPI= yesenable OpenMP= noenable Quad Precision
= no
MPICC
= mpicc
MPIF77= mpif90● configure実行後は，作業ディレクトリやパスを
MPIFC= mpif90
変更しないように（生成されるMakefile中では絶対
MPILIBS=MPIRUN= mpirun
パスが用いられている）．
MPFLAG
=
31

・・・

途中省略

・・・

・・・ 出力の最後部 ・・・/bin/sh ../libtool --mode=link mpif90 -Vaxlib -nofor_main -o spmvtest5 -lm spmvtest5.o L../src -llis -lm
libtool: link: mpif90 -Vaxlib -nofor_main -o spmvtest5 spmvtest5.o L/home/t00003/nfs/lis-1.2.49_mpi/src /home/t00003/nfs/lis-1.2.49_mpi/src/.libs/liblis.a lm
make[1]: Leaving directory `/nfs/all/t00003/lis-1.2.49_mpi/test'
Making all in win32
make[1]: Entering directory `/nfs/all/t00003/lis-1.2.49_mpi/win32'
make[1]: Nothing to be done for `all'.
make[1]: Leaving directory `/nfs/all/t00003/lis-1.2.49_mpi/win32'
make[1]: Entering directory `/nfs/all/t00003/lis-1.2.49_mpi'
make[1]: Nothing to be done for `all-am'.
make[1]: Leaving directory `/nfs/all/t00003/lis-1.2.49_mpi'

最後部でErrorが出ていなければOK32

8. Lis並列(MPI)版のインストール(3/3)

8. ジョブ実行スクリプトの用意(MPI版)
HA8000では，PBS（NQSの一種）でジョブ実行するスクリプトを使用する：

Lis並列(mpi)版のインストールチェック(make check)：$ make check$ cd test$ ls $HOME/sample/qsub_script/● NQSのスクリプトでは左端の“ #@$ ”記号は，これに続くハイ$ cp –p $HOME/sample/qsub_script/qsub_mpi.nqs . フンで始まる文字列がNQSオプ$ cat qsub_mpi.nqs
ションであることを宣言している．
#@$-q lecture・・・ 出力の最後部分 ・・・
checking solvers...
mpiexec: Error: PBS_JOBID not set in environment. Code must be run from a
PBS script, perhaps interactively using "qsub -I".
checking eigensolvers...
mpiexec: Error: PBS_JOBID not set in environment. Code must be run from a
PBS script, perhaps interactively using "qsub -I".

qsub_mpi.nqs

checking SAAMG preconditioner...
mpiexec: Error: PBS_JOBID not set in environment. Code must be run from a
PBS script, perhaps interactively using "qsub -I".
FAIL: test.sh================================
MPIプログラムはqsubを用いな
1 of 1 test failed
いとシステム的に実行できない
Please report to devel@ssisc.org================================
のでErrorになって当然
make[2]: *** [check-TESTS] Error 1
make[2]: Leaving directory `/nfs/all/t00003/lis-1.2.49_mpi/test'
make[1]: *** [check-am] Error 2
make[1]: Leaving directory `/nfs/all/t00003/lis-1.2.49_mpi/test'
make: *** [check-recursive] Error 1$

qsub_omp.nqs

qsub_seq.nqs#@$-N 4#@$-lm 2GB#@$-lT 0:01:00#@$-o OUT#@$-e ERR#@$-s /bin/bash# 標準出力のファイル名# 標準エラー出力のファイル名# スクリプト内のシェルを指定

cd $PBS_O_WORKDIR
pwd# スクリプト内もカレントディレクトリへ移動# pwdは本来は不要だが，確認のため● 赤文字の項目をチェック（NQSオプション“-N”でノード数を指定している）● MPIジョブは mpirun コマンドの引数に実行プログラムを記述
する．

mpirun ./test1 testmat.mtx 0 sol1.txt res1.txt
mpirun ./test1 testmat.mtx 0 sol2.txt res2.txt -i 3 -p 1
#mpirun ./test1 ./mat/add32.mtx 2 sol.txt res.txt -i 5 -p 2
#mpirun ./test1 ./mat/add32.mtx ./mat/add32_rhs1.dat sol.txt res.txt -i 5 -p 2
#mpirun ./test1 ./mat/add32.rua 2 sol.txt res.txt -i 5 -p 2
#mpirun ./test1 ./mat/nos3.mtx 2 sol.txt res.txt -i 5 -p 2
# see manual subsection 3.4#
20100914 created by itosho (ITOH, Shoji)

33

8. HA8000でのジョブ実行(MPI版)

34

8. 実行したジョブの結果確認(MPI版)

ジョブ実行（qsub）と実行中・実行後の状態確認方法：※ 編集の都合上，出力の一部を折り返し表示している．

$ qsub qsub_mpi.nqs$ cat OUT

Request 607391.batch1 submitted to queue: lecture.

/home/t00003/nfs/lis-1.2.49_mpi/test# 確認のpwdの出力
number of processes = 4
100 x 100 matrix 460 entries
Initial vector x = 0
PRECISION : DOUBLESOLVER: BiCG 2
PRECON: NoneSTORAGE: CRS
lis_solve : normal end

number of processes = 4
100 x 100 matrix 460 entries
Initial vector x = 0
PRECISION : DOUBLESOLVER: CGS 3
PRECON: Jacobi
STORAGE: CRS
lis_solve : normal end

BiCG: iter
= 15 iter_double = 15
iter_quad = 0
BiCG: times= 6.489754e-04
BiCG: p_times = 1.597404e-05 (p_c =9.536743e-07 p_i = 1.502037e-05 )BiCG: i_times = 6.330013e-04
BiCG: residual = 7.396952e-15

CGS: iter= 15 iter_double = 15
iter_quad = 0
CGS: times= 6.649494e-04
CGS: p_times = 3.385544e-05 (p_c =9.536743e-07 p_i = 3.290176e-05 )CGS: i_times = 6.310940e-04
CGS: residual = 1.141642e-14

BiCG: rel_resid = 1.067658e-15

CGS: rel_resid = 1.647818e-15$ qstat# ジョブ実行中の表示例

2010/09/15 (Wed) 17:22:58:
REQUESTS on HA8000 cluster
NQS schedule stop time : 2010/09/25 (Sat) 8:55:00 (Remain: 231h 32m 2s)
REQUEST
NAME
OWNER
QUEUE
PRI NODE E-TIME
MEM
STATE
607391.batch1
qsub_mpt00003lecture63
460s
28GB RUNNING$ qstat# ジョブ終了後の表示例

2010/09/15 (Wed) 17:23:24:
REQUESTS on HA8000 cluster
NQS schedule stop time : 2010/09/25 (Sat) 8:55:00 (Remain: 231h 31m 36s)
REQUEST
NAME
OWNER
QUEUE
PRI NODE E-TIME
MEM
STATE
No requests.

$ type lrt
lrt is aliased to `ls -lrt‘$ lrt
・・・ 省略 ・・・-rw------- 1 t00003 t00-rw------- 1 t00003 t00-rw------- 1 t00003 t00-rw------- 1 t00003 t00-rw------- 1 t00003 t00-rw------- 1 t00003 t00

0 Sep 15 17:22 ERR# 標準エラー出力のファイル
3242 Sep 15 17:22 sol1.txt # 解ベクトル1(solution)
208 Sep 15 17:22 res1.txt # 相対残差1(residual)のノルム収束履歴
3242 Sep 15 17:23 sol2.txt # 解ベクトル2(solution)
208 Sep 15 17:23 res2.txt # 相対残差2(residual)のノルム収束履歴
864 Sep 15 17:23 OUT# 標準出力のファイル
35./test1 testmat.mtx 0 sol1.txt res1.txt# 求解1の結果（デフォルト： -i bicg，-p none）./test1 testmat.mtx 0 sol2.txt res2.txt -i 3 -p 1
# 求解2の結果
36

余談: 当実習に役立つユーティリティ１
cdlis： 他版のtestディレクトリ下へ移動する$ ls

morelis： 他版のtestディレクトリ下のファイルの内容を参照するコマンド$HOME/sample/sh/

cdlis_kernel.sh*

余談: 当実習に役立つユーティリティ２$ type morelis

morelis*morelis is hashed (/home/t00003/sample/sh/morelis)$ pwd$ type cdlis/home/t00003/nfs/lis-1.2.49_omp/test

cdlis is aliased to `source $HOME/sample/sh/cdlis_kernel.sh'$ more OUT$ pwd/home/t00003/nfs/lis-1.2.49_omp/test/home/t00003/sample$ cdlis seq# 引数のseqとは下記の赤文字で記したディレクトリ名の一部/home/t00003/nfs/lis-1.2.49_seq/test# seqの部分がSEQ だったら，cdlisの引数はSEQ$ pwd$ morelis seq OUT/home/t00003/nfs/lis-1.2.49_seq/test$ cdlis omp

number of processes = 1・・・/home/t00003/nfs/lis-1.2.49_seq/test# ompも同様である．

/home/t00003/nfs/lis-1.2.49_omp/test

number of processes = 1
・・・

$ pwd
/home/t00003/nfs/lis-1.2.49_omp/test$ pwd/home/t00003/nfs/lis-1.2.49_omp/test

37

38

余談: test1コマンドは他のディレクトリ
test1はダイナミックリンク： でも利用可能$ cdlis$ mkdir$ cp -pseq$HOME/tmp
test1 qsub_seq.nqs testmat.mtx

前半の内容

test1*

次
後半の内容$HOME/tmp$ cd $HOME/tmp$ ls
qsub_seq.nqs

目

testmat.mtx$ $ qsub qsub_seq.nqs
Request 608801.batch1 submitted to queue: lecture.

$ lrt
total 1800
-rw-r--r-- 1 t00003 t00
19290 Oct 31 2007 testmat.mtx-rwxr-xr-x 1 t00003 t00 1669905 Sep 14 16:54 test1*-rw-r--r-- 1 t00003 t00
512 Sep 14 19:24 qsub_seq.nqs-rw------- 1 t00003 t00
0 Sep 18 15:51 ERR-rw------- 1 t00003 t00
844 Sep 18 15:51 OUT-rw-r--r-- 1 t00003 t00
3242 Sep 18 15:51 sol1.txt-rw-r--r-- 1 t00003 t00
208 Sep 18 15:51 res1.txt-rw-r--r-- 1 t00003 t00
3242 Sep 18 15:51 sol2.txt-rw-r--r-- 1 t00003 t00
208 Sep 18 15:51 res2.txt

0. 準備：Lisとテスト問題の入手【4】9. Lisの概観【41】1. HA8000での実習環境設定【5-7】10. 線形方程式用Lisの概要【42】2. Lisインストールの準備作業【8,9】11. 当実習で用いた評価条件【43】3. Lis逐次版のインストール【10-13】

12. 非定常反復法に関する盲点【44-53】

4. 補足事項【14】5. ジョブ実行（逐次版） 【15-17】

13. PETScのインストールと使用要領【54-56】

6. Lis使用に関する解説【18-22】
7. 入力データ形式について【23】8. Lis並列版の準備・インストール【24-39】39【 】：当資料中のスライド番号40

9. Lis (a Library of Iterative Solvers forlinear systems ) の概観‐ 線形方程式求解アルゴリズムの逐次版，並列版（OpenMP, MPI）.
‐ １１種類のデータ形式をサポート. (CRS, CCS, MSR, DIA, …)‐ 倍精度，４倍精度演算サポート‐ Lis ( 小武守，藤井，西田，長谷川 )21 解法 × 11 前処理 × 3 スケーリング ( “前処理無し”，“スケーリング無し” も含める )標準固有値問題

Ax = λx

Ax = b[ http://www.ssisc.org ]Lis (a Library of Iterative Solvers for linear systems)

大規模実疎行列の係数行列 (A)：
線形方程式

10. 線形方程式用Lisの概要

Component of Lis: [Lis‐1.2.49 ver.]これらを解くためのスケーラブルな反復解法ライブラリ
基本的にC言語向けのライブラリだが，Fortranのインタフェース(API)もあり(† Lis の特徴：国内研究者の成果，または，力を入れている)

反 復 解 法
非定常反復解法
CG, BiCG, CGS, BiCGStab, BiCGStab(l ),GPBiCG†, TFQMR, Orthomin(m), GMRES(m)

開発言語はC言語だが，一部SA‐AMG前処理のみFortran90で記述されている．

BiCGSafe†, CR, BiCR†, CRS†, BiCRStab†,
GPBiCR†, BiCRSafe†, FGMRES(m), IDR(s),MINRES

この講習会では線形方程式の方を取り上げて説明します．

定常反復解法41

11. 当実習で用いた評価条件
x0 = 0-initx_zeros true

rk 2 b 2 ≤ 1.0 ×10−12-tol 1.0e-12

スケーリング（広義の前処理の一種）
None,

D −1 Ax = D −1b.

D −1/ 2 AD −1/ 2 D1/ 2 x = D −1/ 2b

Jacobi, Gauss-Seidel, SOR

42

1. ２つの残差ベクトルについて
・求解アルゴリズム中の残差ベクトル・算出された数値解を元の方程式に代入した真の残差ベクトル
2. 残差ベクトルを用いた収束判定について‐ 収束判定条件:

3. 前処理とスケーリングと収束判定との関係‐ ベクトル r ： アルゴリズム中の残差ベクトル ,参考文献：

rk +1 = rk − α k Apk[1] 伊藤祥司，杉原正顯，線形方程式求解に対するクリロフ部分空間法の前処理系に着目した体系的な特性分析，日本応用数理学会2009年度年会講演予稿集，大阪大学，9月(2009).
[2] 伊藤祥司，杉原正顯 ，双ランチョス系統の前処理付きアルゴリズムの改善，計算工学講演会論文集
Vol.15 (2010年5月）pp.171-174．

（CG法の例）[3] 伊藤祥司，杉原正顯，姫野龍太郎，クリロフ部分空間法に対する前処理方式と収束判定について，情報処理学会論文誌コンピューティングシステム(ACS)，Vol.3, No.2, pp.9–19 (June 2010)．

‐ 数値解が有効かどうかの判定 ： 真の残差ベクトルrˆ = b − Axˆ ,Crout ILU, ILUT, (additive schwarz)

12. 非定常反復法に関する盲点

Lisのオプション指定方法（共にデフォルト値）‐ 反復法の初期ベクトル:

前 処 理
None, (Point) Jacobi, ILU(k), SSOR, Hybrid†,
I+S†, SAINV, SA-AMG†, additive Schwarz

x̂ : 数値解[4] 伊藤祥司，杉原正顯 ，双ランチョス系アルゴリズムの前処理変換について，日本応用数理学会
2010年度年会，東京，9月，2010年．

rˆ 2 b 2[5] 伊藤祥司，杉原正顯 ，クリロフ部分空間法に対する前処理方向とライブラリ実装における注意点，
日本応用数理学会2010年度年会，東京，9月，2010年．
43

44

12. 線形方程式求解のための反復解法とその収束性に関する研究の概要
非定常反復法の代表として共役勾配法(CG: Conjugate Gradient method)

非定常反復法では，残差ベクトルに相当する情報も用いて求解アルゴリズムを構成している（収束判定の計算コスト面で有用である）

適当な初期値 x0

p0 = r0 = b − Ax0

12. 当実習で用いた評価条件（再掲）for k = 0,1,L until rk ≤ ε b⎛⎞log10 ⎜⎜ rbk ⎟⎟⎝⎠解へ収束する様子を残差ベクトルのノルムを用いて評価

相対残差ノルム

end‐ 反復法の初期ベクトル:x0 = 0-initx_zeros true

rk 2 b 2 ≤ 1.0 ×10−12-tol 1.0e-12‐ 収束判定条件:

収束判定

α k = ( pk , rk ) ( pk , Apk )xk +1 = xk + α k pk
rk +1 = rk − α k Apk
残差を表すベクトルβ k = − (rk +1 , Apk ) ( pk , Apk )pk +1 = rk +1 + β k pk

Lisのオプション指定方法（共にデフォルト値）‐ ベクトル r ： アルゴリズム中の残差ベクトル ,rk +1 = rk − α k Apk（CG法の例）‐ 数値解が有効かどうかの判定 ： 真の残差ベクトルrˆ = b − Axˆ ,x̂ : 数値解rˆ 2 b 2

反復回数
45

46

12. 前処理付き共役勾配法(PCG法)

12. main関数での真の残差の算出
test1_org.c test1_mod.c との比較：b$ diff test1_org.c test1_mod.c

55a56
2 用の変数を宣言>
LIS_REAL
b_nrm2;// addition
155a157,158
を算出>
lis_vector_nrm2(b, &b_nrm2);// addition2>
177a181>
printf("%s: rel_resid = %e¥n¥n",solvername,resid/b_nrm2); // addition

PCG法のアルゴリズムbrˆ = b − Axˆ ,x̂ : 数値解rˆ 2 b 2 を算出

いるかが収束性向上のポイント/home/t00003/nfs/lis-1.2.49_seq/test

BiCG: residual = 6.399933e-15
BiCG: rel_resid = 9.237507e-16※ただし，実際には上記の変換は行わず，
アルゴリズムで同値な演算を行う．

K が，どれくらい A を近似して$ cat OUT
・・・

不完全コレスキー分解

前処理の演算：r̂ 2

Lisのデフォルト (test1_org.c)rˆ 2 b 2

当講習会で併用 (test1_mod.c)47

q = K-1r (=L-TL-1r ) に対して，
y=L-1r （前進代入），および，
q=L-Ty （後退代入）を計算する
48

12. CGS法と前処理付きCGS法の例
A ≈ K = KLKR

Ax = bこの変換自体は実施しない

x0 , r0 = b − Ax0 , r0# , β −1 = 0,
k = 0,1,2, L; Do

pk = rk + β k −1 z k −1 ,uk = pk + β k −1 ( z k −1 + β k −1uk −1 ),αk =r0# , rk,#0

r , Auk

不完全LU分解K L−1 AK R−1 K R x = K L−1b

αk =

#
0

−1

r , AK uk

End Doαk =D −1b − D −1 Axˆ

,#0r , Auk

r0# , rk,

r0# , Aˆ uk

xk +1 = xk + α k ( pk + z k ),rk +1 = rk − α k AK ( pk + z k ),rk +1 = rk − α k A( pk + z k ),xk +1 = xk + α k ( pk + z k ),r = r − α Aˆ ( p + z ),

−1

βk =

#
0

r , rk +1#
0

r , rkβk =,End Dor0# , rk +1
r0# , rk,

$HOME/sample/test/test1_mod2.cx0 , r0 = K −1 (b − Ax0 ), r0# , β −1 = 0,
k = 0,1,2, L; Do
pk = rk + β k −1 z k −1 ,αk =r0# , rk#0

−1

r , K Aukr , rk +1
r0# , rk

k

k

,
50

rk +1 = rk − α k K −1 A( pk + z k ),βk =End Do#0r , rkK−1(b − Axk )※ 両側・左・右前処理変換後の前処理付きアルゴリズムは一致する（一致性がある）

x0 , r0 = b − Ax0 , r0# , β −1 = 0,

収束判定

k = 0,1,2, L; Do
pk = rk + β k −1 z k −1 ,≤ε

左側は係数行列の部分のみを変換したアルゴリズム.

,

xk +1 = xk + α k ( pk + z k ),r0# , rk +1

K −1 Ax = K −1b

K −1b

z k = pk − α k K −1 Auk ,

51

k

#
0

12. CGS法の２つの異なる左前処理
付きアルゴリズムと収束判定

uk = pk + β k −1 ( z k −1 + β k −1u k −1 ),→ 後者で修正したものが

βk =

k

スケーリングの場合は，係数行列，右辺項とも，
元とは異なるEnd Do

見るからに
左前処理付きアルゴリズム

test1.cでは，元のAとbを保持，または，あらためて読込んで正確に判定する必要がある.

k +1

End Do49となっており，本来求めるべき情報とは異なるので注意．
※ I+S前処理でもスケーリングを用いているので同様に注意．

αk =xk +1 = xk + α k K ( pk + z k ),−1

D −1 Ax = D −1b

この場合の真の残差：

uk = pk + β k −1 ( z k −1 + β k −1uk −1 ),z k = pk − α k Aˆ uk ,−1

12. Lisにおけるスケーリングの注意点

Lisでは，スケーリングすると，
Aとbのデータは上書きされる．

k = 0,1,2, L; Do
pk = rk + β k −1 z k −1 ,z k = pk − α k Auk ,xk +1 = xk + α k ( pk + z k ),,x0 , r0 = bˆ − Aˆ x0 , r , β −1 = 0,

pk = rk + β k −1 z k −1 ,r0# , rkAˆ x = bˆ#0uk = pk + β k −1 ( z k −1 + β k −1uk −1 ),,rk +1 = rk − α k A( pk + z k ),r0# , rk

D = diag ( A) D −1 Ax = D −1b

k = 0,1,2, L; Do

uk = pk + β k −1 ( z k −1 + β k −1uk −1 ),r0# , rk

この変換自体を実施する

x0 , r0 = b − Ax0 , r0# , β −1 = 0,

k = 0,1,2, L; Do
pk = rk + β k −1 z k −1 ,z k = pk − α k AK uk ,βk =Ax = bx0 , r0 = b − Ax0 , r0# , β −1 = 0,z k = pk − α k Auk ,r0# , rk +1

12. CGS法とスケーリングを施した線形方程式に対するCGS法の例

uk = pk + β k −1 ( z k −1 + β k −1u k −1 ),αk =r0# , rk
r0# , AK −1u k,z k = pk − α k AK −1uk ,右側は３種類の前処理方向の変換で，全て同じア
ルゴリズムとなる.しかも，左前処理にも関わらず，収束判定は下式の通り,
収束判定

b − Axkb≤ε

xk +1 = xk + α k K −1 ( pk + z k ),rk +1 = rk − α k AK −1 ( pk + z k ),βk =End Dor0# , rk +1
r0# , rk,52

12. Lisにおける前処理付き
アルゴリズムの状況

13. PETScインストールの準備作業
PETSc利用環境の準備：$ cd $HOME/nfs
$ tar
$ ls

Lisの非定常反復解法 cg(1)～gmres(9)は，右前処理系であり，アルゴリズムの収束判定も正常．

lis-1.2.49/$ mv

ただし，bicgsafe(13)～minres(22)については，伊藤の方では未確認なので，左前処理系にも関わらず収束判定がb − Axk≤εb$HOME/sample/library/petsc-lite-3.1-p3.tar.gz
lis-1.2.49_mpi/

petsc-3.1-p3

lis-1.2.49_omp/

petsc-lite-3.1-p3

lis-1.2.49_seq/

petsc-3.1-p3/# lite版と明示しておく方が良いと思う＃ PETScでは，Intelコンパイラの環境設定してもconfigureでは，gcc，gfortranが設定されるため，ここではGNUコンパイラを用いることにする．
$ cd petsc-lite-3.1-p3
bin/
conf/

config/configure@include/makefile

projects/
src/

TAGS
tutorials/

zope/$ ls -l configure
lrwxrwxrwx 1 t00003 t00 19 Sep 19 13:35 configure -> config/configure.py*と甘い場合があるかも知れないので注意．
※ 定常反復法（jacobi(10)～sor(12)）と非定常反復法とでは論点が異なるので注意．

zxf

53

●PETScのホームページ
http://www.mcs.anl.gov/petsc/petsc-as/●PETScのソース：
http://www.mcs.anl.gov/petsc/petsc-as/download/index.html●PETScのインストール手順：
http://www.mcs.anl.gov/petsc/petsc-as/documentation/installation.html

13. PETScインストール

54

13. PETScの使用要領
サンプルプログラムからLisのtest1のような機能を自作する：

PETScのインストール：$ ./configure --with-cc=gcc --with-fc=gfortran --download-f-blas-lapack=1--download-mpich=1$ pwd/home/t00003/nfs/petsc-lite-3.1-p3$ ls===============================================================================
Configuring PETSc to compile on your system===============================================================================
・・・ 途中省略 ・・・ （１０分位かかる）
xxx=========================================================================xxx
Configure stage complete. Now build PETSc libraries with:make PETSC_DIR=/nfs/all/t00003/petsc-lite-3.1-p3 PETSC_ARCH=linux-gnu-c-debug allxxx=========================================================================xxx

bin/
conf/
config/configure@configure.log@
externalpackages/

include/linux-gnu-c-debug/
makefile

make.log@
projects/
RDict.log

src/
TAGS
tutorials/

zope/$ ls src/mat/examples/tests
ex100.cex113.c・・・ 途中省略 ・・・
ex104.cex117.cex105f.F ex118.c
ex106.cex119.cex107.cex11.c・・・ 以下省略 ・・・$ make all test

● これら２種類の環境変数が重要・・・ 途中省略 ・・・ （１０分位かかる）
Completed building libraries=========================================
Now to check if the libraries are working do:make PETSC_DIR=/nfs/all/t00003/petsc-lite-3.1-p3 PETSC_ARCH=linux-gnu-c-debug test=========================================
Running test examples to verify correct installation
C/C++ example src/snes/examples/tutorials/ex19 run successfully with 1 MPI process
C/C++ example src/snes/examples/tutorials/ex19 run successfully with 2 MPI processes
Fortran example src/snes/examples/tutorials/ex5f run successfully with 1 MPI process
Completed test examples

ex126f.F

ex21.cex37.cex51.cex63f.F

ex77.cex95.cex12.cex13.cex14.cex15.cex27.cex28.cex2.cex30.cex41.cex42.cex43.cex44.cex55.cex56.cex57.cex58.cex67f.F
ex68.cex6.cex70.cex80.cex81.cex85f.F
ex86.cex9.c
makefile
output/＃ 例えば，Matrix Market形式の行列を入力するプログラム例は ex72.c である．
＃ この「tests」下の ex72.c, makefile をPETSC_DIR下にコピーして make する．
＃ ただし，元々の「makefile」を上書きしないよう，makefile2 とコピーした例で以下説明● その際のスクリプト例：#!/bin/bash
export PETSC_DIR="/home/t00003/nfs/petsc-lite-3.1-p3"
export PETSC_ARCH="linux-gnu-c-debug"
55

make –f makefile2 ex72

● PETScの詳細についてはマニュアル参照56