DD-AVX User s Manual (based on Lis ) Toshiaki Hishinuma a, Akihiro Fujii a, Teruo Tanaka a, Hidehiko Hasegawa b

(1)

DD-AVX User’s Manual

バージョン 1.0.0 (based on Lis 1.4.58)

Toshiaki Hishinuma

a

_{, Akihiro Fujii}

a

_{, Teruo Tanaka}

a

_{, Hidehiko Hasegawa}

b

http://www.slis.tsukuba.ac.jp/∼s1530534/DD-AVX.html a

Kogakuin University,bUniversity of Tsukuba

copyright c⃝Toshiaki Hishinuma, Akihiro Fujii, Teruo Tanaka, Hidehiko Hasegawa

(2)

1 はじめに 2 2 導入 3 2.1 システム要件 . . . . 3 2.2 ライブラリのダウンロード・展開 . . . . 3 2.3 DD-AVXと Lis のマージ . . . . 4 2.4 Configure . . . . 4 2.5 実行ファイルの生成・導入 . . . . 5 2.5.1 実行ファイルの生成 . . . . 5 2.5.2 導入 . . . . 9 2.6 検証 . . . . 10 2.6.1 DD-AVX test1 . . . . 10 2.6.2 test1 . . . . 10 2.6.3 test2 . . . . 11 2.6.4 test2b . . . . 11 2.6.5 test3 . . . . 11 2.6.6 test3b . . . . 11 2.6.7 test4 . . . . 12 2.6.8 test5 . . . . 12 2.6.9 etest1 . . . . 12 2.6.10 etest2 . . . . 12 2.6.11 etest3 . . . . 13 2.6.12 etest4 . . . . 13 2.6.13 etest5 . . . . 13 2.6.14 etest6 . . . . 13 2.6.15 spmvtest1 . . . . 14 2.6.16 spmvtest2 . . . . 14 2.6.17 spmvtest2b . . . . 14 2.6.18 spmvtest3 . . . . 15 2.6.19 spmvtest3b . . . . 15 2.6.20 spmvtest4 . . . . 15 2.6.21 spmvtest5 . . . . 16 2.7 制限事項 . . . . 16 3 基本操作 17 3.1 初期化・終了処理 . . . . 17 3.2 ベクトルの操作 . . . . 18 3.3 行列の操作 . . . . 21 3.4 線型方程式の求解 . . . . 27 3.5 固有値問題の求解 . . . . 32 3.6 プログラムの作成 . . . . 34

(3)

3.7 実行ファイルの生成 . . . . 37 3.8 実行 . . . . 39 4 4倍精度演算 41 4.1 4倍精度演算の使用 . . . . 41 4.2 DD-AVXの追加機能 . . . . 43 4.2.1 BSR形式を用いた高速な疎行列ベクトル積 . . . . 43 4.2.2 疎行列ベクトル積における OpenMP のスケジューリング方式自動チューニング . . 43 5 行列格納形式 45 5.1 Compressed Sparse Row (CSR) . . . . 45

5.1.1 行列の作成 (逐次・マルチスレッド環境) . . . . 45

5.1.2 行列の作成 (マルチプロセス環境) . . . . 46

5.1.3 関連する関数 . . . . 46

5.2 Compressed Sparse Column (CSC) . . . . 47

5.2.2 行列の作成 (マルチプロセス環境) . . . . 48

5.2.3 関連する関数 . . . . 48

5.3 Modified Compressed Sparse Row (MSR) . . . . 49

5.3.1 行列の作成 (逐次・マルチスレッド環境) . . . . 49 5.3.2 行列の作成 (マルチプロセス環境) . . . . 50 5.3.3 関連する関数 . . . . 50 5.4 Diagonal (DIA) . . . . 51 5.4.1 行列の作成 (逐次環境) . . . . 51 5.4.2 行列の作成 (マルチスレッド環境) . . . . 52 5.4.3 行列の作成 (マルチプロセス環境) . . . . 53 5.4.4 関連する関数 . . . . 53

5.5 Ellpack-Itpack Generalized Diagonal (ELL) . . . . 54

5.5.2 行列の作成 (マルチプロセス環境) . . . . 55

5.5.3 関連する関数 . . . . 55

5.6 Jagged Diagonal (JAD) . . . . 56

5.6.1 行列の作成 (逐次環境) . . . . 57

5.6.2 行列の作成 (マルチスレッド環境) . . . . 58

5.6.3 行列の作成 (マルチプロセス環境) . . . . 59

5.6.4 関連する関数 . . . . 59

5.7 Block Sparse Row (BSR) . . . . 61

5.7.2 行列の作成 (マルチプロセス環境) . . . . 62

5.7.3 関連する関数 . . . . 62

5.8 Block Sparse Column (BSC) . . . . 63

(4)

5.8.3 関連する関数 . . . . 64

5.9 Variable Block Row (VBR) . . . . 65

5.9.1 行列の作成 (逐次・マルチスレッド環境) . . . . 66 5.9.2 行列の作成 (マルチプロセス環境) . . . . 67 5.9.3 関連する関数 . . . . 68 5.10 Coordinate (COO) . . . . 69 5.10.1 行列の作成 (逐次・マルチスレッド環境) . . . . 69 5.10.2 行列の作成 (マルチプロセス環境) . . . . 70 5.10.3 関連する関数 . . . . 70 5.11 Dense (DNS) . . . . 71 5.11.1 行列の作成 (逐次・マルチスレッド環境) . . . . 71 5.11.2 行列の作成 (マルチプロセス環境) . . . . 72 5.11.3 関連する関数 . . . . 72 参考文献 73 A ファイル形式 79 A.1 拡張 Matrix Market 形式 . . . . 79

A.2 Harwell-Boeing形式 . . . . 80

A.3 ベクトル用拡張 Matrix Market 形式 . . . . 81

(5)

1 はじめに

