当講習会の目標
１．６０分後には，線形方程式（Ax=b）に対して，１００
種類近くのアルゴリズムで求解できるようになる．

「ライブラリ利用：
高性能プログラミング初級入門」
PartⅡ

２．９０分後には，線形方程式に対して，並列計算による求解ができるようになる．
３．１８０分後には，Oakleaf‐FX(FX)上で２種類の線形方程式求解ライブラリを使えるようになっている上に，反復解法に基づく求解アルゴリズムの特徴と注意点を理解している．
．
こんな夢のようなことを実現するのが当講習会！

2013年3月12日（火）東京大学 情報基盤センター
伊藤祥司

LisとPETScだ！

その鍵となるのが・・・

2

当講習会の方針

当講習会で扱うLisとPETScの比較
赤文字：ライブラリ選択の観点

１．東京大学情報基盤センターのOakleaf‐FXにおいて，
疎行列向き線形方程式求解ライブラリを使用できるようにする（マシンの利用方法，システム環境の違い
や，ライブラリのバージョンの違いなどのため，他のマシンでも同様の手順とは限らないことに注意）．
２．Lisの典型的な使用方法・習得要領について説明するが，ライブラリ記述の文法面等は取り上げない（基本的に詳細についてはマニュアルを参照）．
３．PETScの典型的な習得要領について説明するが，ライブラリ記述の文法面等は取り上げない（基本的に詳細についてはマニュアルを参照）．

実
習

Lis‐1.3.0

PETSc‐3.3‐p4○線形方程式用
反復解法

22種類

○

前処理演算

10種類

10種類程度*

○

演算

実数演算

実数演算，複素数演算

○

実行形態

逐次演算，MPI，OpenMP，（逐次演算，）MPI
Hybrid（MPI+OpenMP）

10種類程度**User Manualの記述

ユーザプログラム言語 C言語，Fortran

C言語，C++，Fortran

その他の特長様々なProfiling機能，
簡単なGraphics機能

4倍精度演算線形方程式反復解法以外の解法，分野
固有値解法
3線形方程式用直接解法，
非線形解法，ODE

