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

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

OMP-3

1

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

OMP-32

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

OMP-3

3

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

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

OMP-35

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

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

OMP-3

6

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

7

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

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

src0

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

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

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

10

Control Data: 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法

各方向の要素数11

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

または “rscgrp=school”

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

12

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

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

export OMP_NUM_THREADS=16./L3-sol

実行ファイル名

13

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

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

JOB_NAME
go.sh

STATUS PROJECTRUNNING gt61

RSCGROUP
lectureSTART_DATE
08/27 12:58:08

ELAPSE
00:00:05

TOKEN NODE:COORD
0.0 1

14

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

15

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

実はこのようにしてDirectiveを直接挿入しても良いのだが・・・（1/2）#pragma omp parallel for private(i,VAL,j)for(i=0; i<N; i++) {VAL = D[i] * W[P][i];for(j=indexL[i]; j<indexL[i+1]; j++) {VAL += AL[j] * W[P][itemL[j]-1];}
for(j=indexU[i]; j<indexU[i+1]; j++) {VAL += AU[j] * W[P][itemU[j]-1];}
W[Q][i] = VAL;}

• スレッド数をプログラムで制御できるようにしてみよう• GPU，メニィコアではこのままの方が良い場合もある17

実はこのようにしてDirectiveを直接挿入しても良いのだが・・・（2/2）
for(ic=0; ic<NCOLORtot; ic++){#pragma omp parallel for private (i,WVAL,j)for(i=COLORindex[ic]; i<COLORindex[ic+1]; i++) {WVAL = W[Z][i];for(j=indexL[i]; j<indexL[i+1]; j++) {WVAL -= AL[j] * W[Z][itemL[j]-1];}
W[Z][i] = WVAL * W[DD][i];}}• スレッド数をプログラムで制御できるようにしてみよう• GPU，メニィコアではこのままの方が良い場合もある

18

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

Main Program#include <stdio.h> ...
intmain(){

double *WK;
int NPL, NPU; ISET, ITR, IER; icel, ic0, i;
double xN, xL, xU; Stime, Etime;if(INPUT()) goto error;
if(POINTER_INIT()) goto error;
if(BOUNDARY_CELL()) goto error;if(CELL_METRICS()) goto error;
if(POI_GEN()) goto error;

error:}ISET = 0;
WK = (double *)malloc(sizeof(double)*ICELTOT);if(WK == NULL) {fprintf(stderr, "Error: %s¥n", strerror(errno));goto error;}
Stime = omp_get_wtime();if(solve_ICCG_mc(ICELTOT, NL, NU, indexL, itemL, indexU, itemU,
D, BFORCE, PHI, AL, AU, NCOLORtot, PEsmpTOT,
SMPindex, SMPindexG, EPSICCG, &ITR, &IER)) goto error;
Etime = omp_get_wtime();for(ic0=0; ic0<ICELTOT; ic0++) {icel = NEWtoOLD[ic0];
WK[icel-1] = PHI[ic0];}
for(icel=0; icel<ICELTOT; icel++) {PHI[icel] = WK[icel];}if(OUTUCD()) goto error;
return 0;
return -1;

19

20

struct.h#ifndef __H_STRUCT#define __H_STRUCT#include <omp.h>

ICELTOT:
Number of meshes (NX x NY x NZ)

N:

int ICELTOT, ICELTOTp, N;
int NX, NY, NZ, NXP1, NYP1, NZP1, IBNODTOT;
int NXc, NYc, NZc;

Number of modes

double DX, DY, DZ, XAREA, YAREA, ZAREA;
double RDX, RDY, RDZ, RDX2, RDY2, RDZ2, R2DX, R2DY, R2DZ;
double *VOLCEL, *VOLNOD, *RVC, *RVN;

Number of meshes in x/y/z directions

int **XYZ, **NEIBcell;
int ZmaxCELtot;
int *BC_INDEX, *BC_NOD;
int *ZmaxCEL;
int **IWKX;
double **FCV;NX,NY,NZ:
NXP1,NYP1,NZP1:
Number of nodes in x/y/z directions

IBNODTOT:= NXP1 x NYP1XYZ[ICELTOT][3]:Location of meshes

int my_rank, PETOT, PEsmpTOT;

NEIBcell[ICELTOT][6]:#endif /* __H_STRUCT */Neighboring meshes

PEsmpTOT:
Number of threads

21#ifndef __H_PCG#define __H_PCG
static int N2 = 256;
int NUmax, NLmax, NCOLORtot, NCOLORk, NU,
NL;
int METHOD, ORDER_METHOD;
double EPSICCG;

pcg.h

double *D, *PHI, *BFORCE;
double *AL, *AU;int *INL, *INU, *COLORindex;
int *indexL, *indexU;
int *SMPindex, *SMPindexG;int *OLDtoNEW, *NEWtoOLD;
int **IAL, **IAU;
int *itemL, *itemU;
int NPL, NPU;
#endif /* __H_PCG */NCOLORtot
COLORindex[NCOLORtot+1]

Total number of colors/levels
Index of number of meshes in each color/level(COLORindex[icol+1]- COLORindex[icol])SMPindex [NCOLORtot*PEsmpTOT+1]
SMPindexG[PEsmpTOT+1]
OLDtoNEW, NEWtoOLD

Reference table before/after renumbering

22

Variables/Arrays for Matrix (1/2)
Name

Type

Content

D[N]R

Diagonal components of the matrix (N= ICELTOT)BFORCE[N]R

RHS vector

PHI[N]

R

Unknown vector

indexL[N+2]indexU[N+2]I# of L/U non-zero off-diag. comp. (CRS)

NPL, NPUI

Total # of L/U non-zero off-diag. comp. (CRS)

itemL[NPL]
itemU[NPU]

I

Column ID of L/U non-zero off-diag. comp. (CRS)

AL[NPL]
AU[NPU]R

L/U non-zero off-diag. comp. (CRS)

Name

Type

Content

NL,NU

I

MAX. # of L/U non-zero off-diag. comp. for each mesh (=6)

INL[N]INU[N]

I# of L/U non-zero off-diag. comp.

IAL[N][NL]
IAU[N][NU]

I

Column ID of L/U non-zero off-diag. comp.

23

Variables/Arrays for Matrix (2/2)
Name

Type

Content

NCOLORtot

I

Input:
reordering method + initial number ofcolors/levels≧2: MC, =0: CM, =-1: RCM, -2≧: CMRCM
Output:Final number of colors/levels

COLORindex[NCOLORtot+1]

I

Number of meshes at each color/level
1D compressed array
Meshes in icolth color/level are stored inthis array from COLORindex[icol]
to COLORindex[icol+1]-1

NEWtoOLD[N]I

Reference array from New to Old numbering

OLDtoNEW[N]I

Reference array from Old to New numbering

PEsmpTOTI

Number of Threads

SMPindex[NCOLORtot*PEsmpTOT+1]

I

Array for OpenMP Operations (for Loops with
Data Dependency)

SMPindexG[PEsmpTOT+1]

I

Array for OpenMP Operations (for Loops
without Data Dependency)

Main Program#include <stdio.h> ...
intmain(){

double *WK;
int NPL, NPU; ISET, ITR, IER; icel, ic0, i;
double xN, xL, xU; Stime, Etime;if(INPUT()) goto error;
if(POINTER_INIT()) goto error;
if(BOUNDARY_CELL()) goto error;if(CELL_METRICS()) goto error;
if(POI_GEN()) goto error;

error:}ISET = 0;
WK = (double *)malloc(sizeof(double)*ICELTOT);if(WK == NULL) {fprintf(stderr, "Error: %s¥n", strerror(errno));goto error;}
Stime = omp_get_wtime();if(solve_ICCG_mc(ICELTOT, NL, NU, indexL, itemL, indexU, itemU,
D, BFORCE, PHI, AL, AU, NCOLORtot, PEsmpTOT,
SMPindex, SMPindexG, EPSICCG, &ITR, &IER)) goto error;
Etime = omp_get_wtime();for(ic0=0; ic0<ICELTOT; ic0++) {icel = NEWtoOLD[ic0];
WK[icel-1] = PHI[ic0];}
for(icel=0; icel<ICELTOT; icel++) {PHI[icel] = WK[icel];}if(OUTUCD()) goto error;
return 0;
return -1;

24

25

input: reading INPUT.DAT#include <stdio.h>; <stdlib.h>; <string.h>; <errno.h>#include "struct_ext.h“; "pcg_ext.h“; "input.h"• PEsmpTOT

extern int
INPUT(void){#define BUF_SIZE 1024

– Thread Number

char line[BUF_SIZE];
char CNTFIL[81];double OMEGA;
FILE *fp11;• NCOLORtot

if((fp11 = fopen("INPUT.DAT", "r")) == NULL) {fprintf(stderr, "Error: %s¥n", strerror(errno));return -1;}
sscanf(line, "%d%d%d", &NX, &NY, &NZ);
sscanf(line, "%d", &METHOD);
sscanf(line, "%le%le%le", &DX, &DY, &DZ);
sscanf(line, "%le", &EPSICCG);
sscanf(line, "%d", &PEsmpTOT);
sscanf(line, "%d", &NCOLORtot);
fclose(fp11);
return 0;}

100 100 100
1.00e-02 5.00e-02 1.00e-02
1.00e-08
16
100

NX/NY/NZ
DX/DY/DZEPSICCGPEsmpTOTNCOLORtot– Reordering Method+ Initial Number ofColors/Levels– ≧2: MC
– =0: CM– =-1: RCM– -2≧: CMRCM

