第12回 お試しアカウント付き
並列プログラミング講習会（試行）「ライブラリ利用：高性能プログラミング初級入門」
東京大学情報基盤センター
スーパーコンピューティング部門

2

並列プログラミング講習会（試行）講習会概略

開催日： 2010年9月27日（月）10：30 - 17：00
2010年9月28日（火）10：30 - 17：00



場所：東京大学情報基盤センター 4階 413遠隔講義室



講習会プ グラ講習会プログラム：

9月27日（月） 10：00 - 10：30 受付 10：30 - 12：30 ノートパソコンの設定、テストプログラムの実行など（演習）トパ
設定
トプ グ ム 実行など（演習）（講師：片桐） 14：00 - 15：00 並列プログラミングの基本（座学）（講師：片桐） 15：15 - 17：00 プログラム実習Ⅰ（BLAS）（演習）（講師：片桐）

9月28日（火） 10：30 - 12：30 プログラミング実習Ⅱ（LAPACK,
プログラミング実習Ⅱ（LAPACK ScaLAPACK）（演習）（講師：片桐） 14：00 - 15：30 プログラミング実習Ⅲ（Lis）（演習）（講師：伊藤） 15：45 - 17：00 線形代数の基礎（座学）（講師：伊藤）3

並列プログラミング講習会（試行）テストプログラム起動

25

並列プログラミング講習会（試行）ＵＮＩＸ備忘録
 emacsの起動：

emacs 編集ファイル名 ^x^s （＾はcontrol） ：テキストの保存 ^x ^c ： 終了
（ ^z で終了すると、スパコンの負荷が上がる。絶対にしないこと。）
 ^g: 訳がわからなくなったとき。
 ^k : カーソルより行末まで消す。
カ
より行末ま 消す
消した行は、一時的に記憶される。
 ^y^ : ^kで消した行を、現在のカーソルの場所にコピーする。
^kで消した行を 現在のカ ソルの場所にコピ する ^s 文字列 : 文字列の箇所まで移動する。
 ^x^ goto-line
t li : 指定した行まで移動する。
指定した行まで移動する

26

並列プログラミング講習会（試行）ＵＮＩＸ備忘録
 rm ファイル名：


ファイル名のファイルを消す。

rm *~ : test.c~ などの、~がついたバックアップファイルを消す。

 ls : 現在いるフォルダの中身を見る。
 cd

フォルダ名： フォルダに移動する。

 c
cd.. : 一つ上のフォルダに移動。
のフォルダに移動。
 cd ~ ：ホームディレクトリに行く。訳がわからなくなったとき。
 cat

ファイル名： ファイル名の中身を見る make : 実行ファイルを作る（Makefile があるところでしか実行できない）
 makeclean : 実行ファイルを消す。
（clean がMakefileで定義されていないと実行できない）

27

並列プログラミング講習会（試行）役に立つコマンド
実行ファイルが使用するメモリ量を見る$ size -f <実行ファイル名>
テキスト領域1.

.text
text + .data
data + .bss
bss などの情報が見れます
HA8000の場合 １ノードでのHA8000の場合、１ノードでの物理メモリ上限は32GB（実際は28GB程度）です

malloc()関数で動的メモリ確保される場合はsizeコマンドではわかりません。
ただし、動的確保されるメモリ量が１ノードで処理できる物理サイズを超えると、システム
がハングアップしますので、ご注意ください。

28

並列プログラミング講習会（試行）

デ タ領域
データ領域
（初期化データ、.data（初期化データ、.data））
（未初期化
デ
データ、.
データ、
.bss
bss））
ヒープ領域
（動的データ）

スタック領域サンプルプログラムの実行
初めての並列プログラムの実行

29

並列プログラミング講習会（試行）サンプルプログラムについて
 Ｃ言語版

lec2010HA＜サンプルプログラム名＞C.tar
 Ｆｏｒｔｒａｎ版
Ｆ ｔ版
lec2010HA＜サンプルプログラム名＞Ｆ.tar
 上記のファイルが置いてある場所/home/z30082

30

並列プログラミング講習会（試行）並列版Helloプログラムをコンパイルしよう
1.

2
2.
3
3.
4.
5.
6.
31

/home/z30082 にある lec2010HAsamplesC.tar を
自分のディレクトリにコピ する
自分のディレクトリにコピーする$ cp /home/z30082/lec2010HAsamplesC.tar ./lec2010HAsamplesC tar を展開する
lec2010HAsamplesC.tar$ tar xvf lec2010HAsamplesC.tar
Hello フォルダに入る$ cd Hello
ピュアＭＰＩ用のＭａｋefileをコピーする
ア用をする$ cp Makefile_pure Makefile
make する$ make
実行ファイル(hello)ができていることを確認する$ ls
並列プログラミング講習会（試行）バッチ処理とは スパコン環境では、インタラクティブ実行（コマンドラインで実行すること）はできません
実行すること）はできません。
 ジョブはバッチ処理で実行します。
バッチキュー

ジョブの依頼

バッチ処理
システムが
ジ ブを
ジョブを取り出す
出す

実行
ユーザ
32

スパコン
並列プログラミング講習会（試行）NQSキューの設定のしかた


バッチ処理は、Network Queuing System(NQS)で管理されています。

 以下、NQSの主要コマンドを説明します。
以下 NQSの主要コマンドを説明します

ジョブの投入： qsub <ジョブスクリプトファイル名> 自分が投入したジョブの状況確認：
自分が投入したジ ブの状況確認 qstat
 混雑度を見る： qstat -b
 投入ジョブの削除： qdel
del <ジョブID>
 実行中のジョブの削除： qdel -k <ジョブID>

 本演習のNQSキュー名：
本演習のNQSキュ 名 tutorialt t i l

最大１５分までの実行を許可
 最大ノード数は４
最大ノ ド数は４

33並列プログラミング講習会（試行）NQSキューの種類 “qstatq-b” と入力するとキューの種類を見ることができます tutoriall :本講習会時間内で使えるキ本講習会時間内で使えるキュー lecture :本講習会時間外で使えるキュー
本講習会時間外で使えるキ34

並列プログラミング講習会（試行）qstat –bb の実行画面例[t00002@ha8000-2 ~]$ qstat -b
NQS schedule stop time : 2010/07/23 (Fri) 8:55:00 (Remain: 40h 36m 23s)
QUEUE NAMESTATUS TOTAL RUNNING RUNLIMIT QUEUED HELDdebug
AVAILBL 004000tutorial
AVAILBL 001000lectureAVAILBL 0
0
1
0
0
0

使える
キュー名

35

現在
キュー
使えるか に投入
されて
いる
ジョブ
の総数

同時
実行
数の
上限

待たされ
ている
ジョブの
数

並列プログラミング講習会（試行）IN-TRANSIT

システムが
落ちるまで
の残り時間。
これ以上
長い実行
時間を指定
したら、実行されない。

JOBスクリプトサンプルの説明（ピュアＭＰＩ）（C言語 Fortran言語共通）（C言語、Fortran言語共通）#@$-q lecture
キュー名： lecture#@$ N 4#@$-N
ノード数： 4#@$-J T16
16プロセス実行
実行#@$-lm 2GB#@$-lT 0:01:00
実行時間制限： １分
#@$-o
@
OUT
標準出力フ イル名 OUT
標準出力ファイル名：OUT#@$-e ERR
cd $PBS_O_WORKDIR
標準エラーファイル名:ERR
exportt HF_PRUNST_THREADNUM=1
HF PRUNST THREADNUM 1
mpirun ./numaran-pure.sh ./hello >> ./hello-pure.out NQSを投入した
ディレクトリ
MPIジョブを16＊4=64 プロセスで実行する。
へ移動
hello.out に出力結果を上書きする（”>>”）。
※上書きしな※上書きしないで新規作成（“＞”）：ファイルが無いことを確認
新規作成（“＞”）
イ が無
とを確認
36

並列プログラミング講習会（試行）numaran-pure.shの中身（ピュアＭＰＩ）（C言語 Fortran言語共通）（C言語、Fortran言語共通）
MPIのランク番号を取得#!/bin/bash
ソケット内に4プロセスを配置
を 置
MYRANK=$MXMPIMYRANK$MXMPI_ID
ID
MYVAL=$(expr $MYRANK / 4)
ソケット番号を計算
NODE=$(expr $MYVAL % 4)numactl --cpunodebind=$NODE --membind=$NODE $@プロセスの割り当てソケットとメモリ確保先のソケットを指定

37

並列プログラミング講習会（試行）各ＣＰＵの内部構成

ピュアＭＰＩの実行状況（ノード内）

L3

ＭＰＩプロセス

L2

L2

L2

L1
Core
#0

L1
Core
#1

L1
Core
#2

L1
Core
#3

Memory #0

Memory #1

AMD QuadCore
Opteron#0

AMD
Quad Core
Opteron#1

South
BridgeMyrinetQ d
AMD QuadCore
Opteron#2

Q d
AMD QuadCore
Opteron#3

South
Bridgeg

MyrinetM
Memory#2
38

L2

M
Memory#3
並列プログラミング講習会（試行）

MyrinetMyrinetGbE
RAID

並列版Helloプログラムを実行しよう（ピュアMPI）
 このサンプルのＪＯＢスクリプトは

hello-pure.bashです。
です 配布のサンプルでは、キュー名が”lecture”になっています $ emacs hello-pure.bash
hello pure bash
で、“lecture” →“tutorial” に変更してください
39

並列プログラミング講習会（試行）並列版Helloプログラムを実行しよう（ピュアMPI）
1.
2.
3
3.
4
4.
5.

40

Helloフォルダ中で以下を実行する$ qsub hello-pure.bash
自分の導入されたジョブを確認する
分 導
認す$ qstat
実行が終了すると、以下のファイルが生成される
hello-pure.out
上記ファイルの中身を見てみる$ cat hello-pure.out“HHelloll parallel
ll l world!
ld!”がが、16プロセス*4ノード=64
16プロセス*4ノ ド 64
表示されていたら成功。
並列プログラミング講習会（試行）NQSに関連するTIPS
NQS ジョブが実行されるとき、およびNQSジョブが終了したとき には電子メールが送られてくる
終了したとき、には電子メールが送られてくる。
このメールを読むときは、以下のようにする。
> mailなお、コマンドは以下のとおりです。











41

番号： メールの表示
番号
メ ルの表示 Ex.
E >２
h : メール一覧表示
s : ファイルへ保存 Ex.
E > s filefild: メールの削除 Ex. > d １
x: コマンドの終了（何もしない）
q: コマンドの終了（mboxへ保存）
ほかのメ ルアドレスに飛ばしたい場合は、ほかのメールアドレスに飛ばしたい場合はha8000-１.cc.u-tokyo.ac.jp で.forward ファイル中にアドレスを記述
並列プログラミング講習会（試行）サンプルプログラムについて
 Ｃ言語版

lec2010HA＜サンプルプログラム名＞C.tar
 Ｆｏｒｔｒａｎ版
Ｆ ｔ版
lec2010HA＜サンプルプログラム名＞Ｆ.tar
 上記のファイルが置いてある場所/home/z30082

42

並列プログラミング講習会（試行）各ＣＰＵの内部構成

ピュアＭＰＩの実行状況（ノード内）

L3

ＭＰＩプロセス

L2

L2

L2

L1
Core
#0

L1
Core
#1

L1
Core
#2

L1
Core
#3

Memory #0

Memory #1

AMD QuadCore
Opteron#0

AMD
Quad Core
Opteron#1

South
BridgeMyrinetQ d
AMD QuadCore
Opteron#2

Q d
AMD QuadCore
Opteron#3

South
Bridgeg

MyrinetM
Memory#2
43

L2

M
Memory#3
並列プログラミング講習会（試行）

MyrinetMyrinetGbE
RAID

並列版Helloプログラムを実行しよう（ピュアMPI）
 このサンプルのＪＯＢスクリプトは

hello-pure.bashです。
です 配布のサンプルでは、キュー名が”lecture”になっています $ emacs hello-pure.bash
hello pure bash
で、“lecture” →“tutorial” に変更してください
44

並列プログラミング講習会（試行）並列版Helloプログラムを実行しよう（ピュアMPI）
1.
2.
3
3.
4
4.
5.

45

Helloフォルダ中で以下を実行する$ qsub hello-pure.bash
自分の導入されたジョブを確認する
分 導
認す$ qstat
実行が終了すると、以下のファイルが生成される
hello-pure.out
上記ファイルの中身を見てみる$ cat hello-pure.out“HHelloll parallel
ll l world!
ld!”がが、16プロセス*4ノード=64
16プロセス*4ノ ド 64
表示されていたら成功。
並列プログラミング講習会（試行）並列版Helloプログラムをコンパイルしよう
1.
2.

3.

4.

46

ハイブリッドＭＰＩ用の Makefile をコピーする。
$ cp Makefile_hy4 Makefile
make する。
す 。
$ make clean$ make
実行ファイル（hello）ができていることを確認する。
$ ls
ＪＯＢスクリプト中（hello-hy4.bash）のキュー名を変更する “l t ” → “tutorial”に変更する。
変更する。“lecture”“t t i l”に変更する$ emacs hello-hy4.bash
並列プログラミング講習会（試行）並列版Helloプログラムを実行しよう（ハイブリッドMPI）
1.
2.
3
3.
4
4.
5.

47

Helloフォルダ中で以下を実行する$ qsub hello-hy4.bash
自分の導入されたジョブを確認する
分 導
認す$ qstat
実行が終了すると、以下のファイルが生成される
hello-hy4.out
上記ファイルの中身を見てみる$ cat hello-hy4.out“HHelloll parallel
ll l world!
ld!”がが、4プロセス*4ノード=16
4プロセス*4ノ ド 16
表示されていたら成功。
並列プログラミング講習会（試行）JOBスクリプトサンプルの説明（ハイブリッドＭＰＩ）（C言語 Fortran言語共通）（C言語、Fortran言語共通）#@$-q lecture
キュー名： lecture#@$ N 4#@$-N
ノード数： 4#@$-J T4
4プロセス実行
実行#@$-lm#@$lm 2GB#@$-lT 0:01:00
実行時間制限： １分#@$-o OUT
標準出力ファイル名：OUT#@$-e ERR
cd $PBS_O_WORKDIR
標準エラーファイル名:ERR
export HF_PRUNST_THREADNUM=4
HF PRUNST THREADNUM=4NQSを投入mpirun ./numaran-4.sh ./hello >> ./hello-hy4.out
したディレクトリへ移動
MPIジョブを16プロセスで実行する。
ジ ブを プ セ
実行する
hello.out に出力結果を上書きする（”>>”）。
※上書きしないで新規作成（“＞”）：ファイルが無いことを確認
48

