OpenMPによるマルチコア・
メニィコア並列プログラミング入門
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-36

実行例：OBCX>$ cd /work/gt00/t00XYZ>$ cp /work/gt00/z30088/omp/multicore-f.tar .
>$ tar xvf multicore-f.tar>$ cd multicore/L3
>$ ls
run
src src0>$ cd src>$ make>$ cd ../run>$ ls L3-sol
L3-sol<modify “INPUT.DAT”><modify “go1.sh”>>$ pjsub go1.sh

srcx

reorder0

OMP-3

7

プログラムのありか on OBCX
• 所在
– <$O-L3>: /work/gt00/t00XYZ/multicore/L3– <$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

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

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

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

OMP-39制御データ（INPUT.DAT）
128 128 128
1.00e-00 1.00e-00 1.00e-00
1.0e-08
24
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法

各方向の要素数ジョブスクリプト：go1.sh• /work/gt00/t00XXX/multicore/L3/run/go1.sh• スケジューラへの指令 ＋ シェルスクリプト#!/bin/sh#PJM -N “test1”#PJM -L rscgrp=tutorial#PJM -L node=1#PJM --omp thread=24#PJM -L elapse=00:15:00#PJM -g gt00
#PJM -j#PJM -e err#PJM -o test1.lst

ジョブ名称（省略可）
実行キュー名
ノード数（原則＝1）
スレッド数（1-56，原則1-28）
実行時間
グループ名（財布）
エラー出力ファイル
標準出力ファイル

export KMP_AFFINITY=granularity=fine,compact./L3-sol
プログラム実行
export KMP_AFFINITY=granularity=fine,compact
各スレッドがSocket#0の0番から始まる各コアに順番に割り当てられる

10

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

11

OMP-312

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

OMP-3

13

実はこのようにして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-3

14

実はこのようにして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-315

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

OMP-3

16

プログラムの構成
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-317

module STRUCTICELTOT
要素数（NX×NY×NZ）

module STRUCT
use omp_lib
include 'precision.inc'

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

&

&
&

&

!C
!C-- BOUNDARYs
integer (kind=kint) :: ZmaxCELtot
integer (kind=kint), dimension(:), allocatable :: BC_INDEX, BC_NOD
integer (kind=kint), dimension(:), allocatable :: ZmaxCEL!C!C-- WORK
integer (kind=kint), dimension(:,:), allocatable :: IWKX
real(kind=kreal),
dimension(:,:), allocatable :: FCV

NX,NY,NZ
x,y,z方向要素数
NXP1,NYP1,NZP1
x,y,z方向節点数
IBNODTOT
NXP1×NYP1XYZ(ICELTOT,3)
要素座標
NEIBcell(ICELTOT,6)
隣接要素

integer (kind=kint) :: PEsmpTOT
end module STRUCT

PEsmpTOT
スレッド数

OMP-318

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

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

OMP-320

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

内 容

NEWtoOLD(N)I 入力時にはOrdering手法（≧2：MC，=0：
CM，=-1：RCM,-2≧：CMRCM），
最終的には色数，レベル数が入る
I 各色，レベルに含まれる要素数の一次元圧縮配列，
COLORindex(icol-1)+1から
COLORindex(icol)までの要素がicol番
目の色（レベル）に含まれる。
I 新番号⇒旧番号への参照配列

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

PEsmpTOTI スレッド数

SMPindex(0:NCOLORtot*PEsmpTOT)

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

COLORindex(0:NCOLORtot)

SMPindexG(0:PEsmpTOT)OMP-321

プログラムの構成
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-322

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
24
100• PEsmpTOT– OpenMPスレッド数• NCOLORtot– 色数– 「＝0」の場合はCM– 「＝-1」の場合はRCM– 「≦-2」の場合はCMRCM
NX/NY/NZ
DX/DY/DZEPSICCGPEsmpTOTNCOLORtot

OMP-323

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-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

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.d0
D = 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)26

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)27

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-328

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-329

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===

SMPindex:
for preconditioning
各色内の要素数：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-330

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スレッドの例• 同じ「色」に属する要素は独立⇒並列計算可能• 色の順番に並び替え

31

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) …

Color (ic)

