An Optimized GPU Implementation of a Multicolor Point-Implicit - - PowerPoint PPT Presentation

An Optimized GPU Implementation of a Multicolor Point-Implicit Linear Solver

Outline

• Motivation
• Problem and Challenge
• Approach
• Performance Results
• Future Work and Conclusion

Motivation

FUN3D

• FUN3D is a suite of computational fluid dynamics software developed

at the NASA Langley Research Center to solve the Navier-Stokes equations for a broad range of aerodynamics applications

• For high-energy flows, traditional perfect gas assumptions are invalid

and the system must be further expanded to accommodate finite-rate chemistry models. For these reasons, accurate and efficient simulations of complex aerodynamic flows are challenging and require significant computational resources.

FUN3D

• FUN3D uses an implicit time-integration strategy on unstructured grid

that requires that requires solution of large tightly-coupled system of block-sparse linear equations at each time step

• Computing solutions to these linear systems accounts for a significant

fraction of the overall runtime in virtually all FUN3D simulations.

Problem and Challenges

Linear Solver

• Implicit scheme results in linear systems of equations:

𝑩𝒚 = 𝒄, 𝑩 is a 𝑜 × 𝑜 block matrix for 𝑜 grid points; block is of size 𝑜𝑐 × 𝑜𝑐

• Matrix 𝑩 is segregated into two separate matrices:

𝑩 ≡ 𝑬 + 𝑷, where 𝑬 and 𝑷 represent the diagonal and off-diagonal blocks of 𝑩

• Arrays for storing 𝑬: 𝐸 𝑜𝑐, 𝑜𝑐, 𝑜

Prior to performing each linear solve, each diagonal block 𝐸𝑗 is decomposed in-place into lower and upper triangular matrices 𝑀𝑗 𝑏𝑜𝑒 𝑉𝑗, 1 ≤ 𝑗 ≤ 𝑜

Sparse Block-Matrix Storage Format for 𝑷

Multicolor Solver

• Grid points are colored such that no two adjacent points are assigned the

same color.

• Unknowns of the same color carry no data dependency and may be

updated in parallel.

• Updates of unknowns of each color use the latest updated values at grid

points corresponding to other colors.

• The overall process may be repeated using several outer sweeps over the

entire system.

Distribution of non-zeros in the sample matrix

Solver Challenges

• With an average of few off-diagonal blocks per row, the arithmetic intensity of the

computation is quite low (≈ 0.5 for 983633 sized matrix) – memory bound on GPU

• The number of rows associated with a color can vary significantly. Consequently the

amount of parallelism available for different colors varies significantly

• To support strong scalability, the single node performance for light workload should also

be good

• The solver uses indirect memory addressing
• The implementation must support a broad range of block sizes

Solver Computation

Solver Computation – MPI Implementation

subroutine fun3d_solver do color = 1, nc ! Process rows of each color MPI Communication: update portion of x array end do end subroutine fun3d_solver

This computation is identical to the sequential computation

Approach

Solver Computation

Sparse block-matrix (A) vector operation Multiple triangular solve using dense blocks of matrix D

𝒓 𝑬 𝒚

𝑔𝑝𝑠𝑥𝑏𝑠𝑒 & 𝑐𝑏𝑑𝑙𝑥𝑏𝑠𝑒

Naive Approach

Assign a thread to process a row block. As all rows can be processed independently the parallelism is determined by number of row blocks.

𝒚 𝑷 𝑷𝒚

Solver Computation

Use sparse BLAS: cusparseSbsrmv Use batched dense BLAS: cublasStrsmBatched

Sparse and Dense Blas Formulation

• Advantage: makes code portable with performance on different

architectures

• Issues: For block sizes of interest BLAS routines do not perform well.
• Solution: Write optimized BLAS for block sizes of interest

Overall Structure of Optimized Solver

Approach for sparse matrix vector (for nb = 1 to 16)

Approach for sparse matrix – vector (for nb = 1 to 16)

Approach for sparse matrix – vector (for nb = 1 to 16)

Approach for sparse matrix – vector (for nb = 17 to 64)

Assign warps to four consecutive non-zero blocks in a row. Once the processing of the four blocks is complete, move to the next set of four blocks in the same row. Once the processing is complete for the whole row, aggregate the partial sum being held by four warps using the shared memory.

Approach for triangular solve (for nb <= 32)

Process columns of lower triangular matrix of D from left to right using a single warp Amount of parallelism available to the warp decreases as we move from left to right Need for broadcast a value from previous column processing (b(j1)) – use of shuffle instruction.

Approach for triangular solve (for nb <= 32)

All threads

• f the warp

execute this line

First prefetch

• utside the

loop Prefetches inside the loop

Prefetching

Approach for triangular solve (for nb > 32)

Performance Results

Platform and Testbed

• We ran all our experiments on a Sandy Bridge Host Processors with a

NVIDIA K40 GPU (NASA Pleiades).

• The codes were written in CUDA Fortran, compiled and executed using PGI

Fortran Compiler 15.10, and CUDA Toolkit 7.5.

• Performance testing was done on matrix sparsity structures that occurs
• ften in FUN3D application. Using a sparsity structure we generated

matrices with different block sizes. The blocks were populated using random values. To keep the memory requirement manageable, we only work with a small matrix for large block sizes.

Future Work and Conclusion

• It is possible to get benefit out of emerging architecture but requires

restructuring and designing new algorithms to map to the underlying architecture.

• Explore new sparse storage tuned for GPUs
• Explore NVIDIA PASCAL
• Device memory bandwidth: 720 GB/s