当講習会の目標「ライブラリ利用：
高性能プログラミング初級入門」
PartⅡ

１．６０分後には，線形方程式に対して，１００種類
近くのアルゴリズムで求解できるようになる．

2011年12月22日（木）３．１８０分後には，反復解法に基づく求解アルゴ
リズムの特徴と注意点を理解している．

２．９０分後には，線形方程式に対して，並列計算による求解ができるようになる．

東京大学 情報基盤センター
伊藤祥司
itosho@cc.u‐tokyo.ac.jp

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

Lisだ！
2

目
前半の内容

次

0. 準備：Lisとテスト問題他の入手◆ Lis [ http://www.ssisc.org/ lis/ ]後半の内容

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】
13. test1 処理フローの概要【54,55】

4. 補足事項【14】5. ジョブ実行（逐次版） 【15-17】
6. Lis使用に関する解説【18-22】

14. PETScのインストールと使用要領【56-58】7. 入力データ形式について【23】8. Lis並列版の準備・インストール【24-39】【 】：当資料中のスライド番号3→ Downloadの「Version 1.2.xx」からソースコード，ユーザマニュアル
をダウンロード※ ソースコードは次画面の「ダウンロードページ」→ 「Software License Agreement」へと進む．
ユーザ登録する／しないは各自の自由※ 当資料，当講習会では，lis‐1.2.49 版を前提として説明，実習します．
◆ テスト問題 Matrix Market [ http://math.nist.gov/MatrixMarket ]→ BrowseやSearchの項目のリンクをクリック◆ SESNA : Survey and Evaluation System for Numerical Algorithms [ http://sesna.jp ]求解アルゴリズムの体系的性能評価システム→ Matrix Market の一部問題を Lisの一部を用いて求解した
結果情報のDBと求解結果の評価システム
4

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
add20.mtx
add20_rhs1.mtx

sh/

test/

library

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

add32_rhs2.mtx
add32.rua
fs_183_1.mtx
fs_183_1.rua

setup.env# JOB実行用スクリプト
qsub_seq.nqs# 今回の実習で役立つユーティリティ# Lis実行で使用するmainルーチン

test1_mod.ctest# 実習環境のセットアップ# このような環境を設定$ 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

5

1. HA8000での実習環境設定(3/3)
$ 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/{各自の利用者番号}

nfs

6

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

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

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

morelis*＃ 以後，{各自の利用者番号}を 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

qsub_script$ source setup.env# テスト用入力データ（係数行列，右辺項データ）
memplus.mtx
nos3.rsa
memplus_rhs1.mtx
memplus.rua
nos3.mtx● コマンド実行する部分は太字で表示してます（左端の$はUNIXのコマンドプロンプト）．
● 画面表示は，多少，文字サイズを小さくしてます（当テキスト編集の都合上）．
● 特に注意の無い限り，# 記号から右部分の青文字はコメントです．

$ ln –smatpetsc-lite-3.1-p3.tar.gz
petsc-manual-3.1.pdf

sample/qsub_script:
qsub_mpi.nqs qsub_omp.nqs

sample/test:
test1_mod2.c$ cd sample$ ls# 今回の実習で使用するライブラリ等

add20.rua
add32.mtx
add32_rhs1.dat
add32_rhs1.mtx

sample/sh:
cdlis_kernel.sh*.

setup.env

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

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

mpich-mx-intel11/opt/intel/Compiler/11.0/074/bin/intel64/icc/opt/intel/Compiler/11.0/074/bin/intel64/ifort# 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$HOME/nfs/lis-1.2.49_seq/ 下：
include/config.h，lis.h，lisf.h
src/.libs/
liblis.a# main関数の入替え

rm: remove symbolic link `test1.c'? y
rm: remove regular file `test1.o'? y$ 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 1#@$-lm 2GB#@$-lT 0:01:00#@$-o OUT#@$-e ERR#@$-s /bin/bash# 標準出力のファイル名# 標準エラー出力のファイル名# スクリプト内のシェルを指定

cd $PBS_O_WORKDIR
pwd# スクリプト内もカレントディレクトリへ移動する# pwdは本来は不要だが，確認のため講習会後のキュー名．本日はtutorial (sampleは変更済み)

5. HA8000でのジョブ実行(逐次版)$ qsub qsub_seq.nqs● スクリプトの拡張子は任意

Request 606870.batch1 submitted to queue: lecture.

リクエストID（ジョブ固有の番号）$ qstat# ジョブ実行中の表示例● NQSのスクリプトでは左端の“ #@$ ”記号は，これに続くハイフンで始まる文字列がNQSオプション
であることを宣言している．

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● 赤文字の項目をチェック（NQS
オプション“-N”でノード数を指定している）

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.

./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$ qstat# ジョブ終了後の表示例$ 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どちらでも良い[ ]内はデフォルト値

CG‐i {cg | 1}なし

‐p {none | 0}

BiCG

‐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] 線型方程式系解法，他‐ell [2]解法の l の値

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]