Thread (ip)

Numberingcolor=1

1 2 3 4 5 6 7 8color=2

1 2 3 4 5 6 7 8color=3

1 2 3 4 5 6 7 8color=4

1 2 3 4 5 6 7 8color=5

1 2 3 4 5 6 7 8

ip1= (ic-1)*PEsmpTOT + ip

ip1+1

ParallelAccessing

32

SMPindex
COLORindex(ic )= 100
COLORindex(ic+1)= 200
PEsmpTOT= 8
nn1= 200-100 = 100
num= 100 / 8 = 12 (12.5)
nr = 100 – 12*8=
4
ip0= (ic-1)*PEsmpTOT (ic: starting at 1)
SMPindex(ip0 )= 100
SMPindex(ip0+1)= 113 (13 elements in the 1st thread)
SMPindex(ip0+2)= 126 (13 elements in the 2nd thread)
SMPindex(ip0+3)= 139 (13 elements in the 3rd thread)
SMPindex(ip0+4)= 152 (13 elements in the 4th thread)
SMPindex(ip0+5)= 164 (12 elements in the 5th thread)
SMPindex(ip0+6)= 176 (12 elements in the 6th thread)
SMPindex(ip0+7)= 188 (12 elements in the 7th thread)
SMPindex(ip0+8)= 200 (12 elements in the 8th thread)OMP-333

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-334

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

35!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

poi_gen（5/9）
これ以降は新しい
番号付けを使用

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

配列の宣言