26

cell_metrics#include <stdio.h> ...
extern int
CELL_METRICS(void){
double V0, RV0;
int i;VOLCEL =(double *)allocate_vector(sizeof(double),ICELTOT);
RVC=(double *)allocate_vector(sizeof(double),ICELTOT);XAREA = DY * DZ;YAREA = DZ * DX;ZAREA = DX * DY;

XAREA

RDX = 1.0 / DX;
RDY = 1.0 / DY;
RDZ = 1.0 / DZ;
RDX2 = 1.0 / (pow(DX, 2.0));
RDY2 = 1.0 / (pow(DY, 2.0));
RDZ2 = 1.0 / (pow(DZ, 2.0));
R2DX = 1.0 / (0.5 * DX);
R2DY = 1.0 / (0.5 * DY);
R2DZ = 1.0 / (0.5 * DZ);
V0 = DX * DY * DZ;
RV0 = 1.0 / V0;for(i=0; i<ICELTOT; i++) {VOLCEL[i] = V0;
RVC[i]= RV0;}
return 0; }

DZ
DY
z

DX

y
x

Main Program#include <stdio.h> ...
intmain(){

double *WK;
int NPL, NPU; ISET, ITR, IER; icel, ic0, i;
double xN, xL, xU; Stime, Etime;if(INPUT()) goto error;
if(POINTER_INIT()) goto error;
if(BOUNDARY_CELL()) goto error;if(CELL_METRICS()) goto error;
if(POI_GEN()) goto error;

error:}ISET = 0;
WK = (double *)malloc(sizeof(double)*ICELTOT);if(WK == NULL) {fprintf(stderr, "Error: %s¥n", strerror(errno));goto error;}
Stime = omp_get_wtime();if(solve_ICCG_mc(ICELTOT, NL, NU, indexL, itemL, indexU, itemU,
D, BFORCE, PHI, AL, AU, NCOLORtot, PEsmpTOT,
SMPindex, SMPindexG, EPSICCG, &ITR, &IER)) goto error;
Etime = omp_get_wtime();for(ic0=0; ic0<ICELTOT; ic0++) {icel = NEWtoOLD[ic0];
WK[icel-1] = PHI[ic0];}
for(icel=0; icel<ICELTOT; icel++) {PHI[icel] = WK[icel];}if(OUTUCD()) goto error;
return 0;
return -1;

27

#include "allocate.h"extern int
POI_GEN(void){ int nn;
int ic0, icN1, icN2, icN3, icN4, icN5, icN6;
int i, j, k, ib, ic, ip, icel, icou, icol, icouG;
int ii, jj, kk, nn1, num, nr, j0, j1;
double coef, VOL0, S1t, E1t;
int isL, ieL, isU, ieU;
NL=6; NU= 6;
IAL = (int **)allocate_matrix(sizeof(int),ICELTOT,NL);IAU = (int **)allocate_matrix(sizeof(int),ICELTOT,NU);
BFORCE = (double *)allocate_vector(sizeof(double),ICELTOT);D= (double *)allocate_vector(sizeof(double),ICELTOT);
PHI= (double *)allocate_vector(sizeof(double),ICELTOT);INL = (int *)allocate_vector(sizeof(int),ICELTOT);INU = (int *)allocate_vector(sizeof(int),ICELTOT);
for (i = 0; i <ICELTOT ; i++) {BFORCE[i]=0.0;
D[i] =0.0; PHI[i]=0.0;
INL[i] = 0; INU[i] = 0;
for(j=0;j<6;j++){
IAL[i][j]=0; IAU[i][j]=0;}}
for (i = 0; i <=ICELTOT ; i++) {indexL[i] = 0; indexU[i] = 0;}28

poi_gen (1/9)/**************************
allocate matrix
allocate.c**************************/
void** allocate_matrix(int size,int m,int n){void **aa;
int i;
if ( ( aa=(void ** )malloc( m * sizeof(void*) ) ) == NULL ) {fprintf(stdout,"Error:Memory does not enough! aa in matrix ¥n");exit(1);}
if ( ( aa[0]=(void * )malloc( m * n * size ) ) == NULL ) {fprintf(stdout,"Error:Memory does not enough! in matrix ¥n");exit(1);}for(i=1;i<m;i++) aa[i]=(char*)aa[i-1]+size*n;
return aa;}

for(icel=0; icel<ICELTOT; icel++) {icN1 = NEIBcell[icel][0];icN2 = NEIBcell[icel][1];
icN3 = NEIBcell[icel][2];
icN4 = NEIBcell[icel][3];icN5 = NEIBcell[icel][4];
icN6 = NEIBcell[icel][5];29

poi_gen (2/9)
NEIBcell(icel,6)
NEIBcell[icel][5]if(icN5 != 0) {icou = INL[icel] + 1;
IAL[icel][icou-1] = icN5;
INL[icel]= icou;}if(icN3 != 0) {icou = INL[icel] + 1;
IAL[icel][icou-1] = icN3;
INL[icel]= icou;}if(icN1 != 0) {icou = INL[icel] + 1;
IAL[icel][icou-1] = icN3;
INL[icel]= icou;}if(icN2 != 0) {icou = INU[icel] + 1;
IAU[icel][icou-1] = icN2;
INU[icel]= icou;}if(icN4 != 0) {icou = INU[icel] + 1;
IAU[icel][icou-1] = icN4;
INU[icel]= icou;}}if(icN6 != 0) {icou = INU[icel] + 1;
IAU[icel][icou-1] = icN6;
INU[icel]= icou;}NEIBcell(icel,4)NEIBcell[icel][3]NEIBcell[icel][0]NEIBcell(icel,1)z

NEIBcell(icel,2)
NEIBcell[icel][1]NEIBcell(icel,3)
NEIBcell[icel][2]yNEIBcell[icel][4]NEIBcell(icel,5)x

Lower Triangular Part
NEIBcell[icel][4]= icel – NX*NY + 1
NEIBcell[icel][2]= icel – NX
+1NEIBcell[icel][0]= icel – 1+1“icel” starts at 0

“IAL” starts at 1

12 13 14 15

13 14 15 16

8

9

10 11

9

10 11 12

4

5

6

7

5

6

7

8

0

1

2

3

1

2

3

4

for(icel=0; icel<ICELTOT; icel++) {icN1 = NEIBcell[icel][0];icN2 = NEIBcell[icel][1];
icN3 = NEIBcell[icel][2];
icN4 = NEIBcell[icel][3];icN5 = NEIBcell[icel][4];
icN6 = NEIBcell[icel][5];if(icN5 != 0) {icou = INL[icel] + 1;
IAL[icel][icou-1] = icN5;
INL[icel]= icou;}if(icN3 != 0) {icou = INL[icel] + 1;
IAL[icel][icou-1] = icN3;
INL[icel]= icou;}if(icN1 != 0) {icou = INL[icel] + 1;
IAL[icel][icou-1] = icN3;
INL[icel]= icou;}if(icN2 != 0) {icou = INU[icel] + 1;
IAU[icel][icou-1] = icN2;
INU[icel]= icou;}if(icN4 != 0) {icou = INU[icel] + 1;
IAU[icel][icou-1] = icN4;
INU[icel]= icou;}}if(icN6 != 0) {icou = INU[icel] + 1;
IAU[icel][icou-1] = icN6;
INU[icel]= icou;}30

poi_gen (2/9)
NEIBcell(icel,6)
NEIBcell[icel][5]NEIBcell(icel,4)NEIBcell[icel][3]NEIBcell[icel][0]NEIBcell(icel,1)z

NEIBcell(icel,2)
NEIBcell[icel][1]NEIBcell(icel,3)
NEIBcell[icel][2]yNEIBcell[icel][4]NEIBcell(icel,5)xUpper Triangular Part
NEIBcell[icel][1]= icel + 1+1
NEIBcell[icel][3]= icel + NX
+1NEIBcell[icel][5]= icel + NX*NY + 1“icel” starts at 0

“IAU” starts at 1

12 13 14 15

13 14 15 16

8

9

10 11

9

10 11 12

4

5

6

7

5

6

7

8

0

1

2

3

1

2

3

4

N111:

poi_gen(3/9)

fprintf(stderr, "¥n¥nYou have%8d elements¥n", ICELTOT);
fprintf(stderr, "How many colors do you need ?¥n");fprintf(stderr, " #COLOR must be more than 2 and¥n");fprintf(stderr, " #COLOR must not be more than%8d¥n", ICELTOT);
fprintf(stderr, " if #COLOR= 0 then CM ordering¥n");fprintf(stderr, " if #COLOR=-1 then RCM ordering¥n");fprintf(stderr, " if #COLOR<-1 then CMRCM ordering¥n");fprintf(stderr, "=>¥n");fscanf(stdin, "%d", &NCOLORtot);if(NCOLORtot == 1 && NCOLORtot > ICELTOT) goto N111;
OLDtoNEW = (int *)calloc(ICELTOT, sizeof(int));if(OLDtoNEW == NULL) {fprintf(stderr, "Error: %s¥n", strerror(errno));return -1;}
NEWtoOLD = (int *)calloc(ICELTOT, sizeof(int));
if(NEWtoOLD == NULL) {fprintf(stderr, "Error: %s¥n", strerror(errno));return -1;}
COLORindex = (int *)calloc(ICELTOT+1, sizeof(int));if(COLORindex == NULL) {fprintf(stderr, "Error: %s¥n", strerror(errno));return -1;}
if(NCOLORtot > 0) {MC(ICELTOT, NL, NU, INL, IAL, INU, IAU,
&NCOLORtot, COLORindex, NEWtoOLD, OLDtoNEW);} else if(NCOLORtot == 0) {CM(ICELTOT, NL, NU, INL, IAL, INU, IAU,
&NCOLORtot, COLORindex, NEWtoOLD, OLDtoNEW);} else if(NCOLORtot ==-1) {RCM(ICELTOT, NL, NU, INL, IAL, INU, IAU,
&NCOLORtot, COLORindex, NEWtoOLD, OLDtoNEW);} else if(NCOLORtot <-1) {CMRCM(ICELTOT, NL, NU, INL, IAL, INU, IAU,
&NCOLORtot, COLORindex, NEWtoOLD, OLDtoNEW);}

Reordering
NCOLORtot > 1: Multicolor
NCOLORtot = 0: CM
NCOLORtot =-1: RCM
NCOLORtot <-1: CM-RCM}
31