スレッド数（起動プロ
セス数を含む）の指定

並列プログラミング講習会（試行）numaran-4.shの中身（ピュアＭＰＩ）（C言語 Fortran言語共通）（C言語、Fortran言語共通）
MPIのランク番号を取得#!/bin/bash
ソケット内に1プロセスを配置
MYRANK=$MXMPIMYRANK$MXMPI_ID
ID
ソケット番号を計算
NODE=$(expr $MYRANK % 4)numactl --cpunodebind=$NODE --membind=$NODE $@プロセスの割り当てソケットとメモリ確保先のソケットを
リ確保先のソケットを指定

49

並列プログラミング講習会（試行）各ＣＰＵの内部構成

ハイブリッドＭＰＩの実行状況（ノード内）

L3

ＭＰＩプロセス
スレッド

L2

L2

L2

L1
Core
#0

L1
Core
#1

L1
Core
#2

L1
Core
#3

Memory #0

Memory #1

AMD QuadCore
Opteron#0

AMD
Quad Core
Opteron#1

South
BridgeMyrinetQ d
AMD QuadCore
Opteron#2

Q d
AMD QuadCore
Opteron#3

South
Bridgeg

MyrinetM
Memory#2
50

L2M
Memory#3
並列プログラミング講習会（試行）

MyrinetMyrinetGbE
RAID

並列版Helloプログラムの説明（C言語）#include <stdio.h>#include <mpi<mpi.h>h>

このプログラムは、全ＰＥで起動される