OMP-336!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-337係数の計算：並列に実施可能
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-338!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-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)
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)(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∆y

RDZ =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-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∆y

RDZ =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

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-343

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-344!$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-345

プログラムの構成
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-346

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-347

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-348

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

−1



2
−1
d i =  aii −aik ⋅ d k  = lii
k =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-349

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

−1



2
−1
d i =  aii −aik ⋅ d k  = lii
k =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-350

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-351

行列ベクトル積
依存性が無い⇒独立に計算可能⇒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-352

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-353

内積：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-354

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-355

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-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）(M ){z} = (LDL ){z} = {r}T(L ){z} = {r}前進代入
Forward Substitution(DL ){z} = {z}T

後退代入
Backward SubstitutionOMP-357

前進代入：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===

58

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===

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-360!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-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-3• L2-solへのOpenMPの実装• 実行例• 最適化＋演習

63

OMP-364

計算結果• 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+0365

66

FX10, 1-node/16-cores, 1003（●：MC，△：RCM，-：CM-RCM）
Iterations

Iterations

350

300

MC
RCM

250

CM-RCM

Number of Incompatible Nodes

400

200

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

1.E+02RCM

ICP#

CM-RCM

1.E+01
1.E+00

1

10

12.00

100

1000

10000

1

10.00

MC

Time (solver)RCM

sec./iteration

6.00
4.00
2.00

100

1000

10000

COLOR#

RCM

3.00E-02

CM-RCM

8.00

10

4.00E-02

COLOR#
MC

sec.

1.E+06

CM-RCMTime/
Iteration

2.00E-02

1.00E-02

0.00E+00

0.00
1

10

100

COLOR#

1000

10000

1

10

100

COLOR#

1000

10000

67
500

450

Iterations

OBCX,
1-soket/24-cores,
1283

MC
RCM
CM-RCM

Iterations

400

350

（●：MC，△：RCM，-：CMRCM）

300

250
1

10

100

1000

Color

5.00

MC

RCM

1.50E-02

CM-RCMMC

RCM

CM-RCM

4.00

sec./iterationsec.

1.00E-02

3.00

2.00

5.00E-03

1.00

Time/Iteration

Time (solver)
0.00E+00

0.00
1

10

100

Color

1000

1

10

100

Color

1000

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

68

OMP-E69• マルチコア版コードの実行• 更なる最適化

OMP-1

コンパイル・実行>$ cd /work/gt00/t00XYZ>$ cd multicore/L3/src>$ make>$ ls ../run/L3-sol
L3-sol>$ cd ../run>$ qsub go1.sh

70OMP-371

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

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

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

OMP-372

制御データ（INPUT.DAT）
128 128 128
1.00e-00 1.00e-00 1.00e-00
1.0e-08
24-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-373

go1.sh#!/bin/sh#PJM -N "test"#PJM -L rscgrp=tutorial#PJM -L node=1#PJM --omp thread=24(= PEsmpTOT)#PJM -L elapse=00:15:00#PJM -g gt00
#PJM -j#PJM -e err#PJM -o test1.lst
export KMP_AFFINITY=granularity=fine,compact./L3-sol

74
500

450

Iterations

OBCX,
1-soket/24-cores,
1283

MC
RCM
CM-RCM

Iterations

400

350

（●：MC，△：RCM，-：CMRCM）

300

250
1

10

100

1000

Color

5.00

MC

RCM

1.50E-02

CM-RCMMC

RCM

CM-RCM

4.00

sec./iterationsec.

1.00E-02

3.00

2.00

5.00E-03

1.00

Time/Iteration

Time (solver)
0.00E+00

0.00
1

10

100

Color

1000

1

10

100

Color

1000

OMP-E75• マルチコア版コードの実行• 更なる最適化– その1： OpenMP Statement– その2： Sequential ReorderingOMP-E76

前進代入：現状の並列化（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」でスレッド（～28/56）の生成，消滅が発生– 色ごとにこの部分を通る– 多少のオーバーヘッドがある• 色数が増えるとオーバーヘッドが増すOMP-E77

前進代入：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

78

プログラム類（src0）% cd /work/gt00/t00XYZ% cd multicore/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

srcx

OMP-379

go0.sh#!/bin/sh#PJM -N "test"#PJM -L rscgrp=tutorial#PJM -L node=1#PJM --omp thread=24(= PEsmpTOT)#PJM -L elapse=00:15:00#PJM -g gt00
#PJM -j#PJM -e err#PJM -o test1.lst
export KMP_AFFINITY=granularity=fine,compact./L3-sol0

OMP-E80

プログラム類（srcx）% cd /work/gt00/t00XYZ% cd multicore/L3% ls
run reorder0 src% cd srcx% make% cd ../run% ls L3-solx
L3-solx<modify “INPUT.DAT”><modify “gox.sh”>% pjsub gox.sh

src0

srcx

OMP-381

gox.sh#!/bin/sh#PJM -N "test"#PJM -L rscgrp=tutorial#PJM -L node=1#PJM --omp thread=24(= PEsmpTOT)#PJM -L elapse=00:15:00#PJM -g gt00
#PJM -j#PJM -e err#PJM -o test1.lst
export KMP_AFFINITY=granularity=fine,compact./L3-solx

OMP-382

private(ic,ip,ip1,i,VAL,k)src0とsrcxの !$omp doparallel
ic= 1, NCOLORtot!$omp dodo ip= 1, PEsmpTOT違い（1/2）
ip1= (ic-1)*PEsmpTOT + ip
do i= SMPindex(ip1-1)+1, SMPindex(ip1)
VAL= D(i)IC分解do k= indexL(i-1)+1, indexL(i)srcxはL2-solにOpenMPをそのまま適用したのに近い
両者の性能差はほとんどない

VAL= VAL - (AL(k)**2) * W(itemL(k),DD)
enddo
W(i,DD)= 1.d0/VAL
enddo
enddo
enddo
!$omp end parallel

src0!$omp parallel private(ic,ip1,ip2,i,VAL,k)do ic= 1, NCOLORtot
ip1= SMPindex((ic-1)*PEsmpTOT) + 1
ip2= SMPindex((ic-1)*PEsmpTOT + PEsmpTOT)!$omp dodo i= ip1, ip2
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 parallelsrcx

OMP-383

private(ic,ip,ip1,i,VAL,k)src0とsrcxの !$omp doparallel
ic= 1, NCOLORtot!$omp dodo ip= 1, PEsmpTOT違い（1a/2）
ip1= (ic-1)*PEsmpTOT + ip
do i= SMPindex(ip1-1)+1, SMPindex(ip1)
VAL= D(i)IC分解do k= indexL(i-1)+1, indexL(i)srcxはL2-solにOpenMPをそのまま適用したのに近い
両者の性能差はほとんどない

VAL= VAL - (AL(k)**2) * W(itemL(k),DD)
enddo
W(i,DD)= 1.d0/VAL
enddo
enddo
enddo
!$omp end parallel

src0!$omp parallel private(ic,ip1,ip2,i,VAL,k)do ic= 1, NCOLORtot!$omp dodo i= COLORindex(ic-1)+1, COLORindex(ic)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 parallelsrcx

OMP-384

do private(ip,i,VAL,k)src0とsrcxの !$omp parallel
do ip= 1, PEsmpTOT
do i = SMPindexG(ip-1)+1, SMPindexG(ip)
VAL= D(i)*W(i,P)
違い（2/2）
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)srcxはL2-solにOpenMPをそのまま適用したのに近い
両者の性能差はほとんどない