fprintf(stderr, "¥n# TOTAL COLOR number%8d¥n", NCOLORtot);
return 0;

SMPindex = (int *) allocate_vector(sizeof(int),
NCOLORtot*PEsmpTOT+1);
memset(SMPindex, 0,
sizeof(int)*(NCOLORtot*PEsmpTOT+1));
for(ic=1; ic<=NCOLORtot; ic++) {nn1 = COLORindex[ic] - COLORindex[ic-1];
num = nn1 / PEsmpTOT;
nr = nn1 - PEsmpTOT * num;for(ip=1; ip<=PEsmpTOT; ip++) {if(ip <= nr) {SMPindex[(ic-1)*PEsmpTOT+ip] = num + 1;
} else {SMPindex[(ic-1)*PEsmpTOT+ip] = num;}}}for(ic=1; ic<=NCOLORtot; ic++) {for(ip=1; ip<=PEsmpTOT; ip++) {j1 = (ic-1) * PEsmpTOT + ip;
j0 = j1 - 1;
SMPindex[j1] += SMPindex[j0];}}32

poi_gen (4/9)
SMPindex:
for preconditioning
各色内の要素数：COLORindex[ic]-COLORindex[ic-1]
同じ色内の要素は依存性が無いため，並列に計算可能 ⇒ OpenMP適用
これを更に「PEsmpTOT」で割って「SMPindex」に割り当てる。
前処理で使用SMPindexG = (int *) allocate_vector(sizeof(int),for(ic=0; ic<NCOLORtot; ic++) {PEsmpTOT+1);#pragma omp parallel for ...
memset(SMPindexG, 0, sizeof(int)*(PEsmpTOT+1));
nn = ICELTOT / PEsmpTOT;
nr = ICELTOT - nn * PEsmpTOT;
for(ip=1; ip<=PEsmpTOT; ip++) {SMPindexG[ip] = nn;
if(ip <= nr) {SMPindexG[ip] += 1;}}
for(ip=1; ip<=PEsmpTOT; ip++) {}SMPindexG[ip] += SMPindexG[ip-1];}for(ip=0; ip<PEsmpTOT; ip++) {ip1 = ic * PEsmpTOT + ip;
for(i=SMPindex[ip1]; i<SMPindex[ip1+1]; i++) {
(...)
}
}

33

SMPindex：
前処理向け

for(ic=0; ic<NCOLORtot; ic++) {#pragma omp parallel for ...
for(ip=0; ip<PEsmpTOT; ip++) {ip1 = ic * PEsmpTOT + ip;
for(i=SMPindex[ip1]; i<SMPindex[ip1+1]; i++) {
(...)
}
}
}

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

SMPindex = (int *) allocate_vector(sizeof(int),
NCOLORtot*PEsmpTOT+1);
memset(SMPindex, 0,
sizeof(int)*(NCOLORtot*PEsmpTOT+1));34

poi_gen (4/9)

for(ic=1; ic<=NCOLORtot; ic++) {nn1 = COLORindex[ic] - COLORindex[ic-1];
num = nn1 / PEsmpTOT;
nr = nn1 - PEsmpTOT * num;for(ip=1; ip<=PEsmpTOT; ip++) {if(ip <= nr) {SMPindex[(ic-1)*PEsmpTOT+ip] = num + 1;
} else {SMPindex[(ic-1)*PEsmpTOT+ip] = num;}}}
#pragma omp parallel for ...
for(ic=1; ic<=NCOLORtot; ic++) {for(ip=1; ip<=PEsmpTOT; ip++) {j1 = (ic-1) * PEsmpTOT + ip;
j0 = j1 - 1;
SMPindex[j1] += SMPindex[j0];}}for(ip=0; ip<PEsmpTOT; ip++) {for(i=SMPindexG[ip]; i<SMPindexG[ip]; i++) {
(...)
}
}

SMPindexG = (int *) allocate_vector(sizeof(int),
PEsmpTOT+1);
memset(SMPindexG, 0, sizeof(int)*(PEsmpTOT+1));
nn = ICELTOT / PEsmpTOT;
nr = ICELTOT - nn * PEsmpTOT;
for(ip=1; ip<=PEsmpTOT; ip++) {SMPindexG[ip] = nn;
if(ip <= nr) {SMPindexG[ip] += 1;}}
for(ip=1; ip<=PEsmpTOT; ip++) {SMPindexG[ip] += SMPindexG[ip-1];}

