Document Sample

Automatic Performance Tuning of Sparse-Matrix-Vector-Multiplication (SpMV) and Iterative Sparse Solvers James Demmel www.cs.berkeley.edu/~demmel/cs267_Spr09 Outline • Motivation for Automatic Performance Tuning • Results for sparse matrix kernels – Sparse Matrix Vector Multiplication (SpMV) – Sequential and Multicore results • OSKI = Optimized Sparse Kernel Interface • Tuning Entire Sparse Solvers Berkeley Benchmarking and OPtimization (BeBOP) • Prof. Katherine Yelick • Current members – Kaushik Datta, Mark Hoemmen, Marghoob Mohiyuddin, Shoaib Kamil, Rajesh Nishtala, Vasily Volkov, Sam Williams, … • Previous members – Hormozd Gahvari, Eun-Jim Im, Ankit Jain, Rich Vuduc, many undergrads, … • Many results here from current, previous students • bebop.cs.berkeley.edu Automatic Performance Tuning • Goal: Let machine do hard work of writing fast code • What does tuning of dense BLAS, FFTs, signal processing, have in common? – Can do the tuning off-line: once per architecture, algorithm – Can take as much time as necessary (hours, a week…) – At run-time, algorithm choice may depend only on few parameters (matrix dimensions, size of FFT, etc.) • Can’t always do tuning off-line – Algorithm and implementation may strongly depend on data only known at run-time – Ex: Sparse matrix nonzero pattern determines both best data structure and implementation of Sparse-matrix-vector- multiplication (SpMV) – Part of search for best algorithm just be done (very quickly!) at run-time Source: Accelerator Cavity Design Problem (Ko via Husbands) Linear Programming Matrix … A Sparse Matrix You Encounter Every Day SpMV with Compressed Sparse Row (CSR) Storage Matrix-vector multiply kernel: y(i) y(i) + A(i,j)*x(j) for each row i Only 2 flops per for k=ptr[i] to ptr[i+1] do 2 mem_refs: y[i] = y[i] + val[k]*x[ind[k]] Limited by getting data from memory Example: The Difficulty of Tuning • n = 21200 • nnz = 1.5 M • kernel: SpMV • Source: NASA structural analysis problem Example: The Difficulty of Tuning • n = 21200 • nnz = 1.5 M • kernel: SpMV • Source: NASA structural analysis problem • 8x8 dense substructure: exploit this to limit #mem_refs Speedups on Itanium 2: The Need for Search Mflop/s Best: 4x2 Reference Mflop/s Register Profile: Itanium 2 1190 Mflop/s 190 Mflop/s Power3 - 17% 252 Mflop/s Power4 - 16% 820 Mflop/s Register Profiles: IBM and Intel IA-64 122 Mflop/s 459 Mflop/s Itanium 1 - 8% 247 Mflop/s Itanium 2 - 33% 1.2 Gflop/s 107 Mflop/s 190 Mflop/s Ultra 2i - 11% 72 Mflop/s Ultra 3 - 5% 90 Mflop/s Register Profiles: Sun and Intel x86 35 Mflop/s 50 Mflop/s Pentium III - 21% 108 Mflop/s Pentium III-M - 15% 122 Mflop/s 42 Mflop/s 58 Mflop/s Another example of tuning challenges • More complicated non-zero structure in general • N = 16614 • NNZ = 1.1M Zoom in to top corner • More complicated non-zero structure in general • N = 16614 • NNZ = 1.1M 3x3 blocks look natural, but… • More complicated non-zero structure in general • Example: 3x3 blocking – Logical grid of 3x3 cells • But would lead to lots of ―fill-in‖ Extra Work Can Improve Efficiency! • More complicated non-zero structure in general • Example: 3x3 blocking – Logical grid of 3x3 cells – Fill-in explicit zeros – Unroll 3x3 block multiplies – ―Fill ratio‖ = 1.5 • On Pentium III: 1.5x speedup! – Actual mflop rate 1.52 = 2.25 higher Automatic Register Block Size Selection • Selecting the r x c block size – Off-line benchmark • Precompute Mflops(r,c) using dense A for each r x c • Once per machine/architecture – Run-time “search” • Sample A to estimate Fill(r,c) for each r x c – Run-time heuristic model • Choose r, c to minimize time ~ Fill(r,c) / Mflops(r,c) Accurate and Efficient Adaptive Fill Estimation • Idea: Sample matrix – Fraction of matrix to sample: s [0,1] – Control cost = O(s * nnz ) by controlling s • Search at run-time: the constant matters! • Control s automatically by computing statistical confidence intervals, by monitoring variance • Cost of tuning – Lower bound: convert matrix in 5 to 40 unblocked SpMVs – Heuristic: 1 to 11 SpMVs • Tuning only useful when we do many SpMVs – Common case, eg in sparse solvers Accuracy of the Tuning Heuristics (1/4) See p. 375 of Vuduc’s thesis for matrices NOTE: “Fair” flops used (ops on explicit zeros not counted as “work”) Accuracy of the Tuning Heuristics (2/4) DGEMV Accuracy of the Tuning Heuristics (2/4) Upper Bounds on Performance for blocked SpMV • P = (flops) / (time) – Flops = 2 * nnz(A) • Upper bound on speed: Two main assumptions – 1. Count memory ops only (streaming) – 2. Count only compulsory, capacity misses: ignore conflicts • Account for line sizes • Account for matrix size and nnz • Charge minimum access ―latency‖ ai at Li cache & amem – e.g., Saavedra-Barrera and PMaC MAPS benchmarks Time a i Hits i a mem Hits mem i 1 a1 Loads (a i 1 a i ) Misses i (a mem a ) Misses i 1 Example: Bounds on Itanium 2 Example: Bounds on Itanium 2 Example: Bounds on Itanium 2 Summary of Other Performance Optimizations • Optimizations for SpMV – Register blocking (RB): up to 4x over CSR – Variable block splitting: 2.1x over CSR, 1.8x over RB – Diagonals: 2x over CSR – Reordering to create dense structure + splitting: 2x over CSR – Symmetry: 2.8x over CSR, 2.6x over RB – Cache blocking: 2.8x over CSR – Multiple vectors (SpMM): 7x over CSR – And combinations… • Sparse triangular solve – Hybrid sparse/dense data structure: 1.8x over CSR • Higher-level kernels – A·AT·x, AT·A·x: 4x over CSR, 1.8x over RB – A2·x: 2x over CSR, 1.5x over RB – [A·x, A2·x, A3·x, .. , Ak·x] …. more to say later Potential Impact on Applications: Omega3P • Application: accelerator cavity design [Ko] • Relevant optimization techniques – Symmetric storage – Register blocking – Reordering, to create more dense blocks • Reverse Cuthill-McKee ordering to reduce bandwidth – Do Breadth-First-Search, number nodes in reverse order visited • Traveling Salesman Problem-based ordering to create blocks – Nodes = columns of A – Weights(u, v) = no. of nz u, v have in common – Tour = ordering of columns – Choose maximum weight tour – See [Pinar & Heath ’97] • 2.1x speedup on IBM Power 4 Source: Accelerator Cavity Design Problem (Ko via Husbands) Post-RCM Reordering 100x100 Submatrix Along Diagonal “Microscopic” Effect of RCM Reordering Before: Green + Red After: Green + Blue “Microscopic” Effect of Combined RCM+TSP Reordering Before: Green + Red After: Green + Blue (Omega3P) Optimized Sparse Kernel Interface - OSKI • Provides sparse kernels automatically tuned for user‘s matrix & machine – BLAS-style functionality: SpMV, Ax & ATy, TrSV – Hides complexity of run-time tuning – Includes new, faster locality-aware kernels: ATAx, Akx • Faster than standard implementations – Up to 4x faster matvec, 1.8x trisolve, 4x ATA*x • For ―advanced‖ users & solver library writers – Available as stand-alone library (OSKI 1.0.1h, 6/07) – Available as PETSc extension (OSKI-PETSc .1d, 3/06) – Bebop.cs.berkeley.edu/oski How the OSKI Tunes (Overview) Library Install-Time (offline) Application Run-Time 1. Build for Workload Target 2. Benchmark Matrix from program Arch. monitoring History Generated 1. Evaluate Benchmark Heuristic code Models data models variants 2. Select To user: Data Struct. Matrix handle & Code for kernel calls Extensibility: Advanced users may write & dynamically add “Code variants” and “Heuristic models” to system. How to Call OSKI: Basic Usage • May gradually migrate existing apps – Step 1: ―Wrap‖ existing data structures – Step 2: Make BLAS-like kernel calls int* ptr = …, *ind = …; double* val = …; /* Matrix, in CSR format */ double* x = …, *y = …; /* Let x and y be two dense vectors */ /* Compute y = ·y + a·A·x, 500 times */ for( i = 0; i < 500; i++ ) my_matmult( ptr, ind, val, a, x, , y ); How to Call OSKI: Basic Usage • May gradually migrate existing apps – Step 1: ―Wrap‖ existing data structures – Step 2: Make BLAS-like kernel calls int* ptr = …, *ind = …; double* val = …; /* Matrix, in CSR format */ double* x = …, *y = …; /* Let x and y be two dense vectors */ /* Step 1: Create OSKI wrappers around this data */ oski_matrix_t A_tunable = oski_CreateMatCSR(ptr, ind, val, num_rows, num_cols, SHARE_INPUTMAT, …); oski_vecview_t x_view = oski_CreateVecView(x, num_cols, UNIT_STRIDE); oski_vecview_t y_view = oski_CreateVecView(y, num_rows, UNIT_STRIDE); /* Compute y = ·y + a·A·x, 500 times */ for( i = 0; i < 500; i++ ) my_matmult( ptr, ind, val, a, x, , y ); How to Call OSKI: Basic Usage • May gradually migrate existing apps – Step 1: ―Wrap‖ existing data structures – Step 2: Make BLAS-like kernel calls int* ptr = …, *ind = …; double* val = …; /* Matrix, in CSR format */ double* x = …, *y = …; /* Let x and y be two dense vectors */ /* Step 1: Create OSKI wrappers around this data */ oski_matrix_t A_tunable = oski_CreateMatCSR(ptr, ind, val, num_rows, num_cols, SHARE_INPUTMAT, …); oski_vecview_t x_view = oski_CreateVecView(x, num_cols, UNIT_STRIDE); oski_vecview_t y_view = oski_CreateVecView(y, num_rows, UNIT_STRIDE); /* Compute y = ·y + a·A·x, 500 times */ for( i = 0; i < 500; i++ ) oski_MatMult(A_tunable, OP_NORMAL, a, x_view, , y_view);/* Step 2 */ How to Call OSKI: Tune with Explicit Hints • User calls ―tune‖ routine – May provide explicit tuning hints (OPTIONAL) oski_matrix_t A_tunable = oski_CreateMatCSR( … ); /* … */ /* Tell OSKI we will call SpMV 500 times (workload hint) */ oski_SetHintMatMult(A_tunable, OP_NORMAL, a, x_view, , y_view, 500); /* Tell OSKI we think the matrix has 8x8 blocks (structural hint) */ oski_SetHint(A_tunable, HINT_SINGLE_BLOCKSIZE, 8, 8); oski_TuneMat(A_tunable); /* Ask OSKI to tune */ for( i = 0; i < 500; i++ ) oski_MatMult(A_tunable, OP_NORMAL, a, x_view, , y_view); How the User Calls OSKI: Implicit Tuning • Ask library to infer workload – Library profiles all kernel calls – May periodically re-tune oski_matrix_t A_tunable = oski_CreateMatCSR( … ); /* … */ for( i = 0; i < 500; i++ ) { oski_MatMult(A_tunable, OP_NORMAL, a, x_view, , y_view); oski_TuneMat(A_tunable); /* Ask OSKI to tune */ } Multicore SMPs Used for Tuning SpMV Intel Xeon E5345 (Clovertown) AMD Opteron 2356 (Barcelona) Sun T2+ T5140 (Victoria Falls) IBM QS20 Cell Blade Multicore SMPs with Conventional cache-based memory hierarchy Intel Xeon E5345 (Clovertown) AMD Opteron 2356 (Barcelona) Sun T2+ T5140 (Victoria Falls) IBM QS20 Cell Blade Multicore SMPs with local store-based memory hierarchy Intel Xeon E5345 (Clovertown) AMD Opteron 2356 (Barcelona) Sun T2+ T5140 (Victoria Falls) IBM QS20 Cell Blade Multicore SMPs with CMT = Chip-MultiThreading Intel Xeon E5345 (Clovertown) AMD Opteron 2356 (Barcelona) Sun T2+ T5140 (Victoria Falls) IBM QS20 Cell Blade Multicore SMPs: Number of threads Intel Xeon E5345 (Clovertown) AMD Opteron 2356 (Barcelona) 8 threads 8 threads Sun T2+ T5140 (Victoria Falls) IBM QS20 Cell Blade 128 threads 16* threads *SPEs only Multicore SMPs: peak double precision flops Intel Xeon E5345 (Clovertown) AMD Opteron 2356 (Barcelona) 75 GFlop/s 74 Gflop/s Sun T2+ T5140 (Victoria Falls) IBM QS20 Cell Blade 19 GFlop/s 29* GFlop/s *SPEs only Multicore SMPs: total DRAM bandwidth Intel Xeon E5345 (Clovertown) AMD Opteron 2356 (Barcelona) 21 GB/s (read) 10 GB/s (write) 21 GB/s Sun T2+ T5140 (Victoria Falls) IBM QS20 Cell Blade 42 GB/s (read) 21 GB/s (write) 51 GB/s *SPEs only Multicore SMPs with Non-Uniform Memory Access - NUMA Intel Xeon E5345 (Clovertown) AMD Opteron 2356 (Barcelona) Sun T2+ T5140 (Victoria Falls) IBM QS20 Cell Blade *SPEs only Set of 14 test matrices • All bigger than the caches of our SMPs 2K x 2K Dense matrix stored in sparse format Dense Well Structured (sorted by nonzeros/row) FEM / FEM / Wind FEM / FEM / Protein QCD Economics Epidemiology Spheres Cantilever Tunnel Harbor Ship Poorly Structured hodgepodge FEM / Circuit webbase Accelerator Extreme Aspect Ratio (linear programming) LP SpMV Performance: Naive parallelization • Out-of-the box SpMV performance on a suite of 14 matrices • Scalability isn’t great: Compare to # threads 8 8 128 16 Naïve Pthreads Naïve SpMV Performance: NUMA and Software Prefetching NUMA-aware allocation is essential on NUMA SMPs. Explicit software prefetching can boost bandwidth and change cache replacement policies used exhaustive search SpMV Performance: ―Matrix Compression‖ Compression includes register blocking other formats smaller indices Use heuristic rather than search SpMV Performance: cache and TLB blocking +Cache/LS/TLB Blocking +Matrix Compression +SW Prefetching +NUMA/Affinity Naïve Pthreads Naïve SpMV Performance: Architecture specific optimizations +Cache/LS/TLB Blocking +Matrix Compression +SW Prefetching +NUMA/Affinity Naïve Pthreads Naïve SpMV Performance: max speedup • Fully auto-tuned SpMV performance across the suite of matrices • Included SPE/local store 2.7x 4.0x • optimized version Why do some optimizations work better on some architectures? +Cache/LS/TLB Blocking 2.9x 35x +Matrix Compression +SW Prefetching +NUMA/Affinity Naïve Pthreads Naïve Avoiding Communication in Sparse Linear Algebra • k-steps of typical iterative solver for Ax=b or Ax=λx – Does k SpMVs with starting vector (eg with b, if solving Ax=b) – Finds ―best‖ solution among all linear combinations of these k+1 vectors – Many such ―Krylov Subspace Methods‖ • Conjugate Gradients, GMRES, Lanczos, Arnoldi, … • Goal: minimize communication in Krylov Subspace Methods – Assume matrix ―well-partitioned,‖ with modest surface-to-volume ratio – Parallel implementation • Conventional: O(k log p) messages, because k calls to SpMV • New: O(log p) messages - optimal – Serial implementation • Conventional: O(k) moves of data from slow to fast memory • New: O(1) moves of data – optimal • Lots of speed up possible (modeled and measured) – Price: some redundant computation Locally Dependent Entries for [x,Ax], A tridiagonal, 2 processors Proc 1 Proc 2 A8x A7x A6x A5x A4x A3x A2x Ax x Can be computed without communication Locally Dependent Entries for [x,Ax,A2x], A tridiagonal, 2 processors Proc 1 Proc 2 A8x A7x A6x A5x A4x A3x A2x Ax x Can be computed without communication Locally Dependent Entries for [x,Ax,…,A3x], A tridiagonal, 2 processors Proc 1 Proc 2 A8x A7x A6x A5x A4x A3x A2x Ax x Can be computed without communication Locally Dependent Entries for [x,Ax,…,A4x], A tridiagonal, 2 processors Proc 1 Proc 2 A8x A7x A6x A5x A4x A3x A2x Ax x Can be computed without communication Locally Dependent Entries for [x,Ax,…,A8x], A tridiagonal, 2 processors Proc 1 Proc 2 A8x A7x A6x A5x A4x A3x A2x Ax x Can be computed without communication k=8 fold reuse of A Remotely Dependent Entries for [x,Ax,…,A8x], A tridiagonal, 2 processors Proc 1 Proc 2 A8x A7x A6x A5x A4x A3x A2x Ax x One message to get data needed to compute remotely dependent entries, not k=8 Minimizes number of messages = latency cost Price: redundant work “surface/volume ratio” Remotely Dependent Entries for [x,Ax,A2x,A3x], A irregular, multiple processors Sequential [x,Ax,…,A4x], with memory hierarchy v One read of matrix from slow memory, not k=4 Minimizes words moved = bandwidth cost No redundant work Performance results on 8-Core Clovertown Optimizing Communication Complexity of Sparse Solvers • Example: GMRES for Ax=b on ―2D Mesh‖ – x lives on n-by-n mesh – Partitioned on p½ -by- p½ grid of p processors – A has ―5 point stencil‖ (Laplacian) • (Ax)(i,j) = linear_combination(x(i,j), x(i,j±1), x(i±1,j)) – Ex: 18-by-18 mesh on 3-by-3 grid of 9 processors Minimizing Communication of GMRES • What is the cost = (#flops, #words, #mess) of k steps of standard GMRES? GMRES, ver.1: for i=1 to k w = A * v(i-1) n/p½ MGS(w, v(0),…,v(i-1)) update v(i), H endfor solve LSQ problem with H n/p½ • Cost(A * v) = k * (9n2 /p, 4n / p½ , 4 ) • Cost(MGS) = k2/2 * ( 4n2 /p , log p , log p ) • Total cost ~ Cost( A * v ) + Cost (MGS) • Can we reduce the latency? Minimizing Communication of GMRES • Cost(GMRES, ver.1) = Cost(A*v) + Cost(MGS) = ( 9kn2 /p, 4kn / p½ , 4k ) + ( 2k2n2 /p , k2 log p / 2 , k2 log p / 2 ) • How much latency cost from A*v can you avoid? Almost all GMRES, ver. 2: W = [ v, Av, A2v, … , Akv ] [Q,R] = MGS(W) Build H from R, solve LSQ problem k=3 • Cost(W) = ( ~ same, ~ same , 8 ) • Latency cost independent of k – optimal • Cost (MGS) unchanged • Can we reduce the latency more? Minimizing Communication of GMRES • Cost(GMRES, ver. 2) = Cost(W) + Cost(MGS) = ( 9kn2 /p, 4kn / p½ , 8 ) + ( 2k2n2 /p , k2 log p / 2 , k2 log p / 2 ) • How much latency cost from MGS can you avoid? Almost all GMRES, ver. 3: W = [ v, Av, A2v, … , Akv ] *Q,R+ = TSQR(W) … “Tall Skinny QR” Build H from R, solve LSQ problem W1 R1 R12 W = W2 R2 R1234 W3 R3 R34 W4 R4 • Cost(TSQR) = ( ~ same, ~ same , log p ) • Latency cost independent of s - optimal Minimizing Communication of GMRES • Cost(GMRES, ver. 2) = Cost(W) + Cost(MGS) = ( 9kn2 /p, 4kn / p½ , 8 ) + ( 2k2n2 /p , k2 log p / 2 , k2 log p / 2 ) • How much latency cost from MGS can you avoid? Almost all GMRES, ver. 3: W = [ v, Av, A2v, … , Akv ] *Q,R+ = TSQR(W) … “Tall Skinny QR” Build H from R, solve LSQ problem W1 R1 R12 W = W2 R2 R1234 W3 R3 R34 W4 R4 • Cost(TSQR) = ( ~ same, ~ same , log p ) • Oops Minimizing Communication of GMRES • Cost(GMRES, ver. 2) = Cost(W) + Cost(MGS) = ( 9kn2 /p, 4kn / p½ , 8 ) + ( 2k2n2 /p , k2 log p / 2 , k2 log p / 2 ) • How much latency cost from MGS can you avoid? Almost all GMRES, ver. 3: W = [ v, Av, A2v, … , Akv ] *Q,R+ = TSQR(W) … “Tall Skinny QR” Build H from R, solve LSQ problem W1 R1 R12 W = W2 R2 R1234 W3 R3 R34 W4 R4 • Cost(TSQR) = ( ~ same, ~ same , log p ) • Oops – W from power method, precision lost! Minimizing Communication of GMRES • Cost(GMRES, ver. 3) = Cost(W) + Cost(TSQR) = ( 9kn2 /p, 4kn / p½ , 8 ) + ( 2k2n2 /p , k2 log p / 2 , log p ) • Latency cost independent of k, just log p – optimal • Oops – W from power method, so precision lost – What to do? • Use a different polynomial basis • Not Monomial basis W = [v, Av, A2v, …+, instead … • Newton Basis WN = [v, (A – θ1 I)v , (A – θ2 I)(A – θ1 I)v, …+ or • Chebyshev Basis WC = [v, T1(v), T2(v), …+ Speed ups on 8-core Clovertown Conclusions • Fast code must minimize communication – Especially for sparse matrix computations because communication dominates • Generating fast code for a single SpMV – Design space of possible algorithms must be searched at run-time, when sparse matrix available – Design space should be searched automatically • Biggest speedups from minimizing communication in an entire sparse solver – Many more opportunities to minimize communication in multiple SpMVs than in one – Requires transforming entire algorithm – Lots of open problems EXTRA SLIDES Quick-and-dirty Parallelism: OSKI-PETSc • Extend PETSc‘s distributed memory SpMV (MATMPIAIJ) • PETSc p0 – Each process stores diag (all-local) and off-diag submatrices p1 • OSKI-PETSc: – Add OSKI wrappers p2 – Each submatrix tuned independently p3 OSKI-PETSc Proof-of-Concept Results • Matrix 1: Accelerator cavity design (R. Lee @ SLAC) – N ~ 1 M, ~40 M non-zeros – 2x2 dense block substructure – Symmetric • Matrix 2: Linear programming (Italian Railways) – Short-and-fat: 4k x 1M, ~11M non-zeros – Highly unstructured – Big speedup from cache-blocking: no native PETSc format • Evaluation machine: Xeon cluster – Peak: 4.8 Gflop/s per node Accelerator Cavity Matrix OSKI-PETSc Performance: Accel. Cavity Linear Programming Matrix … OSKI-PETSc Performance: LP Matrix Performance Results • Measured Multicore (Clovertown) speedups up to 6.4x • Measured/Modeled sequential OOC speedup up to 3x • Modeled parallel Petascale speedup up to 6.9x • Modeled parallel Grid speedup up to 22x • Sequential speedup due to bandwidth, works for many problem sizes • Parallel speedup due to latency, works for smaller problems on many processors Speedups on Intel Clovertown (8 core) Extensions • Other Krylov methods – Arnoldi, CG, Lanczos, … • Preconditioning – Solve MAx=Mb where preconditioning matrix M chosen to make system ―easier‖ • M approximates A-1 somehow, but cheaply, to accelerate convergence – Cheap as long as contributions from ―distant‖ parts of the system can be compressed • Sparsity • Low rank Design Space for [x,Ax,…,Akx] (1/3) • Mathematical Operation – Keep all vectors • Krylov Subspace Methods – Improving conditioning of basis • W = [x, p1(A)x, p2(A)x,…,pk(A)x] • pi(A) = degree i polynomial chosen to reduce cond(W) – Preconditioning (Ay=b MAy=Mb) • [x,Ax,MAx,AMAx,MAMAx,…,(MA)kx] – Keep last vector Akx only • Jacobi, Gauss Seidel Design Space for [x,Ax,…,Akx] (2/3) • Representation of sparse A – Zero pattern may be explicit or implicit – Nonzero entries may be explicit or implicit • Implicit save memory, communication Explicit pattern Implicit pattern Explicit nonzeros General sparse matrix Image segmentation Implicit nonzeros Laplacian(graph), ―Stencil matrix‖ for graph partitioning Ex: tridiag(-1,2,-1) • Representation of dense preconditioners M – Low rank off-diagonal blocks (semiseparable) Design Space for [x,Ax,…,Akx] (3/3) • Parallel implementation – From simple indexing, with redundant flops surface/volume ratio – To complicated indexing, with fewer redundant flops • Sequential implementation – Depends on whether vectors fit in fast memory • Reordering rows, columns of A – Important in parallel and sequential cases – Can be reduced to pair of Traveling Salesmen Problems • Plus all the optimizations for one SpMV! Summary • Communication-Avoiding Linear Algebra (CALA) • Lots of related work – Some going back to 1960‘s – Reports discuss this comprehensively, not here • Our contributions – Several new algorithms, improvements on old ones – Unifying parallel and sequential approaches to avoiding communication – Time for these algorithms has come, because of growing communication costs – Why avoid communication just for linear algebra motifs? Possible Class Projects • Come to BEBOP meetings (T 9 – 10:30, 606 Soda) • Incorporate multicore optimizations into OSKI • Experiment with SpMV on GPU – Which optimizations are most effective? • Try to speed up particular matrices of interest – Data mining • Experiment with new [x,Ax,…,Akx] kernel – GPU, multicore, distributed memory – On matrices of interest – Experiment with new solvers using this kernel Extra Slides Tuning Higher Level Algorithms • So far we have tuned a single sparse matrix kernel • y = AT*A*x motivated by higher level algorithm (SVD) • What can we do by extending tuning to a higher level? • Consider Krylov subspace methods for Ax=b, Ax = lx • Conjugate Gradients (CG), GMRES, Lanczos, … • Inner loop does y=A*x, dot products, saxpys, scalar ops • Inner loop costs at least O(1) messages • k iterations cost at least O(k) messages • Our goal: show how to do k iterations with O(1) messages • Possible payoff – make Krylov subspace methods much faster on machines with slow networks • Memory bandwidth improvements too (not discussed) • Obstacles: numerical stability, preconditioning, … Krylov Subspace Methods for Solving Ax=b • Compute a basis for a subspace V by doing y = A*x k times • Find “best” solution in that Krylov subspace V • Given starting vector x1, V spanned by x2 = A*x1, x3 = A*x2 , … , xk = A*xk-1 • GMRES finds an orthogonal basis of V by Gram-Schmidt, so it actually does a different set of SpMVs than in last bullet Example: Standard GMRES r = b - A*x1, = length(r), v1 = r / … length(r) = sqrt(S ri2 ) for m= 1 to k do w = A*vm … at least O(1) messages for i = 1 to m do … Gram-Schmidt him = dotproduct(vi , w ) … at least O(1) messages, or log(p) w = w – h im * vi end for hm+1,m = length(w) … at least O(1) messages, or log(p) vm+1 = w / hm+1,m end for find y minimizing length( Hk * y – e1 ) … small, local problem new x = x1 + Vk * y … Vk = [v1 , v2 , … , vk ] O(k2), or O(k2 log p), messages altogether Latency-Avoiding GMRES (1) r = b - A*x1, = length(r), v1 = r / … O(log p) messages Wk+1 = [v1 , A * v1 , A2 * v1 , … , Ak * v1 ] … O(1) messages [Q, R] = qr(Wk+1) … QR decomposition, O(log p) messages Hk = R(:, 2:k+1) * (R(1:k,1:k))-1 … small, local problem find y minimizing length( Hk * y – e1 ) … small, local problem new x = x1 + Qk * y … local problem O(log p) messages altogether Independent of k Latency-Avoiding GMRES (2) [Q, R] = qr(Wk+1) … QR decomposition, O(log p) messages Easy, but not so stable way to do it: X(myproc) = Wk+1T(myproc) * Wk+1 (myproc) … local computation Y = sum_reduction(X(myproc)) … O(log p) messages … Y = Wk+1T* Wk+1 R = (cholesky(Y))T … small, local computation Q(myproc) = Wk+1 (myproc) * R-1 … local computation Numerical example (1) Diagonal matrix with n=1000, Aii from 1 down to 10-5 Instability as k grows, after many iterations Numerical Example (2) Partial remedy: restarting periodically (every 120 iterations) Other remedies: high precision, different basis than [v , A * v , … , Ak * v ] Operation Counts for [Ax,A2x,A3x,…,Akx] on p procs Problem Per-proc cost Standard Optimized 1D mesh #messages 2k 2 (tridiagonal) #words sent 2k 2k #flops 5kn/p 5kn/p + 5k2 memory (k+4)n/p (k+4)n/p + 8k 3D mesh #messages 26k 26 27 pt stencil #words sent 6kn2p-2/3 + 12knp-1/3 6kn2p-2/3 + 12k2np-1/3 + O(k) + O(k3) #flops 53kn3/p 53kn3/p + O(k2n2p-2/3) memory (k+28)n3/p+ 6n2p-2/3 (k+28)n3/p + 168kn2p-2/3 + O(np-1/3) + O(k2np-1/3) Summary and Future Work • Dense – LAPACK – ScaLAPACK • Communication primitives • Sparse – Kernels, Stencils – Higher level algorithms • All of the above on new architectures – Vector, SMPs, multicore, Cell, … • High level support for tuning – Specification language – Integration into compilers Extra Slides A Sparse Matrix You Encounter Every Day Who am I? I am a Big Repository Of useful And useless Facts alike. Who am I? (Hint: Not your e-mail inbox.) What about the Google Matrix? • Google approach – Approx. once a month: rank all pages using connectivity structure • Find dominant eigenvector of a matrix – At query-time: return list of pages ordered by rank • Matrix: A = aG + (1-a)(1/n)uuT – Markov model: Surfer follows link with probability a, jumps to a random page with probability 1-a – G is n x n connectivity matrix [n billions] • gij is non-zero if page i links to page j • Normalized so each column sums to 1 • Very sparse: about 7—8 non-zeros per row (power law dist.) – u is a vector of all 1 values – Steady-state probability xi of landing on page i is solution to x = Ax • Approximate x by power method: x = Akx0 – In practice, k 25 Current Work • Public software release • Impact on library designs: Sparse BLAS, Trilinos, PETSc, … • Integration in large-scale applications – DOE: Accelerator design; plasma physics – Geophysical simulation based on Block Lanczos (ATA*X; LBL) • Systematic heuristics for data structure selection? • Evaluation of emerging architectures – Revisiting vector micros • Other sparse kernels – Matrix triple products, Ak*x • Parallelism • Sparse benchmarks (with UTK) [Gahvari & Hoemmen] • Automatic tuning of MPI collective ops [Nishtala, et al.] Summary of High-Level Themes • ―Kernel-centric‖ optimization – Vs. basic block, trace, path optimization, for instance – Aggressive use of domain-specific knowledge • Performance bounds modeling – Evaluating software quality – Architectural characterizations and consequences • Empirical search – Hybrid on-line/run-time models • Statistical performance models – Exploit information from sampling, measuring Related Work • My bibliography: 337 entries so far • Sample area 1: Code generation – Generative & generic programming – Sparse compilers – Domain-specific generators • Sample area 2: Empirical search-based tuning – Kernel-centric • linear algebra, signal processing, sorting, MPI, … – Compiler-centric • profiling + FDO, iterative compilation, superoptimizers, self- tuning compilers, continuous program optimization Future Directions (A Bag of Flaky Ideas) • Composable code generators and search spaces • New application domains – PageRank: multilevel block algorithms for topic-sensitive search? • New kernels: cryptokernels – rich mathematical structure germane to performance; lots of hardware • New tuning environments – Parallel, Grid, ―whole systems‖ • Statistical models of application performance – Statistical learning of concise parametric models from traces for architectural evaluation – Compiler/automatic derivation of parametric models Possible Future Work • Different Architectures – New FP instruction sets (SSE2) – SMP / multicore platforms – Vector architectures • Different Kernels • Higher Level Algorithms – Parallel versions of kenels, with optimized communication – Block algorithms (eg Lanczos) – XBLAS (extra precision) • Produce Benchmarks – Augment HPCC Benchmark • Make it possible to combine optimizations easily for any kernel • Related tuning activities (LAPACK & ScaLAPACK) Review of Tuning by Illustration (Extra Slides) Splitting for Variable Blocks and Diagonals • Decompose A = A1 + A2 + … At – Detect ―canonical‖ structures (sampling) – Split – Tune each Ai – Improve performance and save storage • New data structures – Unaligned block CSR • Relax alignment in rows & columns – Row-segmented diagonals Example: Variable Block Row (Matrix #12) 2.1x over CSR 1.8x over RB Example: Row-Segmented Diagonals 2x over CSR Mixed Diagonal and Block Structure Summary • Automated block size selection – Empirical modeling and search – Register blocking for SpMV, triangular solve, ATA*x • Not fully automated – Given a matrix, select splittings and transformations • Lots of combinatorial problems – TSP reordering to create dense blocks (Pinar ‘97; Moon, et al. ‘04) Extra Slides A Sparse Matrix You Encounter Every Day Who am I? I am a Big Repository Of useful And useless Facts alike. Who am I? (Hint: Not your e-mail inbox.) Problem Context • Sparse kernels abound – Models of buildings, cars, bridges, economies, … – Google PageRank algorithm • Historical trends – Sparse matrix-vector multiply (SpMV): 10% of peak – 2x faster with ―hand-tuning‖ – Tuning becoming more difficult over time – Promise of automatic tuning: PHiPAC/ATLAS, FFTW, … • Challenges to high-performance – Not dense linear algebra! • Complex data structures: indirect, irregular memory access • Performance depends strongly on run-time inputs Key Questions, Ideas, Conclusions • How to tune basic sparse kernels automatically? – Empirical modeling and search • Up to 4x speedups for SpMV • 1.8x for triangular solve • 4x for ATA*x; 2x for A2*x • 7x for multiple vectors • What are the fundamental limits on performance? – Kernel-, machine-, and matrix-specific upper bounds • Achieve 75% or more for SpMV, limiting low-level tuning • Consequences for architecture? • General techniques for empirical search-based tuning? – Statistical models of performance Road Map • Sparse matrix-vector multiply (SpMV) in a nutshell • Historical trends and the need for search • Automatic tuning techniques • Upper bounds on performance • Statistical models of performance Compressed Sparse Row (CSR) Storage Matrix-vector multiply kernel: y(i) y(i) + A(i,j)*x(j) for each row i for k=ptr[i] to ptr[i+1] do y[i] = y[i] + val[k]*x[ind[k]] Road Map • Sparse matrix-vector multiply (SpMV) in a nutshell • Historical trends and the need for search • Automatic tuning techniques • Upper bounds on performance • Statistical models of performance Historical Trends in SpMV Performance • The Data – Uniprocessor SpMV performance since 1987 – ―Untuned‖ and ―Tuned‖ implementations – Cache-based superscalar micros; some vectors – LINPACK SpMV Historical Trends: Mflop/s Example: The Difficulty of Tuning • n = 21216 • nnz = 1.5 M • kernel: SpMV • Source: NASA structural analysis problem Still More Surprises • More complicated non-zero structure in general Still More Surprises • More complicated non-zero structure in general • Example: 3x3 blocking – Logical grid of 3x3 cells Historical Trends: Mixed News • Observations – Good news: Moore‘s law like behavior – Bad news: ―Untuned‖ is 10% peak or less, worsening – Good news: ―Tuned‖ roughly 2x better today, and improving – Bad news: Tuning is complex – (Not really news: SpMV is not LINPACK) • Questions – Application: Automatic tuning? – Architect: What machines are good for SpMV? Road Map • Sparse matrix-vector multiply (SpMV) in a nutshell • Historical trends and the need for search • Automatic tuning techniques – SpMV [SC‘02; IJHPCA ‘04b] – Sparse triangular solve (SpTS) [ICS/POHLL ‘02] – ATA*x [ICCS/WoPLA ‘03] • Upper bounds on performance • Statistical models of performance SPARSITY: Framework for Tuning SpMV • SPARSITY: Automatic tuning for SpMV [Im & Yelick ‘99] – General approach • Identify and generate implementation space • Search space using empirical models & experiments – Prototype library and heuristic for choosing register block size • Also: cache-level blocking, multiple vectors • What‘s new? – New block size selection heuristic • Within 10% of optimal — replaces previous version – Expanded implementation space • Variable block splitting, diagonals, combinations – New kernels: sparse triangular solve, ATA*x, Ar*x Automatic Register Block Size Selection • Selecting the r x c block size – Off-line benchmark: characterize the machine • Precompute Mflops(r,c) using dense matrix for each r x c • Once per machine/architecture – Run-time “search”: characterize the matrix • Sample A to estimate Fill(r,c) for each r x c – Run-time heuristic model • Choose r, c to maximize Mflops(r,c) / Fill(r,c) • Run-time costs – Up to ~40 SpMVs (empirical worst case) DGEMV Accuracy of the Tuning Heuristics (1/4) NOTE: “Fair” flops used (ops on explicit zeros not counted as “work”) DGEMV Accuracy of the Tuning Heuristics (2/4) DGEMV Accuracy of the Tuning Heuristics (3/4) DGEMV Accuracy of the Tuning Heuristics (4/4) Road Map • Sparse matrix-vector multiply (SpMV) in a nutshell • Historical trends and the need for search • Automatic tuning techniques • Upper bounds on performance – SC’02 • Statistical models of performance Motivation for Upper Bounds Model • Questions – Speedups are good, but what is the speed limit? • Independent of instruction scheduling, selection – What machines are ―good‖ for SpMV? Upper Bounds on Performance: Blocked SpMV • P = (flops) / (time) – Flops = 2 * nnz(A) • Lower bound on time: Two main assumptions – 1. Count memory ops only (streaming) – 2. Count only compulsory, capacity misses: ignore conflicts • Account for line sizes • Account for matrix size and nnz • Charge min access ―latency‖ ai at Li cache & amem – e.g., Saavedra-Barrera and PMaC MAPS benchmarks Time a i Hits i a mem Hits mem i 1 a1 Loads (a i 1 a i ) Misses i (a mem a ) Misses i 1 Example: Bounds on Itanium 2 Example: Bounds on Itanium 2 Example: Bounds on Itanium 2 Fraction of Upper Bound Across Platforms Achieved Performance and Machine Balance • Machine balance [Callahan ‘88; McCalpin ‘95] – Balance = Peak Flop Rate / Bandwidth (flops / double) • Ideal balance for mat-vec: 2 flops / double – For SpMV, even less Time a1 Loads (a i 1 a i ) Misses i (a mem a ) Misses i • SpMV ~ streaming – 1 / (avg load time to stream 1 array) ~ (bandwidth) – ―Sustained‖ balance = peak flops / model bandwidth Where Does the Time Go? Time a i Hits i a mem Hits mem i 1 • Most time assigned to memory • Caches ―disappear‖ when line sizes are equal – Strictly increasing line sizes Execution Time Breakdown: Matrix 40 Speedups with Increasing Line Size Summary: Performance Upper Bounds • What is the best we can do for SpMV? – Limits to low-level tuning of blocked implementations – Refinements? • What machines are good for SpMV? – Partial answer: balance characterization • Architectural consequences? – Example: Strictly increasing line sizes Road Map • Sparse matrix-vector multiply (SpMV) in a nutshell • Historical trends and the need for search • Automatic tuning techniques • Upper bounds on performance • Tuning other sparse kernels • Statistical models of performance – FDO ’00; IJHPCA ’04a Statistical Models for Automatic Tuning • Idea 1: Statistical criterion for stopping a search – A general search model • Generate implementation • Measure performance • Repeat – Stop when probability of being within e of optimal falls below threshold • Can estimate distribution on-line • Idea 2: Statistical performance models – Problem: Choose 1 among m implementations at run-time – Sample performance off-line, build statistical model Example: Select a Matmul Implementation Example: Support Vector Classification Road Map • Sparse matrix-vector multiply (SpMV) in a nutshell • Historical trends and the need for search • Automatic tuning techniques • Upper bounds on performance • Tuning other sparse kernels • Statistical models of performance • Summary and Future Work Summary of High-Level Themes • ―Kernel-centric‖ optimization – Vs. basic block, trace, path optimization, for instance – Aggressive use of domain-specific knowledge • Performance bounds modeling – Evaluating software quality – Architectural characterizations and consequences • Empirical search – Hybrid on-line/run-time models • Statistical performance models – Exploit information from sampling, measuring Related Work • My bibliography: 337 entries so far • Sample area 1: Code generation – Generative & generic programming – Sparse compilers – Domain-specific generators • Sample area 2: Empirical search-based tuning – Kernel-centric • linear algebra, signal processing, sorting, MPI, … – Compiler-centric • profiling + FDO, iterative compilation, superoptimizers, self- tuning compilers, continuous program optimization Future Directions (A Bag of Flaky Ideas) • Composable code generators and search spaces • New application domains – PageRank: multilevel block algorithms for topic-sensitive search? • New kernels: cryptokernels – rich mathematical structure germane to performance; lots of hardware • New tuning environments – Parallel, Grid, ―whole systems‖ • Statistical models of application performance – Statistical learning of concise parametric models from traces for architectural evaluation – Compiler/automatic derivation of parametric models Acknowledgements • Super-advisors: Jim and Kathy • Undergraduate R.A.s: Attila, Ben, Jen, Jin, Michael, Rajesh, Shoaib, Sriram, Tuyet-Linh • See pages xvi—xvii of dissertation. TSP-based Reordering: Before (Pinar ‘97; Moon, et al ‗04) TSP-based Reordering: After (Pinar ‘97; Moon, et al ‗04) Up to 2x speedups over CSR Example: L2 Misses on Itanium 2 Misses measured using PAPI [Browne ’00] Example: Distribution of Blocked Non-Zeros Sparse/Dense Partitioning for SpTS • Partition L into sparse (L1,L2) and dense LD: L1 x1 b1 L 2 LD x2 b2 • Perform SpTS in three steps: (1) L1 x1 b1 (2) ˆ b2 b2 L2 x1 (3) L x b ˆ D 2 2 • Sparsity optimizations for (1)—(2); DTRSV for (3) • Tuning parameters: block size, size of dense triangle SpTS Performance: Power3 Summary of SpTS and AAT*x Results • SpTS — Similar to SpMV – 1.8x speedups; limited benefit from low-level tuning • AATx, ATAx – Cache interleaving only: up to 1.6x speedups – Reg + cache: up to 4x speedups • 1.8x speedup over register only – Similar heuristic; same accuracy (~ 10% optimal) – Further from upper bounds: 60—80% • Opportunity for better low-level tuning a la PHiPAC/ATLAS • Matrix triple products? Ak*x? – Preliminary work Register Blocking: Speedup Register Blocking: Performance Register Blocking: Fraction of Peak Example: Confidence Interval Estimation Costs of Tuning Splitting + UBCSR: Pentium III Splitting + UBCSR: Power4 Splitting+UBCSR Storage: Power4 Example: Variable Block Row (Matrix #13) Preliminary Results (Matrix Set 2): Itanium 2 Web/IR Dense FEM FEM (var) Bio Econ Stat LP Multiple Vector Performance What about the Google Matrix? • Google approach – Approx. once a month: rank all pages using connectivity structure • Find dominant eigenvector of a matrix – At query-time: return list of pages ordered by rank • Matrix: A = aG + (1-a)(1/n)uuT – Markov model: Surfer follows link with probability a, jumps to a random page with probability 1-a – G is n x n connectivity matrix [n 3 billion] • gij is non-zero if page i links to page j • Normalized so each column sums to 1 • Very sparse: about 7—8 non-zeros per row (power law dist.) – u is a vector of all 1 values – Steady-state probability xi of landing on page i is solution to x = Ax • Approximate x by power method: x = Akx0 – In practice, k 25 MAPS Benchmark Example: Power4 MAPS Benchmark Example: Itanium 2 Saavedra-Barrera Example: Ultra 2i Execution Time Breakdown (PAPI): Matrix 40 Preliminary Results (Matrix Set 1): Itanium 2 Dense FEM FEM (var) Assorted LP Tuning Sparse Triangular Solve (SpTS) • Compute x=L-1*b where L sparse lower triangular, x & b dense • L from sparse LU has rich dense substructure – Dense trailing triangle can account for 20—90% of matrix non-zeros • SpTS optimizations – Split into sparse trapezoid and dense trailing triangle – Use tuned dense BLAS (DTRSV) on dense triangle – Use Sparsity register blocking on sparse part • Tuning parameters – Size of dense trailing triangle – Register block size Sparse Kernels and Optimizations • Kernels – Sparse matrix-vector multiply (SpMV): y=A*x – Sparse triangular solve (SpTS): x=T-1*b – y=AAT*x, y=ATA*x – Powers (y=Ak*x), sparse triple-product (R*A*RT), … • Optimization techniques (implementation space) – Register blocking – Cache blocking – Multiple dense vectors (x) – A has special structure (e.g., symmetric, banded, …) – Hybrid data structures (e.g., splitting, switch-to-dense, …) – Matrix reordering • How and when do we search? – Off-line: Benchmark implementations – Run-time: Estimate matrix properties, evaluate performance models based on benchmark data Cache Blocked SpMV on LSI Matrix: Ultra 2i A 10k x 255k 3.7M non-zeros Baseline: 16 Mflop/s Best block size & performance: 16k x 64k 28 Mflop/s Cache Blocking on LSI Matrix: Pentium 4 A 10k x 255k 3.7M non-zeros Baseline: 44 Mflop/s Best block size & performance: 16k x 16k 210 Mflop/s Cache Blocked SpMV on LSI Matrix: Itanium A 10k x 255k 3.7M non-zeros Baseline: 25 Mflop/s Best block size & performance: 16k x 32k 72 Mflop/s Cache Blocked SpMV on LSI Matrix: Itanium 2 A 10k x 255k 3.7M non-zeros Baseline: 170 Mflop/s Best block size & performance: 16k x 65k 275 Mflop/s Inter-Iteration Sparse Tiling (1/3) x1 t1 y1 • [Strout, et al., ‗01] • Let A be 6x6 tridiagonal • Consider y=A2x x2 t2 y2 – t=Ax, y=At • Nodes: vector elements x3 t3 y3 • Edges: matrix elements aij x4 t4 y4 x5 t5 y5 Inter-Iteration Sparse Tiling (2/3) x1 t1 y1 • [Strout, et al., ‗01] • Let A be 6x6 tridiagonal • Consider y=A2x x2 t2 y2 – t=Ax, y=At • Nodes: vector elements x3 t3 y3 • Edges: matrix elements aij • Orange = everything needed x4 t4 y4 to compute y1 – Reuse a11, a12 x5 t5 y5 Inter-Iteration Sparse Tiling (3/3) x1 t1 y1 • [Strout, et al., ‗01] • Let A be 6x6 tridiagonal • Consider y=A2x x2 t2 y2 – t=Ax, y=At • Nodes: vector elements x3 t3 y3 • Edges: matrix elements aij • Orange = everything needed x4 t4 y4 to compute y1 – Reuse a11, a12 x5 t5 y5 • Grey = y2, y3 – Reuse a23, a33, a43 Inter-Iteration Sparse Tiling: Issues x1 t1 y1 • Tile sizes (colored regions) grow with no. of iterations and increasing out-degree x2 t2 y2 – G likely to have a few nodes with high out-degree (e.g., Yahoo) x3 t3 y3 • Mathematical tricks to limit tile size? – Judicious dropping of edges x4 t4 y4 [Ng‘01] x5 t5 y5 Summary and Questions • Need to understand matrix structure and machine – BeBOP: suite of techniques to deal with different sparse structures and architectures • Google matrix problem – Established techniques within an iteration – Ideas for inter-iteration optimizations – Mathematical structure of problem may help • Questions – Structure of G? – What are the computational bottlenecks? – Enabling future computations? • E.g., topic-sensitive PageRank multiple vector version [Haveliwala ‘02] – See www.cs.berkeley.edu/~richie/bebop/intel/google for more info, including more complete Itanium 2 results. Exploiting Matrix Structure • Symmetry (numerical or structural) – Reuse matrix entries – Can combine with register blocking, multiple vectors, … • Matrix splitting – Split the matrix, e.g., into r x c and 1 x 1 – No fill overhead • Large matrices with random structure – E.g., Latent Semantic Indexing (LSI) matrices – Technique: cache blocking • Store matrix as 2i x 2j sparse submatrices • Effective when x vector is large • Currently, search to find fastest size Symmetric SpMV Performance: Pentium 4 SpMV with Split Matrices: Ultra 2i Cache Blocking on Random Matrices: Itanium Speedup on four banded random matrices. Sparse Kernels and Optimizations • Kernels – Sparse matrix-vector multiply (SpMV): y=A*x – Sparse triangular solve (SpTS): x=T-1*b – y=AAT*x, y=ATA*x – Powers (y=Ak*x), sparse triple-product (R*A*RT), … • Optimization techniques (implementation space) – Register blocking – Cache blocking – Multiple dense vectors (x) – A has special structure (e.g., symmetric, banded, …) – Hybrid data structures (e.g., splitting, switch-to-dense, …) – Matrix reordering • How and when do we search? – Off-line: Benchmark implementations – Run-time: Estimate matrix properties, evaluate performance models based on benchmark data Register Blocked SpMV: Pentium III Register Blocked SpMV: Ultra 2i Register Blocked SpMV: Power3 Register Blocked SpMV: Itanium Possible Optimization Techniques • Within an iteration, i.e., computing (G+uuT)*x once – Cache block G*x • On linear programming matrices and matrices with random structure (e.g., LSI), 1.5—4x speedups • Best block size is matrix and machine dependent – Reordering and/or splitting of G to separate dense structure (rows, columns, blocks) • Between iterations, e.g., (G+uuT)2x – (G+uuT)2x = G2x + (Gu)uTx + u(uTG)x + u(uTu)uTx • Compute Gu, uTG, uTu once for all iterations • G2x: Inter-iteration tiling to read G only once Multiple Vector Performance: Itanium Sparse Kernels and Optimizations • Kernels – Sparse matrix-vector multiply (SpMV): y=A*x – Sparse triangular solve (SpTS): x=T-1*b – y=AAT*x, y=ATA*x – Powers (y=Ak*x), sparse triple-product (R*A*RT), … • Optimization techniques (implementation space) – Register blocking – Cache blocking – Multiple dense vectors (x) – A has special structure (e.g., symmetric, banded, …) – Hybrid data structures (e.g., splitting, switch-to-dense, …) – Matrix reordering • How and when do we search? – Off-line: Benchmark implementations – Run-time: Estimate matrix properties, evaluate performance models based on benchmark data SpTS Performance: Itanium (See POHLL ‘02 workshop paper, at ICS ‘02.) Sparse Kernels and Optimizations • Kernels – Sparse matrix-vector multiply (SpMV): y=A*x – Sparse triangular solve (SpTS): x=T-1*b – y=AAT*x, y=ATA*x – Powers (y=Ak*x), sparse triple-product (R*A*RT), … • Optimization techniques (implementation space) – Register blocking – Cache blocking – Multiple dense vectors (x) – A has special structure (e.g., symmetric, banded, …) – Hybrid data structures (e.g., splitting, switch-to-dense, …) – Matrix reordering • How and when do we search? – Off-line: Benchmark implementations – Run-time: Estimate matrix properties, evaluate performance models based on benchmark data Optimizing AAT*x • Kernel: y=AAT*x, where A is sparse, x & y dense – Arises in linear programming, computation of SVD – Conventional implementation: compute z=AT*x, y=A*z • Elements of A can be reused: a1 T n y a1 an x ak ( ak x ) T aT k 1 n • When ak represent blocks of columns, can apply register blocking. Optimized AAT*x Performance: Pentium III Current Directions • Applying new optimizations – Other split data structures (variable block, diagonal, …) – Matrix reordering to create block structure – Structural symmetry • New kernels (triple product RART, powers Ak, …) • Tuning parameter selection • Building an automatically tuned sparse matrix library – Extending the Sparse BLAS – Leverage existing sparse compilers as code generation infrastructure – More thoughts on this topic tomorrow Related Work • Automatic performance tuning systems – PHiPAC [Bilmes, et al., ‘97], ATLAS [Whaley & Dongarra ‘98] – FFTW [Frigo & Johnson ‘98], SPIRAL [Pueschel, et al., ‘00], UHFFT [Mirkovic and Johnsson ‘00] – MPI collective operations [Vadhiyar & Dongarra ‘01] • Code generation – FLAME [Gunnels & van de Geijn, ‘01] – Sparse compilers: [Bik ‘99], Bernoulli [Pingali, et al., ‘97] – Generic programming: Blitz++ [Veldhuizen ‘98], MTL [Siek & Lumsdaine ‘98], GMCL [Czarnecki, et al. ‘98], … • Sparse performance modeling – [Temam & Jalby ‘92], [White & Saddayappan ‘97], [Navarro, et al., ‘96], [Heras, et al., ‘99], [Fraguela, et al., ‘99], … More Related Work • Compiler analysis, models – CROPS [Carter, Ferrante, et al.]; Serial sparse tiling [Strout ‘01] – TUNE [Chatterjee, et al.] – Iterative compilation [O‘Boyle, et al., ‘98] – Broadway compiler [Guyer & Lin, ‘99] – [Brewer ‘95], ADAPT [Voss ‘00] • Sparse BLAS interfaces – BLAST Forum (Chapter 3) – NIST Sparse BLAS [Remington & Pozo ‘94]; SparseLib++ – SPARSKIT [Saad ‘94] – Parallel Sparse BLAS [Fillipone, et al. ‘96] Context: Creating High-Performance Libraries • Application performance dominated by a few computational kernels • Today: Kernels hand-tuned by vendor or user • Performance tuning challenges – Performance is a complicated function of kernel, architecture, compiler, and workload – Tedious and time-consuming • Successful automated approaches – Dense linear algebra: ATLAS/PHiPAC – Signal processing: FFTW/SPIRAL/UHFFT Cache Blocked SpMV on LSI Matrix: Itanium Sustainable Memory Bandwidth Multiple Vector Performance: Pentium 4 Multiple Vector Performance: Itanium Multiple Vector Performance: Pentium 4 Optimized AAT*x Performance: Ultra 2i Optimized AAT*x Performance: Pentium 4 High Precision GEMV (XBLAS) High Precision Algorithms (XBLAS) • Double-double (High precision word represented as pair of doubles) – Many variations on these algorithms; we currently use Bailey‘s • Exploiting Extra-wide Registers – Suppose s(1) , … , s(n) have f-bit fractions, SUM has F>f bit fraction – Consider following algorithm for S = Si=1,n s(i) • Sort so that |s(1)| |s(2)| … |s(n)| • SUM = 0, for i = 1 to n SUM = SUM + s(i), end for, sum = SUM – Theorem (D., Hida) Suppose F<2f (less than double precision) • If n 2F-f + 1, then error 1.5 ulps • If n = 2F-f + 2, then error 22f-F ulps (can be 1) • If n 2F-f + 3, then error can be arbitrary (S 0 but sum = 0 ) – Examples • s(i) double (f=53), SUM double extended (F=64) – accurate if n 211 + 1 = 2049 • Dot product of single precision x(i) and y(i) – s(i) = x(i)*y(i) (f=2*24=48), SUM double extended (F=64) – accurate if n 216 + 1 = 65537 More Extra Slides Awards • Best Paper, Intern. Conf. Parallel Processing, 2004 – ―Performance models for evaluation and automatic performance tuning of symmetric sparse matrix-vector multiply‖ • Best Student Paper, Intern. Conf. Supercomputing, Workshop on Performance Optimization via High-Level Languages and Libraries, 2003 – Best Student Presentation too, to Richard Vuduc – ―Automatic performance tuning and analysis of sparse triangular solve‖ • Finalist, Best Student Paper, Supercomputing 2002 – To Richard Vuduc – ―Performance Optimization and Bounds for Sparse Matrix-vector Multiply‖ • Best Presentation Prize, MICRO-33: 3rd ACM Workshop on Feedback-Directed Dynamic Optimization, 2000 – To Richard Vuduc – ―Statistical Modeling of Feedback Data in an Automatic Tuning System‖ Accuracy of the Tuning Heuristics (4/4) DGEMV Can Match DGEMV Performance Fraction of Upper Bound Across Platforms Achieved Performance and Machine Balance • Machine balance [Callahan ‘88; McCalpin ‘95] – Balance = Peak Flop Rate / Bandwidth (flops / double) – Lower is better (I.e. can hope to get higher fraction of peak flop rate) • Ideal balance for mat-vec: 2 flops / double – For SpMV, even less Where Does the Time Go? Time a i Hits i a mem Hits mem i 1 • Most time assigned to memory • Caches ―disappear‖ when line sizes are equal – Strictly increasing line sizes Execution Time Breakdown: Matrix 40 Execution Time Breakdown (PAPI): Matrix 40 Speedups with Increasing Line Size Summary: Performance Upper Bounds • What is the best we can do for SpMV? – Limits to low-level tuning of blocked implementations – Refinements? • What machines are good for SpMV? – Partial answer: balance characterization • Architectural consequences? – Help to have strictly increasing line sizes Evaluating algorithms and machines for SpMV • Metrics – Speedups – Mflop/s (―fair‖ flops) – Fraction of peak • Questions – Speedups are good, but what is ―the best?‖ • Independent of instruction scheduling, selection • Can SpMV be further improved or not? – What machines are ―good‖ for SpMV? Statistical Models for Automatic Tuning • Idea 1: Statistical criterion for stopping a search – A general search model • Generate implementation • Measure performance • Repeat – Stop when probability of being within e of optimal falls below threshold • Can estimate distribution on-line • Idea 2: Statistical performance models – Problem: Choose 1 among m implementations at run-time – Sample performance off-line, build statistical model SpMV Historical Trends: Fraction of Peak Motivation for Automatic Performance Tuning of SpMV • Historical trends – Sparse matrix-vector multiply (SpMV): 10% of peak or less • Performance depends on machine, kernel, matrix – Matrix known at run-time – Best data structure + implementation can be surprising • Our approach: empirical performance modeling and algorithm search Sparse CALA Summary • Optimizing one SpMV too limiting • Take k steps of Krylov subspace method – GMRES, CG, Lanczos, Arnoldi – Assume matrix ―well-partitioned,‖ with modest surface-to-volume ratio – Parallel implementation • Conventional: O(k log p) messages • CALA: O(log p) messages - optimal – Serial implementation • Conventional: O(k) moves of data from slow to fast memory • CALA: O(1) moves of data – optimal – Can incorporate some preconditioners • Need to be able to ―compress‖ interactions between distant i, j • Hierarchical, semiseparable matrices … – Lots of speed up possible (measured and modeled) • Price: some redundant computation Motivation for Automatic Performance Tuning • Writing high performance software is hard – Make programming easier while getting high speed • Ideal: program in your favorite high level language (Matlab, Python, PETSc…) and get a high fraction of peak performance • Reality: Best algorithm (and its implementation) can depend strongly on the problem, computer architecture, compiler,… – Best choice can depend on knowing a lot of applied mathematics and computer science • How much of this can we teach? • How much of this can we automate? Examples of Automatic Performance Tuning (1) • Dense BLAS – Sequential – PHiPAC (UCB), then ATLAS (UTK) (used in Matlab) – math-atlas.sourceforge.net/ • Fast Fourier Transform (FFT) & variations – Sequential and Parallel – FFTW (MIT) – www.fftw.org • Digital Signal Processing – SPIRAL: www.spiral.net (CMU) • Communication Collectives (UCB, UTK) • Rose (LLNL), Bernoulli (Cornell), Telescoping Languages (Rice), … • More projects, conferences, government reports, … Potential Impact on Applications: T3P • Application: accelerator design [Ko] • 80% of time spent in SpMV • Relevant optimization techniques – Symmetric storage – Register blocking • On Single Processor Itanium 2 – 1.68x speedup • 532 Mflops, or 15% of 3.6 GFlop peak – 4.4x speedup with multiple (8) vectors • 1380 Mflops, or 38% of peak Examples of Automatic Performance Tuning (2) • What do dense BLAS, FFTs, signal processing, MPI reductions have in common? – Can do the tuning off-line: once per architecture, algorithm – Can take as much time as necessary (hours, a week…) – At run-time, algorithm choice may depend only on few parameters • Matrix dimension, size of FFT, etc. Ultra 2i - 9% 63 Mflop/s Ultra 3 - 5% 109 Mflop/s SpMV Performance (Matrix #2): Generation 2 35 Mflop/s 53 Mflop/s Pentium III - 19% 96 Mflop/s Pentium III-M - 15% 120 Mflop/s 42 Mflop/s 58 Mflop/s Power3 - 13% 195 Mflop/s Power4 - 14% 703 Mflop/s SpMV Performance (Matrix #2): Generation 1 100 Mflop/s 469 Mflop/s Itanium 1 - 7% 225 Mflop/s Itanium 2 - 31% 1.1 Gflop/s 103 Mflop/s 276 Mflop/s Taking advantage of block structure in SpMV • Bottleneck is time to get matrix from memory • Don‘t store each nonzero with index, instead store each nonzero r-by-c block with index – Storage drops by up to 2x, if rc >> 1, all 32-bit quantities – Time to fetch matrix from memory decreases • Change both data structure and algorithm – Need to pick r and c – Need to change algorithm accordingly • In example, is r=c=8 best choice? – Minimizes storage, so looks like a good idea… Example: L2 Misses on Itanium 2 Misses measured using PAPI [Browne ’00] Example: Sparse Triangular Factor • Raefsky4 (structural problem) + SuperLU + colmmd • N=19779, nnz=12.6 M Dense trailing triangle: dim=2268, 20% of total nz Can be as high as 90+%! 1.8x over CSR Cache Optimizations for AAT*x • Cache-level: Interleave multiplication by A, AT – Only fetch A from memory once a1T n … x AA x a1 an ai (ai x) T T T an … i 1 ―axpy‖ dot product • Register-level: aiT to be rc block row, or diag row Example: Combining Optimizations (1/2) • Register blocking, symmetry, multiple (k) vectors – Three low-level tuning parameters: r, c, v X k v * r c += Y A Example: Combining Optimizations (2/2) • Register blocking, symmetry, and multiple vectors [Ben Lee @ UCB] – Symmetric, blocked, 1 vector • Up to 2.6x over nonsymmetric, blocked, 1 vector – Symmetric, blocked, k vectors • Up to 2.1x over nonsymmetric, blocked, k vecs. • Up to 7.3x over nonsymmetric, nonblocked, 1, vector – Symmetric Storage: up to 64.7% savings How the OSKI Tunes (Overview) • At library build/install-time – Pre-generate and compile code variants into dynamic libraries – Collect benchmark data • Measures and records speed of possible sparse data structure and code variants on target architecture – Installation process uses standard, portable GNU AutoTools • At run-time – Library ―tunes‖ using heuristic models • Models analyze user‘s matrix & benchmark data to choose optimized data structure and code – Non-trivial tuning cost: up to ~40 mat-vecs • Library limits the time it spends tuning based on estimated workload – provided by user or inferred by library • User may reduce cost by saving tuning results for application on future runs with same or similar matrix Optimizations in OSKI, so far • Fully automatic heuristics for – Sparse matrix-vector multiply • Register-level blocking • Register-level blocking + symmetry + multiple vectors • Cache-level blocking – Sparse triangular solve with register-level blocking and ―switch-to-dense‖ optimization – Sparse ATA*x with register-level blocking • User may select other optimizations manually – Diagonal storage optimizations, reordering, splitting; tiled matrix powers kernel (Ak*x) – All available in dynamic libraries – Accessible via high-level embedded script language • ―Plug-in‖ extensibility – Very advanced users may write their own heuristics, create new data structures/code variants and dynamically add them to the system Performance Results • Measured – Sequential/OOC speedup up to 3x • Modeled – Sequential/multicore speedup up to 2.5x – Parallel/Petascale speedup up to 6.9x – Parallel/Grid speedup up to 22x • See bebop.cs.berkeley.edu/#pubs

DOCUMENT INFO

Shared By:

Categories:

Tags:

Stats:

views: | 2 |

posted: | 9/16/2011 |

language: | English |

pages: | 261 |

OTHER DOCS BY pengxiang

Docstoc is the premier online destination to start and grow small businesses. It hosts the best quality and widest selection of professional documents (over 20 million) and resources including expert videos, articles and productivity tools to make every small business better.

Search or Browse for any specific document or resource you need for your business. Or explore our curated resources for Starting a Business, Growing a Business or for Professional Development.

Feel free to Contact Us with any questions you might have.