科学技術計算のためのマルチコアプログラミング入門
第Ⅲ部：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

プログラムのありか on Oakleaf-FX• 所在– <$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-37

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

src0

src

% cd src% make% cd ../run% ls L3-sol
L3-sol% <modify “INPUT.DAT”>% <modify “go1.sh”>% pjsub go1.shshort (Fortranのみ)

OMP-38

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

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

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

OMP-39制御データ（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-310

go1.sh#!/bin/sh#PJM -L "node=1"#PJM -L "elapse=00:10:00"#PJM -L "rscgrp=lecture"#PJM -g "gt61"#PJM -j#PJM -o "arcm.lst"
export OMP_NUM_THREADS=16./L3-solPEsmpTOTと一致させるOMP-311• L2-solへのOpenMPの実装• Hitachi SR11000/J2での実行– 計算結果– CM-RCMオーダリング• FX10での実行

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• スレッド数をプログラムで制御できるようにしてみよう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• スレッド数をプログラムで制御できるようにしてみようOMP-315

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

OMP-3

16

プログラムの構成（1/2）
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-317プログラムの構成（2/2）
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-318

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

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

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

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

OMP-322

変数表（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-323

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

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

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

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

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

配列の宣言

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

IAL(j,i):
AL(j,i):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

IAL(j,i):
AL(j,i):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-364• L2-solへのOpenMPの実装• Hitachi SR11000/J2での実行– 計算結果– CM-RCMオーダリング• FX10での実行

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

計算時間（MC=2）
以降，「linear solver」の計算時間のみ問題とする
matrix assembling： poi-genの後半：並列化
linear solver：CG

30.0
matrix assebmling
linear solver

sec.

20.0

10.0

0.0
1

2

4
thread#816

OMP-367

スケーラビリティ
反復回数：MC（2色）：333回，RCM（298レベル）：224回
20.0

30.0

MC=2

MC=2
RCM(298)sec.

Speed-Up

20.0

RCM(298)

16.0

12.0

8.0

10.0

4.0

0.0

0.0
1

2

4
thread#

8

16

1

2

4
thread#• 反復回数ではRCMの方が少ないのに，コア数が増えると，計算時間が逆転する• MC=2は良好なスケーラビリティ（16コアで15.5）

8

16

OMP-368

原 因• 色数が多い（298対2）⇒同期オーバーヘッド• スレッドごとの負荷不均衡– RCMでは要素数の少ない「レベル」が必ず存在する
29 22 16 11

7

4

2

1

37 30 23 17 12

8

5

3

44 38 31 24 18 139650 45 39 32 25 19 14 1055 51 46 40 33 26 20 1559 56 52 47 41 34 27 21
62 60 57 53 48 42 35 28
64 63 61 58 54 49 43 36

OMP-369

原因はこのようなループにある
色数増加⇒同期オーバーヘッド増加
必ずある「色」の計算が終わってから次の「色」に行く
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-370

MCとRCMの比較• MC– 並列性高い，負荷分散も良い（そのように設定されている）– 特に色数少ないと反復回数多い– 色数を増やす，反復回数減るが同期オーバーヘッドの影響• RCM– 収束は早いが，レベル数が多く，同期オーバーヘッドの影響受けやすく，コア数が増えると不利– 負荷分散もいま一つ• 反復回数少なくて同期オーバーヘッドの影響が少ない
方法が無いものか?
– 色数が少なくて，かつ反復回数が少ないという都合の良い方法

OMP-371

16コアにおける結果（●：MC，△：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.

sec./iteration

2.50

2.00

1.50
1.E+00

1.E+01

1.E+02COLOR#

1.E+038.0E-03

6.0E-03

4.0E-03
1.E+00

1.E+01

1.E+02COLOR#

1.E+03

OMP-372