全要素数を「PEsmpTOT」で割って「SMPindexG」に割り当てる。
内積，行列ベクトル積，DAXPYで使用
これを使用すれば，実は，「poi_gen(2/9」の部分も並列化可能「poi_gen(5/9」以降では実際に使用

35

SMPindexG#pragma omp parallel for ...
for(ip=0; ip<PEsmpTOT; ip++) {for(i=SMPindexG[ip]; i<SMPindexG[ip]; i++) {
(...)
}
}

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

36

poi_gen (5/9)

indexL =(int *)allocate_vector(sizeof(int),ICELTOT+1);
indexU =(int *)allocate_vector(sizeof(int),ICELTOT+1);

これ以降は新しい
番号付けを使用

for(i=0; i<ICELTOT; i++){indexL[i+1]=indexL[i]+INL[i];indexU[i+1]=indexU[i]+INU[i];}NPL = indexL[ICELTOT];NPU = indexU[ICELTOT];

Name

itemL = (int *)allocate_vector(sizeof(int),NPL);
itemU = (int *)allocate_vector(sizeof(int),NPU);AL =(double *)allocate_vector(sizeof(double),NPL);AU =(double *)allocate_vector(sizeof(double),NPU);
memset(itemL, 0, sizeof(int)*NPL);
memset(itemU, 0, sizeof(int)*NPU);
memset(AL, 0.0, sizeof(double)*NPL);
memset(AU, 0.0, sizeof(double)*NPU);
for(i=0; i<ICELTOT; i++){for(k=0;k<INL[i];k++){kk= k + indexL[i];itemL[kk]= IAL[i][k];}for(k=0;k<INU[i];k++){kk= k + indexU[i];itemU[kk]= IAU[i][k];}}free(INL); free(INU);
free(IAL); free(IAU);

Type

Content

D[N]R

Diagonal components of the matrix (N= ICELTOT)BFORCE[N]R

RHS vector

PHI[N]

R

Unknown vector

indexL[N+1]indexU[N+1]I# of L/U non-zero off-diag. comp. (CRS)

NPL, NPUI

Total # of L/U non-zero off-diag. comp. (CRS)

itemL[NPL]
itemU[NPU]

I

Column ID of L/U non-zero off-diag. comp. (CRS)

AL[NPL]
AU[NPU]R

L/U non-zero off-diag. comp. (CRS)“itemL” /”itemU”
start at 1
13 14 15 169

10 11 12

5

6

7

8

1

2

3

4

for (i=0; i<N; i++) {q[i]= D[i] * p[i];for (j=indexL[i]; j<indexL[i+1]; j++) {q[i] += AL[j] * p[itemL[j]-1];}
for (j=indexU[i]; j<indexU[i+1]; j++) {q[i] += AU[j] * p[itemU[j]-1];
}
}

37

S1t = omp_get_wtime();#pragma omp parallel for private(ip, icel, ic0, icN1, icN2, icN3, icN4, icN5, icN6, coef, j,
ii, jj, kk, isL, ieL, isU, ieU)
for(ip=0; ip<PEsmpTOT; ip++) {for(icel=SMPindexG[ip]; icel<SMPindexG[ip+1]; icel++) {ic0 = NEWtoOLD[icel];
icN1 = NEIBcell[ic0-1][0];icN2 = NEIBcell[ic0-1][1];
icN3 = NEIBcell[ic0-1][2];
icN4 = NEIBcell[ic0-1][3];icN5 = NEIBcell[ic0-1][4];
icN6 = NEIBcell[ic0-1][5];VOL0 = VOLCEL[ic0];
isL =indexL[icelisU =indexU[icel];];ieL =indexL[icel+1];
ieU =indexU[icel+1];

if(icN5 != 0) {icN5 = OLDtoNEW[icN5-1];
coef = RDZ * ZAREA;
D[icel] -= coef;
if(icN5-1 < icel) {for(j=isL; j<ieL; j++) {if(itemL[j] == icN5) {AL[j] = coef;
break;}}} else {for(j=isU; j<ieU; j++) {if(itemU[j] == icN5) {AU[j] = coef;
break;
}
}
}

...

}

icel: New ID
ic0: Old ID

poi_gen (6/9)
新しい番号付けを使用
neib[icel ][0 ]  icel
xneib[icel ][1]  icel
xneib[icel ][ 2]  icelyneib[icel ][3]  icelyneib[icel ][ 4]  icel
zneib[icel ][5]  icelzyz 
yz 
zx 
zx 
xy 
xy  f icel xyz  038係数の計算：並列に実施可能
SMPindexG を使用
private宣言に注意#pragma omp parallel for private(ip, icel, ic0, icN1, icN2, icN3, icN4, icN5, icN6, coef, j, ii,
jj, kk, isL, ieL, isU, ieU)
for(ip=0; ip<PEsmpTOT; ip++) {for(icel=SMPindexG[ip]; icel<SMPindexG[ip+1]; icel++) {ic0 = NEWtoOLD[icel];
icN1 = NEIBcell[ic0-1][0];39S1t = omp_get_wtime();#pragma omp parallel for private(ip, icel, ic0, icN1, icN2, icN3, icN4, icN5, icN6, coef, j,
ii, jj, kk, isL, ieL, isU, ieU)
for(ip=0; ip<PEsmpTOT; ip++) {for(icel=SMPindexG[ip]; icel<SMPindexG[ip+1]; icel++) {ic0 = NEWtoOLD[icel];
icN1 = NEIBcell[ic0-1][0];icN2 = NEIBcell[ic0-1][1];
icN3 = NEIBcell[ic0-1][2];
icN4 = NEIBcell[ic0-1][3];icN5 = NEIBcell[ic0-1][4];
icN6 = NEIBcell[ic0-1][5];VOL0 = VOLCEL[ic0];
isL =indexL[icelisU =indexU[icel];];ieL =indexL[icel+1];
ieU =indexU[icel+1];

if(icN5 != 0) {icN5 = OLDtoNEW[icN5-1];
coef = RDZ * ZAREA;
D[icel] -= coef;
if(icN5-1 < icel) {for(j=isL; j<ieL; j++) {if(itemL[j] == icN5) {AL[j] = coef;
break;}}} else {for(j=isU; j<ieU; j++) {if(itemU[j] == icN5) {AU[j] = coef;
break;
}
}
}

...

}

icel: New ID
ic0: Old ID

poi_gen (6/9)
New numbering appliedneib[icel ][0 ]  icel
xneib[icel ][1]  icel
xneib[icel ][ 2]  icelyneib[icel ][3]  icelyneib[icel ][ 4]  icel
zneib[icel ][5]  icelzyz 
yz 
zx 
zx 
xy 
xy  f icel xyz  0

40

S1t = omp_get_wtime();#pragma omp parallel for private(ip, icel, ic0, icN1, icN2, icN3, icN4, icN5, icN6, coef, j,
ii, jj, kk, isL, ieL, isU, ieU)
for(ip=0; ip<PEsmpTOT; ip++) {for(icel=SMPindexG[ip]; icel<SMPindexG[ip+1]; icel++) {ic0 = NEWtoOLD[icel];
icN1 = NEIBcell[ic0-1][0];icN2 = NEIBcell[ic0-1][1];
icN3 = NEIBcell[ic0-1][2];
icN4 = NEIBcell[ic0-1][3];icN5 = NEIBcell[ic0-1][4];
icN6 = NEIBcell[ic0-1][5];VOL0 = VOLCEL[ic0];
isL =indexL[icelisU =indexU[icel];];ieL =indexL[icel+1];
ieU =indexU[icel+1];

if(icN5 != 0) {icN5 = OLDtoNEW[icN5-1];
coef = RDZ * ZAREA;
D[icel] -= coef;
if(icN5-1 < icel) {for(j=isL; j<ieL; j++) {if(itemL[j] == icN5) {AL[j] = coef;
break;}}} else {for(j=isU; j<ieU; j++) {if(itemU[j] == icN5) {AU[j] = coef;
break;
}
}
}

...

}

poi_gen (6/9)
New numbering appliedneib[icel ][0 ]  icel
xneib[icel ][1]  icel
xneib[icel ][ 2]  icelyneib[icel ][3]  icelyneib[icel ][ 4]  icel
zneib[icel ][5]  icelzyz 
yz 
zx 
zx 
xy 
xy  f icel xyz  0

41

S1t = omp_get_wtime();#pragma omp parallel for private(ip, icel, ic0, icN1, icN2, icN3, icN4, icN5, icN6, coef, j,
ii, jj, kk, isL, ieL, isU, ieU)
for(ip=0; ip<PEsmpTOT; ip++) {for(icel=SMPindexG[ip]; icel<SMPindexG[ip+1]; icel++) {ic0 = NEWtoOLD[icel];
icN1 = NEIBcell[ic0-1][0];icN2 = NEIBcell[ic0-1][1];
icN3 = NEIBcell[ic0-1][2];
icN4 = NEIBcell[ic0-1][3];icN5 = NEIBcell[ic0-1][4];
icN6 = NEIBcell[ic0-1][5];VOL0 = VOLCEL[ic0];
isL =indexL[icelisU =indexU[icel];];if(icN5-1 < icel) {for(j=isL; j<ieL; j++) {if(itemL[j] == icN5) {AL[j] = coef;
break;}}} else {for(j=isU; j<ieU; j++) {if(itemU[j] == icN5) {AU[j] = coef;
break;
}
}
}

...

}

New numbering appliedneib[icel ][0 ]  icel
x

ieL =indexL[icel+1];
ieU =indexU[icel+1];

if(icN5 != 0) {icN5 = OLDtoNEW[icN5-1];
coef = RDZ * ZAREA;
D[icel] -= coef;

poi_gen (6/9)

neib[icel ][1]  icel
x1
z
ZAREA  xyneib[icel ][ 2]  icel

icN5 < icel
Lower Part

y

RDZ 

y

neib[icel ][3]  icel
neib[icel ][ 4]  icel
zneib[icel ][5]  icelzyz 
yz 
zx 
zx 
xy 
xy  f icel xyz  0

42

S1t = omp_get_wtime();#pragma omp parallel for private(ip, icel, ic0, icN1, icN2, icN3, icN4, icN5, icN6, coef, j,
ii, jj, kk, isL, ieL, isU, ieU)
for(ip=0; ip<PEsmpTOT; ip++) {for(icel=SMPindexG[ip]; icel<SMPindexG[ip+1]; icel++) {ic0 = NEWtoOLD[icel];
icN1 = NEIBcell[ic0-1][0];icN2 = NEIBcell[ic0-1][1];
icN3 = NEIBcell[ic0-1][2];
icN4 = NEIBcell[ic0-1][3];icN5 = NEIBcell[ic0-1][4];
icN6 = NEIBcell[ic0-1][5];VOL0 = VOLCEL[ic0];
isL =indexL[icelisU =indexU[icel];];if(icN5-1 < icel) {for(j=isL; j<ieL; j++) {if(itemL[j] == icN5) {AL[j] = coef;
break;}}} else {for(j=isU; j<ieU; j++) {if(itemU[j] == icN5) {AU[j] = coef;
break;
}
}
}

...

}

New numbering appliedneib[icel ][0 ]  icel
x

ieL =indexL[icel+1];
ieU =indexU[icel+1];

if(icN5 != 0) {icN5 = OLDtoNEW[icN5-1];
coef = RDZ * ZAREA;
D[icel] -= coef;

poi_gen (6/9)

1z
ZAREA  xyRDZ neib[icel ][1]  icel
xneib[icel ][ 2]  icelyneib[icel ][3]  icelyneib[icel ][ 4]  icel
icN5 > icel
Upper Partzneib[icel ][5]  icelzyz 
yz 
zx 
zx 
xy 
xy  f icel xyz  0

if(icN3 != 0) {icN3 = OLDtoNEW[icN3-1];
coef = RDY * YAREA;
D[icel] -= coef;}

if(icN3-1 < icel) {for(j=isL; j<ieL; j++) {if(itemL[j] == icN3) {AL[j] = coef;
break; }}} else {for(j=isU; j<ieU; j++) {if(itemU[j] == icN3) {AU[j] = coef;
break; }}}if(icN1 != 0) {icN1 = OLDtoNEW[icN1-1];
coef = RDX * XAREA;
D[icel] -= coef;}

if(icN1-1 < icel) {for(j=isL; j<ieL; j++) {if(itemL[j] == icN1) {AL[j] = coef;
break;}}} else {for(j=isU; j<ieU; j++) {if(itemU[j] == icN1) {AU[j] = coef;
break;}}}

43

poi_gen (7/9)
neib[icel ][0 ]  icel
xneib[icel ][1]  icel
xneib[icel ][ 2]  icelyneib[icel ][3]  icelyneib[icel ][ 4]  icel
zneib[icel ][5]  icelzyz 
yz 
zx 
zx 
xy 
xy  f icel xyz  0

if(icN2 != 0) {icN2 = OLDtoNEW[icN2-1];
coef = RDX * XAREA;
D[icel] -= coef;}

if(icN2-1 < icel) {for(j=isL; j<ieL; j++) {if(itemL[j] == icN2) {AL[j] = coef;
break;}}} else {for(j=isU; j<ieU; j++) {if(itemU[j] == icN2) {AU[j] = coef;
break;}}}

if(icN4 != 0) {icN4 = OLDtoNEW[icN4-1];
coef = RDY * YAREA;
D[icel] -= coef;}

if(icN4-1 < icel) {for(j=isL; j<ieL; j++) {if(itemL[j] == icN4) {AL[j] = coef;
break; }}} else {for(j=isU; j<ieU; j++) {if(itemU[j] == icN4) {AU[j] = coef;
break; }}}44

poi_gen (8/9)
neib[icel ][0 ]  icel
xneib[icel ][1]  icel
xneib[icel ][ 2]  icelyneib[icel ][3]  icelyneib[icel ][ 4]  icel
zneib[icel ][5]  icelzyz 
yz 
zx 
zx 
xy 
xy  f icel xyz  0#pragma omp parallel for private(ip, icel, ic0, icN1, icN2, icN3, icN4, icN5, icN6,
coef, j, ii, jj, kk, isL, ieL, isU, ieU)…

45

poi_gen (9/9)

if(icN6 != 0) {icN6 = OLDtoNEW[icN5-1];
coef = RDZ * ZAREA;
D[icel] -= coef;
if(icN6-1 < icel) {for(j=isL; j<ieL; j++) {if(itemL[j] == icN6) {AL[j] = coef;
break;}}} else {for(j=isU; j<ieU; j++) {if(itemU[j] == icN6) {AU[j] = coef;
break;}}}}
ii = XYZ[ic0][0];jj = XYZ[ic0][1];kk = XYZ[ic0][2];BFORCE
using original
mesh ID

BFORCE[icel]= -(double)(ii+jj+kk) * VOL0;ii,jj,kk,VOL0:
private}neib[icel ][0 ]  icel
xneib[icel ][1]  icel
xneib[icel ][ 2]  icelyneib[icel ][3]  icelyneib[icel ][ 4]  icel
zneib[icel ][5]  icelzyz 
yz 
zx 
zx 
xy 
xy  f icel xyz  0Main Program#include <stdio.h> ...
intmain(){

double *WK;
int NPL, NPU; ISET, ITR, IER; icel, ic0, i;
double xN, xL, xU; Stime, Etime;if(INPUT()) goto error;
if(POINTER_INIT()) goto error;
if(BOUNDARY_CELL()) goto error;if(CELL_METRICS()) goto error;
if(POI_GEN()) goto error;

error:}ISET = 0;
WK = (double *)malloc(sizeof(double)*ICELTOT);if(WK == NULL) {fprintf(stderr, "Error: %s¥n", strerror(errno));goto error;}
Stime = omp_get_wtime();if(solve_ICCG_mc(ICELTOT, NL, NU, indexL, itemL, indexU, itemU,
D, BFORCE, PHI, AL, AU, NCOLORtot, PEsmpTOT,
SMPindex, SMPindexG, EPSICCG, &ITR, &IER)) goto error;
Etime = omp_get_wtime();for(ic0=0; ic0<ICELTOT; ic0++) {icel = NEWtoOLD[ic0];
WK[icel-1] = PHI[ic0];}
for(icel=0; icel<ICELTOT; icel++) {PHI[icel] = WK[icel];}if(OUTUCD()) goto error;
return 0;
return -1;

46

47

solve_ICCG_mc (1/6)#include <stdio.h>#include <stdlib.h>#include <string.h>#include <errno.h>#include <math.h> etc.
#include "solver_ICCG.h"
extern int
solve_ICCG_mc(int N, int NL, int NU, int *indexL, int *itemL, int *indexU,
int *itemU,
double *D, double *B, double *X, double *AL, double *AU,
int NCOLORtot, int *COLORindex,
int PEsmpTOT, int *SMPindex, int *SMPindexG,
double EPS, int *ITR, int *IER){
double **W;double VAL, BNRM2, WVAL, SW, RHO, BETA, RHO1, C1, DNRM2, ALPHA, ERR;
int i, j, ic, ip, L, ip1;
int R = 0;
int Z = 1;
int Q = 1;
int P = 2;
int DD = 3;

48

W = (double **)malloc(sizeof(double *)*4);
if(W == NULL) {fprintf(stderr, "Error: %s¥n", strerror(errno));return -1;}
for(i=0; i<4; i++) {W[i] = (double *)malloc(sizeof(double)*N);if(W[i] == NULL) {fprintf(stderr, "Error: %s¥n",strerror(errno)); return -1;}}

solve_ICCG_mc(2/6)#pragma omp parallel for private (ip, i)for(ip=0; ip<PEsmpTOT; ip++) {for(i=SMPindexG[ip]; i<SMPindexG[ip+1]; i++) {X[i] = 0.0;
W[1][i] = 0.0;
W[2][i] = 0.0;
W[3][i] = 0.0;}}
for(ic=0; ic<NCOLORtot; ic++) {#pragma omp parallel for private (ip, ip1, i, VAL, j)for(ip=0; ip<PEsmpTOT; ip++) {ip1 = ic * PEsmpTOT + ip;
for(i=SMPindex[ip1]; i<SMPindex[ip1+1]; i++) {VAL = D[i];for(j=indexL[i]; j<indexL[i+1]; j++) {VAL-= AL[j]*AL[j]*W[DD][itemL[j] - 1];}
W[DD][i] = 1.0 / VAL;}}}

Incomplete “Modified”
Cholesky
Factorization

49

Incomplete “Modified” Cholesky
Factorization
1

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



W[DD][i]:
D[i]:

di
aii

itemL[j]:AL[j]:kaik

for(i=0; i<N; i++) {VAL = D[i];for(j=indexL[i]; j<indexL[i+1]; j++) {VAL-= AL[j]*AL[j]*W[DD][itemL[j] - 1];}
W[DD][i] = 1.0 / VAL;}

50

Incomplete “Modified” Cholesky
Factorization: Parallel Version
1

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



W[DD][i]:
D[i]:

di
aii

itemL[j]:AL[j]:kaik

for(ic=0; ic<NCOLORtot; ic++) {#pragma omp parallel for private (ip,ip1,i,VAL,j)for(ip=0; ip<PEsmpTOT; ip++) {ip1 = ic * PEsmpTOT + ip;
for(i=SMPindex[ip1]; i<SMPindex[ip1+1]; i++) {VAL = D[i];for(j=indexL[i]; j<indexL[i+1]; j++) {VAL-= AL[j]*AL[j]*W[DD][itemL[j] - 1];}
W[DD][i] = 1.0 / VAL;}}}

privateに注意。

51solve_ICCG_mc (3/6)#pragma omp parallel for private (ip,i,VAL,j)for(ip=0; ip<PEsmpTOT; ip++) {for(i=SMPindexG[ip]; i<SMPindexG[ip+1]; i++) {VAL = D[i] * X[i];for(j=indexL[i]; j<indexL[i+1]; j++) {VAL += AL[j] * X[itemL[j]-1];}
for(j=indexU[i]; j<indexU[i+1]; j++) {VAL += AU[j] * X[itemU[j]-1];}}
W[R][i] = B[i] - VAL;}BNRM2 = 0.0;
#pragma omp parallel for private (ip, i)
reduction (+:BNRM2)
for(ip=0; ip<PEsmpTOT; ip++) {for(i=SMPindexG[ip]; i<SMPindexG[ip+1]; i++) {BNRM2 += B[i]*B[i];}}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

52

Mat-Vec
NO Data Dependency: SMPindexG#pragma omp parallel for private (ip,i,VAL,j)for(ip=0; ip<PEsmpTOT; ip++) {for(i=SMPindexG[ip]; i<SMPindexG[ip+1]; i++) {VAL = D[i] * X[i];for(j=indexL[i]; j<indexL[i+1]; j++) {VAL += AL[j] * X[itemL[j]-1];}
for(j=indexU[i]; j<indexU[i+1]; j++) {VAL += AU[j] * X[itemU[j]-1];}}
W[R][i] = B[i] - VAL;}