“DD-AVX”は，高速な倍々精度反復解法に向けた倍々精度演算ライブラリである．反復解法ライブラリ Lis[1](http://www.ssisc.org/lis/) と組み合わせて使用することを想定している，Lis に含まれる， • 倍々精度ベクトル演算 • 倍精度疎行列と倍々精度ベクトルの積 (転置含む)

を SIMD 拡張命令 AVX, AVX2 を用いて高速化したルーチン群 [86] からなり，Lis とマージすることで AVX,

AVX2を用いた倍々精度反復解法を解くことを想定している．また，我々が明らかにした以下の 3 つの倍々精度の高速化技法も追加機能として実装した． • AVX に適した疎行列の格納形式 BSR[87](AVX 対応は 4x1 のみ) • OpenMP のスレッド負荷分散スケジューリング方式 balanced[88] • OpenMP のスケジューリング方式の自動選択 [88] これらは．2.6.1 DD-AVX test1 と，4 章の 4 倍精度演算に利用方法を記載した．

また，各関数の使い方は別紙の Function Reference を参照のこと．DD-AVX によって追加された 2 つの関数は，Funtion Reference の 4.9 節に記載した．

本ライブラリはマージ後の使い方は Lis と同じであるように設計したため，Lis を用いて作られたプログラムは，変更せずに利用することが可能である．

(6)

2 導入

2.1 システム要件

本章では，AVX，AVX2 を用いた動作確認環境について述べる．

DD-AVXの導入には C コンパイラが必要である. また，DD-AVX は，UNIX, Linux オペレーティングシ

ステム以外は動作未確認である.

並列計算環境では, OpenMP ライブラリまたは MPI ライブラリを使用する.

SIMD拡張命令 AVX，AVX2 に対応しているかは，導入環境の “/proc/cpuinfo” およびコンパイラの対応

状況の確認が必要である．

表 1 に DD-AVX の AVX での動作確認環境を示す．

表 1: DD-AVX の動作確認環境 (AVX)

CPU (Code name) C Compiler OS

Intel Core i7 2600 K (Sandy Bridge) Intel C/C++ Compiler 12.0.3, 13.1.3, 15.0.0 Fedora 16 gcc 4.6.3, 4.8.3

Intel Core i7 3770 K (Ivy Bridge) Intel C/C++ Compiler 12.0.3, 13.1.3, 15.0.0 Cent OS 6.4 gcc 4.6.3, 4.8.3

Intel Core i7 4770 (Haswell) Intel C/C++ Compiler 12.0.3, 13.1.3, 15.0.0 Cent OS 6.4 gcc 4.6.3, 4.8.3

Intel Core i7 4690 S (Haswell Refresh) Intel C/C++ Compiler 12.0.3, 13.1.3, 15.0.0 Cygwin on Windows 7 gcc 4.6.3, 4.8.3

表 2 に DD-AVX の AVX2 での動作確認環境を示す．

表 2: DD-AVX の動作確認環境 (AVX2)

CPU (Code name) C Compiler OS

Intel Core i7 4770 (Haswell) Intel C/C++ Compiler 12.0.3, 13.1.3, 15.0.0 Cent OS 6.4 gcc 4.8.3

Intel Core i7 4690 S (Haswell Refresh) Intel C/C++ Compiler 12.0.3, 13.1.3, 15.0.0 Cygwin on Windows 7 gcc 4.8.3

2.2 ライブラリのダウンロード・展開

本章では，Lis, DD-AVX を Web ページからダウンロードし，展開する方法を述べる．

本ライブラリは，反復解法ライブラリ Lis(ver.1.4.58) に対し，DD-AVX ライブラリをマージして使用することを想定している．

はじめに，Lis と DD-AVX を同一のフォルダにダウンロードする．ダウンロードをコマンドライン上において行うためには，

(7)

• Lis1.4.58 のダウンロード > wget http://www.ssisc.org/lis/dl/lis-1.4.58.tar.gz • DD-AVX のダウンロード > wget http://www.slis.tsukuba.ac.jp/~s1530534/DD-AVX-files/DD-AVX1.2.0.tar.gz を用いる．次に，これらのアーカイブの展開を行うために，コマンドライン上において， • Lis-1.4.58 の展開

> gunzip -c lis-1.4.58.tar.gz | tar -xv

• DD-AVX の展開

> gunzip -c DD-AVX1.0.0.tar.gz | tar -xv

を行う．これにより，以下の 3 つが同様のフォルダに展開される． 1. DD-AVX/ 2. Lis-1.4.58/ 3. DD-AVX-merge.sh また，これ以降 tar.gz ファイルは必要ないので，不要であれば削除しても構わない．ディレクトリ DD-AVX は，図 1 に示す DD-AVX のサブディレクトリから成る．

2.3 DD-AVX

と Lis のマージ

本章では，展開した lis-1.4.58 と DD-AVX のマージ方法について述べる.

マージには，DD-AVX-merge.sh を用いる．このスクリプトの引数は [lis のディレクトリ], [DD-AVX のディレクトリ], [マージ後の出力ディレクトリ] である．

以下のコマンドを実行し，Lis と DD-AVX をマージする． > sh ./DD-AVX-merge.sh lis-1.4.58/ [output_dir]

これにより, lis-1.4.58 と DD-AVX をマージした [output dir] が得られる.

2.4 Configure

Configureと，そのオプションについて述べる．

“–enable-avx”，“–enable-avx2” が DD-AVX から新たに追加されたオプションである．

ディレクトリ [output dir] において次のコマンドを実行し, Configure を行う．これにより，ユーザが定義した環境向けの Makefile が生成される．

(8)

> ./configure --prefix=<install-dir> [option]とする．

例えば，OpenMP と MPI のハイブリッド並列で 4 倍精度反復解法を AVX2 を用いて，mpicc のオプションに-static intel を付けたいとき，以下のようにする．

> ./configure --prefix=<install-dir> --enable-omp --enable-mpi --enable-quad --enable-avx2 CFLAGS=-static\ intel

本ライブラリを使うためのオプションを表 3 に示す．表 4 に TARGET として指定できる主な計算機環境を示す.

なお，SSE2, AVX, AVX2 は倍々精度精度演算 (“–enable quad”) でしか利用できない．SIMD の設定は， Configure時に設定し，以降は再度 Configure するまで変更できない．

2.5 実行ファイルの生成・導入

本章では，実行ファイルの生成を行い，ライブラリを導入する方法を述べる． 2.5.1 実行ファイルの生成ディレクトリ [output dir] において次のコマンドを入力し, 実行ファイルを生成する. > make 実行ファイルが正常に生成されたかどうかを確認するには, のコマンドを入力し, ディレクトリ [putput dir]/test

(9)

表 3: Configure 時のオプション一覧 ( ./configure --help からも参照可能) --enable-omp OpenMPライブラリを使用 --enable-mpi MPIライブラリを使用 --enable-fortran FORTRAN 77互換インタフェースを使用 --enable-f90 Fortran 90互換インタフェースを使用 --enable-saamg SA-AMG前処理を使用 --enable-quad 倍々精度演算を使用

--enable-longdouble long double型 4 倍精度演算を使用

--enable-longlong 64ビット整数型を使用

--enable-sse2 SIMD拡張命令 SSE2 を使用

--enable-avx SIMD拡張命令 AVX を使用

--enable-avx2 SIMD拡張命令 AVX2 を使用

--enable-debug デバッグモードを使用 --enable-shared 動的リンクを使用 --enable-gprof プロファイラを使用 --prefix=<install-dir> 導入先を指定 TARGET=<target> 計算機環境を指定 CC=<c_compiler> Cコンパイラを指定 CFLAGS=<c_flags> Cコンパイラオプションを指定 F77=<f77_compiler> FORTRAN 77コンパイラを指定 F77FLAGS=<f77_flags> FORTRAN 77コンパイラオプションを指定 FC=<f90_compiler> Fortran 90コンパイラを指定 FCFLAGS=<f90_flags> Fortran90コンパイラオプションを指定 LDFLAGS=<ld_flags> リンクオプションを指定

(10)

表 4: TARGET の一覧 (詳細は lis-($VERSION)/configure.in を参照)

<target> 等価なオプション

cray_xt3_cross ./configure CC=cc FC=ftn CFLAGS="-O3 -B -fastsse -tp k8-64"

FCFLAGS="-O3 -fastsse -tp k8-64 -Mpreprocess" FCLDFLAGS="-Mnomain" ac_cv_sizeof_void_p=8 cross_compiling=yes

ax_f77_mangling="lower case, no underscore, extra underscore" fujitsu_fx10_cross ./configure 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" hitachi_sr16k ./configure CC=cc FC=f90 CFLAGS="-Os -noparallel"

FCFLAGS="-Oss -noparallel" FCLDFLAGS="-lf90s" ac_cv_sizeof_void_p=8

ax_f77_mangling="lower case, underscore, no extra underscore" ibm_bgl_cross ./configure CC=blrts_xlc FC=blrts_xlf90

CFLAGS="-O3 -qarch=440d -qtune=440 -qstrict"

FCFLAGS="-O3 -qarch=440d -qtune=440 -qsuffix=cpp=F90" ac_cv_sizeof_void_p=4 cross_compiling=yes

ax_f77_mangling="lower case, no underscore, no extra underscore" nec_sx9_cross ./configure CC=sxmpic++ FC=sxmpif90 AR=sxar RANLIB=true

ac_cv_sizeof_void_p=8 ax_vector_machine=yes cross_compiling=yes ax_f77_mangling="lower case, no underscore, extra underscore" DD-AVX_cross ./configure CC=icc --enable-omp --enable-mpi --enable-avx DD-AVX2_cross ./configure CC=icc --enable-omp --enable-mpi --enable-avx2

DD-AVX_MPI_cross ./configure CC=icc --enable-omp --enable-mpi --enable-avx --enable-mpi DD-AVX2_MPI_cross ./configure CC=icc --enable-omp --enable-mpi --enable-avx2 --enable-mpi

(11)

に生成された実行ファイルを用いて検証を行う. > make check

このコマンドでは, Matrix Market 形式のファイル test/testmat.mtx から行列, ベクトルデータを読み込 み, BiCG 法を用いて線型方程式 Ax = b の解を求める. 以下に SGI Altix 3700 上での実行結果を示す. な おオプション--enable-omp と--enable-mpi は組み合わせて使用することができる.

既定

matrix size = 100 x 100 (460 nonzero entries) initial vector x = 0

precision : double solver : BiCG 2 precon : none

conv_cond : ||b-Ax||_2 <= 1.0e-12 * ||b-Ax_0||_2 storage : CSR

lis_solve : normal end

BiCG: number of iterations = 15 (double = 15, quad = 0) BiCG: elapsed time = 5.178690e-03 sec.

BiCG: preconditioner = 1.277685e-03 sec. BiCG: matrix creation = 1.254797e-03 sec. BiCG: linear solver = 3.901005e-03 sec. BiCG: relative residual = 6.327297e-15

--enable-omp

max number of threads = 32 number of threads = 2

(12)

--enable-mpi

number of processes = 2

2.5.2 導入ディレクトリ [output dir] において次のコマンドを入力し, 導入先のディレクトリにファイルを複製する. > make install これにより, ディレクトリ<install-dir>に以下のファイルが複製される. ($INSTALLDIR) +bin

| +lsolve esolve hpcg_kernel hpcg_spmvtest spmvtest* +include

| +lis_config.h lis.h lisf.h +lib

| +liblis.a +share

+doc/lis examples/lis man

lis config.hはライブラリを生成する際に, また lis.h は C, lisf.h は Fortran でライブラリを使用する際に必要なヘッダファイルである. liblis.a は生成されたライブラリファイルである. ライブラリファイルが正常に導入されたかどうかを確認するには, ディレクトリ lis-($VERSION) において次のコマンドを入力し, ディレクトリ examples/lis に生成された実行ファイルを用いて検証を行う.

> make installcheck

examples/lis下の test1, etest5, test3b, spmvtest3b は, lsolve, esolve, hpcg kernel, hpcg spmvtest の別名で ($INSTALLDIR)/bin に複製される. examples/lis/spmvtest*も, それぞれ ($INSTALLDIR)/bin に複製される.

(13)

2.6 検証

2.6.1 DD-AVX test1

サンプルプログラムは，[output dir]/test にある．

DD-AVX test1.cは，DD-AVX から追加されたサンプルプログラムで，Lis の test1.c を DD-AVX 向けにしたものである．このプログラムは，4 倍精度での利用を前提としている．

実行するには，ディレクトリ [output dir]/test において

> ./DD-AVX_test1 matrix_filename rhs_setting solution_filename rhistory_filename [options] とすることで，matrix filename から行列データを読み込み, 線型方程式 Ax = b を options で指定され

た解法で解く. このとき，格納形式はブロックサイズ 4x1 の BCRS 形式になる．

また，コメントアウトを外せば，OpenMP のスケジューリング方式を自動で選択する機能も使用できる．解を拡張 Matrix Market 形式で solution filename に, 残差履歴を PLAIN 形式で rhistory filename に書き出す (付録 A を参照). 入力可能な行列データ形式は Matrix Market 形式, 拡張 Matrix Market 形式である. rhs setting には

0 行列データファイルに含まれる右辺ベクトルを用いる

1 b = (1, . . . , 1)T _を用いる

2 b = A× (1, . . . , 1)T _を用いる

rhs filename 右辺ベクトルのファイル名

のいずれかを指定できる. rhs filename は PLAIN 形式, Matrix Market 形式に対応する. test1f.F は test1.cの Fortran 版である.

2.6.2 test1

ディレクトリ [output dir]/test において

> ./test1 matrix_filename rhs_setting solution_filename rhistory_filename [options] と入力すると, matrix filename から行列データを読み込み, 線型方程式 Ax = b を options で指定され た解法で解く. また, 解を拡張 Matrix Market 形式で solution filename に, 残差履歴を PLAIN 形式で rhistory filenameに書き出す (付録 A を参照). 入力可能な行列データ形式は Matrix Market 形式, 拡張 Matrix Market形式である. rhs setting には

0 行列データファイルに含まれる右辺ベクトルを用いる

1 b = (1, . . . , 1)T _を用いる

2 b = A× (1, . . . , 1)T _を用いる

rhs filename 右辺ベクトルのファイル名

のいずれかを指定できる. rhs filename は PLAIN 形式, Matrix Market 形式に対応する. test1f.F は test1.cの Fortran 版である.

(14)

2.6.3 test2

> test2 m n matrix_type solution_filename rhistory_filename [options]

と入力すると, 2 次元 Poisson 方程式を 5 点中心差分で離散化して得られる次数 mn の 5 重対角行列を係数と する線型方程式 Ax = b を, matrix_type で指定された行列格納形式, options で指定された解法で解く. ま た, 解を拡張 Matrix Market 形式で solution filename に, 残差履歴を PLAIN 形式で rhistory filename に書き出す. ただし, 右辺ベクトル b は線型方程式 Ax = b の解ベクトルの値がすべて 1 となるよう設定さ れる. m, n は各次元の格子点数である. test2f.F90 は test2.c の Fortran90 版である.

2.6.4 test2b

> test2b m n matrix_type solution_filename rhistory_filename [options]

と入力すると, 2 次元 Poisson 方程式を 9 点中心差分で離散化して得られる次数 mn の 9 重対角行列を係数と する線型方程式 Ax = b を, matrix_type で指定された行列格納形式, options で指定された解法で解く. ま た, 解を拡張 Matrix Market 形式で solution filename に, 残差履歴を PLAIN 形式で rhistory filename に書き出す. ただし, 右辺ベクトル b は線型方程式 Ax = b の解ベクトルの値がすべて 1 となるよう設定さ れる. m, n は各次元の格子点数である.

2.6.5 test3

> test3 l m n matrix_type solution_filename rhistory_filename [options]

と入力すると, 3 次元 Poisson 方程式を 7 点中心差分で離散化して得られる次数 lmn の 7 重対角行列を係数と する線型方程式 Ax = b を, matrix_type で指定された行列格納形式, options で指定された解法で解く. ま た, 解を拡張 Matrix Market 形式で solution filename に, 残差履歴を PLAIN 形式で rhistory filename に書き出す. ただし, 右辺ベクトル b は線型方程式 Ax = b の解ベクトルの値がすべて 1 となるよう設定さ れる. l, m, n は各次元の格子点数である.

2.6.6 test3b

> test3b l m n matrix_type solution_filename rhistory_filename [options]

と入力すると, 3 次元 Poisson 方程式を 27 点中心差分で離散化して得られる次数 lmn の 27 重対角行列を係数 とする線型方程式 Ax = b を, matrix_type で指定された行列格納形式, options で指定された解法で解く. ま た, 解を拡張 Matrix Market 形式で solution filename に, 残差履歴を PLAIN 形式で rhistory filename に書き出す. ただし, 右辺ベクトル b は線型方程式 Ax = b の解ベクトルの値がすべて 1 となるよう設定さ れる. l, m, n は各次元の格子点数である.

(15)

2.6.7 test4 線型方程式 Ax = b を指定された解法で解き, 解を標準出力に書き出す. 行列 A は次数 12 の 3 重対角行列 A =          2 −1 −1 2 −1 . ._. . ._. . ._. −1 2 −1 −1 2         

である. 右辺ベクトル b は解 x がすべて 1 となるよう設定される. test4f.F は test4.c の Fortran 版で ある.

2.6.8 test5

ディレクトリ [output dir]/test において > test5 n gamma [options]

と入力すると, 線型方程式 Ax = b を指定された解法で解く. 行列 A は次数 n の Toepliz 行列 A =            2 1 0 2 1 γ 0 2 1 . ._. . ._. . ._. . ._. γ 0 2 1 γ 0 2            である. 右辺ベクトル b は解 x がすべて 1 となるよう設定される. 2.6.9 etest1 ディレクトリ [output dir]/test において

> etest1 matrix_filename evector_filename rhistory_filename [options]

と入力すると, matrix filename から行列データを読み込み, 固有値問題 Ax = λx を options で指定さ れた解法で解いて, 指定された固有値を標準出力に書き出す. また, 対応する固有ベクトルを拡張 Matrix Market形式で evector filename に, 残差履歴を PLAIN 形式で rhistory filename に書き出す. 入力可能な行列データ形式は Matrix Market 形式である. etest1f.F は etest1.c の Fortran 版である.

2.6.10 etest2

> etest2 m n matrix_type evector_filename rhistory_filename [options]

と入力すると, 2 次元 Helmholtz 方程式を 5 点中心差分で離散化して得られる次数 mn の 5 重対角行列に関 する固有値問題 Ax = λx を, matrix_type で指定された行列格納形式, options で指定された解法で解き, 指定された固有値を標準出力に書き出す. また, 対応する固有ベクトルを evector filename に, 残差履歴を rhistory filename に書き出す. m, n は各次元の格子点数である.

(16)

2.6.11 etest3

> etest3 l m n matrix_type evector_filename rhistory_filename [options]

と入力すると, 3 次元 Helmholtz 方程式を 7 点中心差分で離散化して得られる次数 lmn の 7 重対角行列に 関する固有値問題 Ax = λx を, matrix_type で指定された行列格納形式, options で指定された解法で解 き, 指定された固有値を標準出力に書き出す. また, 対応する固有ベクトルを拡張 Matrix Market 形式で evector filenameに, 残差履歴を PLAIN 形式で rhistory filename に書き出す. l, m, n は各次元の格子点数である. 2.6.12 etest4 ディレクトリ [output dir]/test において > etest4 n [options] と入力すると, 固有値問題 Ax = λx を指定された解法で解き, 指定された固有値を標準出力に書き出す. 行 列 A は次数 n の 3 重対角行列 A =          2 −1 −1 2 −1 . ._. . ._. . ._. −1 2 −1 −1 2          である. etest4f.F は etest4.c の Fortran 版である.

2.6.13 etest5

> etest5 matrix_filename evalues_filename evectors_filename residuals_filename iters_filename [options]

と入力すると, matrix filename から行列データを読み込み, 固有値問題 Ax = λx を options で指定された解 法で解いて, 指定された固有値を標準出力に書き出す. また, 指定された個数の固有値を evalues filename に, 対応する固有ベクトル, 残差ノルム及び反復回数を evectors filename, residuals filename 及び iters filenameに拡張 Matrix Market 形式で書き出す. 入力可能な行列データ形式は Matrix Market 形式である.

2.6.14 etest6

> etest6 l m n matrix_type evalues_filename evectors_filename residuals_filename iters_filename [options]

と入力すると, 3 次元 Helmholtz 方程式を 7 点中心差分で離散化して得られる次数 lmn の 7 重対角行列に関す る固有値問題 Ax = λx を, matrix_type で指定された行列格納形式, options で指定された解法で解き, 指 定された固有値を標準出力に書き出す. また, 指定された個数の固有値を evalues filename に, 対応する固

(17)

有ベクトル, 残差ノルム及び反復回数を evectors filename, residuals filename 及び iters filename に拡張 Matrix Market 形式で書き出す. l, m, n は各次元の格子点数である.

2.6.15 spmvtest1

ディレクトリ [output dir]/test において > spmvtest1 n iter [matrix_type]

と入力すると, 1 次元 Laplace 演算子を 3 点中心差分で離散化して得られる次数 n の 3 重対角係数行列 A =          2 −1 −1 2 −1 . ._. . ._. . ._. −1 2 −1 −1 2         

と (1, . . . , 1)T _{との積を iter で指定された回数実行し, FLOPS 値を算出する. 必要なら matrix type に}

より, 0 実行可能なすべての行列格納形式について測定する 1-11 行列格納形式の番号のいずれかを指定する. 2.6.16 spmvtest2 ディレクトリ [output dir]/test において > spmvtest2 m n iter [matrix_type]

と入力すると, 2 次元 Laplace 演算子を 5 点中心差分で離散化して得られる次数 mn の 5 重対角係数行列とベ クトル (1, . . . , 1)T _{との積を iter で指定された回数実行し, FLOPS 値を算出する. 必要なら matrix type}

により, 0 実行可能なすべての行列格納形式について測定する 1-11 行列格納形式の番号のいずれかを指定する. m, n は各次元の格子点数である. 2.6.17 spmvtest2b ディレクトリ [output dir]/test において > spmvtest2b m n iter [matrix_type]

と入力すると, 2 次元 Laplace 演算子を 9 点中心差分で離散化して得られる次数 mn の 9 重対角係数行列とベ クトル (1, . . . , 1)T _{との積を iter で指定された回数実行し, FLOPS 値を算出する. 必要なら matrix type}

により,

(18)

1-11 行列格納形式の番号のいずれかを指定する. m, n は各次元の格子点数である.

2.6.18 spmvtest3

> spmvtest3 l m n iter [matrix_type]

と入力すると, 3 次元 Laplace 演算子を 7 点中心差分で離散化して得られる次数 lmn の 7 重対角係数行列と ベクトル (1, . . . , 1)T _{との積を iter で指定された回数実行し, FLOPS 値を算出する. 必要なら matrix type}

により, 0 実行可能なすべての行列格納形式について測定する 1-11 行列格納形式の番号のいずれかを指定する. l, m, n は各次元の格子点数である. 2.6.19 spmvtest3b ディレクトリ [output dir]/test において

> spmvtest3b l m n iter [matrix_type]

と入力すると, 3 次元 Laplace 演算子を 27 点中心差分で離散化して得られる次数 lmn の 27 重対角係数行列と ベクトル (1, . . . , 1)T _{との積を iter で指定された回数実行し, FLOPS 値を算出する. 必要なら matrix type}

により, 0 実行可能なすべての行列格納形式について測定する 1-11 行列格納形式の番号のいずれかを指定する. l, m, n は各次元の格子点数である. 2.6.20 spmvtest4 ディレクトリ [output dir]/test において

> spmvtest4 matrix_filename_list iter [block]

と入力すると, matrix filename list の示す行列データファイルリストから行列データを読み込み, 各行列 とベクトル (1, . . . , 1)T との積を実行可能な行列格納形式について iter で指定された回数実行し, FLOPS 値を算出する. 入力可能な行列データ形式は Matrix Market 形式である. 必要なら block により, BSR, BSC 形式のブロックサイズを指定する.

(19)

2.6.21 spmvtest5

> spmvtest5 matrix_filename matrix_type iter [block]

と入力すると, matrix filename の示す行列データファイルから行列データを読み込み, 行列とベクトル (1, . . . , 1)T _{との積を行列格納形式 matrix type について iter で指定された回数実行し, FLOPS 値を算出}

する. 入力可能な行列データ形式は Matrix Market 形式である. 必要なら block により, BSR, BSC 形式のブロックサイズを指定する.

2.7 制限事項

現バージョンには以下の制限がある. • 行列格納形式 – VBR形式はマルチプロセス環境では使用できない. – CSR形式以外の格納形式は SA-AMG 前処理では使用できない. – マルチプロセス環境において必要な配列を直接定義する場合は, • double-double 型 4 倍精度演算

– 線型方程式解法のうち, Jacobi, Gauss-Seidel, SOR, IDR(s) 法では使用できない.

– 固有値解法のうち, CG, CR 法では使用できない.

– Hybrid前処理での内部反復解法のうち, Jacobi, Gauss-Seidel, SOR 法では使用できない.

– I+S, SA-AMG前処理では使用できない. • long double 型 4 倍精度演算 – Fortranインタフェースでは使用できない. – SA-AMG前処理では使用できない. • 前処理 – Jacobi, SSOR以外の前処理が選択され, かつ行列 A が CSR 形式でない場合, 前処理作成時に CSR形式の行列 A が作成される. – 非対称線型方程式解法として BiCG 法が選択された場合, SA-AMG 前処理は使用できない. – SA-AMG前処理はマルチスレッド計算には対応していない. – SAINV前処理の前処理行列作成部分は逐次実行される. • DD-AVX に関する制約 – 倍精度演算では SIMD 拡張命令は使用できない (自動的に Scalar 命令で実行される). – 倍々精度演算では CSR, BSR 以外は使用できない. – BCRS形式ではブロックサイズ 4x1 以外は AVX,AVX2 で使用できない (自動的に Scalar 命令で実行される). – OpenMPのスケジューリング方式に関する追加事項は倍精度では使用できない.

(20)

3 基本操作

本節では, ライブラリの使用方法について述べる. プログラムでは, 以下の処理を行う必要がある. • 初期化処理 • 行列の作成 • ベクトルの作成 • ソルバ (解法の情報を格納する構造体) の作成 • 行列, ベクトルへの値の代入 • 解法の設定 • 求解 • 終了処理 また, プログラムの先頭には以下のコンパイラ指示文を記述しなければならない. • C #include "lis.h"

• Fortran #include "lisf.h"

lis.h, lisf.hは, 導入時に ($INSTALLDIR)/include 下に格納される.

3.1 初期化・終了処理

初期化, 終了処理は以下のように記述する. 初期化処理はプログラムの最初に, 終了処理は最後に実行しなければならない. C 1: #include "lis.h"

2: LIS_INT main(LIS_INT argc, char* argv[]) 3: { 4: lis_initialize(&argc, &argv); 5: ... 6: lis_finalize(); 7: } Fortran 1: #include "lisf.h" 2: call lis_initialize(ierr) 3: ... 4: call lis_finalize(ierr) 初期化処理初期化処理を行うには, 関数

(21)

• Fortran subroutine lis_initialize(LIS_INTEGER ierr)

を用いる. この関数は, MPI の初期化, コマンドライン引数の取得等の初期化処理を行う.

LIS INTの既定値 int は, プリプロセッサマクロ LONGLONG が定義された場合には long long int に, また LIS INTEGER の既定値 integer は, プリプロセッサマクロ LONGLONG が定義された場合には integer*8 に置き換えられる.

終了処理

終了処理を行うには, 関数

• C LIS_INT lis_finalize()

• Fortran subroutine lis_finalize(LIS_INTEGER ierr) を用いる.

3.2 ベクトルの操作

ベクトル v の次数を global n とする. ベクトル v を nprocs 個のプロセスで行ブロック分割する場合の各 部分ベクトルの行数を local n とする. global n が nprocs で割り切れる場合は local n = global n / nprocs となる. 例えば, ベクトル v を (3.1) 式のように 2 プロセスで行ブロック分割する場合, global n と local n はそれぞれ 4 と 2 となる. v =       0 1 2 3       PE0 PE1 (3.1) (3.1)式のベクトル v を作成する場合, 逐次, マルチスレッド環境ではベクトル v そのものを, マルチプロ セス環境では各プロセスにプロセス数で行ブロック分割した部分ベクトルを作成する. ベクトル v を作成するプログラムは以下のように記述する. ただし, マルチプロセス環境のプロセス数は 2とする. C (逐次・マルチスレッド環境) 1: LIS_INT i,n; 2: LIS_VECTOR v; 3: n = 4; 4: lis_vector_create(0,&v); 5: lis_vector_set_size(v,0,n); /* or lis_vector_set_size(v,n,0); */ 6: 7: for(i=0;i<n;i++) 8: { 9: lis_vector_set_value(LIS_INS_VALUE,i,(double)i,v); 10: }

(22)

C (マルチプロセス環境)

1: LIS_INT i,n,is,ie; /* or LIS_INT i,ln,is,ie; */ 2: LIS_VECTOR v; 3: n = 4; /* ln = 2; */ 4: lis_vector_create(MPI_COMM_WORLD,&v); 5: lis_vector_set_size(v,0,n); /* lis_vector_set_size(v,ln,0); */ 6: lis_vector_get_range(v,&is,&ie); 7: for(i=is;i<ie;i++) 8: { 9: lis_vector_set_value(LIS_INS_VALUE,i,(double)i,v); 10: } Fortran (逐次・マルチスレッド環境) 1: LIS_INTEGER i,n 2: LIS_VECTOR v 3: n = 4 4: call lis_vector_create(0,v,ierr) 5: call lis_vector_set_size(v,0,n,ierr) 6: 7: do i=1,n 9: call lis_vector_set_value(LIS_INS_VALUE,i,DBLE(i),v,ierr) 10: enddo Fortran (マルチプロセス環境) 1: LIS_INTEGER i,n,is,ie 2: LIS_VECTOR v 3: n = 4 4: call lis_vector_create(MPI_COMM_WORLD,v,ierr) 5: call lis_vector_set_size(v,0,n,ierr) 6: call lis_vector_get_range(v,is,ie,ierr) 7: do i=is,ie-1 8: call lis_vector_set_value(LIS_INS_VALUE,i,DBLE(i),v,ierr); 9: enddo 変数宣言第 2 行のように LIS_VECTOR v; と宣言する. ベクトルの作成 ベクトル v の作成には, 関数

• C LIS_INT lis_vector_create(LIS_Comm comm, LIS_VECTOR *v)

• Fortran subroutine lis_vector_create(LIS_Comm comm, LIS_VECTOR v, LIS_INTEGER ierr) を用いる. comm には MPI コミュニケータを指定する. 逐次, マルチスレッド環境では comm の値は無視される.

次数の設定

(23)

• C LIS_INTEGER lis_vector_set_size(LIS_VECTOR v, LIS_INT local_n, LIS_INT global_n)

• Fortran subroutine lis_vector_set_size(LIS_VECTIR v, LIS_INTEGER local_n, LIS_INTEGER global_n, LIS_INTEGER ierr)

を用いる. local n か global n のどちらか一方を与えなければならない.

逐次, マルチスレッド環境では, local n は global n に等しい. したがって, lis_vector_set_size(v,n,0) と lis_vector_set_size(v,0,n) は, いずれも次数 n のベクトルを作成する. マルチプロセス環境においては, lis_vector_set_size(v,n,0) は各プロセス上に次数 n の部分ベクト ルを作成する. 一方, lis_vector_set_size(v,0,n) は各プロセス p 上に次数 mpの部分ベクトルを作成す る. mpはライブラリ側で決定される. 値の代入 ベクトル v の第 i 行に値を代入するには, 関数

• C LIS_INT lis_vector_set_value(LIS_INT flag, LIS_INT i, LIS_SCALAR value, LIS_VECTOR v)

• Fortran subroutine lis_vector_set_value(LIS_INTEGER flag, LIS_INTEGER i, LIS_SCALAR value, LIS_VECTOR v, LIS_INTEGER ierr)

を用いる. マルチプロセス環境では, 部分ベクトルの第 i 行ではなく，全体ベクトルの第 i 行を指定する. flagには

LIS INS VALUE 挿入: v[i] = value

LIS ADD VALUE 加算代入: v[i] = v[i] + value のどちらかを指定する.

ベクトルの複製

既存のベクトルと同じ情報を持つベクトルを作成するには, 関数

• C LIS_INT lis_vector_duplicate(LIS_VECTOR vin, LIS_VECTOR *vout) • Fortran subroutine lis_vector_duplicate(LIS_VECTOR vin, LIS_VECTOR vout,

LIS_INTEGER ierr)

を用いる. 第 1 引数 LIS_VECTOR vin は LIS_MATRIX を指定することも可能である. この関数はベクトルの要素の値は複製しない. 値も複製する場合は，この関数の後に

• C LIS_INT lis_vector_copy(LIS_VECTOR vsrc, LIS_VECTOR vdst)

• Fortran subroutine lis_vector_copy(LIS_VECTOR vsrc, LIS_VECTOR vdst, LIS_INTEGER ierr) を呼び出す.

ベクトルの破棄

(24)

• C LIS_INT lis_vector_destroy(LIS_VECTOR v)

• Fortran subroutine lis_vector_destroy(LIS_VECTOR v, LIS_INTEGER ierr) を用いる.

3.3 行列の操作

係数行列 A の次数を global n × global n とする. 行列 A を nprocs 個のプロセスで行ブロック分割する 場合の各ブロックの行数を local n とする. global n が nprocs で割り切れる場合は local n = global n / nprocsとなる. 例えば, 行列 A を (3.2) 式のように 2 個のプロセスで行ブロック分割する場合, global n と local nはそれぞれ 4 と 2 となる. A =       2 1 1 2 1 1 2 1 1 2       PE0 PE1 (3.2) 目的の格納形式の行列を作成するには以下の 3 つの方法がある. 方法 1: ライブラリ関数を用いて目的の格納形式の配列を定義する場合 (3.2)式の行列 A を CSR 形式で作成する場合, 逐次, マルチスレッド環境では行列 A そのものを, マルチプ ロセス環境では各プロセスにプロセス数で行ブロック分割した部分行列を作成する. 行列 A を CSR 形式で作成するプログラムは以下のように記述する. ただし, マルチプロセス環境のプロ セス数は 2 とする. C (逐次・マルチスレッド環境) 1: LIS_INT i,n; 2: LIS_MATRIX A; 3: n = 4; 4: lis_matrix_create(0,&A); 5: lis_matrix_set_size(A,0,n); /* or lis_matrix_set_size(A,n,0); */ 6: for(i=0;i<n;i++) {

7: if( i>0 ) lis_matrix_set_value(LIS_INS_VALUE,i,i-1,1.0,A); 8: if( i<n-1 ) lis_matrix_set_value(LIS_INS_VALUE,i,i+1,1.0,A); 9: lis_matrix_set_value(LIS_INS_VALUE,i,i,2.0,A);

10: }

11: lis_matrix_set_type(A,LIS_MATRIX_CSR); 12: lis_matrix_assemble(A);

(25)

C (マルチプロセス環境) 1: LIS_INT i,n,gn,is,ie; 2: LIS_MATRIX A; 3: gn = 4; /* or n=2 */ 4: lis_matrix_create(MPI_COMM_WORLD,&A); 5: lis_matrix_set_size(A,0,gn); /* lis_matrix_set_size(A,n,0); */ 6: lis_matrix_get_size(A,&n,&gn); 7: lis_matrix_get_range(A,&is,&ie); 8: for(i=is;i<ie;i++) {

9: if( i>0 ) lis_matrix_set_value(LIS_INS_VALUE,i,i-1,1.0,A); 10: if( i<gn-1 ) lis_matrix_set_value(LIS_INS_VALUE,i,i+1,1.0,A); 11: lis_matrix_set_value(LIS_INS_VALUE,i,i,2.0,A); 12: } 13: lis_matrix_set_type(A,LIS_MATRIX_CSR); 14: lis_matrix_assemble(A); Fortran (逐次・マルチスレッド環境) 1: LIS_INTEGER i,n 2: LIS_MATRIX A 3: n = 4 4: call lis_matrix_create(0,A,ierr) 5: call lis_matrix_set_size(A,0,n,ierr) 6: do i=1,n

7: if( i>1 ) call lis_matrix_set_value(LIS_INS_VALUE,i,i-1,1.0d0,A,ierr) 8: if( i<n ) call lis_matrix_set_value(LIS_INS_VALUE,i,i+1,1.0d0,A,ierr) 9: call lis_matrix_set_value(LIS_INS_VALUE,i,i,2.0d0,A,ierr) 10: enddo 11: call lis_matrix_set_type(A,LIS_MATRIX_CSR,ierr) 12: call lis_matrix_assemble(A,ierr) Fortran (マルチプロセス環境) 1: LIS_INTEGER i,n,gn,is,ie 2: LIS_MATRIX A 3: gn = 4 4: call lis_matrix_create(MPI_COMM_WORLD,A,ierr) 5: call lis_matrix_set_size(A,0,gn,ierr) 6: call lis_matrix_get_size(A,n,gn,ierr) 7: call lis_matrix_get_range(A,is,ie,ierr) 8: do i=is,ie-1

9: if( i>1 ) call lis_matrix_set_value(LIS_INS_VALUE,i,i-1,1.0d0,A,ierr) 10: if( i<gn ) call lis_matrix_set_value(LIS_INS_VALUE,i,i+1,1.0d0,A,ierr) 11: call lis_matrix_set_value(LIS_INS_VALUE,i,i,2.0d0,A,ierr) 12: enddo 13: call lis_matrix_set_type(A,LIS_MATRIX_CSR,ierr) 14: call lis_matrix_assemble(A,ierr) 変数宣言第 2 行のように LIS_MATRIX A; と宣言する. 行列の作成

(26)

行列 A の作成には, 関数

• C LIS_INT lis_matrix_create(LIS_Comm comm, LIS_MATRIX *A)

• Fortran subroutine lis_matrix_create(LIS_Comm comm, LIS_MATRIX A, LIS_INTEGER ierr) を用いる. comm には MPI コミュニケータを指定する. 逐次, マルチスレッド環境では, comm の値は無視される.

次数の設定

次数の設定には, 関数

• C LIS_INT lis_matrix_set_size(LIS_MATRIX A, LIS_INT local_n, LIS_INT global_n) • Fortran subroutine lis_matrix_set_size(LIS_MATRIX A, LIS_INTEGER local_n,

LIS_INTEGER global_n, LIS_INTEGER ierr)

を用いる. local n か global n のどちらか一方を与えなければならない.

逐次, マルチスレッド環境では, local n は global n に等しい. したがって, lis_matrix_set_size(A,n,0) と lis_matrix_set_size(A,0,n) は, いずれも次数 n× n の行列を作成する. マルチプロセス環境においては, lis_matrix_set_size(A,n,0) は各プロセス上に次数 n× N の部分行 列を作成する. N は n の総和である. 一方, lis_matrix_set_size(A,0,n) は各プロセス p 上に次数 mp× n の部分行列を作成する. mpはライブラリ側で決定される. 値の代入 行列 A の第 i 行第 j 列に値を代入するには, 関数

• C LIS_INT lis_matrix_set_value(LIS_INT flag, LIS_INT i, LIS_INT j, LIS_SCALAR value, LIS_MATRIX A)

• Fortran subroutine lis_matrix_set_value(LIS_INTEGER flag, LIS_INTEGER i, LIS_INTEGER j, LIS_SCALAR value, LIS_MATRIX A, LIS_INTEGER ierr)

を用いる. マルチプロセス環境では, 全体行列の第 i 行第 j 列を指定する. flag には LIS INS VALUE 挿入: A[i, j] = value

LIS ADD VALUE 加算代入: A[i, j] = A[i, j] + value のどちらかを指定する.

行列格納形式の設定

行列の格納形式を設定するには, 関数

• C LIS_INT lis_matrix_set_type(LIS_MATRIX A, LIS_INT matrix_type)

• Fortran subroutine lis_matrix_set_type(LIS_MATRIX A, LIS_INTEGER matrix_type, LIS_INTEGER ierr)

(27)

を用いる. 行列作成時の A の matrix_type は LIS_MATRIX_CSR である. 以下に対応する格納形式を示す.

格納形式 matrix_type

Compressed Sparse Row (CSR) {LIS_MATRIX_CSR|1} Compressed Sparse Column (CSC) {LIS_MATRIX_CSC|2} Modified Compressed Sparse Row (MSR) {LIS_MATRIX_MSR|3}

Diagonal (DIA) {LIS_MATRIX_DIA|4}

Ellpack-Itpack Generalized Diagonal (ELL) {LIS_MATRIX_ELL|5}

Jagged Diagonal (JAD) {LIS_MATRIX_JAD|6}

Block Sparse Row (BSR) {LIS_MATRIX_BSR|7}

Block Sparse Column (BSC) {LIS_MATRIX_BSC|8} Variable Block Row (VBR) {LIS_MATRIX_VBR|9}

Coordinate (COO) {LIS_MATRIX_COO|10}

Dense (DNS) {LIS_MATRIX_DNS|11}

行列の組み立て

行列の要素と格納形式を設定した後, 関数

• C LIS_INT lis_matrix_assemble(LIS_MATRIX A)

• Fortran subroutine lis_matrix_assemble(LIS_MATRIX A, LIS_INTEGER ierr)

を呼び出す. lis_matrix_assemble は lis_matrix_set_type で指定された格納形式に組み立てられる.

行列の破棄

不要になった行列をメモリから破棄するには,

• C LIS_INT lis_matrix_destroy(LIS_MATRIX A)

• Fortran subroutine lis_matrix_destroy(LIS_MATRIX A, LIS_INTEGER ierr) を用いる. 方法 2: 目的の格納形式の配列を直接定義する場合 (3.2)式の行列 A を CSR 形式で作成する場合, 逐次, マルチスレッド環境では行列 A そのものを, マルチプ ロセス環境では各プロセスにプロセス数で行ブロック分割した部分行列を作成する. 行列 A を CSR 形式で作成するプログラムは以下のように記述する. ただし, マルチプロセス環境のプロ セス数は 2 とする.

(28)

C (逐次・マルチスレッド環境) 1: LIS_INT i,k,n,nnz; 2: LIS_INT *ptr,*index; 3: LIS_SCALAR *value; 4: LIS_MATRIX A; 5: n = 4; nnz = 10; k = 0; 6: lis_matrix_malloc_csr(n,nnz,&ptr,&index,&value); 7: lis_matrix_create(0,&A); 8: lis_matrix_set_size(A,0,n); /* or lis_matrix_set_size(A,n,0); */ 9: 10: for(i=0;i<n;i++) 11: {

12: if( i>0 ) {index[k] = i-1; value[k] = 1; k++;} 13: index[k] = i; value[k] = 2; k++;

14: if( i<n-1 ) {index[k] = i+1; value[k] = 1; k++;} 15: ptr[i+1] = k; 16: } 17: ptr[0] = 0; 18: lis_matrix_set_csr(nnz,ptr,index,value,A); 19: lis_matrix_assemble(A); C (マルチプロセス環境) 1: LIS_INT i,k,n,nnz,is,ie; 2: LIS_INT *ptr,*index; 3: LIS_SCALAR *value; 4: LIS_MATRIX A; 5: n = 2; nnz = 5; k = 0; 6: lis_matrix_malloc_csr(n,nnz,&ptr,&index,&value); 7: lis_matrix_create(MPI_COMM_WORLD,&A); 8: lis_matrix_set_size(A,n,0); 9: lis_matrix_get_range(A,&is,&ie); 10: for(i=is;i<ie;i++) 11: {

12: if( i>0 ) {index[k] = i-1; value[k] = 1; k++;} 13: index[k] = i; value[k] = 2; k++;

14: if( i<n-1 ) {index[k] = i+1; value[k] = 1; k++;} 15: ptr[i-is+1] = k; 16: } 17: ptr[0] = 0; 18: lis_matrix_set_csr(nnz,ptr,index,value,A); 19: lis_matrix_assemble(A); 配列の関連付け CSR形式の配列をライブラリが扱えるよう行列 A に関連付けるには, 関数

• C LIS_INT lis_matrix_set_csr(LIS_INT nnz, LIS_INT row[], LIS_INT index[], LIS_SCALAR value[], LIS_MATRIX A)

• Fortran subroutine lis_matrix_set_csr(LIS_INTEGER nnz, LIS_INTEGER row(),

LIS_INTEGER index(), LIS_SCALAR value(), LIS_MATRIX A, LIS_INTEGER ierr)

(29)

方法 3: 外部ファイルから行列, ベクトルデータを読み込む場合 外部ファイルから (3.2) 式の行列 A を CSR 形式で読み込む場合, プログラムは以下のように記述する. C (逐次・マルチスレッド・マルチプロセス環境) 1: LIS_MATRIX A; 3: lis_matrix_create(LIS_COMM_WORLD,&A); 6: lis_matrix_set_type(A,LIS_MATRIX_CSR); 7: lis_input_matrix(A,"matvec.mtx"); Fortran (逐次・マルチスレッド・マルチプロセス環境) 1: LIS_MATRIX A 3: call lis_matrix_create(LIS_COMM_WORLD,A,ierr) 6: call lis_matrix_set_type(A,LIS_MATRIX_CSR,ierr) 7: call lis_input_matrix(A,’matvec.mtx’,ierr)

Matrix Market形式による外部ファイル matvec.mtx の記述例を以下に示す.

%%MatrixMarket matrix coordinate real general 4 4 10 1 0 1 2 1.0e+00 1 1 2.0e+00 2 3 1.0e+00 2 1 1.0e+00 2 2 2.0e+00 3 4 1.0e+00 3 2 1.0e+00 3 3 2.0e+00 4 4 2.0e+00 4 3 1.0e+00 外部ファイルから (3.2) 式の行列 A を CSR 形式で, また (3.1) 式のベクトル b を読み込む場合のプログラ ムは以下のように記述する. C (逐次・マルチスレッド・マルチプロセス環境) 1: LIS_MATRIX A; 2: LIS_VECTOR b,x; 3: lis_matrix_create(LIS_COMM_WORLD,&A); 4: lis_vector_create(LIS_COMM_WORLD,&b); 5: lis_vector_create(LIS_COMM_WORLD,&x); 6: lis_matrix_set_type(A,LIS_MATRIX_CSR); 7: lis_input(A,b,x,"matvec.mtx"); Fortran (逐次・マルチスレッド・マルチプロセス環境) 1: LIS_MATRIX A 2: LIS_VECTOR b,x 3: call lis_matrix_create(LIS_COMM_WORLD,A,ierr) 4: call lis_vector_create(LIS_COMM_WORLD,b,ierr) 5: call lis_vector_create(LIS_COMM_WORLD,x,ierr) 6: call lis_matrix_set_type(A,LIS_MATRIX_CSR,ierr) 7: call lis_input(A,b,x,’matvec.mtx’,ierr)

(30)

%%MatrixMarket matrix coordinate real general 4 4 10 1 0 1 2 1.0e+00 1 1 2.0e+00 2 3 1.0e+00 2 1 1.0e+00 2 2 2.0e+00 3 4 1.0e+00 3 2 1.0e+00 3 3 2.0e+00 4 4 2.0e+00 4 3 1.0e+00 1 0.0e+00 2 1.0e+00 3 2.0e+00 4 3.0e+00 外部ファイルからの読み込み 外部ファイルから行列 A のデータを読み込むには, 関数

• C LIS_INT lis_input_matrix(LIS_MATRIX A, char *filename) • Fortran subroutine lis_input_matrix(LIS_MATRIX A,

character filename, LIS_INTEGER ierr)

を用いる. filename にはファイルパスを指定する. 対応するファイル形式は以下の通りである (ファイル形式については付録 A を参照).

• Matrix Market 形式 • Harwell-Boeing 形式

外部ファイルから行列 A とベクトル b, x のデータを読み込むには, 関数

• C LIS_INT lis_input(LIS_MATRIX A, LIS_VECTOR b, LIS_VECTOR x, char *filename) • Fortran subroutine lis_input(LIS_MATRIX A, LIS_VECTOR b, LIS_VECTOR x,

character filename, LIS_INTEGER ierr)

を用いる. filename にはファイルパスを指定する. 対応するファイル形式は以下の通りである (ファイル形式については付録 A を参照). • 拡張 Matrix Market 形式 • Harwell-Boeing 形式

3.4 線型方程式の求解

線型方程式 Ax = b を指定された解法で解く場合, プログラムは以下のように記述する.

(31)

C (逐次・マルチスレッド・マルチプロセス環境) 1: LIS_MATRIX A; 2: LIS_VECTOR b,x; 3: LIS_SOLVER solver; 4: 5: /* 行列とベクトルの作成 */ 6: 7: lis_solver_create(&solver);

8: lis_solver_set_option("-i bicg -p none",solver); 9: lis_solver_set_option("-tol 1.0e-12",solver); 10: lis_solve(A,b,x,solver); Fortran (逐次・マルチスレッド・マルチプロセス環境) 1: LIS_MATRIX A 2: LIS_VECTOR b,x 3: LIS_SOLVER solver 4: 5: /* 行列とベクトルの作成 */ 6: 7: call lis_solver_create(solver,ierr)

8: call lis_solver_set_option(’-i bicg -p none’,solver,ierr) 9: call lis_solver_set_option(’-tol 1.0e-12’,solver,ierr) 10: call lis_solve(A,b,x,solver,ierr)

ソルバの作成

ソルバ (線型方程式解法の情報を格納する構造体) を作成するには, 関数 • C LIS_INT lis_solver_create(LIS_SOLVER *solver)

• Fortran subroutine lis_solver_create(LIS_SOLVER solver, LIS_INTEGER ierr) を用いる.

オプションの設定

線型方程式解法をソルバに設定するには, 関数

• C LIS_INT lis_solver_set_option(char *text, LIS_SOLVER solver)

• Fortran subroutine lis_solver_set_option(character text, LIS_SOLVER solver, LIS_INTEGER ierr)

または

• C LIS_INT lis_solver_set_optionC(LIS_SOLVER solver)

• Fortran subroutine lis_solver_set_optionC(LIS_SOLVER solver, LIS_INTEGER ierr) を用いる. lis_solver_set_optionC は, ユーザプログラム実行時にコマンドラインで指定されたオプションをソルバに設定する関数である.

以下に指定可能なコマンドラインオプションを示す. -i {cg|1}は-i cg または-i 1 を意味する. -maxiter [1000]は, -maxiter の既定値が 1000 であることを意味する.

(32)

線型方程式解法に関するオプション (既定値: -i bicg) 線型方程式解法オプション補助オプション CG -i {cg|1} BiCG -i {bicg|2} CGS -i {cgs|3} BiCGSTAB -i {bicgstab|4}

BiCGSTAB(l) -i {bicgstabl|5} -ell [2] 次数 l

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) BiCGSafe -i {bicgsafe|13} CR -i {cr|14} BiCR -i {bicr|15} CRS -i {crs|16} BiCRSTAB -i {bicrstab|17} GPBiCR -i {gpbicr|18} BiCRSafe -i {bicrsafe|19}