VAL= VAL + AU(k)*W(itemU(k),P)
enddo
W(i,Q)= VAL
enddo
enddo

src0

!$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

srcx

OMP-385

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

Iterations

Iterations

400

350

300

250
1

10

100

Color

1000

OMP-386

ICCG法計算時間（24スレッド）：CM-RCM
色数が増えるとsrcが（やや）遅い：スレッド生成・消滅の影響
src0/srcxの差はほぼ無い，色数多いと全体的に不安定
2.40

src

2.30

src0

srcx

2.20

sec.

2.10
2.00
1.90
1.80
1.70
1.60
1

10

100

Color

1000

OMP-E87• マルチコア版コードの実行• 更なる最適化– その1： OpenMP Statement– その2： Sequential ReorderingOMP-388

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

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-390

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

データ再配置: 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

92

番号付けがinconsistentになる可能性がある• Coalesced– icel1 > icel2, therefore, icel2= itemL[k], whereindexL[icel1]≦k＜indexL[icel1+1]• Sequential– icel1n < icel2n, but still icel2n= itemL[k],where indexL[icel1n]≦k＜indexL[icel1n+2]
color=1

color=2

color=3

color=4

color=512345678 12345678 12345678 12345678 12345678
icel2

icel1

11111 22222 33333 44444 55555 66666 77777 88888
icel1nicel2nOMP-393

データ再配置：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-394

データ再配置：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-395

データ再配置：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-396

データ再配置: Sequential Reordering
Initial Vector

Coalesced
GPUはこちらがお勧め

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

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

Sequential

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-E

97

プログラム類（reorder0）% cd /work/gt00/t00XYZ% cd multicore/L3% ls
run reorder0 src% cd reorder0% make% cd ../run% ls L3-rsol0
L3-rsol0<modify “INPUT.DAT”><modify “gor.sh”>% pjsub gor.sh

src0

srcx

OMP-398

gor.sh#!/bin/sh#PJM -N "test"#PJM -L rscgrp=tutorial#PJM -L node=1#PJM --omp thread=24(= PEsmpTOT)#PJM -L elapse=00:15:00#PJM -g gt00
#PJM -j#PJM -e err#PJM -o test1.lst
export KMP_AFFINITY=granularity=fine,compact./L3-rsol0

OMP-399

制御データ（INPUT.DAT）
128 128 128
1.00e-00 1.00e-00 1.00e-00
1.0e-08
24
-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：前進後退代入と同じ）

100

program MAIN
use STRUCT
use PCGuse solver_ICCG_mc
use solver_ICCG_mc_ft
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

Sequential
Reordering
(1/5)
Main

ISET= 0allocate (WK(ICELTOT))

&
&
&
&
&
&

if (METHOD.eq.0) thencall solve_ICCG_mc( ICELTOT, NPL, NPU, indexL, itemL, indexU, itemU,
D, BFORCE, PHI, AL, AU, NCOLORtot, PEsmpTOT,
SMPindex_new, EPSICCG, ITR, IER)
else
call solve_ICCG_mc_ft( ICELTOT, NPL, NPU, indexL, itemL, indexU, itemU,
D, BFORCE, PHI, AL, AU, NCOLORtot, PEsmpTOT,
SMPindex_new, EPSICCG, ITR, IER)
endif

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

&
&
&
&
&
&

101

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

Coalesced
ic=3

ic=4

ic=5

12345678 12345678 12345678 12345678 12345678

SMPindex_new

Sequential

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