PETScライクな固有値解法ライブラリとしては，SLEPc, the Scalable Library for Eigenvalue Problem
Computations [ http://www.grycap.upv.es/slepc/ ] があるが，PETScの開発グループとは異なる．

4

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

講習内容◆ Lis [ http://www.ssisc.org/ lis/ ]前半の内容→ Downloadの「Version 1.3.xx」からソースコード，ユーザマニュアル
をダウンロード※ 当資料，当講習会では，lis‐1.3.0 版を前提として説明，実習します．
[ http://www.ssisc.org/lis/dl/lis‐1.3.0.tar.gz ]後半の内容

0. 準備，FXでの実習環境設定

9. PETScの概要

1. Lisの概要

10. FXでのPETSc利用の準備

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

11. FXでのPETSc-逐次版利用例

3. Lis逐次版のインストール

12. FXでのPETSc-MPI版利用例

4. ジョブ実行（逐次版）

13. PETSc-複素数演算例題

5. Lis使用に関する解説

14. PETSc使用に関する解説

6. 入力データ形式について

15. PETScテストプログラムの概要

7. Lis並列版のインストール・実行

16. KSPに関する盲点

8. Lisテストプログラムの概要

17. 両ライブラリの習得要領◆ PETSc [ http://http://www.mcs.anl.gov/petsc/ ]→ Downloadのページからソースコード，ユーザマニュアルをダウンロード※ 当資料，当講習会では，Oakleaf‐FXにインストール済みの petsc‐3.3‐p4 版
を前提として説明，実習します（3.0.0版等とは非互換があるので注意）．
[ http://ftp.mcs.anl.gov/pub/petsc/release‐snapshots/petsc‐3.3‐p4.tar.gz ]◆ テスト問題 Matrix Market [ http://math.nist.gov/MatrixMarket ]→ BrowseやSearchの項目のリンクをクリック◆ SESNA : Survey and Evaluation System for Numerical Algorithms [ http://sesna.jp ]求解アルゴリズムの体系的性能評価システム→ Matrix Market の一部問題を Lisの一部を用いて求解した結果情報のDBと求解結果の評価システム

560. 準備：Oakleaf‐FX (FX)のマシン環境

0. 準備：Oakleaf‐FX (FX)のコンパイラ◆ FXのシステム環境と当資料中の記述との対応：◆ FXのコンパイラについて（当資料に関係するもののみ記載）：

inter$ ：インタラクティブノードにおけるコマンドプロンプト（実際は，[${USER}@a10tio297 ${PWD}]$ など）$ （または login$）：ログインノードにおけるコマンドプロンプト（実際は，[${USER}@oakleaf-fx-2 ${PWD}]$ など）
ログインノードのコマンドプロンプトは，
基本的に $ と記述していますが，
インタラクティブジョブ実行する際など，特に区別するときは login$ と記述しています．

今回の講習会ではインタラク
ティブノードは使用不可

MPI
ｺｰﾄﾞ

◆ ジョブ実行用ノード
使用時の注意点：● ジョブ実行用ノード共通
ログインノードにて設定した環境（alias等）は，ジョブ実行用ノードには継承されない（あらためて設定する）．

クロス
コンパイラ

オウン
コンパイラ

Fortran
77/90

frtpx

frt

C言語

fccpx

fcc

C++

FCCpx

FCC

Fortran77/90

mpifrtpx

mpifrtC言語

mpifccpx

mpifccC++

mpiFCCpx

mpiFCC

クロスコンパイラ

オウンコンパイラ

クロスコンパイラ● インタラクティブノードへのログイン
カレントディレクトリは継承される．
● 計算ノードへのバッチジョブ投入
ジョブスクリプトのみが計算ノードへログインするイメージ．Oakleaf-FXのジョブ管理システムではカレント
ディレクトリへログインする（便利な仕様）．

非
MPI
ｺｰﾄﾞ

言語
処理系

「FX10 スーパーコンピュータシステム
Oakleaf-FX 利用手引書 (ver. 2.4)」から

当資料中では「ジョブ実行用ノード」と総称
7

ログインノードにてインタラクティブノード・
計算ノード用の実行モジュールを作成する【通常はこちらを用いれば良い】
当資料の中では，Lisのインストール（サンプルプログラム作成含む），および，PETScのサンプルプログラム作成で用いている．

オウンコンパイラ

ログインノードの負荷が高いときや，コンパイルに相当な時間がかかるときなどに有効（コンパ
イラ自体も実行モジュールである）「FX10 スーパーコンピュータシステム
Oakleaf-FX 利用手引書 (ver. 2.4)」から

当資料中では「ジョブ実行用ノード」と総称
8

0. FXでの実習環境設定(1/5)

0. FXでの実習環境設定(2/5)

各自の$HOME下に実習用sampleファイル一式をコピー：

sample:library/mat/script/sample/library:lis-1.3.0.tar.gz

setup.env

manual/sample/library/manual:lis-manual-1.3.0_ja.pdf$ ls sample
library

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

sh/petsc-3.3-p4_samp/mat

scriptsetup.envsh# 実習環境のセットアップ$ source sample/setup.env

ログインの度に実行する（☆１）# 今回の実習で使用するライブラリ等# 今回の実習で使用するライブラリに対応した版のマニュアル

petsc-manual-3.3-p4.pdf

sample/library/petsc-3.3-p4_samp:
ex1.c ex11.c@ ex11_mod.c ex11_org.c
sample/mat:494_bus.mtx
494_bus.rsa
add20.mtx
add20.rua

各自の実習用環境をセットアップ：# ホームディレクトリ($HOME)へ移動# 実習用sampleファイル一式を入手・展開$ cd[Enter]キー$ tar xf /home/t00003/sample.tar$ ls -FR sample# 今回の実習で使用するPETScの例題プログラム等
ex1_mod.c ex23.c makefile

add20_rhs1.mtx
add32.mtx
add32.rua
add32_rhs1.dat

add32_rhs1.mtx
add32_rhs2.mtx
fs_183_1.mtx
fs_183_1.rua

sample/script:
pjsub_hyb.sh* pjsub_mpi.sh*# テスト用入力データ（係数行列，右辺項データ）
memplus.mtx
nos3.rsa
memplus.rua
memplus_rhs1.mtx
nos3.mtx# pjsub用スクリプト

pjsub_omp.sh*pjsub_seq.sh*# このような環境を設定# .bashrc にて指定することも可能$ cat sample/setup.env
alias ls=‘ls -F’alias lrt='ls -lrt'
alias cp='cp -i'
alias mv='mv -i'
alias rm='rm -i‘
PATH=$PATH:$HOME/sample/sh# 講習会向けJOB実行用ツール

sample/sh:isub*● プロンプトと入力コマンドの部分は太字で表示してます．
● 画面表示は，多少，文字サイズを小さくしてます（当テキスト編集の都合上）．
● 特に注意の無い限り，（# 記号から右部分の）青文字はコメントです．

9

0. FXでの実習環境設定(3/5)
当講習会でのジョブ実行方法について（当講習会固有のコマンド）：＃ これは ☆１のsetup.env を実行した上での話$ type isub

isub is /home/t00003/sample/sh/isub

0. FXでの実習環境設定(4/5)※ isubツールは使用者側の責任で使用して下さい．

＃擬似的なインタラクティブ実行にてバッチジョブを投入するコマンド

isub : Interact mode-like SUBmitting tool

両側をダブルクォーテーション

で囲み，コマンド列を記述
Usage:$ isub {-seq|-omp|-mpi|-hyb} [–th n] "command-line" [pjsub-option]
option:-seq : 1 node
1 thread job-omp : 1 node
8 threads job (default)
isubコマンドを実行したディレクトリ上-mpi : 4 nodes 1 thread job (default)にOUTというファイルができる．
-hyb : 4 nodes 8 threads job (default)
すでに同名のファイルが存在してい-th:

当講習会でのジョブ実行方法について（一般的なバッチジョブ方式）：$ ls sample/script
pjsub_hyb.sh*pjsub_mpi.sh*

pjsub_omp.sh*pjsub_seq.sh*$ more sample/script/pjsub_seq.sh＃ここで 「 -bash: type: isub: not found 」 などと表示されたら，「0.FXでの実習環境設定(2/5)」の☆１を実行$ isub10

Change the number of threads
る場合には上書きされるので注意．
n : The number of threads (default 8) ただし，pjsub‐optionで回避可能．

#!/bin/sh#PJM -L "rscgrp=tutorial"#PJM -L "node=1"#PJM -L "elapse=10:00"#PJM -o OUT#PJM -j###### For Lis examples###

上記の全スクリプトにてデフォルトは make checkである．
実習内容に応じ，スクリプト内のコメント記号を外したり，各自でコマンドを書き加える．

make check# ./test1 testmat.mtx 0# ./test1 testmat.mtx 0

このスクリプトにてpjsub実行したディレクトリ上にOUTというファイルができる．
すでに同名のファイルが存在している場合には上書きされるので注意．

sol1.txt
sol2.txt

rsd1.txt
rsd2.txt -i 3 -p 1

・・・ 以下省略 ・・・

output file (default):
OUT : STDOUT + STDERR
Created by S.Itoh 20121218＃上記のスクリプトを実行モジュールがあるディレクトリにコピーし，必要に応じスクリプトを書き換えた上で：$ pjsub pjsub_seq.sh※ -seq オプションは，あくまでも“1ノード，1スレッドジョブ用” という意味です（非MPIジョブ用という意味では無い）．

11

いま，この時点では，これは実行しないとすれば，逐次（seq）実行のバッチジョブが投入される．

12

0. FXでの実習環境設定(5/5)

1. Lisの概要（線形方程式版）[ http://www.ssisc.org/lis ]当講習会でのジョブ実行用ツール・スクリプトの設定内容：

Lis (a Library of Iterative Solvers for linear systems)

・講習会中のキュー名： tutorial（実習用sampleは変更済み）
・講習会後のキュー名： lecture（ただし，当資料作成ではこちらのキューを使用）$ more $HOME/sample/sh/isub#!/bin/sh
MNODES=4
MTHREADS=8#QUEUE="lecture"
QUEUE="tutorial"
OUTPUTFILE="OUT"● isubコマンドのデフォルト設定：
ノード数 = 4 （mpi，hybの場合）
スレッド数 = 8 （ompi，hybの場合）#!/bin/sh# スクリプト内のシェルを指定#PJM -L "rscgrp=tutorial"#PJM -L "node=1"#PJM -L "elapse=10:00"#PJM –o OUT#PJM -j・・・ 以下省略 ・・・

22 解法 × 11 前処理 × 3 スケーリング（“前処理無し”，“スケーリング無し” も含める）‐ １１種類のデータ形式をサポート： CSR(CRS), CSC(CCS), MSR, DIA, …
‐ 倍精度，４倍精度演算サポート・・・ 以下省略 ・・・$ more $HOME/sample/script/pjsub_seq.sh‐ 開発元： JST CREST SSIプロジェクト（小武守，藤井，西田，長谷川 )新バージョンの予定は無さそうだが，随時，マイナー修正が行われている．
様々な計算機へのインストールが容易．
‐ 線形方程式求解アルゴリズムの逐次版，並列版（OpenMP：マルチスレッド環境, MPI：マルチプロセス環境）

Component of Lis: [Lis‐1.3.0 ver.]● スクリプトの拡張子は任意(† Lis の特徴：国内研究者の成果，または，力を入れて研究している)

反 復 解 法● ジョブ実行スクリプトの左端の“ #PJM ”
記号は，これに続くハイフンで始まる文字列がpjsubコマンドのオプションであることを宣言．
● 赤太文字の項目をチェック（pjsubオプション「-L “node=1”」でノード数を指定）● mpi，hybでは，ノード数=4 と設定

前 処 理

非定常反復解法
CG, BiCG, CGS, BiCGStab, BiCGStab(l ),GPBiCG†, TFQMR, Orthomin(m), GMRES(m)BiCGSafe†, CR, BiCR†, CRS†, BiCRStab†,
GPBiCR†, BiCRSafe†, FGMRES(m), IDR(s),MINRES

スケーリング（広義の前処理の一種）None,D 1 Ax  D 1b.

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

定常反復解法● omp，hybでは，スレッド数=8 と設定 13

None, (Point) Jacobi, ILU(k), SSOR, Hybrid†,
I+S†, SAINV, SA-AMG†,
Crout ILU, ILUT, Additive Schwarz

14

Jacobi, Gauss-Seidel, SOR

2. Lisインストールの準備

1. Lisの制限事項（1.3.0版）
前処理

Lis実習環境の準備と逐次(seq)版の準備：– Jacobi, SSOR 以外の前処理が選択され, かつ行列A がCSR 形式でない場合,
前処理作成時にCSR 形式の行列A が作成される.

$ tar zxf $HOME/sample/library/lis-1.3.0.tar.gz# その他の方法の一例として：（ $ gunzip -c $HOME/sample/library/lis-1.3.0.tar.gz | tar xf - ）– 線型方程式解法としてBiCG 法が選択された場合, SA-AMG 前処理は使用できない.
– SA-AMG 前処理はマルチスレッド環境では使用できない.

$ ls
lis-1.3.0/– SAINV 前処理の前処理行列作成部分は逐次実行される.

sample/$ cp –rp lis-1.3.0 lis-1.3.0_seq

4倍精度演算$ cd lis-1.3.0_seq$ ls
AUTHORSINSTALLCOPYING
Makefile.am
ChangeLog Makefile.in– 線型方程式解法のうち, Jacobi, Gauss-Seidel, SOR, IDR(s) 法では使用できない.
– 固有値解法のうち, CG, CR 法では使用できない.
– Hybrid 前処理での内部反復解法のうち, Jacobi, Gauss-Seidel, SOR 法では使用できない.

NEWSREADME
aclocal.m4# Lis逐次(seq)版の環境準備

build.sh*config/configure*configure.indoc/include/

src/
test/
win/

– I+S, SA-AMG 前処理では使用できない.

行列格納形式– VBR 形式はマルチプロセス環境では使用できない.
– マルチプロセス環境において必要な配列を直接定義する場合は, CSR 形式で作成しなければならない. 目的の格納形式を使用するには, lis_matrix_convert を使用してCSR 形式から変換する.

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

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

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

Lis逐次(seq)版のconfigure：

Lis逐次(seq)版のmake：$ pwd$ make/home/t00003/lis-1.3.0_seq$ ./configure --enable-saamg TARGET=fujitsu_fx10_cross # SAAMG前処理を有効にする
Lis開発側で実績のあるマシンの場合は，TARGETにて計算機環境を指定し，最適化可能だが，実績の無いマシンでも，大抵はインストール＆使用可能である．

・・・ 途中省略 ・・・

C Compiler
C Flags
C Libraries
F90 Compiler
F90 Flags
F90 Libraries

= no
= yes
= no
=
= no
= no
= no
= no

= fccpx= -Kfast,ocl,preex -w -DHAVE_CONFIG_H= -lm= frtpx= -Kfast,ocl,preex -Cpp -fs -DZERO_ORIGIN=1
=

途中省略

・・・

・・・ 出力の最後部 ・・・/bin/sh ../libtool --mode=link frtpx-mlcmain=main -o etest5 -lm
etest5.o -L../src llis -lm
libtool: link: frtpx -mlcmain=main -o etest5 etest5.o -L/home/t00003/lis-1.3.0_seq/src/home/t00003/lis-1.3.0_seq/src/.libs/liblis.a -lm
fccpx -DHAVE_CONFIG_H -I. -I../include -I.. -I../include-Kfast,ocl,preex -w -cetest6.c/bin/sh ../libtool --mode=link frtpx-mlcmain=main -o etest6 -lm
etest6.o -L../src llis -lm
libtool: link: frtpx -mlcmain=main -o etest6 etest6.o -L/home/t00003/lis-1.3.0_seq/src/home/t00003/lis-1.3.0_seq/src/.libs/liblis.a -lm
make[1]: ディレクトリ `/home/t00003/lis-1.3.0_seq/test' から出ます
make[1]: ディレクトリ `/home/t00003/lis-1.3.0_seq' に入ります
make[1]: `all-am' に対して行うべき事はありません.
make[1]: ディレクトリ `/home/t00003/lis-1.3.0_seq' から出ます

・・・ 出力の最後部 ・・・
Enable Fortran Interface
Enable SA-AMG Preconditioner
Enable MPI
Enable OpenMP
Enable Quad Precision
Enable Long Long Int
Enable Dynamic Linking
Enable Profiling・・・

赤文字の５項目をチェック

最後部でErrorが出ていなければOK● configure実行後は，作業ディレクトリやパスを変更しないように（生成されるMakefile中では絶対17パスが用いられている）．

18

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

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

Lis逐次(seq)版のインストールチェック１(isubコマンドによるmake check)：

Lis逐次(seq)版のインストールチェック２(isubコマンドによるmake check)：$ isub–seq”make check”＃クロスコンパイルのため，ログインノードでのジョブ実行は不可能

make check＃実行させたいコマンド列の両側をダブルクォーテーションで囲む[INFO] PJM 0000 pjsub Job 437631 submitted.
ジョブID# ジョブ実行中の表示例$ pjstatOakleaf-FX scheduled stop time: 2012/12/28(Fri) 09:00:00 (Remain: 14days 17:03:11)
JOB_ID
437631

JOB_NAME
isub

STATUS PROJECT RSCGROUP START_DATEQUEUED gt00
lecture --/-- --:--:--$ pjstat

ELAPSE
00:00:00

TOKEN NODE:COORD(0.2) 1# ジョブ終了後の表示例

Oakleaf-FX scheduled stop time: 2012/12/28(Fri) 09:00:00 (Remain: 14days 17:03:06)
No unfinished job found.

$ type lrt
lrt is aliased to `ls -lrt‘$ lrt
・・・ 省略 ・・・-rwxr-xr-x 1 t00003 gt00
-rw-r--r-- 1 t00003 gt00
-rw-r--r-- 1 t00003 gt00
drwxr-xr-x 12 t00003 gt00
drwxr-xr-x 4 t00003 gt00
-rw-r--r-- 1 t00003 gt00

76621 12月 13 15:41 2012 config.status*
31775 12月 13 15:41 2012 config.log
48946 12月 13 15:41 2012 Makefile4096 12月 13 15:45 2012 src/
4096 12月 13 15:56 2012 test/
4684 12月 13 15:56 2012 OUT＃ make check の実行結果
19$ cat OUT
・・・ 出力の途中部分 ・・・
make check-TESTS
make[2]: Entering directory `/home/t00003/lis-1.3.0_seq/test'=== Running test test.shchecking linear solvers...
number of processes = 1
matrix size = 100 x 100 (460 nonzero entries)
initial vector x = 0
precision : double
solver: BiCG 2
precon: nonestorage: CSR
lis_solve : normal end

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

BiCG: number of iterations= 15 (double = 15, quad = 0)
BiCG: elapsed time= 4.179478e-04 sec.
BiCG:
preconditioner= 5.984306e-05 sec.
BiCG:
matrix creation= 3.099442e-05 sec.
BiCG:
linear solver= 3.581047e-04 sec.
BiCG: relative residual 2-norm = 1.689581e-16
checking eigensolvers...
・・・
20

3. Lis逐次版のインストール(5/5)
Lis逐次(seq)版のインストールチェック(番外編：pjsubコマンドによるmake check)：

マニュアルには，最後に上記を実行するように記述されているが，これは，正式にシステムにインストールする（例えば，/usr/local/下など）場合に実行すれば良い．
この実習では，既に作成されている “test1 ” というプログラムを用いて説明するため，
特にインストール先を指定しないでも良い．
現状で，生成済みの関連ファイル群は下記のディレクトリ構成である：＃クロスコンパイルのため，ログインノードでのジョブ実行は不可能[INFO] PJM 0000 pjsub Job 437667 submitted.

$ lrt
・・・ 省略 ・・・-rwxr-xr-x 1 t00003 gt00
-rw-r--r-- 1 t00003 gt00
-rw-r--r-- 1 t00003 gt00
drwxr-xr-x 12 t00003 gt00
-rwxr-xr-x 1 t00003 gt00
drwxr-xr-x 4 t00003 gt00
-rw-r--r-- 1 t00003 gt00

make installは実行しなくても利用可能：$ make install$HOME/sample/script/pjsub_seq.sh をコピーして下記のとおり実行すれば良いが，ここでは実行しなくとも良い．
$ pjsub pjsub_seq.sh

3. Lis逐次版のインストール(補足)76621 12月 13 15:41 2012 config.status*
31775 12月 13 15:41 2012 config.log
48946 12月 13 15:41 2012 Makefile4096 12月 13 15:45 2012 src/
628 12月 13 16:21 2012 pjsub_seq.sh*4096 12月 13 16:21 2012 test/
4684 12月 13 16:21 2012 OUT＃ make check の実行結果$HOME/lis-1.3.0_seq/ 下：
include/lis_config.h，lis.h，lisf.h 他
src/.libs/
liblis.a liblis.la@ liblis.lai

当実習で
は不要！

インストール先のデフォルトは /usr/local/ 下であるが，それを変更する場合には，configure のステップまで戻って，下記のオプションを指定する．
$ ./configure --prefix={インストールディレクトリ}$ cat OUT
・・・ 出力の途中部分 ・・・

configureで指定したTARGET=fujitsu_fx10_crossの実体

isubの実行結果のOUTを上書き更新していることにも注目

CC=fccpx FC=frtpx CFLAGS="-Kfast,ocl,preex“ FCFLAGS="-Kfast,ocl,preex -Cpp-fs" FCLDFLAGS="-mlcmain=main“ ac_cv_sizeof_void_p=8 cross_compiling=yes
ax_f77_mangling="lower case, underscore, no extra underscore“
configureコマンドの引数として，これら全ての環境変数を上記のとおり指定したことと等価．
特に，赤文字で強調したコンパイラと cross_compiling=yes では，クロスコンパイリング
22
環境であることにも注目．

21

4. isubコマンドでジョブ実行(逐次版)

4. isubコマンドでジョブ実行(逐次版)

isubコマンドによるジョブ実行１：

isubコマンドによるジョブ実行２：$ cd $HOME/lis-1.3.0_seq/test$ isub -seq ”./test1 testmat.mtx 0 sol2.txt rsｄ2.txt -i 3 -p 1”$ isub -seq ”./test1 testmat.mtx 0 sol1.txt rsd1.txt”./test1 testmat.mtx 0 sol2.txt res2.txt -i 3 -p 1
[INFO] PJM 0000 pjsub Job 437797 submitted.

./test1 testmat.mtx 0 sol1.txt rsd1.txt[INFO] PJM 0000 pjsub Job 437789 submitted.

$ lrt

$ lrt
-rw-r--r-- 1 t00003 gt00
-rw-r--r-- 1 t00003 gt00
-rw-r--r-- 1 t00003 gt00

3242 12月 13 18:08 2012 sol1.txt # 解ベクトル1(solution)
208 12月 13 18:08 2012 rsd1.txt # 相対残差1(residual)のノルム収束履歴# std{out,err}ファイル
498 12月 13 18:08 2012 OUT

3242 12月 13 18:08 2012 sol1.txt
208 12月 13 18:08 2012 rsd1.txt
3242 12月 13 18:14 2012 sol2.txt
208 12月 13 18:14 2012 rsｄ2.txt
493 12月 13 18:14 2012 OUT# 解ベクトル1(solution)# 相対残差1(residual)のノルム収束履歴# 解ベクトル2(solution)# 相対残差2(residual)のノルム収束履歴# std{out,err}ファイル$ cat OUT

$ cat OUT
number of processes = 1
matrix size = 100 x 100 (460 nonzero entries)
initial vector x = 0
precision : double
solver: BiCG 2
precon: nonestorage: CSR
lis_solve : normal end-rw-r--r-- 1 t00003 gt00
-rw-r--r-- 1 t00003 gt00
-rw-r--r-- 1 t00003 gt00
-rw-r--r-- 1 t00003 gt00
-rw-r--r-- 1 t00003 gt00

※ 「Lis逐次版のインストール(4/5)」のスライドと同じ出力（make checkと同じコマンド，同じ
オプションを指定したということ）．

BiCG: number of iterations= 15 (double = 15, quad = 0)
BiCG: elapsed time= 4.179478e-04 sec.
BiCG:
preconditioner= 5.769730e-05 sec.
BiCG:
matrix creation= 3.099442e-05 sec.
BiCG:
linear solver= 3.602505e-04 sec.
BiCG: relative residual 2-norm = 1.689581e-16

23

number of processes = 1
matrix size = 100 x 100 (460 nonzero entries)
initial vector x = 0
precision : double
solver: CGS 3precon: Jacobi
storage: CSR
lis_solve : normal end
CGS: number of iterations= 15 (double = 15, quad = 0)
CGS: elapsed time
= 1.056194e-03 sec.
CGS:
preconditioner= 7.200241e-05 sec.
CGS:
matrix creation= 3.099442e-05 sec.
CGS:
linear solver= 9.841919e-04 sec.
CGS: relative residual 2-norm = 1.004890e-14

24

5. 色々な解法で解いてみよう

5. いま，Lisで何をやったのだろうか$ isub–seq”./test1”＃ 実行結果 （OUTの内容）

Usage: test1 matrix_filename rhs_setting solution_filename residual_filename [options]
係数行列のファイルを指定
右辺項ベクトルファイル，あるいは，自動生成機能（主にアルゴリズムの性能評価用途）のIDを指定
0 : 拡張Matrix Market形式（右辺ベクトルが行列データに含まれている）を用いる
1 : b = (1,1,...,1)T を用いる
2 : b = A ×(1,1,...,1)T を用いる
solution_filename: 解ベクトルのファイル名
residual_filename: （アルゴリズム中のベクトルを用いた）相対残差ベクトルの２ノルムの履歴[options]:
解法と前処理，その他のオプション指定
matrix_filename:rhs_setting:# 参考： マニュアル 「テストプログラム」の「test1」小節$ more testmat.mtx%%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

25

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

解法

選択オプション

補助オプション

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

CG‐i {cg | 1}

BiCG

‐i {bicg | 2}CGS

‐i {cgs | 3}

BiCGSTAB‐i {bicgstab | 4}

BiCGSTAB( l )‐i {bicgstabl | 5}

GPBiCG‐i {gpbicg | 6}

TFQMR

‐i {tfqmr | 7}

Orthomin(m)‐i {orthomin | 8}

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

GMRES(m)‐i {gmres | 9}

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

Jacobi‐i {jacobi | 10}

Gauss‐Seidel‐i {gs | 11}SOR‐i {sor | 12}‐omega [1.9]緩和係数ωの値(0<ω<2)‐ell [2]解法の l の値

● 上記も含め，全22種類の解法が提供されている．

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

前処理

選択オプション

補助オプション

オプション

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

内容なし‐p {none | 0}

‐maxiter [1000]

最大反復回数

Jacobi‐p {jacobi | 1}

‐tol [1.0e‐12]

収束判定基準（相対残差ノルム）

ILU( k )‐p {ilu | 2}

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

残差の画面表示SSOR‐p {ssor | 3}

‐ssor_w [1.0]緩和係数ωの値(0<ω<2)‐scale [0]Hybrid‐p {hybrid | 4}

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

I+S‐p {is | 5}

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

スケーリングの選択．
※ スケーリングされた係数行列，右辺項ベクトルは，元のデータに上書きされて，プログラムで実行される．

SAINV‐p {sainv | 6}

‐sainv_drop [0.05]‐initx_zeros [true]

初期ベクトル x0 の設定（デフォルトは全要素の値が 0）

SA‐AMG‐p {saamg | 7}

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

Crout ILU

‐p {iluc | 8}

‐iluc_drop [0.05]

ドロップ基準，他ILUT‐p {ilut | 9}

‐ilut_drop [0.05]

ドロップ基準，他

Additive Schwarz‐adds true

‐adds_iter [1]

繰り返し回数

ドロップ基準● 上記の補助オプションは一部のみ掲載している．詳細はマニュアル参照．

26‐omp_num_threads [t] OMP実行スレッド数．t は最大スレッド数 *‐storage [0]行列格納形式‐storage_block [2]BSR，BSCのブロックサイズ* FXでは，デフォルト設定にて FLIB_FASTOMP=TRUE となっているため，当オプションによるプログラム
からのスレッド数指定ではエラーとなります．FXで当オプションを使用する場合，あらかじめ，環境変数にて「FLIB_FASTOMP=FALSE」「FLIB_CNTL_BARRIER_ERR=FALSE」を指定して下さい．
27● 上記のオプションの内容の詳細はlisのマニュアル参照．

28

5. 色々なデータで実行してみよう

6. 入力データ形式について$ pwd● Matrix Market(MM)形式は，データの圧縮度は低いものの，行列のインデックスと値との対応を確認し易く，また，多くのアプリケーションでもMM形式をサポートしたりツール提供されているので，初級者や中級者にとっては都合が良いと思われる．

/home/t00003/lis-1.3.0_seq/test$ ls$HOME/sample/

library/

$ ln
$ ls

mat/

-s
-l

script/setup.envsh/● Lisのtest1.cでは，MM形式のデータを読込み，CSR形式に変換する．MPIのデータ分散にも対応している．

$HOME/sample/mat .
mat● 対称行列の場合には，対角要素から下三角の要素のみが格納されている．この場合，ヘッダには，半角スペースに続いて「symmetric」の文字列が記載されている．
→ mat/下の「494_bus.mtx」「nos3.mtx」が該当

lrwxrwxrwx 1 t00003 gt00 23 12月 13 18:38 2012 mat -> /home/t00003/sample/mat/$ ls mat/494_bus.mtx
494_bus.rsa
add20.mtx
add20.rua

add20_rhs1.mtx
add32.mtx
add32.rua
add32_rhs1.dat

add32_rhs1.mtx
add32_rhs2.mtx
fs_183_1.mtx
fs_183_1.rua

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

参考までに，matの後にディレクトリを示す「/」を付け忘れると：$ ls mat# これは，☆１のsetup.env 影響
mat@$ more $HOME/sample/script/pjsub_seq.sh

nos3.rsa● 非対称行列の場合には，「general」の文字列が記載されている．
→ mat/下の「add{20,32}.mtx」「fs_183_1.mtx」 「memplus.mtx」が該当● 拡張子
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
・・・

● 各々をバッチジョブ実行するときは，左端冒頭の「#」（コメントアウト）記号を削除

29

add32.mtx: MMサイトで提供されているMM形式の係数行列データ
add32_rhs1.mtx :〃MM形式の右辺項ベクトルデータ
add32.rua:〃
Harwell-Boeing(HB)形式の係数行列データ
MMサイトで提供されている add32_rhs1.mtx のヘッダとフォーマットでは，Lisで読み取れないので，
当実習では以下の２通りの修正データを用意してます．問題点と修正内容は，各自で確認して下さい．
add20, memplus等も同様ですが，サンプルではadd32のみ修正データを用意しています．
add32_rhs1.dat : add32_rhs1.mtxを元にPLAIN形式にしたファイル
add32_rhs2.mtx : Lisで読み取り可能なベクトル用拡張MM形式のデータ

30

7. Lis並列(Hybrid)版のインストール(1/3)

7. Lis並列版の準備

Lis並列(Hybrid)版のconfigure：
逐次版と同様，ディレクトリごとコピーして並列版の環境を構築する：$ cd$HOME/lis-1.3.0/test/# lis-1.3.0_seq ではないことに注意

$ ln

-s

# matをシンボリックリンク$ pwd/home/t00003/lis-1.3.0_hyb$HOME/sample/mat .

$ ./configure --enable-omp --enable-mpi --enable-saamg TARGET=fujitsu_fx10_cross# ここで omp と mpi の両方を指定する（もちろん，片方だけ指定も可能）# シンボリックリンクの確認$ ls -l mat

・・・ 途中省略 ・・・

lrwxrwxrwx 1 t00003 gt00 23 12月 13 18:42 2012 mat -> /home/t00003/sample/mat/$ cp –rp$ cp –rplis-1.3.0
lis-1.3.0

lis-1.3.0_omp
lis-1.3.0_mpi# LisのOpenMP(omp)版の環境準備# LisのMPI(mpi)版の環境準備$ cp –rplis-1.3.0

lis-1.3.0_hyb# LisのHybrid(mpi+omp)版の環境準備$ cd

・・・ 出力の最後部 ・・・
Enable Fortran Interface
Enable SA-AMG Preconditioner
Enable MPI
Enable OpenMP
Enable Quad Precision
Enable Long Long Int
Enable Dynamic Linking
Enable Profiling（または，「 cd ~ 」や「 cd $HOME 」など）

$ cd

lis-1.3.0_hyb

C Compiler
C Flags
C Libraries
F90 Compiler
F90 Flags
F90 Libraries
MPI Libraries
MPIRUN Script
31

= no
= yes
= yes
= yes
= no
= no
= no
= no

赤文字の8項目をチェック= mpifccpx= -Kfast,ocl,preex -w -DUSE_MPI -Kopenmp -DHAVE_CONFIG_H= -lm= mpifrtpx= -Kfast,ocl,preex -Cpp -fs -DZERO_ORIGIN=1 -DUSE_MPI -Kopenmp==● configure実行後は，作業ディレクトリやパスを= mpiexec

変更しないように（生成されるMakefile中では絶対32パスが用いられている）．

7. Lis並列(Hybrid)版のインストール(2/3)

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

Lis並列(Hybrid)版のmake：Lis並列(Hybrid)版のインストールチェック(make check)：$ isub -hyb ”make check”$ make・・・

途中省略＃Lisのmake checkでは，MPIプログラムであっても，mpiexecコマンドは不要（make checkで呼び出すtest/test.sh 内でmpiexecを記述しているため）$ cat OUT

・・・ 出力の途中部分 ・・・
make check-TESTS
make[2]: Entering directory `/home/t00003/lis-1.3.0_hyb/test'=== Running test test.sh・・・・・・ 出力の最後部 ・・・/bin/sh ../libtool --mode=link mpifrtpx -Kopenmp -mlcmain=main -o etest6 -lm
etest6.o L../src -llis -lm
libtool: link: mpifrtpx -Kopenmp -mlcmain=main -o etest6 etest6.o -L/home/t00003/lis1.3.0_hyb/src /home/t00003/lis-1.3.0_hyb/src/.libs/liblis.a -lm
make[1]: ディレクトリ `/home/t00003/lis-1.3.0_hyb/test' から出ます
make[1]: ディレクトリ `/home/t00003/lis-1.3.0_hyb' に入ります
make[1]: `all-am' に対して行うべき事はありません.
make[1]: ディレクトリ `/home/t00003/lis-1.3.0_hyb' から出ますchecking linear solvers...
number of processes = 2
max number of threads = 16
number of threads = 2
matrix size = 100 x 100 (460 nonzero entries)
initial vector x = 0
precision : double
solver: BiCG 2
precon: nonestorage: CSR
lis_solve : normal end＃Lisの make check では，ノード数（number of processes）=2 とスレッド数（number of threads）=2は固定されている．

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

BiCG: number of iterations= 15 (double = 15, quad = 0)
BiCG: elapsed time= 2.397419e-03 sec.
BiCG:
preconditioner= 1.086006e-04 sec.
BiCG:
matrix creation= 2.899906e-06 sec.
BiCG:
linear solver= 2.288818e-03 sec.
BiCG: relative residual 2-norm = 1.351335e-16

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

34

・・・

7. isubコマンドでジョブ実行(Hybrid版)

7. isubコマンドでジョブ実行(Hybrid版)

isubコマンドによるジョブ実行：

isubコマンドにてpjsubのオプション指定したジョブ実行：$ cd $HOME/lis-1.3.0_hyb/test$ isub -hyb ”mpiexec ./test1 testmat.mtx 0 sol2.txt rsd2.txt” -o OUT1 -L ”node=2”＃Hybrid版はMPIプログラムでもあるので，mpiexecコマンドが必要$ isub -hyb ”mpiexec ./test1 testmat.mtx 0 sol1.txt rsd1.txt”

mpiexec ./test1 testmat.mtx 0 sol2.txt rsd2.txt[INFO] PJM 0000 pjsub Job 438648 submitted.

mpiexec ./test1 testmat.mtx 0 sol1.txt rsd1.txt[INFO] PJM 0000 pjsub Job 437872 submitted.

$ pjstat$ lrt-rw-r--r-- 1 t00003 gt00
-rw-r--r-- 1 t00003 gt00
-rw-r--r-- 1 t00003 gt00

JOB_ID
4386483242 12月 13 19:00 2012 sol1.txt # 解ベクトル1(solution)
208 12月 13 19:00 2012 rsd1.txt # 相対残差1(residual)のノルム収束履歴# std{out,err}ファイル
547 12月 13 19:00 2012 OUT$ cat OUT
number of processes = 4
max number of threads = 16
number of threads = 8
matrix size = 100 x 100 (460 nonzero entries)
initial vector x = 0
precision : double
solver: BiCG 2
precon: nonestorage: CSR
lis_solve : normal end＃ isubコマンドのデフォルトでは，ノード数（number of processes）=4
スレッド数（number of threads）=8と指定している．
「0.FXでの実習環境設定(5/5)」参照

BiCG: number of iterations= 15 (double = 15, quad = 0)
BiCG: elapsed time= 3.740030e-03 sec.
BiCG:
preconditioner= 1.248000e-04 sec.
BiCG:
matrix creation= 2.699904e-06 sec.
BiCG:
linear solver= 3.615230e-03 sec.
BiCG: relative residual 2-norm = 2.840829e-16

35

JOB_NAME
isub

STATUS PROJECT RSCGROUP START_DATE
ELAPSERUNNING gt00
lecture (12/14 15:18:31) 00:00:00

TOKEN NODE:COORD
0.0 2

$ lrt
-rw-r--r-- 1 t00003 gt00
-rw-r--r-- 1 t00003 gt00
-rw-r--r-- 1 t00003 gt00
-rw-r--r-- 1 t00003 gt00
-rw-r--r-- 1 t00003 gt00
-rw-r--r-- 1 t00003 gt00

$ diff OUT OUT1

3242 12月 13 19:00 2012 sol1.txt
208 12月 13 19:00 2012 rsd1.txt
547 12月 13 19:00 2012 OUT
3242 12月 14 15:18 2012 sol2.txt
208 12月 14 15:18 2012 rsd2.txt
547 12月 14 15:18 2012 OUT1

# 解ベクトル1(solution)# 相対残差1(residual)のノルム収束履歴# std{out,err}ファイル# 解ベクトル2(solution)# 相対残差2(residual)のノルム収束履歴# std{out,err}ファイル2● pjsubオプションの指定は，（isubに限らず一般的なpjsub用の）

2c2
ジョブスクリプト内での指定よりも，pjsubコマンドの引数(option)< number of processes = 4
指定の方が優先される．
--> number of processes = 2＃ isubコマンドのデフォルトでは，
14,18c14,18< BiCG: elapsed time= 3.740030e-03 sec.
・ 出力ファイル： OUT< BiCG:
preconditioner= 1.248000e-04 sec.
・ ノード数（number of processes）=4< BiCG:
matrix creation= 2.699904e-06 sec.
< BiCG:
linear solver= 3.615230e-03 sec.
であるが，今回は pjsubコマンドのオプション< BiCG: relative residual 2-norm = 2.840829e-16
として –o OUT1 –L ”node=2” を指定した．
--36
・・・ 以下省略 ・・・

7. isubコマンドでジョブ実行(Hybrid版)

7. Lisのその他の並列(OMP)版(1/2)

isubコマンドにてpjsubのオプション指定したジョブ実行：$ isub -hyb –th 4 ”mpiexec ./test1 testmat.mtx 0 sol3.txt rsd3.txt” -o OUT2

OMP版： configureオプションとisubコマンドの引数(-omp)に注目：$ lrt$ cd $HOME/lis-1.3.0_omp$ ./configure --enable-omp --enable-saamg TARGET=fujitsu_fx10_cross-rw-r--r-- 1 t00003 gt00
-rw-r--r-- 1 t00003 gt00
-rw-r--r-- 1 t00003 gt00
-rw-r--r-- 1 t00003 gt00
-rw-r--r-- 1 t00003 gt00
-rw-r--r-- 1 t00003 gt00
-rw-r--r-- 1 t00003 gt00
-rw-r--r-- 1 t00003 gt00
-rw-r--r-- 1 t00003 gt00

3242 12月 13 19:00 2012 sol1.txt
208 12月 13 19:00 2012 rsd1.txt
547 12月 13 19:00 2012 OUT
3242 12月 14 15:18 2012 sol2.txt
208 12月 14 15:18 2012 rsd2.txt
547 12月 14 15:18 2012 OUT1
3242 12月 16 13:46 2012 sol3.txt # 解ベクトル3(solution)
208 12月 16 13:46 2012 rsd3.txt # 相対残差3(residual)のノルム収束履歴547 12月 16 13:46 2012 OUT2
# std{out,err}ファイル3

C Compiler
C Flags
C Libraries
F90 Compiler
F90 Flags
F90 Libraries$ diff OUT OUT2
4c4
< number of threads = 8--> number of threads = 4
14,18c14,18< BiCG: elapsed time= 3.740030e-03 sec.
< BiCG:
preconditioner= 1.248000e-04 sec.
< BiCG:
matrix creation= 2.699904e-06 sec.
< BiCG:
linear solver= 3.615230e-03 sec.
< BiCG: relative residual 2-norm = 2.840829e-16--・・・ 以下省略 ・・・

Enable Fortran Interface
Enable SA-AMG Preconditioner
Enable MPI
Enable OpenMP
Enable Quad Precision
Enable Long Long Int
Enable Dynamic Linking
Enable Profiling

= no
= yes
= no
= yes
= no
= no
= no
= no

= fccpx= -Kfast,ocl,preex -w -Kopenmp -DHAVE_CONFIG_H= -lm= frtpx= -Kfast,ocl,preex -Cpp -fs -DZERO_ORIGIN=1 -Kopenmp=$ make$ isub –omp ”make check”＃ isubコマンドのデフォルトでは，・ ノード数（number of threads）=8であるが，今回は pjsubコマンドのオプション
として –th 4 を指定した．
37

7. Lisのその他の並列(OMP)版(2/2)

number of processes = 1
max number of threads = 16
number of threads = 2
matrix size = 100 x 100 (460 nonzero entries)
・・・ 中略 ・・・
storage: CSR
lis_solve : normal end＃詳細な確認は次のスライド(2/2)を参照
38

7. Lisのその他の並列(MPI)版(1/2)

OMP版： isubコマンドの引数(-omp)とOUTの内容に注目：

MPI版： configureオプション，isubコマンドの引数(-mpi)，OUTの内容に注目：$ cd $HOME/lis-1.3.0_omp/test$ isub -omp ”./test1 testmat.mtx 0 sol1.txt rsd1.txt”$ cd $HOME/lis-1.3.0_mpi$ ./configure --enable-mpi --enable-saamg TARGET=fujitsu_fx10_crossEnable Fortran Interface
Enable SA-AMG Preconditioner
Enable MPI
Enable OpenMP
Enable Quad Precision
Enable Long Long Int
Enable Dynamic Linking
Enable Profiling$ cat OUT
number of processes = 1
max number of threads = 16
number of threads = 8
matrix size = 100 x 100 (460 nonzero entries)
initial vector x = 0
precision : double
solver: BiCG 2
precon: nonestorage: CSR
lis_solve : normal end

C Compiler
C Flags
C Libraries
F90 Compiler
F90 Flags
F90 Libraries
MPI Libraries
MPIRUN Script

BiCG: number of iterations= 15 (double = 15, quad = 0)
BiCG: elapsed time= 2.825022e-03 sec.
BiCG:
preconditioner= 1.111031e-04 sec.
BiCG:
matrix creation= 9.536743e-07 sec.
BiCG:
linear solver= 2.713919e-03 sec.
BiCG: relative residual 2-norm = 1.533189e-16

= no
= yes
= yes
=
= no
= no
= no
= no

= mpifccpx= -Kfast,ocl,preex -w -DUSE_MPI -DHAVE_CONFIG_H= -lm= mpifrtpx= -Kfast,ocl,preex -Cpp -fs -DZERO_ORIGIN=1 -DUSE_MPI=== mpiexec$ make$ isub –mpi ”make check”

39

number of processes = 2
matrix size = 100 x 100 (460 nonzero entries)
・・・ 中略 ・・・
storage: CSR
lis_solve : normal end＃詳細な確認は次のスライド(2/2)を参照
40

8. test1 処理フローの概要(1)7. Lisのその他の並列(MPI)版(2/2)

LIS_INT main(LIS_INT argc, char* argv[] ) {lis_initialize

MPI版： isubコマンドの引数(-mpi)，OUTの内容に注目：$ cd $HOME/lis-1.3.0_mpi/test$ isub -mpi ”mpiexec ./test1 testmat.mtx 0 sol1.txt rsd1.txt”$ cat OUT
number of processes = 4
matrix size = 100 x 100 (460 nonzero entries)
initial vector x = 0
precision : double
solver: BiCG 2
precon: nonestorage: CSR
lis_solve : normal end

係数行列 A をmalloc
ベクトル b, x をmalloc
A(MMかHB形式), b(ファイル入力の場合), x(初期値)を読込み

lis_matrix_set_type
lis_matrix_convert
右辺項 b の作成処理

lis_solver_get_timeex
lis_solver_get_residualnorm

行列データをCSRへ変換〃
MMとHB以外の
形式は，マニュアル
5節「行列格納形
求解ルーチン
式」のサンプルや，
特にCSRを扱った
タイマーの値をget
例としてtest[2‐4].c
真の残差ノルム算出
も参照

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

lis_finalize();
return 0;

各関数の引数は，実際のtest1.c
ファイルの記述を確認して下さい．

lis_solve(A,b,x,solver);

BiCG: number of iterations= 15 (double = 15, quad = 0)
BiCG: elapsed time= 2.942821e-03 sec.
BiCG:
preconditioner= 3.380026e-05 sec.
BiCG:
matrix creation= 2.900138e-06 sec.
BiCG:
linear solver= 2.909021e-03 sec.
BiCG: relative residual 2-norm = 1.550156e-16

41

}

42

8. プログラム内で係数行列作成する場合： Lisのtest4,5の例

8. test1 処理フローの概要(2)int lis_solve( ) {lis_precon_create
lis_solve_kernel(A, b, x, solver, precon);}lis_matrix_create
lis_vector_create
lis_input(A,b,x,argv[1]);前処理演算の設定
求解のための事前処理

TAx  b b  Ax, x  1.0,1.0, ,1.0

＃MPIプログラム（Hybridも）の場合，mpiexecコマンドの引数として実行
int lis_solve_kernel( ) {test1の引数をチェック$ ./test4 [options]0 2 1 1 2 1




A
  


1 2 1
 0
1 2 

x0を設定
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);

真の残差ベクトルの2ノルム（絶対残差）※2 scaling後の系に対して算出(cf. ※1)}43係数行列Aの次数n=12（プログラム内
で定義）の3重対角行列$ ./test5 n gamma [options]

2 1
0 2 1 0 2 1A   


 0 2

 0
0

0





1

2

係数行列Aの次数nのToeplitz行列44

8. test4 係数行列作成の概要

8. test5 係数行列作成の概要

LIS_INT main(LIS_INT argc, char* argv[] ) {lis_initialize

LIS_INT main(LIS_INT argc, char* argv[] ) {lis_initialize

lis_matrix_create(LIS_COMM_WORLD,&A); 係数行列 A をmalloc
lis_matrix_set_size(A,0,n);
A のglobalな行数サイズ（n）を設定
lis_matrix_get_size(A,&n,&gn);
A のlocal size(n)とglobal size(gn)情報を取得
lis_matrix_get_range(A,&is,&ie);
分割後のAの当該プロセッサ内での領域を与える
for(i=is;i<ie;i++)test4 では配列Aに値を設定して係数行列を作成{
if( i>0 ) lis_matrix_set_value(LIS_INS_VALUE,i,i-1,-1.0,A); i, j=i-1 に -1.0を設定
if( i<gn-1 ) lis_matrix_set_value(LIS_INS_VALUE,i,i+1,-1.0,A); i, j=i+1 に -1.0を設定
lis_matrix_set_value(LIS_INS_VALUE,i,i,2.0,A);i, j=i に 2.0を設定}
lis_matrix_set_type(A,LIS_MATRIX_CSR); 行列Aの格納形式をCSRとして設定する
lis_matrix_assemble(A);行列Aをプログラム内で使用可能にする・・・ 途中省略 ・・・
lis_finalize();
return 0;}45

gn = atoi(argv[1]);gamma = atof(argv[2]);１番目の引数にて行列サイズを設定
２番目の引数にてガンマの値を設定

lis_matrix_create(LIS_COMM_WORLD,&A); 係数行列 A をmalloc
lis_matrix_set_size(A,0,gn);
A のglobalな行数サイズ（gn）を設定
lis_matrix_get_size(A,&n,&gn);
A のlocal size(n)とglobal size(gn)情報を取得
lis_matrix_malloc_csr(n,3*n,&ptr,&index,&value);
lis_matrix_get_range(A,&is,&ie);
分割後のAの当該プロセッサ内での領域を与える
k = 0;
ptr[0] = 0;
for(ii=is;ii<ie;ii++)
test5 では，直接，CSR形式としてデータ格納{
if( ii>1 ) { jj = ii - 2; index[k] = jj; value[k++] = gamma;}
if( ii<gn-1 ) { jj = ii + 1; index[k] = jj; value[k++] = 1.0;}
index[k] = ii; value[k++] = 2.0;
ptr[ii-is+1] = k;}lis_matrix_set_csr(ptr[ie-is],ptr,index,value,A); CSR形式の配列を行列Aに関連付ける
lis_matrix_assemble(A);行列Aをプログラム内で使用可能にする・・・ 以下省略 ・・・
46

講習内容

9. PETScの概要（線形方程式版）[http://www.mcs.anl.gov/petsc/]

前半の内容

PETSc (Portable, Extensible Toolkit for Scientific Computation)

後半の内容

0. 準備，FXでの実習環境設定

9. PETScの概要

1. Lisの概要

10. FXでのPETSc利用の準備

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

11. FXでのPETSc-逐次版利用例

3. Lis逐次版のインストール

12. FXでのPETSc-MPI版利用例

4. ジョブ実行（逐次版）

13. PETSc-複素数演算例題

5. Lis使用に関する解説

14. PETSc使用に関する解説

6. 入力データ形式について

15. PETScテストプログラムの概要

7. Lis並列版のインストール・実行

16. KSPに関する盲点

8. Lisテストプログラムの概要

17. 両ライブラリの習得要領‐ 開発元： ANL : Argonne National Laboratory
異なるバージョン間の非互換が大きいので注意．
インストール（特にconfigure）は少々難．
‐ 線形方程式求解アルゴリズムの逐次版，並列版（MPIのみ）
14 反復解法 × 10 前処理（“前処理無し” も含める） ＋ 2 直接解法
‐ 複数のデータ形式をサポート：
AIJ(CSR形式のこと), Block AIJ, Symmetric Block AIJ, Dense, Matrix‐Free

Component of PETSc: [PETSc‐3.3‐p4 ver.]反 復 解 法前 処 理

非定常反復解法
Richardson, Chebychev, CG, BiCG, GMRES,FGMRES, Deflated GMRES, GCR, BiCGStab,CGS, TFQMR(1), TFQMR(2), CR, LSQR

None, (Point) Jacobi, Block Jacobi,
SOR (and SSOR), SOR with Eisenstat trick,
IC, ILU, Additive Schwarz, Linear solver,
Combination of preconditioner

直接解法
LU, Cholesky
49

10. FXでのPETSc利用の準備(1/2)

10. FXでのPETSc利用の準備(2/2)

PETSc実習環境の準備：

PETSc実習環境の準備：

Oakleaf-fx のPETScインストールディレクトリ：/usr/local/PETSc/3.3-p4/$ module list
Currently Loaded Modulefiles:
1) TCSuite/GM-1.2.1-04$ pwd$ module avail PETSc/home/t00003$ cp -rp ~/sample/library/petsc-3.3-p4_samp$ cd petsc-3.3-p4_samp
$ ls

.

ex1.c

ex23.cex11.c@50

ex11_mod.c

ex11_org.c

ex1_mod.c＃テキスト作成後にTCSuiteの版数が上がっている場合があるが，当講習会の範囲では，特に気にしないで良い（見込み）．

------------------- /usr/share/Modules/modulefiles/apps -------------------PETSc/3.0.0
PETSc-complex/3.3-p4(default)
PETSc/3.3-p4(default)
makefile＃ex1_mod.c と ex11_mod.c 以外は全て oakleaf-fx の/usr/local/PETSc/3.3-p4/src/ksp/ksp/examples/tutorials/の下のファイルをコピーしたもの（ex11_org.cのみ，ファイル名変更している）．

$ module load PETSc$ module list

（ 「 $ module load PETSc/3.3-p4 」でも可 ）Currently Loaded Modulefiles:
1) TCSuite/GM-1.2.1-04
2) PETSc/3.3-p4$ env | grep -i petsc$ ls -l ex11.c
lrwxrwxrwx 1 t00003 gt00 10

2月 19 16:08 2013 ex11.c -> ex11_mod.c＃ex11.cは ex11_mod.cの方を使用する．
PETSc提供のオリジナルex11.c（ex11_org.cのこと）は不完全なコード（FX固有の問題では無い）．

51

PETSC_ARCH=linux
LIBRARY_PATH=/usr/local/PETSc/3.3-p4/lib
FORT90CPX=-I/usr/local/PETSc/3.3-p4/include/finclude -L/usr/local/PETSc/3.3-p4/lib
INCLUDE_PATH=/usr/local/PETSc/3.3-p4/include
FCCpx_ENV=-I/usr/local/PETSc/3.3-p4/include -L/usr/local/PETSc/3.3-p4/lib
PWD=/home/t00003/petsc-3.3-p4_samp_LMFILES_=/usr/share/Modules/modulefiles/utils/TCSuite/GM-1.2.103:/usr/share/Modules/modulefiles/apps/PETSc/3.3-p4
LOADEDMODULES=TCSuite/GM-1.2.1-03:PETSc/3.3-p4
fccpx_ENV=-I/usr/local/PETSc/3.3-p4/include -L/usr/local/PETSc/3.3-p4/lib
PETSC_DIR=/usr/local/PETSc/3.3-p4

52

11. FXでのPETSc‐逐次版利用例(1/3)

11. FXでのPETSc‐逐次版利用例(2/3)

ex1.c ： PETScの逐次版サンプルプログラム

ex1.c ： PETScの逐次版サンプルプログラム$ module list

$ isubCurrently Loaded Modulefiles:
1) TCSuite/GM-1.2.1-04
2) PETSc/3.3-p4–seq”mpiexec ./ex1”–mpi”mpiexec –n 1 ./ex1”

または$ isub$ pwd/home/t00003/petsc-3.3-p4_samp

ex1は逐次版であるが，MPIコードの1ノード用としてコーディングされている．
したがって，isubコマンドでもMPIの1ノードとして実行する．

$ cat OUTKSP Object: 1 MPI processes
type: gmres
GMRES: restart=30, using Classical (unmodified) Gram-Schmidt
Orthogonalization with no iterative refinement
GMRES: happy breakdown tolerance 1e-30
maximum iterations=10000, initial guess is zerotolerances: relative=1e-05, absolute=1e-50, divergence=10000
left preconditioning
using PRECONDITIONED norm type for convergence test
PC Object: 1 MPI processes
type: jacobi
linear system matrix = precond matrix:
Matrix Object:
1 MPI processes
type: seqaij
rows=10, cols=10
total: nonzeros=28, allocated nonzeros=50
total number of mallocs used during MatSetValues calls =0
not using I-node routines$ make ex1
・・・ makeのメッセージ出力 ・・・＃途中で下記のメッセージが表示されるが，特に気にしないでも良い（見込み）．
oakleaf-fx のバージョンが古いということだが，当実習内での不都合は無い（見込み）．
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++The version of PETSc you are using is out-of-date, we recommend updating to the newreleaseAvailable Version: 3.3.5
Installed Version: 3.3.4
http://www.mcs.anl.gov/petsc/download/index.html++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
・・・ 以下省略 ・・・

53

54

11. FXでのPETSc‐逐次版利用例(3/3)

12. FXでのPETSc‐MPI版利用例(1/2)

ex1.c ： PETScの逐次版サンプルプログラム

ex23.c ： PETScのMPI版サンプルプログラム$ isub –seq ”mpiexec ./ex1 –pc_type sor –ksp_type cgs”

$ ls
OUT

$ cat OUTKSP Object: 1 MPI processes
type: cgs
maximum iterations=10000, initial guess is zerotolerances: relative=1e-05, absolute=1e-50, divergence=10000
left preconditioning
using PRECONDITIONED norm type for convergence test
PC Object: 1 MPI processes
type: sor
SOR: type = local_symmetric, iterations = 1, local iterations = 1, omega= 1linear system matrix = precond matrix:
Matrix Object:
1 MPI processes
type: seqaij
rows=10, cols=10
total: nonzeros=28, allocated nonzeros=50
total number of mallocs used during MatSetValues calls =0
not using I-node routines
Norm of error 1.57716e-06, Iterations 6$ isub –seq ”mpiexec ./ex1 –help”＃ サンプルコードでは，このようにしてhelpを参照できる
55

ex1*

ex11.c@ex11_mod.c

ex11_org.c$ module list
Currently Loaded Modulefiles:
1) TCSuite/GM-1.2.1-04
2) PETSc/3.3-p4

ex11.cex1_mod.c

ex23.cmakefile＃実習状況により，表示されるファイル群は異なり
ますが，少なくとも赤文字のファイルは必須．

$ make ex23
・・・ 省略 ・・・$ isub -mpi ”mpiexec ./ex23”$ cat OUT

ex23はMPI並列版であるので，複数ノードで実行．

KSP Object: 4 MPI processes
type: gmres
GMRES: restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization withno iterative refinement
GMRES: happy breakdown tolerance 1e-30
maximum iterations=10000, initial guess is zerotolerances: relative=1e-07, absolute=1e-50, divergence=10000
left preconditioning
using PRECONDITIONED norm type for convergence test
PC Object: 4 MPI processes
type: jacobi
linear system matrix = precond matrix:
Matrix Object:
4 MPI processes
type: mpiaij
rows=10, cols=10
total: nonzeros=28, allocated nonzeros=70
total number of mallocs used during MatSetValues calls =0
not using I-node (on process 0) routines
5612. FXでのPETSc‐MPI版利用例(2/2)
ex23.c ： PETScのMPI版サンプルプログラム

13. PETSc‐複素数演算例題(MPI版1/2)
ex11.c ： PETScの複素数演算サンプルプログラム（MPI版）$ isub -mpi ”mpiexec ./ex23 –pc_type sor –ksp_type cgs”
$ cat OUTKSP Object: 4 MPI processes
type: cgs
maximum iterations=10000, initial guess is zerotolerances: relative=1e-07, absolute=1e-50, divergence=10000
left preconditioning
using PRECONDITIONED norm type for convergence test
PC Object: 4 MPI processes
type: sor
SOR: type = local_symmetric, iterations = 1, local iterations = 1, omega= 1linear system matrix = precond matrix:
Matrix Object:
4 MPI processes
type: mpiaij
rows=10, cols=10
total: nonzeros=28, allocated nonzeros=70
total number of mallocs used during MatSetValues calls =0
not using I-node (on process 0) routines

$ ls
OUT
ex1*

ex1.c
ex11.c@ex11_mod.c
ex11_org.c

ex1_mod.c
ex23*$ ls -l ex11.c
lrwxrwxrwx 1 t00003 gt00 10

ex23.cmakefile＃実習状況により，表示されるファイル群は異なり
ますが，少なくとも赤文字のファイルは必須．

2月 19 16:08 2013 ex11.c -> ex11_mod.c＃PETSc提供のオリジナル ex11.c はコードが不完全のため，$ module list
実行時にアベンドする．ex11_mod.cは伊藤が修正したもの．
Currently Loaded Modulefiles:
1) TCSuite/GM-1.2.1-04
2) PETSc/3.3-p4$ module unload PETSc$ module list
Currently Loaded Modulefiles:
1) TCSuite/GM-1.2.1-04

PETScで複素数演算を行うために，complex版に切替える．
実数演算でもcomplex版（PETScインストールの
cofigureにてcomplexオプションを有効にしている）
で良いのでは無いかと思うのだが，詳細は未確認．

$ module load PETSc-complex$ module list
Currently Loaded Modulefiles:
1) TCSuite/GM-1.2.1-04
2) PETSc-complex/3.3-p4$ make ex11
$ isub -mpi ”mpiexec ./ex11”57

13. PETSc‐複素数演算例題(MPI版2/2)

58

14. 色々な解法で解いてみよう

ex11.c ： PETScの複素数演算サンプルプログラム（MPI版）

解法
デフォルトはGMRES$ isub -mpi ”mpiexec ./ex11”$ cat OUT
1+1i
Norm of error 0.00050668 iterations 23
1+1i
1+1i
1+1i
$ isub -mpi ”mpiexec ./ex11 –pc_type sor –ksp_type cgs”
$ cat OUT
1+1i
Norm of error 8.56439e-05 iterations 28
1+1i
1+1i
1+1i

59

Options
Database Name（‐ksp_typeの引数）

KSPType

Richardson

KSPRICHARDSON

richardson

Chebychev

KSPCHEBYCHEV

chebychev

Conjugate Gradient

KSPCG

cg

BiConjugate Gradient

KSPBICG

bicg

Generalized Minimal Residual

KSPGMRES

gmres

BiCGSTABKSPBCGS

bcgs

Conjugate Gradient Squared

KSPCGS

cgs

Transpose‐Free Quasi‐Minimal Residual (1)KSPTFQMR

tfqmr

Transpose‐Free Quasi‐Minimal Residual (2)KSPTCQMR

tcqmr

Conjugate Residual

KSPCR

cr

Least Squares Method

KSPLSQR

lsqr

Shell for no KSP method

KSPPREONLY

preonly● 上記は3.0.0版のuser manualの一覧であるが，他の解法も利用可能（3.3-p4版のuser manual参照）． 60

14. 色々な前処理を組合せてみよう
前処理
デフォルトはJacobi

主なオプション

Options
Database Name（‐pc_typeの引数）

PCType14. 色々なオプションを試してみよう
内容[ ]内はデフォルト値‐ksp_max_it [100000]

最大反復回数

Jacobi

PCJACOBI

jacobi‐ksp_rtol [1.0e‐5]

収束判定基準（相対残差ノルム）

Block Jacobi

PCBJACOBI

bjacobi‐ksp_atol [1.0e‐50]

収束判定基準（絶対残差ノルム）

SOR (and SSOR)

PCSOR

sor

‐ksp_divtol [1.0e5]

発散状況の判定基準（相対残差ノルム）

SOR with Eisenstat trick

PCEISENSTAT

eisenstat‐ksp_right_pc [left]

アルゴリズムの前処理系（デフォルトは左前処理系）

Incomplete Cholesky

PCICC

icc

‐mat_type [AIJ形式]

Incomplete LU

PCILU

ilu

Additive Schwarz

PCASM

asm

Linear solver

PCKSP

ksp

Combination of preconditioners

PCCOMPOSITE

composite

LU

PCLU

lu

係数行列の格納形式‐mat_type seqaij‐mat_type mpiaij‐mat_type mpirowbs‐mat_type seqdense‐mat_type mpidense‐mat_type seqbaij‐mat_type mpibaij

Cholesky

PCCholeskycholesky

No preconditioning

PCNONEnone

inc main(inc argc, char **args ) {PetscInitialize
解ベクトル x を定義
x のサイズを設定

MatCreate(PETSC_COMM_WORLD,&A);
行列 A を定義
MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);
行列 Aのサイズを設定
MatSetFromOptions(A);-mat_typeオプションの反映（デフォルトはAIJ形式）
MatSetUp(A);このプログラム内で使えるよう行列 A のデータ構造をセットアップ
Assemble matrix */＃＃＃ 処理内容は２つ後のスライドを参照 ＃＃＃KSPCreate(PETSC_COMM_WORLD,&ksp);
使用するKSPを定義
KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);
行列Aを係数行列として設定し，前処理用行列を定義
KSPGetPC(ksp,&pc);
PCSetType(pc,PCJACOBI);
実際に使用する前処理演算を定義
KSPSetTolerances(ksp,1.e-5, PETSC_DEFAULT,PETSC_DEFAULT, PETSC_DEFAULT);
KSPSetFromOptions(ksp);
KSPSolve(ksp,b,x);
求解ルーチン
KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);
計算結果の出力
PetscFinalize();
return 0;}

rk 2  dtol  b 2

逐次版AIJ形式（CSR(CRS)形式）
MPI版AIJ形式（ 〃 ）MPI版rowbs形式
逐次版密行列MPI版密行列
逐次版block AIJ形式（ 〃 ）MPI版block AIJ形式（ 〃 ）62● 上記やその他のオプションの内容の詳細はPETScのuser manual参照．

15. PETScのex1で使用している係数行列と右辺項

15. ex1 処理フローの概要/*

PETScの魅力の一つとして，makeした実行モジュールに上記のようなオプションを付けることでコマンドラ
イン実行できることが挙げられるが，今回の実習ではバッチ実行のみに制限されているため，コマンドライ
ン記述をスクリプトに書き込む必要がある．

● 上記はuser manualの記述であるが，ソースコードを眺めると他の前処理も利用できそうな様子である． 61

VecCreate(PETSC_COMM_WORLD,&x);
VecSetSizes(x,PETSC_DECIDE,n);

rk 2  max rtol  b 2 , atol

Ax  bb  Ax, x  1.0,1.0, ,1.0

T0 2 1 1 2 1




A
  



1
2

1


 0
1 2 
係数行列Aの次数n=10（プログラム内で定義）の3重対角行列6364

15. ex1の係数行列作成の概要
PetscInt
PetscScalarn = 10, col[3];one = 1.0, value[3];15. ex1の解のファイル出力例(1/3)

要するに一旦下表を構成して，Aに値を設定している

i = 0; col[0] = 0; col[1] = 1; value[0] = 2.0; value[1] = -1.0;
MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

1i～

i = n - 1; col[0] = n - 2; col[1] = n - 1;
MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);

0

n-2
n-1

0

1

2.0
{0}

-1.0
{1}

-1.0
{i-1}

2.0
{i}

-1.0
{n-2}

2.0
{n-1}

$ ls
OUT ex1*
makefile

value [ j ]{ col [ j ] }/* Assemble matrix */value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
for (i=1; i<n-1; i++) {col[0] = i-1; col[1] = i; col[2] = i+1;
MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);}

ex1_mod.c ： PETScのex1を伊藤が改造したプログラム（ex1解ベクトル出力版）

2

$ mv
$ mv

ex1.c

ex1.c
ex1

ex11*

ex11.c@ex11_mod.c

ex11_org.c

ex1_mod.c

ex23*

ex23.c＃実習状況により，表示されるファイル群は異なり
ますが，少なくとも赤文字のファイルは必須．

ex1_org.c
ex1_org$ ln –s ex1_mod.c ex1.c
$ ls

-1.0
{i+1}

OUT
ex1.c@

ex11*
ex11.cex1_mod.c
ex1_org*

ex1_org.c
ex23*

ex23.cmakefile$ ls –l ex1.c
lrwxrwxrwx 1 t00003 gt00 9 12月 14 14:13 2012 ex1.c -> ex1_mod.c$ diff※ ex23(MPI並列版)では/* Set exact solution; then compute right-hand-side vector. */ 変数rstart, rendに各ノード
VecSet(u,one);Tのindexを格納して制御．
b  Au , u  1.0,1.0, ,1.0
MatMult(A,u,b);65

ex1_org.c

ex1_mod.c

PETSc も Lis も，このような要領にてサンプル
プログラムに手を加え，徐々に自分用のシミュレーションコードにアレンジしていくのが近道（なのだと思う）．

36a37,39> /* add - begin by itosho */>PetscViewer
view_out;> /* add - end
by itosho */144a148,154> /* add - begin by itosho */> //VecView(x, PETSC_VIEWER_STDOUT_WORLD); // To STDOUT file>
PetscViewerASCIIOpen(PETSC_COMM_WORLD,"Ex1Sol",&view_out);>
VecView(x, view_out);>
PetscViewerDestroy(&view_out);> /* add - end
by itosho */66

15. ex1の解のファイル出力例(2/3)

15. ex1の解のファイル出力例(3/3)$ module list

$ cat OUT

Currently Loaded Modulefiles:
1) TCSuite/GM-1.2.1-04
2) PETSc/3.3-p4

KSP Object: 1 MPI processes
type: gmres
GMRES: restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization withno iterative refinement
GMRES: happy breakdown tolerance 1e-30
maximum iterations=10000, initial guess is zerotolerances: relative=1e-05, absolute=1e-50, divergence=10000
left preconditioning
using PRECONDITIONED norm type for convergence test
PC Object: 1 MPI processes
type: jacobi
linear system matrix = precond matrix:
Matrix Object:
1 MPI processes
type: seqaij
rows=10, cols=10
total: nonzeros=28, allocated nonzeros=50
total number of mallocs used during MatSetValues calls =0
not using I-node routines$ make ex1
・・・ 省略 ・・・＃この例題では，PETSc-complex/3.3-p4であっても実行可能$ isub -seq ”mpiexec ./ex1”$ lrt
・・・ 省略 ・・・
lrwxrwxrwx 1 t00003 gt00
9 12月 14 14:13 2012 ex1.c -> ex1_mod.c-rwxr-xr-x 1 t00003 gt00 19289246 12月 14 14:22 2012 ex1*
-rw-r--r-- 1 t00003 gt00
71 12月 14 14:22 2012 Ex1Sol-rw-r--r-- 1 t00003 gt00
684 12月 14 14:22 2012 OUT$ cat Ex1Sol
Vector Object:Solution 1 MPI processes
type: seq1
1
1
11
参考：「15. PETScのex1で使用している係数行列と右辺項」
1
1
1
1
1
67

68

16. 線形方程式（Ax=b）求解のための反復解法とその収束性に関する評価の概要

16. KSP法に関する盲点

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

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

適当な初期値 x0

p0  r0  b  Ax0

 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

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

end

Ax  b 前処理行列： K  A 残差ベクトル： r  b  Ax
左前処理系（Left-preconditioned system）：

test1_org.c test1_mod.c との比較：

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

b 2 用の変数を宣言

55a56>
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); // additionbrˆ 2 b 2 を算出$ cat OUT

x0 , r0  K 1  b  Ax0  ,r0# , r0  0, e.g., r0#  r0 ,  1  0,
k  0,1, 2,; Do
uk  rk   k 1qk 1 ,

k 

#
0

1

r , K Apk,

qk  uk   k K 1 Apk ,・・・

xk 1  xk   k (uk  qk ),BiCG: rel_resid = 9.237507e-16

70

r̂ 2

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

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

rk 1  rk   k K A(uk  qk ),1

k 
End Do#0r , rk 1
r0# , rk,p k  K 1 pk , uk  K 1uk , qk  K 1qk ,rk  K 1rk , r0#  K T r0#

右側は両側・左・右３種類の前処理
方向の変換で，全て同じアルゴリズ
ムとなる（一致性）．収束判定は下式の通り

b  Axkb

x0 , r0  b  Ax0 ,r0# , r0  0, e.g ., r0#  r0 ,  1  0,k  0,1, 2,; Do
uk  rk   k 1qk 1 ,残差ベクトルpk  uk   k 1 (qk 1   k 1 pk 1 ),rk  b  Axk

残差ベクトル 1

rk  K

収束判定

rk

× b 

k 

 b  Axk 

K 1  b  Axk br0# , rk
r0# , AK 1 pk,

qk  uk   k AK 1 pk ,

rkK 1  b  Axk 

1

1△ K b

伊藤祥司，杉原正顯，姫野龍太郎，
情報処理学会論文誌ACS, 3(2), 2010．

従来版（Conventional）：

1

K Ax  K b○ b 

pk  uk   k 1 (qk 1   k 1 pk 1 ),r0# , rk

1

収束判定 rk/home/t00003/nfs/lis-1.2.49_seq/test

BiCG: residual = 6.399933e-15

反復回数

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

当スライドの事柄は，Lis‐1.2.95（2012/8/15）版での仕様変更に伴い，（2012/12/20以後の）Lis‐1.3.0を用いた講習会では説明不要となった．

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

アルゴリズム中の残差ベクトルを評価すれば，あらためて残差 rk 1  b  Axk 1 を算出せずに済む．

69

16. main関数での真の残差の算出rˆ  b  Axˆ ,log10  rbk 

相対残差ノルム

参考文献：$ diff test1_org.c test1_mod.c

収束判定

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

2. 残差ベクトルを用いた収束判定について
・Lis ： 概ね右前処理系に基づく・PETSc：デフォルトは左前処理系[2] 伊藤祥司，杉原正顯 ，クリロフ部分空間法に対する前処理方向とライブラリ実装における注意点，
日本応用数理学会2010年度年会，東京，9月，2010年．

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

K b

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

xk 1  xk   k K 1 (uk  qk ),rk 1  rk   k AK 1 (uk  qk ),k 
End Dor0# , rk 1
r0# , rk,72

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

16.前処理方向に関する
従来からの定義

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

Ax  b

ただし， orthomin(k)(ID:8)のみ左前処理系でありながら，別の残差ベクトルを用いて技巧的に収束判定している（好ましくない
方法だが，元のSLAPライブラリの設計を踏襲している）．

前処理行列：

左前処理（ l ）1

1

K AxK b
KL  KKR  Ib   と甘い場合があるかも知れないので注意．

16.クリロフ部分空間(KSP)の基本的な議論





1

1

1

xk  x0  c0 K r0  c1 K A K r0    ck K A
0

k

k

k

0

0

0

1

0

k 1

0

~
A ： 前処理変換された係数行列~r0 ： 前処理変換された残差ベクトル

KL ：左前処理行列

KL  IKR ：右前処理行列

I ：単位行列KR  K○残差ベクトルは本来の情報を保持している．

k ( AK 1 ; r0 ) 

0

2

k

AK 1 K x  b

Kxk 1  Kxk  (b  Axk )    K r  x  z , z    A ; r  .
  , A r , , A r  KSP法の反復の度に  A ; r   Span  r , Ar
空間が拡張されていく1

K L1 AK R1 K R x  K L1b

AK 1 K x  b

前処理変換された残差ベクトルk 1

右前処理（ r ）16.右前処理系の例

A  K  R, Kxk 1  Rxk  b, k  0,1,2,
r0  b  Ax0 を用いると，両側前処理（ b ）[1] Saad, Y., Iterative Methods for Sparse Linear Systems (2nd ed.), SIAM, 2003.
[2] Van der Vorst, H. A., Iterative Krylov Methosd for Large Linear Systems, Cambridge, 2003.
[3] Barrett, R., et al., Templates for the Solution of Linear Systems: Building Blocks for Iterative
Methods, (2nd ed.), SIAM, 1994.
[4] 藤野清次，張 紹良，反復法の数理，朝倉書店，1996.
74※ 定常反復法（jacobi(ID:10)～sor(ID:12)）と非定常反復法とでは論点が異なるので注意．
73

xk 1  K 1 ( Rxk  b)  xk  K 1 (b  Axk ).

（残差ベクトル： r  b  Ax ）K KLKR

bicgsafe(ID:13)～minres(ID:22)については，伊藤の方では未確認のため，左前処理系にも関わらず収束判定が

rk

KA

0

Dongarra, Duff, Sorensen, Van der Vorst,
Solving Linear Systems on Vector and Shared
Memory Computers, SIAM, 1991.
75

rk  b  AK 1Kxk
実質的な残差ベクトル：

Span r0 , ( AK 1 )r0 , ( AK 1 ) 2 r0 ,  , ( AK 1 ) k 1 r0

r  b  AK 1Kx

Kxk  Kx0  c0 r0  c1 ( AK 1 )r0  c2 ( AK 1 ) 2 r0    ck 1 ( AK 1 ) k 1 r0
 Kxk  Kx0  {c0 r0  c1 ( AK 1 )r0  c2 ( AK 1 ) 2 r0    ck 1 ( AK 1 ) k 1 r0 }

 Kx0  k AK 1 ; r0

△生成されるKSPに対し得られる解は，前処理行列が作用したものに基づく．
7616.左前処理系の例
K 1 A x  K 1b
xk 1  xk  ( K 1b  K 1 Axk )1k ( K A; r0 ) 

16.右・左前処理系の例のまとめ

×残差ベクトルは本来の情報を保持していない．

前処理系

rk  K 1b  K 1 Axk

Lisに収められている内の主なKSP法（前半のID
右系 の解法）．その他多くのライブラリやアルゴリズ
ム記述もこちらの系

実質的な残差ベクトル：

Span r0 , ( K 1 A)r0 , ( K 1 A) 2 r0 ,  , ( K 1 A) k 1 r0

1r  K (b  Ax )



左系 PETScのデフォルトxk  x0  {c0 r0  c1 ( K 1 A)r0  c2 ( K 1 A) 2 r0    ck 1 ( K 1 A) k 1 r0 }





77

17. 両ライブラリの習得要領● 係数行列をMatrix Market形式で用意し，test1プログラムを活用（test1自体が単独の数値計算ソフトウェアを構成しているとも考えられる）● test4，test5の記述を参考にして各自のプログラムを開発する
test4 : 係数行列Aの全要素に値を代入する例
test5 : 係数行列Aを直接CRS形式で記述する例

PETSc：● src/ksp/ksp/examples/tutorials{ex1,ex23,ex11}等の記述を参考にして各自のプログラムを開発する（PETScはこの要領で独習するスタイル）
ex1 : 逐次版（MPIコードの1ノード用）：係数行列は実行列
ex23 : MPI並列（ex1の並列版）
ex11 : MPI並列：係数行列は複素行列● 係数行列のファイルをMatrix Market形式で利用するプログラムは，
src/mat/examples/tests/{ex72.c,ex78.c,ex32.c}等を参照して各自で開発する．
もしくは，外部ツールの活用もあるが，pestcのバージョンに対応していないと動作しないものも少なく無い．
● PETScの文法に関しては下記ページ参照
http://www.mcs.anl.gov/petsc/documentation/生成されるKSPと解の関係○保持している△前処理行列が作用したものに基づく

rk  b  AK 1Kxk

Kxk  Kx0  k AK 1 ; r0

×保持していない○前処理行列が作用しないものxk  x0  k K 1 A; r0

rk  K 1b  K 1 Axk









このような解析結果や解決策に関する発表：○生成されるKSPに対し得られる解は，前処理行列が作用しないものである．

Lis：

アルゴリズム中の残差ベクトルの情報

実のところ，前処理付きアルゴリズムの構造上，非常に重要なポイントであるにも関わらず，各々の前処理系にてそこそこの効果が得られるため盲点となっており，
もっぱら，前処理演算開発の方に注力されてきているのが現状である．

 xk  x0  {c0 r0  c1 ( K 1 A)r0  c2 ( K 1 A) 2 r0    ck 1 ( K 1 A) k 1 r0 }
 x0  k K 1 A; r0

ライブラリ

79

[1] 伊藤祥司，杉原正顯，クリロフ部分空間法の前処理系に対する新しい解釈，日本応用数理学会「行列・固有値問題の解法とその応用」研究部会，NII，2009年11月．
[2] 伊藤祥司，杉原正顯，双ランチョス系統の前処理付きアルゴリズムの改善，計算工学講演会論文集
Vol.15，pp.171-174，2010年5月．
[3] 伊藤祥司，杉原正顯，様々な前処理付きCGS法に対するクリロフ部分空間に注目した解析，
日本応用数理学会 2013年 研究部会 連合発表会，3月14日，東洋大学（白山），2013年．

78