53

solve_ICCG_mc (3/6)#pragma omp parallel for private (ip,i,VAL,j)for(ip=0; ip<PEsmpTOT; ip++) {for(i=SMPindexG[ip]; i<SMPindexG[ip+1]; i++) {VAL = D[i] * X[i];for(j=indexL[i]; j<indexL[i+1]; j++) {VAL += AL[j] * X[itemL[j]-1];}
for(j=indexU[i]; j<indexU[i+1]; j++) {VAL += AU[j] * X[itemU[j]-1];}}
W[R][i] = B[i] - VAL;}BNRM2 = 0.0;
#pragma omp parallel for private (ip, i)
reduction (+:BNRM2)
for(ip=0; ip<PEsmpTOT; ip++) {for(i=SMPindexG[ip]; i<SMPindexG[ip+1]; i++) {BNRM2 += B[i]*B[i];}}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

54

Dot Products: SMPindexG, reduction
BNRM2 = 0.0;
#pragma omp parallel for private (ip,i)
reduction (+:BNRM2)
for(ip=0; ip<PEsmpTOT; ip++) {for(i=SMPindexG[ip]; i<SMPindexG[ip+1]; i++) {BNRM2 += B[i]*B[i];}}*ITR = N;
for(L=0; L<(*ITR); L++) {#pragma omp parallel for private(ip, i)for(ip=0; ip<PEsmpTOT; ip++) {for(i=SMPindexG[ip]; i<SMPindexG[ip+1]; i++) {W[Z][i] = W[R][i];}}for(ic=0; ic<NCOLORtot; ic++) {#pragma omp parallel for private (ip,ip1,i,VAL,j)for(ip=0; ip<PEsmpTOT; ip++) {ip1 = ic * PEsmpTOT + ip;
for(i=SMPindex[ip1]; i<SMPindex[ip1+1]; i++) {WVAL = W[Z][i];for(j=indexL[i]; j<indexL[i+1]; j++) {WVAL -= AL[j] * W[Z][itemL[j]-1];}
W[Z][i] = WVAL * W[DD][i];}}}for(ic=NCOLORtot-1; ic>=0; ic--){#pragma omp parallel for private (ip,ip1,i,VAL,j)for(ip=0; ip<PEsmpTOT; ip++) {ip1 = ic * PEsmpTOT + ip;
for(i=SMPindex[ip1]; i<SMPindex[ip1+1]; i++) {SW = 0.0;
for(j=indexU[i]; j<indexU[i+1]; j++) {SW += AU[j] * W[Z][itemU[j]-1];}
W[Z][i]= W[Z][i] - W[DD][i] * SW;
}
}
}

