科学技術計算のためのマルチコアプログラミング入門
Fortran編
第Ⅲ部：OpenMPによる並列化+演習

中島研吾
東京大学情報基盤センター

OMP-3

1

OpenMP並列化• L2-solをOpenMPによって並列化する。
– 並列化にあたってはスレッド数を「PEsmpTOT」によってプログラム内で調節できる方法を適用する• 基本方針– 同じ「色」（または「レベル」）内の要素は互いに独立，したがって並列計算（同時処理）が可能

OMP-32

4色，4スレッドの例
初期メッシュ
57 58 59 60 61 62 63 64
49 50 51 52 53 54 55 56
41 42 43 44 45 46 47 48
33 34 35 36 37 38 39 40
25 26 27 28 29 30 31 32
17 18 19 20 21 22 23 24
9 10 11 12 13 14 15 16
1 2 3 4 5 6 7 8

OMP-3

3

4色，4スレッドの例
初期メッシュ
57 58 59 60 61 62 63 64
49 50 51 52 53 54 55 56
41 42 43 44 45 46 47 48
33 34 35 36 37 38 39 40
25 26 27 28 29 30 31 32
17 18 19 20 21 22 23 24
9 10 11 12 13 14 15 16
1 2 3 4 5 6 7 8

OMP-344色，4スレッドの例
色の順に番号付け
45 61 46 62 47 63 48 64
13 29 14 30 15 31 16 32
41 57 42 58 43 59 44 60
9 25 10 26 11 27 12 28
37 53 38 54 39 55 40 56
5 21 6 22 7 23 8 24
33 49 34 50 35 51 36 52
1 17 ２ 18 ３ 19 ４ 20

OMP-35

4色，4スレッドの例
同じ色の要素は独立：並列計算可能
番号順にスレッドに割り当てる
thread #3
thread #2
thread #1
thread #0

45 61 46 62 47 63 48 64
13 29 14 30 15 31 16 32
41 57 42 58 43 59 44 60
9 25 10 26 11 27 12 28
37 53 38 54 39 55 40 56
5 21 6 22 7 23 8 24
33 49 34 50 35 51 36 52
1 17 ２ 18 ３ 19 ４ 20

OMP-3

6

ファイルコピー：FX10>$ cd <$O-TOP>>$ cp /home/ss/aics60/C/multicore-c.tar .
>$ cp /home/ss/aics60/F/multicore-f.tar .
>$ tar xvf multicore-c.tar>$ tar xvf multicore-f.tar>$ cd multicore
以下のディレクトリが出来ていることを確認
L3 stream
これらを以降<$O-L3>, <$O-stream>OMP-3

7

プログラムのありか on FX10
• 所在
– <$O-L3>/src，<$O-L3>/run• コンパイル，実行方法– 本体• cd <$O-L3>/src• make• <$O-L3>/run/L3-sol（実行形式）– コントロールデータ• <$O-L3>/run/INPUT.DAT– 実行用シェル• <$O-L3>/run/go1.shOMP-38

実行例% cd <$O-L3>
% ls
run
src

src0

reorder0% cd src% make% cd ../run% ls L3-sol
L3-sol% <modify “INPUT.DAT”>% <modify “go1.sh”>% pjsub go1.shOMP-39プログラムの実行
プログラム，必要なファイル等

L3-solポアソン方程式
ソルバー

INPUT.DAT
制御ファイルtest.inp
結果ファイル（ParaView）

OMP-310制御データ（INPUT.DAT）
100 100 100
1.00e-00 1.00e-00 1.00e-00
1.0e-08
16
100

NX/NY/NZ
DX/DY/DZEPSICCGPEsmpTOTNCOLORtot

変数名

型

内

容

NX，NY，NZ

整数

DX，DY，DZ

倍精度実数

各要素の3辺の長さ（X，Y，Z）

EPSICCG

倍精度実数

収束判定値

PEsmpTOT

整数

データ分割数NCOLORtot

整数

Ordering手法と色数≧2：MC法（multicolor），色数=0：CM法（Cuthill-Mckee）=-1：RCM法（Reverse Cuthill-Mckee）≦-2：CM-RCM法

各方向の要素数OMP-311

go1.sh#!/bin/sh#PJM -L "node=1"#PJM -L "elapse=05:00"#PJM -L “rscgrp=small”#PJM -j#PJM -o “test.lst"

または “rscgrp=school”

export OMP_NUM_THREADS=16./L3-solPEsmpTOTと一致させる

12

ジョブスクリプト• <$O-L3>/run/go1.sh• スケジューラへの指令 ＋ シェルスクリプト#!/bin/sh#PJM -L “node=1“#PJM -L “elapse=00:05:00“#PJM -L “rscgrp=school“#PJM -j#PJM -o “test.lst“

ノード数
実行時間
実行キュー名
標準出力ファイル名

export OMP_NUM_THREADS=16./L3-sol

実行ファイル名

13

ジョブ投入，確認等• ジョブの投入• ジョブの確認• ジョブの取り消し・強制終了• キューの状態の確認• キューの詳細構成• 実行中のジョブ数• 同時実行・投入可能数

pjsub スクリプト名
pjstatpjdel ジョブID
pjstat --rsc
pjstat --rsc –x
pjstat --rsc –b
pjstat --limit[z30088@oakleaf-fx-6 S2-ref]$ pjstat
Oakleaf-FX scheduled stop time: 2012/09/28(Fri) 09:00:00 (Remain: 31days 20:01:46)
JOB_ID334730

JOB_NAME
go.sh

STATUS PROJECTRUNNING gt61

RSCGROUP
lectureSTART_DATE
08/27 12:58:08

ELAPSE
00:00:05

TOKEN NODE:COORD
0.0 1

OMP-3• L2-solへのOpenMPの実装• 実行例• 最適化＋演習

14

OMP-315

L2-solにOpenMPを適用• ICCGソルバーへの適用を考慮すると• 内積，DAXPY，行列ベクトル積– もともとデータ依存性無し ⇒ straightforwardな適用可能• 前処理（修正不完全コレスキー分解，前進後退代入）– 同じ色内は依存性無し ⇒ 色内では並列化可能

OMP-316

実はこのようにしてDirectiveを直接挿入しても良いのだが・・・（1/2）!$omp parallel do private(i,VAL,k)do i = 1, N
VAL= D(i)*W(i,P)
do k= indexL(i-1)+1, indexL(i)VAL= VAL + AL(k)*W(itemL(k),P)
enddo
do k= indexU(i-1)+1, indexU(i)VAL= VAL + AU(k)*W(itemU(k),P)
enddo
W(i,Q)= VAL
enddo!$omp end parallel do• スレッド数をプログラムで制御できるようにしてみよう• GPU，メニィコアではこのままの方が良い場合もあるOMP-317

実はこのようにしてDirectiveを直接挿入しても良いのだが・・・（2/2）
do icol= 1, NCOLORtot!$omp parallel do private (i, VAL, k)
do i= COLORindex(icol-1)+1, COLORindex(icol)
VAL= D(i)do k= indexL(i-1)+1, indexL(i)VAL= VAL - (AL(k)**2) * DD(itemL(k))enddo
DD(i)= 1.d0/VAL
enddo!$omp end parallel doenddo• スレッド数をプログラムで制御できるようにしてみよう• GPU，メニィコアではこのままの方が良い場合もあるOMP-318

ICCG法の並列化：OpenMP• 内積：OK• DAXPY：OK• 行列ベクトル積：OK• 前処理

OMP-319

プログラムの構成
program MAIN
use STRUCT
use PCGuse solver_ICCG_mc
implicit REAL*8 (A-H,O-Z)real(kind=8), dimension(:), allocatable :: WK
call INPUT
call POINTER_INIT
call BOUNDARY_CELL
call CELL_METRICS
call POI_GEN
PHI=

0.d0

call solve_ICCG_mc( ICELTOT, NPL, NPU, indexL, itemL, indexU, itemU, D,
BFORCE, PHI, AL, AU, NCOLORtot, PEsmpTOT,
SMPindex, SMPindexG, EPSICCG, ITR, IER)&&&

allocate (WK(ICELTOT))
do ic0= 1, ICELTOT
icel= NEWtoOLD(ic0)
WK(icel)= PHI(ic0)
enddo
do icel= 1, ICELTOTPHI(icel)= WK(icel)
enddo
call OUTUCD
stop
end

結果（PHI）をもとの番号付けにもどす&&&OMP-320

module STRUCT
module STRUCT
use omp_lib
include 'precision.inc'

ICELTOT： 要素数
ICELTOTp: = ICELTOT
N：
節点数

!C
!C-- METRICs & FLUX
integer (kind=kint) :: ICELTOT, ICELTOTp, N
integer (kind=kint) :: NX, NY, NZ, NXP1, NYP1, NZP1, IBNODTOT
integer (kind=kint) :: NXc, NYc, NZc
real (kind=kreal) ::DX, DY, DZ, XAREA, YAREA, ZAREA, RDX, RDY, RDZ,
RDX2, RDY2, RDZ2, R2DX, R2DY, R2DZ&&

real (kind=kreal), dimension(:), allocatable ::VOLCEL, VOLNOD, RVC, RVN&

integer (kind=kint), dimension(:,:), allocatable ::XYZ, NEIBcell

&

&
&
&

NX,NY,NZ： x,y,z方向要素数
NXP1,NYP1,NZP1：
x,y,z方向節点数
IBNODTOT： NXP1×NYP1
XYZ(ICELTOT,3)： 要素座標（後述）
NEIBcell(ICELTOT,6):
隣接要素（後述）&!C!C-- BOUNDARYs
integer (kind=kint) :: ZmaxCELtot
integer (kind=kint), dimension(:), allocatable :: BC_INDEX, BC_NOD
integer (kind=kint), dimension(:), allocatable :: ZmaxCEL

境界条件関連： Z=Zmax!C!C-- WORK
integer (kind=kint), dimension(:,:), allocatable :: IWKX
real(kind=kreal),
dimension(:,:), allocatable :: FCV
integer (kind=kint) :: PEsmpTOT
end module STRUCT

PEsmpTOT：スレッド数

OMP-321

module PCG（これまでとの相違点）
module PCG
integer, parameter :: N2= 256
integer :: NUmax, NLmax, NCOLORtot, NCOLORk, NU, NL
integer :: NPL, NPU
integer :: METHOD, ORDER_METHOD
real(kind=8) :: EPSICCG
real(kind=8), dimension(:), allocatable :: D, PHI, BFORCE
real(kind=8), dimension(:), allocatable :: AL, AU
integer, dimension(:), allocatable :: INL, INU, COLORindex
integer, dimension(:), allocatable :: SMPindex, SMPindexG
integer, dimension(:), allocatable :: OLDtoNEW, NEWtoOLD
integer, dimension(:,:), allocatable :: IAL, IAU
integer, dimension(:), allocatable :: indexL, itemL
integer, dimension(:), allocatable :: indexU, itemU
end module PCG

NCOLORtot
COLORindex(0:NCOLORtot)

色数
各色に含まれる要素数のインデックス（COLORindex(icol)-COLORindex(icol-1)）

SMPindex(0:NCOLORtot*PEsmpTOT) スレッド用配列（後述）
SMPindexG(0:PEsmpTOT)
OLDtoNEW, NEWtoOLD

Coloring前後の要素番号対照表変数表（1/2）

OMP-3

配列・変数名

型

内 容

D(N)

R

対角成分，（N：全メッシュ数）

BFORCE(N)R

右辺ベクトルPHI(N)

R

未知数ベクトル

indexL(0:N)

I

各行の非零下三角成分数（CRS）

indexU(0:N)

I

各行の非零上三角成分数（CRS）

NPL

I

非零下三角成分総数（CRS）

NPU

I

非零上三角成分総数（CRS）

itemL(NPL)I

非零下三角成分（列番号）（CRS）

itemU(NPU)I

非零下三角成分（列番号）（CRS）

AL(NPL)R

非零下三角成分（係数）（CRS）

AU(NPL)R

非零上三角成分（係数）（CRS）

NL,NU

I

各行の非零上下三角成分の最大数（ここでは6）

INL(N)

I

各行の非零下三角成分数

INU(N)

I

各行の非零上三角成分数

IAL(NL,N)I

各行の非零下三角成分に対応する列番号

IAU(NU,N)I

各行の非零上三角成分に対応する列番号22

OMP-323

変数表（2/2）
配列・変数名型

内 容

NCOLORtot

I 入力時にはOrdering手法（≧2：MC，=0：
CM，=-1：RCM,-2≧：CMRCM），
最終的には色数，レベル数が入る

COLORindex(0:NCOLORtot)

NEWtoOLD(N)I 各色，レベルに含まれる要素数の一次元圧縮配列，
COLORindex(icol-1)+1から
COLORindex(icol)までの要素がicol番
目の色（レベル）に含まれる。
I 新番号⇒旧番号への参照配列

OLDtoNEW(N)I 旧番号⇒新番号への参照配列

PEsmpTOTI スレッド数

SMPindex(0:NCOLORtot*PEsmpTOT)

I スレッド用補助配列（データ依存性がある
ループに使用）SMPindexG(0:PEsmpTOT)

I スレッド用補助配列（データ依存性が無い
ループに使用）OMP-324

プログラムの構成
program MAIN
use STRUCT
use PCGuse solver_ICCG_mc
implicit REAL*8 (A-H,O-Z)real(kind=8), dimension(:), allocatable :: WK
call INPUT
call POINTER_INIT
call BOUNDARY_CELL
call CELL_METRICS
call POI_GEN
PHI=

0.d0

