Data Parallel Programming in R David Padua Department - - PowerPoint PPT Presentation

Data ¡Parallel ¡Programming ¡ in ¡R ¡

David ¡Padua ¡ ¡ Department ¡of ¡Computer ¡Science ¡ University ¡of ¡Illinois ¡at ¡Urbana-­‑Champaign ¡

Outline ¡

• Parallelism ¡
• Data ¡Parallel ¡Programming ¡and ¡abstrac?ons ¡
• Hierarchically ¡Tiled ¡Arrays ¡
• Future ¡plans ¡

¡2 ¡

• I. ¡Parallelism ¡
• Parallelism ¡is ¡crucial ¡for ¡ ¡

– Con?nued ¡gains ¡in ¡performance ¡ – Maximum ¡performance ¡at ¡any ¡given ¡?me ¡ – Also ¡the ¡most ¡natural ¡way ¡to ¡program ¡for ¡reac?ve ¡ compu?ng ¡ ¡(but ¡not ¡the ¡topic ¡of ¡this ¡presenta?on) ¡

• Main ¡problem ¡with ¡parallelism ¡is ¡produc?vity. ¡
• Need ¡the ¡right ¡languages, ¡libraries ¡and ¡tools ¡

¡3 ¡

• II. ¡Data ¡parallel ¡programming ¡
• In ¡its ¡simplest ¡form ¡is ¡just ¡the ¡execu?on ¡of ¡the ¡

same ¡opera?on ¡on ¡each ¡element ¡of ¡an ¡aggregate ¡ (array, ¡set, ¡database ¡rela?on). ¡

• Sequen?al ¡execu?on ¡across ¡these ¡opera?ons ¡
• Crucial ¡issue ¡is ¡what ¡should ¡these ¡opera?ons ¡

should ¡look ¡like ¡(research ¡problem) ¡

• There ¡are ¡numerous ¡proposals ¡ ¡

– Array ¡opera?ons ¡(Iversion ¡ca. ¡1960) ¡ – MapReduce ¡(Google, ¡ca. ¡2000) ¡ – Galois ¡(Pingali, ¡ca. ¡2000) ¡

¡4 ¡

• II. ¡Data ¡parallel ¡programming ¡

Array ¡Constructs ¡

• Popular ¡among ¡scien?sts ¡and ¡engineers. ¡

– Fortran ¡90 ¡and ¡successors ¡ – MATLAB ¡ – R ¡

• Parallelism ¡not ¡the ¡reason ¡for ¡this ¡nota?on. ¡

¡5 ¡

• II. ¡Data ¡parallel ¡programming ¡

Array ¡Constructs ¡

• Convenient ¡nota?on ¡ ¡

– Compact ¡ – Higher ¡level ¡of ¡abstrac?on ¡

do i=1,n do i=1,n do j=1,n do j=1,n C(i,j)= A(i,j)+B(i,j) S = S + A(i,j) end do end do end do end do C = A + B S += sum(A) ¡

¡6 ¡

• II. ¡Data ¡parallel ¡programming ¡

Array ¡constructs ¡

• Used ¡in ¡the ¡past ¡for ¡parallelism: ¡Illiac ¡IV, ¡

Connec?on ¡machine ¡

• Today: ¡Intel’s ¡Cilk ¡(mainly ¡for ¡

microprocessor ¡vector ¡extensions) ¡

¡7 ¡

• II. ¡Data ¡parallel ¡programming ¡

Benefits ¡-­‑ ¡Programmability ¡

• Data ¡parallel ¡programming ¡is ¡scalable ¡

– Scales ¡with ¡increasing ¡number ¡of ¡processors ¡by ¡increasing ¡the ¡size ¡of ¡ the ¡data ¡

• Data ¡parallel ¡programs ¡using ¡powerful ¡operator ¡resemble ¡

conven?onal, ¡serial ¡programs. ¡

– Parallelism ¡is ¡encapsulated. ¡ – Parallelism ¡is ¡structured ¡

¡ ¡

• Portable ¡

– Can ¡run ¡on ¡any ¡class ¡of ¡machine ¡for ¡which ¡the ¡appropriate ¡operators ¡ are ¡implemented ¡

• Shared/Distributed ¡Memory, ¡Vector ¡Intrinsics, ¡GPUs ¡
• Interoperates ¡with ¡R ¡! ¡