55

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

56

*ITR = N;
for(L=0; L<(*ITR); L++) {#pragma omp parallel for private(ip, i)for(ip=0; ip<PEsmpTOT; ip++) {for(i=SMPindexG[ip]; i<SMPindexG[ip+1]; i++) {W[Z][i] = W[R][i];}}SMPindexfor(ic=0; ic<NCOLORtot; ic++) {#pragma omp parallel for private (ip,ip1,i,VAL,j)for(ip=0; ip<PEsmpTOT; ip++) {ip1 = ic * PEsmpTOT + ip;
for(i=SMPindex[ip1]; i<SMPindex[ip1+1]; i++) {WVAL = W[Z][i];for(j=indexL[i]; j<indexL[i+1]; j++) {WVAL -= AL[j] * W[Z][itemL[j]-1];}
W[Z][i] = WVAL * W[DD][i];}}}for(ic=NCOLORtot-1; ic>=0; ic--){#pragma omp parallel for private (ip,ip1,i,VAL,j)for(ip=0; ip<PEsmpTOT; ip++) {ip1 = ic * PEsmpTOT + ip;
for(i=SMPindex[ip1]; i<SMPindex[ip1+1]; i++) {SW = 0.0;
for(j=indexU[i]; j<indexU[i+1]; j++) {SW += AU[j] * W[Z][itemU[j]-1];}
W[Z][i]= W[Z][i] - W[DD][i] * SW;}}}

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

57

*ITR = N;
for(L=0; L<(*ITR); L++) {#pragma omp parallel for private(ip, i)for(ip=0; ip<PEsmpTOT; ip++) {for(i=SMPindexG[ip]; i<SMPindexG[ip+1]; i++) {W[Z][i] = W[R][i];}}SMPindexfor(ic=0; ic<NCOLORtot; ic++) {#pragma omp parallel for private (ip,ip1,i,VAL,j)for(ip=0; ip<PEsmpTOT; ip++) {ip1 = ic * PEsmpTOT + ip;
for(i=SMPindex[ip1]; i<SMPindex[ip1+1]; i++) {WVAL = W[Z][i];for(j=indexL[i]; j<indexL[i+1]; j++) {WVAL -= AL[j] * W[Z][itemL[j]-1];}
W[Z][i] = WVAL * W[DD][i];}}}for(ic=NCOLORtot-1; ic>=0; ic--){#pragma omp parallel for private (ip,ip1,i,VAL,j)for(ip=0; ip<PEsmpTOT; ip++) {ip1 = ic * PEsmpTOT + ip;
for(i=SMPindex[ip1]; i<SMPindex[ip1+1]; i++) {SW = 0.0;
for(j=indexU[i]; j<indexU[i+1]; j++) {SW += AU[j] * W[Z][itemU[j]-1];}
W[Z][i]= W[Z][i] - W[DD][i] * SW;}}}

solve_ICCG_mc(4/6)
M z  LDL z  r
TL z  r
Forward SubstitutionDL z  z
T

Backward Substitution

58

Forward Substitution: SMPindex
for(ic=0; ic<NCOLORtot; ic++) {#pragma omp parallel for private (ip,ip1,i,VAL,j)for(ip=0; ip<PEsmpTOT; ip++) {ip1 = ic * PEsmpTOT + ip;
for(i=SMPindex[ip1]; i<SMPindex[ip1+1]; i++) {WVAL = W[Z][i];for(j=indexL[i]; j<indexL[i+1]; j++) {WVAL -= AL[j] * W[Z][itemL[j]-1];}
W[Z][i] = WVAL * W[DD][i];}}}/********************************* {p} = {z} if
ITER=0** BETA = RHO / RHO1 otherwise *********************************/
if(L == 0) {#pragma omp parallel for private(ip,i)for(ip=0; ip<PEsmpTOT; ip++) {for(i=SMPindexG[ip]; i<SMPindexG[ip+1]; i++) {W[P][i] = W[Z][i];}}} else {BETA = RHO / RHO1;#pragma omp parallel for private(ip,i)for(ip=0; ip<PEsmpTOT; ip++) {for(i=SMPindexG[ip];]; i++) {W[P][i] = W[Z][i] + BETA * W[P][i];}}}/***************** {q} = [A]{p} *****************/#pragma omp parallel for private(ip, i)for(ip=0; ip<PEsmpTOT; ip++) {for(i=SMPindexG[ip]; i<SMPindexG[ip+1]; i++) {VAL = D[i] * W[P][i];for(j=indexL[i]; j<indexL[i+1]; j++) {VAL += AL[j] * W[P][itemL[j]-1];}
for(j=indexU[i]; j<indexU[i+1]; j++) {VAL += AU[j] * W[P][itemU[j]-1];}
W[Q][i] = VAL;}}

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

/********************************* {p} = {z} if
ITER=0** BETA = RHO / RHO1 otherwise *********************************/
if(L == 0) {#pragma omp parallel for private(ip,i)for(ip=0; ip<PEsmpTOT; ip++) {for(i=SMPindexG[ip]; i<SMPindexG[ip+1]; i++) {W[P][i] = W[Z][i];}}} else {BETA = RHO / RHO1;#pragma omp parallel for private(ip,i)for(ip=0; ip<PEsmpTOT; ip++) {for(i=SMPindexG[ip];]; i++) {W[P][i] = W[Z][i] + BETA * W[P][i];}}}/***************** {q} = [A]{p} *****************/#pragma omp parallel for private(ip,i)for(ip=0; ip<PEsmpTOT; ip++) {for(i=SMPindexG[ip]; i<SMPindexG[ip+1]; i++) {VAL = D[i] * W[P][i];for(j=indexL[i]; j<indexL[i+1]; j++) {VAL += AL[j] * W[P][itemL[j]-1];}
for(j=indexU[i]; j<indexU[i+1]; j++) {VAL += AU[j] * W[P][itemU[j]-1];}
W[Q][i] = VAL;}}

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

61

/************************* ALPHA = RHO / {p}{q} *************************/
C1 = 0.0;#pragma omp parallel for private(ip,i)reduction(+:C1)
for(ip=0; ip<PEsmpTOT; ip++) {for(i=SMPindexG[ip]; i<SMPindexG[ip+1]; i++){C1 += W[P][i] * W[Q][i];}}ALPHA = RHO / C1;/**************************** {x} = {x} + ALPHA * {p} ** {r} = {r} - ALPHA * {q} ****************************/#pragma omp parallel for private(ip,i)for(ip=0; ip<PEsmpTOT; ip++) {for(i=SMPindexG[ip]; i<SMPindexG[ip+1]; i++){X[i]+= ALPHA * W[P][i];W[R][i] -= ALPHA * W[Q][i];}}DNRM2 = 0.0;
#pragma omp parallel for private(ip,i)
reduction(+:DNRM2)
for(ip=0; ip<PEsmpTOT; ip++) {for(i=SMPindexG[ip]; i<SMPindexG[ip+1]; i++){DNRM2 += W[R][i]*W[R][i];}}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

62

/************************* ALPHA = RHO / {p}{q} *************************/
C1 = 0.0;#pragma omp parallel for private(ip,i)reduction(+:C1)
for(ip=0; ip<PEsmpTOT; ip++) {for(i=SMPindexG[ip]; i<SMPindexG[ip+1]; i++){C1 += W[P][i] * W[Q][i];}}ALPHA = RHO / C1;/**************************** {x} = {x} + ALPHA * {p} ** {r} = {r} - ALPHA * {q} ****************************/#pragma omp parallel for private(ip,i)for(ip=0; ip<PEsmpTOT; ip++) {for(i=SMPindexG[ip]; i<SMPindexG[ip+1]; i++){X[i]+= ALPHA * W[P][i];W[R][i] -= ALPHA * W[Q][i];}}DNRM2 = 0.0;
#pragma omp parallel for private(ip,i)
reduction(+:DNRM2)
for(ip=0; ip<PEsmpTOT; ip++) {for(i=SMPindexG[ip]; i<SMPindexG[ip+1]; i++){DNRM2 += W[R][i]*W[R][i];}}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

63

/************************* ALPHA = RHO / {p}{q} *************************/
C1 = 0.0;#pragma omp parallel for private(ip,i)reduction(+:C1)
for(ip=0; ip<PEsmpTOT; ip++) {for(i=SMPindexG[ip]; i<SMPindexG[ip+1]; i++){C1 += W[P][i] * W[Q][i];}}ALPHA = RHO / C1;/**************************** {x} = {x} + ALPHA * {p} ** {r} = {r} - ALPHA * {q} ****************************/#pragma omp parallel for private(ip,i)for(ip=0; ip<PEsmpTOT; ip++) {for(i=SMPindexG[ip]; i<SMPindexG[ip+1]; i++){X[i]+= ALPHA * W[P][i];W[R][i] -= ALPHA * W[Q][i];}}DNRM2 = 0.0;
#pragma omp parallel for private(ip,i)
reduction(+:DNRM2)
for(ip=0; ip<PEsmpTOT; ip++) {for(i=SMPindexG[ip]; i<SMPindexG[ip+1]; i++){DNRM2 += W[R][i]*W[R][i];}}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