call solve_ICCG_mc( ICELTOT, NPL, NPU, indexL, itemL, indexU, itemU, D,
BFORCE, PHI, AL, AU, NCOLORtot, PEsmpTOT,
SMPindex, SMPindexG, EPSICCG, ITR, IER)

&
&
&

&
&
&

OMP-325

input「IPNUT.DAT」の読み込み!C!C***!C*** INPUT
!C***
!C
!C
INPUT CONTROL DATA
!Csubroutine INPUT
use STRUCT
use PCG
implicit REAL*8 (A-H,O-Z)
character*80 CNTFIL!C!C-- CNTL. file
open (11, file='INPUT.DAT', status='unknown')
read (11,*) NX, NY, NZ
read (11,*) DX, DY, DZ
read (11,*) EPSICCG
read (11,*) PEsmpTOT
read (11,*) NCOLORtot
close (11)!C===
returnend

100 100 100
1.00e-02 5.00e-02 1.00e-02
1.00e-08
16
100• PEsmpTOT– OpenMPスレッド数• NCOLORtot– 色数– 「＝0」の場合はCM– 「＝-1」の場合はRCM– 「≦-2」の場合はCMRCM
NX/NY/NZ
DX/DY/DZEPSICCGPEsmpTOTNCOLORtot

OMP-326

cell_metrics!C!C***!C*** CELL_METRICS!C***!C
subroutine CELL_METRICS
use STRUCT
use PCG
implicit REAL*8 (A-H,O-Z)計算に必要な諸パラメータ!C!C-- ALLOCATE
allocate (VOLCEL(ICELTOT))
allocate (RVC(ICELTOT))!C!C-- VOLUME, AREA, PROJECTION etc.
XAREA= DY * DZYAREA= DX * DZZAREA= DX * DY

XAREA
DZ

RDX= 1.d0 / DX
RDY= 1.d0 / DY
RDZ= 1.d0 / DZ
RDX2= 1.d0 / (DX**2)
RDY2= 1.d0 / (DY**2)
RDZ2= 1.d0 / (DZ**2)
R2DX= 1.d0 / (0.50d0*DX)
R2DY= 1.d0 / (0.50d0*DY)
R2DZ= 1.d0 / (0.50d0*DZ)
V0= DX * DY * DZ
RV0= 1.d0/V0
VOLCEL= V0
RVC
= RV0
return
end

DY
z

DX

y
x

OMP-327

プログラムの構成
program MAIN
use STRUCT
use PCGuse solver_ICCG_mc
implicit REAL*8 (A-H,O-Z)real(kind=8), dimension(:), allocatable :: WK
call INPUT
call POINTER_INIT
call BOUNDARY_CELL
call CELL_METRICS
call POI_GEN
PHI=
&
&
&

0.d0

call solve_ICCG_mc( ICELTOT, NPL, NPU, indexL, itemL, indexU, itemU, D,
BFORCE, PHI, AL, AU, NCOLORtot, PEsmpTOT,
SMPindex, SMPindexG, EPSICCG, ITR, IER)&&&OMP-328

poi_gen（1/9）
subroutine POI_GEN
use STRUCT
use PCG
implicit REAL*8 (A-H,O-Z)!C!C-- INIT.
nn = ICELTOT
nnp= ICELTOTp
NU= 6
NL= 6
allocate (BFORCE(nn), D(nn), PHI(nn))
allocate (INL(nn), INU(nn), IAL(NL,nn), IAU(NU,nn))
PHI= 0.d0D= 0.d0
BFORCE= 0.d0
INL= 0
INU= 0
IAL= 0
IAU= 0OMP-3!C!C +--------------+!C | CONNECTIVITY |!C +--------------+!C===
do icel= 1, ICELTOT
icN1= NEIBcell(icel,1)
icN2= NEIBcell(icel,2)
icN3= NEIBcell(icel,3)
icN4= NEIBcell(icel,4)
icN5= NEIBcell(icel,5)
icN6= NEIBcell(icel,6)29

poi_gen（2/9）
NEIBcell(icel,6)
NEIBcell(icel,4)if (icN5.ne.0.and.icN5.le.ICELTOT) then
icou= INL(icel) + 1
IAL(icou,icel)= icN5INL(icel)= icou
endif
if (icN3.ne.0.and.icN3.le.ICELTOT) then
icou= INL(icel) + 1
IAL(icou,icel)= icN3
INL(icel)= icou
endif
if (icN1.ne.0.and.icN1.le.ICELTOT) then
icou= INL(icel) + 1
IAL(icou,icel)= icN1
INL(icel)= icou
endif
if (icN2.ne.0.and.icN2.le.ICELTOT) then
icou= INU(icel) + 1
IAU(icou,icel)= icN2
INU(icel)= icou
endif
if (icN4.ne.0.and.icN4.le.ICELTOT) then
icou= INU(icel) + 1
IAU(icou,icel)= icN4INU(icel)= icou
endif

!C===

if (icN6.ne.0.and.icN6.le.ICELTOT) then
icou= INU(icel) + 1
IAU(icou,icel)= icN6
INU(icel)= icou
endif
enddo

NEIBcell(icel,1)

NEIBcell(icel,2)

NEIBcell(icel,3)
NEIBcell(icel,5)

下三角成分
NEIBcell(icel,5)= icel – NX*NY
NEIBcell(icel,3)= icel – NX
NEIBcell(icel,1)= icel – 1

OMP-3!C!C +--------------+!C | CONNECTIVITY |!C +--------------+!C===
do icel= 1, ICELTOT
icN1= NEIBcell(icel,1)
icN2= NEIBcell(icel,2)
icN3= NEIBcell(icel,3)
icN4= NEIBcell(icel,4)
icN5= NEIBcell(icel,5)
icN6= NEIBcell(icel,6)30

poi_gen（2/9）
NEIBcell(icel,6)
NEIBcell(icel,4)if (icN5.ne.0.and.icN5.le.ICELTOT) then
icou= INL(icel) + 1
IAL(icou,icel)= icN5INL(icel)= icou
endif
if (icN3.ne.0.and.icN3.le.ICELTOT) then
icou= INL(icel) + 1
IAL(icou,icel)= icN3
INL(icel)= icou
endif
if (icN1.ne.0.and.icN1.le.ICELTOT) then
icou= INL(icel) + 1
IAL(icou,icel)= icN1
INL(icel)= icou
endif
if (icN2.ne.0.and.icN2.le.ICELTOT) then
icou= INU(icel) + 1
IAU(icou,icel)= icN2
INU(icel)= icou
endif
if (icN4.ne.0.and.icN4.le.ICELTOT) then
icou= INU(icel) + 1
IAU(icou,icel)= icN4INU(icel)= icou
endif

!C===

if (icN6.ne.0.and.icN6.le.ICELTOT) then
icou= INU(icel) + 1
IAU(icou,icel)= icN6
INU(icel)= icou
endif
enddo

NEIBcell(icel,1)

NEIBcell(icel,2)

NEIBcell(icel,3)
NEIBcell(icel,5)

上三角成分
NEIBcell(icel,2)= icel + 1
NEIBcell(icel,4)= icel + NX
NEIBcell(icel,6)= icel + NX*NYOMP-331

poi_gen（3/9）!C!C +---------------+!C | MULTICOLORING |!C +---------------+!C===
allocate (OLDtoNEW(ICELTOT), NEWtoOLD(ICELTOT))
allocate (COLORindex(0:ICELTOT))111

並べ替えの実施：
NCOLORtot > 1：Multicolor
NCOLORtot = 0：CM
NCOLORtot =-1：RCM
NCOLORtot <-1：CM-RCM

continue
write (*,'(//a,i8,a)') 'You have', ICELTOT, ' elements.'
write (*,'( a)') 'How many colors do you need ?'write (*,'( a)') ' #COLOR must be more than 2 and'write (*,'( a,i8 )') ' #COLOR must not be more than', ICELTOT
write (*,'( a)') 'CM if #COLOR .eq. 0'
write (*,'( a)') ' RCM if #COLOR .eq.-1'
write (*,'( a)') 'CMRCM if #COLOR .le.-2'
write (*,‘( a)’) ‘=>’
if (NCOLORtot.gt.0) thencall MC (ICELTOT, NL, NU, INL, IAL, INU, IAU,&
NCOLORtot, COLORindex, NEWtoOLD, OLDtoNEW)
endif

&

if (NCOLORtot.eq.0) thencall CM (ICELTOT, NL, NU, INL, IAL, INU, IAU,&
NCOLORtot, COLORindex, NEWtoOLD, OLDtoNEW)
endif

&

if (NCOLORtot.eq.-1) thencall RCM (ICELTOT, NL, NU, INL, IAL, INU, IAU,&
NCOLORtot, COLORindex, NEWtoOLD, OLDtoNEW)
endif

&

if (NCOLORtot.lt.-1) thencall CMRCM (ICELTOT, NL, NU, INL, IAL, INU, IAU,&
NCOLORtot, COLORindex, NEWtoOLD, OLDtoNEW)
endif

&

write (*,'(//a,i8,// )') '### FINAL COLOR NUMBER', NCOLORtot

OMP-332

poi_gen（4/9）
allocate (SMPindex(0:PEsmpTOT*NCOLORtot))
SMPindex= 0
do ic= 1, NCOLORtot
nn1= COLORindex(ic) - COLORindex(ic-1)
num= nn1 / PEsmpTOT
nr = nn1 - PEsmpTOT*num
do ip= 1, PEsmpTOT
if (ip.le.nr) then
SMPindex((ic-1)*PEsmpTOT+ip)= num + 1
else
SMPindex((ic-1)*PEsmpTOT+ip)= num
endif
enddo
enddo
do ic= 1, NCOLORtot
do ip= 1, PEsmpTOT
j1= (ic-1)*PEsmpTOT + ip
j0= j1 - 1
SMPindex(j1)= SMPindex(j0) + SMPindex(j1)
enddo
enddo
allocate (SMPindexG(0:PEsmpTOT))
SMPindexG= 0
nn= ICELTOT / PEsmpTOT
nr= ICELTOT - nn*PEsmpTOT
do ip= 1, PEsmpTOT
SMPindexG(ip)= nn
if (ip.le.nr) SMPindexG(ip)= nn + 1
enddo
do ip= 1, PEsmpTOT
SMPindexG(ip)= SMPindexG(ip-1) + SMPindexG(ip)
enddo
!C===

各色内の要素数：COLORindex(ic)-COLORindex(ic-1)
同じ色内の要素は依存性が無いため，並列に計算可能 ⇒ OpenMP適用
これを更に「PEsmpTOT」で割って「SMPindex」に割り当てる。
前処理で使用

do ic= 1, NCOLORtot!$omp parallel do …
do ip= 1, PEsmpTOT
ip1= (ic-1)*PEsmpTOT+ip
do i= SMPindex(ip1-1)+1, SMPindex(ip1)
(…)
enddo
enddo
!omp end parallel do
enddoOMP-333

SMPindex：
前処理向け

do ic= 1, NCOLORtot!$omp parallel do …
do ip= 1, PEsmpTOT
ip1= (ic-1)*PEsmpTOT+ip
do i= SMPindex(ip1-1)+1, SMPindex(ip1)
(…)
enddo
enddo
!omp end parallel do
enddo

Initial Vector

Coloring(5 colors)+Orderingcolor=1

color=1

color=2

color=2

color=3

color=4

color=3

color=4

color=5

color=512345678 12345678 12345678 12345678 12345678• 5色，8スレッドの例• 同じ「色」に属する要素は独立⇒並列計算可能• 色の順番に並び替え

OMP-334

poi_gen（4/9）
allocate (SMPindex(0:PEsmpTOT*NCOLORtot))
SMPindex= 0
do ic= 1, NCOLORtot
nn1= COLORindex(ic) - COLORindex(ic-1)
num= nn1 / PEsmpTOT
nr = nn1 - PEsmpTOT*num
do ip= 1, PEsmpTOT
if (ip.le.nr) then
SMPindex((ic-1)*PEsmpTOT+ip)= num + 1
else
SMPindex((ic-1)*PEsmpTOT+ip)= num
endif
enddo
enddo
do ic= 1, NCOLORtot
do ip= 1, PEsmpTOT
j1= (ic-1)*PEsmpTOT + ip
j0= j1 - 1
SMPindex(j1)= SMPindex(j0) + SMPindex(j1)
enddo
enddo
allocate (SMPindexG(0:PEsmpTOT))
SMPindexG= 0
nn= ICELTOT / PEsmpTOT
nr= ICELTOT - nn*PEsmpTOT
do ip= 1, PEsmpTOT
SMPindexG(ip)= nn
if (ip.le.nr) SMPindexG(ip)= nn + 1
enddo
do ip= 1, PEsmpTOT
SMPindexG(ip)= SMPindexG(ip-1) + SMPindexG(ip)
enddo
!C===