FGMRES(m) -i {fgmres|20} -restart [40] リスタート値 m

IDR(s) -i {idrs|21} -irestart [2] リスタート値 s

(33)

前処理に関するオプション (既定値: -p none)

前処理オプション補助オプション

なし -p {none|0}

Jacobi -p {jacobi|1}

ILU(k) -p {ilu|2} -ilu_fill [0] フィルインレベル k

SSOR -p {ssor|3} -ssor_w [1.0] 緩和係数 ω (0 < ω < 2)

Hybrid -p {hybrid|4} -hybrid_i [sor] 線型方程式解法

-hybrid_maxiter [25] 最大反復回数

-hybrid_tol [1.0e-3] 収束判定基準

-hybrid_w [1.5] SORの緩和係数 ω (0 < ω < 2)

-hybrid_ell [2] BiCGSTAB(l)の次数 l

-hybrid_restart [40] GMRES(m), Orthomin(m)の リスタート値 m

I+S -p {is|5} -is_alpha [1.0] I + αS(m)のパラメータ α

-is_m [3] I + αS(m)_{のパラメータ m}

SAINV -p {sainv|6} -sainv_drop [0.05] ドロップ基準

SA-AMG -p {saamg|7} -saamg_unsym [false] 非対称版の選択

(行列構造は対称とする) -saamg_theta [0.05|0.12] ドロップ基準 a2ij≤ θ