解決策：CM-RCM• RCM＋Cyclic Mulricoloring 〔土肥，襲田，鷲尾他〕• 手順– まずRCMを施す– Cyclic Multicoloring（CM）の色数を決める（Nc）– RCMの1番目，（Nc+1）番目，（2Nc+1）番目・・・のレベルに属する要素を「1」色に分類する– RCMのk番目，（Nc+k）番目，（2Nc+k）番目・・・のレベルに属する要素を「k」色に分類する– 「k」が「Nc」に達して，要素が「1~Nc」で色付けされたら完了• あとはMCのときと同じように，色の順番に再番号付• RCMの各レベルに対して「Nc」のサイクルで再色付けを実施している– もし同じ色の要素の中に依存性が見つかったら，Nc=Nc+1として最初からやり直し（ここは少し原始的）OMP-373

RCM
29 22 16 11

7

4

2

1 1

9 37

30 23 17 12

8

5

3 2

10 44

38 31 24 18 13

9

6 3

11 50

45 39 32 25 19 14 10 4

12 55

51 46 40 33 26 20 15 5

13 59

56 52 47 41 34 27 21 6

14 62

60 57 53 48 42 35 28 7

15 64

63 61 58 54 49 43 36 8

OMP-374

Nc=4, k=1：1,5,9,13レベルを選択
29 22 16 11

7

4

2

1 1

9 37

30 23 17 12

8

5

3 2

10 44

38 31 24 18 13

9

6 3

11 50

45 39 32 25 19 14 10 4

12 55

51 46 40 33 26 20 15 5

13 59

56 52 47 41 34 27 21 6

14 62

60 57 53 48 42 35 28 7

15 64

63 61 58 54 49 43 36 8

OMP-375

Nc=4, k=2：2,6,10,14レベルを選択
29 22 16 11

7

4

2

1 1

9 37

30 23 17 12

8

5

3 2

10 44

38 31 24 18 13

9

6 3

11 50

45 39 32 25 19 14 10 4

12 55

51 46 40 33 26 20 15 5

13 59

56 52 47 41 34 27 21 6

14 62

60 57 53 48 42 35 28 7

15 64

63 61 58 54 49 43 36 8

OMP-376

CM-RCM（Nc=4）：「色」の順番に再並替
53 36 197

54 37 20

25

8

43 26

k=1: 16k=2: 16k=3: 16
k=4: 162

49 33 173

55 38 211

50 34 184552

61 44 27 10 57 40 2369

56 39 2251 35

14 62 45 28 11 58 41 24
31 15 63 46 29 12 59 42
48 32 16 64 47 30 13 60

OMP-377

16コアにおける結果（●：MC，△：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.

sec./iteration

2.50

2.00

1.50
1.E+00

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

1.E+038.0E-03

6.0E-03

4.0E-03
1.E+00

1.E+01

1.E+02COLOR#

1.E+03

OMP-378

スケーラビリティ
反復回数：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
12.0

8.0

10.0
4.0

0.0

0.0
1

2

4

8

16

thread#16 threads
MC(2):
1.83 sec.
CM-RCM(20): 1.79 sec.

1

2

4
thread#816

OMP-379

CM-RCM• 少ない色数（Nc）で良好な収束を得られる• 計算効率も良い• 実行方法– INPUT.DATで「NCOLORtot=-Nc」とする– L2，L3で有効なオプション（既にL2でも使えていた）• 実装は「cmrcm.f」を参照くださいOMP-380• L2-solへのOpenMPの実装• Hitachi SR11000/J2での実行– 計算結果– CM-RCMオーダリング• FX10での実行

OMP-381

計算結果（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

OMP-382

計算結果（SR11000）：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
12.0

8.0

10.0
4.0

0.0

0.0
1

2

4

8

16

thread#16 threads
MC(2):
1.83 sec.
CM-RCM(20): 1.79 sec.

1

2

4
thread#816

OMP-383

計算結果（T2K東大）：106要素
反復回数：MC（2色）：333回，RCM（298レベル）：224回
CM-RCM（Nc=20）：249回
ccNUMA向け最適化
50.0

20.0
MC=2

MC=2

RCM(298)

40.0

16.0

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

Speed-Up

30.0

20.0

10.0

RCM(298)

12.0

8.0

4.0

0.0

0.0
1

2

4

8

16

thread#16 threads
MC(2):
5.48 sec.
CM-RCM(20): 5.71 sec.

1

2

4

8

16

Memorythread#Memory

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

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

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

OMP-386

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