64

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

OMP-365

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

L1
L3

Memory

L3

L1

Core Core
L1

Memory

L3

L1

Core Core
L1

Memory

L3

L1

Core Core
L1

L1

L2

L2

L2

L2

L2

L2

L2

L2

L1

Core CPU

Memory

L1
L3

L1

Core Core

Memory

L1
L3

L1

Core Core

Memory

L1
L3

L1

Core Core

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

Incompatible Point #ITERATIONS

350

300

250

1.E+04

1.E+02

反復回数

ICP#

200
1.E+00

1.E+01

1.E+02

1.E+03

1.E+00
1.E+00

1.E+04

1.E+01COLOR#

1.E+02

1.E+03

1.E+04COLOR#

1.0E-02

3.00

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

2.50
sec.

OMP-3

2.00

8.0E-03

6.0E-03

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

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

1.E+03

4.0E-03
1.E+00

1.E+01

1.E+02COLOR#

1.E+0366

OMP-3

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

550

反復回数

500

Iterations

450
400
350
300

MC

RCM

CM-RCM

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

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

CM-RCM

1.E+00

250
1

10

100

1

1000

COLOR#

25.00

1000

MC
5.00E-02

sec./iteration

20.00

100

COLOR#

6.00E-02

MC
RCM
CM-RCM

10

sec.

15.00

10.00

5.00

RCM
CM-RCM

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

計算時間
0.00

計算時間／反復

0.00E+00
1

10

COLOR#

100

1000

1

10

COLOR#

100

1000

67

68

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

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

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

70

71

Running L3-sol

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

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

72

制御データ（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法

各方向の要素数73

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

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

74

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

20.0
MC=2

MC=2

RCM(298)

16.0

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

Speed-Up

20.0

10.0

12.0

8.0

4.0
0.0

0.0
1

2

4

8

16

1

2

4

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

8

16

thread#Memory

L1

L1

L1

L1

L1

L1

L1

L2
L1 L1

L1

L1

L1

L1

L1

L1

L1

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

75

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

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

77

前進代入：現状の並列化（C）
for(ic=0; ic<NCOLORtot; ic++) {#pragma omp parallel for private (ip, ip1, i, WVAL, j)for(ip=0; ip<PEsmpTOT; ip++) {ip1 = ic * PEsmpTOT + ip;
for(i=SMPindex[ip1]; i<SMPindex[ip1+1]; i++){WVAL = W[Z][i];for(j=indexL[i]; j<indexL[i+1]; j++){WVAL -= AL[j] * W[Z][itemL[j]-1];}
W[Z][i] = WVAL * W[DD][i];}}}• 「#pragma omp parallel」でスレッド（～16）の生成，消滅が発生– 色ごとにこの部分を通る– 多少のオーバーヘッドがある• 色数が増えるとオーバーヘッドが増す

78

前進代入：Overhead削減（C）#pragma omp parallel private (ic, ip, ip1, i, WVAL, j)for(ic=0; ic<NCOLORtot; ic++) {#pragma omp for
for(ip=0; ip<PEsmpTOT; ip++) {ip1 = ic * PEsmpTOT + ip;
for(i=SMPindex[ip1]; i<SMPindex[ip1+1]; i++){WVAL = W[Z][i];for(j=indexL[i]; j<indexL[i+1]; j++){WVAL -= AL[j] * W[Z][itemL[j]-1];}
W[Z][i] = WVAL * W[DD][i];}}}• このようにすることによって，スレッド生成を前進代入に入る前の一回で済ませることができる• 「#pragma omp for」のループが並列化

79

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

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

src0

80

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

L3-sol

L3-sol05.69 sec.

5.44 sec.

6.54 sec.

6.37 sec.

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

82

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

SMPindex：
前処理向け

for(ic=0; ic<NCOLORtot; ic++) {#pragma omp parallel for ...
for(ip=0; ip<PEsmpTOT; ip++) {ip1 = ic * PEsmpTOT + ip;
for(i=SMPindex[ip1]; i<SMPindex[ip1+1]; i++) {
(...)
}
}
}

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

84

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

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

86

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

87

データ再配置：Sequential Reordering
CM-RCM(2), 4-threads
1st-Color■#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

88

データ再配置：Sequential Reordering
CM-RCM(2), 4-threads
2nd-Color■#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

89

データ再配置: Sequential Reordering
Coalesced
Good for GPU

Initial Vector

Coloring(5 colors)+Orderingcolor=1

color=1

color=2

color=2

color=3

color=4

color=3

color=4

color=5

color=512345678 12345678 12345678 12345678 12345678

Sequential

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

color=1

color=1

color=2

color=2

color=3

color=4

color=5

color=3

color=4

color=512345678 12345678 12345678 12345678 12345678

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

OMP-390

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

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

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

変数名

型

内

容

PEsmpTOT

整数

データ分割数NCOLORtot

整数

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

NFLAG

整数

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

整数

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

92

Sequential Reordering
SMPindex = (int *) allocate_vector(sizeof(int), NCOLORtot * PEsmpTOT + 1);
memset(SMPindex, 0, sizeof(int)*(NCOLORtot*PEsmpTOT+1));
for(ic=1; ic<=NCOLORtot; ic++) {nn1 = COLORindex[ic] - COLORindex[ic-1];
num = nn1 / PEsmpTOT;
nr = nn1 - PEsmpTOT * num;for(ip=1; ip<=PEsmpTOT; ip++) {if(ip <= nr) {SMPindex[(ic-1)*PEsmpTOT+ip] = num + 1;
} else {SMPindex[(ic-1)*PEsmpTOT+ip] = num;}}}

SMPindex
ic=1

ic=2

ic=3

ic=5

12345678 12345678 12345678 12345678 12345678

SMPindex_new
12345 12345 12345 12345 12345 12345 12345 12345

SMPindex_new = (int *) allocate_vector(sizeof(int), NCOLORtot * PEsmpTOT + 1);
memset(SMPindex_new, 0, sizeof(int)*(NCOLORtot*PEsmpTOT+1));
for(ic=1; ic<=NCOLORtot; ic++) {for(ip=1; ip<=PEsmpTOT; ip++) {j1 =(ic-1)*PEsmpTOT + ip;
j0 =j1-1;
SMPindex_new[(ip-1)*NCOLORtot+ic] = SMPindex[j1];SMPindex[j1] = SMPindex[j0] + SMPindex[j1];}}for(ip=1; ip<=PEsmpTOT; ip++) {for(ic=1; ic<=NCOLORtot; ic++) {j1 = (ip-1) * NCOLORtot + ic;
j0 = j1 - 1;
SMPindex_new[j1] += SMPindex_new[j0];
}
}

ic=4

行列ベクトル積の計算法#pragma omp parallel for private(ip,i)for(ip=0; ip<PEsmpTOT; ip++) {for(i=SMPindexG[ip]; i<SMPindexG[ip+1]; i++) {VAL = D[i] * W[P][i];for(j=indexL[i]; j<indexL[i+1]; j++) {VAL += AL[j] * W[P][itemL[j]-1];}
for(j=indexU[i]; j<indexU[i+1]; j++) {VAL += AU[j] * W[P][itemU[j]-1];}
W[Q][i] = VAL;}}METHOD=0#pragma omp parallel for private (ip1,i,VAL,j)for(ip=0; ip<PEsmpTOT; ip++) {for(i=SMPindex[ip*NCOLORtot]; i<SMPindex[(ip+1)*NCOLORtot];i++){
VAL = D[i] * W[P][i];for(j=indexL[i]; j<indexL[i+1]; j++) {VAL += AL[j] * W[P][itemL[j]-1];}
for(j=indexU[i]; j<indexU[i+1]; j++) {VAL += AU[j] * W[P][itemU[j]-1];}
W[Q][i] = VAL;}METHOD=1
}

93

94

前進代入の計算法：色ループは外
for(ic=0; ic<NCOLORtot; ic++) {#pragma omp parallel for private (ip,ip1,i,VAL,j)for(ip=0; ip<PEsmpTOT; ip++) {ip1 = ic * PEsmpTOT + ip;
for(i=SMPindex[ip1]; i<SMPindex[ip1+1]; i++) {WVAL = W[Z][i];for(j=indexL[i]; j<indexL[i+1]; j++) {WVAL -= AL[j] * W[Z][itemL[j]-1];}
W[Z][i] = WVAL * W[DD][i];}}}Originalfor(ic=0; ic<NCOLORtot; ic++) {#pragma omp parallel for private (ip1,i,WVAL,j)for(ip=0; ip<PEsmpTOT; ip++) {ip1 = ip * NCOLORtot + ic;
for(i=SMPindex[ip1]; i<SMPindex[ip1+1]; i++) {WVAL = W[Z][i];for(j=indexL[i]; j<indexL[i+1]; j++) {WVAL -= AL[j] * W[Z][itemL[j]-1];}
W[Z][i] = WVAL * W[DD][i];
}
}
}

New

95

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

3
2
1
3
0

0
5
3
7
0

0
0
0
4
0

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

1 3

1 3 0

1 2 5

1 2 5

4 1 3

4 1 3

3 7 4

3 7 4

1 5

1 5 0

(a) CRS(b) ELL

96

Cases: 1283 meshes
並列化向
け色付け
Coloring
手法
src0

Case-1

reorder0 Case-2

CM-RCM

Case-3

ELL

First Touch
Data
Placement

Further
番号付け
Reordering

Coalesced（図 4（a））

NO無し

Sequential（図 4（b））