2_|a

ii||ajj|

(対称|非対称)

Crout ILU -p {iluc|8} -iluc_drop [0.05] ドロップ基準

-iluc_rate [5.0] 最大フィルイン数の倍率

ILUT -p {ilut|9} -ilut_drop [0.05] ドロップ基準

-ilut_rate [5.0] 最大フィルイン数の倍率

(34)

その他のオプションオプション

-maxiter [1000] 最大反復回数

-tol [1.0e-12] 収束判定基準 tol

-tol_w [1.0] 収束判定基準 tolw -print [0] 残差履歴の出力 -print {none|0} 残差履歴を出力しない -print {mem|1} 残差履歴をメモリに保存する -print {out|2} 残差履歴を標準出力に書き出す -print {all|3} 残差履歴をメモリに保存し, 標準出力に書き出す -scale [0] スケーリングの選択. 結果は元の行列, ベクトルに上書きされる -scale {none|0} スケーリングなし

-scale {jacobi|1} Jacobiスケーリング D−1Ax = D−1b (Dは A = (aij)の対角部分) -scale {symm_diag|2}対角スケーリング D−1/2AD−1/2x = D−1/2b (D−1/2は対角要素の値が 1/√aiiである対角行列) -initx_zeros [1] 初期ベクトル x0 -initx_zeros {false|0} 与えられた値を使用 -initx_zeros {true|1} すべての要素の値を 0 にする -conv_cond [0] 収束条件