¡8 ¡

Opera&ons ¡implemented ¡as ¡ ¡ messages ¡if ¡distributed ¡memory ¡ Opera&ons ¡implemented ¡as ¡ ¡ parallel ¡loops ¡in ¡shared ¡memory ¡ Opera&ons ¡implemented ¡with ¡ ¡ vector ¡intrinsics ¡for ¡SIMD ¡

• II. ¡Data ¡parallel ¡programming ¡

Completeness ¡

• Can ¡all ¡problems ¡be ¡solved ¡in ¡the ¡most ¡

efficient ¡manner ¡with ¡data ¡parallel ¡ programming ¡? ¡

• Most ¡? ¡All ¡?? ¡

– Other ¡(lower ¡level)forms ¡must ¡be ¡there ¡in ¡the ¡ same ¡way ¡that ¡we ¡s?ll ¡use ¡assembly ¡language ¡ some?mes. ¡

¡9 ¡

• II. ¡Data ¡parallel ¡programming ¡

Completeness ¡

• Numerical ¡compu?ng ¡ ¡

– William ¡H. ¡Press, ¡Saul ¡A. ¡Teukolsky, ¡William ¡T. ¡Veeerling, ¡and ¡Brian ¡P. ¡

• Flannery. ¡1996. ¡Numerical ¡Recipes ¡in ¡Fortran ¡90 ¡(2nd ¡Ed.): ¡The ¡Art ¡of ¡Parallel ¡

Scien?fic ¡Compu?ng. ¡Cambridge ¡University ¡Press, ¡New ¡York, ¡NY, ¡USA. ¡

• Graph ¡algorithms ¡