6. 色々なデータで実行してみよう
$ ls
$ ln
$ ls

‐storage_block [2]BSR，BSCのブロックサイズ● 上記のオプションの内容の詳細はマニュアル参照．

-s-lqsub_script/

setup.env

sh/

test/

$HOME/sample/mat .
mat$ ls mat/494_bus.mtx
494_bus.rsa
add20.mtx
add20_rhs1.mtx

add20.rua
add32.mtx
add32_rhs1.dat
add32_rhs1.mtx

add32_rhs2.mtx
add32.rua
fs_183_1.mtx
fs_183_1.rua

memplus.mtx
memplus_rhs1.mtx
memplus.rua
nos3.mtx

参考までに，matの後にディレクトリを示す「/」を付け忘れると：$ ls mat# これは，☆１のsetup.env 影響‐omp_num_threads [t] 実行スレッド数．t は最大スレッド数
行列格納形式

mat/

lrwxrwxrwx 1 t00003 t00 23 Sep 14 19:39 mat -> /home/t00003/sample/mat/

初期ベクトル x0 の設定（デフォルトは全要素の値が0）‐storage [0]$HOME/sample/

library/mat@$ cat qsub_seq.nqs

nos3.rsa● 拡張子
mtx : Matrix Market形式
rua, rsa : Harwell-Boeing形式
dat : ベクトルのPLAIN形式(拡張子自体は任意)#./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形式をサポートしたりツール提供されているので，初級者や中級者にとっては都合が良いと思われる．

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● Lisのtest1.cでは，MM形式のデータを読込み，CRS形式に変換する．MPIのデータ分散にも対応している．

$ rm test1.c ; ln -s test1_mod.c test1.c# 逐次版同様に test1.c を入替える

rm: remove symbolic link `test1.c'? y● 対称行列の場合には，対角要素から下三角の要素のみが格納されている．この場合，ヘッダには，半角スペースに続いて「symmetric」の文字列が記載されている．
→ mat/下の「494_bus.mtx」「nos3.mtx」が該当$ 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/下の「add{20,32}.mtx」「fs_183_1.mtx」 「memplus.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(HB)形式の係数行列データ# ついでにmatもシンボリックリンク# LisのOpenMP(omp)版の環境準備# LisのMPI(mpi)版の環境準備$ which icc ifort

MMサイトで提供されている add32_rhs1.mtx のヘッダとフォーマットでは，Lisで読み取れないので，
当実習では以下の２通りの修正データを用意してます．問題点と修正内容は，各自で確認して下さい．
add20, memplus等も同様ですが，サンプルではadd32のみ修正データを用意しています．
add32_rhs1.dat : add32_rhs1.mtxを元にPLAIN形式にしたファイル
add32_rhs2.mtx : Lisで読み取り可能なベクトル用拡張MM形式のデータ# lis-1.2.49_seq ではないことに注意/opt/intel/Compiler/11.0/074/bin/intel64/icc/opt/intel/Compiler/11.0/074/bin/intel64/ifort# Intelコンパイラのiccのパス# Intelコンパイラのifortのパス＃ Intelコンパイラの環境設定されていない場合には：$ source /opt/itc/mpi/mpiswitch.sh mpich-mx-intel11
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)
Lis並列(omp)版のインストールチェック(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/

・・・ 出力の途中部分 ・・・
make check-TESTS
make[2]: Entering directory `/nfs/all/t00003/lis-1.2.49_omp/test'=== Running test test.shqsub_mpi.nqs