-conv_cond {nrm2_r|0} ||b − Ax||2≤ tol ∗ ||b||2

-conv_cond {nrm2_b|1} ||b − Ax||2≤ tol ∗ ||b − Ax0||2

-conv_cond {nrm1_b|2} ||b − Ax||1≤ tolw∗ ||b||1+ tol

-omp_num_threads [t] 実行スレッド数 (tは最大スレッド数) -storage [0] 行列格納形式 -storage_block [2] BSR, BSC形式のブロックサイズ -f [0] 線型方程式解法の精度 -f {double|0} 倍精度 -f {quad|1} 4倍精度求解 線型方程式 Ax = b を解くには, 関数

• C LIS_INT lis_solve(LIS_MATRIX A, LIS_VECTOR b, LIS_VECTOR x, LIS_SOLVER solver)

• Fortran subroutine lis_solve(LIS_MATRIX A, LIS_VECTOR b, LIS_VECTOR x, LIS_SOLVER solver, LIS_INTEGER ierr)

(35)

3.5 固有値問題の求解

固有値問題 Ax = λx を指定の解法で解く場合, プログラムは以下のように記述する. C (逐次・マルチスレッド・マルチプロセス環境) 1: LIS_MATRIX A; 2: LIS_VECTOR x; 3: LIS_REAL evalue; 4: LIS_ESOLVER esolver; 5: 6: /* 行列とベクトルの作成 */ 7: 8: lis_esolver_create(&esolver);

