科学 技 術 計 算 分 科 会

選 出

SS 研 HPC フォーラム 2008 より海外招待講演High-Productivity Languages for
Peta-Scale Computing

Hans P. Zima
Principal Scientist, Jet Propulsion Laboratory,
California Institute of Technology, Professor
Emeritus, University of Vienna, Austria

Contents
1. Introduction

High-Productivity Languages
forPeta-Scale Computing
Hans P. Zima
Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA

2. Emerging Architectures and Applications
3. Towards High Productivity Programming
4. The High Productivity Language Chapel
5. Alternative Language Approaches
6. Issues in Programming Environments

andUniversity of Vienna, Austria
zima@jpl.nasa.gov

7. Concluding Remarks

Fujitsu HPC Forum 2008
Tokyo, Japan, August 27th, 2008

High-Productivity

The Meaning of “High-Productivity”

Programming and Execution Models

High productivity”
 “
productivity” implies three properties:“High
1. human-centric: programming at a high level of abstraction
2. high-performance: providing “abstraction without guilt”
3. reliability
 Raising the level of abstraction is acceptable only iftarget code performance is not significantly reduced

Programming
Model

Semantics

conceptual view ofdata and control

Productivity
Modelrealizations

Execution
Model This relates to a broad range of topics:

0 language design
0 architecture- and application-adaptive compiler technology
0 operating and runtime systems
0 library design and optimization
0 intelligent tool development
0 fault tolerance

Programming S
Language
P
Programming
Language+
Directives

Contents

Library S
P

abstract machine

Command-line S
InterfaceP

S
PMulticore: An Emerging Technology
 The era of faster sequential processors is over——exponential growth
over
over—exponentialof frequency can no longer be maintained

1. Introduction
2. Emerging Architectures and Applications

0 CMOS manufacturing technology approaches physical limits
0 power wall, memory wall, instruction-level parallelism (ILP) wall
0 Moore’s Law still in force: number of transistors on chip increasing Multicore technology provides continued performance growth

3. Towards High Productivity Programming
4. The High Productivity Language Chapel
5. Alternative Language Approaches
6. Issues in Programming Environments
7. Concluding Remarks

0 a multicore chip is a single chip with two or more independent processing units
0 improvements by multiple cores on a chip rather than higher frequency
0 on-chip resource sharing for cost and performance benefits Multicore systems have been produced since 2000

0 IBM Power 4; Sun Niagara; AMD Opteron; Intel Xeon;…
0 Quadcore systems by AMD, Intel0 IBM/Sony/Toshiba: Cell Broadband EnginePower Processor (PPE) and 8 Synergistic PEs (SPEs)
peak 100 GF double precision (IBM Power XCEll 8i)

 1000 cores on a chip possible with 30nm technology
 “Manycore”
Manycore” chips are already emerging …
“Manycore”

Future Multicore Architectures:
From 10s to 100s of Processors on a Chip