qsub_omp.nqs

赤文字の項目をチェック
下線の項目は計測時間なので，一致して無くても良い#@$-q lecture#@$-N 1#@$-lm 2GB#@$-lT 0:01:00#@$-o OUT#@$-e ERR#@$-s /bin/bash● スクリプトの拡張子は任意

qsub_seq.nqs$ cp –p $HOME/sample/qsub_script/qsub_omp.nqs .
$ cat qsub_omp.nqs

checking solvers...

# 標準出力のファイル名# 標準エラー出力のファイル名# スクリプト内のシェルを指定

export OMP_NUM_THREADS=4 # 環境変数によるスレッド数の指定
cd $PBS_O_WORKDIR# スクリプト内もカレントディレクトリへ移動
pwd# pwdは本来は不要だが，確認のため● 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

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の一種）でジョブ実行するスクリプトを使用する：$ make checknumber 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)
Lis並列(mpi)版のインストールチェック(make check)：

HA8000では，PBS（NQSの一種）でジョブ実行するスクリプトを使用する：$ make check$ cd test$ ls $HOME/sample/qsub_script/

・・・ 出力の最後部分 ・・・
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

qsub_omp.nqs● スクリプトの拡張子は任意

qsub_seq.nqs$ cp –p $HOME/sample/qsub_script/qsub_mpi.nqs . ● NQSのスクリプトでは左端の“ #@$ ”記号は，これに続くハイフ$ cat 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で実行するスクリプトに「make check」を記述すれば，インストールチェックできる．

8. ジョブ実行スクリプトの用意(MPI版)

33

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

ンで始まる文字列がNQSオプショ
ンであることを宣言している．

#@$-q lecture#@$-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)

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
208 Sep 15 17:22 res1.txt
3242 Sep 15 17:23 sol2.txt
208 Sep 15 17:23 res2.txt
864 Sep 15 17:23 OUT# 標準エラー出力のファイル# 解ベクトル1(solution)# 相対残差1(residual)のノルム収束履歴# 解ベクトル2(solution)# 相対残差2(residual)のノルム収束履歴# 標準出力のファイル
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

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

test1*

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

39

38

目

9. Lis (a Library of Iterative Solvers forlinear systems ) の概観

次

前半の内容

後半の内容

大規模実疎行列の係数行列 (A)：

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】
13. test1 処理フローの概要【54,55】

Ax  b

4. 補足事項【14】5. ジョブ実行（逐次版） 【15-17】
6. Lis使用に関する解説【18-22】

14. PETScのインストールと使用要領【56-58】7. 入力データ形式について【23】8. Lis並列版の準備・インストール【24-39】標準固有値問題Ax  x

これらを解くためのスケーラブルな反復解法ライブラリ
基本的にC言語向けのライブラリだが，Fortranのインタフェース(API)もあり
開発言語はC言語だが，一部SA‐AMG前処理のみFortran90で記述されている．

【 】：当資料中のスライド番号
この講習会では線形方程式の方を取り上げて説明します．
40

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

41