9: lis_esolver_set_option("-e ii -i bicg -p none",esolver); 10: lis_esolver_set_option("-etol 1.0e-12 -tol 1.0e-12",esolver); 11: lis_esolve(A,x,evalue,esolver); Fortran (逐次・マルチスレッド・マルチプロセス環境) 1: LIS_MATRIX A 2: LIS_VECTOR x 3: LIS_REAL evalue 4: LIS_ESOLVER esolver 5: 6: /* 行列とベクトルの作成 */ 7: 8: call lis_esolver_create(esolver,ierr)

9: call lis_esolver_set_option(’-e ii -i bicg -p none’,esolver,ierr) 10: call lis_esolver_set_option(’-etol 1.0e-12 -tol 1.0e-12’,esolver,ierr) 11: call lis_esolve(A,x,evalue,esolver,ierr)

ソルバの作成

ソルバ (固有値解法の情報を格納する構造体) を作成するには, 関数 • C LIS_INT lis_esolver_create(LIS_ESOLVER *esolver)

• Fortran subroutine lis_esolver_create(LIS_ESOLVER esolver, LIS_INTEGER ierr) を用いる.

オプションの設定

固有値解法をソルバに設定するには, 関数

• C LIS_INT lis_esolver_set_option(char *text, LIS_ESOLVER esolver)