Sequential
Reordering(2/5)
poi_gen-1

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

102

Coalesced!$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) …

Color (ic)

Thread (ip)

Numberingcolor=1

1 2 3 4 5 6 7 8color=2

1 2 3 4 5 6 7 8color=3

1 2 3 4 5 6 7 8color=4

1 2 3 4 5 6 7 8color=5

1 2 3 4 5 6 7 8

ip1= (ic-1)*PEsmpTOT + ip

ip1+1

ParallelAccessing103

Sequential!$omp parallel private(ic,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) …

Thread (ip)

Color (ic)ip1= (ip-1)*NCOLORtot + ic

thread #1

1 2 3 4 5

thread #2

1 2 3 4 5

thread #3

1 2 3 4 5

thread #4

1 2 3 4 5

thread #5

1 2 3 4 5

thread #6

1 2 3 4 5

thread #7

1 2 3 4 5

thread #8

1 2 3 4 5

ip1+1
Numbering

ParallelAccessing104

Sequential
Reordering (3/5)
poi_gen-2

do ip= 1, PEsmpTOT
do ic= 1, NCOLORtot
icNS= SMPindex_new((ip-1)*NCOLORtot + ic-1)
ic01= SMPindex((ic-1)*PEsmpTOT + ip-1) + 1
ic02= SMPindex((ic-1)*PEsmpTOT + ip )
icou= 0do k= ic01, ic02
icel= NEWtoOLD(k)icou= icou + 1
icelN= icNS + icou
OLDtoNEWnew: Original -> Sequential
OLDtoNEWnew(icel )= icelN
NEWtoOLDnew(icelN)= icel
NEWtoOLDnew: Sequential -> Original
enddo-Original: Initial icel
enddo
enddo
-Sequential
icelN!$omp parallel do private (ip,icel,ic0,ik0)
do ip = 1, PEsmpTOT
do icel= SMPindex_new((ip-1)*NCOLORtot)+1, SMPindex_new(ip*NCOLORtot)
ic0 = NEWtoOLDnew(icel)
ik0 = OLDtoNEW(ic0)indexL(icel)= INL(ik0)indexU(icel)= INU(ik0)

Sequential

Coalesced

enddo
enddo
do icel= 1, ICELTOT
indexL(icel)= indexL(icel) + indexL(icel-1)
indexU(icel)= indexU(icel) + indexU(icel-1)
enddo-Original: Initial ic0-Coalesced
ik0-Sequential
icel

105

Sequential Reordering (4/5)
poi_gen-3
NPL= indexL(ICELTOT)
NPU= indexU(ICELTOT)
allocate (itemL(NPL), AL(NPL))
allocate (itemU(NPU), AU(NPU))
if (NFLAG.eq.0) thenitemL= 0
itemU= 0
AL= 0.d0
AU= 0.d0
else!$omp parallel do private (ip,icel,k)do ip= 1, PEsmpTOT
do icel= SMPindex_new((ip-1)*NCOLORtot+1),&
SMPindex_new(ip*NCOLORtot)
do k= indexL(icel-1)+1, indexL(icel)
itemL(k)= 0
AL(k)= 0.d0
enddo
do k= indexU(icel-1)+1, indexU(icel)
itemU(k)= 0
AU(k)= 0.d0
enddo
enddo
enddo
!$omp end parallel do
endif

!$omp parallel do private (ip,icel,ic0,icN1,icN2,icN3,icN4, icN5,icN6) &!$omp&
private (coef,j,ii,jj,kk)&!$omp&
private (ik0,icN10,icN20,icN30,icN40, icN50,icN60)
do ip = 1, PEsmpTOT
do icel= SMPindex_new((ip-1)*NCOLORtot)+1, SMPindex_new(ip*NCOLORtot)
ic0 = NEWtoOLDnew(icel)
ik0 = OLDtoNEW(ic0)
icN10= NEIBcell(ic0,1)
icN20= NEIBcell(ic0,2)
icN30= NEIBcell(ic0,3)
icN40= NEIBcell(ic0,4)
icN50= NEIBcell(ic0,5)
icN60= NEIBcell(ic0,6)
if (icN50.ne.0) then
icN5= OLDtoNEW(icN50)
coef= RDZ * ZAREA
D(icel)= D(icel) - coef
if (icN5.lt.ik0) thendo j= 1, INL(ik0)
if (IAL(j,ik0).eq.icN5) thenitemL(j+indexL(icel-1))= OLDtoNEWnew(icN50)
AL(j+indexL(icel-1))= coef
exit
endif
enddo
else
do j= 1, INU(ik0)
if (IAU(j,ik0).eq.icN5) thenitemU(j+indexU(icel-1))= OLDtoNEWnew(icN50)
AU(j+indexU(icel-1))= coef
exit
endif
enddo
endif
endif