Multicore Systems are Not Just Small
SMPs or MPPs

 Tile64 (Tilera(Tilera Corporation, 2007)

0 64 identical cores, arranged in an 8X8 grid
0 iMesh on-chip network, 27 Tb/sec bandwidth
0 170-300mW per core; 600 MHz – 1 GHz
0 192 GOPS (32 bit)—about 10 GOPS/Watt-chip inter-core bandwidth is much larger than for a
 IntraIntra
inter
Intra-chip
inter-core
typical parallel machine (SMP or MPP)
-chip inter-core latencies are much smaller IntraIntra
inter
Intra-chip
inter-core Kilocore 1025 (Rapport Inc. and IBM, 2008)

0 Power PC and 1024 8-bit processing elements
0 125 MHz per processing element
0 32X32 “stripes” dedicated to different tasks

 Multicore systems can offer lightweight synchronization-based synchronization is unacceptable: transactional
 LockLock
Lock-based

 512-core SING chip (Alchip
512(Alchip Technologies, 2008)
512-corememory and full/empty bits (Cray MTA) are alternatives

0 for GRAPE-DR, a Japanese supercomputer project-In ProcessingProcessing
In-Memory (PIM) technology offers additional
Processing-In-Memory

 80-core 2 TF research chip from Intel (2011)80
80-coremethods for exploitation of locality

0 2D on-chip mesh network for message passing
0 1.01 TF (3.16 GHz); 62W power—16 GOPS/Watt
0 Note: ASCI Red (1996): first machine to reach 1 TF

4,510 Intel Pentium Pro nodes (200 MHz)500 KW for the machine + 500 KW for cooling of the room

From Eniac (1946) …

Top 500 Performance Development
11.7 PF/s

10/

IBM Roadrunner
1.02 PF/

1 Pflop/s
IBM BlueGene/L
NEC Earth Simulator

100 Tflop/s

9.0 TF/SUM
IBM ASCI W hite
1.17 TF/sIntel ASCI Red

10 Tflop/s
59.7 GF/s
Fujitsu 'NW T'

6-8 years#1

1 Tflop/s

My Laptop#500
100 Gflop/s

2008

2007

2006

2005

2004

2003

2002

2001

2000

1999

1998

1997

1996

1995

1994

1 Gflop/s

1993

0.4 GF/s

10 Gflop/s

100 M flop/s

103 OPS…to LANL Roadrunner: Top 500 #1

Applications
 HPC has become the third pillar of science andengineering, in addition to theory andexperiment Traditional application areas include:

Cell Blade

1,026 TF=1015 OPS
The first machine reaching
Peta-scale performance
17 clusters, each with 192 nodes
Each node contains Opteron and 4 Cells
12,960 Cell chips (100 GF double precision)
Each Cell contains a PowerPC and 8 SPEs
6,948 dual-core Opterons
Total: 122,400 cores

0 DNA Analysis
0 Drug Design
0 Medicine
0 Aerospace
0 Manufacturing
0 Weather Forecasting and Climate Research New architectures facilitate new applications:

0 Graph Traversals
0 Dynamic Programming
0…
0 Backtrack Branch & BoundUC Berkeley’s“Dwarfs”

Contents

High-Level Sequential Languages
The designers of the very first high level programming language were
aware that their success depended on the target code performance:

1. Introduction
2. Emerging Architectures and Applications
3. Towards High Productivity Programming
4. The High Productivity Language Chapel
5. Alternative Language Approaches
6. Issues in Programming Environments

John Backus (1957): “… It was our belief that if FORTRAN … were to translate
any reasonable scientific source program into an object program only half
as fast as its hand-coded counterpart, then acceptance of our system
would be in serious danger …”

High-level algorithmic languages became generally
accepted standards for sequential programming since
their advantages outweighed any performance drawbacks

7. Concluding Remarks

For programming of HPC systems
no similar development took place

Programming Paradigm for
MPPs and Clusters: MPI is State-of-the-Art
The MPI Message-Passing Model
Message
Message-Passing
 widely adopted portable standard for full control of communication
communication
 relatively simple execution model can achieve good performance on commodity clusters

Drawbacks of the MPI Model

MPI vs HPF:
An Example for Locality Management (Jacobi Relaxation)

Sequential Codereal, allocatable A(:, : ), B(:, : )…
do while ( .not. converged )do J=1,N
do I=1,N
B(I,J)=0.25(A(I-1,J)+A(I+1,J)+A(I,J-1)+A(I,J+1))

 low-level paradigm: ““thethe assembly language of parallel programming”
low
programming”
low-level lack of separation between algorithm and communication management
management
 complex, difficult-todifficult
to-change communication structuresdifficult-to-change scalability to peta-scale questionablepeta
peta-scale

enddo
enddo
A(1:N,1:N)=B…

dependence pattern

enddo

Alternatives to MPI have been proposed

Parallelization Based on Data Distribution

 automatic vectorization and parallelization
 libraries for one-sided communication (SHMEM, ARMCI, GASNet)one
GASNet)one-sided High Performance Fortran (HPF), PGAS languages, OpenMP,
OpenMP, etc.

Boundary Exchange in Overlap Regions
Pk-1

Pk

Pk+1

Let A and B be partitioned into blocks of columns mapped todifferent processors. All processors can work concurrently on
their local data, but an exchange must take place at segment
boundaries after each iteration…

Message Passing Approach

do while ( .not. converged )do J=1,M ! Number of local columns
do I=1,N

local view of data, local control,
explicit two-sided communication

After iteration:
Data Exchange…

HPF Approach…

global view of data, global control,
compiler-generated communication

initialize MPI
do while (.not. converged)
do J=1,M
do I=1,N
B(I,J) = 0.25 * (A(I-1,J)+A(I+1,J)+
A(I,J-1)+A(I,J+1))
end doend doA(1:N,1:N) = B(1:N,1:N)

local computation

global computation
do while (.not. converged)
do J=1,N
do I=1,N
B(I,J) = 0.25 * (A(I-1,J)+A(I+1,J)+
A(I,J-1)+A(I,J+1))
end doend doA(1:N,1:N) = B(1:N,1:N)

communication

Processor Pk reads:• rightmost column of Pk-1
• leftmost column of Pk+1.
Processor Pk copies:• its leftmost column to Pk-1
• its rightmost column to Pk+1.
halo regionsPs

The Key Idea ofHigh Performance Fortran (HPF)! purely local operation in each iteration:

B(I,J)=0.25(A(I-1,J)+A(I+1,J)+
A(I,J-1)+A(I,J+1))
enddo
enddo
…

…

P1 P2

if (MOD(myrank,2) .eq. 1) thencall MPI_SEND(B(1,1),N,…,myrank-1,..)
call MPI_RCV(A(1,0),N,…,myrank-1,..)
if (myrank .lt. s-1) thencall MPI_SEND(B(1,M),N,…,myrank+1,..)
call MPI_RCV(A(1,M+1),N,…,myrank+1,..)
endif
else …
…

data distribution
processors P(NUMBER_OF_PROCESSORS)
distribute(*,BLOCK) onto P :: A, B

communicationcompiler-generated

K. Kennedy, C. Koelbel,
Koelbel, and H. Zima: The Rise and Fall of High Performance Fortran: An Historical Object
Object Lesson
Proc. History of Programming Languages III (HOPL III), San Diego,
Diego, June 2007

Fortran+MPI Communicationfor 3D 27-point Stencil (NAS MG rprj3)

Example: Sweep Over Unstructured Mesh in HPF!HPF$ PROCESSORS P(NUMBER_OF_PROCESSORS())
TYPE NODE! type of a node in the unstructured grid…
REAL::V1, V2! flow variables
END TYPE NODE

subroutine
subroutine comm3(u,n1,n2,n3,kk)
comm3(u,n1,n2,n3,kk)
use
use caf_intrinsics
caf_intrinsics

doi=1,nm2
buff(i,buff_id) = 0.0D0

buff_id = 3 + dir

enddo

indx = 0if( axis .eq.

else if( dir .eq. +1 ) thendir = +1dobuff(1:buff_len,buff_id+1)[nbr(axis,dir,k)] =
>

do

enddo
enddo
2,i3)

enddo
enddo

integer
integer n1,
n1, n2,n2, n3,n3, kkkk
doubledouble precision
precision u(n1,n2,n3)u(n1,n2,n3)
integer
integer axis
axis

do

i2=2,n2i2=2,n2-1

buff_id = 3 + dir

indx = indx + 1

buff_len = nm2

if(
dead(kk)) )thenif( .not.
.not. dead(kkdead(kk))thendodo axisaxis == 1,
1, 33

do

enddo

else
else
do
do axisaxis == 1,
1, 33
call
call sync_all()sync_all()if( axis .eq.

endif

i1=1,n1buff_len = buff_len + 1
buff(buff_len,buff(buff_len,

dir = +1

endif

do

if( axis .eq.

2 )thenif( dir .eq. -1 )then

buff_id = 2 + dir

i1=1,n1
indx = indx + 1

buff_len = 0u(i1,n2,i3) = buff(indxbuff(indx,, buff_id )
enddo

do

buff(1:buff_len,buff_id+1)[nbr(axis,dir,k)] =i3=2,n3-1
i3=2,n3dobuff(1:buff_len,buff_id)

if( axis .eq.
doi1=1,n1do

enddo

else if( dir .eq. +1 ) then3 )theni2=1,n2do

i1=1,n1
indx = indx + 1

enddo

3 )thenif( dir .eq. -1 )thenif( axis .eq.

buff(buff_len,buff(buff_len, buff_id ) = u( n1n1-1, i2,i3)

enddo
enddo

if( axis .eq.

enddo
endif

i2=2,n2-1
i2=2,n2buff_len = buff_len + 1

u(i1,n2,i3) = buff(indxbuff(indx,, buff_id )

endif
endif

1 )theni3=2,n3-1
i3=2,n3do

indx = indx + 1

u(i1,i2,n3) = buff(indxbuff(indx,, buff_id )

endif

call
call sync_all()sync_all()

enddo
do

do

i2=1,n2doi3=2,n3-1
i3=2,n3do

i1=1,n1

if( axis .eq.
doi1=1,n1

enddo
endif

i1=1,n1buff_len = buff_len + 1

u(i1,1,i3) = buff(indxbuff(indx,, buff_id )buff(buff_len,buff(buff_len, buff_id ) = u( i1,i2,2)buff(buff_len,buff(buff_len
,

enddo

enddo

buff_id )= u( i1,n2i1,n2-1,i3)

buff_id = 3 + dir

enddo
subroutine
subroutine give3(
give3( axis,
axis, dir,dir, u,u, n1,
n1, n2,n2, n3,n3, kk ))use
use caf_intrinsics
caf_intrinsics

buff(1:buff_len,buff_id+1)[nbr(axis,dir,k)] =>buff(1:buff_len,buff_id)

indx = 0

endif

endif
endif

if( axis .eq.
else if( dir .eq. +1 ) thenimplicit
implicit nonenone

if( axis .eq.

do3 )thendo

integer
integer axis,
axis, dir,dir, n1,
n1, n2,n2, n3,n3, k,k, ierrierrdoubledouble precision
precision u(u( n1,
n1, n2,n2, n3
n3 ))

enddo

do

i1=1,n1dou(i1,i2,n3) = buff(indxbuff(indx,, buff_id )u(1,i2,i3) = buff(indxbuff(indx,, buff_id )

endif

endif

enddo

enddo

if( axis .eq.

dir = -1

enddo

do

buff(1:buff_len,buff_id+1)[nbr(axis,dir,k)] =>buff_idbuff_id == 22 ++ dir
dirbuff_lenbuff_len == 00

buff(1:buff_len,buff_id)

else if( dir .eq. +1 ) thenbuff_id = 2 + dir

i1=1,n1do1 )thenenddo

i3=2,n3-1
i3=2,n3dou(i1,i2,1) = buff(indxbuff(indx,, buff_id )doi3=2,n3-1
do i3=2,n3
i3=2,n3-1doi2=2,n2-1
do i2=2,n2i2=2,n2-1

enddo

i2=2,n2-1
i2=2,n2-

endif

buff_len = buff_len + 1
buff(buff_len,buff_id ) = u( 2,enddo

i2,i3)

enddo

buff_lenbuff_len == buff_len
buff_len ++ 11
buff(buff_len,buff_id
buff(buff_len,buff_id )) == u(u( 2,2, i2,i3)
i2,i3)

subroutine take3( axis, dir, u, n1, n2, n3 )use caf_intrinsics

enddo
enddo
enddo
enddo

i1=1,n1u(i1,1,i3) = buff(indxbuff(indx,, buff_id )if( axis .eq.

indx = indx + 1

enddo

end

indx = indx + 1

i2=1,n2
doreturn2 )theni3=2,n3-1
i3=2,n3do

buff_len = 0
do

endif
endif

if(
if( axis
axis .eq.
.eq. 11 )then
)then
if(
1 )thenif( dir
dir .eq.
.eq. --1)theni2=2,n2-1
i2=2,n2indx = indx + 1

enddo

enddo

indx = indx + 1

buff(buff_len,buff(buff_len, buff_id ) = u( i1,i2,n3i1,i2,n3-1)

enddo

enddo

i1=1,n11 )theni3=2,n3-1
i3=2,n3dobuff(buff_len,buff(buff_len, buff_id ) = u( i1,i2,n3i1,i2,n3-1)

i2=1,n2dobuff_len = buff_len + 1

integer
integer i3,
i3, i2,i2, i1,
i1, buff_len,buff_id
buff_len,buff_id

if( axis .eq.

i1=1,n1buff_len = buff_len + 1

i2=1,n2include''cafnpb.h'
cafnpb.h''include 'cafnpb.hinclude''globals.h'
globals.h''include 'globals.h3 )theni2=1,n2do

if( dir .eq. -1 )thendodir = +1

enddo

enddo

enddo

2 )theni3=2,n3-1
i3=2,n3do

indx = indx + 1

buff_len = buff_len + 1

returnreturnendend

if( axis .eq.
do

enddo

endif

endif

returnif( axis .eq.

3 )theni2=1,n2do

endif

i1=1,n1
indx = indx + 1
u(i1,i2,1) = buff(indxbuff(indx,, buff_id )buff(1:buff_len,buff_id+1)[nbr(axis,dir,k)]buff(1:buff_len,buff_id+1)[nbr(axis,dir,k)] ==buff(1:buff_len,buff_id)
buff(1:buff_len,buff_id)

do

end

implicit none

include 'globals.h'globals.h'
'

enddo
enddo

i1=1,n1

endif

buff_len = buff_len + 1
buff(buff_len,, buff_id ) = u( i1,buff(buff_len

subroutine comm1p( axis, u, n1, n2, n3, kk )

2,i3)

enddo

use caf_intrinsics

returnend

enddo

integer axis, dir, n1, n2, n3doi3=2,n3-1
do i3=2,n3
i3=2,n3-1doi2=2,n2-1
do i2=2,n2i2=2,n2-12 )theni3=2,n3-1
i3=2,n3do

include 'cafnpb.h'cafnpb.h''

else
else if(
if( dir
dir .eq.
.eq. +1+1 )) then
then

endif

double precision u( n1, n2, n3 )implicit none

buff_lenbuff_len == buff_len
buff_len ++ 11
buff(buff_len
buff(buff_len,, buff_id
n1-1, i2,i3)
buff(buff_len,
buff_id )) == u(u( n1n1-1,
i2,i3)>>

2 )theni3=2,n3-1
i3=2,n3do

buff_id )= u( i1,n2i1,n2-1,i3)

enddo
enddo

>

enddo
enddo
call
call zero3(u,n1,n2,n3)
zero3(u,n1,n2,n3)
endif
endif

>>

1 )then

enddo
endif

enddo

i3=2,n3-1
i3=2,n3do

call
call sync_all()sync_all()
1, u,
call
call take3(
take3( axis,
axis, --1,
u, n1,
n1, n2,n2, n3n3 ))callcall take3(
take3( axis,
axis, +1,
+1, u,
u, n1,
n1, n2,n2, n3
n3 ))
else
else
call
call comm1p(
comm1p( axis,
axis, u,u, n1,
n1, n2,n2, n3,n3, kk
kk ))
endif
endif
enddo
enddo

i=1,nm2
buff(i,buff_id) = 0.0D0

enddo
do

i2=2,n2-1
i2=2,n2indx = indx + 1
u(n1,i2,i3) = buff(indxbuff(indx,, buff_id )enddo

u(1,i2,i3) = buff(indxbuff(indx,, buff_id )else if( dir .eq. +1 ) theni3=2,n3-1
i3=2,n3do

i3=2,n3-1
i3=2,n3dobuff(1:buff_len,buff_id)

call
call give3(
give3( axis,
axis, +1,
+1, u,
u, n1,
n1, n2,n2, n3,n3, kk
kk ))
call
1, u,
call give3(
give3( axis,
axis, --1,
u, n1,
n1, n2,n2, n3,n3, kkkk ))Chapel 3D NAS MG Stencil rprj3

u(n1,i2,i3) = buff(indxbuff(indx,, buff_id )i1=1,n1buff_len = buff_len + 1
buff(buff_len,buff(buff_len, buff_id ) = u( i1,

include''cafnpb.h'
cafnpb.h''include 'cafnpb.hinclude''globals.h'
globals.h''include 'globals.hif(
if( nprocs
nprocs .ne.
.ne. 1)
1) thenthencall
call sync_all()sync_all()TYPE(NODE), ALLOCATABLE::GRID(:)
REAL, ALLOCATABLE::EDGE(:,2)
INTEGER, ALLOCATABLE::MAP(:)! mapping array!HPF$ DYNAMIC, DISTRIBUTE(BLOCK)::GRID!HPF$ DYNAMIC, DISTRIBUTE(BLOCK,*)::EDGE!HPF$ DISTRIBUTE(BLOCK)::MAP…
! Read parameters; allocate GRID, MAP; initialize GRID, M
CALL GRID_PARTITIONER(GRID,MAP)!HPF$ REDISTRIBUTE GRID(INDIRECT(MAP))
ALLOCATE(EDGE(M,2))! Initialize and realign EDGE with GRID
! Sweep over edges of the grid:!HPF$ INDEPENDENT,ON HOME(EDGE(J,1)),NEW(N1,N2,DELTAV),REDUCTION(V2)DO J=1,M
N1=EDGE(J,1), N2=EDGE(J,2)…
DELTAV=F(V1(N1),V1(N2))
V2(N1)=V2(N1)-DELTAV
V2(N2)=V2(N2)+DELTAV
ENDDO

i3=2,n3-1
i3=2,n3do

implicit
implicit nonenone

if( axis .eq.

integer buff_id, indxdoinclude 'cafnpb.h'cafnpb.h'
'

enddo
enddo
enddo
enddo

integer i3, i2, i1

buff_id = 3 + dir

integer axis, dir, n1, n2, n3buff(1:buff_len,buff_id+1)[nbr(axis,dir,k)]buff(1:buff_len,buff_id+1)[nbr(axis,dir,k)] ==buff(1:buff_len,buff_id)
buff(1:buff_len,buff_id)

indx = 0

double precision u( n1, n2, n3 )3 )theni2=1,n2doinclude 'globals.h'globals.h''

i1=1,n1buff_len = buff_len + 1
buff(buff_len,buff(buff_len, buff_id ) = u( i1,i2,2)

enddo
enddo

if( axis .eq.
endif
endif
endif
endif

1 )thenif( dir .eq. -1 )then

integer i3, i2, i1, buff_len,buff_id

endif

integer i, kk
kk,, indxdodo
if(
if( axis
axis .eq.
.eq. 22 )then
)then
if(
1 )thenif( dir
dir .eq.
.eq. --1)theni3=2,n3-1
i3=2,n3do

dir = -1i=1,nm2
buff(i,4) = buff(i,3)
buff(i,2) = buff(i,1)

i2=2,n2-1
i2=2,n2indx = indx + 1

buff_id = 3 + dir

enddo

buff_len = nm2
dir = -1

※最終ページに
拡大印刷ありProductivity Challenges for Peta-Scale Systems
 Large-scale hierarchical architectural parallelism

0tens of thousands to hundreds of thousands of processors
0component failures may occur frequently

function rprj3(S,R) {const Stencil: domain(3) = [[-1..1, -1..1, -1..1],// 2727-points

param
coeff:: domaindomain(1)(1) == [0..3];
param coeff
coeff:
domain(1)[0..3]; //// forfor 44 uniqueunique weightweight values
valuesparam
domain(3)(3) == [[[-1..1,-1..1, --1..1,
1..1, --1..1];
1..1]; //27-points
param Stencil:
Stencil: domain
domain(3)// 2727-points

w: [0..3]real = (/0.5, 0.25, 0.125, 0.0625/),functionfunction rprj3(S,rprj3(S, R)
R) {{
param
[[coeff]
coeff]] floatparam w:
w: [coeff
float == (/0.5,(/0.5, 0.25,
0.25, 0.125,
0.125, 0.0625/);
0.0625/);// weights
weights

param
param w3d:
w3d: [(i,j,k)[(i,j,k) inin Stencil]
Stencil] float
float== w((i!=0)w((i!=0) ++ (j!=0)(j!=0) ++ (k!=0));(k!=0));

w3d: [(i,j,k) in Stencil] = w((i!=0) + (j!=0) + (k!=0));
const
const SDSD == S.Domain,
S.Domain,
Rstr
Rstr == R.stride;R.stride;

SS == [ijk[[ijk
ijk inin SD]
SD] sumsum reduce
reduce [off[off inin Stencil]
Stencil](w3d(off)
Rstr*off));*off));(w3d(off) ** R(ijk
R(ijk ++ Rstr
Rstr*off));}}

forall ijk in S.domain do
S(ijk)S(ijk) = sum reduce [off in Stencil] (w3d(off) * R(ijk + R.stride*off));
R.stride*off));}

 Extreme non-uniformity in data access Applications: large, complex, and long-lived

0multi-disciplinary, multi-language, multi-paradigm
0dynamic, irregular, and adaptive
0survive many hardware generations Î portability is important How to exploit the parallelism and locality provided by

the architecture?
0automatic parallelization and locality management are not
powerful enough to provide a general efficient solution

0explicit support for control of parallelism and locality must be
provided by the programming model and the language

Parallel Programming Models
 Fragmented Models
0processor-centric view: code written from the viewpoint of single
threads

0local view of data segments

 Single Program Multiple Data (SPMD) Model
0special class of fragmented model
0single program executed in multiple instances Global-view Models
0global view of data and computation
 burden of partitioning shifts to compiler/runtime
 user may guide this process via language constructs

 Locality-aware Models
0features for mapping data and/or control to the architecture

Languages for High Performance Computing
 HPF Language Family

0predecessors: CM-Fortran, Fortran D, Vienna Fortran
0High Performance Fortran (HPF): HPF-1 (1993); HPF-2(1997)
0successors: HPF+, HPF/JA
 OpenMP Partitioned Global Address Space (PGAS) Languages

0Co-Array Fortran
0UPC
0Titanium
 High-Productivity Languages developed in the HPCS Program

0Chapel
0X10
0Fortress
 Domain-Specific Languages and Abstractions

HPCS Languages

Contents

global view of data, global control

 High-Productivity Computing Systems (HPCS) is a DARPA-sponsored

1. Introduction

program for the development of peta-scale architectures (2002-2010)

2. Emerging Architectures and Applications
3. Towards High Productivity Programming
4. The High Productivity Language Chapel
5. Alternative Language Approaches

 HPCS Languages

0 Chapel (Cascade Project, led by Cray Inc.)
0 X10(PERCS Project, led by IBM)
0 Fortress (HERO Project [until 2006], led by Sun Microsystems)
 These are new, memory-managed, object-oriented languages

0global view of data and computation Î generally no distinction
between local and remote data access in the source code

6. Issues in Programming Environments

0support for explicit data and task parallelism
0explicit locality management
0Chapel is unique in that it provides user-defined data distributions

7. Concluding Remarks

Chapel Language Concepts

Example: Jacobi Relaxation in Chapel

http://chapel.cs.washington.edu
 Explicit high-level control of parallelism

const L:[1..p,1..q] locale = reshape(Locales);0 data parallelism

const n= …, epsilon= …;
const DD:domain(2)=[0..n+1,0..n+1] distributed(block,block)on L;
D: subdomain(DD) = [1..n, 1..n];var delta: real;
var A, Temp: [DD] real; /*array declarations over domain DD */domains, arrays, indices: support distributed data aggregates forall loops and iterators: express data parallel computations


0 task parallelism


cobegin statements: specify task parallel computations
synchronization variables, atomic sections

 Explicit high-level control of locality

A(0,1..n) = 1.0;

0 “locales”: abstract units of locality
0 data distributions: map data domains to sets of locales
0 on clauses: map execution components to sets of locales Close relationship to mainstream languages

0 object-oriented
0 type inference and generic programming
0 modules for Programming-in-the-Large

do {forall (i,j) in D { /* parallel iteration over domain D */Temp(i,j) = (A(i-1,j)+A(i+1,j)+A(i,j-1)+A(i,j+1))/4.0;
delta = max reduce abs(A(D) – Temp(D));A(D) = Temp(D);} while (delta > epsilon);
writeln(A);Note: Some of the features discussed in the following have the status of research
proposals and are currently not part of the official Chapel language specification

Task Parallelism in Chapel

Example: Jacobi Relaxation in Chapel
const L:[1..p,1..q] locale = reshape(Locales);const n= …, epsilon= …;
const DD:domain(2)…distributed(block,block on L;
D: subdomain(DD) = [1..n, 1..n];var delta: real;
var A, Temp: [DD] real; Task Creation
cobegin { S1,…Sn}

executes the Si in parallel (i = 1,…n)

Locale Grid L

 Task Localizationon L(i,j) do f(A(i,j))

A(0,1..n) = 1.0;
do {forall (i,j) in D {Temp(i,j) = (A(i-1,j)+A(i+1,j)+A(i,j-1)+A(i,j+1))/4.0;
delta = max reduce abs(A(D) – Temp(D));A(D) = Temp(D);Key Features} while (delta > epsilon);•global view of data/control•explicit parallelism (forall)
writeln(A);•high-level locality control•NO explicit communication•NO local/remote distinction
in source code

executes f(A(i,j) on locale L(i,j)

 Task Synchronization- atomic sections-

sync variables-

single-assignment variables

Chapel’s Framework for
User-Defined Distributions

Aspects of Locality
Locale: an abstract Provides functionality for:0distributing index sets across locales
0arranging data within a locale
0defining specialized distribution libraries

unit of locality

Locale Set

P
M
P
M
P
M

P
M
P
M

P
M
P
M
P
M

P
M
P
M
P
M

 This capability is in its effect similar to function specification

0unstructured meshes
0multi-block problems
0multi-grid problems
0distributed sparse matrices

distributework

distribute
data

work

domainalign data with work
(affinity)align data
domainUser-Defined Distributions:Global Mapping

Locality Control in Chapel: Basic Concepts
 Domain: first class entity
0 components: index set, distribution, associated arrays, iterators/* declaration of distribution classes MyC and MyB: */class MyC: Distribution {const z:int;
const ntl:int;

 Array—Mapping from a Domain to a Set of Variables Framework for User-Defined Distributions: three levels
1. naï
ïve use of a predefined library distribution (block, cyclic, indirect,
rect,…
…)
na
indi
naïve
indirect,…)
2. specification of a distribution by
global mapping: index set Æ localesinterface for the definition of mapping, distribution segments, iterators
system-provided default functionality can be overridden by user

3. specification of a distribution by global mapping andlayout mapping: index set Æ locale data space High-Level Control of Communication
0 user-defined specification of halos; communication assertions

function map(i:index(source)):locale {return Locales(mod(ceil(i/z-1)+1,ntl));}class MyB: Distribution {var bl:int = ...;
function map(i: index(source)):locale {return Locales(ceil(i/bl));
}

2

j
1

1

3
3

2

4
4

5

6

5

6

7
7

8

9
9

8

1

2

3

4

i

2

7

bw = 3 (bandwidth)11

p=4 (number of locales)

12

Distribution—global map:

13

Blocks of bw diagonals are
cyclically mapped to locales

15

4

16
17

1

18

8
9

10

14

3

5

6

Diagonal A/d = { A(i,j) | d=i+j }2/* global mapping for MyC *//* block length *//* global mapping for MyB */}/* use of distribution classes MyC and MyB in declarations: */const D1C: domain(1) distributed(MyC(z=100))=1..n1;
const D1B: domain(1) distributed(MyB) on Locales(1..num_locales/10)=1..n1;var
A1: [D1C] real;var
A2: [D1B] real;

ExampleMatrix-Vector Multiplication (sparse CRS)

Example: Banded Distributiond/* block size *//* number of target locales*/

0 53
0 0
19 0
0 0
0 0
0 0
0 0

0 0 0
0 0
0
0 0 0
0 0
0
0 17 00 0 930 0 0

0 0 0
0 21 00 0 1672 0 00 0 0
0 0 0
0 13 0R0
D0 C0 ________ ____
53
19
17
93

2
1
4
5

Layout:Each diagonal is represented
as a one-dimensional dense
array. Arrays in a locale arereferenced by a pointer array0 0 0 00 23 69 0
27 0 0 11

44
37
0

0 0 190 0 00 64 0

1
2
2
3
3
4
5
5

R01
D01 C01 ________ ____

53
21
19
16
17
72
93
13

27
18
46
57

1
21
2
3
34
4
54
5

D20 C2 R2____ ____ ____

D3 C3____ ____

R3

23
53
19
69
17
27
93
11

44
19
37
64

1
3
4
5

2
3
1
4

1
1
3
5

5
8
5
7

____

Contents

User-Defined Halos User-Defined Specification of halo (ghost cells) Compiler/Runtime System
0 allocates local images of remote data
0 defines mapping between remote objects and their images

1. Introduction
2. Emerging Architectures and Applications
3. Towards High Productivity Programming Halo Management
0 update
0 flush

4. The High Productivity Language Chapel
5. Alternative Language Approaches
distribution
segment

6. Issues in Programming Environments
7. Concluding Remarks

PGAS Language Overview

Domain-Specific Languages

Support for global view of data, but local control

 Partitioned Global Address Space (PGAS) languages

are based on the SPMD model Providing a shared-memory, global view, of data, General-purpose languages are limited in their ability toaccommodate the abstractions of a scientific domain Domain-specific languages provide abstractions tailoredcombined with support for locality

to a specific domain

0 global address space is logically partitioned, with each portion

0narrowing of the semantic gap between the programming

mapped to a processor

0 single-sided shared-memory communication (instead of MPI-stylemessage passing)language and the application domain

0separation of domain expertise from parallelization and resource
management

0 in general, local and remote references distinguished in the source
code

0 implemented via one-sided communication libraries (e.g., GASNet) Local control of execution via processor-centric view
 Main representatives: Co-Array Fortran (CAF), Unified

 Domain-specific knowledge can be used to improve

program analysis and support V&V and fault tolerance.
 Telescoping supports the automatic generation of

domain-specific languages by generating specialized,
optimized versions of libraries

Parallel C (UPC), Titanium

Contents

Issues in Programming Environments
 Legacy Code Migration

1. Introduction
2. Emerging Architectures and Applications

 (Semi) Automatic Tuning

3. Towards High Productivity Programming

0closed loop adaptive control: measurement, decision-making,4. The High Productivity Language Chapel

0information exposure: users, compilers, runtime systems
0learning from experience: databases, data mining, reasoning

5. Alternative Language Approaches
6. Issues in Programming Environments
7. Concluding Remarks

actuationsystems Fault Tolerance
0massive parallelism poses new reliability problems
0fault anticipation, detection, localization, analysis, and recovery

Legacy Code Migration

Offline Performance Tuning Rewriting Legacy Codes
0preservation of intellectual content
0performance portability: exploit new hardware and new algorithms
0code size may preclude complete rewrite: incremental porting Language, compiler, tool, and runtime support
0(semi) automatic tools for migrating code
0translation of performance-critical sections requires highly-

Transformation
systemsourceprogramactuators

Parallelizing
compiler

KnowledgeBase static analysis, using domain and/or architecture-specific knowledge

Expert System
Advisormonitoring
sensors

sophisticated software for automatic adaptation
 reverse engineering of the original program

restructuringcommands

instrumented
targetprogram pattern matching and concept comprehension
 optimizing code generation guided by the target architecture

HPCS system
Application Execution

end oftuning cycle

High-Capability On-Board System:
An Example

Autonomy and Fault Tolerance forHigh-Performance Space-Borne Computing

Fault-Tolerant High-Capability Computational Subsystem

Spacecraft
ControlComputer

High-Performance
Computing System (HPCS)

SystemController(SYSC)

Multi-coreComputeP P PEngineM M M
P P P P
ClusterM M M MP P P P

(SCC)

Mars Sample Return
Neptune Triton
Explorer
Communication
Subsystem(COMM)

Intelligent
Mass Data
Storage(IMDS)

Interface
fabricIntelligent
Processor
In
MemoryData
ServerM M M M

Instrument Interface

Instruments
EuropaExplorer

Titan Explorer…

Europa Astrobiology
Laboratory

Future missions need Autonomy and High-Capability On-Board Computing:this can be accomplished by extending traditional spacecraft architectures

A Framework for Introspection

Fault tolerance for the HPCS component—protecting it
against transient faults—is a specific concern for this model

EARTH

An Introspection Module (IM)Introspection…
Introspection…

Introspection System provides dynamic monitoring, analysis, and feedback,sensors

enabling system to become self-aware:
self-aware and contextcontext
self-aware
context-aware:0monitoring execution behavior
0reasoning about its internal state0changing the system or system state when necessary exploits adaptively the available threads.
.
.

Application

Monitoring

Knowledge
BaseAnalysis

SystemKnowledge

Recovery can be applied to different scenarios, including:0fault tolerance
0performance tuning
0power management
0behavior analysis
0intrusion detection

Inference Engine(SHINE).
.
.

For hierarchically structured architectures,a corresponding hierarchy of introspection
modules will be constructedactuators

DomainKnowledge…

Prognostics

Application
Knowledge

Example
BRD Distribution with CRS Layout

Conclusion
 HPCS languages constitute an important step towards high-productivity

programming for massively parallel peta-scale architectures Acceptance of a new language depends on many criteria, including:

0 functionality and target code performance0 mature, industrial-strength compiler and runtime system technology
0 easy integration/migration of legacy codes
0 familiarity of users with conventional features
0 flexibility to deal with new hardware developments
 Many research challenges remain

0 high-level language features for multi-threading
0 architecture- and application-adaptive compilation and runtime systems that
employ intelligent search strategies (ATLAS-like)

0 intelligent tools and middleware that provide efficient support for program
development, performance tuning, fault tolerance, and power management

0 performance-porting of legacy applicationsclass BRD: Distribution {…
function map(i:index(source)):locale{…}; /* global mapping for dense domain */function GetDistributionSegment(loc:locale):domain(1){…}; /* “box” for loc */…
}class CRS: LocalSegment {const loc: locale = this.getLocale();/* declaration of dense and sparse distribution segment for locale loc: */const locD: domain(2);const locDD: sparse domain(locD) = GetDistributionSegment(loc);…
const LocalDomain: domain(1)=1..nnz;/* local data domain *//* persistent data structures in the local segment: */var cx: [LocalDomain] index(locD(2));/* column index vector */var ro: [l1..u1+1] index(xLocalDomain); /* row vector */…
function define_column_vector(): {[z in LocalDomain] cx(z)=nz2x(z)(2)}function define_row_vector(): {…}…
/* mapping global index to index in local data domain: */functionlayout(i: index(D)): index(LocalDomain) return(x2nz(i))
constructor LocalSegment(){define_column_vector(); define_row_vector(); }}Case Study: Introspection Sensors
for Performance Tuning

Implementation Target Architecture:
Cluster of Cell Broadband Engines

Introspection sensors yield information about theexecution of the application:

CBE-1

Cell Broadband Engine CBE-i..
.
Cluster

PowerPC
Processor
Element

PPE
L1

L2

SystemMemory Hardware Monitors
I/O-level Software Monitoring (at message-passing level) LowLow
message
Low-levelmessage-passingI

0 waiting times for blocking send and receive
0 communication transfer times
0 barrier synchronization times

InterConnection

C

CBE-i

Element Interconnect Bus (EIB)

0 …
…

Network

N

Synergistic
Processor
Elements

..
.

SPE-1

0 accumulators: counting standard events (cache misses, loads, FP ops,..)0 timers: analysis of latencies and stalls
0 programmable watch events for special conditions

. . .

SPE-8

-Level Software Monitoring ((atat the level of a high-level
 HighHigh
high
High-Level
high-level
language)

0 timing for redistribution of a globally distributed collection
0 timing for function invocation, loop, or program region
0 timing for computing a communication schedule (“inspector”)
0 evaluation of assertions and invariants
0 …

CBE-n

Fault tolerance must be applied across all levels of the system hierarchy:
SPE Æ PPE Æ CBE Æ Cluster

Case Study: Introspection Actuators
for Performance Tuning

Berkeley’s 7 +6 Dwarfs
1. Dense Linear Algebra (BLAS, ScaLAPACK,MATLAB)
ScaLAPACK,MATLAB)

Introspection actuators provide mechanisms, data, andcontrol paths for implementing feedback to the application,
depending on results of analysis and prediction:
 Instrumentation and Measurement Retargeting Resource Reallocation

2. Sparse Linear Algebra (SpMV
SpMV,, SuperLU)((SpMV,
SuperLU)
3. Spectral Methods (FFT)4. N-Body Methods (Barnes-Hut, Fast Multipole)(BarnesMultipole)
N-Body(Barnes-Hut,
5. Structured Grids (Cactus, Magneto-Hydrodynamics)
MagnetoMagneto-Hydrodynamics)
6. Unstructured Grids (ABAQUS,FIDAP)
7. Monte Carlo Computational Steering
0 changing the implementation of an application section

8. Combination Logic (Encryption; Cyclic Redundancy Codes—
—CRC)
Codes
Codes—CRC)

 changing a function implementation by choosing a more efficient algorithm

9. Graph Traversal (Quicksort
Quicksort))((Quicksort) changing the implementation of a loop

10. Dynamic Programming changing the distribution of key data structures, with the goal of load balancing

 Program Restructuring and Recompilation (offline)

11. Backtrack and Branch and Bound
12. Construction of graphical models (Bayesian networks, Hidden Markov
Markov Models)
13. Finite State Machines

The Traditional Approach will not Scale Traditional approach based on rad-hard processors and fixedradrad-hard

redundancy (e.g.,Triple Modular Redundancy——TMR)
Redundancy
Redundancy—TMR)’09 Launch)
0 Current Generation (Phoenix and Mars Science Lab ––’09Single BAE Rad 750 Processor
256 MB of DRAM and 2 GB Flash Memory (MSL)

200 MIPS peak, 14 Watts available power (14 MIPS/W)

0 ST8 Honeywell Dependable Multiprocessor
COTS system with Rad 750 controller (100 MIPS) and IBM PowerPC 750FX (1300 MIPS) 120 MIPS/Watt Performance






Fault tolerant architecture Rad-hard processors today lag commercial architectures by a factorRad
Rad-hardof about 100 (and growing) By 2015: a single rad-hard processor may deliver about 1 GF—
radGF—
rad-hardorders of magnitude below requirements
 COTS-based multicore systems will be able to provide the required
COTS
COTS-based

capability, but there are serious issues to be addressed…
addressed…

Introspection versus Traditional V&V
 Introspection
0focuses on execution time monitoring, analysis, recovery
0actual work considers transient and hard faults, not design errors
errors Verification & Validation:
0focuses on design errors
0is applied before actual program execution
 Verification has the goal to prove that a program

conforms to its specification for all legal inputs
 Test proves or disproves correctness of the program forspecific (range of) inputs
 Both verification and test are not complete:complete:0problems may be undecidable or intractable
0tests can prove existence of faults, not their total absence

Key Applications for Future Architectures:
Berkeley’s “Dwarfs”
1. Dense Linear Algebra (BLAS, ScaLAPACK,MATLAB)
ScaLAPACK,MATLAB)
2. Sparse Linear Algebra (SpMV
SpMV,, SuperLU)((SpMV,
SuperLU)
3. Spectral Methods (FFT)4. N-Body Methods (Barnes-Hut, Fast Multipole)(BarnesMultipole)
N-Body(Barnes-Hut,
5. Structured Grids (Cactus, Magneto-Hydrodynamics)
MagnetoMagneto-Hydrodynamics)
6. Unstructured Grids (ABAQUS,FIDAP)

Early Alternatives to MPI
 Automatic Vectorization and Parallelization

0 automatic vectorization (for inner loops) and parallelization (for SMPs)
were successful in limited contexts

0 in general, automatic parallelization is essentially intractable Data parallel languages for MPPs and clusters

0 pioneered by compiler projects at Caltech (Cosmic Cube) and U of Bonn
(SUPERB Fortran parallelizer)

0 key features of data parallel languages

7. Monte Carlo
8. Combination Logic (Encryption; Cyclic Redundancy Codes—
—CRC)
Codes
Codes—CRC)
9. Graph Traversal (Quicksort
Quicksort))((Quicksort)
10. Dynamic Programming
11. Backtrack and Branch and Bound
12. Construction of graphical models (Bayesian networks, Hidden Markov
Markov Models)
13. Finite State Machines

Sensors and Actuators
 Sensors and actuators link the introspection framework to

the application and the environment

 Sensors:
Sensors: provide input to the introspection system

Examples for sensor-provided inputs:
sensorsensor-provided
0 state of a variable, data structure, synchronization object
0 value of an assertion
0 state of a temperature sensor or hardware counter Actuators:
Actuators: provide feedback from the introspection system

Examples for actuator-triggered actions:actuator
actuator-triggered
0 modification of program components (methods and data)
0 modification of sensor/actuator sets (including activation and deactivation)
0 local recovery
0 signaling fault to next higher level in a hierarchical system
0 requesting actions from lower levels in a hierarchical system global name space
 single thread of control
 loosely synchronous parallel computation
 automatic generation of communication

0 key language developments
 IVTRAN (1973) – for the SIMD ILLIAV IV – first language to allow control of datalayout MPP languages: Kali, Fortran D, Vienna Fortran, Connection Machine Fortran High Performance Fortran (HPF) result of a standardization effort

Application-Aware Fault Detection
Assertions: Examples
 Assertions based on general program structures
0 values and value ranges for variables, subscript expressions, pointers
0 sequential and parallel control flow patterns
0 locality and communication assertions
0 independence assertions for data-parallel loops
0 real-time constraints
0 safety and liveness properties
 Domain-specific assertions: exploiting knowledge about:

0 target system: hardware and software
0 application domain libraries: pre- and post conditions, argument constraints
 data structure invariants control constraints
 data representation and distribution knowledge (e.g., CRS for distributed

sparse matrices)
 communication patterns and schedules for parallel constructs

Space Flight Avionics and Microprocessors
History and Outlook-induced transient hardware
 SEUs and MBUs are radiationradiation
radiation-induced

Rad-hard
Rad-hard components
components areare always
alwaysat
at least
least 22 generations
generations behindbehindcommercial
commercial State-of-the-Art
State-of-the-Art

errors, which may corrupt software in multiple ways:

Multi-Core Regime
HMC –HeterogeneousMulti-core

100,000,000

HMC (64/1024SP)
HMC (8/256SP)

10,000,000

HMC (4/64SP)
HMC (2/16SP)

1,000,000

HMC (1/9SP)

Computational Rate (MIPS)

Transient Faults

OASIS/
Hyperion

FPGA s)
+gate
(Core

100000

SAR

10000

1000

TS
CO

100

68040/40
1068030/50
68020/33
80486/25

200 200
333
PPC604/132 150
Pentium II/233
PPC603e/133
Pentium/60 68060/75
PPC601/110
80486/50
PPC601/80
Stardust
Mars Pathfinder AIM(RAD6000)(RAD6000)

Clementine HKP(1750)

1

86

Galileo CDS(1802)
87

88

89

90

Mars Global Surveyor(1750A)

Mars Observer EDF(1750A)
91

92

93

94

in
d) S
-har
R ad
ht (
Flig

Cassini(1750A)

Virtx V

RAD750 III

FPGA (Core Only)

Deep Impact(RadLite750)

Deep Space 1(RAD6000)

80386/33

0.1

X 103

Pentium 4/2530
raPentium 4/2000
PPC7470/1250
re E
PPC7455/1000-Co Pentium II/450Pentium III/450 PPC7441/700
gle
Virtx (2/10M)
266 300
Sin Pentium Pro/150
400 PPC7400/450

SIRTF(RAD6000)

e
Cor
gle-

Mars Pathfinder Rover(80C85)
95

96

97

98

Intel

99

00

01

02

03

04

Launch Year

Motorola 680X0

05

06

07

08

09

10

PowerPC

11

12

13

14

15

Missions

16

17

18

19

Multi-Core

20 21 22 23 24 25

FPGA

0 instruction codes and addresses
0 user data structures
0 synchronization objects
0 protected OS data structures
0 synchronization and communication
 Potential effects include:
0 wrong or illegal instruction codes and addresses
0 wrong user data in registers, cache, or DRAM
0 buffer overflows
0 control flow errors
0 unwarranted exceptions
0 hangs and crashes
0 synchronization and communication faults

Source: Contributions from Dan Katz (LSU), Larry Bergman (JPL), and others

Basic Parallel Architecture Paradigms
Control Unit
CPU…

Distributions: Outline of a Formal Framework

CPU
Interconnection Network

Interconnection NetworkP…

A data distribution

Interconnection Network

node

…

LM

…

node

P
MM

LM

 Let I denote the index set of a domain, and L the index domain for a set locales.

MMδ: I Æ L
is a total function that specifies for each element in I an associated locale

LM Let I1,I2, denote index domains. An alignment from I1 to I2, is a total function

LM

Single-Instruction-Multiple-Data (SIMD)
Hardware implementation of a data parallelmodel of computation

Symmetric Multiprocessors (SMP)
Uniform Memory Access (UMA)

Distributed-Memory Multiprocessors(MPPs)
Non-Uniform Memory Access (NUMA)ILLIAC IV, DAP, BSP, CM-1, CM-2, MasPar MP-1…

Alliant FX, Digital VAX 8800, Sequent Balance,…

Cosmic Cube, Suprenum, Transputers,…

Vector Computers

α: I1 Æ I2
that associates an index in I2, with every index of I1. If I2 has a distribution, δ2,then a distribution, δ1, for I1, is obtained as δ1 = δ2 ο αClusters

Cray 1, Cyber 205, …

 Affinity between distributed data and threads can be formalized in a similar way

In modern multicore-based architectures, such building blocks
may be hierarchically combined in many different configurations

Top 500 Architectures

Domains
 Concept influenced by HPF templates, ZPL regions
 Domains are first-class objects
 Domain components

0index set
0distribution
0set of arrays Index sets are general sets of “names”

0Cartesian products of integer intervals (as in Fortran95 etc.)0sparse subsets of Cartesian products
0sets of object instances, e.g., for graph-based data structures
 Iterators based on domains

Example: Possible Extensions for the CELL
Matrix-Vector Multiply
var A: [1..m,1..n] real;var x: [1..n]real;
var y: [1..m]real;(original)
Chapelversiony = sum reduce(dim=2) forall (i,j) in [1..m,1..n] A(i,j)*x(j);param n_spe = 8;
const SPE:[1..n_spe] locale;ExampleMatrix-Vector Multiply on the CELL: V2/* number of synergistic processors (SPEs) *//* declaration of SPE array */var A: [1..m,1..n] real distributed(block,*) on SPE;var x: [1..n]real replicatedon SPE;var y: [1..m]real distributed(block)on SPE;y = sum reduce(dim=2) forall (i,j) in [1..m,1..n] A(i,j)*x(j);Chapel with
(implicit)
heterogeneous
semantics

param n_spe = 8;
const SPE:[1..n_spe] locale;const PPE: locale

var A: [1..m,1..n] real on PPE linked(AA) distributed(block,*) on SPE;var x: [1..n]real on PPE linked(xx) replicatedon SPE;var y: [1..m]real on PPE linked(yy) distributed(block)on SPE;
AA=A; xx=x;/* copy and distribute A, x to SPEs */yy=sum reduce(dim=2) forall (i,j) in [1..m,1..n] on locale(xx(j)) AA(i,j)*xx(j);y=yy;/* copy yy back to PPE */SPE1

PPE Memory
A1
A2
A3
A4
A5
A6
A7
A8

A

y1
y2
y3
y4
y5
y6
y7
y8

y

SPEk local memory (k=4)

SPE3

SPE4
SPE5

xm

Ak: k-th block of rows
yk: k-th block of elements
xk: k-th element

A

SPE6
SPE7

A1
A2
A3
Ak
A5
A6
A7
A8

y1
y2
y3
yk
y5
y6
y7
y8

x1
x2

y

x

xm

A

A1
A2
A3
A4
A5
A6
A7
A8

y1
y2
y3
y4
y5
y6
y7
y8

SPE8

x1
x2

y

SPE4
x

SPE5

AA

SPE6
xm

SPE7

AA1
AA2
AA3
AAk
AA5
AA6
AA7
AA8

yy1
yy2
yy3
yyk
yy5
yy6
yy7
yy8

xx1
xx2

yy

xx

xxm

SPE8

User-Defined Distributions: Global Mapping(2)/* declaration of distribution class MyC1: */class MyC1: Distribution {const ntl:int;
function map(i:index(source)):locale {return Locales(mod(i-1,ntl)+1);}

SPEk local memory (k=4)

SPE2

PPE Memory

SPE3xChapel/HETMC with
explicit transfers

SPE1

SPE2
x1
x2

/* number of synergistic processors (SPEs) *//* declaration of SPE locale array *//* declaration of PPE locale *//* cyclic(1) *//* number of target locales *//* global mapping for MyC1 *//* set of local iterators : */iterator DistSegIterator(loc: index(target)): index(source) {const N: int = getSource().extent;
const k: int = locale_index(loc);
for i in k..N by ntl { yield(i); }}/* distribution segment : */function GetDistributionSegment(loc: index(target)): Domain {const N: int = getSource().extent;
const k: int = locale_index(loc);
return (k..N by ntl);}

An Approach to Application-Oriented
Introspection-Based Fault Tolerance in the HPCS
-dependent) fault model Approach based on a (mission(mission(mission-dependent)
0 classifies faults (fault types, severity)
0 specifies fault probabilities, depending on environment
0 prescribes recovery actions

 Addressing fault detection, analysis, isolation, recovery
 Exploiting knowledge from different sources
0 automatic generation of assertions based on:







static analysis and profilingproperties of target system hardware and software
application domain (libraries, data structures, data distributions)

0 user-provided assertions and invariants

 Leveraging existing technology}/* use of distribution class MyC1 in declarations: */const D1C1: domain(1) distributed(MyC1()) on Locales(1..4)=1..16;var
A1: [D1C1] real;...

X10 and Fortress: Some Key Properties
 X10 --- the IBM HPCS Language

0object-oriented; serial sublanguage based on Java0an array sublanguage supports the distribution of multidimensional arrays via standard methods

0sequential and parallel iterators, either local or global
0asynchronous activities
 Fortress --- the SUN HPCS Language

0object-oriented, with some relationship to Java
0supports Unicode and conventional mathematical notation:e.g.,y = a sin 2 x + cos 2 x log log x

0strong security model
0support for language “growth” via inclusion of libraries
0by default, arrays are distributed and loops are parallel

0 fixed-redundancy for small critical areas in a program
0 Algorithm-Based Fault Tolerance (ABFT): standard matrix methods
0 integration of high-level generator systems such as CMU’s “SPIRAL”

Co-Array Fortran Extension of Fortran to allow SPMD-style programming Introduces a new type of array dimension (co-array) torefer to the cooperating instances (“images”) of an
SPMD program, making processor boundaries explicit:
integer :: a(n.m) [*]this introduces a shared co-array a with n*m integers
local to each processor image
 Non-local variables can be directly referenced based ona corresponding syntax extension:a(1,:) [p]references the first row of co-array a in processor p
 a barrier provides synchronization between images

An Approach to Application-Oriented
Introspection-Based Fault Tolerance in the HPCS

UPC
 Support for a global address space model for SPMD

parallel programs, in which threads share part of theiraddress space The shared space is logically partitioned into fragments,

each of which is associated with a thread Shared arrays are distributed in block-cyclic fashionparallel loop
 Additional features include special constructs forpointers (private/shared), non-blocking barriers, andcollective operations

Example: PGAS vs. HPCS
Setting up a block-distributed array in Titanium vs. Chapel
Titanium: a dialect of Java that supports distributed multi-dimensional arrays,
iterators, subarrays, and synchronization/communication primitives
Chapel Code Fragment

Titanium Code Fragment// determine parameters of local block:
Point<3> startCell = myBlockPos * numCellsPerBlockSide;
Point<3> endCell = startCell + (numCellsPerBlockSide-[1,1,1]);//create local myBlock array:double [3d] myBlock = new double[startCell:endCell];const D: domain(3) distributed (block)= [l1..u1,l2..u2,l3..u3];…
var A: [D] real;…

//build the distributed structure://declare blocks as 1D-array of references (one element per processor)
blocks.exchange(myBlock);

myBlock Addressing fault detection, analysis, isolation, recovery
 Exploiting knowledge from different sources
0 automatic generation of assertions based on: The upc_forall construct supports work sharing for a

blocks

0 classifies faults (fault types, severity)
0 specifies fault probabilities, depending on environment
0 prescribes recovery actions

among threads

P0-dependent) fault model Approach based on a (mission(mission(mission-dependent)

P1

blocksmyBlock

P2

blocksmyBlock

Source: K.Yelick et al.: Parallel Languages and Compilers: Perspective from the Titanium Experience

static analysis and profilingproperties of target system hardware and software
application domain (libraries, data structures, data distributions)

0 user-provided assertions and invariants

 Leveraging existing technology
0 fixed-redundancy for small critical areas in a program
0 Algorithm-Based Fault Tolerance (ABFT): standard matrix methods
0 integration of high-level generator systems such as CMU’s “SPIRAL”

Fortran+MPI Communicationfor 3D 27-point Stencil (NAS MG rprj3)
subroutine
subroutine comm3(u,n1,n2,n3,kk)
comm3(u,n1,n2,n3,kk)
use
use caf_intrinsics
caf_intrinsics

dou(n1,i2,i3) = buff(indx,buff(indx, buff_id )i3=2,n3i3=2,n3-1
doimplicit
implicit nonenone

buff(buff_len,buff(buff_len, buff_id ) = u( i1,
includecafnpb.h''''cafnpb.h'include 'cafnpb.hincludeglobals.h''''globals.h'include 'globals.hdodir = +1dobuff(1:buff_len,buff_id+1)[nbr(axis,dir,k)] =i2=2,n2i2=2,n2-1

buff_id = 3 + dir

indx = indx + 1

buff_len = nm2

do

enddo

indx = indx + 1
u(n1,i2,i3) = buff(indx,buff(indx, buff_id )i1=1,n1do

dir = +1

endif

if( axis .eq.

enddo2 )thenif( dir .eq. -1 )thenbuff(1:buff_len,buff_id+1)[nbr(axis,dir,k)] =>indx = indx + 1

buff_id = 2 + diru(i1,n2,i3) = buff(indx,buff(indx, buff_id )buff_len = 0

dobuff(1:buff_len,buff_id)

i3=2,n3i3=2,n3-1
doif( axis .eq.
doi1=1,n1

endif

enddo
else if( dir .eq. +1 ) then

enddo
enddo

do

i2=1,n2dodo

i1=1,n1i3=2,n3i3=2,n3-1
dobuff_len = buff_len + 1

u(i1,i2,n3) = buff(indx,buff(indx, buff_id )if( axis .eq.
doi1=1,n12 )thenenddo

i3=2,n3i3=2,n3-1
do

endif

i1=1,n1buff_len = buff_len + 1

u(i1,1,i3) = buff(indx,buff(indx, buff_id )enddo

buff(buff_len,buff(buff_len,

enddo

enddo

buff_id )= u( i1,n2i1,n2-1,i3)

buff(1:buff_len,buff_id+1)[nbr(axis,dir,k)] =>buff(1:buff_len,buff_id)

buff_id = 3 + dir

endif

endif

indx = 0

endif
if( axis .eq.

else if( dir .eq. +1 ) thenimplicit
implicit nonenone

if( axis .eq.

dodou(i1,i2,n3) = buff(indx,buff(indx, buff_id )buff(1:buff_len,buff_id)

endif

dir = -1

if( axis .eq.

i1=1,n1

if( axis .eq.
doindx = indx + 1

indx = indx + 1

enddo

1 )thenenddo

i3=2,n3i3=2,n3-1dou(i1,i2,1) = buff(indx,buff(indx, buff_id )returnendif

i2=2,n2i2=2,n2-1
buff_len = buff_len + 1

enddo

buff(buff_len,buff_id ) = u( 2,enddo

i2,i3)

use caf_intrinsics

if( axis .eq.
do

enddo
subroutine take3( axis, dir, u, n1, n2, n3 )

i1=1,n1u(i1,1,i3) = buff(indx,buff(indx, buff_id )i2=1,n2do2 )theni3=2,n3i3=2,n3-1
dobuff_id = 2 + dir
buff_len = 0

do

endif

do
-1
i3=2,n3
do i3=2,n3i3=2,n3-1do-1i2=2,n2
do i2=2,n2i2=2,n2-1

enddo

endif
endif

endif

returnif( axis .eq.

3 )theni2=1,n2do

i1=1,n1
indx = indx + 1
u(i1,i2,1) = buff(indx,buff(indx, buff_id )

enddo
enddo

do

end

implicit none

include 'globals.h''globals.h'2 )thenenddo

i3=2,n3i3=2,n3-1
do

enddo

i1=1,n1

endif

buff_len = buff_len + 1

include 'cafnpb.h''cafnpb.h'

else
else if(
if( dir
dir .eq.
.eq. +1+1 )) thenthenbuff(buff_len,buff(buff_len, buff_id ) = u( i1,subroutine comm1p( axis, u, n1, n2, n3, kk )

2,i3)

enddo

use caf_intrinsics
endif

double precision u( n1, n2, n3 )implicit none

if( axis .eq.

integer buff_id, indxdoinclude 'cafnpb.h''cafnpb.h'
integer i3, i2, i1

include 'globals.h''globals.h'

buff_id = 3 + dir

integer axis, dir, n1, n2, n3

indx = 0

double precision u( n1, n2, n3 )3 )theni2=1,n2
do

enddo
enddo

i1=1,n1buff_len = buff_len + 1
buff(buff_len,buff(buff_len, buff_id ) = u( i1,i2,2)

enddo
enddo

if( axis .eq.

1 )thenif( dir .eq. -1 )then

integer i3, i2, i1, buff_len,buff_id

endif

integer i, kk,
kk, indxdodo

i3=2,n3i3=2,n3-1
do

dir = -1

indx = indx + 1

i=1,nm2
buff(i,4) = buff(i,3)
buff(i,2) = buff(i,1)

i2=2,n2i2=2,n2-1
buff_id = 3 + dir

enddo

buff_len = nm2
dir = -1

returnend

enddo

integer axis, dir, n1, n2, n3

if(
if( axis
axis .eq.
.eq. 22 )then
)then
if(
1 )thenif( dir
dir .eq.
.eq. --1
)then

enddo

enddo
endif

do
else if( dir .eq. +1 ) then

end

endif
endif
endif
endif

indx = indx + 1
u(1,i2,i3) = buff(indx,buff(indx, buff_id )enddo

buff(1:buff_len,buff_id+1)[nbr(axis,dir,k)] =

endif

>>

i2=2,n2i2=2,n2-1

enddo
enddo

buff(1:buff_len,buff_id+1)[nbr(axis,dir,k)]buff(1:buff_len,buff_id+1)[nbr(axis,dir,k)] ==buff(1:buff_len,buff_id)
buff(1:buff_len,buff_id)1 )theni3=2,n3i3=2,n3-1
do

enddo

i1=1,n1
indx = indx + 1

enddo
enddo

>

buff_lenbuff_len == buff_len
buff_len ++ 11
buff(buff_len,, buff_id-1, i2,i3)
buff(buff_len
n1buff(buff_len,
buff_id )) == u(u( n1n1-1,
i2,i3)
enddo
enddo

do

i1=1,n1buff(buff_len,buff(buff_len, buff_id ) = u( i1,i2,n3i1,i2,n3-1)

i2=1,n2dobuff(buff_len,buff(buff_len, buff_id ) = u( i1,i2,n3i1,i2,n3-1)do-1
i3=2,n3
do i3=2,n3i3=2,n3-1do-1i2=2,n2
do i2=2,n2i2=2,n2-1

if( axis .eq.

buff_len = buff_len + 1

i1=1,n1buff_len = buff_len + 1

integer
integer i3,
i3, i2,i2, i1,
i1, buff_len,buff_id
buff_len,buff_id3 )theni2=1,n2do

i2=1,n2do

integer
integer axis,
axis, dir,dir, n1,
n1, n2,n2, n3,n3, k,k, ierrierrdoubledouble precision
precision u(u( n1,
n1, n2,n2, n3n3 ))buff(1:buff_len,buff_id+1)[nbr(axis,dir,k)]buff(1:buff_len,buff_id+1)[nbr(axis,dir,k)] ==buff(1:buff_len,buff_id)
buff(1:buff_len,buff_id)

do3 )thenif( dir .eq. -1 )then

includecafnpb.h''''cafnpb.h'include 'cafnpb.hincludeglobals.h''''globals.h'include 'globals.h

dir = +1

enddo

enddo

enddo

>>

i1=1,n1
indx = indx + 1

endif

indx = indx + 1

buff(buff_len,buff(buff_len, buff_id ) = u( i1,i2,2)

end
end

buff_lenbuff_len == buff_len
buff_len ++ 11
buff(buff_len,buff_id
buff(buff_len,buff_id )) == u(u( 2,2, i2,i3)
i2,i3)
enddo
enddo

3 )theni2=1,n2
do

enddo

call
call zero3(u,n1,n2,n3)
zero3(u,n1,n2,n3)
endif
endif
returnreturnif(
if( axis
axis .eq.
.eq. 11 )then
)then
if(
1 )thenif( dir
dir .eq.
.eq. --1
)then

do

enddo

3 )thenif( dir .eq. -1 )thenif( axis .eq.

buff(buff_len,buff(buff_len, buff_id ) = u( n1n1-1, i2,i3)

enddo
enddo

if( axis .eq.

endif

i2=2,n2i2=2,n2-1
buff_len = buff_len + 1

u(i1,n2,i3) = buff(indx,buff(indx, buff_id )

endif

enddo

1 )theni3=2,n3i3=2,n3-1
doindx = indx + 1

call
call sync_all()sync_all()call
call sync_all()sync_all()buff_idbuff_id == 22 ++ dir
dirbuff_lenbuff_len == 00

i1=1,n1

enddo

endif
endif

subroutine
subroutine give3(
give3( axis,
axis, dir,dir, u,u, n1,
n1, n2,n2, n3,n3, kk ))use
use caf_intrinsics
caf_intrinsics2 )theni3=2,n3i3=2,n3-1
dobuff_id )= u( i1,n2i1,n2-1,i3)

enddo

call
call sync_all()sync_all()
call
1, u,
call take3(
take3( axis,
axis, --1,
u, n1,
n1, n2,n2, n3n3 ))if( axis .eq.

endifbuff_len = buff_len + 1
buff(buff_len,buff(buff_len,

endif

enddo

i3=2,n3i3=2,n3-1
docall
call give3(
give3( axis,
axis, +1,
+1, u,
u, n1,
n1, n2,n2, n3,n3, kk
kk ))
call
1, u,
call give3(
give3( axis,
axis, --1,
u, n1,
n1, n2,n2, n3,n3, kk
kk ))

enddo

i=1,nm2
buff(i,buff_id) = 0.0D0

enddo

if(
if( nprocs
nprocs .ne.
.ne. 1)
1) thenthencall
call sync_all()sync_all()i2=2,n2i2=2,n2-1

enddo

u(1,i2,i3) = buff(indx,buff(indx, buff_id )do1 )theni3=2,n3i3=2,n3-1
doi3=2,n3i3=2,n3-1

buff(1:buff_len,buff_id)

else if( dir .eq. +1 ) then

enddo
enddo
else
else
do
do axisaxis == 1,
1, 33

indx = 0if( axis .eq.

else if( dir .eq. +1 ) then

enddo

if(
) )thendead(kkif( .not.
.not. dead(kk)dead(kk))thendodo axisaxis == 1,
1, 33

call
call take3(
take3( axis,
axis, +1,
+1, u,
u, n1,
n1, n2,n2, n3
n3 ))
else
else
call
call comm1p(
comm1p( axis,
axis, u,u, n1,
n1, n2,n2, n3,n3, kkkk ))buff_id = 3 + dir

2,i3)

enddo

>

i=1,nm2
buff(i,buff_id) = 0.0D0

enddo

enddo

do
integer
integer n1,
n1, n2,n2, n3,n3, kkkk
doubledouble precision
precision u(n1,n2,n3)u(n1,n2,n3)
integer
integer axis
axis

do

enddo

i1=1,n1buff_len = buff_len + 1

科

学 技 術 計 算 分 科 会

選 出2008 年度 科学技術計算分科会 より

アプリケーションとアーキテクチャ

工学院大学
小柳 義夫

サイエンティフィック・システム研究会

200８年 8 月 27 日

SS研 HPCフォーラム2008「ペタスケール・アプリを支える基盤技術」

アプリケーションとアーキテクチャ
工学院大学
小柳義夫

[アブストラクト]
コンピュータのアーキテクチャとアプリケーションがどのような相互作用を行って来たかについて総括する。
最初の電子式自動計算機 ENIAC は、弾頭計算という特定のアプリのために製作されたが、多くの目的のために利用された。1964 年の IBM 360 シリーズの発表以来「汎用コンピュータ」がコンピュータ設計の基本原理とされ、どのようなアプリケーションにも対応できるという幻想を振りまいた。1976 年の Cray-1 を嚆矢とするスーパーコンピュータの出現は、科学技術計算を目的としたコンピュータの意義を確認させた。その後、マイクロプロセッサの高速化と低価格化により「汎用チップを用いた専用計算機」が可能になったが、マルチコア化のためにデータ転送がボトルネックになりつつある。与えられたコンピュータシステムでいかに速く計算するかという努力は重要であるが、新しいコンピュータはアプリケーション（複数）を想定して設計するべきである。
何にでも速いコンピュータを作るという目標では設計上の判断ができず中途半端なものになってしまうであろう。

[キーワード]アーキテクチャ、アルゴリズム、専用計算機、並列計算機、科学技術計算［概要］
１． はじめにこの講演を依頼されて思い出したのは、現在神戸で建設が進んでいる次世代スーパーコンピュータの概念
設計評価作業部会（2007 年）での議論である。文部科学省の提案によると、これまでの日本の主導的スーパーコンピュータであった、NWT, cp-pacs, 地球シミュレータが専用計算機であったのに対し、これは汎用スーパーコンピュータであり、どんな問題にも対応できるということが強調された。私やＩ先生はこれに強く反対
し、NWT も cp-pacs も地球シミュレータも、特定のアプリを目的として設計されたが、実は汎用的計算機であり、特定のアプリケーションを強く意識して設計したがゆえに、（逆説的に）多くのアプリケーションに役立つ汎用コンピュータになったことを指摘した。私が言うまでもなく、設計とは、有限の資源のもとにおける多くの技
術要素のトレードオフであるから、最初から汎用を設計目標にすると、「あれもほしい、これもほしい」ということになって、中途半端なものになってしまう危険がある。
考えてみると、普通の計算機が汎用で特定目的の計算機が例外、という概念は 1960 年代以降のものであって、最初からそうだったわけではない。チューリングやフォン・ノイマンの計算機モデルは先駆的であったが、初期の計算機は特定目的の計算機であり、しばらくは事務用と科学技術計算用の計算機は別に作られていた。これを統一したＩＢＭの System360 は画期的であり、その後の各社のコンピュータのモデルとなったが、これは直前の状況へのいわばアンチテーゼとして出されたものであった。1980 年代以降には、ベクトルスーパ
ーコンピュータや並列計算機が出現し、状況は変わったにもかかわらず、「コンピュータは汎用であるべきである」という神話が生きているのは残念なことである。
以下では、「アプリを意識した計算機の方が役に立つ」という仮説を証明しようとしたのであるが、実際はな
かなか難しい。

２． 機械式計算器におけるアーキとアプリ
現存する最古の歯車式計算器は、パスカルの作った「パスカリーヌ」(1642)と言われているが、そのアプリ
は父親の税金計算であった。ライプニッツは桁ずらしにより乗除算の出来る計算器を作ったが、そのアプリ
は天文計算であった。バベッジの階差機関(1822)の目的は、天文学や暦だと言われている。機械式計算器
でも、アプリが想定されていることは（考えてみれば当たり前であるが）注目に値する。日本では、1902 年に矢頭良一が発明し製作した「自答算盤」は政府関係に多く売れたそうである。アプリは？

３． 計算機第１世代
真空管式計算機の ENIAC (1946)は配線によってプログラムするもので、その初期の目的は弾道計算であった。その後プログラム内蔵式に改良され、多くの目的に使われた。

４． 計算機第 2 世代（トランジスタ）＝HPC(-3)世代(STRETCH)
1955 年頃からトランジスタ式の計算機が登場したが、10 進演算が主の事務用計算機と、浮動小数演算が主の科学技術用計算機とに分かれていた。私は学生時代に OKITAC5090 を使ったことがあるが、10 進 12桁の計算機であった。浮動小数もあったが、指数部 2 桁、仮数部 10 桁の 10 進表現であった。本来、事務用だったのでしょうか？
このころ IBM は STRETCH の開発をはじめ、1961 年にロスアラモス研究所に納入した。7 台しか製作されなかったが、既存の計算機より 100 倍高速で、スーパーコンピュータの走りと言えよう。アプリは 3 次元の流体
力学計算などであったと思われる。

５． 計算機第 3 世代（メインフレーム）＝HPC(-2)世代(CDC6600, IBM2938)
1964 年ＩＢＭはＩＣを使った 360 シリーズを発表し、科学技術用・事務用を問わず 360 度の全方位に向けたコンピュータを標榜した。汎用コンピュータの歴史がここに始まり、われわれもその呪縛のもとにある。
同じ年、アメリカ原子力委員会は各社に「革新的な計算機」を製作するよう促し、ここから後年の CDC
Star-100、TI ASC、 ILLIAC-IV などが生まれたと言われている。やはりこの年、CDC 社は CDC6600 を発表
し、科学技術用計算機として多数納入された。アメリカ空軍がイリノイ大学と ILLIAC-IV の契約を結んだのもこの年である。
IBM は最初のベクトル計算機 IBM 2938 Array Processor を 1968 年に出荷した。IBM360 の I/O チャンネルに接続するもので、アプリは石油探査であった。「汎用コンピュータ」を提唱した IBM 自身が、このようなアプリ
を明確に意識した計算機を出していることは注目に値する。

６． 計算機 3．5 世代＝HPC(-1)世代(ASC, Star-100, ILLIAC IV, BSP)
IBM は 1970 年に LSI を使った 370 シリーズを発表し、各社もこれに追随したが、このころベクトル計算機が登場しつつあった。TI 社の ASC と、CDC 社の Star-100 である。また、並列計算機としては、ILLIAC IV (1976)
や BSP(1980 年開発中止)などがある。

７． HPC 神代時代(Cray-1, 75APU, IAP, DAP, Cosmic Cube, MPP)
Cray-1 が登場したのは 1976 年であるが、そのアプリは科学技術全般であろう。同時にこのころ、様々な並列計算機も登場した。ICL DAP、CM*、LAU、Cosmic Cube、MPP などである。日本では FACOM 230-75 APU
や、日立等の IAP が登場している。これらはまだ本格的なベクトル計算機とは言えないので、HPC 神代時代と名付けることにする。注目すべきことは、このころ日本でも多くの並列計算機の試作が行われていることである。しかしいずれも実験機であって商品にはなったものはない。アプリは、シミュレーション、画像処理、信号処理などであった。
このうち、寿命の長かったのは、Cosmic Cube→ipsc→Paragon→ASCI Red、pacs→pax→cp-pacs→
pacs-cs や Grape シリーズなどで、いずれもアプリが明確なものである、と私には見える。

８． HPC 第 1 世代(FPS, Cyber205, XMP, S810, V200, SX-2, PAX, iPSC/1, NCUBE/1, CM-1, T)
その後、1980 年代になると、ベクトルコンピュータの全盛期を迎えるが、注目すべきことは、アメリカなどでは、同時に並列コンピュータの商品が山ほど登場していることである。日本では試作機は多いが、あくまで特定目的に限定されていた。

９． 専用シミュレータ（Ising, QCD, MD）
このころから ASIC が利用可能になり、Ising 模型シミュレータ、格子ゲージ模型シミュレータ、古典粒子系シミュレータ、モンテカルロ専用機などいろいろ登場している。

１０．1990 年代の HPC
1993 年に Cray が T3D を、IBM が SP-1 を発表して、超並列分野に参入した。その後、1980 年代に登場していた多くのベンチャー超並列企業の倒産・吸収が続いた。その原因は様々で、この事情の裏にアプリがどう関係しているかは不明である。

１１．ベンチマークと性能評価
コンピュータの性能を評価するためにベンチマーク・プログラムが用いられる。だから、ベンチマーク性能の高いコンピュータがよいコンピュータだと思いがちである。ベンチマークとしては、Mix、MIPS、LFK、Linpack、SPEC などいろいろあるが、いずれも万能ではない。
そもそもベンチマークはアプリケーションとアーキテクチャとの界面である。アプリケーションから見ればそのアプリケーションの実行速度を推定するためのデータである。従って、ユーザ側としては、対象となるアプ
リケーションの実行速度をよく推定できるようなベンチマークが望ましい。逆にアーキテクチャから見れば、アーキテクチャを設計するときの目標である。すべてのアプリケーションに対して最適なシステムを設計することは困難であるので、代表的な負荷を取って設計を考える。この両面（とくに後者）を忘れてはならない。
結局、ベンチマークもアプリケーションと独立なものでなく、性能の絶対評価ではない。

１２．おわりに
歴史的にどこまで証明できたかは問題であるが、コンピュータの歴史においてアプリケーションが重要な論点であったことはご理解いただけたと思う。今後、アーキテクチャが複雑化し、マルチコア・メニーコアが普及すると、かなりクセのあるコンピュータが利用されるようになることが予想される。どんな計算も高速に実行できることは望み薄で、アプリケーションとの相性が問題になる。さらには、モデル化、アルゴリズムまでも再考する必要が生じるであろう。

─企画委員会からの要望─

アーキテクチャと応用アプリの歴史総括• 何を狙ってそのアーキができて、アプリはどの様に対応してきたか。
• 実装技術、その時の課題、ブレークスルー、
応用アプリの特性、高速化技術など。
• 大変難しい要望である。どこまで答えられた
かは自信がないが。

アプリケーションとアーキテクチャ
小柳義夫
工学院大学情報学部

2008/8/27

1

2008/8/27

2

何を言いたかったのか

概要• 日本の次世代スーパーコンピュータ計画「汎用」スーパーである、と強調された。曰く「NWT
やcp-pacsや地球シミュレータは専用スーパーであった」• 「『汎用』では設計が出来ない。」 I先生や私が強調。
• 「普通の計算機が『汎用』で、特定目的の計算機は例外」という考え方は時代的な概念。
• 「汎用」コンピュータはアンチテーゼとして登場• 仮説「アプリを意識した計算機の方が役立つ。」を歴史的に証明したかった。
• それは成功したか？？？？
2008/8/27

3歯車式計算器におけるアーキとアプリ
計算機第１世代（真空管）
計算機第2世代（トランジスタ）＝HPC第(-3)世代(STRETCH)
計算機第3世代（メインフレーム）＝HPC第(-2)世代(CDC6600, IBM2938)
計算機3．5世代＝HPC第(-1)世代(ASC, Star-100,
ILLIAC IV, BSP)
HPC神代時代(Cray-1, 75APU, IAP, DAP, Cosmic
Cube, MPP)
HPC第1世代(FPS, Cyber205, XMP, S810, V200, SX2, PAX, iPSC/1, NCUBE/1, CM-1, T-series)
専用シミュレータ（Ising, QCD, MD）
性能とベンチマーク
コンピュータの設計戦略

•
•
•
•
•
•
•
•
•
•
2008/8/27

4「汎用」コンピュータの神話

歯車式計算器• 歴史的に、コンピュータは目的（アプリ）を想定して設計され、製造されてきた。
• アンチテーゼとしての「汎用」コンピュータの概念• 最初の計算機械：ウィルヘルム・シッカート(1592－
1635)チュービンゲン大学教授（現存せず）• 1642年、パスカル「パスカリーヌ」（加減算）– アプリ：税金計算– チューリング・マシン（計算可能性）– フォン・ノイマン・アーキテクチャ（プログラム内蔵）– IBM System 360 --- 「メインフレーム」の登場– Attack of killer micros• まず、歴史的な「計算器」から
2008/8/27

5• 167４年、ライプニッツ（桁ずらしにより乗除算も）– アプリ：天文計算• 1902年、矢頭良一が発明し製作した「自答算盤」– 販売先：陸軍省、内務省、農事試験場• 1923年、大本寅治郎「虎印計算器」、後のタイガー
計算器
2008/8/27

6

ライプニッツ(1646年-1716年)

チャールズ・バベッジ（Charles Babbage ）• 1792～1871• 1822年に階差機関（difference engine）構想• 微分や積分の記号を発明• 形式言語に当たるものを初めて考案• 2進法を研究• どんな推論も代数計算のように単純で機械的な作業に置き換えることができる（人工知能？）• チューリングの先駆？
2008/8/27– 完成せず。1991年、バベッジの本来の設計に基づいて階差機関が組み立てられ、完全に機能。
– アプリ：天文学、暦• 1830年代、さらに汎用的な解析機関を構想。いかなる数学の関数も計算できる機能を持つ。
– パンチカード（Jacquard Loomで使用）で制御– 世界で初めての「プログラム可能」な計算機– これも実現せず。
7

2008/8/27

8

パンチカードシステム（Punch Card System ）コンピュータの前夜• ハーマン・ホレリスが発明。1890年の米国国勢調査のデータ処理で初めて使用された。
• 余談：Fortranで3HABCのＨは彼に由来。
• 80桁のパンチカード（私も使った）• 1896年、ホレリスはTabulating Machine社を設立（IBM社の源流の一つ）• 1936年、アラン・チューリングが万能計算機械（チューリング
マシン）の論文を発表。計算可能性、unsolvable problem• 1938年、ドイツのコンラート・ツーゼが、機械式の計算機V1（後にZ1と改名）を作成。アプリは暗号か？
• 1939年、ツーゼが演算部がリレー、記憶部が機械式の計算機Z2を作成。
• 1940年、ツーゼが全リレー式のZ3を作成。Z3はプログラム
可能な最初の計算機である。その後Z4も作成。
• 1942年、ジョン・アタナソフとクリフォード・ベリーが電子素子を使って演算処理をする世界初の機械 ABCを作成。
アプリは連立1次方程式。実際には稼働せず。
• 1944年、Harvard Mark I出荷（電気機械式計算器）
アプリは？

2008/8/27

2008/8/27

9

計算機第１世代（真空管）

計算機第２世代（トランジスタ）• 1946年：最初の真空管式コンピュータENIACが完成– プログラムは配線による。後にプログラム内蔵へ。
– 10進法の回路。
– アプリ：弾道計算が初期の目的。その後いろいろ。

• まもなく商用化：UNIVAC-1（1951年）、IBM701（1953年）

2008/8/27

10

11• 1955年頃からトランジスタが利用される• NEAC2201(1958)、 IBM7090(1958)、 CDC1604 (1958)、
UNIVAC1100(1960)、 TOSBAC2100(1959)、
HITAC301(1959)、OKITAC5090(1960)、
FACOM222(1961)など• 証券会社、商事会社、保険会社などに納入• 科学技術計算用：浮動小数演算が主。
事務計算用：十進演算が主。
両者が区別されていた。
2008/8/27

12ＨＰＣ（-3）世代• 1956：IBM starts 7030 project (known as
STRETCH) 100 times faster than IBM 704
using transistors (2進法)– Atomic Energy Commission at Los Alamos.
– 1961: first STRETCH computer to LANL (only 7
built )– アプリ：三次元の流体力学計算など– much of the technology re-surfaces in the later
IBM 7090 and 7094.
2008/8/27

13

2008/8/27

IBM STRETCH (1960)14

http://www03.ibm.com/ibm/history/history/images/ibm_stretch.jpg

第３世代（メインフレーム）

第３世代• 1964年：IBMがICを使った360シリーズを発表• 科学技術用・事務用を問わず360度の全方位に向けたコンピュータを標榜• 汎用コンピュータの歴史が始まる。
• ファミリー概念• 第３世代メインフレームの代表例• UNIVAC1108(1965) 、GE635/645(1965)、
FACOM230(1965)、HITAC 8400(1966)、
NEAC2200(1964)、CDC6400(1966)、
B8500(1967)、NCRのCentury100(1968)、
TOSBAC5600(1971)• 汎用大型計算機の時代– 同じアーキテクチャをもった小型機から大型機まで用意– 新機種開発にあたっても過去の機種との互換性を重視
2008/8/27

15

2008/8/27

ＨＰＣ (-２)世代

ＨＰＣ (-２)世代• 1964：Atomic Energy Commission urges
manufacturers to look at "radical” machine
structures. This leads to CDC Star-100, TI
ASC, and ILLIAC-IV.
• 1964: Control Data Corporation produces
CDC 6600（科学技術用計算機）• 1964: Air Force signs ILLIAC-IV contract
with University of Illinois. (Burroughs andTI)
2008/8/27

16

17• ベクトル計算機の登場– 演算パイプラインによってベクトル演算を高速に実行– 1965年にIBM社のSenzigらによって最初に提案。
– パイプライン方式と並列方式を比較• 石油探査のためIBM 2938 Array Processor(1968)– IBM360のI/Oチャンネルに接続する付加プロセッ
サとして開発2008/8/27

18

Vector registers

2008/8/27

出典：IBM Systems Journal 8, #2, pp. 118-135 (1969)19

2008/8/27

20

1.59 MFlops

2008/8/27

21

2008/8/27

HPC（-１）世代

メインフレーム第3.5世代• ベクトル（パイプライン）計算機• IBMは1970年にLSIを使ったIBM370シリーズを発表し、各社もこれに追随• 日本では、通産省がIBM対抗機の開発を促進するため行政指導を行って企業連合– ASC (1972, Texas Instruments社）• 30 MFlops、pipeline 1本、7機製作• 高速性。拡張性のある設計。アプリは？
• Cray-1の登場により一夜のうちに販売停止– 富士通・日立はMシリーズ（IBM互換）– 三菱・沖はCOSMOシリーズ(IBM非互換）– 日本電気・東芝はACOSシリーズ(IBM非互換）– Star-100(1973, CDC社）STring ARray Computer• CDC7600の流れ。100 Mflopsを目標。
• Memory-to-memory architecture。立ち上がり遅い。
• スカラー速度を犠牲にした。結局性能は出ず。
• 50 MFlops、4機製作。2台はLLNLへ。アプリは？

• いずれも1974年に発表した。
2008/8/27

22

23

2008/8/27

24

ILLIAC IV

HPC（-１）世代
並列計算機– ILLIAC IV (1973, Burroughs社、1976完全稼働）• 並列度64、50 MFlops、1機製作• Illinois大学からNASAへ。アプリ：流体力学など• 並列処理の人材を世界中に（日本にも）送り出した。

– BSP (1974設計開始、1980開発中止、Burroughs社）• 並列とパイプライン技術を併用• 並列度16、理論性能50 MFLops• 17個のmemory unit+ crossbar switch• アプリは？
2008/8/27

25

2008/8/27

26

http://ed-thelen.org/comp-hist/vs-illiac-iv.html

Cray-1

HPC神代時代• Seymour Cray
– CDC (Control Data Corporation) 社において、
CDC6600 (1964, 1 MFlops) およびCDC7600(1969、5 MFLops) を設計– 1972年、CDC8600 計画が社内で拒否される– 同社を退社し、CRI社 (Cray Research Inc.) を設立• CRI社は1976年160 MFlopsの性能をもつ
Cray-1を出荷し、ロスアラモス研究所に納入
2008/8/27

27

2008/8/27

HPC神代時代（Cray-1)

解
HPC神代時代（富士通）
誤
で

• 実装

• FACOM 230-75 APU (1977）– 4ゲートのICを高密度に実装するという画期的技術
でこのような性能を実現したことは驚異– メモリを高速化。クロック12.5ns。スカラーも速かった。

• 自動ベクトル化コンパイラ– 画期的であった。でも、– 性能は十分でなく、使いこなすには技能が必要であった• アプリは？ 科学技術計算全般
2008/8/27

http://www.tfh-wildau.de/l202/team2/cray1.jpg28

29

これ

ま

– 22 MFlops、2機製作、航空技術研究所に納入– ベクトルレジスタを持つベクトル計算機– AP-FORTRANという拡張言語方式– 間接参照や条件文のベクトル演算をサポートしていた点は注目される（CrayではXMP後半から）– アプリ：流体計算など？

• 商業的には成功とは言えないが、日本の最初のベクトル計算機であった。
2008/8/27

30

HPC神代時代（日立）

HPC神代時代（日立）• IAP (Integrated Array Processor)というメイ
ンフレームに対する付加プロセッサ• 間接参照、総和、内積、1次漸化式の自動並列化– アプリ：汎用• HITAC M-180 IAP (1978),• M-200H IAP (1979, 48MFlops)（筆者が利用）• M-280H IAP (1982, 67 MFlops)– 仮想空間上のデータに対してベクトル演算– 性能向上はほどほど（数倍程度）– 高度な自動ベクトル化コンパイラを装備– TSSでも使える
2008/8/27

31– Cray-1がまだ完全にはサポートしていなかった機能を有していたことが特徴である。

• M-280H IAP は、世界で初めて条件付きdoloopを自動ベクトル化できた。
• キャッシュに頼ったベクトル演算には限界がありメインフレーム自体の高速化とともに姿を消した。
2008/8/27

32

HPC神話時代の日本の並列計算機

HPC神代時代（NEC他）• 1979：GMMCS(群馬大)、PACS-9(京大)• 1980：PPA（北大）、KDSS-1(豊田工大）など1983年までに30近く（高橋義造氏による）。
いずれも実験機。商品はなし。
• マイクロプロセッサの登場により可能• 日本電気– ACOS-1000 IAP (1982, 28 Mflops)• 三菱電機– MELCOM COSMO IAP– 詳細は不明– MC6800, Z80, LSI11, i8086など– 専用LSIを用いたものはHAL, QA2など少数• アプリ：シミュレーション、画像、信号処理など
2008/8/27

33

2008/8/27

HPC神代時代の外国の並列計算機• 1980：DAP（専用）, CM*(LSI-11), LAU(DF,i8085+Am2901)• 1981：Manchester DFM（専用）• 1983：Cosmic Cube（i8086+i8087), MPP
（専用）

34

ICL DAP• 1972: paper by Dr Stewart F Reddaway• 1974: prototype 製造開始
• 1979: Queen Mary Collegeに納入– 64x64 1 bit PE (SIMD)– アプリ：？？

括弧内はCPUチップ。

2008/8/27

35

2008/8/27

36

神代のアプリとアーキ

Goodyear MPP• ベクトル計算機は汎用科学技術用、並列計算機は専用的計算機という位置づけ。
• マイクロプロセッサの登場により、アプリに特化したコンピュータが可能になった。
• 超並列機は、欧米では商品開発、日本では実験機のみ。
• アプリが明確なものは寿命が長い(?)• 198?: Goodyear's array processor, a4x256 1-bit processing element (PE)computer.
• 1983: Goodyear’s MPP, a 128x128 2-D
array of 1-bit PEs. 納入：NASA Goddard
Space Flight Center– アプリ：初期は衛星画像解析、その後、開口合成
レーダー、海流、宇宙線、NN,レイトレーシングなど。
2008/8/27

37

38(continued)
1984
1984
1984
1985
1985
1985
1985
1985
1985
1986

1981
1981
1982
1982

FPS-164 (64 bits)
CDC Cyber 205
400MF
Cray XMP-2 Steve Chen 630MF
Cosmic Cube in Caltech, Alliant FX/8 delivered,
HEP installed
1983 HITAC S-810/20
630MF
1983 FACOM VP-200
570MF
1983 Encore, Sequent and TMC founded,
ETA span off from CDC

39

専用シミュレータ（この頃から）• この頃からASICが利用可能になった。
• Ising模型シミュレータ：米UCSB(1982)、オランダ
Delft工科大学(1983)、東大(1988)など• 格子ゲージ模型シミュレータ：イタリアのAPE(1986)、筑波大学のqcdpax(1989)とcp-pacs(1996)、アメリカのGF11(1991), ACP-MAPS, QCDSP(1995)、
QCDOC(2002)など• 古典粒子系シミュレータ：Delft工科大学のDMDP(1988)、東大のGRAPE(1989～)、理研のMD-GRAPE (1999)、富士ゼロックスのMD Engine(1996)など• モンテカルロ専用機：Monte-4(1993原研、vector型）
2008/8/27

2008/8/27

1st Generation (1H of 1980’s)

1st Generation (1H of 1980’s)

2008/8/27– Cosmic Cube→ipsc→Paragon→ASCI Red– pacs, pax→cp-pacs, pacs-cs– Grape series

41

Multiflow founded
Cray XMP-4 1260MF
PAX-64J completed (Tsukuba)
NEC SX-2 1300MF
FPS-264
Convex C1
Cray-2
1952MFIntel iPSC/1, T414, NCUBE/1, Stellar, Ardent…
FACOM VP-400
1140MF
CM-1 shipped, FPS T-series (max 1TF!!)

2008/8/27

40

専用から汎用へ• 汎用指向の格子ゲージシミュレータ– cp-pacs→SR2201/SR8000– QCDSP, QCDOC→BG/L– APE, GF11（いずれもSIMD)は汎用化せず• その他の例（果たして専用だったのか？）– NWT(数値風洞）→VPP500/300/700– 地球シミュレータ→SX-6/7/8

2008/8/27

42

性能とは

1990年代のHPC• IBMとCrayの超並列参入• コンピュータの性能評価：ベンチマーク– 1993 Cray T3D (Alpha chip)– 1993 IBM SP-1 (Power chip)アーキテクチャ• ベンチャー超並列企業の倒産・吸収– 1995– 1996

Convex, KSR, nCUBE, Pyramid
MasPar, TMC

ベンチマーク• アプリケーションとアーキテクチャの関係？？
– 倒産の原因は？アプリが不明確だったから？？
アプリケーション
2008/8/27

43

2008/8/27

ベンチマークコンピュータの設計戦略• Mix, MIPSからLFK, Linpack, Euroben, SPEC,
HPCC, ….
• アプリケーションから見れば：自分のアプリケー
ションの計算速度を推定するためのデータ• アーキテクチャから見れば：アーキテクチャをチューニングする時の目標• ベンチマークの選択：コンピュータの利用目的の選択。
2008/8/27

4445• 設計：有限な資源の中での選択「あれもほしい」「これもほしい」ではだめ• 「汎用」の陥穽（おとしあな）：汎用といっても、暗黙のアプリのクラスを前提にしている。無前提の「汎用」はない。
• 多くの成功したHPCは、特定のアプリを前提に設計されている。結果的に多くのアプリでもよい性能を出すことが多い。
2008/8/27

46