• Fortran subroutine lis_esolver_set_option(character text, LIS_ESOLVER esolver, LIS_INTEGER ierr)

または

• C LIS_INT lis_esolver_set_optionC(LIS_ESOLVER esolver)

(36)

を用いる. lis_esolver_set_optionC は, ユーザプログラム実行時にコマンドラインで指定されたオプションをソルバに設定する関数である. 以下に指定可能なコマンドラインオプションを示す. -e {pi|1}は-e pi または-e 1 を意味する. -emaxiter [1000]は, -emaxiter の既定値が 1000 であることを意味する. 固有値解法に関するオプション (既定値: -e pi) 固有値解法オプション補助オプション Power -e {pi|1} Inverse -e {ii|2} -i [cg] 線型方程式解法

Approximate Inverse -e {aii|3} -i [cg] 線型方程式解法

Rayleigh Quotient -e {rqi|4} -i [cg] 線型方程式解法

Subspace -e {si|5} -ss [2] 部分空間の大きさ -m [0] モード番号 Lanczos -e {li|6} -ss [2] 部分空間の大きさ -m [0] モード番号 CG -e {cg|7} CR -e {cr|8} 前処理に関するオプション (既定値: -p ilu) 前処理オプション補助オプションなし -p {none|0} Jacobi -p {jacobi|1}

ILU(k) -p {ilu|2} -ilu_fill [0] フィルインレベル k

SSOR -p {ssor|3} -ssor_w [1.0] 緩和係数 ω (0 < ω < 2)

Hybrid -p {hybrid|4} -hybrid_i [sor] 線型方程式解法

-hybrid_maxiter [25] 最大反復回数

-hybrid_tol [1.0e-3] 収束判定基準

-hybrid_w [1.5] SORの緩和係数 ω (0 < ω < 2)

-hybrid_ell [2] BiCGSTAB(l)の次数 l