Sequential
Reordering(5/5)
poi_gen-4
icel: Sequential
ic0: Original
ik0: Coalesced
icN50: Original
icN5 : Coalesced
icN5>ik0: Upper (AU)
icN5<ik0: Lower (AL)106

前進代入の計算法：色ループは外Color #2 Thread #1-#(Pe)Color #3 Thread #1-#(Pe)Color #4 Thread #1-#(Pe)Color #Nc Thread #1-#(Pe)Coalesced
Thread #1

Color #1-#(Nc)Thread #2

Color #1-#(Nc)Thread #3

Color #1-#(Nc)Thread #4

Color #1-#(Nc)Thread#(Pe)

Color #1-#(Nc)・・・!$omp parallel private(ic,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

Color #1 Thread #1-#(Pe)・・・!$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
enddo
!$omp end parallel

Sequential

107

108

行列ベクトル積!$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-3109

ICCG法計算時間：CM-RCM
Coalesced（src,src0）と比較してSequential（reorder0）が安定で速い（特に色数が大きい場合）（24-threads）
2.40

src

2.30

src0

reorder0

2.20

sec.

2.10
2.00
1.90
1.80
1.70
1.60
1

10

100

Color

1000

OMP-3

110

Results on OBCX： 1283
24 cores, src: MC（2-colors）：424 iter’s，2.287 sec.
RCM（382-levels）：287 iter’s，2.117 sec.
CM-RCM（Nc=20）：318 iter’s，1.830 sec.
Computation Time

Speed-Up based on RCM (src)with 1 thread

25.00

15.00

MC(2)
RCM
CM-RCM(20)

20.00

MC(2)
RCM
CM-RCM(20)
Ideal

Speed-Up

sec.

10.00
15.00

10.00

5.00
5.00

0.00

0.00
0

4

8

12

16

Color

20

24

28

0

4

8

12

16

Color

20

24

28

OMP-3

111

Results on OBCX： 1283
24 cores, reorder0：MC（2-colors）：424 iter’s，2.287 sec.
RCM（382-levels）：287 iter’s，1.857 sec.
CM-RCM（Nc=20）：318 iter’s，1.754 sec.
Computation Time

Speed-Up based on RCM (src)with 1 thread

25.00

15.00

MC(2)
RCM
CM-RCM(20)

20.00

MC(2)
RCM
CM-RCM(20)
Ideal

Speed-Up

sec.

10.00
15.00

10.00

5.00
5.00

0.00

0.00
0

4

8

12

16

Color

20

24

28

0

4

8

12

16

Color

20

24

28

112

OBCX：1ノード：2CPU（ソケット）
Memory
96 GB

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

2933 MHz×6ch
140.8 GB/sec

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #0: 0th-27th cores

UPI
UPI
UPI

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #1: 28th-55th cores

Ultra Path Interconnect
10.4 GT/sec ×3= 124.8 GB/sec

Socket #0

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

Memory96 GB2933 MHz×6ch
140.8 GB/sec

Socket #1• ここまでは1CPU（1ソケット）を使ってきた• 1ノード（2CPU，56コア）を使ってみよう• NUMAアーキテクチャ（Non-Uniform Memory Access）• メモリは各CPUに搭載されていて独立，異なるCPUのローカルメモリ上のデータをアクセスすることは可能• できるだけローカルメモリ上のデータを使って計算するのが効率的

OMP-3113