!$omp parallel do …
do ip= 1, PEsmpTOT
do i= SMPindexG(ip-1)+1, SMPindexG(ip)
(…)
enddo
enddo
!$omp end parallel do
全要素数を「PEsmpTOT」で割って「SMPindexG」に割り当てる。
内積，行列ベクトル積，DAXPYで使用
これを使用すれば，実は，「poi_gen(2/9」の部分も並列化可能「poi_gen(5/9」以降では実際に使用OMP-335

SMPindexG!$omp parallel do …
do ip= 1, PEsmpTOT
do i= SMPindexG(ip-1)+1, SMPindexG(ip)
(…)
enddo
enddo
!$omp end parallel do

ip=1

ip=1

ip=2

ip=2

ip=3

ip=3

ip=4

ip=5

ip=4

ip=5

ip=6

ip=6

ip=7

ip=7

各スレッドで独立に計算：行列ベクトル積，内積，DAXPY等

ip=8

ip=8

OMP-3

36!C!C-- 1D arraynn = ICELTOT
allocate (indexL(0:nn), indexU(0:nn))
indexL= 0
indexU= 0
do icel= 1, ICELTOT
indexL(icel)= INL(icel)
indexU(icel)= INU(icel)
enddo
do icel= 1, ICELTOT
indexL(icel)= indexL(icel) + indexL(icel-1)
indexU(icel)= indexU(icel) + indexU(icel-1)
enddo
NPL= indexL(ICELTOT)
NPU= indexU(ICELTOT)
allocate (itemL(NPL), AL(NPL))
allocate (itemU(NPU), AU(NPU))!C===

itemL= 0
itemU= 0
AL= 0.d0
AU= 0.d0

配列の宣言

poi_gen（5/9）
これ以降は新しい
番号付けを使用OMP-337!C!C +-----------------------------------+!C | INTERIOR & NEUMANN BOUNDARY CELLs |!C +-----------------------------------+!C===!$omp parallel do private(ip,icel,ic0,icN1,icN2,icN3,icN4,icN5,icN6) &!$omp&
private (VOL0,coef,j,ii,jj,kk)
do ip = 1, PEsmpTOT
do icel= SMPindexG(ip-1)+1, SMPindexG(ip)ic0 = NEWtoOLD(icel)
icN1= NEIBcell(ic0,1)
icN2= NEIBcell(ic0,2)
icN3= NEIBcell(ic0,3)
icN4= NEIBcell(ic0,4)
icN5= NEIBcell(ic0,5)
icN6= NEIBcell(ic0,6)

icel： 新しい番号
ic0: 古い番号

VOL0= VOLCEL (ic0)
if (icN5.ne.0) then
icN5= OLDtoNEW(icN5)
coef= RDZ * ZAREA
D(icel)= D(icel) - coef
if (icN5.lt.icel) thendo j= 1, INL(icel)
if (IAL(j,icel).eq.icN5) thenitemL(j+indexL(icel-1))= icN5
AL(j+indexL(icel-1))= coef
exit
endif
enddo
else
do j= 1, INU(icel)
if (IAU(j,icel).eq.icN5) thenitemU(j+indexU(icel-1))= icN5
AU(j+indexU(icel-1))= coef
exit
endif
enddo
endif
endif

poi_gen（6/9）
新しい番号付けを使用
neib (icel ,1)  icel
xneib (icel , 2)  icel
xneib (icel ,3)  icelyneib (icel , 4)  icelyneib (icel ,5)  icel
zneib (icel , 6 )  icelzyz 
yz 
zx 
zx 
xy 
xy  f icel xyz

OMP-338

係数の計算：並列に実施可能
SMPindexG を使用
private宣言に注意!C!C +-----------------------------------+!C | INTERIOR & NEUMANN BOUNDARY CELLs |!C +-----------------------------------+!C===!$omp parallel do private (ip,icel,ic0,icN1,icN2,icN3,icN4,icN5,icN6)!$omp&
private (VOL0,coef,j,ii,jj,kk)
do ip = 1, PEsmpTOT
do icel= SMPindexG(ip-1)+1, SMPindexG(ip)ic0 = NEWtoOLD(icel)
icN1= NEIBcell(ic0,1)
icN2= NEIBcell(ic0,2)
icN3= NEIBcell(ic0,3)
icN4= NEIBcell(ic0,4)
icN5= NEIBcell(ic0,5)
icN6= NEIBcell(ic0,6)
VOL0= VOLCEL (ic0)…

&OMP-339!C!C +-----------------------------------+!C | INTERIOR & NEUMANN BOUNDARY CELLs |!C +-----------------------------------+!C===!$omp parallel do private(ip,icel,ic0,icN1,icN2,icN3,icN4,icN5,icN6) &!$omp&
private (VOL0,coef,j,ii,jj,kk)
do ip = 1, PEsmpTOT
do icel= SMPindexG(ip-1)+1, SMPindexG(ip)ic0 = NEWtoOLD(icel)
icN1= NEIBcell(ic0,1)
icN2= NEIBcell(ic0,2)
icN3= NEIBcell(ic0,3)
icN4= NEIBcell(ic0,4)
icN5= NEIBcell(ic0,5)
icN6= NEIBcell(ic0,6)

icel： 新しい番号
ic0: 古い番号

VOL0= VOLCEL (ic0)
if (icN5.ne.0) then
icN5= OLDtoNEW(icN5)
coef= RDZ * ZAREA
D(icel)= D(icel) - coef
if (icN5.lt.icel) thendo j= 1, INL(icel)
if (IAL(j,icel).eq.icN5) thenitemL(j+indexL(icel-1))= icN5
AL(j+indexL(icel-1))= coef
exit
endif
enddo
else
do j= 1, INU(icel)
if (IAU(j,icel).eq.icN5) thenitemU(j+indexU(icel-1))= icN5
AU(j+indexU(icel-1))= coef
exit
endif
enddo
endif
endif

poi_gen（6/9）
新しい番号付けを使用
neib (icel ,1)  icel
xneib (icel , 2)  icel
xneib (icel ,3)  icelyneib (icel , 4)  icelyneib (icel ,5)  icel
zneib (icel , 6 )  icelzyz 
yz 
zx 
zx 
xy 
xy  f icel xyz

OMP-340!C!C +-----------------------------------+!C | INTERIOR & NEUMANN BOUNDARY CELLs |!C +-----------------------------------+!C===!$omp parallel do private(ip,icel,ic0,icN1,icN2,icN3,icN4,icN5,icN6) &!$omp&
private (VOL0,coef,j,ii,jj,kk)
do ip = 1, PEsmpTOT
do icel= SMPindexG(ip-1)+1, SMPindexG(ip)ic0 = NEWtoOLD(icel)
icN1= NEIBcell(ic0,1)
icN2= NEIBcell(ic0,2)
icN3= NEIBcell(ic0,3)
icN4= NEIBcell(ic0,4)
icN5= NEIBcell(ic0,5)
icN6= NEIBcell(ic0,6)
VOL0= VOLCEL(ic0)

if (icN5.ne.0) then
icN5= OLDtoNEW(icN5)
coef= RDZ * ZAREA
D(icel)= D(icel) - coef
if (icN5.lt.icel) thendo j= 1, INL(icel)
if (IAL(j,icel).eq.icN5) thenitemL(j+indexL(icel-1))= icN5
AL(j+indexL(icel-1))= coef
exit
endif
enddo
else
do j= 1, INU(icel)
if (IAU(j,icel).eq.icN5) thenitemU(j+indexU(icel-1))= icN5
AU(j+indexU(icel-1))= coef
exit
endif
enddo
endif
endif

poi_gen（6/9）
新しい番号付けを使用
neib (icel ,1)  icel
xneib (icel , 2)  icel
xneib (icel ,3)  icelyneib (icel , 4)  icelyneib (icel ,5)  icel
zneib (icel , 6 )  icelzyz 
yz 
zx 
zx 
xy 
xy  f icel xyz

OMP-341!C!C +-----------------------------------+!C | INTERIOR & NEUMANN BOUNDARY CELLs |!C +-----------------------------------+!C===!$omp parallel do private(ip,icel,ic0,icN1,icN2,icN3,icN4,icN5,icN6) &!$omp&
private (VOL0,coef,j,ii,jj,kk)
do ip = 1, PEsmpTOT
do icel= SMPindexG(ip-1)+1, SMPindexG(ip)ic0 = NEWtoOLD(icel)(ic0)

if (icN5.ne.0) then
icN5= OLDtoNEW(icN5)
coef= RDZ * ZAREA
D(icel)= D(icel) - coef

新しい番号付けを使用
neib (icel ,1)  icel
x

icN1= NEIBcell(ic0,1)
icN2= NEIBcell(ic0,2)
icN3= NEIBcell(ic0,3)
icN4= NEIBcell(ic0,4)
icN5= NEIBcell(ic0,5)
icN6= NEIBcell(ic0,6)
VOL0= VOLCEL

poi_gen（6/9）

neib (icel , 2)  icel
1
z
ZAREA  xyRDZ if (icN5.lt.icel) thendo j= 1, INL(icel)
if (IAL(j,icel).eq.icN5) thenitemL(j+indexL(icel-1))= icN5
AL(j+indexL(icel-1))= coef
exit
endif
icN5がicelより小さ
enddo
ければ下三角成分else
do j= 1, INU(icel)
if (IAU(j,icel).eq.icN5) thenitemU(j+indexU(icel-1))= icN5
AU(j+indexU(icel-1))= coef
exit
endif
enddo
endif
endif

x

neib (icel ,3)  icelyneib (icel , 4)  icelyneib (icel ,5)  icel
zneib (icel , 6 )  icelzyz 
yz 
zx 
zx 
xy 
xy  f icel xyz

OMP-342!C!C +-----------------------------------+!C | INTERIOR & NEUMANN BOUNDARY CELLs |!C +-----------------------------------+!C===!$omp parallel do private(ip,icel,ic0,icN1,icN2,icN3,icN4,icN5,icN6) &!$omp&
private (VOL0,coef,j,ii,jj,kk)
do ip = 1, PEsmpTOT
do icel= SMPindexG(ip-1)+1, SMPindexG(ip)ic0 = NEWtoOLD(icel)(ic0)

if (icN5.ne.0) then
icN5= OLDtoNEW(icN5)
coef= RDZ * ZAREA
D(icel)= D(icel) - coef

新しい番号付けを使用
neib (icel ,1)  icel
x

icN1= NEIBcell(ic0,1)
icN2= NEIBcell(ic0,2)
icN3= NEIBcell(ic0,3)
icN4= NEIBcell(ic0,4)
icN5= NEIBcell(ic0,5)
icN6= NEIBcell(ic0,6)
VOL0= VOLCEL

poi_gen（6/9）

neib (icel , 2)  icel
1
z
ZAREA  xyRDZ if (icN5.lt.icel) thendo j= 1, INL(icel)
if (IAL(j,icel).eq.icN5) thenitemL(j+indexL(icel-1))= icN5
AL(j+indexL(icel-1))= coef
exit
endif
icN5がicelより大きenddo
ければ上三角成分else
do j= 1, INU(icel)
if (IAU(j,icel).eq.icN5) thenitemU(j+indexU(icel-1))= icN5
AU(j+indexU(icel-1))= coef
exit
endif
enddo
endif
endif

x

neib (icel ,3)  icelyneib (icel , 4)  icelyneib (icel ,5)  icel
zneib (icel , 6 )  icelzyz 
yz 
zx 
zx 
xy 
xy  f icel xyz

OMP-343

if (icN3.ne.0) then
icN3= OLDtoNEW(icN3)
coef= RDY * YAREA
D(icel)= D(icel) - coef
if (icN3.lt.icel) thendo j= 1, INL(icel)
if (IAL(j,icel).eq.icN3) thenitemL(j+indexL(icel-1))= icN3
AL(j+indexL(icel-1))= coef
exit
endif
enddo
else
do j= 1, INU(icel)
if (IAU(j,icel).eq.icN3) thenitemU(j+indexU(icel-1))= icN3
AU(j+indexU(icel-1))= coef
exit
endif
enddo
endif
endif
if (icN1.ne.0) then
icN1= OLDtoNEW(icN1)
coef= RDX * XAREA
D(icel)= D(icel) - coef
if (icN1.lt.icel) thendo j= 1, INL(icel)
if (IAL(j,icel).eq.icN1) thenitemL(j+indexL(icel-1))= icN1
AL(j+indexL(icel-1))= coef
exit
endif
enddo
else
do j= 1, INU(icel)
if (IAU(j,icel).eq.icN1) thenitemU(j+indexU(icel-1))= icN1
AU(j+indexU(icel-1))= coef
exit
endif
enddo
endif
endif

poi_gen（7/9）
neib (icel ,1)  icel
xneib (icel , 2)  icel
xneib (icel ,3)  icelyneib (icel , 4)  icelyneib (icel ,5)  icel
zneib (icel , 6 )  icelzyz 
yz 
zx 
zx 
xy 
xy  f icel xyz

OMP-344

if (icN2.ne.0) then
icN2= OLDtoNEW(icN2)
coef= RDX * XAREA
D(icel)= D(icel) - coef
if (icN2.lt.icel) thendo j= 1, INL(icel)
if (IAL(j,icel).eq.icN2) thenitemL(j+indexL(icel-1))= icN2
AL(j+indexL(icel-1))= coef
exit
endif
enddo
else
do j= 1, INU(icel)
if (IAU(j,icel).eq.icN2) thenitemU(j+indexU(icel-1))= icN2
AU(j+indexU(icel-1))= coef
exit
endif
enddo
endif
endif
if (icN4.ne.0) then
icN4= OLDtoNEW(icN4)
coef= RDY * YAREA
D(icel)= D(icel) - coef
if (icN4.lt.icel) thendo j= 1, INL(icel)
if (IAL(j,icel).eq.icN4) thenitemL(j+indexL(icel-1))= icN4
AL(j+indexL(icel-1))= coef
exit
endif
enddo
else
do j= 1, INU(icel)
if (IAU(j,icel).eq.icN4) thenitemU(j+indexU(icel-1))= icN4
AU(j+indexU(icel-1))= coef
exit
endif
enddo
endif
endif

poi_gen（8/9）
neib (icel ,1)  icel
xneib (icel , 2)  icel
xneib (icel ,3)  icelyneib (icel , 4)  icelyneib (icel ,5)  icel
zneib (icel , 6 )  icelzyz 
yz 
zx 
zx 
xy 
xy  f icel xyz

OMP-345!$omp parallel do private(ip,icel,ic0,icN1,icN2,icN3,icN4,icN5,icN6) &!$omp&
private (VOL0,coef,j,ii,jj,kk)…
if (icN6.ne.0) then
icN6= OLDtoNEW(icN6)
coef= RDZ * ZAREA
D(icel)= D(icel) - coef
if (icN6.lt.icel) thendo j= 1, INL(icel)
if (IAL(j,icel).eq.icN6) thenitemL(j+indexL(icel-1))= icN6
AL(j+indexL(icel-1))= coef
exit
endif
enddo
else
do j= 1, INU(icel)
if (IAU(j,icel).eq.icN6) thenitemU(j+indexU(icel-1))= icN6
AU(j+indexU(icel-1))= coef
exit
endif
enddo
endif
endif
ii= XYZ(ic0,1)
もとの座標に従って
jj= XYZ(ic0,2)
BFORCE計算
kk= XYZ(ic0,3)
BFORCE(icel)= -dfloat(ii+jj+kk) * VOL0
enddo
enddo
!$omp end parallel do!C===ii,jj,kk,VOL0はprivate

poi_gen（9/9）
係数の計算（境界面以外）

neib (icel ,1)  icel
xneib (icel , 2)  icel
xneib (icel ,3)  icelyneib (icel , 4)  icelyneib (icel ,5)  icel
zneib (icel , 6 )  icelzyz 
yz 
zx 
zx 
xy 
xy  f icel xyz

OMP-3

46

プログラムの構成
program MAIN
use STRUCT
use PCGuse solver_ICCG_mc
implicit REAL*8 (A-H,O-Z)real(kind=8), dimension(:), allocatable :: WK
call INPUT
call POINTER_INIT
call BOUNDARY_CELL
call CELL_METRICS
call POI_GEN
PHI=
&
&
&

0.d0

call solve_ICCG_mc( ICELTOT, NPL, NPU, indexL, itemL, indexU, itemU, D,
BFORCE, PHI, AL, AU, NCOLORtot, PEsmpTOT,
SMPindex, SMPindexG, EPSICCG, ITR, IER)&&&

この時点で，係数，右辺ベクトル
ともに，「新しい」番号にしたがって
計算，記憶されている。

OMP-347

solve_ICCG_mc（1/6）!C***!C*** module solver_ICCG_mc!C***!
module solver_ICCG_mc
contains!C!C*** solve_ICCG!C
subroutine solve_ICCG_mc&( N, NPL, NPU, indexL, itemL, indexU, itemU, D, B, X,&
AL, AU, NCOLORtot, PEsmpTOT, SMPindex, SMPindexG,&
EPS, ITR, IER)
implicit REAL*8 (A-H,O-Z)
integer :: N, NL, NU, NCOLORtot, PEsmpTOT
real(kind=8), dimension(N)real(kind=8), dimension(N)real(kind=8), dimension(N)

:: D
:: B
:: X

real(kind=8), dimension(NPL) :: ALreal(kind=8), dimension(NPU) :: AU
integer, dimension(0:N):: indexL, indexU
integer, dimension(NPL):: itemL
integer, dimension(NPU):: itemU
integer, dimension(0:NCOLORtot*PEsmpTOT):: SMPindex
integer, dimension(0:PEsmpTOT):: SMPindexG
real(kind=8), dimension(:,:), allocatable :: W
integer, parameter :: R= 1
integer, parameter :: Z= 2
integer, parameter :: Q= 2
integer, parameter :: P= 3
integer, parameter :: DD= 4&&&

OMP-348

solve_ICCG_mc（2/6）!C!C +------+!C | INIT |!C +------+!C===
allocate (W(N,4))!$omp parallel do private(ip,i)do ip= 1, PEsmpTOT
do i = SMPindexG(ip-1)+1, SMPindexG(ip)
X(i) = 0.d0
W(i,2)= 0.0D0
W(i,3)= 0.0D0
W(i,4)= 0.0D0
enddo
enddo
!$omp end parallel dodo ic= 1, NCOLORtot!$omp parallel do private(ip,ip1,i,VAL,k)do ip= 1, PEsmpTOT
ip1= (ic-1)*PEsmpTOT + ip
do i= SMPindex(ip1-1)+1, SMPindex(ip1)
VAL= D(i)do k= indexL(i-1)+1, indexL(i)VAL= VAL - (AL(k)**2) * W(itemL(k),DD)
enddo
W(i,DD)= 1.d0/VAL
enddo
enddo
!$omp end parallel do
enddo

不完全修正
コレスキー分解OMP-349

不完全修正コレスキー分解
1



2
1
d i   aii 
aik  d k   lii
k 1


i 1



W(i,DD):
D(i):

di
aii

itemL(j):AL(j):kaik

do i= 1, N
VAL= D(i)do k= indexL(i-1)+1, indexL(i)VAL= VAL - (AL(k)**2) * W(itemL(k),DD)
enddo
W(i,DD)= 1.d0/VAL
enddoOMP-350

不完全修正コレスキー分解：並列版
1



2
1
d i   aii 
aik  d k   lii
k 1


i 1



W(i,DD):
D(i):

di
aii

itemL(j):AL(j):kaik

do ic= 1, NCOLORtot!$omp parallel do private(ip,ip1,i,VAL,k)do ip= 1, PEsmpTOT
ip1= (ic-1)*PEsmpTOT + ip
do i= SMPindex(ip1-1)+1, SMPindex(ip1)
VAL= D(i)do k= indexL(i-1)+1, indexL(i)VAL= VAL - (AL(k)**2) * W(itemL(k),DD)
enddo
W(i,DD)= 1.d0/VAL
enddo
enddo
!$omp end parallel do
enddo

privateに注意。

OMP-351

solve_ICCG_mc（3/6）!C +-----------------------+!C | {r0}= {b} - [A]{xini} |!C +-----------------------+!C===!$omp parallel do private(ip,i,VAL,k)do ip= 1, PEsmpTOT
do i = SMPindexG(ip-1)+1, SMPindexG(ip)
VAL= D(i)*X(i)do k= indexL(i-1)+1, indexL(i)VAL= VAL + AL(k)*X(itemL(k))
enddo
do k= indexU(i-1)+1, indexU(i)VAL= VAL + AU(k)*X(itemU(k))
enddo
W(i,R)= B(i) - VAL
enddo
enddo
!$omp end parallel do
BNRM2= 0.0D0!$omp parallel do private(ip,i) reduction(+:BNRM2)
do ip= 1, PEsmpTOT
do i = SMPindexG(ip-1)+1, SMPindexG(ip)
BNRM2 = BNRM2 + B(i) **2
enddo
enddo
!$omp end parallel do!C===

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(i)i = i-1/p(i)q(i)x(i)= x(i-1) + ip(i)r(i)= r(i-1) - iq(i)check convergence |r|endOMP-352

行列ベクトル積
依存性が無い⇒独立に計算可能⇒SMPindexG使用!$omp parallel do private(ip,i,VAL,k)do ip= 1, PEsmpTOT
do i = SMPindexG(ip-1)+1, SMPindexG(ip)
VAL= D(i)*X(i)do k= indexL(i-1)+1, indexL(i)VAL= VAL + AL(k)*X(itemL(k))
enddo
do k= indexU(i-1)+1, indexU(i)VAL= VAL + AU(k)*X(itemU(k))
enddo
W(i,R)= B(i) - VAL
enddo
enddo
!$omp end parallel doOMP-353

solve_ICCG_mc（3/6）!C +-----------------------+!C | {r0}= {b} - [A]{xini} |!C +-----------------------+!C===!$omp parallel do private(ip,i,VAL,k)do ip= 1, PEsmpTOT
do i = SMPindexG(ip-1)+1, SMPindexG(ip)
VAL= D(i)*X(i)do k= indexL(i-1)+1, indexL(i)VAL= VAL + AL(k)*X(itemL(k))
enddo
do k= indexU(i-1)+1, indexU(i)VAL= VAL + AU(k)*X(itemU(k))
enddo
W(i,R)= B(i) - VAL
enddo
enddo
!$omp end parallel do
BNRM2= 0.0D0!$omp parallel do private(ip,i) reduction(+:BNRM2)
do ip= 1, PEsmpTOT
do i = SMPindexG(ip-1)+1, SMPindexG(ip)
BNRM2 = BNRM2 + B(i) **2
enddo
enddo
!$omp end parallel do!C===

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(i)i = i-1/p(i)q(i)x(i)= x(i-1) + ip(i)r(i)= r(i-1) - iq(i)check convergence |r|endOMP-354

内積：SMPindexG使用，reduction
BNRM2= 0.0D0!$omp parallel do private(ip,i) reduction(+:BNRM2)
do ip= 1, PEsmpTOT
do i = SMPindexG(ip-1)+1, SMPindexG(ip)
BNRM2 = BNRM2 + B(i) **2
enddo
enddo
!$omp end parallel doOMP-355

ITR= N
do L= 1, ITR!C!C +----------------+!C | {z}= [Minv]{r} |!C +----------------+!C===!$omp parallel do private(ip,i)do ip= 1, PEsmpTOT
do i = SMPindexG(ip-1)+1, SMPindexG(ip)
W(i,Z)= W(i,R)
enddo
enddo
!$omp end parallel dodo ic= 1, NCOLORtot!$omp parallel do private(ip,ip1,i,WVAL,j)do ip= 1, PEsmpTOT
ip1= (ic-1)*PEsmpTOT + ip
do i= SMPindex(ip1-1)+1, SMPindex(ip1)
WVAL= W(i,Z)
do j= 1, INL(i)WVAL= WVAL - AL(j,i) * W(IAL(j,i),Z)
enddo
W(i,Z)= WVAL * W(i,DD)
enddo
enddo
!$omp end parallel do
enddo
do ic= NCOLORtot, 1, -1!$omp parallel do private(ip,ip1,i,SW,j)do ip= 1, PEsmpTOT
ip1= (ic-1)*PEsmpTOT + ip
do i= SMPindex(ip1-1)+1, SMPindex(ip1)
SW = 0.0d0
do j= 1, INU(i)SW= SW + AU(j,i) * W(IAU(j,i),Z)
enddo
W(i,Z)= W(i,Z) - W(i,DD) * SW
enddo
enddo
!$omp end parallel do
enddo
!C===

solve_ICCG_mc（4/6）
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(i)i = i-1/p(i)q(i)x(i)= x(i-1) + ip(i)r(i)= r(i-1) - iq(i)check convergence |r|endOMP-356

ITR= N
do L= 1, ITR!C!C +----------------+ ここでは「SMPindex」を使う!C | {z}= [Minv]{r} |!C +----------------+!C===!$omp parallel do private(ip,i)do ip= 1, PEsmpTOT
do i = SMPindexG(ip-1)+1, SMPindexG(ip)
W(i,Z)= W(i,R)
enddo
enddo
!$omp end parallel dodo ic= 1, NCOLORtot!$omp parallel do private(ip,ip1,i,WVAL,k)do ip= 1, PEsmpTOT
ip1= (ic-1)*PEsmpTOT + ip
do i= SMPindex(ip1-1)+1, SMPindex(ip1)
WVAL= W(i,Z)
do k= indexL(i-1)+1, indexL(i)WVAL= WVAL - AL(k) * W(itemL(k),Z)
enddo
W(i,Z)= WVAL * W(i,DD)
enddo
enddo
!$omp end parallel do
enddo
do ic= NCOLORtot, 1, -1!$omp parallel do private(ip,ip1,i,SW,k)do ip= 1, PEsmpTOT
ip1= (ic-1)*PEsmpTOT + ip
do i= SMPindex(ip1-1)+1, SMPindex(ip1)
SW = 0.0d0
do k= indexU(i-1)+1, indexU(i)SW= SW + AU(k) * W(itemU(k),Z)
enddo
W(i,Z)= W(i,Z) - W(i,DD) * SW
enddo
enddo
!$omp end parallel do
enddo
!C===

solve_ICCG_mc（4/6）
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(i)i = i-1/p(i)q(i)x(i)= x(i-1) + ip(i)r(i)= r(i-1) - iq(i)check convergence |r|endOMP-357

ITR= N
do L= 1, ITR!C!C +----------------+ ここでは「SMPindex」を使う!C | {z}= [Minv]{r} |!C +----------------+!C===!$omp parallel do private(ip,i)do ip= 1, PEsmpTOT
do i = SMPindexG(ip-1)+1, SMPindexG(ip)
W(i,Z)= W(i,R)
enddo
enddo
!$omp end parallel dodo ic= 1, NCOLORtot!$omp parallel do private(ip,ip1,i,WVAL,k)do ip= 1, PEsmpTOT
ip1= (ic-1)*PEsmpTOT + ip
do i= SMPindex(ip1-1)+1, SMPindex(ip1)
WVAL= W(i,Z)
do k= indexL(i-1)+1, indexL(i)WVAL= WVAL - AL(k) * W(itemL(k),Z)
enddo
W(i,Z)= WVAL * W(i,DD)
enddo
enddo
!$omp end parallel do
enddo
do ic= NCOLORtot, 1, -1!$omp parallel do private(ip,ip1,i,SW,k)do ip= 1, PEsmpTOT
ip1= (ic-1)*PEsmpTOT + ip
do i= SMPindex(ip1-1)+1, SMPindex(ip1)
SW = 0.0d0
do k= indexU(i-1)+1, indexU(i)SW= SW + AU(k) * W(itemU(k),Z)
enddo
W(i,Z)= W(i,Z) - W(i,DD) * SW
enddo
enddo
!$omp end parallel do
enddo
!C===

solve_ICCG_mc（4/6）
M z  LDL z  r
TL z  r
前進代入
Forward SubstitutionDL z  z
T

後退代入
Backward SubstitutionOMP-358

前進代入：SMPindex使用
do ic= 1, NCOLORtot!$omp parallel do private(ip,ip1,i,WVAL,k)do ip= 1, PEsmpTOT
ip1= (ic-1)*PEsmpTOT + ip
do i= SMPindex(ip1-1)+1, SMPindex(ip1)
WVAL= W(i,Z)
do k= indexL(i-1)+1, indexL(i)WVAL= WVAL - AL(k) * W(indexL(k),Z)
enddo
W(i,Z)= WVAL * W(i,DD)
enddo
enddo
!$omp end parallel do
enddoOMP-3!C +-----------------------------+!C | {p} = {z} if
ITER=1|!C | BETA= RHO / RHO1 otherwise |!C +-----------------------------+!C===
if ( L.eq.1 ) then!$omp parallel do private(ip,i)do ip= 1, PEsmpTOT
do i = SMPindexG(ip-1)+1, SMPindexG(ip)
W(i,P)= W(i,Z)
enddo
enddo
!$omp end parallel do
else
BETA= RHO / RHO1!$omp parallel do private(ip,i)do ip= 1, PEsmpTOT
do i = SMPindexG(ip-1)+1, SMPindexG(ip)
W(i,P)= W(i,Z) + BETA*W(i,P)
enddo
enddo
!$omp end parallel do
endif
!C===
!C +-------------+!C | {q}= [A]{p} |!C +-------------+!C===!$omp parallel do private(ip,i,VAL,k)do ip= 1, PEsmpTOT
do i = SMPindexG(ip-1)+1, SMPindexG(ip)
VAL= D(i)*W(i,P)
do k= indexL(i-1)+1, indexL(i)VAL= VAL + AL(k)*W(itemL(k),P)
enddo
do k= indexU(i-1)+1, indexU(i)VAL= VAL + AU(k)*W(itemU(k),P)
enddo
W(i,Q)= VAL
enddo
enddo
!$omp end parallel do
!C===

59

solve_ICCG_mc（5/6）
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(i)i = i-1/p(i)q(i)x(i)= x(i-1) + ip(i)r(i)= r(i-1) - iq(i)check convergence |r|endOMP-3!C +-----------------------------+!C | {p} = {z} if
ITER=1|!C | BETA= RHO / RHO1 otherwise |!C +-----------------------------+!C===
if ( L.eq.1 ) then!$omp parallel do private(ip,i)do ip= 1, PEsmpTOT
do i = SMPindexG(ip-1)+1, SMPindexG(ip)
W(i,P)= W(i,Z)
enddo
enddo
!$omp end parallel do
else
BETA= RHO / RHO1!$omp parallel do private(ip,i)do ip= 1, PEsmpTOT
do i = SMPindexG(ip-1)+1, SMPindexG(ip)
W(i,P)= W(i,Z) + BETA*W(i,P)
enddo
enddo
!$omp end parallel do
endif
!C===
!C +-------------+!C | {q}= [A]{p} |!C +-------------+!C===!$omp parallel do private(ip,i,VAL,k)do ip= 1, PEsmpTOT
do i = SMPindexG(ip-1)+1, SMPindexG(ip)
VAL= D(i)*W(i,P)
do k= indexL(i-1)+1, indexL(i)VAL= VAL + AL(k)*W(itemL(k),P)
enddo
do k= indexU(i-1)+1, indexU(i)VAL= VAL + AU(k)*W(itemU(k),P)
enddo
W(i,Q)= VAL
enddo
enddo
!$omp end parallel do
!C===

60

solve_ICCG_mc（5/6）
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(i)i = i-1/p(i)q(i)x(i)= x(i-1) + ip(i)r(i)= r(i-1) - iq(i)check convergence |r|endOMP-361!C!C +---------------------+!C | ALPHA= RHO / {p}{q} |!C +---------------------+!C===
C1= 0.d0!$omp parallel do private(ip,i) reduction(+:C1)
do ip= 1, PEsmpTOT
do i = SMPindexG(ip-1)+1, SMPindexG(ip)C1= C1 + W(i,P)*W(i,Q)
enddo
enddo
!$omp end parallel do!C===

ALPHA= RHO / C1

!C
!C +----------------------+!C | {x}= {x} + ALPHA*{p} |!C | {r}= {r} - ALPHA*{q} |!C +----------------------+!C===!$omp parallel do private(ip,i)do ip= 1, PEsmpTOT
do i = SMPindexG(ip-1)+1, SMPindexG(ip)
X(i) = X(i)+ ALPHA * W(i,P)
W(i,R)= W(i,R) - ALPHA * W(i,Q)
enddo
enddo
!$omp end parallel do
DNRM2= 0.d0!$omp parallel do private(ip,i) reduction(+:DNRM2)
do ip= 1, PEsmpTOT
do i = SMPindexG(ip-1)+1, SMPindexG(ip)
DNRM2= DNRM2 + W(i,R)**2
enddo
enddo
!$omp end parallel do!C===solve_ICCG_mc（6/6）
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(i)i = i-1/p(i)q(i)x(i)= x(i-1) + ip(i)r(i)= r(i-1) - iq(i)check convergence |r|endOMP-362!C!C +---------------------+!C | ALPHA= RHO / {p}{q} |!C +---------------------+!C===
C1= 0.d0!$omp parallel do private(ip,i) reduction(+:C1)
do ip= 1, PEsmpTOT
do i = SMPindexG(ip-1)+1, SMPindexG(ip)C1= C1 + W(i,P)*W(i,Q)
enddo
enddo
!$omp end parallel do!C===

ALPHA= RHO / C1

!C
!C +----------------------+!C | {x}= {x} + ALPHA*{p} |!C | {r}= {r} - ALPHA*{q} |!C +----------------------+!C===!$omp parallel do private(ip,i)do ip= 1, PEsmpTOT
do i = SMPindexG(ip-1)+1, SMPindexG(ip)
X(i) = X(i)+ ALPHA * W(i,P)
W(i,R)= W(i,R) - ALPHA * W(i,Q)
enddo
enddo
!$omp end parallel do
DNRM2= 0.d0!$omp parallel do private(ip,i) reduction(+:DNRM2)
do ip= 1, PEsmpTOT
do i = SMPindexG(ip-1)+1, SMPindexG(ip)
DNRM2= DNRM2 + W(i,R)**2
enddo
enddo
!$omp end parallel do!C===solve_ICCG_mc（6/6）
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(i)i = i-1/p(i)q(i)x(i)= x(i-1) + ip(i)r(i)= r(i-1) - iq(i)check convergence |r|endOMP-363!C!C +---------------------+!C | ALPHA= RHO / {p}{q} |!C +---------------------+!C===
C1= 0.d0!$omp parallel do private(ip,i) reduction(+:C1)
do ip= 1, PEsmpTOT
do i = SMPindexG(ip-1)+1, SMPindexG(ip)C1= C1 + W(i,P)*W(i,Q)
enddo
enddo
!$omp end parallel do!C===

ALPHA= RHO / C1

!C
!C +----------------------+!C | {x}= {x} + ALPHA*{p} |!C | {r}= {r} - ALPHA*{q} |!C +----------------------+!C===!$omp parallel do private(ip,i)do ip= 1, PEsmpTOT
do i = SMPindexG(ip-1)+1, SMPindexG(ip)
X(i) = X(i)+ ALPHA * W(i,P)
W(i,R)= W(i,R) - ALPHA * W(i,Q)
enddo
enddo
!$omp end parallel do
DNRM2= 0.d0!$omp parallel do private(ip,i) reduction(+:DNRM2)
do ip= 1, PEsmpTOT
do i = SMPindexG(ip-1)+1, SMPindexG(ip)
DNRM2= DNRM2 + W(i,R)**2
enddo
enddo
!$omp end parallel do!C===solve_ICCG_mc（6/6）
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(i)i = i-1/p(i)q(i)x(i)= x(i-1) + ip(i)r(i)= r(i-1) - iq(i)check convergence |r|endOMP-3• L2-solへのOpenMPの実装• 実行例• 最適化＋演習

64

OMP-365

計算結果• Hitachi SR11000/J2 1ノード（16コア）• 1003要素
MemoryL3Core Core
L1

L1
L3

Memory

L3

L1

Core Core
L1

Memory

L3

L1

Core Core
L1

Memory

L3

L1

Core Core
L1

L1

L2

L2

L2

L2

L2

L2

L2

L2

L1

Core CPU

Memory

L1
L3

L1

Core Core

Memory

L1
L3

L1

Core Core

Memory

L1
L3

L1

Core Core

MemorySR11000，16コアにおける結果，1003（●：MC，△：RCM，-：CM-RCM）
1.E+06

Incompatible Point #ITERATIONS

350

300

250

1.E+04

1.E+02

反復回数

ICP#

200
1.E+00

1.E+01

1.E+02

1.E+03

1.E+00
1.E+00

1.E+04

1.E+01COLOR#

1.E+02

1.E+03

1.E+04COLOR#

1.0E-02

3.00

計算時間／反復計算時間
sec./iteration

2.50
sec.

OMP-3

2.00

8.0E-03

6.0E-03

挙動お
かしい
1.50
1.E+00

1.E+01
1.E+02COLOR#

1.E+03

4.0E-03
1.E+00

1.E+01

1.E+02COLOR#

1.E+0366

OMP-3

FX10，16コアにおける結果，1283（●：MC，△：RCM，-：CM-RCM）
Number of Incompatible Nodes

550

反復回数

500

Iterations

450
400
350
300

MC

RCM

CM-RCM

1.E+07
1.E+06
1.E+05
1.E+04
1.E+03MC

1.E+02RCM
1.E+01ICP#

CM-RCM

1.E+00

250
1

10

100

1

1000

COLOR#

25.00

1000

MC
5.00E-02

sec./iteration

20.00

100

COLOR#

6.00E-02

MC
RCM
CM-RCM

10

sec.

15.00

10.00

5.00

RCM
CM-RCM

4.00E-02
3.00E-02
2.00E-02
1.00E-02

計算時間
0.00

計算時間／反復

0.00E+00
1

10

COLOR#

100

1000

1

10

COLOR#

100

1000

67

OMP-3• L2-solへのOpenMPの実装• 実行例• 最適化＋演習

68

OMP-E69• マルチコア版コードの実行• 更なる最適化• STREAM• プロファイラ，コンパイルリスト分析等

OMP-1

コンパイル・実行>$ cd <$O-L3>/src>$ make>$ ls ../run/L3-sol
L3-sol>$ cd ../run>$ pjsub go1.sh

70OMP-371

プログラムの実行
プログラム，必要なファイル等

L3-solポアソン方程式
ソルバー

INPUT.DAT
制御ファイルtest.inp
結果ファイル（ParaView）

OMP-372

制御データ（INPUT.DAT）
100 100 100
1.00e-00 1.00e-00 1.00e-00
1.0e-08
16-10

NX/NY/NZ
DX/DY/DZEPSICCGPEsmpTOTNCOLORtot

変数名

型

内

容

NX，NY，NZ

整数

DX，DY，DZ

倍精度実数

各要素の3辺の長さ（X，Y，Z）

EPSICCG

倍精度実数

収束判定値

PEsmpTOT

整数

データ分割数NCOLORtot

整数

Ordering手法と色数≧2：MC法（multicolor），色数=0：CM法（Cuthill-Mckee）=-1：RCM法（Reverse Cuthill-Mckee）≦-2：CM-RCM法

各方向の要素数OMP-173

go1.sh#!/bin/sh#PJM -L "node=1"#PJM -L "elapse=00:10:00"#PJM -L "rscgrp=small"#PJM -j#PJM -o “test.lst“
export OMP_NUM_THREADS=16./L3-sol

標準出力ファイル名
スレッド数，通常 =PEsmpTOTOMP-374

計算結果（FX10@東大）：106要素
反復回数：MC（2色）：333回，RCM（298レベル）：224回
CM-RCM（Nc=20）：249回
30.0

20.0
MC=2

MC=2

RCM(298)

16.0

CM-RCM(20)RCM(298)
CM-RCM(20)sec.

Speed-Up

20.0

10.0

12.0

8.0

4.0
0.0

0.0
1

2

4

8

16

1

2

4

thread#16 threads
MC(2):
2.42 sec.
CM-RCM(20): 2.01 sec.

8

16

thread#Memory

L1

L1

L1

L1

L1

L1

L1

L2
L1 L1

L1

L1

L1

L1

L1

L1

L1

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

75

実習（1）
• 色々なケースでやってみよう– 問題サイズ– スレッド数– 色数，色分け法（MC，RCM，CM-RCM）

OMP-E76• マルチコア版コードの実行• 更なる最適化– その1： OpenMP Statement– その2： Sequential Reordering– その3： ELL• STREAM• プロファイラ，コンパイルリスト分析等

OMP-E77

前進代入：現状の並列化（Fortran）
do ic= 1, NCOLORtot!$omp parallel do private(ip,ip1,i,WVAL,k)do ip= 1, PEsmpTOT
ip1= (ic-1)*PEsmpTOT + ip
do i= SMPindex(ip1-1)+1, SMPindex(ip1)
WVAL= W(i,Z)
do k= indexL(i-1)+1, indexL(i)WVAL= WVAL - AL(k) * W(itemL(k),Z)
enddo
W(i,Z)= WVAL * W(i,DD)
enddo
enddo
!$omp end parallel do
enddo• 「!omp parallel」でスレッド（～16）の生成，消滅が発生– 色ごとにこの部分を通る– 多少のオーバーヘッドがある• 色数が増えるとオーバーヘッドが増すOMP-E78

前進代入：Overhead削減（Fortran）!$omp parallel private(ic,ip,ip1,i,WVAL,k)do ic= 1, NCOLORtot!$omp dodo ip= 1, PEsmpTOT
ip1= (ic-1)*PEsmpTOT + ip
do i= SMPindex(ip1-1)+1, SMPindex(ip1)
WVAL= W(i,Z)
do k= indexL(i-1)+1, indexL(i)WVAL= WVAL - AL(k) * W(itemL(k),Z)
enddo
W(i,Z)= WVAL * W(i,DD)
enddo
enddo
endd
!$omp end parallel• このようにすることによって，スレッド生成を前進代入に入る前の一回で済ませることができる• 「!omp do」のループが並列化

OMP-E

79

プログラム類% cd <$O-L3>% ls
run reorder0

src% cd src0% make% cd ../run% ls L3-sol0
L3-sol0% <modify “INPUT.DAT”>% <modify “go0.sh”>% pjsub go0.sh

src0OMP-E80

計算結果：L3-sol0が速い
N=1283NCOLORtot= -20
CM-RCM (20)318 Iterations
NCOLORtot= -1
RCM (382 levels)
287 Iterations

L3-sol

L3-sol05.69 sec.

5.44 sec.

6.54 sec.

6.37 sec.

OMP-E81• マルチコア版コードの実行• 更なる最適化– その1： OpenMP Statement– その2： Sequential Reordering– その3： ELL• STREAM• プロファイラ，コンパイルリスト分析等

OMP-382

現在のオーダリングの問題
• 色付け
– MC
– RCM
– CM-RCM• 同じ色に属する要素は独立：並列計算可能• 「色」の順番に番号付け• 色内の要素を各スレッドに振り分ける• 同じスレッド（すなわち同じコア）に属する要素は連続の番号ではない– 効率の低下OMP-383

SMPindex：
前処理向け

do ic= 1, NCOLORtot!$omp parallel do …
do ip= 1, PEsmpTOT
ip1= (ic-1)*PEsmpTOT+ip
do i= SMPindex(ip1-1)+1, SMPindex(ip1)
(…)
enddo
enddo
!omp end parallel do
enddo

Initial Vector

Coloring(5 colors)+Orderingcolor=1

color=1

color=2

color=2

color=3

color=4

color=3

color=4

color=5

color=512345678 12345678 12345678 12345678 12345678• 5色，8スレッドの例• 同じ「色」に属する要素は独立⇒並列計算可能• 色の順番に並び替え

OMP-384

データ再配置：Sequential Reordering• 同じスレッドで処理するデータをなるべく連続に配置するように更に並び替え– 効率の向上が期待される• 係数行列等のアドレスが連続になる• 局所性が高まる（2ページあと）• 番号の付け替えによって要素の大小関係は変わるが，上三角，下三角の関係は変えない（もとの計算と反復回数は変わらない）– 従って自分より要素番号が大きいのにIAL（下三角）に含まれていたりするOMP-385

データ再配置: Sequential Reordering
各スレッド上でメモリアクセスが連続となるよう更に並び替え
5 colors, 8 threads
Initial Vector
Coloring(5 colors)+Orderingcolor=1

color=1

color=2

color=2

color=3

color=4

color=5

color=3

color=4

color=5

Coalesced

12345678 12345678 12345678 12345678 12345678

Sequential

11111 22222 33333 44444 55555 66666 77777 88888

OMP-386

データ再配置：Sequential Reordering
CM-RCM(2), 4-threads
スレッド上のデータ連続性：キャッシュ有効利用，プリフェッチが効きやすくなる
45 10 39 5 35 2 33 1
17 46 11 40 6 36 3 34
53 18 47 12 41 7 37 4
24 54 19 48 13 42 8 38
59 25 55 20 49 14 43 9
29 60 26 56 21 50 15 44
63 30 61 27 57 22 51 16
32 64 31 62 28 58 23 52

29 18 15 5 11 2 9 1
33 30 19 16 6 12 3 1045 34 31 20 25 7 13 4
40 46 35 32 21 26 8 1459 49 47 36 41 22 27 17
53 60 50 48 37 42 23 28
63 54 61 51 57 38 43 24
56 64 55 62 52 58 39 44

CM-RCM(2)Sequential Reordering, 4-threads

OMP-387

データ再配置：Sequential Reordering
CM-RCM(2), 4-threads
第1色■#0 thread, ■#1，■#2，■#3
45 10 39 5 35 2 33 1
17 46 11 40 6 36 3 34
53 18 47 12 41 7 37 4
24 54 19 48 13 42 8 38
59 25 55 20 49 14 43 9
29 60 26 56 21 50 15 44
63 30 61 27 57 22 51 16
32 64 31 62 28 58 23 52

29 18 15 5 11 2 9 1
33 30 19 16 6 12 3 1045 34 31 20 25 7 13 4
40 46 35 32 21 26 8 1459 49 47 36 41 22 27 17
53 60 50 48 37 42 23 28
63 54 61 51 57 38 43 24
56 64 55 62 52 58 39 44

CM-RCM(2)Sequential Reordering, 4-threads

OMP-388

データ再配置：Sequential Reordering
CM-RCM(2), 4-threads
第2色■#0 thread, ■#1，■#2，■#3
45 10 39 5 35 2 33 1
17 46 11 40 6 36 3 34
53 18 47 12 41 7 37 4
24 54 19 48 13 42 8 38
59 25 55 20 49 14 43 9
29 60 26 56 21 50 15 44
63 30 61 27 57 22 51 16
32 64 31 62 28 58 23 52

29 18 15 5 11 2 9 1
33 30 19 16 6 12 3 1045 34 31 20 25 7 13 4
40 46 35 32 21 26 8 1459 49 47 36 41 22 27 17
53 60 50 48 37 42 23 28
63 54 61 51 57 38 43 24
56 64 55 62 52 58 39 44

CM-RCM(2)Sequential Reordering, 4-threads

OMP-389

データ再配置: Sequential Reordering
Coalesced
GPUはこちらがお勧め

Initial Vector

Coloring(5 colors)+Orderingcolor=1

color=1

color=2

color=2

color=3

color=4

color=3

color=4

color=5

color=512345678 12345678 12345678 12345678 12345678

Sequential

Initial Vector
Coloring(5 colors)+Ordering
各スレッド上で不連続なメモリ
アクセス（色の順に番号付け）

color=1

color=1

color=2

color=2

color=3

color=4

color=5

color=3

color=4

color=512345678 12345678 12345678 12345678 12345678

11111 22222 33333 44444 55555 66666 77777 88888
スレッド内で連続に番号付け

OMP-390

プログラムのありか• 所在– <$L3>/reorder0• コンパイル，実行方法– 本体• cd <$L3>/reorder0• make• <$L3>/reorder0/L3-rsol0（実行形式）– コントロールデータ• <$L3>/run/INPUT.DAT– 実行用シェル• <$L3>/run/gor.shOMP-391

制御データ（INPUT.DAT）
100 100 100
1.00e-00 1.00e-00 1.00e-00
1.0e-08
16
100
0
0

NX/NY/NZ
DX/DY/DZEPSICC
PEsmpTOTNCOLORtot
NFLAG
METHOD

変数名

型

内

容

PEsmpTOT

整数

データ分割数NCOLORtot

整数

Ordering手法と色数≧2：MC法（multicolor），色数=0：CM法（Cuthill-Mckee）=-1：RCM法（Reverse Cuthill-Mckee）≦-2：CM-RCM法

NFLAG

整数

0：First-Touch無し，1：あり
今回は無関係METHOD

整数

行列ベクトル積のループ構造（0：従来通り，1：前進後退代入と同じ）

OMP-392

再配置：Sequential Reordering
allocate (SMPindex(0:PEsmpTOT*NCOLORtot))
SMPindex= 0
do ic= 1, NCOLORtot
nn1= COLORindex(ic) - COLORindex(ic-1)
num= nn1 / PEsmpTOT
nr = nn1 - PEsmpTOT*num
do ip= 1, PEsmpTOT
if (ip.le.nr) then
SMPindex((ic-1)*PEsmpTOT+ip)= num + 1
else
SMPindex((ic-1)*PEsmpTOT+ip)= num
endif
enddo
enddo

SMPindex
ic=1

ic=2

ic=3

ic=4

ic=5

12345678 12345678 12345678 12345678 12345678

SMPindex_new
12345 12345 12345 12345 12345 12345 12345 12345

allocate (SMPindex_new(0:PEsmpTOT*NCOLORtot))
SMPindex_new(0)= 0
do ic= 1, NCOLORtot
do ip= 1, PEsmpTOT
j1= (ic-1)*PEsmpTOT + ip
j0= j1 - 1
SMPindex_new((ip-1)*NCOLORtot+ic)= SMPindex(j1)SMPindex(j1)= SMPindex(j0) + SMPindex(j1)
enddo
enddo
do ip= 1, PEsmpTOT
do ic= 1, NCOLORtot
j1= (ip-1)*NCOLORtot + ic
j0= j1 - 1
SMPindex_new(j1)= SMPindex_new(j0) + SMPindex_new(j1)
enddo
enddo

OMP-393

行列ベクトル積の計算法!$omp parallel do private(ip,i,VAL,k)do ip= 1, PEsmpTOT
do i = SMPindexG(ip-1)+1, SMPindexG(ip)
VAL= D(i)*W(i,P)
do k= indexL(i-1)+1, indexL(i)VAL= VAL + AL(k)*W(itemL(k),P)
enddo
do k= indexU(i-1)+1, indexU(i)VAL= VAL + AU(k)*W(itemU(k),P)
enddo
W(i,Q)= VAL
enddo
enddo
!$omp end parallel do

METHOD=0!$omp parallel do private(ip,i,VAL,k)do ip= 1, PEsmpTOT
do i= SMPindex((ip-1)*NCOLORtot)+1, SMPindex(ip*NCOLORtot)
VAL= D(i)*W(i,P)
do k= indexL(i-1)+1, indexL(i)VAL= VAL + AL(k)*W(itemL(k),P)
enddo
do k= indexU(i-1)+1, indexU(i)VAL= VAL + AU(k)*W(itemU(k),P)
enddo
W(i,Q)= VAL
enddo
enddo
METHOD=1!$omp end parallel doOMP-3

前進代入の計算法：色ループは外!$omp parallel private(ip,ip1,i,WVAL,k)do ic= 1, NCOLORtot!$omp dodo ip= 1, PEsmpTOT
ip1= (ic-1)*PEsmpTOT + ip
do i= SMPindex(ip1-1)+1, SMPindex(ip1)
WVAL= W(i,Z)
do k= indexL(i-1)+1, indexL(i)WVAL= WVAL - AL(k) * W(itemL(k),Z)
enddo
W(i,Z)= WVAL * W(i,DD)
enddo
enddo
enddo
!$omp end parallel!$omp parallel private(ip,ip1,i,WVAL,k)do ic= 1, NCOLORtot!$omp dodo ip= 1, PEsmpTOT
ip1= (ip-1)*NCOLORtot + ic
do i= SMPindex(ip1-1)+1, SMPindex(ip1)
WVAL= W(i,Z)
do k= indexL(i-1)+1, indexL(i)WVAL= WVAL - AL(k) * W(itemL(k),Z)
enddo
W(i,Z)= WVAL * W(i,DD)
enddo
enddo
enddo
!$omp end parallel

OriginalNew

94

OMP-395

行列格納方法ELL（Ellpack-Itpack）: ループ長固定，プリフェッチ効きやすい
1
1

4

0
1

3
2
1
3
0

0
5
3
7
0

0
0
0
4
0

0
0
0

0
5

1 3

1 3 0

1 2 5

1 2 5

4 1 3

4 1 3

3 7 4

3 7 4

1 5

1 5 0

(a) CRS(b) ELLOMP-3

96

実施ケース
並列化向
け色付け
手法
src0

Case-1

reorder0 Case-2

CM-RCM

Case-3

ELL

番号付けFirst Touch
Data
Placement

Coalesced（図 4（a））

無し

Sequential（図 4（b））

有り

color=1

ELL

Initial Vector

Initial Vector

color=1

CRS

Sequential

Coalesced
Coloring(5 colors)+Ordering

係数行列格納形式

color=2

color=2

color=3

color=4

color=5

color=3

color=4

color=512345678 12345678 12345678 12345678 12345678

Coloring(5 colors)+Orderingcolor=1

color=2

color=3

color=4

color=5

各スレッド上でcolor=1
color=2
color=3
color=4
color=5
不連続なメモリ
アクセス（色の 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8
順に番号付け）

11111 22222 33333 44444 55555 66666 77777 88888
スレッド内で連続に番号付け

97

T2K/Tokyo

Cray XE6 (Hopper)

MemoryMemoryMemory

L2
L1

L3
L2 L2
L1 L1

L2
L1

L2
L1

L3
L2 L2
L1 L1

L2
L1

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

L1
L2

L1 L1
L2 L2
L3

L1
L2

L1
L2

L1 L1
L2 L2
L3

L1
L2

MemoryMemory

L2
L1

L2
L1

L3
L2 L2
L1 L1

L2
L1

L2
L1

L2
L1

L2
L1

L3
L2 L2
L1 L1

L2
L1

L2
L1

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

L1
L2

L1
L2

L1 L1
L2 L2
L3

L1
L2

L1
L2

L1
L2

L1
L2

L1 L1
L2 L2
L3

L1
L2

L1
L2

MemoryMemoryMemoryFujitsu FX10 (Oakleaf-FX)
Memory

L1

L1

L1

L1

L1

L1

L1

L2
L1 L1

L1

L1

L1

L1

L1

L1

L1

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

OMP-398

計算結果：色数～反復回数（CM-RCM）
450

Iterations

400

350

300

250
1

10

COLOR#

100

1000

OMP-399

計算結果：色数～計算時間：FX10
10.00• CASE-1(src0)⇒
CASE-2(reorder0)

6.00

sec.

– 色数が大きくなると多少効果あり– 色数が多くなるほど，色当たり，スレッド当たりの計算量は減少8.00

4.00

Case-1Case-22.00

Case-30.00• CASE-2は色が変わってもスレッド上のデータが連続– First Touchの効果はナシ• ELLの効果は大

1

10

100

COLOR#

Case-1: src0
Case-2: reorder0
Case-3: reorder0 + ELL

1000

OMP-3

100

Fujitsu FX10: CASE-1, CM-RCM(2)L1デマンドミス：25.6%, Mem. Throughput:41.8GB/sec.
前処理部分（前進後退代入）
整数ロードメモリアクセス待ち
浮動小数点ロードキャッシュアクセス待ち
命令フェッチ待ち
1命令コミット[秒]
4.0E+00

浮動小数点ロードメモリアクセス待ち
整数演算待ちバリア同期待ち整数レジスタ書き込み制約

ストア待ち
浮動小数点演算待ち
uOPコミット
2/3命令コミット

整数ロードキャッシュアクセス待ち
分岐命令待ちその他の待ち4命令コミット

3.5E+00

3.0E+00

2.5E+00

2.0E+00

1.5E+00

1.0E+00

5.0E-01

0.0E+00
Thread 0

Thread 1

Thread 2

Thread 3

Thread 4

Thread 5

Thread 6

Thread 7

Thread 8

Thread 9

Thread 10

Thread 11

Thread 12

src0: CRS, Coalesced

Thread 13

Thread 14

Thread 15

OMP-3

101

Fujitsu FX10: CASE-2, CM-RCM(2)25.6%, 41.8GB/sec.
整数ロードメモリアクセス待ち
浮動小数点ロードキャッシュアクセス待ち
命令フェッチ待ち
1命令コミット[秒]
4.0E+00

浮動小数点ロードメモリアクセス待ち
整数演算待ちバリア同期待ち整数レジスタ書き込み制約

ストア待ち
浮動小数点演算待ち
uOPコミット
2/3命令コミット

整数ロードキャッシュアクセス待ち
分岐命令待ちその他の待ち4命令コミット

3.5E+00

3.0E+00

2.5E+00

2.0E+00

1.5E+00

1.0E+00

5.0E-01

0.0E+00
Thread 0

Thread 1

Thread 2

Thread 3

Thread 4

Thread 5

Thread 6

Thread 7

Thread 8

Thread 9

Thread 10

Thread 11

Thread 12

Thread 13

reorder0: CRS, Sequential

Thread 14

Thread 15

OMP-3

102

Fujitsu FX10: CASE-1, CM-RCM(382)
37.7%, 28.7GB/sec.
整数ロードメモリアクセス待ち
浮動小数点ロードキャッシュアクセス待ち
命令フェッチ待ち
1命令コミット[秒]
4.0E+00

浮動小数点ロードメモリアクセス待ち
整数演算待ちバリア同期待ち整数レジスタ書き込み制約

ストア待ち
浮動小数点演算待ち
uOPコミット
2/3命令コミット

整数ロードキャッシュアクセス待ち
分岐命令待ちその他の待ち4命令コミット

3.5E+00

3.0E+00

2.5E+00

2.0E+00

1.5E+00

1.0E+00

5.0E-01

0.0E+00
Thread 0

Thread 1

Thread 2

Thread 3

Thread 4

Thread 5

Thread 6

Thread 7

Thread 8

Thread 9

Thread 10

Thread 11

Thread 12

src0: CRS, Coalesced

Thread 13

Thread 14

Thread 15

OMP-3

103

Fujitsu FX10: CASE-2, CM-RCM(382)
29.3%, 32.6GB/sec.
整数ロードメモリアクセス待ち
浮動小数点ロードキャッシュアクセス待ち
命令フェッチ待ち
1命令コミット[秒]
4.0E+00

浮動小数点ロードメモリアクセス待ち
整数演算待ちバリア同期待ち整数レジスタ書き込み制約

ストア待ち
浮動小数点演算待ち
uOPコミット
2/3命令コミット

整数ロードキャッシュアクセス待ち
分岐命令待ちその他の待ち4命令コミット

3.5E+00

3.0E+00

2.5E+00

2.0E+00

1.5E+00

1.0E+00

5.0E-01

0.0E+00
Thread 0

Thread 1

Thread 2

Thread 3

Thread 4

Thread 5

Thread 6

Thread 7

Thread 8

Thread 9

Thread 10

Thread 11

Thread 12

Thread 13

reorder0: CRS, Sequential

Thread 14

Thread 15

OMP-3104

計算結果のまとめ：Fujitsu FX10
詳細プロファイラ
上段: L1デマンドミス率
下段：メモリースループット
src0
CASE-1CRS+
Coalesced

reorder0
CASE-2
CRS+
Sequential

CASE-3
ELL+
Sequential

25.5 %

25.6 %

5.42 %

CM-RCM(2)41.8 GB/sec. 41.8 GB/sec. 64.0 GB/sec.
37.7 %

29.3 %

16.5 %

CM-RCM(382)28.7 GB/sec. 32.6 GB/sec. 52.2 GB/sec.

OMP-3105

計算結果のまとめ：Fujitsu FX10
詳細プロファイラ
上段：CM-RCM(20),下段：CMRCM(382)
Instructions

Memory Access
SIMD（%）
Throughput（%）

Case-2
CRS

1.83×1011
1.83×1011

7.17
6.90

50.2
44.5

Case-3ELL

6.71×1010
5.96×1010

16.3
16.2

69.8
67.0

Case-1: src0
Case-2: reorder0
Case-3: reorder0 + ELL

OMP-3106

計算結果：色数～計算時間：Cray XE6
30.00• CASE-1(src0)⇒
CASE-2(reorder0)の変化大• CRS⇒ELLの効果はそれほど顕著では無い

Case-2

20.00

sec.

– NUMA向け最適化の効果– First Touchとの合わせ技
ではあるがCase-1Case-3

10.00

0.00
1

10

100

COLOR#

Case-1: src0
Case-2: reorder0
Case-3: reorder0 + ELL

1000

OMP-3107

計算結果のまとめ
CM-RCM(20)
一反復当
計算時
たり計算
間（秒）
時間（秒）
Fujitsu
FX10

Cray
XE6

CM-RCM(382) = RCM
計算時
間（秒）

一反復当
たり計算
時間（秒）

Case-15.44

1.71×10-2

6.37

2.22×10-2

Case-25.76

1.81×10-2

5.78

2.02×10-2

Case-33.90

1.23×10-2

3.61

1.26×10-2

Case-1

19.7

6.26×10-2

18.7

6.52×10-2

Case-27.45

2.34×10-2

7.40

2.58×10-2

Case-37.14

2.25×10-2

6.77

2.36×10-2

Case-1: src0
Case-2: reorder0
Case-3: reorder0 + ELL

OMP-3

108

Fujitsu FX10: CASE-3, CM-RCM(2)5.4%, 64.0GB/sec.
整数ロードメモリアクセス待ち
浮動小数点ロードキャッシュアクセス待ち
命令フェッチ待ち
1命令コミット[秒]
4.0E+00

浮動小数点ロードメモリアクセス待ち
整数演算待ちバリア同期待ち整数レジスタ書き込み制約

ストア待ち
浮動小数点演算待ち
uOPコミット
2/3命令コミット

整数ロードキャッシュアクセス待ち
分岐命令待ちその他の待ち4命令コミット

3.5E+00

3.0E+00

2.5E+00

2.0E+00

1.5E+00

1.0E+00

5.0E-01

0.0E+00
Thread 0

Thread 1

Thread 2

Thread 3

Thread 4

Thread 5

Thread 6

Thread 7

Thread 8

Thread 9

Thread 10

Thread 11

ELL, Sequential

Thread 12

Thread 13

Thread 14

Thread 15

OMP-3

109

Fujitsu FX10: CASE-3, CM-RCM(382)
16.5%, 52.2GB/sec.
整数ロードメモリアクセス待ち
浮動小数点ロードキャッシュアクセス待ち
命令フェッチ待ち
1命令コミット[秒]
4.0E+00

浮動小数点ロードメモリアクセス待ち
整数演算待ちバリア同期待ち整数レジスタ書き込み制約

ストア待ち
浮動小数点演算待ち
uOPコミット
2/3命令コミット

整数ロードキャッシュアクセス待ち
分岐命令待ちその他の待ち4命令コミット

3.5E+00

3.0E+00

2.5E+00

2.0E+00

1.5E+00

1.0E+00

5.0E-01

0.0E+00
Thread 0

Thread 1

Thread 2

Thread 3

Thread 4

Thread 5

Thread 6

Thread 7

Thread 8

Thread 9

Thread 10

Thread 11

ELL, Sequential

Thread 12

Thread 13

Thread 14

Thread 15

OMP-3110• マルチコア版コードの実行• 更なる最適化• STREAM• プロファイラ，コンパイルリスト分析等

OMP-3111

計算結果（FX10@東大）：106要素
反復回数：MC（2色）：333回，RCM（298レベル）：224回
CM-RCM（Nc=20）：249回, L3-sol
30.0

20.0
MC=2

MC=2

RCM(298)

16.0

CM-RCM(20)RCM(298)
CM-RCM(20)sec.

Speed-Up

20.0

10.0

12.0

8.0

4.0
0.0

0.0
1

2

4

8

16

1

2

4

thread#16 threads
MC(2):
2.42 sec.
CM-RCM(20): 2.01 sec.

8

16

thread#Memory

L1

L1

L1

L1

L1

L1

L1

L2
L1 L1

L1

L1

L1

L1

L1

L1

L1

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

OMP-3112

何故16倍にならないか?
• 16スレッドがメモリにアクセスすると，1スレッドの場合と比較して，スレッド当り（コア当り）メモリ性能は低下• 疎行列はmemory-boundなためその傾向がより顕著– 疎行列計算の高速化：研究途上の課題• 問題規模が比較的小さい

113

疎行列・密行列
do i= 1, N
Y(i)= D(i)*X(i)do k= index(i-1)+1, index(i)kk= item(k)Y(i)= Y(i) + AMAT(k)*X(kk)
enddo
enddo

do j= 1, N
Y(j)= 0.d0
do i= 1, N
Y(j)= Y(j) + A(i,j)*X(i)
enddo
enddo

• “X” in RHS– 密行列：連続アクセス，キャッシュ有効利用– 疎行列：連続性は保証されず，キャッシュを有効に活用できず• より「memory-bound」

OMP-3114

GeoFEM Benchmark
ICCG法の性能（固体力学向け）
SR11K/J2

SR16K/M1

T2K

FX10

京

Core #/Node

16

32

16

16

8

Peak Performance(GFLOPS)

147.2

980.5

147.2

236.5

128.0

STREAM Triad (GB/s)

101.0

264.2

20.0

64.7

43.3

B/F

0.686

0.269

0.136

0.274

0.338

GeoFEM (GFLOPS)

19.0

72.7

4.69

16.0

11.0

% to Peak

12.9

7.41

3.18

6.77

8.59

LLC/core (MB)

18.0

4.00

2.00

0.75

0.75

疎行列ソルバー：Memory-BoundOMP-ESTREAM ベンチマーク

http://www.cs.virginia.edu/stream/• メモリバンド幅を測定するベンチマーク– Copy:– Scale:– Add:– Triad:c(i)= a(i)c(i)= s*b(i)c(i)= a(i) + b(i)c(i)= a(i) + s*b(i)---------------------------------------------Double precision appears to have 16 digits of accuracy
Assuming 8 bytes per DOUBLE PRECISION word---------------------------------------------Number of processors =16
Array size =2000000
Offset=
0
The total memory requirement is732.4 MB(
45.8MB/task)
You are running each test 10 times-The *best* time for each test is used
*EXCLUDING* the first and last iterations------------------------------------------------------------------------------------------------------Function
Rate (MB/s) Avg time
Min time Max timeCopy:
18334.1898
0.02800.02790.0280Scale:
18035.1690
0.02840.02840.0285
Add:
18649.4455
0.04120.04120.0413
Triad:
19603.8455
0.03940.03920.0398115

OMP-E116

マイクロプロセッサの動向
CPU性能，メモリバンド幅のギャップ

http://www.cs.virginia.edu/stream/OMP-1

117

実行：OpenMPバージョン>$ cd <$O-stream>>$ pjsub run.shOMP-1118

run.sh#!/bin/sh#PJM -L "rscgrp=school"#PJM -L "node=1"#PJM -L "elapse=10:00"#PJM -jexport PATH=…
export LD_LIBRARY_PATH=…
export PARALLEL=16
export OMP_NUM_THREADS=16

スレッド数（1-16）./stream.out > 16-01.lst 2>&1

出力ファイル名OMP-1

119

Results of Triad

Memory<$O-stream>/stream/*.lst
Peak is 85.3 GB/sec., 75%
Thread #

L1

L1

L1

L1

L1

L1

L1

L2
L1 L1

L1

L1

L1

L1

L1

L1

L1

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

MB/sec.

Speed-up

1

8606.141.00

216918.81

1.97

434170.72

3.97

859505.92

6.91

1664714.32

7.52

C

120

実習（2）
• 実際にやってみよ• スレッド数を変える• 1CPU版，MPI版もある– FORTRAN，C– STREAMのサイト

OMP-3121• マルチコア版コードの実行• 更なる最適化• STREAM• プロファイラ，コンパイルリスト分析等– 利用支援ポータル⇒ドキュメント閲覧⇒プログラ
ム開発支援ツール⇒プロファイラ使用手引書⇒「3章：詳細プロファイラ」

122

デフォルト
>$ cd <$O-L3>/src>$ make>$ ls ../run/L3-sol
L3-sol>$ cd ../run>$ pjsub go1.sh
F90= frtpx
F90OPTFLAGS= -Kfast,openmp -Qt
F90FLAGS =$(F90OPTFLAGS)

Kshortloop>$ cd <$O-L3>/short>$ make>$ ls ../run/L3-sols
L3-sols>$ cd ../run>$ pjsub gos.shコンパイル・実行• -Qt– コンパイルリスト出力– ＊.lst• Cでは「-Qt」は使えません，「-Nsrc」を使ってください。
– 画面に出てしまいますがF90= frtpx
F90OPTFLAGS= -Kfast,openmp -Kshortloop=3-Qt
F90FLAGS =$(F90OPTFLAGS)

123＊.lstの見方

124

SIMD情報125

自動並列化情報126

solver_ICCG_mc.lst (src)
101
102
103
104
105
106
107
108

1
1
1
1
1
1
1
2

!C
!C +----------------+!C | {z}= [Minv]{r} |!C +----------------+
!C===

109
110
111
112
113
114
115
116
117
118
119
120
121
122

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

p
p
p
p

123
124
125
126
127
128
129
130

5
5
5
4
4
3
2
2

p
p
p
p
p
p

p

8v
8v
8v

p
p
p
p

4v
4v
4v

!$omp parallel do private(ip,i)do ip= 1, PEsmpTOT<<< Loop-information Start >>><<< [OPTIMIZATION]
<<<
SIMD
<<<
SOFTWARE PIPELINING<<< Loop-information End >>>do i = SMPindexG(ip-1)+1, SMPindexG(ip)
W(i,Z)= W(i,R)
enddo
enddo
!$omp end parallel do
Stime= omp_get_wtime()call fapp_start ("precond", 1, 1)
do ic= 1, NCOLORtot!$omp parallel do private(ip,ip1,i,WVAL,k)do ip= 1, PEsmpTOT
ip1= (ic-1)*PEsmpTOT + ip
do i= SMPindex(ip1-1)+1, SMPindex(ip1)
WVAL= W(i,Z)<<< Loop-information Start >>><<< [OPTIMIZATION]
<<<
SIMD
<<<
SOFTWARE PIPELINING<<< Loop-information End >>>do k= indexL(i-1)+1, indexL(i)WVAL= WVAL - AL(k) * W(itemL(k),Z)
enddo
W(i,Z)= WVAL * W(i,DD)
enddo
enddo
!$omp end parallel do
enddoOMP-3127

3.5 精密PA可視化機能（Excel形式）（1/3） 測定範囲指定，コンパイル・リンク
call start_collection ("SpMV")!$omp parallel do private(ip,i,VAL,k)do ip= 1, PEsmpTOT
do i= SMPindex((ip-1)*NCOLORtot)+1, SMPindex(ip*NCOLORtot)
VAL= D(i)*W(i,P)
do k= 1, 3
VAL= VAL + AL(k,i)*W(itemL(k,i),P)
enddo
do k= 1, 3
VAL= VAL + AU(k,i)*W(itemU(k,i),P)
enddo
W(i,Q)= VAL
enddo
enddo
!$omp end parallel docall stop_collection ("SpMV")OMP-3128

3.5 精密PA可視化機能（Excel形式）（2/3） データ収集：7回実行
ディレクトリ名：pa1~pa7，-Hpa=1~7#!/bin/sh#PJM -L "node=1"#PJM -L "elapse=00:05:00"#PJM -L "rscgrp=debug"#PJM -g "pz0088"#PJM -j#PJM -o "3.lst"#PJM --mpi "proc=1"
export OMP_NUM_THREADS=16
fapp -C -d pa1 -Ihwm -Hpa=1 ./sol-r3k

OMP-3129

3.5 精密PA可視化機能（Excel形式）（3/3） データ解析：変換＋Excelシート（p.80）
fapppx -A -d pa1 -o output_prof_1.csv -tcsv -Hpa
fapppx -A -d pa2 -o output_prof_2.csv -tcsv -Hpa
fapppx -A -d pa3 -o output_prof_3.csv -tcsv -Hpa
fapppx -A -d pa4 -o output_prof_4.csv -tcsv -Hpa
fapppx -A -d pa5 -o output_prof_5.csv -tcsv -Hpa
fapppx -A -d pa6 -o output_prof_6.csv -tcsv -Hpa
fapppx -A -d pa7 -o output_prof_7.csv -tcsv -Hpa

OMP-3130

まとめ• 「有限体積法から導かれる疎行列を対象としたICCG法」を題材とした，データ配置，reorderingなど，科学技術計算のためのマルチコアプログラミングにおいて重要なアルゴリズムについての講習• 更に理解を深めるための，FX10を利用した実習• オーダリングの効果

OMP-3131

今後の動向• メモリバンド幅と性能のギャップ– BYTE/FLOP，中々縮まらない• マルチコア化，メニーコア化– Intel Xeon/Phi• >105コアのシステム– Exascale：＞108• オーダリング– グラフ情報だけでなく，行列成分の大きさの考慮も必要か?
– 最適な色数の選択：研究課題（特に悪条件問題）• OpenMP+MPIのハイブリッド⇒一つの有力な選択– プロセス内（OpenMP）の最適化が最もcritical• 本講習会の内容が少しでも役に立てば幸いである

132

略 称

FX10

名 称

Fujitsu SPARC64 IX
fx

動作周波数（GHz）

1.848

1.053

2.80

16 (16)60 (240)

10 (20)

16

240

10

DDR3

GDDR5

DDR3

236.5

1,010.9

224.0

32

8

64

85.1

320

59.7

コア数（有効スレッド数）
使用スレッド数
メモリ種別
理論演算性能（GFLOPS）
主記憶容量（GB）
理論メモリ性能（GB/sec.）

MIC

Intel Xeon Phi
Intel Xeon E5-2680
5110P
v2 (Ivy-Bridge-EP)(Knights Corner)

L1:32KB/core
L2:12MB/socket

L1:32KB/core
L2:512KB/core

コンパイルオプション -Kfast, openmp-O3 -openmp mmic -alignarray64byte

キャッシュ構成

IvyB

L1:32KB/core
L2:256KB/core
L3:25MB/socket-O3 -openmp -ipo
-xAVX
-align array32byte

133

Coalesced, Sequential• NUMA向け最適化– First Touch
Data Placement• 計算と同じ順番
で初期化– Sequential
Reordering（再オーダリング）• 色の順番にリオ－ダリング（Coalesced）• 更にスレッド上で番号が連続となるように再番号付け

CoalescedColoring(5 colors)+Orderingcolor=1

color=1

Initial Vector

color=2

color=2

color=3

color=4

color=5

color=3

color=4

color=512345678 12345678 12345678 12345678 12345678

SequentialColoring(5 colors)+Orderingcolor=1

Initial Vector

color=2

color=3

color=4

color=5

各スレッド上でcolor=1
color=2
color=3
color=4
color=5
不連続なメモリ
アクセス（色の 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8
順に番号付け）

11111 22222 33333 44444 55555 66666 77777 88888
スレッド内で連続に番号付け

134

Numberin
g
AR-0

行列格納
形式

その他

CRS
行方向
（Rowwise）

AR-1
AR-2

外側ループ

Coalesced
Sliced ELL

間接メモリアクセス

列方向（Columnwise）
ブロック化あり

AC-1
AC-2
BR-0

CRS

BR-1

行方向
（Rowwise）

Sequential
BC-1
BC-2

Sliced ELL

列方向（Columnwise）
ブロック化

135

ELL：外側ループ：行方向：Row-wise
従来，中島がやってきたやり方：CRS with fixed length
前進代入のループ!$omp parallel
do icol= 1, NCOLORtot!$omp dodo ip = 1, PEsmpTOT
do i= Index(ip-1,icol)+1, Index(ip,icol)
do k= 1, 6
Z(i)= Z(i) - AML(k,i)*Z(IAML(k,i))
enddo
Z(i)= Z(i) / DD(i)
enddo
enddo
enddo
!omp end parallel

k

i

136

ELL：外側ループ：列方向：Column-wise
こちらが本来のELL: Jagged Diagonal
係数行列のアクセス不連続!$omp parallel
do icol= 1, NCOLORtot!$omp dodo ip = 1, PEsmpTOT
do k= 1, 6
do i= Index(ip-1,icol)+1, Index(ip,icol)
Z(i)= Z(i) + AML (i,k)*Z(IAML(i,k))
enddo
enddo
do i= Index(ip-1,icol)+1, Index(ip,icol)
Z(i)= Z(i) / DD(i)
enddo
enddo
enddo
!omp end parallelkiEach Color/Thread

137

ELL：外側ループ：列方向：Column-wise
ブロック化版，各色・スレッドごとに別々のブロックで係数行列を記憶!$omp parallel
do icol= 1, NCOLORtot!$omp dodo ip = 1, PEsmpTOT
blkID= (ip-1)*NCOLORtot + ip
do k= 1, 6
do i= IndexB(ip-1,blkID,icol)+1, &iIndexB(ip ,blkID,icol)
locID= i – IndexB(ip-1,blkID,icol)
Z(i)= Z(i) +&AMLb(locID,k,blkID)* X(IAMLb(locID,k,blkID))
enddo
enddo
do i= IndexB(ip-1,blkID,icol)+1, IndexB(ip ,blkID,icol)
Z(i)= Z(i) / DD(i)
enddo
enddo
enddo
!omp end parallel

k

400

• ELLの場合

380
360

Iterations

色数と計算時間（Solver）

138– AR-1: Coalesced + Row-wise– BR-1: Sequential + Row-wise

340
320
300
280
260
1

10

100

sec.

Color#• 色数が減ると反復回数が減る 12.00
FX10: AR-1 FX10: BR-1
が，同期オーバーヘッドは増MIC: AR-1
MIC: BR-1
10.00
加IvyB: AR-1
IvyB: BR-1
8.00
• MICはこれが顕著
6.00• 計算時間を考慮した最適色
4.00
数はアーキテクチャによって
異なる
2.00
1
10
100
• FX10，IvyBはRCMが良い
Color#

1000

1000

10.00色数と計算時間
100色まで

139

AR-1
8.00

AR-2

sec.

AC-1
AC-2

6.00

• FX10，MICはELLの効果大（CRS遅い：AR-0，BR-0） 4.00
• IvyBは効果少，色数の影響
2.00もほとんど無い
0
10.00

AR-0

IvyB

AR-0

FX10

BR-0
BR-1
BC-1
BC-2

20

40

60

80

100

Color#

10.00

AR-0

MIC

AR-1
8.00

AR-1

AR-2

8.00

AR-2

AC-2

6.00

BR-0

AC-1
sec.

sec.

AC-1

AC-2

6.00

BR-0

BR-1
BC-1

4.00

BR-1
BC-1

4.00

BC-2
2.00

BC-2
2.00

0

20

40

60

Color#

80

100

0

20

40

60

Color#

80

100

140

各手法の結果（Solver計算時間）
7.00

FX10: 30 colors
MIC: 10 colors

6.00

IvyB: 40 colors

sec.

5.00
4.00
3.00
2.00
1.00
0.00
AR-0 AR-1 AR-2 AC-1 AC-2 BR-0 BR-1 BC-1 BC-2
CASE ID

141

各手法の結果（Solver計算時間）:FX10
7.00

FX10: 30 colors
MIC: 10 colors

6.00

IvyB: 40 colors

sec.

5.00
4.00
3.00
2.00
1.00
0.00
AR-0 AR-1 AR-2 AC-1 AC-2 BR-0 BR-1 BC-1 BC-2
CASE ID• ELLの効果大（AR-1, AR-2,
BR-1）• Column-wiseの性能低い（AC-X，
BC-X）：命令数，
L1ミス• ブロック化の効果はある（XC-1⇒XC-2）• CoalescedとSequentialは同じに見える• AR-2悪くなさそう

142

各手法の結果（Solver計算時間）:MIC
7.00

FX10: 30 colors
MIC: 10 colors

6.00

IvyB: 40 colors

sec.

5.00

• ELLの効果大• Column-wiseの性能高い（AC-X，
BC-X）：ベクトル化，SIMD• ブロック化の効
果多少ある（XC-1⇒XC-2）

4.00
3.00
2.00
1.00
0.00
AR-0 AR-1 AR-2 AC-1 AC-2 BR-0 BR-1 BC-1 BC-2
CASE ID

143

各手法の結果（Solver計算時間）:IvyB
7.00

FX10: 30 colors
MIC: 10 colors

6.00

IvyB: 40 colors

sec.

5.00

• ELLの効果小さい• Column-wiseは遅い• Sequential（B）良い• 間接参照効果小（AR-2），AR-1とAR-2は同じ

4.00
3.00
2.00
1.00
0.00
AR-0 AR-1 AR-2 AC-1 AC-2 BR-0 BR-1 BC-1 BC-2
CASE ID