-hybrid_restart [40] GMRES(m), Orthomin(m)の リスタート値 m

I+S -p {is|5} -is_alpha [1.0] I + αS(m)_{のパラメータ α}

-is_m [3] I + αS(m)のパラメータ m

SAINV -p {sainv|6} -sainv_drop [0.05] ドロップ基準

SA-AMG -p {saamg|7} -saamg_unsym [false] 非対称版の選択

(行列構造は対称とする) -saamg_theta [0.05|0.12] ドロップ基準 a2

ij≤ θ2|aii||ajj|

(対称|非対称)

crout ILU -p {iluc|8} -iluc_drop [0.05] ドロップ基準

-iluc_rate [5.0] 最大フィルイン数の倍率

ILUT -p {ilut|9} -ilut_drop [0.05] ドロップ基準

-ilut_rate [5.0] 最大フィルイン数の倍率

(37)

その他のオプションオプション -emaxiter [1000] 最大反復回数 -etol [1.0e-12] 収束判定基準 -eprint [0] 残差履歴の出力 -eprint {none|0} 残差履歴を出力しない -eprint {mem|1} 残差履歴をメモリに保存する -eprint {out|2} 残差履歴を標準出力に書き出す -eprint {all|3} 残差履歴をメモリに保存し, 標準出力に書き出す

-ie [ii] Subspace, Lanczos法の内部で使用する固有値解法の指定

-ie {pi|1} Power (Subspaceのみ)

-ie {ii|2} Inverse

-ie {aii|3} Approximate Inverse -ie {rqi|4} Rayleigh Quotient

-shift [0.0] 固有値のシフト量 -initx_ones [1] 初期ベクトル x0 -initx_ones {false|0} 与えられた値を使用 -initx_ones {true|1} すべての要素の値を 1 にする -omp_num_threads [t] 実行スレッド数 tは最大スレッド数 -estorage [0] 行列格納形式 -estorage_block [2] BSR, BSC形式のブロックサイズ -ef [0] 固有値解法の精度 -ef {double|0} 倍精度 -ef {quad|1} 4倍精度求解 固有値問題 Ax = λx を解くには, 関数

• C LIS_INT lis_esolve(LIS_MATRIX A, LIS_VECTOR x, LIS_REAL evalue, LIS_ESOLVER esolver)

• Fortran subroutine lis_esolve(LIS_MATRIX A, LIS_VECTOR x, LIS_ESOLVER esolver, LIS_INTEGER ierr)

を用いる.

3.6 プログラムの作成

(38)

行列 A は次数 12 の 3 重対角行列 A =          2 −1 −1 2 −1 . ._. . ._. . ._. −1 2 −1 −1 2          である. 右辺ベクトル b は解 x がすべて 1 となるよう設定される. このプログラムはディレクトリ lis-($VERSION)/test にある.

(39)

検証プログラム: test4.c

1: #include <stdio.h> 2: #include "lis.h"

3: main(LIS_INT argc, char *argv[]) 4: { 5: LIS_INT i,n,gn,is,ie,iter; 6: LIS_MATRIX A; 7: LIS_VECTOR b,x,u; 8: LIS_SOLVER solver; 9: n = 12; 10: lis_initialize(&argc,&argv); 11: lis_matrix_create(LIS_COMM_WORLD,&A); 12: lis_matrix_set_size(A,0,n); 13: lis_matrix_get_size(A,&n,&gn) 14: lis_matrix_get_range(A,&is,&ie) 15: for(i=is;i<ie;i++) 16: {

17: if( i>0 ) lis_matrix_set_value(LIS_INS_VALUE,i,i-1,-1.0,A); 18: if( i<gn-1 ) lis_matrix_set_value(LIS_INS_VALUE,i,i+1,-1.0,A); 19: lis_matrix_set_value(LIS_INS_VALUE,i,i,2.0,A); 20: } 21: lis_matrix_set_type(A,LIS_MATRIX_CSR); 22: lis_matrix_assemble(A); 23: 24: lis_vector_duplicate(A,&u); 25: lis_vector_duplicate(A,&b); 26: lis_vector_duplicate(A,&x); 27: lis_vector_set_all(1.0,u); 28: lis_matvec(A,u,b); 29: 30: lis_solver_create(&solver); 31: lis_solver_set_optionC(solver); 32: lis_solve(A,b,x,solver); 33: lis_solver_get_iter(solver,&iter);

34: printf("number of iterations = %d\n",iter); 35: lis_vector_print(x); 36: lis_matrix_destroy(A); 37: lis_vector_destroy(u); 38: lis_vector_destroy(b); 39: lis_vector_destroy(x); 40: lis_solver_destroy(solver); 41: lis_finalize(); 42: return 0; 43: } }

(40)

検証プログラム: test4f.F 1: implicit none 2: 3:#include "lisf.h" 4: 5: LIS_INTEGER i,n,gn,is,ie,iter,ierr 6: LIS_MATRIX A 7: LIS_VECTOR b,x,u 8: LIS_SOLVER solver 9: n = 12 10: call lis_initialize(ierr) 11: call lis_matrix_create(LIS_COMM_WORLD,A,ierr) 12: call lis_matrix_set_size(A,0,n,ierr) 13: call lis_matrix_get_size(A,n,gn,ierr) 14: call lis_matrix_get_range(A,is,ie,ierr) 15: do i=is,ie-1

16: if( i>1 ) call lis_matrix_set_value(LIS_INS_VALUE,i,i-1,-1.0d0,

17: . A,ierr)

18: if( i<gn ) call lis_matrix_set_value(LIS_INS_VALUE,i,i+1,-1.0d0,

19: . A,ierr) 20: call lis_matrix_set_value(LIS_INS_VALUE,i,i,2.0d0,A,ierr) 21: enddo 22: call lis_matrix_set_type(A,LIS_MATRIX_CSR,ierr) 23: call lis_matrix_assemble(A,ierr) 24: 25: call lis_vector_duplicate(A,u,ierr) 26: call lis_vector_duplicate(A,b,ierr) 27: call lis_vector_duplicate(A,x,ierr) 28: call lis_vector_set_all(1.0d0,u,ierr) 29: call lis_matvec(A,u,b,ierr) 30: 31: call lis_solver_create(solver,ierr) 32: call lis_solver_set_optionC(solver,ierr) 33: call lis_solve(A,b,x,solver,ierr) 34: call lis_solver_get_iter(solver,iter,ierr) 35: write(*,*) ’number of iterations = ’,iter 36: call lis_vector_print(x,ierr) 37: call lis_matrix_destroy(A,ierr) 38: call lis_vector_destroy(b,ierr) 39: call lis_vector_destroy(x,ierr) 40: call lis_vector_destroy(u,ierr) 41: call lis_solver_destroy(solver,ierr) 42: call lis_finalize(ierr) 43: 44: stop 45: end

3.7 実行ファイルの生成

test4.cから実行ファイルを生成する方法について述べる. ディレクトリ lis-($VERSION)/test にあ

る検証プログラム test4.c を SGI Altix 3700 上の Intel C/C++ Compiler (icc), Intel Fortran Compiler

(ifort)でコンパイルする場合の例を以下に示す. SA-AMG 前処理には Fortran90 で記述されたコードが

含まれるため, SA-AMG 前処理を使用する場合には Fortran90 コンパイラでリンクしなければならない. また, マルチプロセス環境ではプリプロセッサマクロ USE_MPI が定義されなければならない. 64bit 整数型を

(41)

使用する場合, C プログラムについてはプリプロセッサマクロ_LONGLONG が, また Fortran プログラムについてはプリプロセッサマクロ LONGLONG が定義されなければならない.

逐次環境

コンパイル

> icc -c -I($INSTALLDIR)/include test4.c

リンク

> icc -o test4 test4.o -llis

リンク(--enable-saamg)

> ifort -nofor_main -o test4 test4.o -llis

マルチスレッド環境

コンパイル

> icc -c -openmp -I($INSTALLDIR)/include test4.c

リンク

> icc -openmp -o test4 test4.o -llis

> ifort -nofor_main -openmp -o test4 test4.o -llis

マルチプロセス環境

コンパイル

> icc -c -DUSE_MPI -I($INSTALLDIR)/include test4.c

リンク

> icc -o test4 test4.o -llis -lmpi

> ifort -nofor_main -o test4 test4.o -llis -lmpi

マルチスレッド・マルチプロセス環境

コンパイル

> icc -c -openmp -DUSE_MPI -I($INSTALLDIR)/include test4.c

リンク

> icc -openmp -o test4 test4.o -llis -lmpi

> ifort -nofor_main -openmp -o test4 test4.o -llis -lmpi

次に, test4f.F から実行ファイルを生成する方法について述べる. ディレクトリ lis-($VERSION)/test にある検証プログラム test4f.F を SGI Altix 3700 上の Intel Fortran Compiler (ifort) でコンパイルする場合の例を以下に示す. Fortran のユーザプログラムにはコンパイラ指示文が含まれるため, プリプロセッサを使用するようコンパイラに指示しなければならない. ifort の場合のオプションは-fpp である.

逐次環境

コンパイル

> ifort -c -fpp -I($INSTALLDIR)/include test4f.F

リンク

> ifort -o test4f test4f.o -llis

DD-AVX User s Manual (based on Lis ) Toshiaki Hishinuma a, Akihiro Fujii a, Teruo Tanaka a, Hidehiko Hasegawa b