11. 当実習で用いた評価条件[ http://www.ssisc.org ]Lisのオプション指定方法（共にデフォルト値）

Lis (a Library of Iterative Solvers for linear systems)‐ 反復法の初期ベクトル:‐ 線形方程式求解アルゴリズムの逐次版，並列版（OpenMP, MPI）.
‐ １１種類のデータ形式をサポート. (CRS, CCS, MSR, DIA, …)‐ 倍精度，４倍精度演算サポート‐ Lis ( 小武守，藤井，西田，長谷川 )22 解法 × 11 前処理 × 3 スケーリング ( “前処理無し”，“スケーリング無し” も含める )x0  0-initx_zeros true

rk 2 b 2  1.0 1012-tol 1.0e-12‐ 収束判定条件:‐ ベクトル r ： アルゴリズム中の残差ベクトル ,Component of Lis: [Lis‐1.2.49 ver.]反 復 解 法
非定常反復解法
CG, BiCG, CGS, BiCGStab, BiCGStab(l ),GPBiCG†, TFQMR, Orthomin(m), GMRES(m)BiCGSafe†, CR, BiCR†, CRS†, BiCRStab†,
GPBiCR†, BiCRSafe†, FGMRES(m), IDR(s),MINRES
定常反復解法Jacobi, Gauss-Seidel, SOR

rk 1  rk   k Apk(† Lis の特徴：国内研究者の成果，または，力を入れている)
前 処 理（CG法の例）

None, (Point) Jacobi, ILU(k), SSOR, Hybrid†,
I+S†, SAINV, SA-AMG†,
Crout ILU, ILUT, additive Schwarz‐ 数値解が有効かどうかの判定 ： 真の残差ベクトルrˆ  b  Axˆ ,スケーリング（広義の前処理の一種）
1

D Ax  D b.

None,

D

1

1 / 2x̂ : 数値解rˆ 2 b 2

AD 1/ 2 D1/ 2 x  D 1/ 2b
42

43

12. 線形方程式求解のための反復解法とその収束性に関する評価の概要

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

非定常反復法の代表として共役勾配法(CG: Conjugate Gradient method)

1. ２つの残差ベクトルについて
・求解アルゴリズム中の残差ベクトル・算出された数値解を元の方程式に代入した
真の残差ベクトル

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

適当な初期値 x0

p0  r0  b  Ax0

収束判定

for k  0,1, until rk   b

2. 残差ベクトルを用いた収束判定について
3. 前処理とスケーリングと収束判定との関係
参考文献：[1] 伊藤祥司，杉原正顯，姫野龍太郎，クリロフ部分空間法に対する前処理方式と収束判定について，情報処理学会論文誌コンピューティングシステム(ACS)，Vol.3, No.2, pp.9–19 (June 2010)．
[2] 伊藤祥司，杉原正顯 ，クリロフ部分空間法に対する前処理方向とライブラリ実装における注意点，
日本応用数理学会2010年度年会，東京，9月，2010年．

end

解へ収束する様子を残差ベクトルのノルムを用いて評価

相対残差ノルム 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


log10  rbk 



反復回数

44

12. 当実習で用いた評価条件（再掲）x0  0-initx_zeros true

rk 2 b 2  1.0 1012-tol 1.0e-12

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

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

45

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); // additionb‐ 収束判定条件:rˆ  b  Axˆ ,‐ ベクトル r ： アルゴリズム中の残差ベクトル ,rk 1  rk   k Apk（CG法の例）x̂ : 数値解rˆ 2 b 2 を算出$ cat OUT
/home/t00003/nfs/lis-1.2.49_seq/test‐ 数値解が有効かどうかの判定 ： 真の残差ベクトルrˆ  b  Axˆ ,b$ diff test1_org.c test1_mod.cx̂ : 数値解

・・・
BiCG: residual = 6.399933e-15

rˆ 2 b 2

BiCG: rel_resid = 9.237507e-16
46

r̂ 2

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

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

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

12. CGS法と前処理付きCGS法の例
A  K  KLKR

Ax  bPCG法のアルゴリズム
不完全コレスキー分解

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

k  0,1,2, ; Do
pk  rk   k 1 z k 1 ,k  0,1,2, ; Do
pk  rk   k 1 z k 1 ,uk  pk   k 1 ( z k 1   k 1uk 1 ),uk  pk   k 1 ( z k 1   k 1uk 1 ),#0

