内容に関するご質問はida＠cc.u-tokyo.ac.jp
まで、お願いします。

第140回 お試しアカウント付き
並列プログラミング講習会「ライブラリ利用：科学技術計算の効率化入門」

スパコンと線形計算ライブラリ（BLAS, LAPACK）
東京大学情報基盤センター 特任准教授 伊田 明弘

1

チュートリアルの流れ
1.
2.
3.

現代のスパコンと並列計算
BLAS, LAPACKの説明プログラム実習Ⅰ
・BALS/DGEMMサンプルプログラム実行

4.

プログラム実習Ⅱ
・LAPACK/DGESVサンプルプログラム実行

2

お試しアカウント付き講習会

SMPクラスタ型並列計算システムの概念図■転送速度
Core占有cache

> > >

約10倍

共有cache約10倍

メインメモリ
約10倍

通信網

■容量
・Core占有cache：
数十～数百KB・共有cache：
数十MB
・メインメモリ：
数十GB
3

スーパーコンピュータの変遷■クラスタ型が現代では主流

Share(%)

MPP■クラスタ型計算システムの性能を発揮させるのは大変
Cluster

SMP
SingleProcessor’95

Constellations
’05
’10
’00
Years

Top 500 list of Supercomputers(http://www.top500.org)
スーパーコンピュータのアーキテクチャの変遷

・キャッシュの有効利用
・SIMDの活用
・計算負荷バランス
・データ競合回避
・計算順序保障・ハイパースレッドの活用
・通信量の削減
・効率的な通信パターン
・同期回数の削減
ライブラリを活用することでプログラミングコスト削減
4

SMPクラスタ型並列計算システムの例
クラスタ型並列計算システムの例
SMPクラスタ名
クラスタ名
SMP数
SMP間
全体
データ転送性能
Processor数
理論ピーク演算性能
SMP
メモリ容量
メモリ転送速度
Processor名
Core数
Processor 理論演算性能(Core)
Core占有cache(L2)共有cache(L3)ReedbushReedbush-U
＠東大基盤センタ
420

Oakbridge-CX
＠東大基盤センタ
1368

25Gbit/s × 4

100Gbit/s

21209.6 Gflops
256GB
153.6 GB/sec
Intel Xeon E5 TM
18
33.6 GFlops
256KB
45MB

2
4838.4 GFlops
192GB
281.6 GB/sec
Intel Xeon Platinum 8280 TM

28
86.4 GFlops
25MB
38.5MB■GB:Giga Byteの略, １GB=８Gbit■flops : １秒間に行われる浮動小数点演算の回数
計算機の性能を表す指標によく使用される
5

並列計算の必要性と種類■逐次コードで使える計算資源は、スパコン全体の極々一部
・スパコンは、普通のＣＰＵの集合に過ぎずちっともスーパーではない
・単体ＣＰＵの性能は、普通のパソコンと大差ない■並列計算の種類① 逐次コードの単純並列実行
・小規模計算を大量のパラメータに対して行う
・何も工夫することなく、最高の並列計算効率
・アルゴリズムの工夫やCasheの使い方などが高速化の鍵② 並列コード
・複数のＣＰＵを連携させて、大規模計算を行う
・並列化効率を向上させるため、うまく連携させる工夫が必要
・Casheの使い方など逐次コードで鍵となった技術もやはり重要

6

性能評価指標－台数効果
 台数効果

式： SP  TS / TP TS ：逐次の実行時間、
TP ：P台での実行時間
 S  P のとき、理想的な（ideal）速度向上
P


 並列化効率

式： EP  SP / P 100 [%] 飽和性能速度向上の限界 Saturation、「さちる」
Ｐ
7

お試しアカウント付き講習会

アムダールの法則
 逐次実行時間をK とする。
そのうち、並列化ができる割合を α とする。
 このとき、台数効果は以下のようになる。

SP  K /(K / P  K(1)) 1/( / P  (1)) 上記の式から、たとえ無限大の数のプロセッサを使っても（P→∞）、台数効果は、高々 １／（１－α） である。
（アムダールの法則）


8

並列化が全体の90％だと、無限個のプロセッサを使用しても、
１/(１-0.9) = １0 倍 にしかならない！
→高性能を達成するためには、少しでも並列化効率を上げる
実装をすることがとても重要である
お試しアカウント付き講習会

アムダールの法則の直観例
並列化できない部分（１ブロック） 並列化できる部分（８ブロック）●逐次実行●並列実行（４並列）●並列実行（８並列）

9＝88.8%が並列化可能

9/3=3倍

9/2=4.5倍 ≠ 6倍

お試しアカウント付き講習会

並列コードの種類
ＭＰＩプロセスとOpenMPスレッドを用いた並列化が二大メジャー手段■プロセス並列– MPI （Message Passing Interface））– HPF （High Performance Fortran）■スレッド並列– Open MP• ユーザが並列化指示行を記述– Open ACC– Pｔhread (POSIX スレッド)– CILK– Javaプロセスとスレッドの違い•メモリを意識するかどうかの違い
メモリを意識するかどうかの違い•別メモリは「プロセス」
別メモリは「プロセス」•同一メモリは「スレッド」
同一メモリは「スレッド」

10

スレッド並列とプロセス並列
スレッド並列とプロセス並列■スレッド並列メモリ

・ＳＭＰのように、メモリが多数のＣＰＵで共有されていることが前提

スレッド

スレッド

スレッド

スレッド

ＣＰＵ

・任意の情報へどのスレッドからも直接アクセス可能・通信網の先にある他のＳＭＰが持つ情報へはアクセス不能■プロセス並列メモリ

・メモリが多数のＣＰＵで共有されていたとしてもメモリ空間を分離

プ ロセス

通信

プ ロセス

プ ロセス

プ ロセス

ＣＰＵ

通信

・他のプロセスが持つ情報へは「通信」することでアクセス
・ 「通信」することにより、他のＳＭＰが持つ情報へもアクセス可能
11

MPI＋OpenMPハイブリッドプログラミングモデル■ＳＭＰクラスタを考えるとプロセス並列とスレッド並列の組み合わせが合理的・各MPI-processが、１つ以上,SMP-core数以下の固定数のOpenMP-Threadを持つ・演算の多くはOpenMP-Threadが共有メモリ下で実行する
・各MPI-processは独立したメモリ空間を有し，process間通信でデータをやりとりする
・各MPI-processが何Thread持つのが良いかは計算システムと実装で決まる■SMP-core数が４つの場合の組み合わせ数が４つの場合の組み合わせ
1process1thread(準pure-MPI)

ThreadThreadThreadprocess

Threadprocess

process
process
process

1process4thread
SMP

SMP

SMP

Thread1process2thread

ThreadThreadThreadprocess

ThreadThreadThreadprocess

Thread12

MPIの概要と特徴
 メッセージパッシング用のライブラリ規格の１つ• メッセージパッシングのモデルである• コンパイラの規格、特定のソフトウエアやライブラリを指すものではない

 分散メモリ型並列計算機（クラスタ）に対応可能• １プロセッサにおけるメモリサイズやファイルサイズの制約を打破可能• プロセッサ台数の多い並列システムを用いる実行に向く‐１プロセッサ換算で膨大な実行時間の計算を、短時間で処理可能

 異なるシステムへの移植が容易
API（Application Programming Interface）の標準化

 スケーラビリティ、性能が良い
通信処理をユーザが記述することによるアルゴリズムの最適化が可能

13

OpenMPの概要と特徴
 共有メモリ型並列計算機用のDirectiveの統一規格
 ディレクティヴ（指示行）の形で利用• プログラムを並列に実行するための動作をユーザーが明示• 挿入直後のループが並列化される等• コンパイラがサポートしていなければ，コメントとみなされる

 逐次コードからの並列化が比較的容易• 指示行（ディレクティヴ）を挿入するだけで手軽に並列化が可能• ループ単位など少しずつ段階的に並列化が可能

 OpenMPがMPIより簡単ということはない• OpenMP実行環境は，依存関係，衝突，デッドロック，競合条件，結果としてプログラムが誤った実行につながるような問題に関するチェックは要求されていない• プログラムが正しく実行されるよう構成するのはユーザーの責任14スパコン上でのライブラリ利用■ライブラリ・ソフトウェア使えばスパコンの性能を楽に引き出せる■並列化の各ステージごとに様々なソフトウェアが提供されている① SMPクラスタ対応パッケージソフト
・格子データや定義ファイルを用意するだけで大規模計算が可能
・ユーザは並列化や通信をあまり意識する必要がない
・OpenFOAM，FrontFlow, ppOpen-APPL/BEMなど② ＳＭＰクラスタ対応ライブラリ、クラスタ対応ライブラリ
・通信機能を持つルーチンを呼ぶことで大規模計算を行う
・ライブラリルーチンが行う計算以外の並列計算はユーザが行う・ＭＰＩ，Super-LU， ScaLAPACKなど③ SMP対応ライブラリ、逐次計算用ライブラリ
・逐次計算コードで高速化したい部分の機能を有するルーチンを呼ぶ
・複数ＳＭＰを使いたければ、通信をユーザが行う必要あり
・Intel MKL，富士通SSL2，BLAS， LAPACKなど
15

BLAS, LAPACK, ScaLAPACKの説明

16

お試しアカウント付き講習会

参考文献（１）

BLAS
http://www.netlib.org/blas/
2. LAPACK
http://www.netlib.org/lapack/3. ScaLAPACK
http://www.netlib.org/scalapack/1.

Netlib ：科学技術計算用ソフトウェア・リポジトリ（ＡＴ＆Ｔ、ベル研究所、テネシー大学、
オークリッジ国立研究所が共同管理）
17

お試しアカウント付き講習会

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

18

お試しアカウント付き講習会

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

お試しアカウント付き講習会

BLASサブルーチンの分類
 BLASでは、以下のように分類わけをして、サブルーチンの命名規則を統一
1.
2.
3.
4.

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

ベクトルとベクトルの演算 レベル２ BLAS： 行列とベクトルの演算 レベル３ BLAS： 行列と行列の演算
 なるべく高いレベルのルーチンを使用するのが効果的
20

お試しアカウント付き講習会

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

21

お試しアカウント付き講習会

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

22

お試しアカウント付き講習会

レベル３ BLAS
 行列-行列積などの演算
 例：

C←αAB+βC
 レベル２ BLASより更にデータ再利用を追及
2
 行列-行列積では、行列データ O(n ) に対して演算は O ( n 3 ) なので、データ再利用性が原理的に高い。
 行列積は、アルゴリズムレベルでもブロック化できる。
さらにデータの局所性を高めることができる。

23

お試しアカウント付き講習会

BLASの命名規則
 関数名：XYYYY

データ型
S:単精度、D：倍精度、C：複素、Z：倍精度複素 YYYY： 計算の種類
 X：

 レベル１：

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

お試しアカウント付き講習会

BLASルーチン一覧

25

お試しアカウント付き講習会

インタフェース例：DGEMM （１/４）
 DGEMM(TRANSA,TRANSB, M, N, K, ALPHA, A, LDA, B, LDB,
BETA, C, LDC)
C := alpha*op( A )*op( B ) + beta*C の計算をする op( X ) = X もしくは op( X ) = X’（Xの転置行列）


 引数


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

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

26

お試しアカウント付き講習会

インタフェース例：DGEMM （２/４）


M（入力） - INTEGER
 op( A )



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の転置を入れる。
27

お試しアカウント付き講習会

インタフェース例：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。
 TRANSA = ‘N’ or ‘n’のときは、k×nの配列の要素に行列Bを含まないといけない。そうでないときは、n×kの配列に行列Bの転置を入れる。


LDB（入力） - INTEGER
 行列Bの最初の次元数を入れる。 TRANSA = ‘N’ もしくは ‘n’ なら、LDA

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

28

お試しアカウント付き講習会

インタフェース例：DGEMM （４/４）
BETA （入力） - DOUBLE PRECISION スカラ値

BETAの値を設定する。

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


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

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

お試しアカウント付き講習会

ＢＬＡＳの問題点
 BLASの問題点

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


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

2.


3.

例：レベル１BLASを用いた行列・ベクトル積＜メモリ浪費＞や＜演算性能低下＞の要因に

ソースコードが読みにくくなる
 BLASのインタフェースを熟知しないと、かえって処理が理解できない30

お試しアカウント付き講習会

LAPACKとは
 Linear Algebra PACKage（線形代数パッケージ）
BLASを組み合わせて使用することで、より高度な線形代数計算を実現
 線形方程式系の解法
 固有値問題の解法
 特異値計算 最小二乗法 使いたいものを探すのが大変なほど

数多くのルーチンが用意されている
 HP： http://www.netlib.org/lapack/31

お試しアカウント付き講習会

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

32

お試しアカウント付き講習会

LAPACKルーチン一覧■先頭にデータ型(“s”,“d”,”c”,”z”)のいづれかを付けるとルーチン名になる
bdsdc bdsqr disna gbbrd gbcon gbequ gbequb gbrfs gbrfsx gbsv gbsvx gbsvxx gbtrf
gbtrs gebak gebal gebrd gecon geequ geequb gees geesx geev geevx gehrd gejsv gelqf
gels gelsd gelss gelsy geqlf geqp3 geqpf geqrf geqrfp gerfs gerfsx gerqf gesdd gesv gesvd
gesvj gesvx gesvxx getrf getri getrs ggbak ggbal gges ggesx ggev ggevx ggglm gghrd gglse
ggqrf ggrqf ggsvd ggsvp gtcon gtrfs gtsv gtsvx gttrf gttrs hgeqz hsein hseqr opgtr opmtr
orgbr orghr orglq orgql orgqr orgrq orgtr ormbr ormhr ormlq ormql ormqr ormrq
ormrz ormtr pbcon pbequ pbrfs pbstf pbsv pbsvx pbtrf pbtrs pftrf pftri pftrs pocon
poequ poequb porfs porfsx posv posvx posvxx potrf potri potrs ppcon ppequ pprfs
ppsv ppsvx pptrf pptri pptrs pstrf ptcon pteqr ptrfs ptsv ptsvx pttrf pttrs sbev sbevd
sbevx sbgst sbgv sbgvd sbgvx sbtrd sfrk spcon spev spevd spevx spgst spgv spgvd spgvx
sprfs spsv spsvx sptrd sptrf sptri sptrs stebz stedc stegr stein stemr steqr sterf stev
stevd stevr stevx sycon syequb syev syevd syevr syevx sygst sygv sygvd sygvx syrfs
syrfsx sysv sysvx sysvxx sytrd sytrf sytri sytrs tbcon tbrfs tbtrs tfsm tftri tfttp tfttr
tgevc tgexc tgsen tgsja tgsna tgsyl tpcon tprfs tptri tptrs tpttf tpttr trcon trevc trexc
trrfs trsen trsna trsyl trtri trtrs trttf trttp tzrzf

33

お試しアカウント付き講習会

インタフェース例：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の次元数。 N >= 0。

34

お試しアカウント付き講習会

インタフェース例：DGESV （２/３）
NRHS (入力) – INTEGER 右辺ベクトルの数。行列Bの次元数。 NRHS >= 0。

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



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

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

35

お試しアカウント付き講習会

インタフェース例：DGESV （３/３）
B (入力／出力) – DOUBLE PRECISION, DIMENSION(:,:)
 入力時は、右辺ベクトルの

N×NRHS 行列Bを入れる。
 出力時は、もし、INFO = 0 なら、N×NRHS行列である解行列Xが戻る。


LDB (入力) －INTEGER 配列Bの最初の次元の大きさ。 LDB >= max(1,N)。



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

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

36

お試しアカウント付き講習会

C言語でＢＬＡＳ, LAPACKを使う時の注意■BLAS, LAPACKはFortranで書かれている・ルーチン名が変えれている（DGEMM⇒dgemm_）
・引数はポインタ引き渡し
・配列インデックスは１-オリジン
・２次元配列はColomnMajor

37

ColomnMajor

Row Major

Fortran

C言語
お試しアカウント付き講習会

プログラム実習Ⅰ(BLAS/DGEMM)

38

お試しアカウント付き講習会

サンプルプログラムの説明（BLAS）
 １コア（逐次）実行版です スレッド並列化には、スレッド並列化版をリンクしてコンパイルの上、並列実行数の指定をする必要があります BLASスレッドでは、OpenMPのスレッド数指定法と同じ方法で、実行するスレッド数が指定できます ジョブスクリプト内に使用スレッド情報を記載#PJM --omp thread=10

39

お試しアカウント付き講習会

実習を行うあたっての注意点( BLAS DGEMM)
 サンプルプログラムのファイル名

Mat-Mat-BLAS-OBCX.tar
 ジョブスクリプトファイルmat-mat-blas.bashのキュー名を、
lecture から tutorial に変更してから
pjsub してください lecture : 実習時間外のキュー（同時実行数１）
 tutorial : 実習時間内のキュー（同時実行数4+）40

お試しアカウント付き講習会

中BLAS DGEMMのサンプルプログラムの実行（Ｃ言語版/Fortran言語版共通）
 $ cd /work/gt00/t/work/gt00/tXXXXX (XXXXXは数値)// 実行してWork領域に移動する
 以下のコマンドを実行する$ cp /work/gt00/z30107/Mat-Mat-BLAS-OBCX.tar ./$ tar xvf Mat-Mat-BLAS-OBCX.tar$ cd Mat-Mat-BLAS$ cd C //C言語の人$ cd F //Fortran言語の人$ make$ pjsub mat-mat-blas.bash
 実行が終了したら、以下を実行する$ cat mat-mat-blas.bash.oXXXXXX (XXXXXXは数値)
41

お試しアカウント付き講習会

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

N = 9100
Mat-Mat time = 190.748147 [sec.]
7901.214364 [MFLOPS]
OK!

42

お試しアカウント付き講習会

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

N =
9100
Mat-Mat time[sec.] = 103.923363756388
MFLOPS = 14502.4366198134
OK!

43

お試しアカウント付き講習会

サンプルプログラムの説明（Ｃ言語） #define

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

お試しアカウント付き講習会

Fortran言語のサンプルプログラムの注意
 行列サイズNNの宣言は、以下のファイルにあります。

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

integer NN
parameter (NN=91００)

45

お試しアカウント付き講習会

時間計測方法（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();t0 = t2 - t1;
ierr = MPI_Reduce(&t0, &t_w, 1,
MPI_DOUBLE,MPI_MAX, 0,
MPI_COMM_WORLD);

46

各プロセッサーで、ｔ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 MPI_BARRIER(MPI_COMM_WORLD, ierr)t2 = MPI_WTIME(ierr)
t0 = t2 - t1
call MPI_REDUCE(t0, t_w, 1,&
MPI_DOUBLE_PRECISION,&MPI_MAX, 0, MPI_COMM_WORLD, ierr)
47

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

ＯBCXでのBLAS呼び出し
 インテルコンパイラから、BLAS (Intel MKLライブラリ)を呼び出す場合、以下のオプションを付けます。
BLASが逐次(1コア実行）
mpiifort -mkl <プログラム名> ：Fｏｒｔｒａｎ言語
mpiicc -mkl <プログラム名> ：Ｃ言語
BLASがスレッド実行
mpiifort –mkl=parallel -qopenmp <プログラム名>:Ｆｏｒｔｒａｎ言語
mpiicc –mkl=parallel -qopenmp <プログラム名>：Ｃ言語
48

お試しアカウント付き講習会

演習課題１（BLAS）
1.

MyMatMat関数の行列積コードを、BLAS のDEGMM
ルーチンを呼び出して高速化してください コードで行列-行列積の実行部をコメントアウト
C版：ファイルmat-mat-blas.cの100行目
Fortran版：ファイルmat-mat-blas.fの110行目
 DEGMMの引数の並びに注意して記述

C版：ファイルmat-mat-blas.cの106行目
Fortran版：ファイルmat-mat-blas.fの116行目
 C言語は、以下に注意してください 関数名： DGEMM_
 全ての引数が ポインタ引き渡し
49

お試しアカウント付き講習会

演習課題１のヒント
 倍精度演算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, ALPHA, A, N, B, N, BETA, C, N )Aが転置しているか
Bが転置しているか

alpha
の値
Aの
アドレス

Mの大きさ
Nの大きさ

Aの１次元目
の要素数

Bの
アドレス

Bの１次元目
の要素数

Kの大きさ

50

beta
の値

お試しアカウント付き講習会

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

演習課題２（BLAS）
2-1 スレッド並列コードの実行
スレッド並列版のDGEMMを呼び出し、いろいろなスレッド数で実行してください
ヒント1：３ページ前のスライドを見て、ファイル「Makefile」を開き、「FLAGS=」の右辺を書き換える
ヒント2：ファイル「mat-mat-blas.bash」を開き、「thread=1」の右辺を変更してジョブを投入する

2-2 台数効果による性能向上の評価
逐次実行に比べ、どれだけ高速化されたか、性能評価をしてください
51

お試しアカウント付き講習会

BLAS DGEMM回答

52

お試しアカウント付き講習会

BLAS DGEMM演習問題１の回答（C言語版）

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

お試しアカウント付き講習会

BLAS DGEMM演習問題１の回答（Fotran言語版）
double precision ALPHA,BETA
ALPHA=1.0d0
BETA=1.0d0
CALL DGEMM('N', 'N', n, n, n, ALPHA,&
A, n, B, n, BETA, C, n)

54

お試しアカウント付き講習会

プログラム実習Ⅱ(LAPACK)

55

お試しアカウント付き講習会

内容
 サンプルプログラムの説明
 LAPACK/DGESVサンプルプログラム実行（C言語、Fortran言語）

56

お試しアカウント付き講習会

サンプルプログラムで計算する問題■格子状に並んだ粒子間の熱伝導問題■詳細資料
http://nkl.cc.u-tokyo.ac.jp/FEMintro/
57

お試しアカウント付き講習会

サンプルプログラム・シリアル版の計算手順
 制御データ入力
 粒子生成 係数マトリクス計算

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


 出力

58

MicroAVS用
お試しアカウント付き講習会

サンプルプログラムの説明
 誤差ERRの計算は，10×10×10の問題に特化されています．
それ以外の問題サイズでは機能しません．
 test_org.c、test_org.fは、ライブラリを用いず、手製の反復解法（CG法）を用いて、線形方程式系の求解を行うコードです。

59

お試しアカウント付き講習会

サンプルプログラムの実行（LAPACK/DGESV）

60

お試しアカウント付き講習会

サンプルプログラムの説明（ LAPACK版）
 １コア（逐次）実行版です スレッド並列化版は、スレッド並列化版をリンクしてコンパイルの上、並列実行数の指定をする必要があります。
 LAPACKはスレッド並列化のみ対応しているため、複数ノードを通信しつつ利用することはできません ノード内のみＬＡＰＡＣＫを使い、ノード間はＭＰＩで並列化するような使い方は可能です たとえば、領域分割法を利用し、ノード内のみＬＡＰＡＣＫで連立一次方程式を解く場合

61

お試しアカウント付き講習会

実習を行うあたっての注意点(LAPACK dgesv)
 サンプルプログラムのファイル名

lecLAPACK-OBCX.tar
 ジョブスクリプトファイルgo.sh中のキュー名をlecture から tutorial に変更してから
pjsub してください。
 lecture : 実習時間外のキュー（同時実行数１）
 tutorial : 実習時間内のキュー（同時実行数4+）62

お試しアカウント付き講習会

LAPACK dgesvのサンプルプログラムの実行（Ｃ言語版/Fortran言語版共通）
 /work/gt00/t00***に移動する
 以下のコマンドを実行する$ cp /work/gt00/z30107/lecLAPACK-OBCX.tar ./$ tar xvf lecLAPACK-OBCX.tar$ cd sphere-LAPACK$ cd C //C言語の人 $ cd F //Fortran言語の人$ make$ pjsub go.bash
 実行が終了したら、以下を実行する$ 63cat go.bash.oXXXXXX(XXXXXXは数字)
お試しアカウント付き講習会

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

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

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

お試しアカウント付き講習会

LAPACK dgesvのサンプルプログラムの実行（Fortran言語版）
 以下のような結果が見えれば成功
TIME[sec] = 6.053596735000610E-008
MFLOPS = 11070575549.3297
991 -2.272792E+00
992 -2.276715E+00
993 -2.288410E+00
994 -2.307670E+00
995 -2.334166E+00
996 -2.367479E+00
997 -2.407125E+00
998 -2.452584E+00

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

999 -2.503330E+00
1000 -2.558846E+00
ERR = 627.423511391423
65

お試しアカウント付き講習会

OBCXでのLAPACK呼び出し
 インテルコンパイラから、LAPACK (Intel MKLライブラリ)を呼び出す場合、Makefileに以下のオプションを付けますLAPACKが逐次(1コア実行）
mpiifort -mkl <プログラム名> ：Fｏｒｔｒａｎ言語
mpiicc -mkl <プログラム名> ：Ｃ言語
LAPACKがスレッド実行
mpiifort –mkl=parallel -qopenmp <プログラム名>:Ｆｏｒｔｒａｎ言語
mpiicc –mkl=parallel -qopenmp <プログラム名>：Ｃ言語66

お試しアカウント付き講習会

演習課題（LAPACK）
1.

test.c もしくは test.f のメイン関数（手続き）の連立一次方程式の求解を行う部分（test.c:159行目、test.f : 114行目）に、
LAPACK dgesvルーチンを呼び出してください 行列：ATAN, サイズ：N, 右辺ベクトル：RHS, ピボット:PIV, 引数の並びに注意してください C言語は、dgesv のFortran手続きの呼び出しになります。以下に注意してください 関数名：

dgesv_

 全ての引数がポインタ引き渡し
 誤差（ERR）が十分小さいことを確認してください。
67

お試しアカウント付き講習会

演習課題1 LAPACK dgesv回答

68

お試しアカウント付き講習会

演習課題 LAPCK dgesvの回答（C言語版）

dgesv_(&nn, &inc, amat2, &nn,piv, rhs, &nn, &info);69

お試しアカウント付き講習会

演習課題 LAPCK dgesvの回答（Fotran言語版）

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

70

お試しアカウント付き講習会

アンケートへ回答のお願い
 第140回講習会「科学技術計算の効率化入門」アンケートQRコード71

お疲れさまでした72

お試しアカウント付き講習会