– Aydın ¡Buluç ¡and ¡John ¡R ¡Gilbert. ¡2011. ¡The ¡Combinatorial ¡BLAS: ¡design, ¡ implementa&on, ¡and ¡applica&ons. ¡Int. ¡J. ¡High ¡Perform. ¡Comput. ¡Appl. ¡25, ¡4 ¡ (November ¡2011), ¡496-­‑509. ¡ – Keshav ¡Pingali, ¡Donald ¡Nguyen, ¡Milind ¡Kulkarni, ¡Mar?n ¡Burtscher, ¡M. ¡Amber ¡ Hassaan, ¡Rashid ¡Kaleem, ¡Tsung-­‑Hsien ¡Lee, ¡Andrew ¡Lenharth, ¡Roman ¡ Manevich, ¡Mario ¡Méndez-­‑Lojo, ¡Dimitrios ¡Prountzos, ¡and ¡Xin ¡Sui. ¡2011. ¡The ¡ tao ¡of ¡parallelism ¡in ¡algorithms. ¡In ¡Proceedings ¡of ¡the ¡32nd ¡ACM ¡SIGPLAN ¡ conference ¡on ¡Programming ¡language ¡design ¡and ¡implementa?on ¡(PLDI ¡'11). ¡ ACM, ¡New ¡York, ¡NY, ¡USA, ¡12-­‑25. ¡ ¡

• Database ¡algorithms ¡

– Anand ¡Rajaraman ¡and ¡Jeff ¡Ullman. ¡Mining ¡of ¡Massive ¡Datasets ¡.Cambridge ¡ University ¡Press, ¡2011. ¡

• … ¡

¡10 ¡

• II. ¡Data ¡parallel ¡programming ¡

Transla?ng ¡to ¡SPMD ¡

• SPMD ¡is ¡the ¡nota?on ¡of ¡

choice ¡for ¡distribute ¡ memory ¡machines ¡(and ¡ GPUs). ¡

• Easy ¡to ¡convert ¡from ¡

array ¡nota?on ¡to ¡SPMD ¡

• form. ¡
• This ¡is ¡an ¡op?miza?on. ¡

¡11 ¡ real ¡a, ¡b, ¡x(1000) ¡ a=sin ¡(b) ¡ x(:)=x(:)+a ¡ real ¡a, ¡b, ¡x(1000/p) ¡ ¡ /* ¡a ¡and ¡b ¡are ¡replicated ¡*/ ¡ a=sin(b) ¡ x(:)=x(:)+a ¡

• III. ¡Hierarchically ¡Tiled ¡Arrays: ¡ ¡

Our ¡data ¡parallel ¡nota?on ¡for ¡array ¡ computa?ons ¡

¡12 ¡

• Recognizes ¡the ¡importance ¡of ¡blocking/?ling ¡

for ¡locality ¡and ¡parallel ¡programming. ¡ ¡

• Makes ¡?les ¡first ¡class ¡objects. ¡ ¡

– Referenced ¡explicitly. ¡ ¡ – Manipulated ¡using ¡array ¡opera?ons ¡such ¡as ¡ reduc?ons, ¡gather, ¡etc.. ¡

• G. ¡Bikshandi, ¡J. ¡Guo, ¡D. ¡Hoeflinger, ¡G. ¡Almasi, ¡B. ¡Fraguela, ¡M. ¡Garzarán, ¡D. ¡Padua, ¡ ¡

and ¡C. ¡von ¡Praun. ¡Programming ¡for ¡ ¡Parallelism ¡and ¡Locality ¡with ¡Hierarchically ¡Tiled. ¡ ¡ PPoPP, ¡March ¡2006. ¡ ¡

Hierarchically ¡Tiled ¡Arrays ¡

¡13 ¡

Locality Locality Distributed

Hierarchically ¡Tiled ¡Arrays ¡

¡14 ¡

Vector Addressing

a(1:2,1:2) a(1,1:4) In general, a(vn)

¡15 ¡

HTA Addressing

¡16 ¡

HTA Addressing

h{1,1:2} (hta) h{2,1} (array)

hierarchical ¡17 ¡

HTA Addressing

h(3,4) ↔ scalar

flattened

h{1,1:2} (hta) h{2,1} (array)

hierarchical ¡18 ¡

HTA Addressing

h{1:2,2}(1:2,2) ↔ hta

hybrid

h(3,4) ↔ scalar

flattened

h{1,1:2} (hta) h{2,1} (array)

hierarchical ¡19 ¡

h(3,4) ↔ scalar

flattened

HTA Addressing

h{ i + j == 3} ↔ hta

logical indexing

h{1:2,2}(1:2,2) ↔ hta

hybrid

h{1,1:2} (hta) h{2,1} (array)

hierarchical ¡20 ¡

repmat(h, [1, 3]) circshift( h, [0, -1] ) transpose(h)

Higher ¡level ¡opera?ons ¡

¡21 ¡

Higher ¡Level ¡Opera?ons ¡

• Many ¡operators ¡part ¡of ¡the ¡library ¡

– reduce, ¡circular ¡shiv, ¡replicate, ¡transpose, ¡etc ¡

• A ¡map ¡opera?on ¡(hmap) ¡

– Applies ¡user ¡defined ¡operators ¡to ¡each ¡?le ¡of ¡the ¡ HTA ¡

• And ¡corresponding ¡?les ¡if ¡mul?ple ¡HTAs ¡are ¡passed ¡as ¡

input ¡

– Applica?on ¡of ¡operator ¡occurs ¡in ¡parallel ¡across ¡ ?les ¡

¡

¡22 ¡

User ¡Defined ¡Opera?ons ¡

HTA ¡X(3,3)[10] ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡// ¡3x3 ¡?les ¡of ¡10 ¡elements ¡ HTA ¡Y(3,3)[10] ¡ … ¡ hmap( ¡F(), ¡X, ¡Y ¡) ¡ ¡ F(HTA ¡x, ¡HTA ¡y) ¡{ ¡ ¡ ¡ ¡ ¡ ¡y ¡[i] ¡= ¡x[i] ¡* ¡x[i] ¡-­‑ ¡3 ¡ } ¡

¡23 ¡

F(…) F(…) F(…) F(…) F(…) F(…) F(…) F(…) F(…)

X Y

Cannon's Matrix Multiplication

A00 A01 A10 A20 A11 A21 A22 A12 A02 B01 B10 B20 B11 B21 B22 B12 B02 B00 A00 B00 A01 B11 A02 B22 A12 B21 A11 B10 A10 B02 A22 B20 A20 B01 A21 B12 A00 B00 A01 B11 A02 B22 A12 B21 A11 B10 A10 B02 A22 B20 A20 B01 A21 B12

initial skew shift-multiply-add

¡24 ¡

Cannon's Matrix Multiplication

%Main loop for i = 1:n c = c + a * b; a = circshift( a, [0, -1] ); b = circshift( b, [-1, 0] ); end

¡25 ¡

FT

¡26 ¡

Advantages ¡of ¡?ling ¡as ¡a ¡first ¡class ¡object ¡for ¡

• p?miza?on ¡
• HTAs ¡have ¡been ¡implemented ¡as ¡C++ ¡and ¡

MATLAB ¡libraries. ¡

– For ¡shared ¡and ¡distributed ¡memory ¡machines. ¡

• Dense ¡and ¡sparse ¡versions ¡
• Implemented ¡several ¡benchmark ¡suites. ¡
• Performance ¡is ¡compe??ve ¡with ¡OpenMP, ¡ ¡MPI, ¡

and ¡TBB ¡counterparts ¡

• Furthermore, ¡the ¡HTA ¡nota?on ¡produces ¡code ¡

more ¡readable ¡than ¡other ¡nota?ons. ¡It ¡ significantly ¡reduces ¡number ¡of ¡lines ¡of ¡code. ¡ ¡

¡27 ¡

EP CG MG FT LU Lines of code

Lines ¡of ¡Code. ¡HTA ¡vs. ¡MPI ¡

¡ ¡Advantages ¡of ¡?ling ¡as ¡a ¡first ¡class ¡object ¡ ¡ ¡

¡28 ¡

Performance ¡Results ¡

¡29 ¡

1 2 4 8 16 32 64 128 1 2 3 4 5 6 7 8

Number of processors Running time in seconds (log2)

G5 NAS G5 HTA gen G5 HTA opt Itanium NAS Itanium HTA opt 1 2 4 8 16 32 64 128 3 4 5 6 7 8 9 10

Number of processors Running time in seconds (log2)

G5 NAS G5 HTA gen G5 HTA opt Itanium NAS Itanium HTA opt 1 2 4 8 16 32 64 128 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5

Number of processors Running time in seconds (log2)

G5 NAS G5 HTA gen G5 HTA opt Itanium NAS Itanium HTA opt 1 2 4 8 16 32 64 128 3 4 5 6 7 8 9 10 11 12 13

Number of processors Running time in seconds (log2)

G5 NAS G5 HTA gen G5 HTA opt Itanium NAS Itanium HTA opt

MG ¡ IS ¡ FT ¡ CG ¡

Plans ¡

¡

• We ¡are ¡considering ¡extending ¡ROpt ¡(our ¡

extended ¡R ¡interpreter) ¡with ¡HTA ¡opera?ons ¡ and ¡enable ¡execu?on ¡on ¡distributed ¡memory ¡

• machines. ¡

Work ¡with ¡Peng ¡Wu ¡(IBM ¡Research) ¡and ¡Haichuan ¡Wang ¡(Illinois) ¡

¡30 ¡

ROpt ¡

• Extends ¡R ¡byte ¡code ¡

interpreter ¡using ¡ specializa?on. ¡

– New ¡Op ¡codes ¡

• For ¡specific ¡data ¡types ¡
• For ¡frequently ¡occurring ¡

code ¡sequences ¡

– Simpler ¡data ¡ representa?on ¡ – Automa?c ¡op?miza?on ¡

¡31 ¡

Preliminary ¡results ¡with ¡ROpt ¡

• On ¡a ¡set ¡of ¡kernels ¡containing ¡mainly ¡scalar ¡opera?ons ¡ ¡

– ROpt ¡delivers ¡average ¡speedup ¡of ¡3.56 ¡over ¡the ¡bytecode ¡ R ¡interpreter ¡ ¡ – ROpt ¡is ¡12 ¡?mes ¡slower ¡than ¡C ¡ – The ¡bytecode ¡interpreter ¡speedup ¡over ¡the ¡original ¡R ¡ interpreter ¡is ¡2.6. ¡

• On ¡the ¡shootout ¡benchmarks ¡

– 2.07 ¡speedup ¡over ¡the ¡bytecode ¡interpreter ¡but ¡100 ¡?mes ¡ slower ¡than ¡C. ¡ – The ¡bytecode ¡interpreter ¡is ¡2.5 ¡?mes ¡faster ¡than ¡the ¡

• riginal ¡interpreter. ¡

¡32 ¡