r , rkK L1 AK R1 K R x  K L1bx0 , r0  b  Ax0 , r0# ,  1  0,x0 , r0  b  Ax0 , r0# ,  1  0,

※ただし，実際には上記の変換は行わず，
アルゴリズムで同値な演算を行う．

不完全LU分解k 

r0# , rk,

K が，どれくらい A を近似して

k 

いるかが収束性向上のポイント

z k  pk   k Auk ,z k  pk   k AK 1uk ,xk 1  xk   k ( pk  z k ),xk 1  xk   k K 1 ( pk  z k ),rk 1  rk   k A( pk  z k ),rk 1  rk   k AK 1 ( pk  z k ),前処理の演算：q = K-1r (=L-TL-1r ) に対して，

k 

y=L-1r （前進代入），および，
q=L-Ty （後退代入）を計算する
48

12. CGS法とスケーリングを施した線形方程式に対するCGS法の例
Ax  bこの変換自体を実施する

x0 , r0  b  Ax0 , r0# ,  1  0,
k  0,1,2, ; Do
pk  rk   k 1 z k 1 ,uk  pk   k 1 ( z k 1   k 1uk 1 ),k 

r0# , rk,

r0# , AukD  diag ( A) D 1 Ax  D 1br0# , Aˆ uk

z k  pk   k Auk ,z k  pk   k Aˆ uk ,xk 1  xk   k ( pk  z k ),xk 1  xk   k ( pk  z k ),r  r   Aˆ ( p  z ),rk 1  rk   k A( pk  z k ),

k 
End Dor0# , rk 1
r0# , rk

,

k 1

k 
End Dokk

r0# , rk 1
r0# , rkk

r , rk 1#0

r , rk,k 

End Dor0# , AK 1uk

r0# , rk 1
r0# , rk,

End Do49● 対角スケーリング：D 1 Ax  D 1b

Aˆ x  bˆ

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

uk  pk   k 1 ( z k 1   k 1uk 1 ),,#0● Jacobiスケーリング：k  0,1,2, ; Do
pk  rk   k 1 z k 1 ,r0# , rk

r0# , Auk

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

x0 , r0  bˆ  Aˆ x0 , r0# ,  1  0,

k 

,

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

例えば，Jacobiスケーリングの
ときの真の残差ベクトル：D 1b  D 1 Axˆ
となっており，本来求めるべき情報とは異なるので注意．
※ I+S前処理でもスケーリングを用いているので同様に注意．

D 1/ 2 AD 1/ 2 D1/ 2 x  D 1/ 2b
test1.cでは，元のAとbを保持，または，あらためて読込んで正確に判定する必要がある.
→ 後者で修正したものが$HOME/sample/test/test1_mod2.c「rel_resid2=」の出力が追加される．

51

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

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

収束判定

rk

x0 , r0  K 1 (b  Ax0 ), r0# ,  1  0,× b

k  0,1,2, ; Do
pk  rk   k 1 z k 1 ,k 

r , rk#0

1

r , K Auk

,

,

K

1



 b  Axk 
1

K b



x0 , r0  b  Ax0 , r0# ,  1  0,
k  0,1,2, ; Do
pk  rk   k 1 z k 1 ,uk  pk   k 1 ( z k 1   k 1u k 1 ),k 

左側は係数行列の部分のみを変換したアルゴリズム.
右側は３種類の前処理方向の変換で，全て同じア
ルゴリズムとなる.しかも，左前処理にも関わらず，収束判定は下式の通り

rk 1  rk   k K 1 A( pk  z k ),r0# , rkb

1

xk 1  xk   k ( pk  z k ),k K 1  b  Axk 

rkz k  pk   k K 1 Auk ,r0# , rk 1

△ K b

uk  pk   k 1 ( z k 1   k 1u k 1 ),#0※ 両側・左・右前処理変換後の前処理付きアルゴリズムは一致する（一致性がある）

