# LongBench v2 / 66ee4287821e116aacb21258

task_id: 03cff3a8-48a6-51fd-8e5b-96519acd6dee
task_key: train--66ee4287821e116aacb21258
task_revision_id: 1

{"choice_A":"They use the same techniques, such as the scaling and solving methods, to apply FP16 in multigrid preconditioner. However, their evaluations were on different architectures with different implementations. The performance speedups of the work on CPU were more prominent than the work on GPU. But the conclusions were universal and not architecture-specific.","choice_B":"Their scaling methods were different. The article \"Three-precision algebraic multigrid on GPUs\" showed the design and implementation of unstructured AMG for GPU. The guidelines and algorithms proposed in the other article were also applicable in unstructured scenarios and on GPU, but the evaluations were based on structured-specific multigrid on CPU.","choice_C":"They are based on different types of multigrids. The article \"Three-precision algebraic multigrid on GPUs\" focused on unstructured AMG that is suitable for general problems. Still, it did not discuss details about its implementations, such as where the scaling is located in the whole multigrid procedure. The other article provided a comprehensive investigation including guidelines, algorithms, and implementations, but those discussions were restricted to structured-specific multigrid.","choice_D":"They were on different hardware platforms, and had different software designs for multigrid algorithms. The article \"Three-precision algebraic multigrid on GPUs\" focused on unstructured multigrid on GPU, while the other article only showed algorithms, experiments and discussions on structured-specific multigrid on CPU. The conclusions on CPU were difficult to port to GPU.","context":"FP16 Acceleration in Structured Multigrid Preconditioner for\nABSTRACT\nHalf-precision hardware support is now almost ubiquitous. In con-\ntrast to its active use in AI, half-precision is less commonly em-\nployed in scientific and engineering computing. The valuable propo-\nsition of accelerating scientific computing applications using half-\nprecision prompted this study. Focusing on solving sparse linear\nsystems in scientific computing, we explore the technique of uti-\nlizing FP16 in multigrid preconditioners. Based on observations of\nsparse matrix formats, numerical features of scientific applications,\nand the performance characteristics of multigrid, this study formu-\nlates four guidelines for FP16 utilization in multigrid. The proposed\nalgorithm demonstrates how to avoid FP16 overflow through scal-\ning. A setup-then-scale strategy prevents FP16’s limited accuracy\nand narrow range from interfering with the multigrid’s numerical\nproperties. Another strategy, recover-and-rescale on the fly, re-\nduces the memory footprint of hotspot kernels. The extra precision-\nconversion overhead in mix-precision kernels is addressed by the\ntransformation of storage formats and SIMD implementation. Two\nablation experiments validate the effectiveness of our algorithm\nand parallel kernel implementation on ARM and X86 architectures.\nWe further evaluate three idealized and five real-world problems to\ndemonstrate the advantage of utilizing FP16 in a multigrid precon-\nditioner. The average speedups are approximately 2.75x and 1.95x\nin preconditioner and end-to-end workflow, respectively.\nCCS CONCEPTS\n• Computing methodologies →Parallel algorithms; • Mathe-\nmatics of computing →Solvers.\nThis work is licensed under a Creative Commons Attribution International\n4.0 License.\nICPP ’24, August 12–15, 2024, Gotland, Sweden\n© 2024 Copyright held by the owner/author(s).\nACM ISBN 979-8-4007-1793-2/24/08\nhttps://doi.org/10.1145/3673038.3673040\nKEYWORDS\nmultigrid, sparse matrix, structured grid, preconditioner\nACM Reference Format:\nYi Zong, Peinan Yu, Haopeng Huang, and Wei Xue. 2024. FP16 Acceleration\nin Structured Multigrid Preconditioner for Real-World Applications. In\nThe 53rd International Conference on Parallel Processing (ICPP ’24), August\n12–15, 2024, Gotland, Sweden. ACM, New York, NY, USA, 11 pages. https:\n//doi.org/10.1145/3673038.3673040\n1\nINTRODUCTION\nIn recent years, there has been a significant increase in hardware\nsupport for half-precision floating-point arithmetic [11]. This has\nbeen primarily driven by the growing demand for more efficient\nand high-performance computing. Having witnessed the pivotal\nshift towards lower precision in the fields of artificial intelligence\n(AI) and machine learning (ML) [20], it is of particular interest\nto leverage the benefits of half-precision to achieve significant\nperformance gains in scientific and industrial applications.\nThese applications, such as numerical weather prediction [21],\nocean and groundwater modeling [12], radiation hydrodynamics\nsimulation [34], and solid mechanics [14], primarily rely on numeri-\ncal simulations, especially numerical solutions to partial differential\nequation (PDE) problems. The essential part of simulations is solv-\ning a linear system 𝐴𝑥= 𝑏, which usually imposes stricter accuracy\nrequirement than AI, and are more sensitive to precision. Most\nscientific applications currently use FP64 in the solver, in which\na lower-precision preconditioner is a common accelerating tech-\nnique [1, 10]. Multigrid (MG) is a method of optimal computational\ncomplexity 𝑂(𝑁) in solving large-scale sparse linear systems [32],\nand thus widely used as the preconditioner [28].\nThe opportunity to accelerate multigrid by half-precision lies in\nits multi-level framework. Its hierarchical grids could be regarded as\na series of iterative refinement processes [19], while half-precision\niterative refinement itself has been validated in highly ideal exper-\nimental conditions, such as HPL-MxP [3]. Unlike other one-level\npreconditioners, such as incomplete lower-upper (ILU) factorization\nand symmetric Gauss-Seidel (SymGS), multigrid preconditioners\nconsume more time in solving linear systems. Its dominant role\nprovides a higher upper bound of end-to-end (E2E, i.e., the entire\n52\n\n\nICPP ’24, August 12–15, 2024, Gotland, Sweden\nYi Zong, Peinan Yu, Haopeng Huang, and Wei Xue\nprocess of solving the linear system) speedups by half-precision\naccording to Amdahl’s Law, as will be seen in our experiments.\nHowever, real-world applications are complicated. The numer-\nical distributions of nonzero entries in six matrices discretized\nfrom real-world problems are shown in Figure 1. There are large\nspans of nonzero magnitudes, and wide gaps among them. The\nranges are already beyond the compass of FP16 (IEEE 754 Stan-\ndard). An even worse problem is that accuracy may be lost near the\nupper bound of FP16, introducing further convergence issues. Brain-\nFloat16 (BF16) [13] is another half-precision format with the same\nrange as FP32. However, with fewer bits to represent significands,\nits worse accuracy is worrying to use in scientific computations.\nMore complicated challenges arise also from the multi-level\nframework. Matrices and vectors on different levels could be as-\nsigned as different precisions. There are 9𝑛(three kinds of precision\nfor matrices multiplied by three kinds for vectors give 9 combina-\ntions on a specific level) possible combinations of mix-precision in a\n𝑛-level multigrid, which makes searching for the best combination\ndifficult. Knowing how to find a good-performance combination in\ntheory is of high significance. Meanwhile, the numerical hazards\n(i.e., overflow and underflow) and accuracy loss of half-precision\nshould be treated more carefully since the connection between\ntwo consecutive levels in algebraic multigrids is determined by\na triple-matrix product. Predicting the error propagation along a\nchain of triple-matrix products is difficult. In addition to the above\nalgorithmic issues, implementation matters as the extra precision-\nconversion overhead may slow down half-precision speed. The\nmix-precision kernels require instruction-level optimization even if\nthe numerical stability is guaranteed. These issues will be addressed\nin the following from a top-down perspective.\n1e-18\n1e-15\n1e-12\n1e-09\n1e-06\n0.001\n1    \n1000 \n1e+06\n1e+09\n10-1\n100\n101\n102\npercent / %\nrhd\noil\nweather\nIEEE 754 fp16 range\nrhd-3T\noil-4C\nsolid-3D\nFigure 1: Numerical distributions of nonzero entries in six\nreal-world matrices. Detailed descriptions refer to section 6.1.\nIn this article, we have only investigated the performance of\nFP16 defined in the IEEE 754 standard. We have opted for FP16\nbecause it provides higher accuracy and is better suited for scientific\ncomputing applications, and has more extensive instruction set\nsupport than BF16 on X86 and ARM processors. Specifically, this\npaper makes the following contributions.\n• We analyze the potential gains and risks of applying FP16\nin multigrids, and summarize four guidelines. The good-\nperformance combinations of precisions are presented.\n• We propose an algorithm to adapt multigrid to the use of\nFP16. The two strategies, setup-then-scale and recover-and-\nrescale on the fly, avoid the negative influence of FP16’s\nlimited range and accuracy loss.\n• We present a high-performance parallel implementation of\nmix-precision kernels to hide the precision-conversion in-\nstruction overhead and achieve E2E speedup for ARM and\nX86 architectures.\n• We evaluate three ideal problems and five real-world prob-\nlems from four application fields. The average speedups of\nMG preconditioner reach 2.7x and 2.8x on ARM and X86,\nrespectively, while the average speedups of the entire work-\nflow are 1.9x and 2.0x.\n2\nRELATED WORK\nThe background of multigrid is first revisited. An overview of multi-\ngrid is shown in Figure 2. Based on hierarchical grids, multigrid\neliminates high-frequency errors on the finer grid by smoothers, cal-\nculates and restricts residuals to the next coarser grid, and computes\ncoarse-grid corrections recursively. Depending on how coarser\ngrids are generated, there are two types of multigrid methods. Geo-\nmetric multigrid (GMG) rediscretizes the PDE of interest in coarser\nresolutions [32]. Because rediscretization requires users’ explicit\ninvolvement with application-specific knowledge, GMG is unavail-\nable in mainstream libraries. A more black-box and user-friendly\nsolution is algebraic multigrid (AMG), which only utilizes infor-\nmation from the assembled finest-level matrix 𝐴ℎto generate all\ncoarser grids automatically.\nSmoother 𝑨𝟐𝒉𝒖𝟐𝒉= 𝒇𝟐𝒉\nResidual 𝒓𝟐𝒉= 𝒇𝟐𝒉−𝑨𝟐𝒉𝒖𝟐𝒉\nSmoother 𝑨𝟐𝒉𝒖𝟐𝒉= 𝒇𝟐𝒉\nError correction 𝒖𝟐𝒉+= 𝒆𝟐𝒉 \nSmoother 𝑨𝟒𝒉𝒖𝟒𝒉= 𝒇𝟒𝒉\nResidual 𝒓𝟒𝒉= 𝒇𝟒𝒉−𝑨𝟒𝒉𝒖𝟒𝒉\nSmoother 𝑨𝒉𝒖𝒉= 𝒇𝒉\nResidual 𝒓𝒉= 𝒇𝒉−𝑨𝒉𝒖𝒉\nSmoother 𝑨𝟒𝒉𝒖𝟒𝒉= 𝒇𝟒𝒉\nError correction 𝒖𝟒𝒉+= 𝒆𝟒𝒉 \nSmoother 𝑨𝒉𝒖𝒉= 𝒇𝒉\nError correction 𝒖𝒉+= 𝒆𝒉 \nRestrict 𝒇𝟖𝒉= 𝑹𝟒𝒉𝒓𝟒𝒉\nRestrict 𝒇𝟒𝒉= 𝑹𝟐𝒉𝒓𝟐𝒉\nRestrict 𝒇𝟐𝒉= 𝑹𝒉𝒓𝒉\nInterp 𝒆𝟒𝒉= 𝑷𝟖𝒉𝒖𝟖𝒉\nInterp 𝒆𝟐𝒉= 𝑷𝟒𝒉𝒖𝟒𝒉\nInterp 𝒆𝒉= 𝑷𝟐𝒉𝒖𝟐𝒉\n𝜴𝒉: given 𝑨𝒉\n𝜴𝟐𝒉: 𝑨𝟐𝒉= 𝑹𝒉𝑨𝒉𝑷𝟐𝒉\nComponent\nCoarsening\nSmoother\nRestrict & Interpolate\nResidual\nError Correction\nKernel Type\nSpGEMM\nSpMV/SpTRSV\nStencil/SpMV\nSpMV\nAXPY\n𝜴𝟒𝒉: 𝑨𝟒𝒉= 𝑹𝟐𝒉𝑨𝟐𝒉𝑷𝟒𝒉\nFigure 2: Multigrid overview. V-Cycle is in the solve phase.\nMultigrid is set up once and solved iteratively. Thus, the total\ntime 𝑇tot of employing a multigrid in iterative solvers is\n𝑇tot = 𝑇setup + #iter · 𝑇single\n(1)\nwhere 𝑇setup is setup time, 𝑇single is single-iteration solve time, and\n#iter is number of iterations.\nThe setup phase is denoted by green arrows in Figure 2. Multi-\ngrid constructs a series of coarser grids Ω2ℎ, Ω4ℎ, ..., based on the\nfinest-level grid Ωℎ, where the superscripts denote grid spacings of\nlevels. The essential process of setup is the triple-matrix-product in\ncomputing the coarser operator (i.e., matrix) 𝐴2𝑘ℎbased on 𝐴2𝑘−1ℎ\nof the finer level.\nIn the solve phase, denoted by blue arrows, multigrid executes\na V-Cycle [32], which starts from the finest level, traverses to the\n53\n\n\nFP16 Acceleration in Structured Multigrid Preconditioner for Real-World Applications\nICPP ’24, August 12–15, 2024, Gotland, Sweden\ncoarsest, and then reverses. On each level, multigrid invokes a\nsmoother to solve 𝐴∗ℎ𝑢∗ℎ= 𝑓∗ℎapproximately, computes residual\n𝑟∗ℎin the downward pass, and computes error correction to update\nsolution 𝑢∗ℎ, and then invokes the smoother again in the upward\npass. Restrictions of residuals 𝑟∗ℎand interpolations of errors 𝑒∗ℎ\noccur in downward and upward pass, respectively.\nHypre [15], MueLu [26], and AmgX [22] are popular AMG li-\nbraries, but still lack mix-precision support. The practice of mix-\nprecision multigrid is ahead of the theory. Precisions lower than\nFP64 reduce the memory volumes to accelerate sparse linear solvers\nthat are typically memory-bounded. As listed in Table 1, most of the\nprevious works mainly focused on using FP32, simply by changing\nthe precision of the multigrid preconditioner from FP64 to FP32 in\ntheir workflows. GMG is likely to obtain more promising end-to-\nend speedups because its matrix-free format could be compressed\ninto 50% of its original memory volume of FP64 by FP32, and 25% by\nFP16. Meanwhile, AMG for unstructured grid problems needs extra\ninteger indices to store the matrix, and the indices usually cannot\nbe compressed. Only the part of floating point data of unstructured\nAMG could leverage the benefits of lower-precision. Their practical\nexperience showed that using FP32 is safe and efficient for most\nscenarios since the number of iterations (#iter) of mix-FP32/FP64\nin their problems kept close to the original full-FP64 results.\nTable 1: Summary of mix-precision multigrid preconditioner.\n’Scale?\" means whether its method includes scaling to avoid\noverflow, and ’N/N’ for not needed. ’P.C.’ stands for precondi-\ntioner.\nRef.\nType\nScale?\nP.C. precision\nP.C. Speedup\nE2E Speedup\n[9]\nGMG\nN/N\nFP32\n~2.0x\n~1.7x\n[5]\nAMG\nN/N\nFP32\n1.1x~1.5x\nunclear\n[27]\nAMG\nN/N\nFP32\nunclear\n1.19x\n[8]\nGMG\nN/N\nFP32\n1.9x\n1.6x\n[35]\nGMG\nN/N\nFP32\n2.0x\n1.18x\n[33]\nAMG\nYes\nFP16/FP32\nunclear\n1.05x~1.35x\nOurs\nAMG\nYes\nFP16/FP32\n2.75x\n1.95x\nThe explorations of utilizing half-precision in multigrids are less\nactive than FP32, probably due to its unsafe behaviors of underflow,\noverflow, or accuracy loss. A model problem, isotropic Poisson with\nhomogeneous Dirichlet boundary condition, was tested with FP16\nin [24], where both storage and computation were in FP16. How-\never, its methods did not include out-of-range treatment necessary\nfor real-world matrices. A more recent work developed a three-\nprecision algebraic multigrids [33] as part of the Ginkgo library,\nwhich enabled arbitrary combinations of precisions (FP64, FP32,\nand FP16) on different levels. Eight problems from SuiteSparse ma-\ntrix collection [4] were tested, and E2E speedups of V-cycle mainly\nwere less than 1.2x.\nStructMG [40] is an algebraic multigrid preconditioner for struc-\ntured grid problems and has faster speed and better scalability than\nthe SOTA library hypre. The experiments in this article will be\nbased on StructMG, because other mainstream multigrid libraries\nare incapable of mix-precision configuration. Our guidelines and\nalgorithms do NOT make assumptions about the background\nproblems and can be applied to other libraries. Furthermore,\nour implementation is applicable to hypre’s SMG, PFMG,\nand SysPFMG that are designed for structure-grid problems.\nIt is worth mentioning that our guidelines and algorithms\nare also applicable to unstructured multigrid. But our subse-\nquent analysis will reveal that unstructured multigrid is not\nsuitable for harnessing the mix-precision advantages.\n3\nFP16 UTILIZATION GUIDELINES\nLower-precision accelerates each single iteration of multigrid at the\ncost of poorer convergence (i.e., faster 𝑇single but with more #iter\nin Equation 1). The previous works have shown that the full-FP32\nmultigrid as a preconditioner in FP64 iterative methods is safe and\nefficient since the #iter is kept close to its full-FP64 counterpart.\nSituations become more complicated with the more range-limited\nFP16. The basic rule of optimization is to make a balance between\nthe decreasing 𝑇single and the probably increasing #iter.\n3.1\nEagerly convert matrices to FP16.\nFor performance concern, in solving the linear system 𝐴𝑥= 𝑏,\nmatrix 𝐴is the hotspot of FP16-compressed storage, because its\nmemory access volume is often as several times large as the vectors\n𝑏and 𝑥. Given a 𝑚× 𝑚matrix 𝐴, the two vectors 𝑏and 𝑥require\n2𝑚elements, and the percentage of the matrix is\npercent𝐴=\n𝑛𝑛𝑧(𝐴)\n𝑛𝑛𝑧(𝐴) + 2𝑚× 100%\n(2)\nWe make a statistic covering the matrices whose 𝑚> 104 in\nSuiteSparse, and find that 85% of the matrices have percent𝐴> 0.7,\nand 60% of the matrices have percent𝐴> 0.8. If the integer ar-\nray indices of CSR (compressed sparse rows) format are consid-\nered, percent𝐴would be more dominant. For structured matrices,\npercent𝐴are 0.78, 0.88, and 0.90 when their nonzero patterns are\n3d7, 3d19, and 3d27, respectively.\nThe above observations still hold in the context of multigrid,\nand furthermore, the coarser-level matrices 𝐴2𝑘ℎ(𝑘≥1) would\nhave more significant percents than that of the finest-level matrix\n𝐴ℎ. This is because the Galerkin coarsening makes the nonzero\npatterns of coarser-level matrices expand [6]. Thus, the top priority\nand most beneficial task is reducing the lower bounds of memory\naccess volumes of matrices by FP16.\n3.2\nPrefer structured matrix formats that do not\nuse per-element index arrays.\nSparse linear solvers are typically memory-bounded programs. The\nupper bound of lower-precision speedup could be estimated based\non the reduction of memory access volume. Since only the floating-\npoint data of an algebraic multigrid could be compressed, struc-\ntured grid problems are more favored by FP16. In structured-grid-\nspecific multigrids, such as StructMG [40] and hypre’s SMG, PFMG,\nSysPFMG [16], the matrices are stored in structured-grid-diagonal\n(SG-DIA) format [18] to avoid extra integer indices arrays that\nare necessary for compressed sparse formats like CSR, CSC (com-\npressed sparse columns), and COO (coordinate format) to represent\nunstructured grids. We refer to a matrix as a structured matrix if\nits nonzero pattern fits the SG-DIA format, which happens when\nthe PDE is discretized on a structured grid.\n54\n\n\nICPP ’24, August 12–15, 2024, Gotland, Sweden\nYi Zong, Peinan Yu, Haopeng Huang, and Wei Xue\nThe upper bounds of lower-precision speedup in the precondi-\ntioner of SG-DIA and CSR are listed in Table 2. For unstructured\ngrid problems with CSR format, each nonzero entry requires a\nfloating-point number and an integer index in storage. Additionally,\nthe row pointer (also referred to as row map), an integer array in the\nlength of the number of rows plus 1, is amortized to each nonzero\nentry by 𝛿= (𝑚+ 1)/𝑛𝑛𝑧(𝐴) in Table 2. Their upper bounds are re-\nmarkably lower than structured grid problems with SG-DIA format,\nespecially when the problem size is too large that int64 is necessary\nfor indices. The E2E speedup would further decrease, which accords\nwith the previous experiments reported in Table 1. Without loss of\ngenerality, unstructured grid problems with CSC or COO format\ncan get results similar to CSR ones.\nTable 2: Estimated upper bound of speedup based on mini-\nmal memory access volume of matrix. ’P.C.’ stands for pre-\nconditioner. 𝛿= 15% in average for 2216 square matrices in\nSuiteSparse.\nBytes per Nonzero\nUpper Bound of P.C. Speedup\nFP64\nFP32\nFP16\nFP64/FP32\nFP32/FP16\nFP64/FP16\nSG-\n8\n4\n2\n2\n2\n4\nDIA\nCSR\n12+4𝛿\n8+4𝛿\n6+4𝛿\n12+4𝛿\n8+4𝛿<1.5\n8+4𝛿\n6+4𝛿<1.3\n12+4𝛿\n6+4𝛿<2\nint32\nCSR\n16+8𝛿\n12+8𝛿\n10+8𝛿\n16+8𝛿\n12+8𝛿<1.3\n12+8𝛿\n10+8𝛿<1.2\n16+8𝛿\n10+8𝛿<1.6\nint64\n3.3\nUse FP16 at the finest possible level.\nMatrices and vectors on different levels may have different pre-\ncisions. The number of possible combinations of precision grows\nexponentially with the number of levels. However, we will demon-\nstrate that only very limited choices are worthwhile based on the\ndistinct feature of multigrid.\nGrid complexity 𝐶G and operator complexity 𝐶O [30] are two\nmetrics to evaluate the computational overhead of a multigrid\n𝐶G =\nÍ\n𝑙𝑛𝑙\n𝑛0\nand\n𝐶O =\nÍ\n𝑙𝑍𝑙\n𝑍0\n(3)\nwhere 𝑛𝑙and 𝑍𝑙denote the number of unknowns and nonzero\nentries, respectively, on level 𝑙. 𝑙= 0 corresponds to the finest level.\nWe have evaluated 10 example problems from MFEM [31] that\nare suitable for multigrid on 11 kinds of input meshes 1 and collected\nthe 𝐶G and 𝐶O in these 60 cases (one problem on one mesh is one\ncase). These representative cases cover most of the applications\nusing MG as a preconditioner. The statistics are shown in Figure 3.\nAn obvious situation is that most cases have low complexities. The\ncumulative frequency curves verify that 𝐶G and 𝐶O are lower than\n1.2 and 1.5, respectively, in 80% of the cases, and lower than 1.15\nand 1.22, respectively, in 60% of the cases.\n1Specifically, ex1p(Laplace), ex4p(Grad-div), ex5p(Darcy), ex6p(Laplace with AMR),\nex8p(DPG for Laplace), ex14p(DG Diffusion), ex26p(Multigrid Preconditioner) are\ntested with meshes including star.mesh, beam-tri.mesh, escher-p2.mesh, mobius-\nstrip.mesh, square-disc-p3.mesh, and rt-2d-q3.mesh. Ex2p(Linear Elasticity), ex17p(DG\nLinear Elasticity) and ex21p(AMR for linear elasticity) are tested with meshes including\nbeam-tri.mesh, beam-quad.mesh, beam-tet.mesh, beam-hex.mesh, beam-wedge.mesh,\nand beam-quad-nurbs.mesh. Refer to https://mfem.org/examples for more detailed\ndescription.\nTheir low complexities result from aggressive coarsening [36]\nin AMG, and high-dimensional coarsening in StructMG. Lower\n𝐶G and 𝐶O benefit the overall performance of multigrid, and thus,\naggressive coarsening has been recommended in the best practice\nguide [6] by the hypre team. The outlier (𝐶G > 1.5 and 𝐶O > 2)\nin Figure 3 are the cases whose default settings do not include\naggressive coarsening.\n𝐶G\n𝐶O\nFigure 3: Complexity statistics for 60 MFEM cases.\nThis indicates that the aggregated overhead of all coarser grids\nis not as significant as half the importance of the finest grid. The\ndominant role of the finest grid could also be observed in previous\nworks on performance modeling [7]. Therefore, the finer level is\naccelerated by FP16, the more performance gains will be obtained.\nFP16 should be utilized to reduce the memory access volumes on\nlevels as fine as possible, which is in contrast to [33] that proposed\n’DP-SP-HP’(i.e., the first level’s matrix uses double precision, the\nsecond level’s matrix uses single precision, and the other levels’\nmatrices use half precision). Despite the massive number of possible\ncombinations of precisions in multigrid, only those who apply FP16\non the finer levels are worthwhile from a performance perspective.\nBased on this observation, multigrid should allow switching the\nprecision from FP16 to higher precisions since a specific level.\n3.4\nAvoid vectors in FP16.\nFor scientific safety concerns, vectors are not suitable to be stored\nor computed in FP16 because overflow could happen out of control.\nIn solving a linear system 𝐴𝑥= 𝑏, the given matrix 𝐴is static while\nthe unknown vector 𝑥changes dynamically. It is difficult to predict\nwhich element of 𝑥may overflow sometime. However, the cost of\noverflow is intolerable. Even if a single element overflows to ’inf’,\nit immediately propagates to ’NaN’, crashing the entire program.\nThe above four guidelines instruct the following algorithmic\ndesigns. Vectors should be kept in FP32, while FP16 compresses the\nmemory volumes of matrices as fine-level as possible. Aiming for\nbetter leverage of FP16, structured grids should be considered to\ndiscretize the PDE of interest.\n4\nALGORITHMS\nThis section involves three kinds of precision. Variables of different\nprecision in Algorithm 1, 2, 3 are marked in different colors. The\n:::::::::::::::\niterative precision marked in red with waved underline is the\ncomputation and storage precision of iterative solvers, usually as\nFP64 or FP32, which is determined by users’ applications. The\ncomputation precision of preconditioners is marked in blue\n55\n\n\nFP16 Acceleration in Structured Multigrid Preconditioner for Real-World Applications\nICPP ’24, August 12–15, 2024, Gotland, Sweden\nwith straight underline, usually as FP32. The ............................\nstorage precision of\n......................\npreconditioners is marked in green with dotted underline, usually\nas FP16.\n4.1\nSetup-then-scale\nThe strategy in the setup phase is setup-then-scale. As shown in\nLine 1-3 in Algorithms 1, a normal setup process is completed first\nin high precision, which includes Galerkin coarsening on all levels.\nAlgorithm 1: MG_setup_for_FP16\nInput: matrix :\n𝐴\nOutput: matrices ....\n𝐴0, ....\n𝐴1, ..., ....\n𝐴𝐿on hierarchical grids and\nsmoothers ...\n𝑆0, ...\n𝑆1, ..., ....\n𝑆𝐿, and 𝑄0, 𝑄1, ..., 𝑄𝐿if needed\n1 for 𝑖= 0, 1, . . . , 𝐿−1 do\n2\n𝐴𝑖+1\n:: ←𝑅𝑖\n:𝐴𝑖\n: 𝑃𝑖+1\n:: // Galerkin coarsening\n3 end\n4 for 𝑖= 0, 1, . . . , 𝐿do\n5\nif need to scale then // truncation after scaling\n6\n::\n𝑄𝑖←\n1\n𝐺𝑖extract_diagonals( :\n𝐴𝑖);\n7\n𝐴𝑖\n: ←𝑄𝑖\n:\n−1/2𝐴𝑖\n: 𝑄𝑖\n:\n−1/2;\n8\n𝐴𝑖\n... ←𝐴𝑖\n: ;\n9\n𝑄𝑖←𝑄𝑖\n: ;\n10\nelse // direct truncation\n11\n𝐴𝑖\n... ←𝐴𝑖\n:\n12\nend\n13\n𝑆𝑖\n.. ←smoother_setup(𝐴𝑖\n: ) ;\n// setup smoothers\n14 end\nThe special treatment to adapt to FP16 is the need to scale branch\nfrom Line 5 to 9. The purpose of scaling is to transform nonzero\nvalues of the matrix ::\n𝐴𝑖into FP16 range, where 𝐺𝑖is a constant\nparameter. We have the following theorem to ensure the safety of\nFP16.\nTheorem 4.1. Truncation to FP16 would not cause unsafe ’inf’ if\n𝐺𝑖is appropriately chosen.\nProof. Since scaling operations are independent across different\nlevels, level indices are ommited in the proof. Given the matrix on a\nspecific level 𝐴(𝑎𝑖𝑗) ∈R𝑁×𝑁, the diagonal matrix 𝑄in Algorithm 1\nis\n𝑄= 1\n𝐺diag(𝑎11,𝑎22, ...,𝑎𝑁𝑁)\nThe element at 𝑖-row and 𝑗-column after multiplication of\n𝑄−1/2𝐴𝑄−1/2 will be\n𝐺\n𝑎𝑖𝑗\n𝑎1/2\n𝑖𝑖𝑎1/2\n𝑗𝑗\nLet 𝑆:= FP16_MAX denote the upper bound of FP16. The abso-\nlute sign of |∗| is omitted because representation range of floating-\npoint numbers consists of a positive half and a negative half that\nare completely symmetric around zero. The following condition\nneeds to be satisfied for all 𝑖, 𝑗to avoid overflow:\n𝐺\n𝑎𝑖𝑗\n𝑎1/2\n𝑖𝑖𝑎1/2\n𝑗𝑗\n< 𝑆\nWhen 𝑎𝑖𝑗is large, it requires 𝐺to be small. Therefore, the maximal\n𝐺that avoids overflow is\n𝐺max = 𝑆max\n𝑖,𝑗\n\n\n\n\n\n\n\n𝑎1/2\n𝑖𝑖𝑎1/2\n𝑗𝑗\n𝑎𝑖𝑗\n\n\n\n\n\n\n\nAs long as a 𝐺< 𝐺max is chosen, overflow of FP16 is avoided in\nthe truncation after scaling. Since that 𝑄is stored in preconditioner\ncomputation precision (usually FP32, but not FP16 in any cases),\nand will not overflow due to a small 𝐺. Notice that the square-root\noperation requires all diagonal entries of the original matrix 𝐴\nare positive, which is naturally included by M-matrix property. It\ndoes not make new assumptions because the convergent theory of\nmultigrid itself requires 𝐴to be an M-matrix [32].\n□\nMatrices ::\n𝐴𝑖are truncated to storage precision in Line 8 and\ndiagonal matrices ::\n𝑄𝑖are truncated to computation precision in\nLine 9. If scaling is unnecessary, direct truncation in Line 11 follows\nthe standard setup process. Based on previous observations, to\navoid overflow, whether there is a need to scale or not depends on\nif there exists::\n𝐴𝑖’s values greater than FP16_MAX.\nFinally, high-precision matrices are used to setup correspond-\ning smoothers. Data in smoothers, such as the factorized lower\nand upper triangular matrices ˜\n𝐿, ˜\n𝑈in ILU, are calculated in itera-\ntive precision followed by truncation to storage precision. Notice\nthat high-precision :\n𝐴𝑖and ::\n𝑄𝑖will no longer be used in multigrid\npreconditioner after the setup phase.\n4.2\nRecover-and-rescale on the fly during solve\nAll the optimizations in this article focus on preconditioners, so\nnothing special is applied to iterative solvers. Due to space lim-\nitations, we use a stationary iteration to illustrate how an FP16-\naccelerated preconditioner is employed. Other iterative algorithms,\nsuch as Conjugated Gradient (CG) and Generalized Minimal Resid-\nual (GMRES), have the same way of invoking the preconditioner.\nAs shown in Line 3 in Algorithm 2, residual is calculated by\nhigh-precision matrices and vectors in iterative solvers. The pre-\nconditioner accepts a lower-precision residual vector as input and\nreturns an error vector of the same precision. Explicit precision\ntransitions are in Line 4 and 6.\nAlgorithm 2: Stationary iterative method.\nInput: matrix :\n𝐴, right-hand-side :\n𝑏, initial solution :\n𝑥\nOutput: approximated solution :\n𝑥\n1 Initialize preconditioner: MG_setup_for_FP16(\n:\n𝐴);\n2 while not converged do\n3\n:\n𝑟←:\n𝑏- ::\n𝐴𝑥;\n4\n𝑟←:\n𝑟;\n// truncate residual\n5\n𝑒←MG_solve_with_FP16(𝑟) ;\n// apply multigrid\n6\n:\n𝑒←𝑒;\n// recover error\n7\n:\n𝑥←:\n𝑥+ :\n𝑒;\n8 end\nThere is nothing in iterative precision throughout the V-Cycle\nof multigrid, as shown in Algorithm 3. The basic idea is that matri-\nces’ nonzero entries are stored in FP16 and thus need recovering,\n56\n\n\nICPP ’24, August 12–15, 2024, Gotland, Sweden\nYi Zong, Peinan Yu, Haopeng Huang, and Wei Xue\nwhile vectors’ values throughout the solve phase remain unchanged\nas FP32. Therefore, matrices must recover and rescale whenever\nsparse matrix-vector product (SpMV) and sparse triangular solve\n(SpTRSV) are performed. The way of rescaling is shown in Line 7\nof Algorithm 3. Since computation performs in FP32, the diagonal\nmatrices 𝑄𝑖are used to rescale values of ...\n𝐴𝑖.\nNote that the recover-and-rescale is on the fly. Matrices data of\nFP16 is retrieved, transformed into FP32, and used to calculate only\nwhen needed. Data of FP32 is not explicitly maintained; otherwise,\nthe advantage of accessing fewer volumes of memory no longer\nexists. Rescaling in smoother_solve at Line 4 and 17 is similar if\n..\n𝑆𝑖was setup using a scaled ::\n𝐴𝑖previously. SymGS and ILU could\nbe configured as the smoothers of different levels by users, and a\ndirect solver for the coarsest level.\nAlgorithm 3: MG_solve_with_FP16\nInput: right-hand-size 𝑏\nOutput: approximated solution 𝑥\n1 Initialize: 𝑓0 ←𝑏;\n2 for 𝑖= 0, 1, ..., 𝐿−1 do // forward part of V-Cycle\n3\nfor 𝑗= 0, 1, ..., 𝜈1 −1 do // pre-smoothing 𝜈1 times\n4\n𝑢𝑖←smoother_solve( ..\n𝑆𝑖, 𝑓𝑖, 𝑢𝑖);\n5\nend\n6\nif scaled in setup then // compute residual\n7\n𝑟𝑖←𝑓𝑖−𝑄𝑖1/2 ...\n𝐴𝑖𝑄𝑖1/2 𝑢𝑖;\n8\nelse\n9\n𝑟𝑖←𝑓𝑖- ...\n𝐴𝑖𝑢𝑖;\n10\nend\n11\nif 𝑖< 𝐿−1 then // restrict residual to next level\n12\n𝑓𝑖+1 ←Restrict(𝑟𝑖);\n13\nend\n14 end\n15 for 𝑖= 𝐿−1, ..., 1, 0 do // backward parf of V-Cycle\n16\nfor 𝑗= 0, 1, ..., 𝜈2 −1 do // post-smoothing 𝜈2 times\n17\n𝑢𝑖←smoother_solve( ....\n𝑆𝑖𝑇, 𝑓𝑖, 𝑢𝑖) ;\n18\nend\n19\nif 𝑖> 0 then // interpolate error to next finer grid\n20\n𝑢𝑖−1 ←𝑢𝑖−1 + Interpolate(𝑢𝑖);\n21\nend\n22 end\n23 𝑥←𝑢0 ;\n4.3\nPractical Remarks\nFurther improvement of replacing vector with FP16 is limited for\nperformance concerns. In contrast, overflow risks in vectors’ val-\nues are potential for scientific accuracy concern, as discussed in\nSection 3.1 and 3.4.\nOur method introduces the additional overhead of accessing 𝑄𝑖\nin the solve phase. A natural counterpart would be \"scale-then-\nsetup\", which first scales the problem matrix :\n𝐴, then passes it to\nmultigrid to setup matrices ::\n𝐴𝑖on all levels, and finally truncates\nthem to ...\n𝐴𝑖. Though the scale-then-setup strategy is free of 𝑄𝑖inside\nV-Cycle, our setup-then-scale in Section 4.1 has two-fold advan-\ntages. First, our setup-then-scale does not require users’ explicit\ninvolvement to scale the original problem matrix and keeps multi-\ngrid as a \"black-box\" tool. Most importantly, it does not interfere\nwith the original setup phase, especially the chain of triple-matrix\nproducts. Otherwise, scale-then-setup negatively influences the\nchain because only floating-point numbers in the FP16 range could\nbe used even if triple-matrix products are calculated with FP64,\nand even worse, may still incur overflow or underflow. In practice\nof our various problems, the overhead of diagonal matrices 𝑄𝑖is\ncost-efficient compared to the increase of #iter, because 𝑄𝑖only\noccupies the same amount of memory as a vector. The advantage\nof setup-then-scale over scale-then-setup will be verified in the\nablation experiments in Section 7.1.\nUnlike overflow, which inevitably crashes the program by pro-\nducing NaN values, underflow is not as catastrophic. However,\nit can sometimes slow down convergence speed and, in extreme\ncases, lead to non-convergence. Underflow occurs due to the\ntriple-matrix products during the setup phase, where multipli-\ncation causes the absolute values of matrices at coarser levels\nto become increasingly smaller. A simple but effective method\nto avoid the adverse effects of underflow is to switch back to\nhigher precision such as FP32. Specifically, the latest version of\nStructMG [40] introduces a tunable parameter called shift_levid.\nFrom the shift_levid level to the coarsest level, matrices are stored\nin the computation precision of preconditioners rather than in the\n......................\nstorage precision. Switching to higher precision at coarser levels\ndoes not introduce significant overhead, because the finer levels\nare the primary hotspots in multigrid methods, as analysed in Sec-\ntion 3.3.\n5\nINSTRUCTION-LEVEL OPTIMIZATION\nBased on the FP16 utilization guidelines and the algorithms, the data\nof matrices and vectors in multigrid are stored in FP16 and FP32,\nrespectively. SpMV and SpTRSV kernels become mix-precision,\nwhich introduces the additional overhead of precision conversion.\nA representative profiling [37] indicated that SymGS smoother (a\nspecialized form of SpTRSV) accounts for 78% of the entire execu-\ntion time of the HPCG benchmark, while SPMV has the second-\nlargest contribution of 20%. Therefore, the overhead could not be\nneglected, as will be shown by the results of Section 7.2, and requires\ninstruction-level optimization.\n5.1\nVectorization: From AOS to SOA\nThe implementation is based on structured-grid-specific multigrids,\nsuch as StructMG, and hypre’s SMG, PFMG, and SysPFMG. The\nmatrices are stored in array-of-structure (AOS) format in these\nmultigrids, as shown in Figure 4. The nonzero entries with the\nsame subscripts (e.g., (𝑖, 𝑗)) are stored contiguously in the mem-\nory. The superscripts, 0~3, of the matrix entries indicate which\nneighboring position corresponds to this matrix entry. For this\nexample, 0 corresponds to (𝑖−1, 𝑗−1), and 1, 2, 3 correspond\nto (𝑖−1, 𝑗), (𝑖−1, 𝑗+ 1), (𝑖, 𝑗−1), respectively. Alternatively, the\nstructure-of-array (SOA) format shuffles the entries and stores those\nwith the same superscripts contiguously.\nIn full-FP32 situations, AOS format is usually good enough to\nfully utilize the memory bandwidth, which requires only one load\n(ldr) instruction for each 4-byte entry to prepare for multiplica-\ntion. However, when the matrix is stored in FP16, AOS requires\n57\n\n\nFP16 Acceleration in Structured Multigrid Preconditioner for Real-World Applications\nICPP ’24, August 12–15, 2024, Gotland, Sweden\nAOS (in FP16)\nSOA (in FP16)\nScalar instruction for AOS\n𝑨𝑖,𝑗\n0 ←fcvt ldr 𝑨𝑖,𝑗\n0\n𝑨𝑖,𝑗\n1 ←fcvt ldr 𝑨𝑖,𝑗\n1\n𝑨𝑖,𝑗\n2 ←fcvt ldr 𝑨𝑖,𝑗\n2\n𝑨𝑖,𝑗\n3 ←fcvt ldr 𝑨𝑖,𝑗\n3\n𝒚𝑖𝑗←𝑨𝑖,𝑗\n0 ∗𝒙𝑖−1 𝑗−1\n+𝑨𝑖,𝑗\n1 ∗𝒙𝑖−1 𝑗\n+𝑨𝑖,𝑗\n2 ∗𝒙𝑖−1 𝑗+ 1\n+𝑨𝑖,𝑗\n3 ∗𝒙𝑖\n𝑗−1\nSIMD for SOA\n𝑨𝑖,𝑗:(𝑗+3)\n0\n←fcvt ldr 𝑨𝑖,𝑗:(𝑗+3)\n0\n𝑨𝑖,𝑗:(𝑗+3)\n1\n←fcvt ldr 𝑨𝑖,𝑗:(𝑗+3)\n1\n𝑨𝑖,𝑗:(𝑗+3)\n2\n←fcvt ldr 𝑨𝑖,𝑗:(𝑗+3)\n2\n𝑨𝑖,𝑗:(𝑗+3)\n3\n←fcvt ldr 𝑨𝑖,𝑗:(𝑗+3)\n3\n𝒚𝑖𝑗: (𝑗+ 3) ←𝑨𝑖,𝑗:(𝑗+3)\n0\n∗𝒙𝑖−1\n𝑗−1 : (𝑗+ 2)\n+𝑨𝑖,𝑗:(𝑗+3)\n1\n∗𝒙𝑖−1\n𝑗\n∶(𝑗+ 3)\n+𝑨𝑖,𝑗:(𝑗+3)\n2\n∗𝒙𝑖−1\n𝑗+ 1 : (𝑗+ 4)\n+𝑨𝑖,𝑗:(𝑗+3)\n3\n∗𝒙𝑖\n𝑗−1 : (𝑗+ 2)\n𝑨𝑖,𝑗\n0\n𝑨𝑖,𝑗\n1\n𝑨𝑖,𝑗\n2\n𝑨𝑖,𝑗\n3\n𝑨𝑖,𝑗+1\n0\n𝑨𝑖,𝑗+1\n1\n𝑨𝑖,𝑗+1\n2\n𝑨𝑖,𝑗+1\n3\n𝑨𝑖,𝑗+2\n0\n𝑨𝑖,𝑗+2\n1\n𝑨𝑖,𝑗+2\n2\n𝑨𝑖,𝑗+2\n3\n𝑨𝑖,𝑗+3\n0\n𝑨𝑖,𝑗+3\n1\n𝑨𝑖,𝑗+3\n2\n𝑨𝑖,𝑗+3\n3\n𝑨𝑖,𝑗\n0\n𝑨𝑖,𝑗+1\n0\n𝑨𝑖,𝑗+2\n0\n𝑨𝑖,𝑗+3\n0\n𝑨𝑖,𝑗\n1\n𝑨𝑖,𝑗+1\n1\n𝑨𝑖,𝑗+2\n1\n𝑨𝑖,𝑗+3\n1\n𝑨𝑖,𝑗\n2\n𝑨𝑖,𝑗+1\n2\n𝑨𝑖,𝑗+2\n2\n𝑨𝑖,𝑗+2\n3\n𝑨𝑖,𝑗\n3\n𝑨𝑖,𝑗+1\n3\n𝑨𝑖,𝑗+2\n3\n𝑨𝑖,𝑗+3\n3\nFigure 4: Illustration of SpMV in AOS and SOA format.\nGreen and blue data are in FP16 and FP32, respectively. The\ndashed box emphasizes the precision-conversion instruc-\ntions needed to prepare for FP32 multiplication. SOA is SIMD-\nfriendly. Computations of 4 elements (𝑗: (𝑗+ 3) written in\nMatlab style) are vectorized.\none load, and additionally one floating-point convert (fcvt) instruc-\ntion for each 2-byte entry. The arithmetic intensity during data\npreparation is now as 4 times high as the full-FP32 situation. A\ndecrease of bandwidth efficiency2 will be observed for AOS in the\nmix-FP16/FP32 situation. Therefore, SOA format is necessary to\namortize the precision-conversion overhead via vectorization. As\nshown in the dashed box in Figure 4, it requires only one ldr, and\none fcvt instruction for every four 2-byte entries if the SIMD length\nis 128-bit. The vectorization of SpMV is straightforward because\ncomputations of different elements are independent. SpTRSV re-\nquires more sophisticated parallel strategy [39, 40].\n6\nEXPERIMENTAL SETUP\nThis section describes the problems, solvers settings and machines\nconfigurations.\n6.1\nProblems\nOur selection of test problems aims to cover different domains and\ncharacteristics as extensively as possible. Laplace27 is an idealized\nbenchmark problem in performance evaluation and modeling, such\nin hypre’s reports [17] and HPCG [39]. Laplace27∗108 is an ideal-\nized problem that we construct by multiplying the coefficients of\nlaplace27 by 108 to observe the influence of out-of-range of FP16.\nRhd and rhd-3T are from radiation hydrodynamics [34]. \"3T\" means\nthree temperatures (radiation, electron, and ion). Oil and oil-4C are\nfrom petroleum reservoir simulation. Settings of SPE1 and SPE10\nbenchmarks [23] are combined to generate larger cases via Open-\nCAEPoro [38]. \"4C\" means four components (oil, water, gas, and\ndissolved gas in live oil). Weather is from atmospheric dynamics,\nprovided by the dynamic core of GRAPES-MESO [25], the national\nNWP system of the Chinese Meteorological Administration [2].\nThe largest case of 637M #dof is a 2km resolution of Chinese re-\ngional forecasting in Dec 2018. Solid-3D cases discretized from the\n2The bandwidth efficiency is the measured bandwidth of the kernel over the bandwidth\nof the stream benchmark, which can be regarded as the evaluation of architecture\nefficiency of the kernel. The measured bandwidth refers to the minimal theoretical\nmemory volume to access divided by the measured kernel time.\nweak form of linear elasticity problem in solid mechanics [29], are\ngenerated by ourselves. \"3D\" means three displacements associated\nwith each element. Data of all problems are available online3.\nThese problems’ characteristics are listed in Table 3. The basic\ninformation includes PDE types4, nonzero patterns (’Pattern’ field)5,\ntotal degrees of freedom (’#dof’ field), total numbers of nonzero\nentries (’#nnz’ field).\nTheir numerical features are also considered, including whether\nthe problem is discretized from a real-world application, whether\nits range is out-of-range of FP16, and the distance from FP16 (’Dist.’\nfield) if it is out. The ranges of different problems and FP16 are\ndisplayed in Figure 1. Anisotropy (also referred to as multi-scale\nproperty [34]) is an important metric indicating that the charac-\nteristics of the linear system vary based on the direction in which\nthey are measured. The more anisotropic a linear system is, the\nmore difficult it is for the solver to converge on it, and the more\nchallenging it is for FP16 to accelerate. The qualitative statements\nare in the ’Aniso.’ field of Table 3, and more detailed visualiza-\ntion can be found in Figure 5. Laplace27 and laplace27∗108 are\nfully isotropic and have constant coefficients. Rhd-3T is highly\nanisotropic due to non-smooth coefficients and multi-physics cou-\npling in three-temperature equations. Rhd is relatively isotropic\nafter decoupling from the rhd-3T system. Oil-4C and oil are both\nhighly anisotropic due to inhomogeneous permeability. The strong\nanisotropy of the weather problem comes from irregular earth to-\npography and nonuniform latitudinal spacing. Solid-3D is relatively\nisotropic and has homogeneous coefficients. The condition numbers\nare also included in Table 3. The condition number of weather is\nevaluated based on a much smaller-size matrix of 9.95M discretized\nfrom the same problem, because the original size of 637M is too\nlarge to compute the condition number in a reasonable time.\n(i) Ranges statistics\n(ii) Multi-scale measurements statisti\nAnisotropic\nIsotropic\nFigure 5: Statistics of multi-scale metric [34] in six problems.\n6.2\nSolvers\nIterative Krylov [28] solvers, and the precision corresponding to\nSection 4 in different problems are in Table 3. Users’ applications\ndetermine the iterative precisions (in red). The detailed configu-\nration of StructMG as a preconditioner can be found online3. All\nproblems have low 𝐶G and 𝐶O, which matches our observation and\nguideline in Section 3.3.\n3https://zenodo.org/records/10023590. Data of laplace27 and laplace27∗108 are not\nincluded because they can be easily constructed by readers.\n4Scalar PDE means only one unknown is associated with each grid element, while\nvector PDE has multiple knowns.\n53d15 and 3d19 expand to 3d27 on coarser grids in StructMG and hypre’s SMG, PFMG,\nand SysPFMG.\n58\n\n\nICPP ’24, August 12–15, 2024, Gotland, Sweden\nYi Zong, Peinan Yu, Haopeng Huang, and Wei Xue\nTable 3: Test problems characteristics. ’M’ for million and ’B’ for billion in ’#dof’ and ’#nnz’ fields. More detailed statistics of\n’Dist.’ and ’Aniso.’ fields refer to Figure 1 and 5.\nProblem\nBasic Information\nNumerical Features\nSolver Information\nPDE\nPattern\n#dof\n#nnz\nReal-world?\nOut-of-FP16?\nDist.\nAniso.\nCond.\nPrecision\nSolver\n𝐶G\n𝐶O\nlaplace27\nscalar\n3d27\n16.8 M\n453 M\nNo\nNo\nNone\n3e+03\n:::\nFP64/ FP32/ ......\nFP16\nCG\n1.14\n1.14\nlaplace27∗108\nscalar\n3d27\n16.8 M\n453 M\nNo\nYes\nFar\nNone\n3e+03\n:::\nFP64/ FP32/ ......\nFP16\nCG\n1.14\n1.14\nrhd\nscalar\n3d7\n2.10 M\n14.7 M\nYes\nYes\nFar\nLow\n1e+08\n:::\nFP64/ FP32/ ......\nFP16\nCG\n1.14\n1.14\noil\nscalar\n3d7\n31.5 M\n220.2 M\nYes\nNo\nHigh\n1e+04\n:::\nFP64/ FP32/ ......\nFP16\nGMRES\n1.14\n1.14\nweather\nscalar\n3d19\n637 M\n12.1 B\nYes\nYes\nNear\nHigh\n1e+05\n:::\nFP64/ FP32/ ......\nFP16\nGMRES\n1.31\n1.44\nrhd-3T\nvector\n3d7\n6.30 M\n52.4 M\nYes\nYes\nFar\nHigh\n1e+15\n:::\nFP64/ FP32/ ......\nFP16\nCG\n1.14\n1.14\noil-4C\nvector\n3d7\n31.5 M\n880 M\nYes\nYes\nNear\nHigh\n1e+05\n:::\nFP64/ FP32/ ......\nFP16\nGMRES\n1.14\n1.14\nsolid-3D\nvector\n3d15\n11.8 M\n531 M\nNo\nYes\nFar\nLow\n1e+07\n:::\nFP64/ FP32/ ......\nFP16\nCG\n1.14\n1.26\n6.3\nMachines\nExperiments are evaluated on ARM and X86 platforms, as shown\nin Table 4. Similar to [40], the best result at a specific degree of\nparallelism is reported from tests of various MPI/OpenMP ratios\nwith load-balance process partitions. 1:1 (MPI-only), 1:2, 1:4, 1:8,\n1:16, and 1:32 are tested when each NUMA has 32 available cores.\n1:1, 1:2, 1:3, 1:5, 1:6, 1:10, 1:15, and 1:30 are tested when only 30\ncores are available because ARM’s McKernel mechanism reserves\ntwo cores for OS to reduce system noise for large-scale tests.\nTable 4: Machines Configurations.\nSystem\nARM\nX86\nProcessor\nKunpeng 920-6426\nAMD EPYC-7H12\nFrequency\n2.60 GHz\n2.60~3.30 GHz\nCores per node\n128 (64 per socket)\n128 (64 per socket)\nL1/L2/L3 per core\n64 KB/512 KB/1 MB\n32 KB/512 KB/4 MB\nStream Triad BW\n138 GB/s\n100 GB/s\nMemory per Node\n512 GB DDR4-2933\n256 GB DDR4-3200\nMax Nodes\n64\n64\nNetwork\n100Gbps InfiniBand\n100Gbps InfiniBand\nMPI/Compiler\nOpenMPI-4.1.4/gcc-9.3.0\nIntel-OneAPI-2021.6\n7\nRESULTS AND ANALYSIS\nThis section will present the results in a local-to-global perspective.\nA controlled variables experiment first verifies the algorithmic\neffect. Kernel performance is then measured to demonstrate the\nnecessity of hiding precision-conversion overhead. The above two\neffects will be combined to see the overall speedups in a single\nprocessor. Finally, strong scalability tests are presented.\n7.1\nAlgorithmic Effect\nDescending curves of residual norm display how fast a precon-\nditioned solver could converge to the desired solution, as shown\nin Figure 6. Laplace27, laplace27∗108, weather, rhd, and rhd-3T\nare analyzed in this ablation experiment, because they have dis-\ntinct numerical features in different dimensions (i.e., out-of-FP16,\ndistance, and anisotropy) that are labeled in the upper-right cor-\nner of sub-figures. A controlled variable comparison could obtain\nrepresentative conclusions.\nFive combinations of precisions and strategies are evaluated.\nFull64 is the baseline whose curves descend at the highest rate in\nall problems. K64P32D32 is the FP64 iterative solver preconditioned\nNumber of Iteration\nRelative Residual Norm: ||r||2/||b||2\n2\n4\n6\n8\n10\n12\n10-9\n10-6\n10-3\n100\n(a) laplace27\n2\n4\n6\n8\n10\n12\n10-9\n10-6\n10-3\n100\n(b) laplace27*108\n5\n10\n15\n10-9\n10-6\n10-3\n100\n(c) weather\n10\n20\n30\n40\n50\n60\n10-9\n10-6\n10-3\n100\n103\n(d) rhd\n20\n40\n60\n80\n100\n10-9\n10-6\n10-3\n100\n(e) rhd-3T\nFull64\nK64P32D32\nK64P32D16-none\nK64P32D16-scale-setup\nK64P32D16-setup-scale\nin FP16\nisotropic\nfar from FP16\nisotropic\nfar from FP16\nweakly-anisotropic\nnear FP16\nhighly-anisotropic\nfar from FP16\nhighly-anisotropic\nFigure 6: Descending curves of relative residual norm in five\nrepresentative problems. ’K’, ’P’, and ’D’ in the legends stand\nfor iterative precision, computation precision of precondi-\ntioners, and storage precision of preconditioners in Section 4,\nrespectively. The yellow curve of ’K64P32D16-none’ only oc-\ncurs in (a) because it fails to solve the other problems. The\nother 4 curves coincide in (b), and all 5 curves coincide in (a).\nThe oscillations in (c) are due to false convergence of GMRES.\nby an FP32 multigrid, commonly used in previous works. There\nare three strategies based on K64P32D16 where FP16 is used in\nstorage in the multigrid preconditioner. ’None’ means no scaling\nis applied, and thus would lead to NaN in all problems except\nlaplace27. ’Scale-setup’ and ’setup-scale’ correspond to the scale-\nthen-setup and setup-then-scale in Section 4.1, respectively. All five\ncombinations perform nearly identically in the idealized benchmark\nof laplace27, showing completely overlapping residual reduction\ncurves in Figure 6(a). All four combinations, except for ’K64P32D16-\nnone’, perform almost identically in the idealized benchmark of\n59\n\n\nFP16 Acceleration in Structured Multigrid Preconditioner for Real-World Applications\nICPP ’24, August 12–15, 2024, Gotland, Sweden\nlaplace27*∗108, as shown in Figure 6(b). The advantages of setup-\nthen-scale over scale-then-setup are demonstrated in three real-\nworld problems. Even in the weather6 problem whose values are\nnear the range of FP16, setup-then-scale shows faster-descending\ncurves and results in 11 iterations to converge to ||𝑟||2/||𝑏||2 <1e-\n10. In comparison, scale-then-setup needs 15 iterations to reach the\nsame convergent threshold. Their differences are further magnified\nand showcased in rhd and rhd-3T problems, where the values are\nfar from the range of FP16. Scale-then-setup could not converge in\nthese two problems.\n7.2\nKernel Optimization Effect\nSection 5.1 has discussed the extra overhead of precision-conversion\ninstruction, for which the SIMD-friendly SOA format should be\nused to obtain the expected speedups. The two most essential ker-\nnels involving matrices of FP16 in V-cycle are SpMV and SpTRSV,\noccupying over 80% of multigrid time [37]. We evaluated these\ntwo kernels in StructMG on matrices of different patterns. The\nbaselines are full-FP32 kernels of AOS format without precision-\nconversion overhead (denoted by ’MG-fp32/fp32’ in Figure 7). Our\nbaseline SpTRSV and SpMV implementations are ~3.5x and ~1.8x\nfaster than ARM Performance Library (ARMPL, latest version 23.10)\non ARM architecture, and ~2.2x and ~1.2x faster than MKL (version\n2022.1.0) on X86, respectively. The time of symbolic analysis has\nbeen excluded for ARMPL and MKL. More importantly, our mix-\nprecision SpMV and SpTRSV of SOA format with SIMD (denoted\nby ’MG-fp16/fp32(opt)’) could further reduce the kernel times. The\ntime reduction is proportional to the reduction of memory volumes\nneeded to access in the kernel. The higher the ratio of volumes the\nmatrix occupies, the higher the speedup would be. The theoretical\nmaximum reachable speedups (denoted by ’Max-fp16/fp32’) are\nbased on the reductions of memory volumes needed to access (i.e.,\nthe memory volume of full-FP32 divided by that of mix-precision),\nwhich represents the upper-bounds of the mix-precision perfor-\nmance. As shown in Figure 7, our optimized implementations could\nreach similar speedups with the maximum ones. On the other hand,\nas we have expected, performance degradation is observed for mix-\nprecision kernels of AOS format (denoted by ’MG-fp16/fp32(naive)’),\nwhich are straightforward extensions of the baseline counterparts.\nThe degradation is more pronounced for SpTRSV, which requires\nsophisticated parallelization. The experiments in Figure 7 utilized\nall cores within a single NUMA (32 cores on ARM and 64 cores on\nX86) with multi-threading and speedups are geometrically averaged\nover problem sizes of 2563, 2883, 3203, 3523 and 3843.\n7.3\nEnd-to-end Improvement\nOur algorithm aims to minimize the increase of #iter when FP16\nis utilized, and the efficient implementation aims to reduce 𝑇single\nin Equation (1). Combining the two sides leads to an overall high-\nperformance solver preconditioned by FP16-accelerated multigrid.\nFigure 8 shows the improvement on a single ARM processor\nafter accelerating the MG preconditioner by FP16 in the entire\nfull-FP64 workflow. All times have been normalized according to\nthe total time of the Full64 solver. The additional overhead of the\n6The iterative precision of weather is increased to FP64 in this subsection 7.1 to collect\na longer history of descending residual for analysis.\n0\n0.2\n0.4\n0.6\n0.8\n1\n1.2\n1.4\n1.6\n1.8\n2\n3d4 3d10 3d14\n3d7 3d19 3d27\n3d4 3d10 3d14\n3d7 3d19 3d27\nSpeedup\nMax-fp16/fp32\nMG-fp16/fp32(opt)\nMG-fp16/fp32(naïve)\nMG-fp32/fp32\nARMPL@ARM or MKL@X86\nSpMV\nSpTRSV\nSpMV\nSpTRSV\nARM\nX86\nFigure 7: Ablation experiment of kernel optimization effect.\nSpeedups are over MG-fp32/fp32 (i.e., the best implementa-\ntion of full-FP32 precision).\nlaplace27\nlaplace27*108\nrhd\noil\nweather\nrhd-3T\noil-4C\nsolid-3D\n0\n0.2\n0.4\n0.6\n0.8\n1\nNormalized Time\n11\n11\n55\n20\n8\n59\n15\n15\n11\n11\n65\n20\n12\n81\n15\n16\n3.7x\n3.2x\n1.9x\n2.7x\n1.8x\n1.8x\n3.8x\n3.4x\nsetup overhead\nMG preconditioner\nother\nFigure 8: Performance of solving the linear systems on a\nsingle ARM processor. Left column: Full64. Right column:\nK64P32D16 with the setup-then-scale strategy. The #iters to\nconverge are displayed on the top of columns. Preconditioner\nspeedups are indicated within the orange column.\nlaplace27\nlaplace27*108\nrhd\noil\nweather\nrhd-3T\noil-4C\nsolid-3D\n0\n0.2\n0.4\n0.6\n0.8\n1\nNormalized Time\n11\n11\n51\n20\n8\n59\n15\n16\n11\n11\n61\n20\n12\n77\n15\n16\n3.4x\n3.2x\n2.0x\n3.0x\n1.8x\n2.3x\n3.5x\n3.7x\nsetup overhead\nMG preconditioner\nother\nFigure 9: Performance on a single X86 processor. The legends\nand meanings are the same as Figure 8.\nmix-precision setup is shown in blue columns, which shows that\nthe setup-then-scale strategy introduces only limited overhead\nduring setup. MG preconditioner speedups are indicated in orange\ncolumns, and are case-dependent. For the most idealized benchmark\nproblem, laplace27, the speedup 3.70x approaches the upper-bound\nof 4.0x shown in Table 2, because the nonzero pattern of 3d27\nrequires 27 times more memory than a vector. Laplace27∗108 needs\n60\n\n\nICPP ’24, August 12–15, 2024, Gotland, Sweden\nYi Zong, Peinan Yu, Haopeng Huang, and Wei Xue\nan additional diagonal matrix to scale and rescale, thus gaining a\nslightly lower speedup. Oil’s matrix of the 3d7 pattern occupies\na smaller proportion of memory than 3d27 and obtains a smaller\nspeedup. Meanwhile, vector PDE problems (especially oil-4C and\nsolid-3D) are more favored by FP16 because each nonzero entry\nin matrices is a small dense matrix of 𝑟× 𝑟where 𝑟is the number\nof variables located at each element. It is worth noting that the\nincreases of #iter in rhd, rhd-3T, and weather slow down their\nspeeds.\nThe E2E (i.e., the sum of setup overhead, MG preconditioner and\nother) speedup of the entire workflow could reach 2.39x, 2.21x,\n1.73x, 1.74x, 1.92x, 1.78x, 2.32x, 2.45x for these eight problems\non a single ARM processor, respectively. The results are similar on\na single X86 processor, as shown in Figure 9.\n7.4\nScalability Test\nStrong scalability is evaluated for the above problems, as shown in\nFigure 10. The parallel efficiencies on ARM of the mix-precision\nsolver could reach 96%, 89%, 63%, 99%, 98%, 71%, 93%, 62% of its\nfull-iterative-precision counterpart in terms of total time in these\nproblems, respectively. The scaling behaviors on X86 are similar.\nMix-precision solvers could maintain nearly perfect scalability in\nmedium and large-size problems. Their strong scalability will not\nsurpass that of the full-iterative-precision ones because using FP16\nin storage accelerates the computation part. In contrast, after op-\ntimization, the communication part becomes more dominant in\nE2E time. Another reason for degraded scalability in small size\nproblems is the underutilization of SIMD when there are too few\n#dof per core. The extra precision-conversion overhead inhibits the\nperformance in particularly small problems (such as the problems\nof rhd, rhd-3T and solid-3D). The expected speedups estimated by\nthe reduction of memory volumes could be observed when #dof per\ncore is large enough. In most cases, leveraging FP16 acceleration at\nthe cost of scalability is worthwhile.\n8\nDISCUSSION\nIn our experiments, both pre-smoothing and post-smoothing are\napplied once in the V-cycle. Increasing the number of smoothings\nmakes the multigrid method more time-consuming, which results in\na more significant E2E speedup when accelerated by half-precision.\nWhile some previous studies listed in Table 1 applied smoothings\nmore than once, this article maintains the numbers of smoothing\nas 1 throughout the experiments. This decision is based on the\nobservation that additional smoothings are generally less efficient\nin reducing time-to-solution for most problems.\nAnother half-precision format, BF16, is not discussed in this\narticle. This format will be worth exploring when instruction set\nsupport for BF16 becomes more extensive on CPUs. BF16 does not\nrequire scaling to avoid overflow because its range is the same as\nFP32. However, its accuracy is even worse than FP16. In our pre-\nliminary evaluation on GPUs, the #iter of using FP16 as the storage\nprecision of preconditioners is always fewer than or equal to that\nof BF16. An noticeable gap is observed that FP16 and BF16 increase\nthe #iter by 19% and 59%, respectively, compared to the full-FP64\nbaseline in the rhd problem. Therefore, FP16 appears to be more\nsuitable than BF16 for scientific computations that require strict\n# of cores\n64\n128\n256\n512\n1024\n10-1\n100\n(a) laplace27 total time\n64\n128\n256\n512\n1024\n100\n(b) laplace27*108 total time\n64\n128\n256\n512\n1024\n2048\n10-1\n100\n(c) rhd total time\n120\n240\n480\n960\n1920\n3840\n10-1\n100\n(d) oil total time\n240\n480\n960\n1920\n3840\n7680\n10-1\n100\n(e) weather total time\n64\n128\n256\n512\n1024\n2048\n10-1\n100\n(f) rhd-3T total time\n120\n240\n480\n960\n1920\n3840\n10-1\n100\n(g) oil-4C total time\n120\n240\n480\n960\n1920\n3840\n10-1\n100\n(h) solid-3D total time\nFull*-ARM\nMix16-ARM\nFull*-X86\nMix16-X86\nperfect scaling\nFigure 10: Results of strong scalability tests. The units of\ntotal times are all in seconds. ’Full*’ precision stands for the\noriginal workflow of the iterative precision, while ’Mix16’\nmeans the multigrid preconditioner is configured as the mix-\nprecision in Table 3.\nnumerical accuracy. BF16 may need more sophisticated accuracy\ncompensation techniques to leverage its advantage in reducing\nmemory access volume.\nThe transformation from AOS to SOA format discussed in Sec-\ntion 5.1 extends seamlessly to GPU implementations. Optimizing\nmemory bandwidth efficiency on GPUs necessitates coalesced mem-\nory access among threads within a warp. Thus, SOA proves more\nadvantageous even in scenarios utilizing full-FP32 precision.The\nfurther precision reduction of FP32 to FP16 for matrices is straight-\nforward based on SOA. Despite new computing hardware innova-\ntions on GPUs, represented by tensor cores, they are difficult to\nleverage significantly for sparse linear algebra kernels that are typi-\ncally memory-bounded. Similar to CPUs, our GPU implementation\nhas achieved near 100% memory bandwidth efficiency, approaching\nthe theoretical performance limit.\n9\nCONCLUSION\nMultigrid and half-precision match with each other. Multigrid pro-\nvides the numerical tolerance to lower precision due to its multi-\nlevel framework [19]. In the meantime, half-precision acceleration\nobtains much more significant end-to-end speedup because of multi-\ngrid’s dominance in runtime.\nIn this article, we investigated the use of FP16 in multigrid precon-\nditioners. Potential gains and risks are analyzed based on the matrix\n61\n\n\nFP16 Acceleration in Structured Multigrid Preconditioner for Real-World Applications\nICPP ’24, August 12–15, 2024, Gotland, Sweden\nformat and characteristics of multigrid. A complete algorithm is\nproposed with proof to avoid FP16 overflow. Our setup-then-scale\nstrategy is lightweight in setup and successfully prevents FP16 from\ninterfering with the triple-matrix products. Eight problems (3 ideal-\nized, 5 from real-world applications) are evaluated to demonstrate\nthe effectiveness of our algorithms and implementations. Geometric\naverage speedups of 2.7x and 2.8x could be obtained in MG precon-\nditioner on ARM and X86, respectively, which contributes average\nspeedups of 1.9x and 2.0x in the entire workflow. The maximal\nspeedup 3.8x in preconditioner has approached the upper-bound of\n4.0x and results in a maximal E2E speedup of 2.5x. Our algorithm\nmakes no assumptions about the background problems and can\nbe applied to other multigrids. The kernel optimization technique\ncan be directly ported to hypre’s structured-grid-specific multigrids\nsuch as SMG, PFMG, and SysPFMG.\nIt is worth mentioning that our guidelines and algorithms\nare also applicable to unstructured multigrid where matri-\nces are usually stored in CSR format. But it is difficult for\nunstructured problems to obtain comparable speedups. The\nreasons for their inadequate performance stem from two aspects.\nOn one hand, the memory volumes of the extra integer arrays in\nCSR format could not be reduced by mix-precision storage, as in-\ndicated in Table 2. On the other hand, the indirect memory access\nand inability to leverage vectorization in CSR format further de-\ncrease the bandwidth efficiency, as discussed in Section 5.1. These\ndeficiencies significantly offset the programming efforts in utilizing\nFP16 in unstructured multigrid.\nACKNOWLEDGMENTS\nWe would like to thank Xinliang Wang, Qin Wang, Zhaohui Ding,\nand others for their valuable suggestions. This work is supported by\nthe National Natural Science Foundation of China (NO.U2242210).\nWei Xue (Email: xuewei@tsinghua.edu.cn) is the corresponding\nauthor of this paper.\nREFERENCES\n[1] Ahmad Abdelfattah and et al. 2021. A survey of numerical linear algebra meth-\nods utilizing mixed-precision arithmetic. The International Journal of High\nPerformance Computing Applications 35 (Mar 2021).\n[2] China Meteorological Administration. 2016. GRAPES Numerical Weather Predic-\ntion System. Retrieved July 7, 2023 from https://www.cma.gov.cn/2011xwzx/\n2011xqxxw/2011xqxyw/202110/t20211030_4079298.html\n[3] Innovative Computing Laboratory at University of Tennessee. 2023. HPL-MXP\nmixed-precision benchmark. Retrieved March 3, 2023 from https://hpl-mxp.org/\n[4] Timothy A. Davis and Yifan Hu. 2011. The University of Florida Sparse Matrix\nCollection. ACM Trans. Math. Softw. 38, 1, Article 1 (Dec 2011), 25 pages.\n[5] Maximilian Emans and Albert van der Meer. 2010. Mixed-precision AMG as\nlinear equation solver for definite systems. Procedia Computer Science 1, 1 (2010),\n175–183. https://doi.org/10.1016/j.procs.2010.04.020 ICCS 2010.\n[6] Robert D. Falgout and Jacob B. Schroder. 2014. Non-Galerkin Coarse Grids for\nAlgebraic Multigrid. SIAM J. Sci. Comput. 36, 3 (Jan 2014), C309–C334.\n[7] Hormozd Gahvari and et al. 2012.\nModeling the Performance of an Alge-\nbraic Multigrid Cycle Using Hybrid MPI/OpenMP. In 2012 41st International\nConference on Parallel Processing. 128–137.\n[8] S. L. Glimberg and et al. 2013. A Fast GPU-Accelerated Mixed-Precision Strategy\nfor Fully Nonlinear Water Wave Computations. In Numerical Mathematics and\nAdvanced Applications 2011. Berlin, Heidelberg, 645–652.\n[9] Dominik Goddeke and Robert Strzodka. 2011. Cyclic Reduction Tridiagonal\nSolvers on GPUs Applied to Mixed-Precision Multigrid. IEEE Trans. Parallel\nDistrib. Syst. 22, 1 (jan 2011), 22–32. https://doi.org/10.1109/TPDS.2010.61\n[10] Nicholas J. Higham and Theo Mary. 2022. Mixed precision algorithms in numeri-\ncal linear algebra. Acta Numerica 31 (2022), 347–414.\n[11] Nhut-Minh Ho and et al. 2017. Exploiting half precision arithmetic in Nvidia\nGPUs. In 2017 IEEE High Performance Extreme Computing Conference (HPEC).\n[12] X. Huang and et al. 2016. P-CSI v1.0, an accelerated barotropic solver for the\nhigh-resolution ocean model component in the Community Earth System Model\nv2.0. Geoscientific Model Development 9, 11 (2016), 4209–4225.\n[13] Intel. 2018. BFLOAT16 - hardware numerics definition. Retrieved Nov 30, 2023\nfrom https://www.intel.com/content/dam/develop/external/us/en/documents/\nbf16-hardware-numerics-definition-white-paper.pdf\n[14] Carlo Janna, Andrea Comerlati, and Giuseppe Gambolati. 2009. A Comparison\nof Projective and Direct Solvers for Finite Elements in Elastostatics. Adv. Eng.\nSoftw. 40, 8 (aug 2009), 675–685. https://doi.org/10.1016/j.advengsoft.2008.11.010\n[15] Lawrence Livermore National Lab. 2023. Documentation for hypre. Retrieved\nMarch 3, 2023 from https://hypre.readthedocs.io/en/latest\n[16] Lawrence Livermore National Lab. 2023. Structured multigrid in HYPRE. Re-\ntrieved March 3, 2023 from https://hypre.readthedocs.io/en/latest/solvers-smg-\npfmg.html\n[17] Ruipeng Li and Ulrike Meier Yang. 2021. Performance Evaluation of hypre Solvers.\n(Feb 2021). https://doi.org/10.2172/1764323\n[18] Daniel Lowell and et al. 2013. Stencil-Aware GPU Optimization of Iterative\nSolvers. SIAM Journal on Scientific Computing 35, 5 (2013), S209–S228.\n[19] Stephen F. McCormick, Joseph Benzaken, and Rasmus Tamstorf. 2020. Algebraic\nerror analysis for mixed-precision multigrid solvers. SIAM J. Sci. Comput. 43\n(2020), S392–S419.\n[20] Paulius Micikevicius, Sharan Narang, Jonah Alben, Gregory Diamos, Erich\nElsen, David Garcia, Boris Ginsburg, Michael Houston, Oleksii Kuchaiev, Ganesh\nVenkatesh, and Hao Wu. 2018. Mixed Precision Training. arXiv:1710.03740 [cs.AI]\n[21] Eike H. Müller and et al. 2014. Massively parallel solvers for elliptic partial\ndifferential equations in numerical weather and climate prediction. Quarterly\nJournal of the Royal Meteorological Society 140, 685 (2014), 2608–2624.\n[22] M. Naumov and et al. 2015. AmgX: A Library for GPU Accelerated Algebraic\nMultigrid and Preconditioned Iterative Methods. SIAM Journal on Scientific\nComputing 37, 5 (2015), S602–S626.\n[23] Society of Petroleum Engineers. 2023. SPE Comparative Solution Project. Re-\ntrieved March 3, 2023 from https://www.spe.org/web/csp/datasets/set02.htm\n[24] Kyaw Linn Oo and Andreas Vogel. 2020. Accelerating Geometric Multigrid Pre-\nconditioning with Half-Precision Arithmetic on GPUs. arXiv:2007.07539 [cs.MS]\n[25] China Meteorological News Press. 2014. An Introduction of GRAPES. Retrieved\nJuly 7, 2023 from https://www.cma.gov.cn/en/NewsReleases/MetInstruments/\n201403/t20140327_241784.html\n[26] Trilinos project. 2023. MueLu.\nRetrieved Nov 26, 2023 from https://trilinos.\ngithub.io/muelu.html\n[27] Christian Richter and et al. 2014.\nGPU-accelerated mixed precision alge-\nbraic multigrid preconditioners for discrete elliptic field problems. In 9th IET\nInternational Conference on Computation in Electromagnetics. 1–2.\n[28] Yousef Saad. 2003. Iterative Methods for Sparse Linear Systems (second ed.).\nSociety for Industrial and Applied Mathematics.\n[29] Martin H. Sadd. 2005. Elasticity: Theory, Applications, and Numerics. Academic\nPress. https://doi.org/10.1016/B978-0-12-605811-6.X5000-3\n[30] K. Stuben. 2000. Algebraic Multigrid (AMG) : An Introduction With Applications.\n[31] MFEM team. 2023. MFEM examples.\nRetrieved Nov 26, 2023 from https:\n//mfem.org/examples\n[32] Ulrich Trottenberg and et al. 2001.\nMultigrid. Academic Press, San Diego,\nCalifornia, USA.\n[33] Yu-Hsiang Mike Tsai and et al. 2023. Three-precision algebraic multigrid on\nGPUs. Future Generations Computer Systems 149 (12 2023).\n[34] Xiaowen Xu and et al. 2017. Algebraic interface-based coarsening AMG precon-\nditioner for multi-scale sparse matrices with applications to radiation hydrody-\nnamics computation. Numerical Linear Algebra with Applications 24, 2 (2017),\ne2078. https://doi.org/10.1002/nla.2078\n[35] Takateru Yamagishi and et al. 2016. GPU Acceleration of a Non-Hydrostatic\nOcean Model with a Multigrid Poisson/Helmholtz Solver. Procedia Comput. Sci.\n80, C (jun 2016), 1658–1669. https://doi.org/10.1016/j.procs.2016.05.502\n[36] Ulrike Meier Yang. 2010. On long-range interpolation operators for aggressive\ncoarsening. Numerical Linear Algebra with Applications 17, 2-3 (2010), 453–472.\n[37] Xiaojian Yang and et al. 2023. Optimizing Multi-Grid Computation and Paral-\nlelization on Multi-Cores. In Proceedings of the 37th International Conference\non Supercomputing (ICS ’23). 227–239. https://doi.org/10.1145/3577193.3593726\n[38] Chensong Zhang and et al. 2023. OpenCAEPoro. Retrieved March 3, 2023 from\nhttps://github.com/OpenCAEPlus/OpenCAEPoro/tree/main/examples/spe10\n[39] Qianchao Zhu and et al. 2021.\nEnabling and Scaling the HPCG Bench-\nmark on the Newest Generation Sunway Supercomputer with 42 Million Het-\nerogeneous Cores. In Proceedings of the International Conference for High\nPerformance Computing, Networking, Storage and Analysis (SC ’21). ACM, Ar-\nticle 57, 13 pages. https://doi.org/10.1145/3458817.3476158\n[40] Yi Zong and et al. 2024. POSTER: StructMG: A Fast and Scalable Structured\nMultigrid. In Proceedings of the 29th ACM SIGPLAN Annual Symposium on\nPrinciples and Practice of Parallel Programming, PPoPP 2024. ACM, 478–480.\nhttps://doi.org/10.1145/3627535.3638482\n62\n\n\nFuture Generation Computer Systems 149 (2023) 280–293\nContents lists available at ScienceDirect\nFuture Generation Computer Systems\njournal homepage: www.elsevier.com/locate/fgcs\nThree-precision algebraic multigrid on GPUs\nYu-Hsiang Mike Tsai a,∗, Natalie Beams b, Hartwig Anzt b,a\na Karlsruhe Institute of Technology, Kaiserstraße 12, 76131, Karlsruhe, Germany\nb Innovative Computing Laboratory, University of Tennessee, 1122 Volunteer Blvd, 37996 TN, Knoxville, USA\na r t i c l e\ni n f o\nArticle history:\nReceived 16 February 2023\nReceived in revised form 15 July 2023\nAccepted 19 July 2023\nAvailable online 23 July 2023\nKeywords:\nAlgebraic multigrid\nMixed precision\nPortability\nGPUs\na b s t r a c t\nRecent research has demonstrated that using low precision inside some levels of an algebraic multigrid\n(AMG) solver can improve performance without negatively impacting the AMG quality. In this paper,\nwe build upon previous research and implement an AMG that can use double, single, and half precision\nfor the distinct multigrid levels. The implementation is platform-portable across GPU architectures\nfrom AMD, Intel, and NVIDIA. In an experimental analysis, we demonstrate that the use of half precision\ncan be a viable option in multigrid. We evaluate the performance of different AMG configurations and\ndemonstrate that mixed precision AMG can provide runtime savings compared to a double precision\nAMG.\n© 2023 Elsevier B.V. All rights reserved.\n1. Introduction\nFinite element simulations drive a large portion of computer-\nassisted research and development. The principle behind finite\nelement simulations is the discretization of the weak form of the\ngoverning partial differential equation(s) on a mesh and the com-\nputation of solution approximations through local basis functions\non each mesh element. Computing these approximations requires\nthe solution of a linear system that reflects the contributions\nof the local basis functions to the global solution. Given the\nproperties of the finite element discretization, the local solutions\nare only directly coupled to their neighboring mesh elements,\nresulting in linear systems that are typically large and sparse.\nFor one-dimensional finite element simulations, direct sparse\nlinear solvers based on Gaussian elimination are generally a good\nchoice, but for two- and three-dimensional simulations, the large\nfill-in occurring during the factorization makes direct solvers\nunattractive. Instead, iterative solvers generating a sequence of\nsuccessively-better solution approximations are often preferred.\nA successful iterative solver must find a balance between how\nquickly each iteration can be performed and how many iterations\nare required to reach the desired accuracy tolerance. A simple\niterative process may be very cheap to apply on a per-iteration\nbasis but do a poor job in terms of reducing some (or all) frequen-\ncies of the error with each iteration. This motivates the concept\nof multigrid methods [1–4], which use a hierarchy of successively\nsmaller problems (coarse grids), related through operators that\nrestrict from fine to coarse and prolongate back. These different\nlevels of refinement allow the same iterative solvers to target\n∗Corresponding author.\nE-mail address:\nyu-hsiang.tsai@kit.edu (Y.-H.M. Tsai).\ndifferent frequencies of the total error at the same time, based\non how the error is represented on each grid. Traditionally, one\ndistinguishes two classes of multigrid methods: those that use the\ngeometric mesh information to derive the hierarchy of grids are\ncalled ‘‘Geometric Multigrid Methods’’ (GMG); those that derive\nthe hierarchy of grids exclusively from the large sparse matrix\nare called ‘‘Algebraic Multigrid Methods’’ (AMG). AMG methods\nare particularly effective for discretizations of elliptic PDEs on un-\nstructured grids and problems that lack an underlying geometric\nmesh.\nMixed precision methods for numerical computing have a\nlong history and are still an active area of research [5]. With\nthe widespread use of GPU accelerators in scientific computing,\nmixed precision methods exploiting their lower (single, half, or\nother formats like bfloat16) precision capabilities while retaining\nsuitable accuracy are of particular interest recently [6]. When\nusing mixed precision in the context of solving linear systems,\nperhaps the most common choice is a form of iterative refine-\nment (IR). Originally proposed as a way to reduce accumulated\nround-off errors [7,8], it is a framework in which the solution is\nsuccessively refined through corrections found by solving a linear\nsystem with the residual of the previous iteration as the right-\nhand side. These ‘‘inner’’ solves – taking place inside the outer\niterative process – can often be done at a lower precision than\nthe calculation of the residual vectors and updates to the solution.\nThis offers the opportunity to accelerate the time to solution, as\nhas been demonstrated on both CPUs [9,10] and GPUs [11,12].\nPrevious mixed precision multigrid work on GPUs has gen-\nerally focused on using a reduced or mixed precision MG cycle\nin the inner solve of a higher-precision IR process [13–15]. This\nstrategy is given careful theoretical consideration by Tamstorf\net al. in [16,17]. In most cases, the multigrid itself is not mixed\nhttps://doi.org/10.1016/j.future.2023.07.024\n0167-739X/© 2023 Elsevier B.V. All rights reserved.\n\n\nY.-H.M. Tsai, N. Beams and H. Anzt\nFuture Generation Computer Systems 149 (2023) 280–293\nprecision, meaning that the entire MG solver uses the lower\nprecision, rather than combining different precisions on the dif-\nferent levels. Oo and Vogel [15] consider several configurations of\nmixed precision GMG for the inner solve of IR, with up to three\ndifferent precisions on different levels, including half precision.\nThey try both ‘‘directions’’ for the lower precision use: on the\ncoarsest grids, with higher precisions on the finer grids, and\nvice versa. The former is similar in spirit to our three-precision\nconfigurations. While not a GPU implementation, a notable work\nof mixed precision MG apart from an IR framework is Buttari\net al. [18], where the authors deploy single precision and block\nlow rank factorization to decrease the coarse grid solution time\nof a very large, distributed GMG solver. In our current work,\nwe also do not use lower or mixed precision multigrid inside\niterative refinement, opting instead to test its use as a ‘‘drop\nin’’ replacement for any other double precision preconditioner\nin a conjugate gradient (CG) solver. We also briefly consider the\nperformance of mixed precision AMG as a stand-alone solver in\nSection 4.4.\nCompared to existing literature, this paper presents the fol-\nlowing novel contributions:\n1. We evaluate the use half precision within AMG with re-\nspect to preconditioner accuracy and CG convergence;\n2. We demonstrate that using half precision for some AMG\ncomponents can be a viable option even if aiming for\nhigh-accuracy solutions;\n3. We compare the performance of mixed precision AMG on\nAMD MI250X, Intel PVC, and NVIDIA H100 GPUs.\nThe rest of the paper is structured as follows: Following some\nbrief background in Section 2, details of the design of Ginkgo’s\nAMG, showcasing its flexibility in terms of mixed precision con-\nfiguration, are given in Section 3. Section 4 evaluates the effec-\ntiveness and performance of the mixed precision AMG when used\nas a preconditioner inside a CG iterative solver on NVIDIA, AMD,\nand Intel GPUs. We then briefly discuss using mixed precision\nAMG as a standalone solver and the potential of mixed precision\nin other types of multigrid cycles. The findings are summarized\nin Section 5.\n2. Background on AMG\nAlgebraic multigrid (AMG) [2,4] is a popular choice for solving\nor preconditioning linear problems originating from finite ele-\nment discretizations. Unlike geometric multigrid (GMG), which\nrelies on using information about the underlying geometric mesh,\nan AMG solver is constructed directly from the sparse system ma-\ntrix. Similarly to GMG, AMG builds a hierarchy of consecutively-\ncoarser grids and computes error correction terms on the coarser\ngrids to improve the solution on finer grids. After creating the\nhierarchy, we can use it for a multigrid method like the V-\nCycle in Algorithm 1. It restricts the residual on a fine grid to\na coarser grid, then uses the coarser grid to obtain an error\ncorrection that is prolongated back to the finer grid to update the\nsolution approximation. These correction computations generally\nentail a few iterations of an iterative method, called a ‘‘smoother’’\nbecause it acts to smooth the high-frequency errors on the scale\nof that grid, while the coarsest grid may opt for a direct solve of\nthe restricted problem, which is much smaller than the original\nmatrix.\nOur previous work [19] developed an AMG implementation\nas part of the Ginkgo library that allows the user to employ\nsingle precision for coarse (i.e., not the finest) multigrid levels.\nWe now present an extension of this work, enabling half precision\ncomputations. This results in the first AMG implementation that\ncan employ double, single, and half precision computations for\nAlgorithm 1 V-cycle multigrid method. We use blue, red, and\nbrown colors, respectively, to indicate the precision for: ma-\ntrices (A), working vectors on the fine level (r), and working\nvectors on the next coarsest level (g, e). The presmoothers and\npostsmoothers also use the working precision for computation\nand matrix precision for the system matrix.\n1: procedure Vcycle(A, x, b)\n2:\nx = PreSmooth(x, b)\n3:\nr = b - Ax\n4:\ng = Restrict(r)\n5:\ne = zero\n6:\nVcycle(Coarse, e, g)\n7:\nx += Prolong(e)\n8:\nx = PostSmooth(x, b)\n9: end procedure\nFig. 1. The MultigridLevel class with its components.\nthe distinct multigrid levels on AMD, Intel, and NVIDIA GPUs.\nWe demonstrated in our first paper that Ginkgo AMG can be\ncompetitive with existing open source AMG implementations –\nnamely, NVIDIA’s AmgX [20] and Lawrence Livermore National\nLaboratory’s HYPRE [21] – on an NVIDIA GPU. In this work,\nwe focus on the evaluation of different precision configurations,\nincluding the use of half precision, within Ginkgo AMG.\n3. Design of the flexible and platform-portable AMG\nThe design of the Ginkgo AMG is driven by three main goals:\nflexibility, performance, and platform portability.\nIn Ginkgo, we define a MultigridLevel class, visualized in\nFig. 1, that contains a fine grid matrix from which it constructs\nthe coarse grid matrix (C) via the coarsening algorithm, as well as\nthe restriction (R) and prolongation (P) operators. The fine matrix\n(F) is the given matrix or the coarse matrix from the finer level,\nbut with the storage precision or format potentially altered by\nthe MultigridLevel according to the algorithm requirement or\nsettings. The prolongation operator (P) is an n × m matrix built\nfrom F, and the restriction operator (R) is an m×n matrix defined\nas PT . The coarse matrix (C) is an m×m matrix formed as C = RFP,\ni.e., standard Galerkin coarsening. The MultigridLevel class is\nvisualized with its key components in Fig. 1.\nGinkgo’s AMG implementation allows the use of different\nprecision formats on different levels of the multigrid hierarchy,\nresulting in a mixed precision AMG. The precision conversion\nhappens on-the-fly in the restriction and prolongation operations.\nThe gray and red portions of Listing 1 show a standard AMG\nV-cycle with a max level depth of 10, a smoother sm that is\nused for all smoothing operations, a MultigridLevel mg_lvl,\nand a coarse level solver, coarsest_solver. The number of\nsmoothing sweeps is a parameter of the smoother object sm.\nIn the configuration shown in the gray and green parts of\nListing 1, we enable mixed precision by adding two Multigrid-\nLevels and two smoothers in the configuration list, with _f\n281\n\n\nY.-H.M. Tsai, N. Beams and H. Anzt\nFuture Generation Computer Systems 149 (2023) 280–293\n1 multigrid::build()\n2\n.with_max_levels(10u)\n/ /\ne q u a l\nt o\nNVIDIA / AMGX\n11\nmax\nl e v e l s\n3\n.with_min_coarse_row(64u)\n4 DP\n.with_pre_smoother(sm)\n5 |\n.with_mg_level(mg_lvl)\n6 DP\n.with_coarest_solver(coarest_solver)\n7 MP\n.with_pre_smoother(sm, sm_f)\n8 ||\n.with_mg_level(mg_lvl , mg_lvl_f)\n9 ||\n.with_coarest_solver(coarest_solver_f)\n10 ||\n.with_level_selector(\n11 ||\n[](const size_type level, const LinOp*)\n12 ||\n-> size_type {\n13 ||\nreturn level >= 1 ? 1 : 0;\n14 MP\n})\nListing 1:\nConfiguration of a Ginkgo Multigrid object.\nLines with a red background are used when configuring for\ndouble precision (DP), while the green background indicates\nconfiguration for mixed precision (MP).\nindicating ‘‘float’’ or single precision. We also need to config-\nure the level_selector to describe the desired scheme. Here,\nwhen the level is larger than or equal to 1 (that is, all levels except\nthe finest grid), we use the second pair (mg_lvl_f, sm_f).\nWhen the level is less than 1, we use the first pair (mg_lvl, sm).\nTaken together, this configuration generates a mixed precision\nMultigrid where only the finest level is using double precision,\nand all other levels use single precision. This particular mixed\nprecision Multigrid configuration allows for smooth integration\nas a preconditioner into an iterative solver using double precision,\nas the input and output vectors, as well as the original matrix,\nremain in double precision.\nAs smoother applications in the form of vector operations\nare relatively cheap, the runtime of an AMG cycle is gener-\nally dominated by the residual computations that involve sparse\nmatrix–vector multiplications (SpMVs). We implement several\ncommon optimization techniques to reduce the overhead of these\noperations. Since memory allocation on the GPU is known to be\nexpensive [22,23], we use a workspace for allocating the opera-\ntor components and intermediate operations. The workspace is\navailable over the complete lifetime of the operator.\nWe also reduce the overall number of residual evaluations\nwhen possible: we skip residual computation if an initial guess\nis zero, because we know the residual will be equal to the right-\nhand side vector. When using AMG as a preconditioner, we only\ncompute the explicit residual if the user requests it, e.g. for mon-\nitoring convergence. A few optimizations regarding the handling\nof residual vectors are specific to the machinery of the Ginkgo\nlibrary itself. For example, if we need the residual in the AMG\nsolver, and an internal component already computed it, this resid-\nual is ensured to be accessible from outside the component. We\nalso ‘‘split’’ the iteration termination check within Ginkgo, such\nthat reaching an iteration limit terminates the algorithm before\nthe residual is computed for convergence checking purposes.\nTo enable both platform and performance portability, we im-\nplement the Ginkgo AMG using a backend model as described\nin [24,25], where we complement an algorithm skeleton invoking\na sequence of subroutines with backends containing the corre-\nsponding subroutines as heavily optimized GPU kernels in the\nvendor-native programming languages. Specifically, we imple-\nment CUDA kernels for NVIDIA GPUs, DPC++ kernels for Intel\nGPUs, and HIP kernels for AMD GPUs. Instead of having three\ncomplete stand-alone AMG implementations for the distinct GPU\narchitectures with the corresponding kernel sets, we use C++\nruntime polymorphism for automatically selecting and invoking\nthe suitable kernels when executing the AMG algorithm. This\nFig. 2. Meshes used for MFEM diffusion experiments. Left: L-shape mesh with 7\nlevels of uniform refinement (49,152 elements); Right: Beam mesh with 3 levels\nof uniform refinement (4,096 elements).\nallows the deployment of the AMG solver without having to\nmaintain different variants for different hardware architectures.\nBy doing so, we keep the cutting-edge features from vendors’\nofficial compilers without waiting for another compiler to adopt\nthe new features.\n4. Experimental evaluation\nWe consider a set of benchmark problems from the SuiteS-\nparse Matrix Collection [26], as well as two diffusion problems\nexported from the MFEM finite element library. MFEM [27,28] is a\npopular open-source finite element library with support for high-\norder meshes and basis functions, among many other features.\nOur exported problems come from MFEM’s ‘‘example 1’’, solving\na standard diffusion problem −∇· (c∇u) = 1, where c is a given\ncoefficient. We use homogeneous Dirichlet boundary conditions.\nTwo of MFEM’s provided meshes are tested; they are shown in\nFig. 2. For the ‘‘L-shape’’ mesh, a constant coefficient of c =\n1 is used, while the ‘‘beam’’ mesh uses a piecewise constant\ncoefficient with a jump from 1 to 0.1 at the midpoint of the\nlength of the beam. All tests use standard third order tensor-\nproduct basis functions on the Legendre–Gauss–Lobatto nodes\nand MFEM’s default choices for quadrature points based on the\norder of basis functions. A summary of all matrices is provided\nin\nTable 1, including their sizes and the range of the absolute\nvalues of all their nonzeros.\nWe evaluate Ginkgo’s AMG implementation on GPU architec-\ntures from AMD, Intel, and NVIDIA. The GPUs and corresponding\ncompilers are listed along with some key characteristics in Ta-\nble 2. Note that the AMD MI250X has two graphics compute dies\n(GCDs), but the GCDs are seen individually by the system. In the\nfollowing, we always consider just one GCD. Similarly, one Intel\nPVC GPU contains two tiles, but we run exclusively on one tile.\nFor Intel GPUs, to our best knowledge, Ginkgo is currently the\nonly library providing an AMG implementation with all portions\nof the calculation taking place on the GPU.\nThe experimental results are arranged as follows: first, we\ndiscuss a simple model for the potential performance gains of\nmixed precision AMG based on the performance of sparse matrix–\nvector multiplication (SpMV) operations and data storage. Next,\nmixed precision AMG (V-cycle) is used as a preconditioner inside\na CG solver for our selected benchmark problems, and we analyze\nthe effects on convergence rate and time per iteration. We discuss\nthe challenges of using half precision in AMG and strategies\nfor mitigation. Finally, we briefly consider AMG as a standalone\nsolver and the potential for mixed precision in other multigrid\ncycles, like the W-cycle.\n4.1. SpMV performance as proxy for AMG\nFor an idea of the potential performance benefits of mixed\nprecision AMG on our target architectures, we analyze the per-\nformance of the sparse matrix–vector operation (SpMV), which is\n282\n\n\nY.-H.M. Tsai, N. Beams and H. Anzt\nFuture Generation Computer Systems 149 (2023) 280–293\nTable 1\nMatrix characteristics for the selected problems.\nProblem\nSize\nElements\nabs. value range\nMFEM\nbeam (-o3 -l3)\n120,625\n14,070,001\n8.9×10−7–1.0\nL-shape (-o3 -l7)\n443,905\n11,066,881\n4.0×10−3–6.0\n2cubes_sphere\n101,492\n1,647,264\n6.7×10−15–2.5×1010\nthermal2\n1,228,045\n8,580,313\n1.7×10−7–4.9\nSuiteSparse\ncage14\n1,505,785\n27,130,349\n0.011–0.94\ncage13\n445,315\n7,479,343\n0.012–0.93\noffshore\n259,789\n4,242,673\n7.2×10−21–7.5×1014\ntmt_sym\n726,713\n5,080,961\n8.5×10−14–19\nTable 2\nGPU characteristics.\nGPU\nPeak Perf. (DP)\nPeak Perf. (SP)\nMem. size\nBandwidth\nCompiler\nAMD MI250Xa (1 GCD)\n24 TFLOP/s\n24 TFLOP/s\n64 GB\n1.6 TB/s\nHIP 5.3\nNVIDIA H100b (PCIE)\n26 TFLOP/s\n51 TFLOP/s\n80 GB\n2.0 TB/s\nCUDA 12.0\nIntel PVCc (1 Tile)\n22.8 TFLOP/s\n22.8 TFLOP/s\n64 GB\n1.6 TB/s\nDPC++ 2023.1\nNote: PVC is not released as of this writing, so the official characteristics may be changed when released.\naFrom Frontier, Oak Ridge National Laboratory, USA.\nbFrom BwUniCluster 2.0, bwHPC, Germany.\ncFrom Sunspot, Argonne National Laboratory, USA.\nused in residual computations. We expect the residual computa-\ntions on each level to be the major limiting factor of performance;\nin comparison, the solution phase of our smoothers, as well as the\nrestriction and prolongation operations, are simpler. We extract\nthe original matrices and the coarse matrices from levels 1–\n10 in an 11-level multigrid. After 2 warmup applications, we\naverage the results from 10 SpMV evaluations. The SpMV is a\nmemory-bound operation, with performance tied to the amount\nof memory traffic required rather than the number of floating\npoint operations performed. The memory involved in a SpMV\noperation of an n × n CSR matrix with nnz non-zero elements\ncan be written as:\n(nnz + n + 1) × I + (nnz + 2n) × V.\nI and V are the size, in bytes, of the IndexType and ValueType\nused for the CSR matrix storage. When nnz ≫n, the storage is\napproximated as nnz ×(I +V). Thus, in terms of reduced memory\ntraffic, the speedup of a lower precision ValueType V2 over higher\nprecision V1 can be estimated as\nnnz(I + V1)\nnnz(I + V2) = I + V1\nI + V2\n.\nFor example, the speedup of single precision SpMV over double\nprecision SpMV would be 1.5x, and the speedup of half precision\nSpMV over double precision SpMV would be 2x.\nGinkgo has recently introduced native half precision support.\nThe classical CSR SpMV kernel in Ginkgo is shown in Listing 2.\nIn Fig. 4 we visualize the speedup of single and half precision\nSpMVs over double precision on the H100 GPU for the matrices\nof the distinct multigrid levels for each of the test matrices.\nThe matrix size decreases for coarser multigrid levels, numbered\nwith higher level numbers, as shown in Fig. 3. For the classical\nCSR implementation, the performance difference between the\nsingle and half precision SpMVs is negligible on H100. A possible\nexplanation is that NVIDIA GPUs are optimized for each thread\naccessing at least 4 bytes of memory, which forms a 128 byte\ncache line. To adjust for this, we consider a ‘‘packed’’ version of\nCSR SpMV. The implementation, shown in Listing 3, processes\ntwo half precision values in one memory access. With the packed\nvariant, the performance of the half precision SpMV generally\noutperforms that of single precision on H100 in Fig. 5, especially\nfor larger cases, which is in the left side in the plots. We note\nthat speedup values can exceed expectations due to cache effects.\nFigs. 6 and 7 combine the performance data and arrange the\n1\ntemplate < i n t\nsubwarp_size >\n2 classical_csr_spmv(row_ptrs , col_idxs , val, b, c)\n{\n3\n/ /\nc r e a t e\na\nsubwarp\nw i t h\nt h e\nsubwarp\ns i z e\n4\nauto subwarp_tile = ...;\n5\n/ /\ng e t\nt h e\ni n d e x\no f\nt h r e a d\ni n\na\nsubwarp\n6\nconst\nauto subid = ...;\n7\n/ /\ng e t\nt h e\np r o c e s s i n g\nrow\ni n d e x\n8\nauto row = ...;\n9\nconst\nauto ind_end = row_ptrs[row + 1];\n10\nValueType temp_val = zero<ValueType >();\n11\n/ /\neach\nt h r e a d\na c c u m u l a t e s\nt h e\nr e s u l t ;\ns w i t c h\nn e x t\nwindow\nby\ns u b w a r p _ s i z e\n12\nf o r ( auto ind = row_ptrs[row] + subid; ind <\nind_end; ind += subwarp_size) {\n13\ntemp_val += val[ind] * b[col_idxs[ind]];\n14\n}\n15\n/ /\nuse\ns h u f f l e\nt o\ng e t\nt h e\nsummation\no f\na l l\nt h r e a d s\ni n\na\nsubwarp\n16\nauto subwarp_result = ...;\n17\ni f\n(subid == 0) {\n18\n/ /\nw r i t e\nt h e\nr e s u l t\nt o\nmemory\n19\nc[row] = subwarp_result;\n20\n}\n21 }\nListing 2: Classical CSR SpMV as implemented in Ginkgo\nFig. 3. The ratio of the number of stored elements in each level of a multigrid\nhierarchy to that of the finest level (level 0), using Parallel Graph Match [20]\nwith size 2 for aggregation.\n283\n\n\nY.-H.M. Tsai, N. Beams and H. Anzt\nFuture Generation Computer Systems 149 (2023) 280–293\n1\ntemplate < i n t\nsubwarp_size >\n2 classical_csr_spmv_pack(row_ptrs , col_idxs , val,\nb, c) {\n3\n/ /\nc r e a t e\na\nsubwarp\nw i t h\nt h e\nsubwarp\ns i z e\n4\nauto subwarp_tile = ...;\n5\n/ /\ng e t\nt h e\ni n d e x\no f\nt h r e a d\ni n\na\nsubwarp\n6\nconst\nauto subid = ...;\n7\n/ /\ng e t\nt h e\np r o c e s s i n g\nrow\ni n d e x\n8\nauto row = ...;\n9\nconst\nauto ind_end = row_ptrs[row + 1];\n10\nValueType temp_val = zero<ValueType >();\n11\n/ /\neach\nt h r e a d\na c c u m u l a t e s\nt h e\nr e s u l t ;\ns w i t c h\nn e x t\nwindow\nby\ns u b w a r p _ s i z e ∗2\n12\n/ /\neach\nt h r e a d\nh a n d l e s\ntwo\ne l e m e n t s\nf o r\nh a l f\n13\nf o r ( auto ind = row_ptrs[row] + subid * 2;\nind < ind_end; ind += subwarp_size * 2)\n{\n14\ntemp_val += val[ind] * b[col_idxs[ind]];\n15\ni f\n(ind + 1 < ind_end) {\n16\n/ /\nAdd\nt h e\nn e x t\ne l e m e n t\ni f\na v a i l a b l e\n17\ntemp_val += val[ind+1] * b[col_idxs[\nind+1]];\n18\n}\n19\n}\n20\n/ /\nuse\ns h u f f l e\nt o\ng e t\nt h e\nsummation\no f\na l l\nt h r e a d s\ni n\na\nsubwarp\n21\nauto subwarp_result = ...;\n22\ni f\n(subid == 0) {\n23\n/ /\nw r i t e\nt h e\nr e s u l t\nt o\nmemory\n24\nc[row] = subwarp_result;\n25\n}\n26 }\nListing 3:\nClassical CSR SpMV: packed variant with altered\nmemory access\nFig. 4. H100 speedup of single and half precision SpMVs compared to double\nprecision for the matrices produced on each level of AMG.\nmatrices according to increasing nonzero count. Fig. 7 shows clear\nimprovement for the packed half variant compared to Fig. 6. For\nboth versions, we see noticeable speedup for matrices containing\nmore than 2e5 elements. The additional condition in line 15\nof Listing 3 is only known at runtime, and is not based solely\non the index of the thread, which poses a challenge for the\ncompiler. This additional complexity compared to the classical\nSpMV kernel could be a factor in the cases where the packed\nvariant performs worse than single precision, particularly since\nthis happens for the smaller matrices. With less total memory\nmovement involved, there is less potential speedup to cover for\nthe effects of the extra conditional statement.\nWe repeat the same experiments on one GCD of an AMD\nMI250X GPU. Here, the original CSR implementation achieves\nhigher performance in half precision than in single precision,\nFig. 5. H100 speedup of single and half precision SpMVs compared to double\nprecision for the matrices produced on each level of AMG. The half precision\nSpMV uses the ‘‘packed half’’ variant.\nFig. 6. H100 speedup of single (◦) and half precision (×) SpMVs for all matrices\nin the AMG hierarchies, arranged by total number of nonzeros.\nFig. 7. H100 speedup of single (◦) and half precision (×) SpMVs for all matrices\nin the AMG hierarchies, arranged by total number of nonzeros. The half precision\nSpMV uses the ‘‘packed half’’ variant.\nshown in Fig. 8, and the packed variant does not improve per-\nformance compared to the original version. Similar to the NVIDIA\n284\n\n\nY.-H.M. Tsai, N. Beams and H. Anzt\nFuture Generation Computer Systems 149 (2023) 280–293\nFig. 8. MI250X speedup of single and half precision SpMVs compared to double\nprecision for the matrices produced on each level of AMG.\nFig. 9. MI250X speedup of single (◦) and half precision (×) SpMVs for all\nmatrices in the AMG hierarchies, arranged by total number of nonzeros.\nH100 experiments, we observe significant performance benefits\nonly for matrices containing more than 2e5 elements in Fig. 9.\nLike the H100, using the packed variant for half precision\ngives better performance than the original version for CSR SpMV\non one tile of PVC in Figs. 10 and 11. The finest matrix of\ncage14 can reach around 6x speedup for half precision compared\nto double precision. This may be related to Ginkgo’s choice of\nsubwarp_size for half precision compared to double precision.\nIn Ginkgo, the subwarp size (1, 16, or 32) is chosen according\nto the average number of non-zeros per rows in the matrix for\nPVC. Intel PVC only supports 16 and 32 for subgroup (analogous\nto a warp in NVIDIA terminology) size, and Intel does not support\nthe same sub-subgroup features as subwarp in CUDA currently.\nFor the SpMV kernel on PVC, Ginkgo uses the subgroup as a\nsubwarp safely because there is no communication out of the\nsubgroup/subwarp. Size 1 is still available for PVC because it\ncorresponds to one thread per row. Because threads handle two\nmatrix elements each in half precision, the kernel can assign a\nsmaller subwarp size in half precision when the average num-\nber of non-zeros per row is less than 32, and this smaller size\nmay perform particularly well for some matrices. However, half\nprecision is slower than single precision for several matrices.\nMoreover, the single precision speedup is quite close to 1. In this\ncase, the kernel may use too large of a subwarp size for the short\nrows such that we do not fully utilize the GPU. The CUDA and\nFig. 10. PVC speedup of single and half precision SpMVs compared to double\nprecision for the matrices produced on each level of AMG. The half precision\nSpMV uses the ‘‘packed half’’ variant.\nFig. 11. PVC speedup of single (◦) and half precision (×) SpMVs for all matrices\nin the AMG hierarchies, arranged by total number of nonzeros. The half precision\nSpMV uses the ‘‘packed half’’ variant.\nHIP backends may select subwarp sizes of 2, 4, or 8, while PVC is\nlimited to choosing between 1 and 16 for short rows, due to the\nlimitations mentioned above.\n4.1.1. Extension to mixed precision AMG\nWe can compute storage approximations for an entire AMG\ncycle based on the number of SpMVs on each level. In our im-\nplementation, coarsening usually aggregates two nodes together,\nwhile the exact compression ratio depends on the sparsity pat-\ntern. Fig. 3 reports the relative nonzero count of the matrices in\nthe distinct AMG levels in comparison to a compression ratio of\n2. Accumulating over an AMG hierarchy of N + 1 levels, we have\nN\n∑\ni=0\nCi\n1\n2i {(nnz + n + 1) × I + (nnz + 2n) × Vi}\n≈\nN\n∑\ni=0\nCi\nnnz × (I + Vi)\n2i\n,\nwhere Vi is the size of the precision on level i and Ci is a constant\ndetermined by the total number of SpMVs performed on that\nlevel. For example, in the Jacobi smoother V-cycle configuration\n285\n\n\nY.-H.M. Tsai, N. Beams and H. Anzt\nFuture Generation Computer Systems 149 (2023) 280–293\nbelow, we have Ci = 2 for all but the coarsest level, coming from\none pre- and one post-smoother application on those levels. From\nhere, we can compute an estimate for the potential speedup of\nthe mixed precision configuration as\n∑N\ni=0 Ci\n(I+V)\n2i\n∑N\ni=0 Ci\n(I+Vi)\n2i\n.\nExtending the work of [19] to also include the use of half\nprecision, we begin with four mixed precision settings:\n1. (DP): all levels use double precision.\n2. (DP-SP): The first level uses double precision, and all\ncoarser levels use single precision.\n3. (DP-SP-HP): The first level uses double precision, the sec-\nond level uses single precision, and all coarser levels use\nhalf precision.\n4. (DP-HP): The first level uses double precision, and all\ncoarser levels use half precision.\nWe call these settings ‘‘uniform level’’, meaning that the ma-\ntrix and vector working precision (marked in blue and red in\nAlgorithm 1) are always the same on a particular level. For the\n(DP-SP) and (DP-HP) configurations described here, the SpMV-\nbased speedup estimate would predict a speedup of 1.2x and 1.3x\nover double precision AMG, respectively.\n4.2. Challenges of half precision in AMG\nWe first consider a straightforward continuation of the AMG\nexperiments in [19]. For defining the restriction operators, the\nAMG implementation in Ginkgo uses parallel graph match (PGM),\nwhich was introduced by Naumov et al. [20] as a GPU-based\nalgorithm for deriving a coarse approximation through explo-\nration of the graph representation of a matrix. It is a type of\naggregation method, in which nodes in the fine grid are combined\nto form a single coarse grid node, and can efficiently operate on\nsparse matrices stored in the CSR format. In Ginkgo, we deploy\nPGM with deterministic aggregation of size 2, meaning PGM will\nalways try to aggregate two adjacent nodes.\nThe maximum number of multigrid levels is 11, with a min-\nimum of 64 rows in the coarsest matrix. We set the stopping\ncriterion as implicit relative residual norm reduction of 10−12 or\nmaximum of 700 iterations. The pre-/post-smoothing is weighted\nscalar Jacobi with a weight of 0.9, i.e., xi+1 = xi + 0.9D−1(b −Axi)\nwhere xi is the solution at iteration i, D is the diagonal matrix of\nA, and b is the right-hand side of the linear system being solved.\nThe same relaxation is used on the coarse grid problem, but with\nfour relaxation sweeps instead of one.\nWe refer to Algorithm 1 to describe the different precision\nusage in the implementation. The right-hand side b and solution\nx are from other solvers or user inputs, so their precision is\ndetermined outside of the multigrid cycle in Algorithm 1. Residual\nr (red) is based on the current level’s working precision. After\nrestricting r, we store the result g and zero initial guess e (brown)\nin the working precision of the next coarsest level. Precision\nconversion is part of the restriction and prolongation operations.\nAs mentioned previously, the matrix precision (blue) uses the\nsame as the current level’s working precision in ‘‘uniform levels’’.\nWe collect all data by performing 2 warmup applications,\nfollowed by 5 evaluation applications, of the AMG-preconditioned\nCG solver; the times of the 5 evaluation runs are averaged. The\nresults are summarized in Table 3. In (DP-SP), we see similar\ntrends as reported in [19], where the solver retains the same con-\nvergence behavior but has some improvement in performance. In\n(DP-SP-HP) and (DP-HP), some experiments fail to converge and\nsome matrices see delayed convergence. (DP-SP-HP) for cage13\nFig. 12. H100 AMG-preconditioned CG with Jacobi smoothers: beam problem\niteration counts as a function of the first level to switch to single or half\nprecision.\nand cage14 have the same number of iterations as (DP) and\n(DP-SP), so we can expect some performance benefit. In the\nbeam problem, the solver using half precision exhibits a conver-\ngence delay that cannot be compensated by faster computations.\nTable 1 shows that the cage13 and cage14 matrices are the\nonly ones whose non-zero values fall entirely in the range that\ncan be represented by half precision (which is approximately\n6×10−5–6.55×104).\n4.2.1. Effects of delayed precision changes in the AMG hierarchy\nFrom Table 3, the effect that the use of half precision can have\non convergence is clear. Using the beam problem as an example,\nwe plot heatmaps for the number of iterations and the total time\nwith different precision settings, delaying (in terms of distance\nfrom the finest level) the switch to single precision, half precision,\nor both. Each table has two axes:\n• x-axis (column): The level when we change to half precision,\nand\n• y-axis (row): The level when we change to single precision.\nEach entry (x, y) means the setting is the following:\n1. use double precision from the finest level (0) to the (x−1)\nlevel;\n2. use single precision from the x level to (y−1) level;\n3. use half precision from the y level to the coarsest level.\nThe last column (‘‘no’’) only changes to single precision, and the\nlast row (‘‘no’’) only changes to half precision (skipping single\nprecision). Also, the interaction of these two (no, no) is the purely\ndouble precision setting (DP). The coarsest solver always uses\nthe same precision as the last level. In terms of the previously-\nused configurations, (no, no) = (DP), (1, no) = (DP-SP), (no, 1) =\n(DP-HP), and (1, 2) = (DP-SP-HP). We currently do not consider\nconfigurations where we switch to higher precision for coarser\nlevels, though Ginkgo’s configuration flexibility, demonstrated\nby these tests, would make this simple to consider for future\nnumerical behavior analysis. As the iteration counts in Fig. 12\nmatch for all configurations in the two right-most columns, the\nfastest configurations are those in the upper right corner, which\nswitch to single precision sooner than those on the bottom right;\nsee\nFig. 13. Notably, these are the configurations that avoid or\nbarely use half precision: any speedup from using half precision\nearlier in the multigrid cycle is not enough to make up for the\nextra iterations required for convergence.\n4.2.2. Obstacles to convergence\nBased on these initial experiments, we identified three issues\naffecting the AMG-preconditioned CG convergence for the matri-\nces in Table 1 when using half precision in a ‘‘uniform precision\nlevels’’ configuration:\n286\n\n\nY.-H.M. Tsai, N. Beams and H. Anzt\nFuture Generation Computer Systems 149 (2023) 280–293\nTable 3\nPerformance of CG preconditioned with an AMG V-cycle, scalar Jacobi with uniform level precision configurations, on H100. The ‘‘packed half’’ SpMV implementation\nwas used.\nProblem\nGinkgo’s AMG (DP)\nGinkgo’s AMG (DP SP)\nGinkgo’s AMG (DP SP HP)\nGinkgo’s AMG (DP HP)\nres. norm\n#iter\nTime [ms]\nres. norm\n#iter\nTime [ms]\nres. norm\n#iter\nTime [ms]\nres. norm\n#iter\nTime [ms]\n2cubes_sphere\n6.56151e−09\n20\n14.0444\n6.56151e−09\n20\n13.676\nNaN\n700\n430.653\nNaN\n700\n443.579\ncage13\n3.68272e−10\n11\n15.206\n3.68271e−10\n11\n13.750\n4.23527e−10\n11\n13.892\n6.146e−10\n15\n17.779\ncage14\n3.41273e−10\n10\n33.494\n3.41273e−10\n10\n30.196\n3.72638e−10\n10\n29.433\n7.02332e−10\n13\n35.629\noffshore\n1594.67\n700\n679.344\n1325.81\n700\n608.188\nNaN\n700\n629.504\nNaN\n700\n586.523\nthermal2\n2.18369e−06\n349\n477.996\n2.37802e−06\n425\n536.342\nNaN\n700\n853.791\nNaN\n700\n803.019\ntmt_sym\n6.94616e−05\n359\n360.928\n7.41395e−05\n401\n372.497\nNaN\n700\n645.689\nNaN\n700\n599.100\nbeam-pw-sv0.1-o-3-l-3\n3.05273e−15\n44\n45.011\n3.05462e−15\n44\n42.417\n5.52572e−15\n86\n79.725\n7.26702e−15\n127\n119.101\nl-shape-const-o-3-l-7\n4.63552e−14\n160\n199.818\n4.78055e−14\n171\n191.208\n5.39855e−10\n700\n783.019\n3.59387e−07\n700\n758.162\nFig. 13. H100 AMG-preconditioned CG with Jacobi smoothers: beam problem\ntotal solve time as a function of the first level to switch to single or half\nprecision.\n• The matrix values are out of the range of half precision:\n2cubes_sphere and offshore.\nWhen applying the matrix with a zero vector, Inf ∗0 = NaN.\n• The application of the Jacobi smoother is out of the range of\nhalf-precision: thermal2 and tmt_sym.\nFor some diagonal values that are within the range of half\nprecision, the inverted values are not.\n• The residual is too small in coarse levels: L-shape.\nThis occurs for each level of multigrid where the residual\npassed to the coarser level after the restriction is zero in\nhalf precision due to underflow. This means the right-hand\nside of the coarser level becomes zero, giving the trivial so-\nlution on that level, and thus no contribution to the overall\ncorrection. L-shape does not encounter the NaN issue, but\nthe convergence is slower than the other settings. This issue\nmay materialize after a few iterations, rather than from the\nbeginning.\nWe utilize two mitigation strategies. First, row/column scaling\nas given in Algorithm 2.5 of [29], detailed in Algorithm 2, is\napplied to the 2cubes_sphere and offshore matrices to ensure\nall values are in the half precision range. To directly compare\nthe convergence of higher precisions and mixed precision, we\nperform the scaling on these matrices for all precision configu-\nrations, not just those involving half. We also add a scaling of\n1/2 to the restriction matrix in the aggregation method to protect\nagainst values exceeding the allowable range after the merging\nstep. With scaling in Algorithm 2, the (DP-SP-HP) configuration\non 2cubes_sphere can solve the problem without incurring NaNs.\nSecond, we employ the common mixed precision technique\nof decoupling the working vector precision from the matrix stor-\nage precision, analogous to computing IR residuals in a higher\nprecision to combat roundoff errors. Using a higher precision\nfor the vectors helps avoid zero residuals and Jacobi smoother\napplication overflow, but we retain most of the benefit from half\nprecision in SpMV as the matrix accounts for the bulk of the\nmemory movement. This also allows us to use a lower precision\nto store the finest level matrix. This will incur some extra setup\ncost and memory usage, as it requires copying and converting\nAlgorithm 2 Symmetry-preserving row and column equilibration\n(one iteration of Algorithm 2.5 in [29]).\n1: procedure Scale(A)\n2:\nfor all i in 0...n −1 do\n3:\nLet ri = max(abs(A(i, :)))1/2\n4:\nLet ci = max(abs(A(:, i)))1/2\n5:\nend for\n6:\nfor all row, col in A do\n7:\nA(row, col) = A(row,col)\nrrowccol\n8:\nend for\n9: end procedure\nthe original double precision matrix to store in half precision, but\nincreases the potential speedup for each application of the AMG\ncycle.\nReturning to Algorithm 1, we now also have ‘‘non-uniform\nlevels’’, where the matrix precision can be chosen separately from\nthe vector/working precision. The arithmetic operations in the\nresidual computation will always use the highest precision format\namong input, output, and matrix precisions and store the result\nin the output precision. In our configurations, the working vector\nprecision always uses more bits than the matrix precision when\nthey differ, so the precision used in the arithmetic operations is\nalways the working vector precision.\nWe consider the following non-uniform configurations:\n• The vector precision (also the arithmetic precision in our\ncases):\n1. (DP-SP): the first level’s vector uses double precision,\nbut the other levels’ vectors use single precision.\n2. (DP): all vectors use double precision.\n• The matrix precision:\n1. (SP): all matrices in multigrid use single precision.\n2. (HP): all matrices in multigrid use half precision.\n3. (DP-SP-HP): the first level’s matrix uses double pre-\ncision, the second level’s matrix uses single precision,\nand the other levels’ matrices use half precision.\nThus, we have 6 possible combinations from these options. In\nthe following discussion and figures, we use the notation (Work-\ning Precision, Matrix Precision) to represent a specific precision\nconfiguration. The uniform levels from previous experiments will\nnow be represented with the same precision shorthand repeated,\ne.g. (DP-HP, DP-HP).\n4.2.3. Additional configuration options to improve convergence\nEven in double precision, some of our test problems are dif-\nficult for the simple AMG configuration used in the previous\nexperiments; the offshore problem did not converge within the\niteration limit. To address this issue, as well as to compare mixed\n287\n\n\nY.-H.M. Tsai, N. Beams and H. Anzt\nFuture Generation Computer Systems 149 (2023) 280–293\nprecision performance in AMG cycles using more complicated\nsmoothers, we add two additional smoother configuration op-\ntions. The first is block Jacobi, with all non-block-size settings\nidentical to the scalar Jacobi smoother. We use 32 as the max-\nimum block size for block Jacobi smoothers (see [30] for details\nof Ginkgo’s block Jacobi implementation). Block Jacobi also did\nnot succeed in achieving convergence for the offshore problem,\nso we will omit the configurations based on scalar and block\nJacobi for this problem. We also test a configuration using an ℓ1-\nJacobi Chebyshev smoother, combining ℓ1-Jacobi from [31] with\nthe Chebyshev iteration, as in [32]. For this smoother, we set 2\niterations each for pre-/post-smoothing, as well as 2 iterations\nfor the coarse solver. As described in [32], when using ℓ1-Jacobi\nas the inner solver of Chebyshev on SPD (symmetric positive\ndefinite) matrices, all eigenvalues are in (0, 1]. For the two non-\nSPD matrices – cage13 and cage14 – we omit this smoother\nconfiguration in the following tests. This smoother is able to\nimprove the (DP) convergence of the offshore problem such that\nit converges before reaching the iteration limit of 700 iterations.\n4.3. Comparison of uniform and non-uniform precision configura-\ntions\nOn the H100, the convergence benefit of the non-uniform con-\nfiguration is clear. (DP-SP, HP) usually achieves the best speedup\nin total solve time in Figs. 14 to 16. For thermal2 and tmt_sym,\nwe need more iterations with single precision working vectors\n(Figs. 17 to 19), which affects performance. With double precision\nworking vectors on every level, we can keep the same iteration\nas the all-double multigrid setup, so (DP, HP) is the best for\nthe thermal2 and tmt_sym cases. (DP-SP, HP) in tmt_sym still\nshows some speedup compared to the full double settings be-\ncause the speedup per iteration can compensate for the increase\nin iterations. In Fig. 19, (DP-SP-HP, DP-SP-HP) and (DP-HP, DP-\nHP) configurations do not converge because ℓ1-Jacobi’s addition\nof the absolute values of all off-diagonal entries to the diagonal\nvalue causes overflow in half precision. We can get up to 1.35x\nspeedup in the preconditioned CG with Jacobi smoothers, up\nto 1.27x speedup with block Jacobi smoothers, and up to 1.24x\nspeedup with ℓ1-Jacobi Chebyshev smoothers. Note that if no\nspeedup was achieved by a particular configuration, it does not\nhave a bar for that problem. We embed additional information\ninto speedup figures: ‘‘NaN (×)’’ indicates the result residual\nnorm is NaN, ‘‘More iter (◦)’’ indicates the configuration requires\nmore iterations than the full double precision configuration, ‘‘No\nspeedup per iter (△)’’ indicates the time per iteration is not faster\nthan the full double precision configuration, and ‘‘Not converged\n(▽)’’ indicates the configuration does not converge.\nThe performance trends are very similar when executing on\nthe MI250X in Figs. 20 to 22. (DP-SP, HP) is usually the fastest\noption with Jacobi and blockJacobi smoothers. For cage14, (DP-\nSP, HP) performs the best, and can get up to 1.13x speedup\nwith Jacobi smoothers and up to 1.16x speedup with blockJacobi\nsmoothers. In Fig. 22, we get 1.25x speedup from (DP-SP, SP)\nfor the beam problem, but around 1.05x–1.10x speedup from\n(DP-SP, HP) or (DP, HP) in other cases. We expect less speedup\nfrom mixed precision than for H100 because of the smaller\nperformance\ndifferences\nin\nthe\nSpMV\nanalysis\nfrom\nFigs. 7 and 9.\nIn Figs. 23 to 25, we present the results obtained from running\non Intel’s PVC GPU (one tile). The performance trends are different\nfrom those we observed for the MI250X and the H100. We get\naround 1.15x speedup for the tmt_sym problem with (DP, DP-SP-\nHP) and for cage13 with (DP-SP, HP), while seeing almost 1.25x\nspeedup for the cage14 (DP-SP, HP) case with Jacobi and Block-\nJacobi smoothers. In ℓ1-Jacobi-Chebyshev, only (DP-SP, DP-SP) in\nFig. 14. Speedup in total solve time for CG with AMG V-cycle preconditioning,\nscalar Jacobi configuration. Results are for H100 with the packed half SpMV\nvariant. * denotes the matrix was scaled prior to solving.\nFig. 15. Speedup in total solve time for CG with AMG V-cycle preconditioning,\nblock Jacobi configuration. Results are for H100 with the packed half SpMV\nvariant. * denotes the matrix was scaled prior to solving.\nFig. 16. Speedup in total solve time for CG with AMG V-cycle preconditioning,\nℓ1-Jacobi-Chebyshev configuration. Results are for H100 with the packed half\nSpMV variant. * denotes the matrix was scaled prior to solving.\n288\n\n\nY.-H.M. Tsai, N. Beams and H. Anzt\nFuture Generation Computer Systems 149 (2023) 280–293\nFig. 17. Total iterations for CG with AMG V-cycle preconditioning, scalar Jacobi\nconfiguration. Results are shown for H100 with the packed half SpMV variant.\n* denotes the matrix was scaled prior to solving.\nFig. 18. Total iterations for CG with AMG V-cycle preconditioning, block Jacobi\nconfiguration. Results are shown for H100 with the packed half SpMV variant.\n* denotes the matrix was scaled prior to solving.\nFig. 19. Total iterations for CG with AMG V-cycle preconditioning, ℓ1-Jacobi-\nChebyshev configuration. Results are shown for H100 with the packed half SpMV\nvariant. * denotes the matrix was scaled prior to solving.\nbeam and (DP-SP, HP) for L-shape get close to 1.2x speedup. The\nother configurations does not show much speedup with mixed\nprecision on PVC, which is expected from the SpMV performance\nin Figs. 10 and 11.\nFinally, in Table 4, we collect the fastest (in total solve time)\nconfiguration for each matrix. The benefit of the non-uniform\nFig. 20. Speedup in total solve time for CG with AMG V-cycle preconditioning,\nscalar Jacobi configuration. Results are for one GCD of MI250X. * denotes the\nmatrix was scaled prior to solving.\nFig. 21. Speedup in total solve time for CG with AMG V-cycle preconditioning,\nblock Jacobi configuration. Results are for one GCD of MI250X. * denotes the\nmatrix was scaled prior to solving.\nFig. 22. Speedup in total solve time for CG with AMG V-cycle preconditioning,\nℓ1-Jacobi-Chebyshev configuration. Results are for one GCD of MI250X. * denotes\nthe matrix was scaled prior to solving.\n289\n\n\nY.-H.M. Tsai, N. Beams and H. Anzt\nFuture Generation Computer Systems 149 (2023) 280–293\nFig. 23. Speedup in total solve time for CG with AMG V-cycle preconditioning,\nscalar Jacobi configuration. Results are for one tile of PVC with the packed half\nvariant SpMV. * denotes the matrix was scaled prior to solving.\nFig. 24. Speedup in total solve time for CG with AMG V-cycle preconditioning,\nblock Jacobi configuration. Results are for one tile of PVC with the packed half\nvariant SpMV. * denotes the matrix was scaled prior to solving.\nFig. 25. Speedup in total solve time for CG with AMG V-cycle preconditioning,\nℓ1-Jacobi-Chebyshev configuration. Results are for one tile of PVC with the\npacked half variant SpMV. * denotes the matrix was scaled prior to solving.\nFig. 26. Speedup in total solve time for a standalone V-cycle AMG solver (scalar\nJacobi configuration) on H100, using packed half SpMV. * denotes the matrix was\nscaled prior to solving.\nmixed precision configurations is clear across all three architec-\ntures. Notably, for each matrix, the same smoother was the best\nchoice for all three architectures, and it was the scalar Jacobi\nsmoother for all but the offshore matrix. In terms of specific\nprecision configurations, MI250X and H100 always align, yet PVC\nis only the same for the cage14 and L-shape problems. The\n‘‘winning’’ configuration for every problem uses half precision\nmatrices on at least some levels, demonstrating that while the\nuse of half precision presents more challenges than single preci-\nsion, both in terms of convergence and efficient implementation,\nsolutions like matrix scaling, the use of higher working precision,\nand the packed SpMV variant can help address the issues that\nmay arise.\n4.4. Standalone multigrid solver\nWe also demonstrate Ginkgo’s AMG as a standalone solver on\nH100 in Figs. 26 and 27. From the success of the non-uniform\nsettings with half precision in the preconditioner experiments,\nwe use the Jacobi smoother with the same parameters as before,\nexcept with the absolute residual norm stopping criterion set\nto 1e−9. There are four cases that converge in double preci-\nsion settings (DP, DP): 2cubes_sphere, cage13, cage14, and beam.\nThe (DP-HP, DP-HP) configuration is generally poor for conver-\ngence. This is expected from the previously-described challenges\nof using half precision: the 2cubes_sphere and beam problems\ndo not converge, and cage13 and cage14 require more itera-\ntions than the other configurations. However, the non-uniform\nsettings with higher working precision preserve convergence. A\nnon-uniform configuration with half precision matrices performs\nthe best for all but the 2cubes_sphere problem. Overall, the stan-\ndalone multigrid shows higher speedup than the preconditioned\nCG because there is no additional double precision work outside\nof the multigrid. The cage14 problem attains up to 1.45x speedup\nwith (DP-SP, HP) configuration. For these parameters, using the\nsmoother in double precision only at the finest level is sufficient\nto preserve convergence and final accuracy without an outer CG\nsolver or iterative refinement.\n4.5. Other multigrid cycles\nThe standalone multigrid solver’s greater potential for speedup\nwith mixed precision is due to, in essence, reducing the ratio of\n290\n\n\nY.-H.M. Tsai, N. Beams and H. Anzt\nFuture Generation Computer Systems 149 (2023) 280–293\nTable 4\nThe fastest configuration for AMG-preconditioned CG, in terms of total solve time, across all smoothers and all\nmixed precision settings.\nSmoother\nH100\nMI250X(1GCD)\nPVC(1tile)\n2cubes_sphere*\nJacobi\n(DP-SP, HP)\n(DP-SP, HP)\n(DP-SP, DP-SP)\ncage13\nJacobi\n(DP-SP, HP)\n(DP-SP, HP)\n(DP, HP)\ncage14\nJacobi\n(DP-SP, HP)\n(DP-SP, HP)\n(DP-SP, HP)\noffshore*\nℓ1-Jacobi-Chebyshev\n(DP-SP, HP)\n(DP-SP, HP)\n(DP, HP)\nthermal2\nJacobi\n(DP, HP)\n(DP, HP)\n(DP, DP-SP-HP)\ntmt_sym\nJacobi\n(DP, HP)\n(DP, HP)\n(DP, DP-SP-HP)\nbeam-pw-sv0.1-o-3-l-3\nJacobi\n(DP-SP, HP)\n(DP-SP, HP)\n(DP-SP, DP-SP)\nl-shape-const-o-3-l-7\nJacobi\n(DP, HP)\n(DP, HP)\n(DP, HP)\nThe matrices marked with * were scaled to avoid exceeding the representation range of half precision.\nFig. 27. Total iterations for a standalone V-cycle AMG solve (scalar Jacobi\nconfiguration) on H100, using packed half SpMV. * denotes the matrix was scaled\nprior to solving.\ndouble precision use to lower precision use. To raise the ‘‘maxi-\nmum speedup ceiling’’ while still using AMG as a preconditioner\ninside a higher-precision linear solver like CG, the use of cycles\nother than the V-cycle is an option. As we move through the grid\nhierarchy from finest to coarsest, each level contributes less to\nthe total overall runtime and memory movement requirements\nof the AMG as the matrix sizes decrease. Thus, coarse levels\nwhich could achieve speedup ‘‘individually’’, i.e., relative to their\ncorresponding levels in a fully double precision multigrid, will\nhave a limited impact on the speedup of the full cycle. If we visit\nthe coarse levels more often during one full cycle of multigrid,\nresulting in more operations on the coarse matrices, we decrease\nthe ratio of double precision to lower precision for one full\niteration and thereby increase the potential benefit of the mixed\nprecision AMG. F- and W-cycles, for example, spend more time\non the coarse levels during each multigrid cycle. We try the same\nCG experiments on a small W-cycle with 4 levels. In Figs. 28 to 30,\nwe can get up to 1.6x speedup with Jacobi smoothers, 1.45x\nspeedup with block Jacobi smoothers, and 1.4x speedup with ℓ1-\nJacobi Chebyshev smoothers on H100. As with the V-cycle, when\nseparating the precision of working vectors and matrices on a\nlevel, we can use lower precision for the matrix on the finest level\nto enable even higher speedups.\n5. Conclusion\nIn this paper, we implement and evaluate a mixed precision\nalgebraic multigrid (AMG) method that allows the use of double\nprecision, single precision, and half precision. We demonstrate\nthat when using AMG as a preconditioner inside an iterative\nsolver, some linear systems allow for using single precision or\neven half precision on the coarser multigrid levels without im-\npacting the solution quality. We use techniques such as scaling\nFig. 28. Speedup in total solve time for CG with AMG W-cycle preconditioning,\nscalar Jacobi configuration with 4 levels. Results are for H100 with the packed\nhalf SpMV variant. * denotes the matrix was scaled prior to solving.\nFig. 29. Speedup in total solve time for CG with AMG W-cycle preconditioning,\nblock Jacobi configuration with 4 levels. Results are for H100 with the packed\nhalf SpMV variant. * denotes the matrix was scaled prior to solving.\nand decoupling the precision formats for the working vectors and\nmatrices on the same level to address convergence issues caused\nby the use of half precision. When using higher precision in work-\ning vectors, we found more mixed precision cases that maintain\nthe same number of iterations as the full double precision settings\n291\n\n\nY.-H.M. Tsai, N. Beams and H. Anzt\nFuture Generation Computer Systems 149 (2023) 280–293\nFig. 30. Speedup in total solve time for CG with AMG W-cycle preconditioning,\nℓ1-Jacobi-Chebyshev configuration with 4 levels. Results are for H100 with the\npacked half SpMV variant. * denotes the matrix was scaled prior to solving.\nwhile taking less time. Performance evaluations were completed\non AMD, Intel, and NVIDIA GPUs.\nCRediT authorship contribution statement\nYu-Hsiang Mike Tsai: Conceptualization, Methodology, Writ-\ning – original draft, Investigation, Formal analysis, Visualization.\nNatalie Beams: Formal analysis, Resources, Writing – review &\nediting. Hartwig Anzt: Conceptualization, Supervision, Writing –\nreview & editing.\nDeclaration of competing interest\nThe authors declare that they have no known competing finan-\ncial interests or personal relationships that could have appeared\nto influence the work reported in this paper.\nData availability\nData will be made available on request.\nAcknowledgments\nThe authors acknowledge support by the state of Baden-\nWürttemberg through bwHPC. This work was supported by the\n‘‘Impuls und Vernetzungsfond’’ of the Helmholtz Association un-\nder grant VH-NG-1241, and the US Exascale Computing Project\n(17-SC-20-SC), a collaborative effort of the U.S. Department of\nEnergy Office of Science and the National Nuclear Security Ad-\nministration. This research used resources of the Oak Ridge Lead-\nership Computing Facility at the Oak Ridge National Laboratory,\nwhich is supported by the Office of Science of the U.S. Depart-\nment of Energy under Contract No. DE-AC05-00OR22725, and\nresources of the Argonne Leadership Computing Facility, which\nis a DOE Office of Science User Facility supported under Contract\nDE-AC02-06CH11357.\nReferences\n[1] A. Brandt, Multi-level adaptive solutions to boundary-value problems,\nMath. Comp. 31 (138) (1977) 333–390, http://dx.doi.org/10.1090/S0025-\n5718-1977-0431719-X.\n[2] J.W. Ruge, K. Stüben, 4. Algebraic Multigrid, in: Multigrid Methods, in:\nFrontiers in Applied Mathematics, Society for Industrial and Applied Math-\nematics, 1987, pp. 73–130, http://dx.doi.org/10.1137/1.9781611971057.\nch4.\n[3] P. Wesseling, C.W. Oosterlee, Geometric multigrid with applications to\ncomputational fluid dynamics, Numerical Analysis 2000. Vol. VII: Partial\nDifferential Equations, J. Comput. Appl. Math. 128 (1) (2001) 311–334,\nhttp://dx.doi.org/10.1016/S0377-0427(00)00517-3.\n[4] R. Falgout, An introduction to algebraic multigrid, Comput. Sci. Eng. 8 (6)\n(2006) 24–33, http://dx.doi.org/10.1109/MCSE.2006.105.\n[5] N.J. Higham, T. Mary, Mixed precision algorithms in numerical lin-\near algebra, Acta Numer. 31 (2022) 347–414, http://dx.doi.org/10.1017/\nS0962492922000022.\n[6] A. Abdelfattah, H. Anzt, E.G. Boman, E. Carson, T. Cojean, J. Dongarra, A.\nFox, M. Gates, N.J. Higham, X.S. Li, et al., A survey of numerical linear\nalgebra methods utilizing mixed-precision arithmetic, Int. J. High Perform.\nComput. Appl. 35 (4) (2021) 344–369.\n[7] J.H. Wilkinson, Rounding Errors in Algebraic Processes, in: Prentice-Hall\nseries in automatic computation, Prentice-Hall, Englewood Cliffs, N.J, 1964.\n[8] C.B. Moler, Iterative Refinement in Floating Point, J. ACM 14 (2) (1967)\n316–321, http://dx.doi.org/10.1145/321386.321394.\n[9] J. Langou, J. Langou, P. Luszczek, J. Kurzak, A. Buttari, J. Dongarra, Exploiting\nthe Performance of 32 bit Floating Point Arithmetic in Obtaining 64 bit\nAccuracy (Revisiting Iterative Refinement for Linear Systems), in: SC ’06:\nProceedings of the 2006 ACM/IEEE Conference on Supercomputing, 2006,\np. 50, http://dx.doi.org/10.1109/SC.2006.30.\n[10] E. Carson, N.J. Higham, Accelerating the Solution of Linear Systems by\nIterative Refinement in Three Precisions, SIAM J. Sci. Comput. 40 (2) (2018)\nA817–A847, http://dx.doi.org/10.1137/17M1140819.\n[11] A. Haidar, S. Tomov, J. Dongarra, N.J. Higham, Harnessing GPU Tensor\nCores for Fast FP16 Arithmetic to Speed up Mixed-Precision Iterative\nRefinement Solvers, in: SC18: International Conference for High Perfor-\nmance Computing, Networking, Storage and Analysis, 2018, pp. 603–613,\nhttp://dx.doi.org/10.1109/SC.2018.00050.\n[12] J.A.\nLoe,\nC.A.\nGlusa,\nI.\nYamazaki,\nE.G.\nBoman,\nS.\nRajamanickam,\nExperimental Evaluation of Multiprecision Strategies for GMRES on GPUs.\n[13] D. Göddeke, R. Strzodka, Cyclic reduction tridiagonal solvers on GPUs\napplied to mixed-precision multigrid, IEEE Trans. Parallel Distrib. Syst. 22\n(1) (2010) 22–32.\n[14] Y. Sumiyoshi, A. Fujii, A. Nukada, T. Tanaka, Mixed-precision AMG method\nfor many core accelerators, in: Proceedings of the 21st European MPI Users’\nGroup Meeting, 2014, pp. 127–132.\n[15] K.L. Oo, A. Vogel, Accelerating Geometric Multigrid Preconditioning with\nHalf-Precision Arithmetic on GPUs, 2020, arXiv:2007.07539 [cs].\n[16] S.F. McCormick, J. Benzaken, R. Tamstorf, Algebraic error analysis for\nmixed-precision multigrid solvers, SIAM J. Sci. Comput. 43 (5) (2021)\nS392–S419.\n[17] R. Tamstorf, J. Benzaken, S.F. McCormick, Discretization-Error-Accurate\nMixed-Precision Multigrid Solvers, SIAM J. Sci. Comput. 43 (5) (2021)\nS420–S447, http://dx.doi.org/10.1137/20M1349230.\n[18] A. Buttari, M. Huber, P. Leleux, T. Mary, U. Rüde, B. Wohlmuth, Block\nlow-rank single precision coarse grid solvers for extreme scale multigrid\nmethods, Numer. Linear Algebra Appl. 29 (1) (2022) e2407, http://dx.doi.\norg/10.1002/nla.2407.\n[19] Y.-H.M. Tsai, N. Beams, H. Anzt, Mixed Precision Algebraic Multigrid\non GPUs, in: R. Wyrzykowski, J. Dongarra, E. Deelman, K. Karczewski\n(Eds.), Parallel Processing and Applied Mathematics, in: Lecture Notes\nin Computer Science, Springer International Publishing, Cham, 2023,\npp. 113–125, http://dx.doi.org/10.1007/978-3-031-30442-2_9.\n[20] M. Naumov, M. Arsaev, P. Castonguay, J. Cohen, J. Demouth, J. Eaton, S.\nLayton, N. Markovskiy, I. Reguly, N. Sakharnykh, et al., AmgX: A library for\nGPU accelerated algebraic multigrid and preconditioned iterative methods,\nSIAM J. Sci. Comput. 37 (5) (2015) S602–S626.\n[21] U.M. Yang, et al., BoomerAMG: a parallel algebraic multigrid solver and\npreconditioner, Appl. Numer. Math. 41 (1) (2002) 155–177.\n[22] NVIDIA, CUDA Best Practices, 2023, https://docs.nvidia.com/cuda/cuda-c-\nbest-practices-guide/index.html. (Accessed 14 July 2023).\n[23] AMD, rocBLAS Contributor’s Guide, 2023, https://rocblas.readthedocs.io/en/\nrocm-5.3.0/Contributors_Guide.html, 2023-07-14.\n292\n\n\nY.-H.M. Tsai, N. Beams and H. Anzt\nFuture Generation Computer Systems 149 (2023) 280–293\n[24] T. Cojean, Y.-H.M. Tsai, H. Anzt, Ginkgo—A math library designed for\nplatform portability, Parallel Comput. 111 (2022) 102902, http://dx.doi.org/\n10.1016/j.parco.2022.102902, URL https://www.sciencedirect.com/science/\narticle/pii/S0167819122000096.\n[25] H. Anzt, T. Cojean, G. Flegar, F. Göbel, T. Grützmacher, P. Nayak, T. Ribizel,\nY.M. Tsai, E.S. Quintana-Ortí, Ginkgo: A Modern Linear Operator Algebra\nFramework for High Performance Computing, ACM Trans. Math. Software\n48 (1) (2022) 2:1–2:33, http://dx.doi.org/10.1145/3480935.\n[26] T.A. Davis, Y. Hu, The University of Florida sparse matrix collection, ACM\nTrans. Math. Softw. 38 (1) (2011) 1–25.\n[27] R. Anderson, J. Andrej, A. Barker, J. Bramwell, J.-S. Camier, J.C.V. Dobrev, Y.\nDudouit, A. Fisher, T. Kolev, W. Pazner, M. Stowell, V. Tomov, I. Akkerman,\nJ. Dahm, D. Medina, S. Zampini, MFEM: A modular finite element methods\nlibrary, Comput. Math. Appl. 81 (2021) 42–74, http://dx.doi.org/10.1016/j.\ncamwa.2020.06.009.\n[28] MFEM: Modular finite element methods [Software], 2017, http://dx.doi.\norg/10.11578/dc.20171025.1248, mfem.org.\n[29] N.J. Higham, S. Pranesh, M. Zounon, Squeezing a Matrix into Half Precision,\nwith an Application to Solving Linear Systems, SIAM J. Sci. Comput. 41 (4)\n(2019) A2536–A2551, http://dx.doi.org/10.1137/18M1229511.\n[30] G. Flegar, H. Anzt, T. Cojean, E.S. Quintana-Ortí, Adaptive Precision Block-\nJacobi for High Performance Preconditioning in the Ginkgo Linear Algebra\nSoftware, ACM Trans. Math. Softw. 47 (2) (2021) http://dx.doi.org/10.1145/\n3441850.\n[31] A.H. Baker, R.D. Falgout, T.V. Kolev, U.M. Yang, Multigrid Smoothers for\nUltraparallel Computing, SIAM J. Sci. Comput. 33 (5) (2011) 2864–2887,\nhttp://dx.doi.org/10.1137/100798806, URL https://epubs.siam.org/doi/abs/\n10.1137/100798806.\n[32] A. El Haman Abdeselam, A. Napov, Y. Notay, Porting an aggregation-\nbased algebraic multigrid method to GPUs, ETNA - Electron. Trans. Numer.\nAnal. 55 (2022) 687–705, http://dx.doi.org/10.1553/etna_vol55s687, URL\nhttps://hw.oeaw.ac.at?arp=0x003da4b8.\nYu-Hsiang Mike Tsai is a Ph.D. student supervised by\nHartwig Anzt in Karlsruhe Institute of Technology.\nNatalie Beams is a research scientist at the Innovative\nComputing Laboratory at the University of Tennessee.\nHartwig Anzt is the Director of the Innovative Com-\nputing Lab and Professor in the Electrical Engineering\nand Computer Science Department of the University of\nTennessee. He also holds a Senior Research Scientist\nPosition at Steinbuch Centre for Computing at the\nKarlsruhe Institute of Technology.\n293","difficulty":"hard","domain":"Multi-Document QA","length":"short","question":"These two articles both focused on mix-precision acceleration, especially involving FP16, in multigrid preconditioners. What are the differences between them?","sub_domain":"Academic"}

Source: https://huggingface.co/datasets/zai-org/LongBench-v2

initial import

Posting: /agents

GET /api/v1/write?intent=publish&task_id=03cff3a8-48a6-51fd-8e5b-96519acd6dee&body={url_encoded_text}&agent_name={optional_name}&nonce={optional_random_id}