データをローカルメモリに置く方法（1/2）• NUMAアーキテクチャでは，プログラムにおいて変数や配列を宣言した時点では，物理的メモリ上に記憶領域は確保されず，ある変数を最初にアクセスしたコア（の属するソケット）のローカルメモリ上に，その変数の記憶領域（ページ）が確保される。
• これをFirst Touch Data Placement（First Touch）と呼び，配列の初期化手順により得られる性能が大幅に変化する場合があるため，注意が必要である。
• 例えばある配列を初期化する場合，特に指定しなければ0
番のソケットで初期化が行われるため，記憶領域は0番ソ
ケットのローカルメモリ上に確保される。

OMP-3114

データをローカルメモリに置く方法（2/2）• したがって，他のソケットでこの配列のデータをアクセスする場合には，必ず0番ソケットのメモリにアクセスする必要があるため，高い性能を得ることは困難である。
• 配列の初期化を，実際の計算の手順にしたがってOpenMPを使って並列に実施すれば，実際に計算を担当
するソケットのメモリにその配列の担当部分の記憶領域が確保され，より効率的に計算を実施することができる。
• 1ソケットしか使用しない場合はこのような配慮は不要– OpenMP/MPIハイブリッドで1CPU（ソケット）当りに1つのMPIプ
ロセスを使用する場合も同様– numactl

OMP-3

115

First Touch Data Placement“Patterns for Parallel Programming” Mattson, T.G. et al.
To reduce memory traffic in the system, it is important to keep thedata close to the PEs that will work with the data (e.g. NUMA control).
On NUMA computers, this corresponds to making sure the pages ofmemory are allocated and “owned” by the PEs that will be working
with the data contained in the page.
The most common NUMA page-placement algorithm is the“first touch” algorithm, in which the PE first referencing a region ofmemory will have the page holding that memory assigned to it.
A very common technique in OpenMP program is to initialize
data in parallel using the same loop schedule as will be used laterin the computations.

116

Local/Remote Memory

L2
L1

MemoryMemoryMemoryMemory

L3

L3

L3

L3

L2
L1

L2
L1

L2
L1

L2
L1

L2
L1

L2
L1

L2
L1

L2
L1

L2
L1

L2
L1

L2
L1

L2
L1

L2
L1

L2
L1

L2
L1

Core Core Core Core

Core Core Core Core

Core Core Core Core

Core Core Core Core

Core Core Core Core

Core Core Core Core

Core Core Core Core

Core Core Core Core

L1
L2

L1
L2

L1
L2

L1
L2

L1
L2

L1
L2

L1
L2

L1
L2

L1
L2

L1
L2

L1
L2

L1
L2

L1
L2

L1
L2

L1
L2

L3

L3

L3

L3

MemoryMemoryMemoryMemoryLocal Memory

Remote Memory

L1
L2

OMP-3117

制御データ（INPUT.DAT）
128 128 128
1.00e-00 1.00e-00 1.00e-00
1.0e-0848-50
0
0

NX/NY/NZ
DX/DY/DZEPSICC
PEsmpTOTNCOLORtot
NFLAG (0 or 1)
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-3118

go2.sh: reorder0#!/bin/sh#PJM -N "test"#PJM -L rscgrp=tutorial#PJM -L node=1#PJM --omp thread=48(= PEsmpTOT)#PJM -L elapse=00:15:00#PJM -g gt00
#PJM -j#PJM -e err#PJM -o test1.lst
export KMP_AFFINITY=granularity=fine,compact./L3-rsol0

OMP-3

119

Array Initialization: NFLAG=0/1 (1/3)
poi_gen.f!C!C-- ARRAY init.
Pages are allocated at theif (NFLAG.eq.0) then
BFORCE= 0.d0
local memory of the master
PHI= 0.d0threadD= 0.d0
OLDtoNEWnew= 0
NEWtoOLDnew= 0
else!$omp parallel do private (ip,icel)
do ip= 1, PEsmpTOT
do icel= SMPindex_new((ip-1)*NCOLORtot+1),&
SMPindex_new(ip*NCOLORtot)
BFORCE(icel)= 0.d0
PHI
(icel)= 0.d0
Pages are allocated at theD(icel)= 0.d0local memory of each thread
OLDtoNEWnew(icel)= 0
NEWtoOLDnew(icel)= 0
enddo
enddo
!$omp end parallel doendif