K 1 Ax  K 1b

収束判定

End Dork○ b b  Axkb

r0# , rk

Lisの非定常反復解法 cg(1)～gmres(9)は，右前処理系であり，アルゴリズムの収束判定も正常．
ただし，bicgsafe(13)～minres(22)については，伊藤の方では未確認なので，左前処理系にも関わらず収束判定が,rk

z k  pk   k AK 1u k ,b

r0# , AK 1uk

xk 1  xk   k K 1 ( pk  z k ),と甘い場合があるかも知れないので注意．

rk 1  rk   k AK 1 ( pk  z k ),

k 

r0# , rk 1
r0# , rk,

End Do※ 定常反復法（jacobi(10)～sor(12)）と非定常反復法とでは論点が異なるので注意．
52

53

13. test1 処理フローの概要(1)13. test1 処理フローの概要(2)

int main(int argc, char* argv[] ) {lis_initialize
lis_matrix_create
lis_vector_create
lis_input(A,b,x,argv[1]);係数行列 A をmalloc
ベクトル b, x をmalloc
A(MMとHB形式), b(ファイル入力の場合),
x(初期値)を読み込み

lis_matrix_set_type
lis_matrix_convert
右辺項 b の作成処理

行列データをCRSへ変換〃
MMとHB以外の

lis_solve(A,b,x,solver);

求解ルーチン

lis_solver_get_timeex
lis_solver_get_residualnorm

タイマーの値をget
真の残差ノルム算出

lis_output_vector(x,LIS_FMT_MM,argv[3]);lis_solver_output_rhistory(solver, argv[4]);解を出力
残差履歴を出力

int lis_solve( ) {lis_precon_create
lis_solve_kernel(A, b, x, solver, precon);}

前処理演算の設定
求解のための事前処理

int lis_solve_kernel( ) {test1の引数をチェック
x0を設定

形式は，マニュア
ル5節「行列格納形式」のサンプル
を参照

if(I+S前処理選択) スケーリング処理
else (スケーリング選択) スケーリング処理
lis_solver_executey=D^{-1/2}x

I+S前処理は要スケーリング
CG法ではSYMM_DIAG※1 scalingでは A, b のデータを上書き求解（例えば，lis_gmresを呼び出す）「SYMM_DIAG(対角scaling)」&&「I+S前処理以外」のとき実行

lis_vector_nrm2(t,&nrm2);

lis_finalize();
return 0;}

真の残差ベクトルの2ノルム（絶対残差）※2 scaling後の系に対して算出(cf. ※1)

}
54

55

14. PETScインストールの準備作業
PETSc利用環境の準備：

PETScのインストール：$ cd $HOME/nfs
$ tar
$ ls

zxf

lis-1.2.49/$ mv$ ./configure --with-cc=gcc --with-fc=gfortran --download-f-blas-lapack=1--download-mpich=1$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/===============================================================================
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

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

$ make all test

$ cd petsc-lite-3.1-p3
bin/
conf/

config/configure@include/makefile

projects/
src/

TAGS
tutorials/● これら２種類の環境変数が重要・・・ 途中省略 ・・・ （１０分位かかる）
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

zope/$ ls -l configure
lrwxrwxrwx 1 t00003 t00 19 Sep 19 13:35 configure -> config/configure.py*●PETScのホームページ
http://www.mcs.anl.gov/petsc/●PETScのソース：
http://www.mcs.anl.gov/petsc/download/index.html●PETScのインストール手順：
http://www.mcs.anl.gov/petsc/documentation/installation.html

56

14. PETScの使用要領
サンプルプログラムからLisのtest1のような機能を自作する：$ pwd/home/t00003/nfs/petsc-lite-3.1-p3

$ ls
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・・・ 以下省略 ・・・

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"
make –f makefile2 ex72

14. PETScインストール● PETScの詳細についてはマニュアル参照58

57