void main(int argc, char* argv[]) {int
int

MPIの初期化

myid, numprocs;
ierr, rc;

自分のID番号を取得：各ＰＥで値は異なる
各
値 異なる

ierr = MPI_Init(&argc, &argv);
ierr = MPI_Comm_rank(MPI_COMM_WORLD, &myid);ierr = MPI_Comm_size(MPI_COMM_WORLD,MPI Comm size(MPI COMM WORLD, &numprocs);
printf("Hello parallel world! Myid:%d ¥n", myid);rc = MPI_Finalize();exit(0);}51

MPIの終了全体のプロセッサ台数を取得：各ＰＥで値は同じ（演習環境では
６４、もしくは１６）並列プログラミング講習会（試行）並列版Helloプログラムの説明（Fortran言語）
このプログラムは、全ＰＥで起動される

program main
include 'mpif'mpif.h'h'common /mpienv/myid,numprocs
integer myid,
myid numprocs
integer ierr

MPIの初期化
自分のID番号を取得：各ＰＥで値は異なる

call MPIMPI_INIT(ierr)
INIT(ierr)call MPI_COMM_RANK(MPI_COMM_WORLD, myid, ierr)call MPI_COMM_SIZE(MPI_COMM_WORLD, numprocs, ierr)print *, "Hello parallel world! Myid:", myid
call MPI_FINALIZE(ierr)
( )
stop
end
52

MPIの終了全体のプロセッサ台数を取得：各ＰＥで値は同じ
各ＰＥで値は同じ（演習環境では
６４、もしくは１６）並列プログラミング講習会（試行）時間計測方法（C言語）
double t0, t1, t2, t_w;..
ierr = MPI_Barrier(MPI_COMM_WORLD);
t1 = MPI_Wtime();

バリア同期後時間を習得し保存＜ここに測定したいプログラムを書く＞
ierr = MPI_Barrier(MPI_COMM_WORLD);
t2 = MPI_Wtime();
MPI Wtime();t0 = t2 - t1;
ierr = MPI_Reduce(&t0,
MPI Reduce(&t0 &t
&t_w,
w 1,
1
MPI_DOUBLE,MPI_MAX, 0,
MPI_COMM_WORLD);

53

各プロセッサーで、ｔ0の値は異なる。
この場合は、最も遅いものの場合は 最も遅 も
値をプロセッサ０番が受け取る

並列プログラミング講習会（試行）時間計測方法（Fortran言語）
double precision t0, t1, t2, t_w
double precision MPI_WTIME..
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
t1 = MPI_WTIME(ierr)

バリア同期後時間を習得し保存＜ここに測定したいプログラムを書く＞
call MPIMPI_BARRIER(MPI_COMM_WORLD,BARRIER(MPI COMM WORLD ierr)t2 = MPI_WTIME(ierr)
t0 = t2 - t1
call MPI_REDUCE(t0, t_w, 1,&
MPI_DOUBLE_PRECISION,&MPI_MAX,MPI MAX 0,
0 MPI_COMM_WORLD,MPI COMM WORLD ierr)54

各プロセッサーで、ｔ0の値は異なる。
この場合は、最も遅いも場合最も遅 ものの値をプロセッサ０番が受け取る

並列プログラミング講習会（試行）任意のプロセス数での実行方法（東大センタ限定） 以下のように、
以下のように、*.bashファイルを修正してください.bashファイルを修正してください ２プロセス／ノードでの実行#@$-JJ T2#@$
 ４プロセス／ノードでの実行#@$-J T4
 pプロセス／ノードでの実行#@$-J Tp
 実習環境では、１～１６プロセス／ノードまで
プ

しかないことに注意
55

並列プログラミング講習会（試行）サンプルプログラムの説明










Hello/


並列版H ll プログラム
並列版Helloプログラム



hello-pure.bash, hello-hy4.bash : NQSジョブスクリプトファイル

Cpi/


円周率計算プログラム

cpi-pure.bash NQSジョブスクリプトファイル

Wa1/


逐次転送方式による総和演算

wa1-pure.bash NQSジョブスクリプトファイル

Wa2/


二分木通信方式による総和演算

wa2-pure.bashp
NQSジョブスクリプトファイル
Q

Cpi_m/


円周率計算プログラムに時間計測ルーチンを追加したもの

cpi_m-pure.bashib h NQSジョブスクリプトファイル
NQSジ ブスクリプトフ イル

56並列プログラミング講習会（試行）演習課題
逐次転送方式のプログラムを実行

1.


Wa1 のプログラム

二分木通信方式のプログラムを実行2.


Wa2のプログラム

時間計測プログラムを実行

3.


Cpi_mのプログラム

プロセス数を変化させて、サンプルプログラムを実行
Helloプログラムを、以下のように改良

4.
5.


57

MPI_Sendを用いて、PE0からＣｈａｒ型のデータ
を用 て、
から
型のデ タ
“Hello World!!”を、その他のPEに送信する
その他のPEでは、MPI_Recvで受信して表示する
並列プログラミング講習会（試行）総和演算プログラム（逐次転送方式）


各プロセスが所有するデータを、全プロセスで加算し、あるプ セス１つが結果を所有する演算を考える。
あるプロセス１つが結果を所有する演算を考える。

 素朴な方法（逐次転送方式）（0番でなければ）左隣のプロセスからデータを受信する；（
番でなければ）左隣のらデ タを受信する；
左隣のプロセスからデータが来ていたら；

1.
2.
1.
2.
3.
4.

受信する;
＜自分のデータ＞と＜受信データ＞を加算する;
自デデ
算（63番でなければ）右隣のプロセスに＜2の加算した結果を＞送信する;
処理を終了する;

実装上の注意

左隣りとは、（myid-１）のIDをもつプロセス

右隣りとは、（myid+１）のIDをもつプロセス58

myid=0のプロセスは、左隣りはないので、受信しない

myid=p-１のプロセスは、右隣りはないので、送信しない
y p のプは、右隣りはな ので、送信しな
並列プログラミング講習会（試行）バケツリレー方式による加算
所有データ

所有データ

０

１

０

所有データ

所有データ

２

３

１

CPU0

CPU１
CPU
１
送信

０ ＋ １ ＝ １

３
CPU２
CPU
２

送信

CPU３
CPU
３
送信

１ ＋ ２ ＝ ３

３ ＋ ３ ＝ ６

最終結果

59

並列プログラミング講習会（試行）１対１通信利用例（逐次転送方式 C言語）
（逐次転送方式、C言語）
受信用システム配列の確保
受信用シ
テ 配列 確保
void main(int argc,
argc char
char* argv[]) {MPI_Status istatus;
自分より一つ少ない….
ID番号（myid-1）から、
ydsendbuf = myid;
double型データ１つをdrecvbuf = 0.0;
受信しdrecvbuf変数に
if (myid != 0) {ierr = MPI_Recv(&drecvbuf,(, 1,, MPI_DOUBLE,, myid-1,y, 0,, 代入
MPI_COMM_WORLD, &istatus);}
自分より一つ多い
dsendbuf = dsendbuf + drecvbuf;
ID番号（myid+1）にID番号（myid+1）に、if (myid != nprocs-1) {ierr = MPI_Send(&dsendbuf, 1, MPI_DOUBLE, myid+1, 0, dsendbuf変数に入っているdouble型データ
MPI_COMM_WORLD);
１つを送信}
if (myid == nprocs-1) printf ("Total = %4.2lf ¥n", dsendbuf);….
}
60

並列プログラミング講習会（試行）１対１通信利用例（逐次転送方式 Fortran言語）（逐次転送方式、Fortran言語）
program main
受信用システム配列の確保
受信用シ
テ 配列 確保
integerege istatus(MPI
s a us(_S
STATUS
US_S
SIZE))….
dsendbuf = myid
自分より一つ少ない
drecvbuf = 0.0
ID番号（myid-1）から、
yif ((myid
id .ne. 0) thenthdouble型データ１つを
call MPI_RECV(drecvbuf, 1, MPI_DOUBLE_PRECISION,
受信しdrecvbuf変数に&
myid-1, 0, MPI_COMM_WORLD, istatus, ierr)endif
代入
dsendbuf = dsendbuf + drecvbuf
if (myid .ne. numprocs-1) thencall MPI_SEND(dsendbuf, 1, MPI_DOUBLE_PRECISION,&
myid+1,
id+1 0,0 MPI_COMM_WORLD,MPI COMM WORLD ierr)i )自分より一つ多い
endif
ID番号（myid+1）に、if (myid .eq. numprocs-1) thendsendbuf変数に変数
print *,, "Total = ",, dsendbuf
pendif
入っているdouble型….
データ１つを送信
stop
end
61

並列プログラミング講習会（試行）総和演算プログラム（二分木通信方式）
 二分木通信方式

k = １;for (i=0; i < log2(nprocs); i++)iff ( (myid( d & k) == k)

1.
2.
3.




else

4.



62

(myid – k)番 プロセス からデータを受信；
自分のデータと自分のデタと、受信デ
受信データを加算する；
タを加算する；k = k * 2;
(myid + k)番 プロセス に、データを転送する；
処理を終了する；

並列プログラミング講習会（試行）総和演算プログラム（二分木通信方式）
３段目＝ｌog2（８）段目段目
g （ ）段目

3

1

3

２段目
0

1段目

1

2

7
5

3

4

7

5

6

7

7
7

3

0
63

5

3

1
1

2

3

4

7
5

並列プログラミング講習会（試行）67総和演算プログラム（二分木通信方式）
 実装上の工夫

要点： プロセス番号の２進数表記の情報を利用する
要点プ
番号
進数表記 情報を利用する
 第i段において、受信するプロセスの条件は、以下で書ける：
myid & k が k と一致と 致 ここで、k = 2^(i-１) 。
 つまり、プロセス番号の２進数表記で右からi番目のビットが立っている

プロセスが、送信することにする

また、送信元のプロセス番号は、以下で書ける：
myid + k つまり 、通信が成立するPE番号の間隔は２＾(i-１)

64←二分木なので送信プ セ に い は 上記の逆が成り立
送信プロセスについては、上記の逆が成り立つ。

並列プログラミング講習会（試行）総和演算プログラム（二分木通信方式）
 逐次転送方式の通信回数

明らかに
明らかに、nprocs－１
１回

 二分木通信方式の通信回数

見積もりの前提
 各段で行われる通信は、完全に並列で行われる（通信の衝突は発生しない）

段数の分の通信回数となる
 つまり、log2(nprocs)
g ( p
)回


 両者の通信回数の比較

プロセッサ台数が増すと、通信回数の差（＝実行時間）がッサ台数が増すと、通信回数の差（ 実行時間）がとても大きくなる １０２４構成では、１０２３回 対 １０回！




65

でも、必ずしも二分木通信方式がよいとは限らない（通信衝突の多発）
並列プログラミング講習会（試行）並列プログラミングの基礎（座学）
東京大学情報基盤センター 特任准教授 片桐孝洋

２０１０年９月２７日（月）１４：００－１５：００
66

並列プログラミング講習会（試行）本講義の流れ
1.
2.
3.
4.
5
5.
6.
7
7.
8.

67

東大スーパーコンピュータの概略
並列プログラミングの基礎
性能評価指標基礎的なMPI関数
データ分散方式
デ
タ分散方式
ベクトルどうしの演算
ベクトル 行列積
ベクトル-行列積
リダクション演算

並列プログラミング講習会（試行）東大ス パ コンピュ タの概略
東大スーパーコンピュータの概略

68

並列プログラミング講習会（試行）東京大学情報基盤センター スパコン
HITACHI SR11000 model J2

T2K東大（HA8000クラスタシステム）

Total Peak p
performance: 18.8 TFLOPS
Total number of nodes: 128
Total memory: 16384 GB
Peak performance per node : 147.2 GFLOPSMain memory per node
: 128 GBDisk capacity: 94.2 TB
IBM POWER5+ 2.3GHz

Total Peak p
performance: 140 TFLOPS
Total number of nodes: 952
Total memory: 32000 GB
Peak performance per node : 147.2 GFLOPSMain memory per node
: 32 GBGB, 128 GB
Disk capacity: 1 PB
AMD Quad Core Opteron 2.3GHz
ノ ド製品名 HITACHI HA8000-tc/RS425
ノード製品名：HITACHI
HA8000 tc/RS42569

並列プログラミング講習会（試行）69T2Kオープンスパコン東大版
HA8000-tc/RS425クラスタ群
HA8000tc/RS425クラスタ群
全体仕様

ノード単体仕様CPU
ソ ケ ッ ト 数 (CPU数)
理論演算性能主記憶容量
ディスク

AMD社製Quad Core Opteron 2.3 GHz
4 (16)

総理論演算性能
総主記憶容量
計算ノード数

140.1344 TFlops
31.25 TByte
948

147.2 GFlopsl
32 GByte (936ノード)、128GByte(16ノード)
250 GByte

ログインノード数4

計算ノード群タイプB
プ (36 + 16 ノード)
＆ログインノード群(4ノード)
計算

外部部接続ルータ
タ

ログイン
ログイン
ノ ド
ノード
ログイン
ノード
ログイン
ノード
ノード

70

ノード
計算
計算
ノード
計算
ノード
ノード

ネ ッ ト ワ ー ク 性 能 （ 両 タイプA(512,128ノード)：5 GByte/sec (Myrinet-10G x 4)
方 向 物 理 デ ー タ 転 送 性 タイプB(256,36,16)： 2.5 GByte/sec (Myrinet-10G x 2)
能）ファイルシステム容量
1 PByteOS

RedHat Enterprise Linux 5

バッチジョブシステム

現有バッチジョブシステムと同様の機能が提供されます

コンパイラ

日立製作所社製最適化Fortran, C/C++ (OpenMP有)通信ライブラリ

MPI通信ライブラリ(MPICH-MX)

数値計算ライブラリ

MSL2
MATRIX/MPP
MATRIX/MPP/SSS
BLAS, LAPACK, ScaLAPACK
Gaussian03

アプリケーション

計算ノード群 タイプB
(256ノード)計算
ノード 計算
ノード
計算
計算
計算
ノード
計算
ノード
ノ 計算
ノード
ド 計算
ノード ノードノード

計算ノード群 タイプA
(128ノード)
計算
ノード
計算
計算
ノード
計算
ノード
ド
ノード

並列プログラミング講習会（試行）ストレージシステム
計算ノード群 タイプA
(512ノード) 計算
1 PByte
計算
計算
ノード
計算
ノード
ノード
ノード
計算
計算
計算
計算
ノード
計算
計算
ノード
計算
計算
ノード
ノ
ノード
ノ
計算
ド
計算
ノド ド
計算
ノ
ノード
ド ノード
ノード
計算
ノード
ノード
ノード
ノード
ノード
70

タイプＡ
ノード数
ネットワーク性能

512, 128 （２つの群から構成）
5 GB/sec (片方向）F ll bisection
Full
bi
ti
b d idth
bandwidth
CPU
AMD 社 製 Quad Core Opteron
2.3 GHzソ ケ ッ ト 数 (CPU 4 (16)数)
理論演算性能147.2 GFlops
主記憶容量
32 GByte
ディスク
250 GByte
512ノード71

並列プログラミング講習会（試行）128ノード

71

タイプB
ノード数

256, 36, 16 （3つの群から構成）

ネットワーク性能
ネットワ
ク性能
CPU

2.5 GB/sec (片方向）Full bisection bandwidth
AMD社製Quad Core Opteron 2.3 GHz

ソケット数(CPU数)

4 (16)

理論演算性能主記憶容量

147.2 GFlopsy(16ノードは128GByte)y32 GByte

ディスク

250 GByte
256ノード

36ノード

72

並列プログラミング講習会（試行）なお、ログインノード
４台の構成は３６ノードと同様です

16ノード

72

T2KオープンスパコンCPU緒元（AMD Qu
（
Quad
d Co
Core
e OpteOpterono 2.3GHz）
2 3G ）：東大、京大、筑波大に設置のスパコン
項目

値

L1キャッシュ64 Kbytes (命令、データ、双方は分離)
2 Way Associativity（ライトバック、3サイクル）
キャッシュライン：64 bytes、LRU置換
Full Associativity
2Mbyte Page:8エントリ, 4Kbyte Page:32エントリ
512 Kbytes (2300 MHz)
2048 Kbytes
3 Way
3 Way (整数、アドレス生成、浮動小数点(整数 アドレス生成 浮動小数点 )MMX, SSE, SSE2
乱発行乱終了(Out-of-order) 整数、浮動小数点
 レ ガ シ ー モ ー ド ： 汎 用 32bit: 8 個 、 128bit-XMM:8
128bit XMM:8 個 、64bit MMX:8個（x87互換用:8個と同一）
 64bit モ ー ド ： 汎 用 64bit:16 個 、 128bit-XMM:16 個 、64bit MMX:8個（x87互換用:8個と同MMX:8個（x87互換用:8個と同一））

L1データTLB
L2キャッシュL3キャッシュ命令デコード演算実行
SIMD命令実行
命令実行
レジスタ

システムバス731000 MHz

並列プログラミング講習会（試行）各ＣＰＵの内部構成

物理ノードの構成（タイプＡ群）

L3
L2

L2

L2

L1
Core
#0

L1
Core
#1

L1
Core
#2

L1
Core
#3

Memory #0

Memory #1

AMD QuadCore
Opteron#0

AMD
Quad Core
Opteron#1

South
BridgeMyrinetQ d
AMD QuadCore
Opteron#2

Q d
AMD QuadCore
Opteron#3

South
Bridgeg

MyrinetM
Memory#2
74

L2

M
Memory#3
並列プログラミング講習会（試行）

MyrinetMyrinetGbE
RAID

ノード構成（東大、タイプＡ群）

各ＣＰＵの内部構成
L3

アクセス
速い

共有メモリ
でアクセス
時間が
不均一
（ccNUMA）

L2

L2

L2

L1
Core
#0

L1
Core
#1

L1
Core
#2

L1
Core
#3

Memory #0

Memory #1

AMD QuadCore
Opteron#0

AMD
Quad Core
Opteron#1

South
BridgeMyrinetQ d
AMD QuadCore
Opteron#3

South
Bridgeg

MyrinetMyrinet

アクセス
遅い

Q d
AMD QuadCore
Opteron#2

M
Memory#2
75

L2M
Memory#3
並列プログラミング講習会（試行）

MyrinetGbE
RAID

ピュアＭＰＩ：MPIプロセス数 ＝ コア数
ハイブリッドＭＰＩ：（１つの例）
MPIプロセス数：４、OpenMPスレッド数：４

各ＣＰＵの内部構成
L3
L2

L2

L2

L2

L1
Core
#0

L1
Core
#1

L1
Core
#2

L1
Core
#3

Memory #0

Memory #1

AMD QuadCore
Opteron#0
OpenMPp

AMD
Quad Core
Opteron#1

South
BridgeMyrinetQ d
AMD QuadCore
Opteron#3

South
Bridgeg

Myrinet

各スレッド
のローカル
メモリ
アクセス
の推進

スレッド
MPI
プロセス

Q d
AMD QuadCore
Opteron#2

M
Memory#2
76

M
Memory#3
並列プログラミング講習会（試行）

MyrinetMyrinetGbE
RAID

スパコンでの全体メモリ構成図
コア
０

コア
１

コア
２

コア
３

Ｌ１

Ｌ１

Ｌ１

Ｌ１

Ｌ２

Ｌ２

Ｌ２

Ｌ２

…

コア
１２

コア
１３

コア
１４

コア
１５

Ｌ１

Ｌ１

Ｌ１

Ｌ１

Ｌ２

Ｌ２

Ｌ２

Ｌ２

Ｌ３

…

コア
０

コア
１

コア
２

コア
３

Ｌ１

Ｌ１

Ｌ１

Ｌ１

Ｌ２

Ｌ２

Ｌ２

Ｌ２

Ｌ３

コア
１

コア
２

コア
３

Ｌ１

Ｌ１

Ｌ１

Ｌ１

Ｌ２

Ｌ２

Ｌ２

Ｌ２

…

Ｌ３

コア
１４

コア
１５

Ｌ１

Ｌ１

Ｌ１

Ｌ１

Ｌ２

Ｌ２

Ｌ２

Ｌ２

Ｌ３
メインメモリ

コア
１２

コア
１３

コア
１４

コア
１５

コア
０

コア
１

コア
２

コア
３

Ｌ１

Ｌ１

Ｌ１

Ｌ１

Ｌ１

Ｌ１

Ｌ１

Ｌ１

Ｌ２

Ｌ２

Ｌ２

Ｌ２

Ｌ２

Ｌ２

Ｌ２

Ｌ２

Ｌ３

… …

…

Ｌ３

メインメモリ

77

コア
１３

Ｌ３

メインメモリ

コア
０

…

コア
１２

コア
１２

コア
１３

コア
１４

コア
１５

Ｌ１

Ｌ１

Ｌ１

Ｌ１

Ｌ２

Ｌ２

Ｌ２

Ｌ２

Ｌ３
メインメモリ

メモリが多段に階層化（L1、L2、L3、分散メモリ）

並列プログラミング講習会（試行）高速ネットワーク（５Ｇバイト／秒
×双方向）（タイプＡ群）

Ｔ２Ｋオープンスパコンの詳細情報
 以下のページをご参照ください
以下のジを 参照くださ 利用申請方法
 運営体系
 料金体系
 利用の手引などがご覧になれます。
などがご覧になれます
http://www.cc.u-tokyo.ac.jp/service/ha8000/
78

並列プログラミング講習会（試行）東大情報基盤センタHA8000クラスタシステムの料金表（2008年6月1日制定 ） パーソナルコース（年間）

コース１： １１０，０００円 ： ４ノード
 コース２： １５０，０００円 ： ８ノード
 コース３：コ ス３ ２２０，０００円
２２０ ０００円 ：１６ノード
１６ノ ド
 コース４： ３００，０００円 ：３２ノード
 コ
コース５：
ス５： ４４０，０００円 ：６４ノ
：６４ノードド 月に1度、ひとつ上のコースで実行可能。コース５は、２５６ノード。
 専用キュー
専用キ

１，０００，０００円 ：８ノード

 ノード固定

１，５００，０００円 ：８ノード

 以上は、「大学等」料金79

並列プログラミング講習会（試行）並列プログラミングの基礎

80

並列プログラミング講習会（試行）並列プログラミングとは何か？
 逐次実行のプログラム（実行時間T）を、ｐ台の計算機を使って、T / ｐ にすること.
T
T / ｐ
 素人考えでは自明。
 実際は、できるかどうかは、対象処理の内容
実際は できるかどうかは 対象処理の内容（アルゴリズム）で 大きく 難しさが違う
アルゴリズム上、絶対に並列化できない部分の存在
 通信のためのオーバヘッドの存在 通信立ち上がり時間
 データ転送時間81

並列プログラミング講習会（試行）並列と並行 並列（Parallel）

物理的に並列（時間的に独立）
 ある時間に実行されるものは多数
T 並行（Concurrent）

論理的に並列（時間的に依存）
 ある時間に実行されるものは１つ（＝１プロセッサで実行）T

時分割多重、疑似並列
 OSによるプロセス実行スケジューリング（ラウンドロビン方式）
82

並列プログラミング講習会（試行）並列計算機の分類
 Michael J. Flynn教授（スタンフォード大）の分類（１９６６）
 単一命令・単一データ流（SISD, Single Instruction Single Data Stream）
 単一命令・複数データ流（SIMD, Single Instruction Multiple Data Stream）
 複数命令・単一データ流（MISD, Multiple Instruction Single Data Stream）
 複数命令・複数データ流（MIMD, Multiple Instruction Multiple Data Stream）

83

並列プログラミング講習会（試行）並列計算機のメモリ型による分類
1.

共有メモリ型（SMP、
Symmetric Multiprocessor）

2.

分散メモリ型（メッセージパッシング）

33.

分散共有メモリ型（DSM、
Distributed Shared Memory）
84

並列プログラミング講習会（試行）並列計算機のメモリ型による分類

4.

共有・非対称メモリ型（ccNUMA、
Cache Coherent NonUniform Memory Access）

85

並列プログラミング講習会（試行）並列プログラミングのモデル
 実際の並列プログラムの挙動はMIMD
 アルゴリズムを考えるときは＜SIMDが基本＞
 複雑な挙動は理解できないので

86

並列プログラミング講習会（試行）並列プログラミングのモデル
 MIMD上での並列プログラミングのモデル
1.

SPMD（Single Program Multiple Data）
 １つの共通のプログラムが、並列処理開始時に、
全プロセッサ上で起動する
 MPI（バージョン１）のモデル

2.

Master / Worker（Master/ Slave））（
 １つのプロセス（Master）が、複数のプロセス（Worker）を管理（生成 消去）する
管理（生成、消去）する。

87

並列プログラミング講習会（試行）並列プログラムの種類
 マルチプロセス

MPI （Message（
Passing Interface））
 HPF （High Performance Fortran）


自動並列化F t コンパイラ
自動並列化Fortranコンパイラ

ユーザがデータ分割方法を明示的に記述

 マルチスレッド
ルチスレッドプロセスとスレッドの違い•メモリを意識するかどうかの違い•別メモリは、「プロセス」、」•同一メモリは、「スレッド」

Pｔhread (POSIX スレッド)
 Solaris Thread ((Sun Solaris OS用)用) NT thread (Windows NT系、Windows95以降)


スレッドの Fork（分離） と Join（融合） を明示的に記述

 JavaJ

言語仕様としてスレッドを規定

 Ope
Open MP

88

ユーザが並列化指示行を記述
並列プログラミング講習会（試行）並列処理の実行形態（１） データ並列

データを分割することで並列化する。
 データの操作（＝演算）は同一となる。
 データ並列の例：行列－行列積


 1 2 3   9 8 7   1 * 9  2 * 6  3 * 3 1 * 8  2 * 5  3 * 2 1 * 7  2 * 4  3 *1  
 ＝  4 * 9  5 * 6  6 * 3 4 * 8  5 * 5  6 * 2 4 * 7  5 * 4  6 *1 4 5 6   6 5 4   7 * 9  8 * 6  9 * 3 7 * 8  8 * 5  9 * 2 7 * 7  8 * 4  9 *1 
7 8 9 3 2 1 

 

●並列化

全CPUで共有

CPU０  1

2 3  9 8 7 
 1 * 9  2 * 6  3 * 3 1 * 8  2 * 5  3 * 2 1 * 7  2 * 4  3 *1  
 ＝  4 * 9  5 * 6  6 * 3 4 * 8  5 * 5  6 * 2 4 * 7  5 * 4  6 *1CPU１  4 5 6   6 5 4  7 * 9  8 * 6  9 * 3 7 * 8  8 * 5  9 * 2 7 * 7  8 * 4  9 *1 
CPU２  7 8 9   3 2 1   並列に計算：初期データは異なるが演算は同一並列に計算 初期デ タは異なるが演算は同
89

並列プログラミング講習会（試行）並列処理の実行形態（２） タスク並列タスク（ジョブ）を分割することで並列化する。
タスク（ジブ）を分割する とで並列化する
 データの操作（＝演算）は異なるかもしれない。
 タスク並列の例：カレーを作る
タスク並列の例 カレ を作る

仕事１：野菜を切る
 仕事２：肉を切る
 仕事３：水を沸騰させる
 仕事４：野菜・肉を入れて煮込む
仕事 カレ
を入れる
 仕事５：カレールゥを入れる


●並列化

仕事１
仕事２

仕事４

仕事３
仕事
90

仕事５

時間
並列プログラミング講習会（試行）性能評価指標
並列化の尺度

91

並列プログラミング講習会（試行）性能評価指標－台数効果
 台数効果SP  TS / TP (0  S p)
TP ：P台での実行時間
 TS ：逐次の実行時間、
 P台用いてSP  P のとき、理想的な（ideal）速度向上
 P台用いてSP  P のとき、スーパリニア・スピードアップ式： 主な原因は、並列化により、データアクセスが局所化されて、主な並
デ タ
が 所 されキャッシュヒット率が向上することによる高速化
 並列化効率

式： EP  SP / P 100 (0  Ep) [%] 飽和性能速度向上の限界 Saturation、「さちる」
Saturation 「さちる」


92

Ｐ
並列プログラミング講習会（試行）アムダールの法則
K とする。
そのうち 並列化ができる割合を α とする。
そのうち、並列化ができる割合をとする
 このとき、台数効果は以下のようになる。
 逐次実行時間をSP  K /(K / P  K(1)) 1/( / P  (1))  1/((1/ P 1) 1) 上記の式から、たとえ無限大の数のプロセッサを使っても（Pも（P→∞））、台数効果は、高々台数効果は 高々 １／（１１／（１－α）
α） である。
（アムダールの法則）
93

全体の９０％が並列化できたとしても、無限大の数のプロセッ
サをつかっても、 １/(１-0.9) = １0 倍 にしかならない！
→高性能を達成するためには、少しでも並列化効率を上げる
実装をすることがとても重要である
並列プログラミング講習会（試行）アムダールの法則の直観例
並列化できない部分（１ブロック） 並列化できる部分（８ブロック）●逐次実行＝88.8%が並列化可能%が並列化可能●並列実行（４並列）

9/3 3倍
9/3=3倍●並列実行（８並列）

9/2=4 5倍 ≠ 6倍
9/2=4.5倍

94

並列プログラミング講習会（試行）１ 基本演算１.
 逐次処理では、「データ構造」が重要
 並列処理においては、「データ分散方法」が重要

になる！
1.

2.
3
3.
4.

各PEの「演算負荷」を均等にする
ロード・バランシング： 並列処理の基本操作の一つ

粒度調整
各PEの「利用メモリ量」を均等にする
演算に伴う通信時間を短縮する
各PEの「データ・アクセスパターン」を高速な方式にする（＝逐次処理におけるデータ構造と同じ）

 行列データの分散方法＜次元レベル＞： １次元分散方式、２次元分散方式
 ＜分割レベル＞： ブロック分割方式、サイクリック（循環）分割方式
95

並列プログラミング講習会（試行）１.1.1
N/4行

PE=0

N/4行

PE 1
PE=1

N/4行

PE=2

N/4行

PE=3

１次元分散
•（行方向)（行方向) ブロック分割方式•(Block, *) 分散方式

N列
1行•(行方向) サイクリック分割方式•(Cyclic, *) 分散方式

2行•(行方向)ブロック・サイクリック分割方式•(Cyclic(2), *) 分散方式

この例の「２」： ＜ブロック幅＞とよぶ
96

並列プログラミング講習会（試行）2次元分散

１.1.2
N/2

N/2

N/2

PE=0 PE=1

N/2

PE=2 PE=3•ブロック・ブロック分割方式•(Block, Block)分散方式•サイクリック・サイクリック分割方式•(Cyclic, Cyclic)分散方式

0

1

0

1

0

1

0

1

2

3

2

3

2

3

2

3

0

1

0

1

0

1

0

1

2

3

2

3

2

3

2

3

0

1

0

1

0

1

0

1

2

3

2

3

2

3

2

3

0

0

1

1

0

0

1

1

0

1

0

1

0

1

0

1

0

0

1

1

0

0

1

1

2

3

2

3

2

3

2

3

2

2

3

3

2

2

3

3

2

2

3

3

2

2

3

3

0

0

1

1

0

0

1

1

0

0

1

1

0

0

1

1

2

2

3

3

2

2

3

3

2

2

3

3

2

2

3

3

•二次元ブロック・サイクリック分割方式•(Cyclic(2), Cyclic(2))分散方式

並列プログラミング講習会（試行）97１ 2 ベクトルどうしの演算
１.2 以下の演算z  a x y

ここで、αはスカラ、ｚ、ｘ、ｙ はベクトル

 どのようなデータ分散方式でも並列処理が可能

ただし、スカラ α は全PEで所有する。
 ベクトルはO(ｎ)のメモリ領域が必要なのに対し、スカラはO(１)のメモリ領域で大丈夫。
O(１)のメモリ領域で大丈夫＝→スカラメモリ領域は無視可能
 計算量：O(N/P)
 あまり面白くないｚ98

並列プログラミング講習会（試行）

＋

α

ｘ

ｙ

１ 3 行列とベクトルの積１.3
 ＜行方式＞と＜列方式＞がある。
＜デ タ分散方式＞と＜方式＞組のみ合わせがあり 少し面白い＜データ分散方式＞と＜方式＞組のみ合わせがあり、少し面白い
①②

①②

①②

①

①
②

=

…

…

…

…

②
=

①
②

for(i=0;i<n;i++){y[i]=0.0;
for(j=0;j<n;j++){y[i][i] += a[i][j]*x[j];[i][j] [j]}}for(j=0; j<n; j++) y[j]=0.0;
for(j=0; j<n; j++) {for (i=0; i<n; i++) {y[i][i] += a[i][j]*x[j];[i][j] [j]}}＜行方式＞： 自然な実装
C言語向き
99

＜列方式＞： Fortran言語向き
並列プログラミング講習会（試行）１.3 行列とベクトルの積＜行方式の場合＞＜行方向分散方式＞ ：行方式に向く分散方式
PE=0

PE=0

PE=1

PE=1

=

PE=2
PE=3

PE=2
PE=3

右辺 クトルを MPI_Allgather関数
右辺ベクトルを
MPI Allgather関数を利用し、全PEで所有する各PE内で行列ベクトル積を行う
各PE内で行列
クトル積を行う＜列方向分散方式＞ ：ベクトルの要素すべてがほしいときに向く=
各PE内で行列 クトル積
各PE内で行列-ベクトル積
を行う
100

+

+

+

MPI_Reduce関数で総和を求める
MPI
Reduce関数で総和を求める（※ある１PEにベクトルすべてが集まる）
並列プログラミング講習会（試行）１.3 行列とベクトルの積＜列方式の場合＞＜行方向分散方式＞ ：無駄が多く使われない
PE=0

PE=0

PE=1

PE=1

PE=2

= + + + = PE=2

PE=3

PE=3

右辺 クトルを MPI_Allgather関数
右辺ベクトルを
MPI Allgather関数を利用して、全PEで所有する結果をMPI_Reduce関数により結果をMPI
Reduce関数により総和を求める＜列方向分散方式＞ ：列方式に向く分散方式=
各PE内で行列 クトル積
各PE内で行列-ベクトル積
を行う
101

+

+

+

MPI_Reduce関数で総和を求める
MPI
Reduce関数で総和を求める（※ある１PEにベクトルすべてが集まる）
並列プログラミング講習会（試行）

参考文献
1.
2.

3.
4.
5.

102

MPI並列プログラミング、P.パチェコ 著 ／ 秋葉 博 訳
並列プログラミング虎の巻MPI版、青山幸也 著、
理化学研究所情報基盤センタ（ http://accc.riken.jp/HPC/training/text.html ）Message Passing Interface Forum（ http://www.mpi-forum.org/ ）MPI-Jメーリングリスト（http://phase.hpcc.jp/phase/mpi-j/ml/）
並列コンピュータ工学、富田眞治著、昭晃堂（１９９６）

並列プログラミング講習会（試行）プログラム実習Ⅰ(BLAS)（演習）
東京大学情報基盤センター 特任准教授 片桐孝洋

２０１０年９月２７日（月）１５：１５－１７：００
103

並列プログラミング講習会（試行）演習内容
 BLASとはとは GOTO BLASとは
 LAPACKとは ScaLAPACKとは BLASの利用法と実習（DGEMM）

104

並列プログラミング講習会（試行）１ 8 BLASとPBLAS
１.8
 BLAS（Basic Linear Algebra Subprograms、

基本線形代数副プログラム集）
 線形代数計算で用いられる、基本演算を標準化
線形代数計算で用いられる 基本演算を標準化（API化）したもの。
 普通は、密行列用の線形代数計算用の基本演算
普通は 密行列用の線形代数計算用の基本演算の副プログラムを指す。
 疎行列の基本演算用の＜スパースBLAS＞というも
疎行列 基本演算用
パと うものあるが、まだ定着していない。

105

並列プログラミング講習会（試行）１ 8 BLASとPBLAS
１.8
 高性能なライブラリの＜作成の手間＞と、＜プログラム再利用性＞を高める目的で提案
 BLAS演算の性能改善を、個々のユーザが、個々のプログラムで独立して行うのは、ソフトウエア開発効率が悪い
 ＜工学的＞でない
 性能を改善するチュ
性能を改善するチューニングは
ニングは、経験のないユ
経験のないユーザザには無理 ＜職人芸＞ ！＝ ＜科学、工学＞＜科学 工学＞
 BLASは、数学ソフトウエアにおける＜ソフトウエア工学＞のはしり
106

並列プログラミング講習会（試行）１ 8 BLASとPBLAS
１.8
 BLASでは、以下のように分類わけをして、類

サブルーチンの命名規則を統一
1.
2.
3.
4.

演算対象のベクトルや行列の型（整数型、実数型、複素型）
行列形状（対称行列、三重対角行列）
データ格納形式（帯行列を二次元に圧縮）演算結果が何か（行列、ベクトル） 演算性能から、以下の３つに演算を分類
演算性能
、以
演算を分類
 レベル１ BLAS：

ベクトルとベクトルの演算 レベル２ BLAS： 行列とベクトルの演算 レベル３ BLAS： 行列と行列の演算
107

並列プログラミング講習会（試行）１.8 BLASとPBLAS
 レベル１BLAS ベクトル内積、ベクトル定数倍の加算、など

y ← αx+ y データの読み出し回数、演算回数がほほ同じ
 データの再利用（キャッシュに乗ったデータの再利用によるデータアクセス時間の短縮）がほとんどできない 実装による性能向上が、あまり期待できない ほとんど、計算機ハードウエアの演算性能
 レベル１BLASのみで演算を実装すると、演算が本来持っているデ タ再利用性がなくなる
るデータ再利用性がなくなる
 例：行列-ベクトル積を、レベル１BLASで実装
 例：

108

並列プログラミング講習会（試行）１.8 BLASとPBLAS
BLAS
 行列-ベクトル積などの演算

 レベル２

 例： y ← α A x + β y 前進/後退代入演算、T x = y （Tは三角行列）をxについて解く演算、を含む
 レベル１BLASのみの実装よる、データ再利用性の喪失を回避する目的で提案
 行列とベクトルデータに対して、データの再利用性あり デ
データアクセス時間をタアクセス時間を、実装法により短縮可能
実装法により短縮可能 （実装法により）性能向上がレベル１BLASに比べしやすい（が十分でない）
109

並列プログラミング講習会（試行）

１8
１.8

BLASとPBLAS

BLAS
 行列-行列積などの演算

 レベル３C←αAB+βC
 共有記憶型の並列ベクトル計算機では、レベル２ BLASでも性能向上が達成できない。
 例
例：

 並列化により１PE当たりのデータ量が減少する。
 より大規模な演算をとり扱わないと、再利用の効果がない。

 行列-行列積では、行列データ
行列 行列積 は 行列デ タ

2

対O(n ) に対して演算は O ( n 3 ) なので、データ再利用性が原理的に高い。
 行列積は、アルゴリズムレベルでもブロック化できる。
行列積は アルゴリズムレベルでもブロック化できる
さらにデータの局所性を高めることができる。
110

並列プログラミング講習会（試行）典型的なBLASの性能
性能[FLOPS]

理論性能の限界

BLAS３

BLAS２

キャッシュ
の効果

BLAS1

行列サイズ
111

並列プログラミング講習会（試行）BLAS利用例
 倍精度演算BLAS3
C := alpha*op( A )*op( B ) + beta*C
A: M*K; B:K*N; C:M*N;
CALL DGEMM( ‘N’,‘N’ ‘N’,‘N’ n, n, n, ALPHA,
ALPHA A
A, N
N, B,
B N,
N BETA,
BETA C,C N)
Aが転置しているか
転置
Bが転置しているか
Mの大きさ
Nの大きさ

alpha
の値
Aの
アドレス
Aの１次元目
の要素数

Bの
アドレス

Bの１次元目の要素数

Kの大きさ

引数が多すぎ！ とても使いにくい！
112

beta
の値

並列プログラミング講習会（試行）

Cのアドレス
Cの１次元目
の要素数

１8
１.8

BLASとPBLAS PBLAS
 並列版のBLAS

 BLASとほぼ同じインタフェースをもつ
 BLAS利用者が、容易に移行できる ライブラリ再利用の目的で開発
ライ ラリ再利用 目的 開発（リン・パック） 固有値計算用ライブラリ：EISPACK （アイス・パック）
 これらを統合したライブラリ：LAPACK （エル・エー・パック）
 連立一次方程式用ライブラリ：LINPACK アルゴリズムレベルでブロック化して、レベル３

BLASを利用するアルゴリズムを新規開発
 LAPACKを分散メモリ並列化：ScaLAPACK（スカラ・パック） 並列版BLASとして、PBLASを利用
並列版BLASとして PBLASを利用
113

並列プログラミング講習会（試行）BLASの機能詳細
http://www.netlib.org/blas/pg 命名規則： 関数名：XYYYY 詳細はHP:
 X：

データ型
デ
タ型
S:単精度、D：倍精度、C：複素、Z：倍精度複素 YYYY：YYYY 計算の種類
 レベル１：

例：AXPY：ベクトルをスカラー倍して加算
例
AXPY ベクトルをスカラ 倍して加算 レベル２：
例 GEMV 一般行列とベクトルの積
例：GEMV:
般行列とベクトルの積 レベル３：
例：GEMM:一般行列どうしの積
114

並列プログラミング講習会（試行）インタフェース例：DGEMM（１/４）

 DGEMM(TRANSA,TRANSB,(TRANSA
TRANSB M,
M N,
N K
K, ALPHA,
ALPHA A,A LDA,LDA B,B LDB,LDB
BETA, C, LDC)C :== alpha*op(
l h * ( A )*op()* ( B ) + beta*C
b t *C の計算をする op( X ) = X もしくは op( X ) = X’（Xの転置行列）


 引数


TRANSA（入力） ‐ CHARACTER*１
 TRANSA は op( A ) の操作を指定する。以下の文字列を指定。
の操作を指定する 以下の文字列を指定
 TRANSA = ‘N’ もしくは 'n', op( A ) = A
 TRANSA = ‘T’ もしくは 't', op(p( A ) = A' TRANSA = 'C' or 'c', op( A ) = A‘

TRANSB（入力） ‐ CHARACTER*１
 TRANSB は op( B ) の操作を指定する。以下同様。

115

並列プログラミング講習会（試行）インタフェース例：DGEMM


M（入力） - INTEGER
 op( A )

と 行列 Cの行の大きさを指定する。
行 大きさを指定する

N（入力） - INTEGER
 op(( B )



（２/４）

と 行列 Cの列の大きさを指定する。
Cの列の大きさを指定する

K（入力） - INTEGER
 op( A ) の列の大きさ、およびの列の大きさ および op( B ) の行の大きさを指定する。



ALPHA（入力） - DOUBLE PRECISION スカラ値

ALPHAの値を設定する。

A（入力） - DOUBLE PRECISION 行列Aの配列。大きさは ( LDA, ka )で、ka はTRANSA = ‘N’ or ‘n’のときはk。 そうでないときは、m。
そう な ときは
 TRANSA = ‘N’ or ‘n’のときは、m×kの配列の要素に行列Aを含まないといけない。そうでないときは、k mの配列に行列Aの転置を入れる。
いけない。そうでないときは、k×mの配列に行列Aの転置を入れる。
116

並列プログラミング講習会（試行）インタフェース例：DGEMM


（３/４）

LDA（入力） - INTEGER
 行列Aの最初の次元数を入れる。
行列
最初 次 数を入れる TRANSA = ‘N’ もしくはも くは ‘n’ なら、LDAならは max( １, m )でなくてはならない。そうでないなら、 LDA はmax( １, k )でなくてはならない。


B（入力） - DOUBLE PRECISION
 行列Bの配列。大きさは ( LDB, kb )で、kb はTRANSA = ‘N’ or ‘n’のときはn。 そうでないときは、k。
そうでないときは k TRANSA = ‘N’ or ‘n’のときは、k×nの配列の要素に行列Bを含まないといけない。そうでないときは、n×kの配列に行列Bの転置を入れる。
けな 。そう な
き 、
配列 行列
転置を入れる。
LDB（入力） - INTEGER
 行列Bの最初の次元数を入れる。 TRANSA = ‘N’ もしくは ‘n’ なら、LDA

は max( １, k )でなくてはならない。そうでないなら、 LDB はmax( １, n )でなくてはならない。

117

並列プログラミング講習会（試行）インタフェース例：DGEMM

BETA （入力） - DOUBLE PRECISION スカラ値
カ 値



（４/４）

BETAの値を設定する。
値を設定するC （入力／出力） - DOUBLE PRECISION
 行列Cの配列。
行列Cの配列 入力時、m × nの配列に行列Cを入れる。配列には、BETAが0でない限
り、行列Cを入れる。この場合は、Cの入力は必要ない。
 出力時、この配列に、 m × n行列の演算結果( alpha*op( A )*op( B ) + beta*C )が上書きされて戻る。


LDC （入力） - INTEGER
 行列Cの最初の次元数を入れる。LDC は118

max(１, m )でなくてはならない。

並列プログラミング講習会（試行）ＢＬＡＳの問題点
 BLASの問題点

BLASやPBLASを用いると、データ再利用性や並列性が低下するかもしれない1.


インタフェースに合わせるため、無駄な処理
インタフェースに合わせるため
無駄な処理（配列への代入等）が必要になる場合も

2
2.


3.

例：レベル１BLASにおける行列-ベクトル積＜メモリ浪費＞や＜演算性能低下＞の要因に

ソースコードが読みにくくなる
 BLASのインタフェースを熟知しないと、かえって処理が理解できない まあ、再利用性・性能とのトレードオフでしょうが
まあ 再利用性 性能とのトレ ドオ でし うが
119

並列プログラミング講習会（試行）GOTO BLASとは 後藤和茂 氏により開発された、ソースコードが無償入手可能な、高性能ＢＬＡＳの実装（ライブラリ）
償
高
ブ
 特徴
マルチコア対応がなされている
 多くのコモディティハードウエア上の実装に特化 Intel Nehalem and Atom systems
 VIA Nanoprocessor
 AMD Shanghai
Sh h i andd Istanbul
It b l等
 テキサス大学先進計算センター（ＴＡＣＣ）で、
テキサス大学先進計算センター（ＴＡＣＣ）で

GOTO BLAS2として、ソースコードを配布している
120

ＨＰ：http://www tacc utexas edu/tacc projects/gotoblas2/
ＨＰ：http://www.tacc.utexas.edu/tacc-projects/gotoblas2/
並列プログラミング講習会（試行）LAPACK 密行列に対する、連立一次方程式の解法、および固有値の解法の“標準”アルゴリズムル チンを“標準”アルゴリズムルーチンを無償で提供
 その道の大学の専門家が集結
カリフォルニア大バークレー校：
James Demmel教授
 テネシー大ノックスビル校：
Jack Dongarra教授
 HP

http://www netlib org/lapack/
http://www.netlib.org/lapack/121

並列プログラミング講習会（試行）LAPACKの命名規則
 命名規則： 関数名：XYYZZZ X： データ型
S:単精度、D：倍精度、C：複素、Z：倍精度複素 YY： 行列の型
BD: 二重対角、DI：対角、GB：一般帯行列、GE：一般行列、
HE:複素エルミート、HP：複素エルミート圧縮形式、SY：対称
複素
ミ ト
複素
ミ ト圧縮形式
対称行列、….
 ZZZ：ZZZ 計算の種類
TRF: 行列の分解、TRS：行列の分解を使う、CON：条件数の計算、RFS：計算解の誤差範囲を計算、TRI：三重対角行列の分解、EQU：スケーリングの計算、…

122

並列プログラミング講習会（試行）インタフェース例：DGESV （１/３）
 DGESV(N, NRHS, A, LDA, IPIVOT, B, LDB, INFO)A X = B の解の行列Xを計算をする
 A * X = B、ここで A はN×N行列で、 X と B は N×NRHS行列とする。
 行交換の部分枢軸選択付きのLU分解
行交換 部分枢軸選択付き
分解 でA を A = P * L * U と分解する。ここで、P は交換行列、L は下三角行列、Uは上三角行列である。
である 分解されたA は、連立一次方程式A * X = Bを解くのに使われる。


 引数


N (入力) - INTEGER
 線形方程式の数。行列Aの次元数。

123

N >= 0。

並列プログラミング講習会（試行）インタフェース例：DGESV （２/３）
NRHS (入力) – INTEGER 右辺ベクトルの数。行列Bの次元数。 NRHS >= 0。

A (入力／出力) – DOUBLE PRECISION, DIMENSION(:,:)
 入力時は、N×Nの行列Aの係数を入れる。
入力時は N×Nの行列Aの係数を入れる= P*L*Uを圧縮して出力する。
Lの対角要素は１であるので、収納されていない。

 出力時は、Aから分解された行列LとU



LDA (入力) – INTEGER 配列Aの最初の次元の大きさ。 LDA>= max(1,N)。

IPIVOT (出力) － DOUBLE PRECISION, DIMENSION(:)
 交換行列Aを構成する枢軸のインデックス。 行列のi行がIPIVOT(i)行と交換されている。
換されている

124

並列プログラミング講習会（試行）インタフェース例：DGESV （３/３）
B (入力／出力) – DOUBLE PRECISION, DIMENSION(:,:)
N×NRHS 行列Bを入れる。
 出力時は、もし、INFO = 0 なら、N×NRHS行列である解行列Xが戻る。
 入力時は、右辺ベクトルのLDB (入力) －INTEGERINTEGER 配列Bの最初の次元の大きさ。 LDB >= max(1,N)。



INFO (出力) ーINTEGERINTEGER = 0: 正常終了
 < 0: もし INFO = -i

なら i-th 行の引数の値がおかしい。
 > 0: もし INFO = i なら U(i,i) が厳密に0である。 分解は終わるが、Uの分解は特異なため、解は計算されない。

125

並列プログラミング講習会（試行）ScaLAPACK
 密行列に対する、連立一次方程式の解法、および固有値の解法の“標準”
びアルゴリズムルーチンの並列化版を無償で提供
 ユーザインタフェースはLAPACKに＜類似＞ ソフトウエアの＜階層化＞がされている
内部ルーチンはLAPACKを利用 並列インタフェースはBLACS
 データ分散方式に、２次元ブロック・サイクリック分散方式を採用
 HP: http://www.netlib.org/scalapack/126

並列プログラミング講習会（試行）ScaLAPACKソフトウエア構成図（参照）http://www.netlib.org/scalapack/poster.html
127

並列プログラミング講習会（試行）BLACSとPBLAS BLACS
ScaLAPACK中で使われる通信機能を関数化したもの。
 通信ライブラリは、MPI、PVM、各社が提供する通信ライブラリを想定し、ScaLAPACK内でコード修正せずに使うことを目的とする いわゆる、通信ライブラリのラッパー的役割でScaLAPACK内で利用

現在、MPIがデファクトになったため、MPIで構築された
現在
MPIがデフ クトにな たため MPIで構築された
BLACSのみ、現実的に利用されている。
 なので、ScaLAPACKはMPIでコンパイルし、起動して利用するなので ScaLAPACKはMPIでコンパイルし 起動して利用する

 PBLASBLACSを用いてBLASと同等な機能を提供する関数群
 並列版BLASといってよい。


128

並列プログラミング講習会（試行）ScaLAPACKの命名規則
 原則：

LAPACKの関数名の頭に“P”を付けたもの
 そのほか、BLACS、PBLAS、データ分散を
そのほか BLACS PBLAS データ分散を制御するためのScaLAPACK用関数がある。

129

並列プログラミング講習会（試行）インタフェース例：PDGESV （１/４）
 PDGESV( N, NRHS, A, IA, JA, DESCA, IPIV, B, IB, JB, DESCB,INFO )
sub(A) X = sub(B) の解の行列Xを計算をする
解 行列 を計算をする
 ここで sub(A) はN×N行列を分散したA(IA:IA+N-１, JA:JA+N-１)の行列
 X と B は N×NRHS行列を分散したB(IB:IB+N-１, JB:JB+NRHS-１)の行列
 行交換の部分枢軸選択付きのLU分解 でsub(A) をsub(A) = P * L * U と分解する。ここで、P は交換行列、
L は下三角行列、Uは上三角行列である。
 分解されたsub(A) は、連立一次方程式sub(A) * X = sub(B)を解くのに使われる。
解くのに使われる

130

並列プログラミング講習会（試行）インタフェース例：PDGESV （２/４）
N (大域入力) – INTEGER
 線形方程式の数。行列Aの次元数。
線形方程式の数 行列Aの次元数

N >= 0。
0

NRHS (大域入力) – INTEGER 右辺ベクトルの数。行列Bの次元数。 NRHS >
>= 0。



A (局所入力／出力) – DOUBLE PRECISION, DIMENSION(:,:)
 入力時は、N×Nの行列Aの局所化された係数を配列A(LLD_A, LOCc( JA+N-１))を入れる。
 出力時は、Aから分解された行列LとU = P*L*Uを圧縮して出力する。
Lの対角要素は１であるので 収納されていない
Lの対角要素は１であるので、収納されていない。

IA(大域入力) －INTEGER ：sub(A)の最初の行のインデックス
 JA(大域入力) －INTEGERINTEGER ：sub(A)の最初の列のインデックス
 DESCA (大域かつ局所入力) – INTEGER 分散された配列Aの記述子。
分散さ
列
記述 。
131

並列プログラミング講習会（試行）インタフェース例：PDGESV （３/４）


IPIVOT (局所出力) － DOUBLE PRECISION, DIMENSION(:)
 交換行列Aを構成する枢軸のインデックス。
交換行列Aを構成する枢軸のインデックス 行列のi行がIPIVOT(i)行と交換されている。分散された配列( LOCr(M_A)+MB_A )として戻る。


B ((局所入力／出力)) – DOUBLE PRECISION,, DIMENSION(:,:)(,)
N×NRHSの行列Bの分散されたものを(LLD_B, LOCc(JB+NRHS-１))に入れる。
 出力時は、もし、INFO
出力時は もし INFO = 0 なら、N×NRHS行列である解行列Xが、なら N NRHS行列である解行列Xが行列Bと同様の分散された状態で戻る。
 入力時は、右辺ベクトルのIB(大域入力) －INTEGERINTEGER sub(B)の最初の行のインデックス



JJB(大域入力)() －INTEGER sub(B)の最初の列のインデックス



DESCB (大域かつ局所入力) – INTEGER
 分散された配列Bの記述子。

132

並列プログラミング講習会（試行）インタフェース例：PDGESV （３/４）


INFO (大域出力) ーINTEGER
 = 0: 正常終了
 < 0:
 もし i番目の要素が配列で、i番目の要素が配列で そのj要素の値がおかしいなら、
そのj要素の値がおかしいなら

INFO = -(i*100+j)となる。
 もしi番目の要素がスカラーで、かつ、その値がおかしいなら、
INFO = -iとなる。
となる
 > 0: もし INFO = Kのとき U(IA+K-１, JA+K-１) が厳密に0である。
分解は完了するが、分解されたＵは厳密に特異なので、解は計算できない。

133

並列プログラミング講習会（試行）サンプルプログラムの実行（BLAS DGEMM）

134

並列プログラミング講習会（試行）ＵＮＩＸ防忘録
 emacsの起動：

emacs 編集ファイル名 ^x^s （＾はcontrol） ：テキストの保存 ^x ^c ： 終了
（ ^z で終了すると、スパコンの負荷が上がる。絶対にしないこと。）
 ^g^: 訳がわからなくなったとき。
訳がわからなくな たとき ^k : カーソルより行末まで消す。消した行は一時記憶される。
 ^y : ^kで消した行を、現在のカーソルの場所にコピーする。
 ^s 文字列 : 文字列の箇所まで移動する。
文字列の箇所まで移動する
 ^x goto-line : 指定した行まで移動する。
135

並列プログラミング講習会（試行）ＵＮＩＸ防忘録
 rm ファイル名：


ファイル名のファイルを消す。

rm *~ : test.c~ などの、~がついたバックアップファイルを消す。

 ls : 現在いるフォルダの中身を見る。
 cd

フォルダ名： フォルダに移動する。

 c
cd.. : 一つ上のフォルダに移動。
のフォルダに移動。
 cd ~ ：ルートディレクトリに行く。訳がわからなくなったとき。
 cat

ファイル名： ファイル名の中身を見る make : 実行ファイルを作る （Makefile があるところでしか実行できない）
 make136

clean : 実行ファイルを消す。
並列プログラミング講習会（試行）BLAS DGEMMサンプルプログラムの注意点
 C言語版、Fortran言語版のファイル名（共通）lecBLAS.tar
 ジョブスクリプトファイルmat-mat-blas.bash
ジョブスクリプトファイルmat mat blas bash

中のキュー名をllectureから tutorial
i l に変更してから
qsub してください。
 lecuter : 実習時間外のキュー（同時実行数１）
 tutorial : 実習時間内のキ
実習時間内のキュー（同時実行数4+）（同時実行数4 ）137

並列プログラミング講習会（試行）BLAS DGEMMのサンプルプログラムの実行（Ｃ言語版）
 以下のコマンドを実行する$ cp /home/z30082/lecBLAS.tar ./$ tar xvf lecBLAS.tar$ cd Mat-Mat-BLAS$ cd C$ make$ 実行が終了したら、以下を実行する
実行が終了したら 以下を実行する$ cat mat-mat-blas.out
138

並列プログラミング講習会（試行）BLAS DGEMMのサンプルプログラムの実行（Ｃ言語版）
 以下のような結果が見えれば成功

N = 1000
Mat-Mat time = 0.885865 [sec.]
2257 680346 [MFLOPS]
2257.680346
OK!

139

並列プログラミング講習会（試行）BLAS DGEMMのサンプルプログラムの実行（Fortran言語版）
 以下のコマンドを実行する$ cp /home/z30082/lecBLAS.tar ./$ tar xvf lecBLAS.tar$ cd Mat-Mat-BLAS$ cd F$ make$ qsub mat-mat-blas.bash
mat mat blas bash
 実行が終了したら、以下を実行する
実行が終了したら 以下を実行する$ cat mat-mat-blas.out
140

並列プログラミング講習会（試行）BLAS DGEMMのサンプルプログラムの実行（Fortran言語版）
 以下のような結果が見えれば成功

N =
1000
Mat-Mat time[sec.] = 0.853669448999999969
MFLOPS = 2342
2342.82718831429702
82718831429702
OK!

141

並列プログラミング講習会（試行）サンプルプログラムの説明（Ｃ言語） #define

N １０００の、数字を変更すると、行列サイズが変更できます #define DEBUG １「１」にすると、行列-行列積の演算結果が検が証証できます。
す。
 MyMatMat関数の仕様
 Ｄｏｕｂｌｅ型N×N行列ＡとＢの行列積をおこない、
Ｄ ｂｌ 型N×N行列ＡとＢの行列積をおこない

Ｄｏｕｂｌｅ型Ｎ×Ｎ行列Ｃにその結果が入ります142

並列プログラミング講習会（試行）Fortran言語のサンプルプログラムの注意
 行列サイズNNの宣言は、以下のファイルにあります。

mat-mat-blas inc
mat-mat-blas.inc 行列サイズ変数が、ＮＮとなっています。
行
ズ変数が
な
ます

integerg NNparameter (NN=１０００)

143

並列プログラミング講習会（試行）GOTO BLAS呼び出しオプション
 日立コンパイラから、GOTO BLASを呼び出す

場合、以下のオプションを付けます。
C言語
cc <プログラム名> -L/opt/itc/lib
L/opt/itc/lib -lgoto
lgoto
Fortran言語
f90 <プログラム名> -L/opt/itc/lib -lgoto※任意数のスレッド並列版（ノード内）を利用する場合は -lgoto_smp
１６スレッド限定版（ノード内）を利用する場合は -lgoto_smp16
を lgotoの代わりに付けてください
を、-lgotoの代わりに付けてください
144

並列プログラミング講習会（試行）サンプルプログラムの説（BLAS）
 １コア（逐次）実行版です スレッド並列化版は、対応するスレッド並列化版をリンクしてコンパイルの上、並列実行数の指定をする必要があります。
 GOTO BLASでは、OpenMPのスレッド数指定法と同じBLASでは OpenMPのスレッド数指定法と同じ方法で、実行するスレッド数が指定できます OMP_NUM_THREADS環境変数に実行スレッド数を入れる
環境変数 実行ッ 数を入れる

145

並列プログラミング講習会（試行）演習課題１（BLAS）
1.

MyMatMat関数（手続き）の行列積コードを、
y
BLAS DEGMMルーチンの呼び出しにより高速化してください 引数の並びに注意してください。
 C言語は、DGEMM
言語はのFortran手続き
手続きを呼び出します。以下に注意してください。
 関数名： DGEMM_
 全ての引数が ポインタ引き渡し

146

並列プログラミング講習会（試行）演習課題２（BLAS）
2.

147

スレッド並列版のDGEMMを呼び出し、いろいろなスレッド数で実行してください。
また逐次実行に比べ どれだけ高速化されまた逐次実行に比べ、どれだけ高速化されたか、性能評価をしてください。

並列プログラミング講習会（試行）BLAS DGEMM回答
回答

148

並列プログラミング講習会（試行）BLAS DGEMMの回答（C言語版）

double ALPHA,, BETA;;
char TEX[1] = {'N'};
ALPHA=1.0;
ALPHA=1
0;
BETA=1.0;
dgemm_(&TEX, &TEX, &n, &n, &n, &ALPHA,
A &n,
A,
& B,
B &n,
& &BETA,&BETA C,C &n);& )149

並列プログラミング講習会（試行）BLAS DGEMMの回答（Fotran言語版）

double pprecision ALPHA,BETA,
ALPHA=1.0d0
BETA=1 0d0
BETA=1.0d0
CALL DGEMM('N',('N', n, n, n, ALPHA,&
A, n, B, n, BETA, C, n)

150

並列プログラミング講習会（試行）参考文献（１）

BLAS
http://www.netlib.org/blas/
2. LAPACK
http://www.netlib.org/lapack/
ttp://www. et b.o g/ apac /3. ScaLAPACK
http://www.netlib.org/scalapack/
//
/
/
4 スパ
4.
スパースBLAS
スBLAS
http://math.nist.gov/spblas/
1.

151

並列プログラミング講習会（試行）参考文献（２）
1.
2.

3.
4.
5.
6.

152

MPI並列プログラミング、P.パチェコ 著 ／ 秋葉 博 訳
並列プログラミング虎の巻MPI版、青山幸也 著、
理化学研究所情報基盤センタ（ http://accc.riken.jp/HPC/training/text.html ）Message Passing Interface Forum（ http://www.mpi-forum.org/ ）MPI-Jメーリングリスト（http://phase.hpcc.jp/phase/mpi-j/ml/）
並列コンピュータ工学、富田眞治著、昭晃堂（１９９６）
並列数値処理 ー高速化と性能向上のためにー、
金田康正 編著、コロナ社（２０１０）
並列プログラミング講習会（試行）プログラム実習Ⅱ(LAPACK, ScaLAPACK)（演習）
東京大学情報基盤センター 特任准教授 片桐孝洋

２０１０年９月２８日（火）１０：３０－１２：３０
153

並列プログラミング講習会（試行）演習内容
 粒子間熱伝導問題の説明（座学）
 LAPACKの利用法と演習（C言語、Fortran言語）（C言語
Fortran言語） ScaLAPACKの利用法と演習（Fortran言語のみ）

154

並列プログラミング講習会（試行）粒子間熱伝導問題
 東京大学情報基盤センター

中島研吾教授から提供頂いた、PPTとサンプルプログラムを利用 詳細資料 http://nkl.cc.u-tokyo.ac.jp/09e/CE23/CE23.pdf

155

並列プログラミング講習会（試行）問題設定（1/5）/
 空間上に規則正しく，X,Y,Z方向にNX,NY,NZ個ずつ等間隔（DX,DY,DZ）に配置された粒子の集合体を考える。

156

並列プログラミング講習会（試行）問題設定（2/5）
粒子iに関する熱流束の/釣り合い
DEL 
RADI
RADI maxj

AREAT j  Ti   HCONVi  SURF  T0  Ti   QVOL V  0
DELij / CONDij

i

Y

X
157

並列プログラミング講習会（試行）問題設定（3/5）
粒子間熱伝導DEL 
RADIRADI max

i

AREA
T j  Ti   HCONVi  SURF  T0  Ti   QVOL V  0
DELij / CONDij
粒子iを中心とする，半径RADImaxの球（円）

i
Y

X
158

並列プログラミング講習会（試行）問題設定（3/5）
粒子間熱伝導DEL 
RADIRADI max

i

AREA
T j  Ti   HCONVi  SURF  T0  Ti   QVOL V  0
DELij / CONDij
粒子iを中心とする，半径RADImaxの球（円）

i
Y

X
159

並列プログラミング講習会（試行）

粒子間の距離DELijがRADImaxのよりも小さいも さ
粒子群とのみ熱伝導がある。

問題設定（3/5）
粒子間熱伝導DEL  RADI max


i

AREA
T j  Ti   HCONVi  SURF  T0  Ti   QVOL V  0
DELij / CONDij
断面積＝AREA
断面積
AREA

粒子iを中心とする，半径RADImaxの球（円）

熱伝導率
CONDiji

粒子間の距離DELijがRADImaxのよりも小さい
粒 群と粒子群とのみ熱伝導が熱伝導がある。
長さDEL，断面積AREA，
熱伝導率CONDijの管で連結

160
並列プログラミング講習会（試行）熱伝導率：CONDij
DEL 
RADI
RADI max

i

AREA
T j  Ti   HCONVi  SURF  T0  Ti   QVOL V  0
DELij / CONDijｊ

COND 0
CONDij 
DEL20
min 10 ,10





DELij
i

COND0 「基準」熱伝導率

161

並列プログラミング講習会（試行）問題設定（4/5）
対流熱伝達DEL  RADI max


i

AREA
T j  Ti   HCONVi  SURF  T0  Ti   QVOL V  0
DELij / CONDij
X=0
対流熱伝達率HCONVi= HCONV if X=0=0
otherwise
i

Y

雰囲気温度
T0
X

162

熱交換面積
SURF

並列プログラミング講習会（試行）問題設定（5/5）
体積発熱
DEL 
RADI
RADI max

i

AREA
T j  Ti   HCONVi  SURF  T0  Ti   QVOL V  0
DELij / CONDij

体積発熱量

QVOL

発熱体積

V

V  C1  VOL  C2  VOL  DELQ
DELQ  X i  Yi  Z i
2

2

2

C1, C2 定数VOL
各粒子の「基準」体積
DELQ 粒子中心と原点との距離163

並列プログラミング講習会（試行）粒子 に関する りあ
粒子iに関するつりあい
DEL  RADI max


j

AREA
T j  Ti   HCONVi  SURF  T0  Ti   QQVOL V  0
DELij / CONDij  DEL RADI max

AREA
AREA
j  DEL / COND T j    j DEL / COND Ti

 
ij
ij 
ij
ij

 HCONVi  SURF  T0  HCONVi  SURF  Ti  QVOL  V  0
DEL  RADI max

DEL  RADI max  DEL RADI max


AREA
AREA
 HCONVi  SURF Ti   
Tj 
 
DELij / CONDij
j
j


 DELij / CONDij 

AMAT(i,i)（対角成分）( )AMAT(i,j)（非対角成分）( ,j)（非対角成分）

  HCONVi  SURF  T0  QVOL  V
RHS（右辺）164

並列プログラミング講習会（試行）係数行列 AMAT
 影響範囲 RADImax が小さければ，係数行列は疎

局所的な効果 影響範囲が大きければ，係数行列は密

非ゼ 成分 割合が大きくなる
非ゼロ成分の割合が大きくなるi

165

並列プログラミング講習会（試行）疎行列と密行列： {q}=[A]{p}{q}{p}疎行列：差分法，有限要素法，有限体積法
do i= 1, N
q(i)= D(i)*p(i)
d k
do
k= i
index(i-1)+1,
d (i 1) 1 i
index(i)d (i)q(i)= q(i) + AMATs(k)*p(item(k))
enddo
enddo

密行列 スペクトル法 分子動力学 境界要素法
密行列：スペクトル法，分子動力学，境界要素法

166

do i= 1, N
q(i)= 0.d0
0 d0
do j= 1, N
q(i)= q(i) + AMATd(i,j)*p(j)
enddo
dd
enddo
並列プログラミング講習会（試行）係数行列 AMAT
 影響範囲 RADImax が小さければ，係数行列は疎

局所的な効果 影響範囲が大きければ，係数行列は密影響範囲 大きければ，係数行列は密

非ゼロ成分の割合が大きくなる

 影響範囲が大きくなる場合のことを考えて，プログラムの中数( ,j) 」
では係数行列を「密（AMATd(i,j)）」として扱う。

167

並列プログラミング講習会（試行） 問題設定
 1CPU用プログラムの解説

168

並列プログラミング講習会（試行）sps: ファイル，データ

inp.1
inp
1
（名称固定）
制御データ

169

sphere.fld(名称固定)

計算結果（形状）

spheredata(名称固定)

計算結果
（温度）

sps

並列プログラミング講習会（試行）制御データ：inp1
10 10 101.e0 1.e0 1.e0
1.e0 10.e0 1.e0 1.e0
1.e24 0.e0 10.e0
8.0e0
1.0e0 0.1e0

170

NX, NY, NZ
DX, DY, DZ
VOL, AREA, QVOL, COND0HCONV,, T0,, SURF
RADImax
C1, C2並列プログラミング講習会（試行）計算結果（形状）：sphere.fld*.fldの形式であれば名称は任意
fldの形式であれば名称は任意# AVS field file
必須
ndim=
ndim
3
三次元モデルであることを示す
次元 デルである とを示す
dim1=
10
NX
dim2=
10
NY
dim3=
10
NZ
nspace=3veclen=1data= float
field= uniform
label= temperature
variable 1 file=./spheredata/
filetype=ascii 温度ファイル名
温度フ イル名

171

並列プログラミング講習会（試行）計算結果（温度）：spheredata
計算結果（形状）からの参照が正しければ名称は任意

1.188402E-24
5.388751E+00
9.485117E+00
1.311638E+01
1 630137E+01
1.630137E+01…

open (22, file='spheredata', status='unknown')
do i= 1, NX*XY*NZ
write (22,'(1pe16.6)') PHI(i)
enddo

172

並列プログラミング講習会（試行）シリアル版の計算手順
 制御データ入力
 粒子生成粒 生成
 係数マトリクス計算

熱伝導 対流熱伝達 発熱
 CG法による連立一次方程式求解（オリジナル版
ジ版 test_org.f, test_org.c）
密行列
 点ヤコビ前処理
ビ


 出力

173

MicroAVS用
並列プログラミング講習会（試行）test.f：シリアル版（1/10）
初期設定
implicit REAL*8 (A-H,O-Z)real(kind=8) :: VOL, AREA, QVOL, COND0, COND, RADImaxreal(kind=8) :: HCONV, T0, SURF, DEL, DEL0, coef1, coef2
real(kind=8), dimension(:) , allocatable :: XC, YC, ZC
real(kind=8), dimension(:,:), allocatable :: AMAT
real(kind=8), dimension(:) , allocatable :: RHS
real(kind=8) dimension(:) , allocatable :: PHIreal(kind=8),real(kind=8), dimension(:,:), allocatable :: W
integer :: NX, NY, NZ, N
integerg:: R, Z, P, Q, DD!C!C +------+!C | INIT |!C +------+!C===!C
open (11, file= 'inp1', status='unknown')
read (11,*) NX, NY, NZ
read (11,*) DX, DY, DZ
read ((11,*), ) VOL,, AREA,, Q
QVOL,, COND0
read (11,*) HCONV, T0, SURF
read (11,*) RADImax
read (11,*) coef1, coef2
close (11)
N= NX*NY*NZ
174

並列プログラミング講習会（試行）test.f：シリアル版（2/10）
粒子中心の座標（XC YC ZC）
粒子中心の座標（XC,YC,ZC）
allocate (XC(N)(XC(N), YC(N)YC(N), ZC(N))
icou= 0do k= 1, NZ
do j= 1, NY
do i= 1, NX
icou= icou + 1
XC(icou)= dfloat(i-1)*DX
YC(icou)= dfloat(j-1)*DY
ZC(icou)= dfloat(k-1)*DZ
enddo
enddo
enddo
!C===

粒子の通し番号 ii （1~N）
ii=(k-1)*NX*NY + (j-1)*NX + i

175

並列プログラミング講習会（試行）test.f：シリアル版（3/10）
マトリクス生成：熱伝導部分!C!C +--------+++!C | MATRIX |!C +--------+!C===
allocate (AMAT(N,N), RHS(N))

 DEL RADI max
AREA HCONVi  SURF Ti 
 
/
DEL
COND
j


ij
ij
DEL  RADI max AREAj  DEL / COND T j 

ij
ij

  HCONVi  SURF  T0  QVOL  V

AMAT= 0.d0
RHS = 0.d0
do i= 1, N
do j= 1
1, N
if (j.ne.i) then
DEL= dsqrt((XC(i)-XC(j))**2 + (YC(i)-YC(j))**2&&+ (ZC(i)-ZC(j))**2)
if (DEL.le.RADImax) then
COND= COND0/(10.d0**dmin1(DEL,20.d0))
coef= COND*AREA / DEL
AMAT(i,j)= coef
COND 0
AMAT(i,i)= AMAT(i,i) - coef
CONDijj 
endif
min 10 DEL ,10 20
endif
enddo
enddo



176

並列プログラミング講習会（試行）test.f：シリアル版（3/10）
マトリクス生成：熱伝導部分!C!C +--------+++!C | MATRIX |!C +--------+!C===
allocate (AMAT(N,N), RHS(N))

 DEL RADI max
AREA HCONVi  SURF Ti 
 
/
DEL
COND
j


ij
ij
DEL  RADI max AREAj  DEL / COND T j 

ij
ij

  HCONVi  SURF  T0  QVOL  V

AMAT= 0.d0
RHS = 0.d0
do i= 1, N
do j= 1
1, N
if (j.ne.i) then
DEL= dsqrt((XC(i)-XC(j))**2 + (YC(i)-YC(j))**2&+ (ZC(i)-ZC(j))**2)
if (DEL.le.RADImax) then
COND= COND0/(10.d0**dmin1(DEL,20.d0))
coef= COND*AREA / DEL
AMAT(i,j)= coef
AMAT(i,i)= AMAT(i,i) - coef 対角項
endif
endif
enddo
enddo

177

並列プログラミング講習会（試行）&test.f：シリアル版（4/10）
マトリクス生成：体積発熱
do i= 1, N
DELQ= dsqrt(XC(i)**2 + YC(i)**2 + ZC(i)**2)
RHS(i)= -QVOL*(coef1*VOL + coef2*DELQ*VOL)
enddo
i= 1
do k= 1, NZ
do j= 1, NY
ic= (k-1)*NX*NY + (j-1)*NX + i
AMAT(ic,ic)= -HCONV*SURF+ AMAT(ic,ic)
RHS (ic)= -HCONV*SURF*T0 + RHS (ic)
enddo
enddo
!C===

 DEL RADI max

AREA
HCONV
SURF




Ti 
j DEL / COND
i


ij
ij
DEL  RADI max AREAj  DEL / COND T j 


ij
ij
  HCONVi  SURF  T0  QVOL  V
178

V  C1 VOL  C2  VOL  DELQ
DELQ  X i  Yi  Z i

並列プログラミング講習会（試行）

2

2

2

test.f：シリアル版（4/10）
マトリクス生成：対流熱伝達
do i= 1, N
DELQ= dsqrt(XC(i)**2 + YC(i)**2 + ZC(i)**2)
RHS(i)= -QVOL*(coef1*VOL + coef2*DELQ*VOL)
enddo
i= 1
do k= 1, NZ
do j= 1, NY
ic= (k-1)*NX*NY + (j-1)*NX + i
AMAT(ic,ic)= -HCONV*SURF+ AMAT(ic,ic)
RHS (ic)= -HCONV*SURF*T0 + RHS (ic)
enddo
enddo
!C===

 DEL RADI max
AREA HCONVi  SURF Ti 
 
DEL
/
COND
j

ij
ij

DEL  RADI max AREATj  DEL / COND j 


ij
ij
  HCONVi  SURF  T0  QVOL  V
179

並列プログラミング講習会（試行）前処理付き共役勾配法
Preconditioned Conjugate Gradient Method （CG）
Compute r(0)= b-[A]x(0)for i= 1, 2, …
solve [M]z(i-1)= r(i-1)
i-1= r(i-1) z(i-1)
if i=1p(1)= z(0)elsei-1= i-1/i-2p(i)= z(i-1) + i-1 p(i-1)
endifq(i)= [A]p
[A] (i)i = i-1/p(i)q(i)x(i)= x(i-1) + ip(i)(i-1)1) -  q(i)r(i)= r(iicheck convergence |r|
end

180

前処理：対角スケーリング
x(i) ：ベクトル
i ：スカラー

並列プログラミング講習会（試行）test_org.f：シリアル版（5/10）
CG法①!C!C +---------------+!C | CG iterations |!C +---------------+!C===
EPS= 1.d-08
allocate (W(N,4), PHI(N))
W = 0
0.d0
d0
PHI= 0.d0
R = 1
Z = 2
Q = 2
P = 3
DD= 4
do i= 1, N
W(i DD) 1.0D0
W(i,DD)=
1 0D0 / AMAT(i,i)
AMAT(i i)
enddo

W(i,1)= W(i,R) ⇒ {r}W(i,2)= W(i,Z) ⇒ {z}W(i,2)= W(i,Q) ⇒ {q}W(i,3)= W(i,P) ⇒ {p}W(i,4)= W(i,DD) ⇒ 1/DIAG
181

Compute r(0)= b-[A]x(0)for i= 1, 2, …
solve [M]z(i-1)= r(i-1)
i-1= r(i-1) z(i-1)
if i=1p(1)= z(0)elsei-1= i-1/i-2p(i)= z(i-1) + i-1 p(i-1)
endifq(i)= [A]p
[A] (i)i = i-1/p(i)q(i)x(i)= x(i-1) + ip(i)(i-1)1) -  q(i)r(i)= r(iicheck convergence |r|end

並列プログラミング講習会（試行）test_org.f：シリアル版（6/10）
CG法②!C!C-- {r0}= {b} - [A]{xini} |!C
do i= 1, N
W(i,R) = RHS(i)do j= 1, N
W(i,R) = W(i,R) - AMAT(i,j)*PHI(j)
enddo
enddo
BNRM2= 0
0.0D0
0D0
do i= 1, N
BNRM2 = BNRM2
enddo+ RHS(i)**2!C******************************************
do iter= 1, N
!C!C-- {z}= [Minv]{r}do i= 1, N
W(i,Z)= W(i,DD) * W(i,R)
enddo

182

Compute r(0)= b-[A]x(0)for i= 1, 2, …
solve [M]z(i-1)= r(i-1)
i-1= r(i-1) z(i-1)
if i=1p(1)= z(0)elsei-1= i-1/i-2p(i)= z(i-1) + i-1 p(i-1)
endifq(i)= [A]p
[A] (i)i = i-1/p(i)q(i)x(i)= x(i-1) + ip(i)(i-1)1) -  q(i)r(i)= r(iicheck convergence |r|end

並列プログラミング講習会（試行）test_org.f：シリアル版（7/10）
CG法③!C!C-- RHO= {r}{z}!C
RHO= 0.d0
do i= 1, N
RHO= RHO + W(i,R)*W(i,Z)
enddo!C!C-- {p} = {z} if!C
BETA= RHO / RHO1

ITER=1otherwise

if ( iter.eq.1 ) thendo i= 1, N
W(i,P)= W(i,Z)
enddo
else
BETA= RHO / RHO1
do i= 1, N
W(i P) W(iW(i,P)=
W(i,Z)
Z) + BETA*W(i,P)
BETA*W(i P)
enddo
endif

183

Compute r(0)= b-[A]x(0)for i= 1, 2, …
solve [M]z(i-1)= r(i-1)
i-1= r(i-1) z(i-1)
if i=1p(1)= z(0)elsei-1= i-1// i-2
p(i)= z(i-1) + i-1 p(i-1)
endifq(i)= [A]p
[A] (i)i = i-1/p(i)q(i)x(i)= x(i-1) + ip(i)(i-1)1) -  q(i)r(i)= r(iicheck convergence |r|end

並列プログラミング講習会（試行）test_org.f：シリアル版（8/10）
CG法④!C!C-- {q}= [A]{p}do i= 1, N
W(i,Q) = 0.d0
do j= 1, N
W(i,Q) = W(i,Q) + AMAT(i,j)*W(j,P)
enddo
enddo
!C
!C-- ALPHA= RHO / {p}{q}p qC1= 0.d0
do i= 1, N
C1= C1 + W(i,P)*W(i,Q)
enddo
ALPHA= RHO / C1
!C!C-- {{x}=} {{x}} + ALPHA*{p}{p}!C{r}= {r} - ALPHA*{q}do i= 1, N
PHI(i) = PHI(i) + ALPHA * W(i,P)
W (i(i,R)=R) W(i
W(i,R)R) - ALPHA * W(i,Q)
W(i Q)
enddo
184

Compute r(0)= b-[A]x(0)for i= 1, 2, …
solve [M]z(i-1)= r(i-1)
i-1= r(i-1) z(i-1)
if i=1p(1)= z(0)elsei-1= i-1/i-2p(i)= z(i-1) + i-1 p(i-1)
endifq(i)= [A]p
[A] (i)i = i-1/p(i)q(i)x(i)= x(i-1) + ip(i)(i-1)1) -  q(i)r(i)= r(iicheck convergence |r|end

並列プログラミング講習会（試行）test_org.f：シリアル版（9/10）
CG法⑤
DNRM2 = 0.0
do i= 1
1, N
DNRM2= DNRM2 + W(i,R)**2
enddo
RESID= dsqrt(DNRM2/BNRM2)
write (*, 1000) iter, RESID
format (i5, 1pe16.6)

1000

if ( RESIDRESID.le.EPS)
le EPS) goto 900
RHO1 = RHO
enddo!C******************************************
IER = 1900 continue
!C
!C===

185

Compute r(0)= b-[A]x(0)for i= 1, 2, …
solve [M]z(i-1)= r(i-1)
i-1= r(i-1) z(i-1)
if i=1p(1)= z(0)elsei-1= i-1/i-2p(i)= z(i-1) + i-1 p(i-1)
endifq(i)= [A]p
[A] (i)i = i-1/p(i)q(i)x(i)= x(i-1) + ip(i)(i-1)1) -  q(i)r(i)= r(iicheck convergence |r|end

並列プログラミング講習会（試行）test_org.f：シリアル版（10/10）
結果出力!C!C +--------+!C | OUTPUT |!C +--------+
!C===
N1= 1
N3= 3
!C
!C-- MESH
open (22, file='sphere.fld', status='unknown')
write (22,'(a)')'# AVS field file'write (22(22,'(a(a,i5)
i5)')) 'ndim='ndim , N3
write (22,'(a,i5)') 'dim1=', NX
write (22,'(a,i5)') 'dim2=', NY
write (22,'(a,i5)') 'dim3=', NZ
write (22,'(a,i5)') 'nspace=', N3
write (22,'(a,i5)') 'veclen=', N1
write (22,'(a,i5)') 'data= float'
write (22,'(a,i5)') 'field= uniform'
write (22,'(a,i5)') 'label= temperature'
write (22(22,'(a(a,i5)
i5)')) 'variablevariable 1 file=
file=./spheredata/spheredata filetype=ascii'
filetype=ascii
close (22)!C
!C-- RESULTS
open (22, file='spheredata', status='unknown')
do i= 1, N
write (22,'(1pe16.6)') PHI(i)186enddo
並列プログラミング講習会（試行）!C===サンプルプログラムの実行（LAPACK dgesv）

187

並列プログラミング講習会（試行）LAPACK dgesvサンプルプログラムの注意点
 C言語版、Fortran言語版のファイル名（共通）lecLAPACK.tar
 ジョブスクリプトファイルgo.sh
ジョブスクリプトファイルgo sh

中のキュ 名を中のキュー名をlecture から tutorial に変更してから
qsub
b してください。
してください lecuter : 実習時間外のキュー（同時実行数１）
 tutorial : 実習時間内のキュー（同時実行数4+）

188

並列プログラミング講習会（試行）PGFコンパイラ(ACMLライブラリ利用のため)のコマンド
 Ｃ言語、Ｆｏｒｔｒａｎ言語ともに、ＰＧＦコンパイラを利用するため、以下のコマンドを実行してください。

source /opt/itc/mpi/mpiswitch.sh mpich-mx-pgi
なお、サンプルプログラムのファイルに上記
なお、サンプルプログラムのファイルに上記コマンドがあります。以下を入力してください。

$ cat pgf
pgf.shsh189

並列プログラミング講習会（試行）LAPACK dgesvのサンプルプログラムの実行（Ｃ言語版）
 以下のコマンドを実行する$ cp /home/z30082/lecLAPACK.tar ./$ tar xvf lecLAPACK.tar$ cd sphere-LAPACK$ cd C$ make$ qsub go.sh
go sh 実行が終了したら、以下を実行する
実行が終了したら 以下を実行する$ cat test.out
190

並列プログラミング講習会（試行）LAPACK dgesvのサンプルプログラムの実行（Ｃ言語版）
 以下のような結果が見えれば成功
time = 0.000000 [sec.]
inf [MFLOPS]
990 -2.272792e+00
991 -2.276715e+00
992 -2.288410e+00
2.288410e 00
993 -2.307670e+00
994 -2.334166e+00
995 -2.367479e+00
996 -2.407125e+00
997 -2.452584e+00-2 452584e+00

連立一次方程式求解部分が実装されて無いので、
装
時間が0秒、
誤差ERRが大きくなっていますが、正常な動作です。

998 -2.503330e+00
999 -2.558846e+00
err = 6.274235e+02
191

並列プログラミング講習会（試行）LAPACK dgesvのサンプルプログラムの実行（Fortran言語版）
 以下のコマンドを実行する$ cp /home/z30082/lecLAPACK.tar ./$ tar xvf lecLAPACK.tar$ cd sphere-LAPACK$ cd F$ make$ qsub go.sh
go sh 実行が終了したら、以下を実行する
実行が終了したら 以下を実行する$ cat test.out
192

並列プログラミング講習会（試行）LAPACK dgesvのサンプルプログラムの実行（Fortran言語版）
 以下のような結果が見えれば成功
TIME[sec] = 9.9999999999999995E-007
MFLOPS =670168000.0000000

991 -2.272792E+00
992 -2.276715E+00
993 -2.288410E+00
2.288410E 00
994 -2.307670E+00
995 -2.334166E+00
996 -2.367479E+00
997 -2.407125E+00
998 -2.452584E+00-2 452584E+00

連立一次方程式求解部分が実装されて無いので、時間が極端に短く、誤差ERRが大きくなっていますが、正常な動作です。

999 -2.503330E+00
1000 -2.558846E+00
EPS =
193

627.4235113914224
並列プログラミング講習会（試行）サンプルプログラムの説（LAPACK）
 １コア（逐次）実行版です スレッド並列化版は、スレッド並列化版をリンクしてコンパイルの上、並列実行数の指定をする
、
実 数
定 す
必要があります。
 LAPACKはスレッド並列化のみ対応しているため、複数ノードを通信しつつ利用することはできません 通信しつつ利用するには、分散版であるScaLAPACKを通信しつつ利用するには 分散版であるScaLAPACKを

利用します
 ノ
ノード内のみＬＡＰＡＣＫを使いド内のみＬＡＰＡＣＫを使い、ノノード間はＭＰＩで並列化
ド間はＭＰＩで並列化するような使い方は可能です たとえば、領域分割法を利用し、ノード内のみＬＡＰＡＣＫで連立一次方程式を解く場合
194

並列プログラミング講習会（試行）ACML LAPACK呼び出しオプション
 以下のオプションを付けます。

C言語

pgcc <プログラム名>-lacmllacml -lpgftnrtl
lpgftnrtl -pgf90libs
pgf90libs※C言語からFortranのACMLライブラリが提供するLAPACK手続き
を呼び出すため f90ランタイムライブラリの指定が必要ですを呼び出すため、f90ランタイムライブラリの指定が必要ですFortran言語
Fortran言語

pgf95 <プログラム名> -lacml
195

並列プログラミング講習会（試行）演習課題1（LAPACK）
1.

test.c もしくは test.f のメイン関数（手続き）の連立一次方程式求解部分を、LAPACK dgesv
り、高速ルーチンの呼び出しにより、高速化してください 引数の並びに注意してください。
 C言語は、dgesv
C言語は dgesv のFortran手続きの呼び出しになります。以下に注意してください。
 関数名：
関数名ddgesv_dgesv_

 全ての引数がポインタ引き渡し

 デバッグとして、誤差（ERR）が十分小さいことを確認してください。
196

並列プログラミング講習会（試行）演習課題2（LAPACK）
CG法による解法ルーチンtest_org.cgもしくはtest_org.f をコンパイルして、
LAPACK版と実行速度を比べてください。
LAPACK版と実行速度を比べてください
問題サイズを大きくして、LAPACK版、CG法版の実行速度を比べてください。

2.

3.

197

並列プログラミング講習会（試行）LAPACK dgesv回答
g
回答

198

並列プログラミング講習会（試行）LAPCK dgesvの回答（C言語版）

dgesv_(&nn,d(& &inc,&i amat2,
t2 &&nn,
piv rhs,piv,
rhs &nn,&nn &info);

199

並列プログラミング講習会（試行）LAPCK dgesvの回答（Fotran言語版）

call DGESV(N, INC, AMAT, N,&
PIV
PIV, RHS,
RHS N,N INFO)200

並列プログラミング講習会（試行）サンプルプログラムの実行（ScaLAPACK pdgesv）

201

並列プログラミング講習会（試行）ScaLAPACK pdgesvサンプルプログラムの注意点
 Fortran言語版のファイル名

lecLAPACK.tar※ScaLAPACKはC言語版のサンプルはありません ジョブスクリプトファイルgo.sh

中のキュー名をlecture から tutorial に変更してから
qsub してください。
 lecuter : 実習時間外のキュー（同時実行数１）
 tutorial : 実習時間内のキュー（同時実行数4+）202

並列プログラミング講習会（試行）ScaLAPACKにおけるプロセッサ
ScaLAPACKにおけるプロセッサー・グリッド
グリッド
 ＭＰＩで起動する総プロセス数は、プロセッサー・グリッドと呼ばれる２次元プロセッサ構造で管理されます。
ばプ
造
８×１構成 例： ８プロセス

２×４構成P(0,0)
P(1,0)

203

P(0,1)
P(1,1)

P(0,2)
P(1,2)

P(0,3)
P(1,3)

４×２構成P(0 0)
P(0,0)

P(0 1)
P(0,1)

P(1,0)

P(1,1)

P(2,0)

P(2,1)

P(3,0)

P(3,1)

並列プログラミング講習会（試行）P(0,0)
P(1,0)

P(4,0)

P(2,0)

P(5,0)

P(3,0)

P(6 0)
P(6,0)
P(7,0)( )ScaLAPACKにおけるデータ分散
 ブロック幅ＭＢの２次元ブロックサイクリック分散となります 例：

８プロセス（2×4プロセッサー・グリッド）、ＭＢ=２

9

P(0,0)

P(0,1)

P(1,0)

P(1,1)

P(2,0)

P(2,1)

P(3,0)

P(3,1)

9

大域ビュー（行列Ａ）
大域ビ
（行列Ａ）
204

局所ビュー(sub(Ａ))
局所ビ
並列プログラミング講習会（試行）コンテキストの定義
コンテキスト：プロセッサ・グリッドに関する情報体
BLACS_GET( ICONTXT, WHAT, VAL )ＢＬＡＣＳの内部情報を取得する
ICONTXT: 入力、整数型、
入力 整数型
WHATで指定するコンテキストに設定を試みる値（引数）。
（引数が）無い時は無視される。
ＷＨＡＴ：入力、整数型、コンテキストに設定する内容
WHAT = 0 : デフォルトのシステムコンテキスト
WHAT = 1 : BLACS メッセージの ID 幅
WHAT = 2 : コンパイルされるBLACS デバックレベル
WHAT = 10: ICONTXTにより制御されるBLACSコンテキストを定義するために用いられるシステムコンテキストの参照
WHAT = 11: 現在使っているマルチリング構造のリング数
WHAT = 12: 現在使っている一般化した木構造の枝数

ＶＡＬ 出力 整数型 ンテキスト情報
ＶＡＬ：出力、整数型、コンテキスト情報
205

並列プログラミング講習会（試行）BLACSグリッドの定義
BLACS_GRIDINIT(BLACSGRIDINIT( ICONTXT, ORDER, NPROW, NPCOL )ICONTXT : 入力／出力、整数型
コンテキストORDER：入力、文字型＊１
プロセスをどのようにＢＬＡＣＳのグリッドに割り当てるか指定‘R’ : 行方向の自然なオーダリングを使う。
行方向の自然なオ ダリングを使う
‘C’ : 列方向の自然なオーダリングを使う。
そうでないなら：行方向の自然なオ ダリングを使う。
そうでないなら：行方向の自然なオーダリングを使う。

206

並列プログラミング講習会（試行）BLACSグリッドの情報入手
BLACS_GRIDINFOBLACSGRIDINFO(ICONTXT, NPROW, NPCOL, MYPROW, MYPCOL )ICONTXT : 入力、整数型
コンテキストNPROW：出力、整数型
プロセッサ グリッドの行数
プロセッサ・グリッドの行数
NPＣＯＬ：出力、整数型
プロセッサ グリッドの列数
プロセッサ・グリッドの列数
ＭＹPROW：出力、整数型
プロセッサ・グリッドにおける自分の行方向の認識番号
ＭＹPＣＯＬ：出力、整数型
プロセッサ・グリッドにおける自分の列方向の認識番号
207

並列プログラミング講習会（試行）データ分散を容易にする関数
PDELSET( A, IA, JA, DESCA, ALPHA )大域配列の添え字から、局所配列の適する場所にデータを収納する
Ａ：局所出力、倍精度型
収納すべき局所配列
ＩＡ：大域入力 整数型
ＩＡ：大域入力、整数型
大域的配列における１次元目の添え字
Ｊ 大域入力、整数型
ＪＡ：大域入力、整数型
大域的配列における２次元目の添え字
ＤＥＳＣＡ：大域かつ局所入力、整数型配列
局所配列Ａの記述子
所
述ＡＬＰＨＡ：局所入力、倍精度型
代入すべき値
208

並列プログラミング講習会（試行）データ分散を容易にする関数:例
大域配列ＡＭＡＴを、適切な局所行列Ａに代入、
代

DO J = 1, N
DO I = 1, N
CALL PDELSET( A, I, J, DESCA, AMAT(I,J))ENDDOCALL PDELSET( B, J, 1, DESCB, RHS(J))
ENDDO
大域ベクトルＲＨＳを、適切な局所ベクトルＢに代入

209

並列プログラミング講習会（試行）PDELSET関数による方法の注意点
 PDELSET関数を用いて行列を分散する方法は容易ですが、サンプルプログラムのように、大域行列Aを全プロセスで所サ
プ プ グ ム よう
大域行列 を全プ所有していると、メモリに関するスケーラビリティが無くなります。


つまり、実行できる行列サイズNが、１ノードのメモリ総量でまり 実行 きる行列サイズ が ノ ド メ リ総量決まる。実行ノード数を増加しても、Nを大きくできない。

 本格的な並列プログラムを書くためには

PDELSET関数を用いず
 局所ビューの観点で、直接sub(A)に値を代入する
局所ビの観点で 直接 b(A)に値を代入する
プログラムを書く必要があります。
 そのためには、ブロック・サイクリック分散方式の特徴をそのためには ブロ ク サイクリ ク分散方式の特徴を理解して、行列の添え字を計算する必要があります 本講習会のMPI基礎編で、関連技術の演習があります。
本講習会のMPI基礎編で 関連技術の演習があります
210

並列プログラミング講習会（試行）ScaLAPACK pdgesvのサンプルプログラムの実行（Fortran言語版）
 以下のコマンドを実行する$ cp /home/z30082/lecScaLAPACK.tar ./$ tar xvf lecScaLAPACK.tar$ cd sphere-ScaLAPACK$ cd F$ make$ qsub go.sh
go sh 実行が終了したら、以下を実行する
実行が終了したら 以下を実行する$ cat test.out
211

並列プログラミング講習会（試行）ScaLAPACK pdgesvのサンプルプログラムの実行（Fortran言語版）
 以下のような結果が見えれば成功
TIME[sec] = 0.309944152832031250E-005
MFLOPS = 351361040.384000003
ScaLAPACK Example Program -- Sep, 2010
Solvingg Ax=b where A is a 1000 byy 1000 matrix witha block size of 32
Running on 64 processes, where the process grid is 8 by 8
INFO code returned by PDGESV = 0
EPS = 627.423511391422494
連立一次方程式求解部分が実装されて無いので、時間が極端に短く、誤差ERRが大きくなっていますが、正常な動作です。
常な動作 す212

並列プログラミング講習会（試行）サンプルプログラムの説明(ScaLAPACK) ピュアＭＰＩ版です スレッド並列を利用したハイブリッドＭＰＩ版を利用するには、スレッド並列版をリンクしてコンパイルし、スレッド並列数を指定す指定する必要がありますがあす ６４プロセス実行（プロセッサー・グリッド８×８）、

問題サイズ1000、ブロック幅32に特化されています プロセス数、プロセッサー・グリッド、問題サイズ、ブロック幅を変更する場合、PARAMETER文の定数を変更する必要
があります

213

並列プログラミング講習会（試行）ScaLAPACK呼び出しオプション
 以下のオプションを付けます。

Fortran言語

mpif90f90 <プログラム名>
プ グ ム名 -noparallel-lscalapacklscalapack_sc
sc -llapackllapack_sc
sc-lblacsF77_sc -lblacsBASE_sc-lblacsC_sc
lbl
C-L/opt/itc/lib
L/ t/it /lib -lgoto
l t※要素並列化されていない（スレッド並列化されていない）ＬＡＰACKとＢＬＡＳを用いたScaLAPACKの呼び出し手順です。
※ＢＬＡＳは、GOTO BLASを用います。
※MPI通信を用いる版です※MPI通信を用いる版です。
214

並列プログラミング講習会（試行）演習課題1（ScaLAPACK）
1.

test.f のメイン手続きの連立一次方程式求解部
分を、ScaLAPACK pdgesv手続きの呼び出しにより、高速り、高速化してください 引数の並びに注意してください。
 デバッグとして、誤差（ERR）が十分小さいことを確認
デバッグとして 誤差（ERR）が十分小さいことを確認してください。

215

並列プログラミング講習会（試行）演習課題2（ScaLAPACK）
問題サイズを大きくして実行速度を比べてください。
並列実行数を64並列から変更して実行してください。

2.
3
3.
プロセッサー・グリッドを８×８から変更する必要があります。
必要があります上記の２、３ともに、データ分散（行列の確保のサイズ）の変更（コードの修正）が必須ですので 注意してください
ですので、注意してください。

216

並列プログラミング講習会（試行）演習課題３（ScaLAPACK）
大域配列matに大域的な係数を代入せず、
局所配列 A = sub(mat) に直接、各プロセスで必要な係数を代入する
プログラムを作成し、メモリに関するス
ケーラビリティがあるように、プログラムビ
プ グ
を改良してください。

4.

大域配列matの確保は不要となります

217

並列プログラミング講習会（試行）ScaLAPACK pdgesv回答
p g
回答

218

並列プログラミング講習会（試行）LAPCK dgesvの回答（Fortran言語版）

CALL PDGESV(&
NN, NRHS,&
A
A, IA
IA, JA,
JA DESCA,DESCA IPIV,
IPIV
&
B
B, IB
IB, JB,JB DESCB,DESCB&INFO )
219

並列プログラミング講習会（試行）お疲れさまでした220

並列プログラミング講習会（試行）