YES有り

color=1

ELL

Initial Vector

Initial Vector

color=1

CRS

Sequential

Coalesced
Coloring(5 colors)+Ordering

係数行列格納MatrixStorage
Format形式

color=2

color=2

color=3

color=4

color=5

color=3

color=4

color=512345678 12345678 12345678 12345678 12345678

Coloring(5 colors)+Orderingcolor=1

color=2

color=3

color=4

color=5

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

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

97

T2K/Tokyo

Cray XE6 (Hopper)

MemoryMemoryMemory

L2
L1

L3
L2 L2
L1 L1

L2
L1

L2
L1

L3
L2 L2
L1 L1

L2
L1

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

L1
L2

L1 L1
L2 L2
L3

L1
L2

L1
L2

L1 L1
L2 L2
L3

L1
L2

MemoryMemory

L2
L1

L2
L1

L3
L2 L2
L1 L1

L2
L1

L2
L1

L2
L1

L2
L1

L3
L2 L2
L1 L1

L2
L1

L2
L1

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

L1
L2

L1
L2

L1 L1
L2 L2
L3

L1
L2

L1
L2

L1
L2

L1
L2

L1 L1
L2 L2
L3

L1
L2

L1
L2

MemoryMemoryMemoryFujitsu FX10 (Oakleaf-FX)
Memory

L1

L1

L1

L1

L1

L1

L1

L2
L1 L1

L1

L1

L1

L1

L1

L1

L1

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

OMP-398

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

Iterations

400

350

300

250
1

10

COLOR#

100

1000

OMP-399

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

6.00

sec.

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

4.00

Case-1Case-22.00

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

1

10

100

COLOR#

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

1000

OMP-3

100

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

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

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

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

3.5E+00

3.0E+00

2.5E+00

2.0E+00

1.5E+00

1.0E+00

5.0E-01

0.0E+00
Thread 0

Thread 1

Thread 2

Thread 3

Thread 4

Thread 5

Thread 6

Thread 7

Thread 8

Thread 9

Thread 10

Thread 11

Thread 12

src0: CRS, Coalesced

Thread 13

Thread 14

Thread 15

OMP-3

101

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

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

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

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

3.5E+00

3.0E+00

2.5E+00

2.0E+00

1.5E+00

1.0E+00

5.0E-01

0.0E+00
Thread 0

Thread 1

Thread 2

Thread 3

Thread 4

Thread 5

Thread 6

Thread 7

Thread 8

Thread 9

Thread 10

Thread 11

Thread 12

Thread 13

reorder0: CRS, Sequential

Thread 14

Thread 15

OMP-3

102

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

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

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

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

3.5E+00

3.0E+00

2.5E+00

2.0E+00

1.5E+00

1.0E+00

5.0E-01

0.0E+00
Thread 0

Thread 1

Thread 2

Thread 3

Thread 4

Thread 5

Thread 6

Thread 7

Thread 8

Thread 9

Thread 10

Thread 11

Thread 12

src0: CRS, Coalesced

Thread 13

Thread 14

Thread 15

OMP-3

103

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

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

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

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

3.5E+00

3.0E+00

2.5E+00

2.0E+00

1.5E+00

1.0E+00

5.0E-01

0.0E+00
Thread 0

Thread 1

Thread 2

Thread 3

Thread 4

Thread 5

Thread 6

Thread 7

Thread 8

Thread 9

Thread 10

Thread 11

Thread 12

Thread 13

reorder0: CRS, Sequential

Thread 14

Thread 15

OMP-3104

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

reorder0
CASE-2
CRS+
Sequential

CASE-3
ELL+
Sequential

25.5 %

25.6 %

5.42 %

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

29.3 %

16.5 %

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

OMP-3105

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

Memory Access
SIMD（%）
Throughput（%）

Case-2
CRS

1.83×1011
1.83×1011

7.17
6.90

50.2
44.5

Case-3ELL

6.71×1010
5.96×1010

16.3
16.2

69.8
67.0

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

OMP-3106

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

Case-2

20.00

sec.

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

10.00

0.00
1

10

100

COLOR#

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

1000

OMP-3107

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

Cray
XE6

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

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

Case-15.44

1.71×10-2

6.37

2.22×10-2

Case-25.76

1.81×10-2

5.78

2.02×10-2

Case-33.90

1.23×10-2

3.61

1.26×10-2

Case-1

19.7

6.26×10-2

18.7

6.52×10-2

Case-27.45

2.34×10-2

7.40

2.58×10-2

Case-37.14

2.25×10-2

6.77

2.36×10-2

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

OMP-3

108

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

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

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

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

3.5E+00

3.0E+00

2.5E+00

2.0E+00

1.5E+00

1.0E+00

5.0E-01

0.0E+00
Thread 0

Thread 1

Thread 2

Thread 3

Thread 4

Thread 5

Thread 6

Thread 7

Thread 8

Thread 9

Thread 10

Thread 11

ELL, Sequential

Thread 12

Thread 13

Thread 14

Thread 15

OMP-3

109

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

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

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

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

3.5E+00

3.0E+00

2.5E+00

2.0E+00

1.5E+00

1.0E+00

5.0E-01

0.0E+00
Thread 0

Thread 1

Thread 2

Thread 3

Thread 4

Thread 5

Thread 6

Thread 7

Thread 8

Thread 9

Thread 10

Thread 11

ELL, Sequential

Thread 12

Thread 13

Thread 14

Thread 15

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

111

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

20.0
MC=2

MC=2

RCM(298)

16.0

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

Speed-Up

20.0

10.0

12.0

8.0

4.0
0.0

0.0
1

2

4

8

16

1

2

4

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

8

16

thread#Memory

L1

L1

L1

L1

L1

L1

L1

L2
L1 L1

L1

L1

L1

L1

L1

L1

L1

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

112

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

113

疎行列・密行列for(i=0;i<N;i++){Y[i] = Diag[i] * X[i];for(k=Index[i];k<Index[i+1];k++){Y[i] += AMat[k]*X[Item[k]];}}for(j=0;j<N;i++){Y[j]= 0.0;
for (i=0;i<N;i++){Y[j]+= A[j][i]*X[i];}}• “X” in RHS– 密行列：連続アクセス，キャッシュ有効利用– 疎行列：連続性は保証されず，キャッシュを有効に活用できず• より「memory-bound」

114

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

115

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

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

http://www.cs.virginia.edu/stream/

117

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

118

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

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

出力ファイル名

119

Results of Triad

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

L1

L1

L1

L1

L1

L1

L1

L2
L1 L1

L1

L1

L1

L1

L1

L1

L1

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

MB/sec.

Speed-up

1

8606.141.00

216918.81

1.97

434170.72

3.97

859505.92

6.91

1664714.32

7.52

C

120

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

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

122

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

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

123＊.lstの見方

124

SIMD情報125

自動並列化情報126

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

1
1
1
1
1
1
1
2

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

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

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

p
p
p
p

123
124
125
126
127
128
129
130

5
5
5
4
4
3
2
2

p
p
p
p
p
p

p

8v
8v
8v

p
p
p
p

4v
4v
4v

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

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

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

OMP-3129

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

OMP-3130

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

OMP-3131

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

132

略 称

FX10

名 称

Fujitsu SPARC64 IX
fx

動作周波数（GHz）

1.848

1.053

2.80

16 (16)60 (240)

10 (20)

16

240

10

DDR3

GDDR5

DDR3

236.5

1,010.9

224.0

32

8

64

85.1

320

59.7

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

MIC

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

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

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

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

キャッシュ構成

IvyB

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

133

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

CoalescedColoring(5 colors)+Orderingcolor=1

color=1

Initial Vector

color=2

color=2

color=3

color=4

color=5

color=3

color=4

color=512345678 12345678 12345678 12345678 12345678

SequentialColoring(5 colors)+Orderingcolor=1

Initial Vector

color=2

color=3

color=4

color=5

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

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

134

Numberin
g
AR-0

行列格納
形式

その他

CRS
行方向
（Rowwise）

AR-1
AR-2

外側ループ

Coalesced
Sliced ELL

間接メモリアクセス

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

AC-1
AC-2
BR-0

CRS

BR-1

行方向
（Rowwise）

Sequential
BC-1
BC-2

Sliced ELL

列方向（Columnwise）
ブロック化

135

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

k

i

136

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

137

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

k

400

• ELLの場合

380
360

Iterations

色数と計算時間（Solver）

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

340
320
300
280
260
1

10

100

sec.

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

1000

1000

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

139

AR-1
8.00

AR-2

sec.

AC-1
AC-2

6.00

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

AR-0

IvyB

AR-0

FX10

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

20

40

60

80

100

Color#

10.00

AR-0

MIC

AR-1
8.00

AR-1

AR-2

8.00

AR-2

AC-2

6.00

BR-0

AC-1
sec.

sec.

AC-1

AC-2

6.00

BR-0

BR-1
BC-1

4.00

BR-1
BC-1

4.00

BC-2
2.00

BC-2
2.00

0

20

40

60

Color#

80

100

0

20

40

60

Color#

80

100

140

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

FX10: 30 colors
MIC: 10 colors

6.00

IvyB: 40 colors

sec.

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

141

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

FX10: 30 colors
MIC: 10 colors

6.00

IvyB: 40 colors

sec.

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

142

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

FX10: 30 colors
MIC: 10 colors

6.00

IvyB: 40 colors

sec.

5.00

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

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

143

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

FX10: 30 colors
MIC: 10 colors

6.00

IvyB: 40 colors

sec.

5.00

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

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