A very common technique in OpenMP program for optimization is toinitialize data in parallel using the same loop schedule as will be usedlater in the computations.

OMP-3120

Array Initialization: NFLAG=0/1 (2/3)
if (NFLAG.eq.0) thendo icel= 1, ICELTOT
Pages are allocated at theindexL(icel)= 0
local memory of the masterindexU(icel)= 0
enddo
threadelse!$omp parallel do private (ip,icel,k)do ip= 1, PEsmpTOT
do icel= SMPindex_new((ip-1)*NCOLORtot+1),&
SMPindex_new(ip*NCOLORtot)
indexL(icel)= 0
indexU(icel)= 0
Pages are allocated at theenddo
local memory of each thread
enddo
endif

A very common technique in OpenMP program for optimization is toinitialize data in parallel using the same loop schedule as will be usedlater in the computations.

OMP-3

121

Array Initialization: NFLAG=0/1 (3/3)
if (NFLAG.eq.0) thenitemL= 0
Pages are allocated at the
itemU= 0
local memory of the master
AL= 0.d0
AU= 0.d0
threadelse!$omp parallel do private (ip,icel,k)do ip= 1, PEsmpTOT
do icel= SMPindex_new((ip-1)*NCOLORtot+1),&
SMPindex_new(ip*NCOLORtot)
do k= indexL(icel-1)+1, indexL(icel)
itemL(k)= 0
AL(k)= 0.d0
enddo
do k= indexU(icel-1)+1, indexU(icel)
itemU(k)= 0
AU(k)= 0.d0
enddo
Pages are allocated at theenddo
local memory of each thread
enddo!$omp end parallel doendif

A very common technique in OpenMP program for optimization is toinitialize data in parallel using the same loop schedule as will be usedlater in the computations.

122

Sequential Reordering (4/5)
poi_gen-3
NPL= indexL(ICELTOT)
NPU= indexU(ICELTOT)
allocate (itemL(NPL), AL(NPL))
allocate (itemU(NPU), AU(NPU))
if (NFLAG.eq.0) thenitemL= 0
itemU= 0
Pages are allocated at the
AL= 0.d0
local memory of the master
AU= 0.d0
else
thread!$omp parallel do private (ip,icel,k)do ip= 1, PEsmpTOT
do icel= SMPindex_new((ip-1)*NCOLORtot+1),&
SMPindex_new(ip*NCOLORtot)
do k= indexL(icel-1)+1, indexL(icel)
itemL(k)= 0
AL(k)= 0.d0
Pages are allocated at theenddo
do k= indexU(icel-1)+1, indexU(icel)
local memory of each thread
itemU(k)= 0
AU(k)= 0.d0
enddo
enddo
enddo
!$omp end parallel do
endifOMP-3123

計算結果：reoder0, L3-rsol0
128 128 128
1.00e-0 1.00e-0 1.00e0
1.0e-0848-1
0
0
Thread #
24
28
48
56

NX/NY/NZ
DX/DY/DZOMEGA, EPSICCG
PEsmpTOTNCOLORtot
NFLAG
METHOD

NFLAG
0
0
0

RCM
1.736
1.859
1.550

CM-RCM(50)
1.670
1.764
1.275

1
0
1

1.274
1.681
1.321

0.891
1.297
0.819

OMP-3124

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

OMP-3125

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

ICCG法計算時間：CM-RCM, reorder0
全ての計算を単精度・倍精度で実施した場合の比較（倍精度⇒単精度）反復回数：20-25%増加，
計算時間：30-40%減少
1.40

Iterations

1.30

ICCG Time

Single/Double

1.20
1.10
1.00
0.90
0.80
0.70
0.60
0.50
1

10

100

Color

1000

OMP-3

127

この後• UIDは一ヶ月有効• 何かあったらいつでも気軽に質問してください• 利用方法等に関する質問は相談窓口ではなく中島へ• 一ヶ月間に限りZoomでの個別のセッション（一人最大45分～60分程度）にも対応しますので，ご連絡ください。
• プログラム類は研究等にご利用いただいて結構です。
商用プログラムの一部になる場合，これを元にしたプログラムでビジネスをやる場合のみご相談ください。
• アンケートもよろしく