Pith. sign in

REVIEW 1 major objections 1 minor 34 references

Adaptive Sketching Based Construction of H2 Matrices on GPUs

T0 review · 1 major / 1 minor · reviewed 2026-08-15 · deepseek-v4-flash

Pith's one-line read A bottom-up adaptive sketching algorithm constructs H2 matrices in linear time and memory, and runs on GPUs with up to 1000x speedups over existing methods.

desk verdict Strong GPU engineering for H2 construction, but the adaptive accuracy guarantee is heuristic and needs proof or explicit framing as such. read the letter →

arxiv 2506.16759 v1 pith:VXGKX7JZ submitted 2025-06-20 cs.MS

classification cs.MS
keywords H2matrixadaptivesketchingbottom-upconstructionGPUbatchedkernelsinterpolativedecompositionnestedbasislinearcomplexityrandomizednumericalalgebra
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

The paper sets out to show that an $\mathcal{H}^2$ matrix — a memory-efficient hierarchical representation of a data-sparse matrix whose row and column bases are nested across levels — can be constructed from scratch in time and memory proportional to the matrix dimension $N$, given only a black-box sketching routine (a fast matrix–vector product with random vectors) and a routine to evaluate matrix entries. Earlier sketching-based constructions worked top-down and consumed many random vectors, especially for three-dimensional problems, which made them slow and memory-hungry. The proposed algorithm works bottom-up: at each level of the cluster tree it strips off the contributions of already-known dense and coupling blocks from the sketched samples, applies a row interpolative decomposition to the residual samples to select skeleton rows, and propagates the reduced samples upward, adding fresh random vectors only when a local QR-based test says the node has not yet converged. On a GPU the same algorithm runs with batched kernels that handle variable-size blocks, and the paper reports up to $13\times$ speedup over its own CPU version, up to $1000\times$ speedup over an existing GPU top-down construction, and up to $660\times$ speedup over an existing CPU sketching-based $\mathcal{H}$ construction. If these claims hold, $\mathcal{H}^2$ construction ceases to be the bottleneck it has been for large kernel matrices, integral equation operators, and low-rank updates inside sparse solvers.

What carries the argument

The load-bearing object is the bottom-up skeletonization sweep together with the sample-propagation identity $Y^l = W^l Y^{l-1}$, where $W^l$ is a block unit-triangular matrix assembled from the interpolation matrices of interpolative decompositions at level $l$. This identity means that once the leaf level has subtracted its dense blocks, every higher level's sample matrix is obtained by simple block operations on the previous level's samples, so the black-box operator $K_{\mathrm{blk}}$ is invoked only for the initial sketch and for adaptive additions. The second essential piece is the adaptive convergence test: for each node, a QR factorization of the local residual sample matrix $Y_{\mathrm{loc}}^\tau$ is computed, and the node is declared converged when the smallest diagonal entry of the triangular factor falls below an absolute threshold $\epsilon_{\mathrm{abs}}$; samples are added in blocks until every node at the level passes. On the GPU side, the machinery is batching: all operations are marshaled into a few kernel launches per level, with the block-sparse row matrix multiplication split into at most $C_{\mathrm{sp}}$ batched non-uniform matrix multiplications to avoid atomics.

What would settle it

Run the adaptive algorithm on a matrix whose off-diagonal blocks need noticeably higher ranks at coarser levels than at the leaves — for example a three-dimensional integral-equation kernel with strong admissibility and a target tolerance of $10^{-6}$ — and record, at every node, whether the local QR-based test reports convergence while the final measured relative error $\|K_{\mathrm{comp}} - K\|/\|K\|$ exceeds the target; the claimed adaptive guarantee falls if the local test passes everywhere yet the global error is out of tolerance.

Watch

Extended reading notes

Core claim

The central claim, stated on the paper's own terms, is that the nested-basis $\mathcal{H}^2$ format can be built in linear complexity by a bottom-up sweep that never re-applies the black-box sketching operator after the initial samples are taken. At the leaf level the algorithm forms $Y_{\mathrm{loc}}^\tau = Y^\tau - \sum_{b\in \mathcal{N}_\tau} D_{\tau,b}\Omega_b$, removing the directly evaluated dense inadmissible blocks so that only admissible interactions remain in the sample matrix; a row interpolative decomposition of $Y_{\mathrm{loc}}^\tau$ then picks skeleton indices $\tilde{I}_\tau$ and defines the leaf basis. At every higher level, the samples and random vectors of a node are assembled from its children and the residual after subtracting the coupling-block contributions is again decomposed, yielding the transfer matrices $E_{\tau_1}, E_{\tau_2}$ and the skeleton set, with coupling matrices $B_{\tau,b} = K(\tilde{I}_\tau,\tilde{I}_b)$ evaluated directly. The key identity $Y^l = W^l Y^{l-1}$, with $W^l$ block unit-triangular, is what lets the sweep propagate samples without further sketching. The paper asserts that for a fixed tolerance the sample count remains $O(1)$ — 256 in the reported experiments — making the whole construction an $O(N)$ algorithm when the rank is $r = O(\log 1/\epsilon)$.

Load-bearing premise

The accuracy guarantee rests on the assumption that a node's approximation is converged once the smallest diagonal entry of the QR factor of its local sample matrix falls below an absolute threshold, and that checking this locally at every node, level by level, keeps the global relative error of the assembled $\mathcal{H}^2$ matrix within the target tolerance; the paper does not prove this and notes that its simple error compensation does not fully account for errors accumulated while sweeping up the tree.

Editorial extensions

If this is right

  • For a fixed tolerance, construction costs $O(N)$ time and memory, so $\mathcal{H}^2$ representations of kernel and integral-equation matrices become practical at sizes where top-down sample counts would exhaust memory.
  • The sample count stops growing with $N$ — 256 in the reported runs instead of thousands to tens of thousands — which removes the dominant cost of the black-box sketching operator and is the main source of the reported speedups.
  • Because the construction is adaptive, users no longer need to guess the rank in advance; the algorithm adds sample blocks until the local convergence test is met.
  • The same CPU/GPU code base with batched variable-size kernels can be used for low-rank updates of existing $\mathcal{H}^2$ matrices, a step that appears in multifrontal factorizations and Schur-complement updates.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • If the local QR convergence test is sound, the same adaptive criterion could be applied to other hierarchical formats (HODLR, HSS, fully black-box bottom-up skeletonization), potentially giving them the same sample-count reduction.
  • A practitioner comparing methods may want to normalize by achieved accuracy: the baselines in the paper were constructed to $10^{-5}$ while the proposed algorithm targeted $10^{-6}$, so a time-to-accuracy comparison would put the speedup ratios in context.
  • The batched variable-size marshaling pattern is reusable beyond $\mathcal{H}^2$: the same kernel design should accelerate other data-sparse formats and sparse multifrontal front updates on GPUs.
  • The paper's admitted lack of a formal error analysis for the adaptive tolerance suggests a concrete research target: derive an a priori bound linking the per-node QR diagonal threshold to the global relative error, which would let users set $\epsilon_{\mathrm{abs}}$ with confidence.
Share X Bluesky LinkedIn Reddit HN

Signed reviews

No signed human review yet.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

1 major / 1 minor

Summary. The paper presents a bottom-up sketching-based algorithm for constructing H2 matrices with linear complexity, together with a GPU implementation based on batched dense linear algebra and batched entry evaluation. The algorithm requires a black-box matrix-vector product and an entry-evaluation routine, extends an earlier HSS construction method to strongly admissible H2 matrices, and includes an adaptive sampling mechanism that adds random vectors until a local QR-based convergence test is satisfied. Numerical experiments on covariance matrices, volume integral equation operators, and H2 update problems report speedups up to 13x over a CPU version, up to 1000x over a GPU top-down H2 construction in H2Opus, and up to 660x over ButterflyPACK, with memory scaling consistent with O(N). The paper claims to be the first GPU implementation of the bottom-up sketching-based H2 construction class.

Significance. If the correctness and complexity claims are fully substantiated, this is a valuable contribution: it demonstrates that bottom-up, strongly admissible H2 construction can be made adaptive, sample-efficient, and GPU-resident, and the batched implementation of non-uniform block-sparse operations is a nontrivial engineering achievement. The stated reduction in the number of random vectors relative to top-down methods is a concrete and useful advance. However, the central accuracy guarantee currently rests on an unproven local convergence test, and the headline speedup numbers are obtained from comparisons at different target tolerances; both issues need to be resolved before the performance claims can be accepted at face value.

major comments (1)
  1. [Section IV-B / Section V-C] The performance analysis in Section IV-B claims that the batched operations require only L=O(log N) kernel launches and that kernel-launch overhead is negligible, but the numerical breakdown in Fig. 7 shows that the convergence test (dominated by batched QR) takes up a significant fraction of GPU runtime, especially for smaller problems. The claim that the kernel-launch cost is negligible appears to be based on asymptotic counts rather than measured overhead, and the profiling data suggest that batched QR efficiency is a practical bottleneck. The authors should reconcile the asymptotic kernel-launch argument with the profiled time distribution, or restrict the claim accordingly.
minor comments (1)
  1. [Section I] The phrase 'It's also worth mentioning' should be changed to 'It is also worth mentioning' for formal style.

Circularity Check

0 steps flagged · score 0.0 of 10

No significant circularity: the H2 construction is an independent compression pipeline; deferred error analysis and a self-referential validation metric are caveats, not circular reductions.

full rationale

Algorithm 1 takes a black-box sketching operator K_blk(·) and an entry-evaluation function as inputs, and arrives at the H2 representation via interpolative decompositions of sketched sample blocks. No parameter is fitted to the reported accuracy, and no equation in the derivation is defined in terms of the output H2 matrix. The O(N) complexity claim is argued directly from the operation counts in Algorithm 1 (O(r^2 N) sampling and O(rN) entry evaluations under a fixed tolerance), not from a self-citation. The adaptive convergence test in Section III-B is a local QR-based criterion whose global error behavior is explicitly deferred to future work, and Section V-D admits that the error compensation 'does not fully account for the approximation errors as we sweep up the tree.' These are acknowledged limitations in the accuracy guarantee, but they are not circular: the algorithm does not define its reported error as the stopping criterion, and the experiments report measured errors below the target threshold. The relative error is measured against the provided sampler K_blk, and for covariance/IE experiments that sampler is itself an H2Opus H2 approximation. This makes the validation somewhat self-referential, but it is a stated experimental choice, not a load-bearing derivation; the construction is still an independent compression of the given oracle and is benchmarked against external packages (ButterflyPACK, STRUMPACK, H2Opus). Overall, no claimed result reduces by construction to its own inputs, so the paper is not circular.

Assumptions & free parameters 3 free parameters · 5 assumptions · 0 invented entities

The central claim rests on the black-box matvec oracle and cheap entry evaluation (domain assumptions), the standard randomized ID machinery (standard math), and an unproven adaptive convergence heuristic that the authors explicitly defer for future work. The supplied free parameters are hand-chosen algorithm settings, not fitted constants.

free parameters (3)
  • Leaf size = 64, 128, 256 in experiments
    Cluster tree leaf size set by hand in Section V.A; it trades memory against runtime and rank.
  • Initial sample block size = 256 (main runs), 32 (adaptive runs)
    Number of random vectors used for the first sketch in Algorithm 1; chosen by hand, not fitted to data.
  • Admissibility parameter eta = 0.7
    Strong admissibility threshold used for all reported runs in Section V.B; affects the matrix partitioning and sparsity constant.
assumptions (5)
  • domain assumption The target matrix admits a strongly-admissible H2 representation with block ranks r = O(log 1/epsilon) and bounded sparsity constant C_sp.
    Section III assumes a hierarchical block partitioning allowing low-rank compression is already computed; this is the standard low-rank assumption for H2 matrices.
  • domain assumption The black-box operator K_blk computes Y = K Omega in O(Nd) time and entry evaluation is cheap.
    Section III-A defines the input requirements; the linear-complexity claim depends on these oracle costs.
  • standard math Randomized sketching with Gaussian-like matrices provides a good column space with high probability, and column-pivoted QR gives a valid interpolative decomposition.
    Section II-B uses column ID via pivoted QR; the randomized range-finder guarantees are imported from the literature without restatement.
  • ad hoc to paper The QR-diagonal convergence test at each node guarantees the global H2 approximation error stays below the target tolerance.
    Section III-B introduces the test; the authors leave error analysis to future work and Section V-D notes the error compensation is not fully accounted for as samples move up the tree.
  • domain assumption The matrix is symmetric and real-valued.
    Section II-A states this simplification and notes extension to nonsymmetric or complex matrices is straightforward.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Adaptive Sketching Based Construction of H2 Matrices on GPUs." pith.science (2026). https://pith.science/paper/VXGKX7JZ

@misc{pith2026250616759,
  author       = {Pith},
  title        = {Pith review of: Adaptive Sketching Based Construction of H2 Matrices on GPUs},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/VXGKX7JZ}},
  note         = {Machine review of arXiv:2506.16759}
}
abstract

We develop a novel linear-complexity bottom-up sketching-based algorithm for constructing a $H^2$ matrix, and present its high performance GPU implementation. The construction algorithm requires both a black-box sketching operator and an entry evaluation function. The novelty of our GPU approach centers around the design and implementation of the above two operations in batched mode on GPU with accommodation for variable-size data structures in a batch. The batch algorithms minimize the number of kernel launches and maximize the GPU throughput. When applied to covariance matrices, volume IE matrices and $H^2$ update operations, our proposed GPU implementation achieves up to $13\times$ speedup over our CPU implementation, and up to $1000\times$ speedup over an existing GPU implementation of the top-down sketching-based algorithm from the H2Opus library. It also achieves a $660\times$ speedup over an existing sketching-based $H$ construction algorithm from the ButterflyPACK library. Our work represents the first GPU implementation of the class of bottom-up sketching-based $H^2$ construction algorithms.

Figures

Figures reproduced from arXiv: 2506.16759 by the authors.

Figure 2
Figure 2. The matrix tree for the hierarchical matrix in Fig. 1 [PITH_FULL_IMAGE:figures/full_fig_p002_2.png] view at source ↗
Figure 4
Figure 4. (a)-(b) Block partitioning of a hierarchical matrix for a 3D problem of size [PITH_FULL_IMAGE:figures/full_fig_p005_4.png] view at source ↗
Figure 5
Figure 5. The time of the CPU and GPU implementations of Algorithm 1 for the covariance and IE matrices as well as the [PITH_FULL_IMAGE:figures/full_fig_p008_5.png] view at source ↗
Figures from the paper (2 more)
Figure 6
Figure 6. Figure 6: (a) The memory of Algorithm 1 for the covariance and IE matrices. (b) The memory of Algorithm 1 and a few other [PITH_FULL_IMAGE:figures/full_fig_p008_6.png]
Figure 7
Figure 7. Figure 7: A breakdown of the construction time by percentage of time taken by each phase on (a) CPU and (b) GPU for varying [PITH_FULL_IMAGE:figures/full_fig_p009_7.png]

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

34 extracted references · 33 canonical work pages

  1. [1]

    A sparse matrix arithmetic based onH-matrices. part i: Introduction to-matrices,

    W. Hackbusch, “A sparse matrix arithmetic based onH-matrices. part i: Introduction to-matrices,”Computing, vol. 62, no. 2, pp. 89–108, 1999

  2. [2]

    Existence ofH-matrix approximants to the inverse fe-matrix of elliptic operators withl ∞-coefficients,

    M. Bebendorf and W. Hackbusch, “Existence ofH-matrix approximants to the inverse fe-matrix of elliptic operators withl ∞-coefficients,” Numerische Mathematik, vol. 95, pp. 1–28, 2003

  3. [3]

    Introduction to hierar- chical matrices with applications,

    S. B ¨orm, L. Grasedyck, and W. Hackbusch, “Introduction to hierar- chical matrices with applications,”Engineering analysis with boundary elements, vol. 27, no. 5, pp. 405–422, 2003

  4. [4]

    A fast direct solver for integral equations on locally refined boundary discretizations and its application to multiphase flow simulations,

    Y . Zhang, A. Gillman, and S. Veerapaneni, “A fast direct solver for integral equations on locally refined boundary discretizations and its application to multiphase flow simulations,”Advances in Computational Mathematics, vol. 48, no. 5, p. 63, 2022

  5. [5]

    Hierarchical interpolative factorization for self green’s function in 3d modified poisson-boltzmann equations,

    Y . Tu, Z. Xu, and H. Yang, “Hierarchical interpolative factorization for self green’s function in 3d modified poisson-boltzmann equations,” Communications on Applied Mathematics and Computation, pp. 1–26, 2024

  6. [6]

    Randomized sparse direct solvers,

    J. Xia, “Randomized sparse direct solvers,”SIAM Journal on Matrix Analysis and Applications, vol. 34, no. 1, pp. 197–227, 2013

  7. [7]

    A robust parallel preconditioner for indefinite systems using hierarchical matrices and randomized sampling,

    P. Ghysels, S. L. Xiaoye, C. Gorman, and F.-H. Rouet, “A robust parallel preconditioner for indefinite systems using hierarchical matrices and randomized sampling,” in2017 IEEE International Parallel and Distributed Processing Symposium (IPDPS). IEEE, 2017, pp. 897– 906

  8. [8]

    Scalable and memory-efficient kernel ridge regression,

    G. Ch ´avez, Y . Liu, P. Ghysels, X. S. Li, and E. Rebrova, “Scalable and memory-efficient kernel ridge regression,” in2020 IEEE International parallel and distributed processing symposium (IPDPS). IEEE, 2020, pp. 956–965

Show all 34 references
  1. [9]

    Fast direct methods for Gaussian processes,

    S. Ambikasaran, D. Foreman-Mackey, L. Greengard, D. W. Hogg, and M. O’Neil, “Fast direct methods for Gaussian processes,”IEEE transactions on pattern analysis and machine intelligence, vol. 38, no. 2, pp. 252–265, 2015

  2. [10]

    A superfast algorithm for Toeplitz systems of linear equations,

    S. Chandrasekaran, M. Gu, X. Sun, J. Xia, and J. Zhu, “A superfast algorithm for Toeplitz systems of linear equations,”SIAM Journal on Matrix Analysis and Applications, vol. 29, no. 4, pp. 1247–1266, 2008

  3. [11]

    AnO(NlogN)fast direct solver for partial hierarchically semi-separable matrices: with application to radial basis function interpolation,

    S. Ambikasaran and E. Darve, “AnO(NlogN)fast direct solver for partial hierarchically semi-separable matrices: with application to radial basis function interpolation,”Journal of Scientific Computing, vol. 57, pp. 477–501, 2013

  4. [12]

    A fast solver for HSS representations via sparse matrices,

    S. Chandrasekaran, P. Dewilde, M. Gu, W. Lyons, and T. Pals, “A fast solver for HSS representations via sparse matrices,”SIAM Journal on Matrix Analysis and Applications, vol. 29, no. 1, pp. 67–81, 2007

  5. [13]

    A direct solver with O(N)complexity for integral equations on one-dimensional domains,

    A. Gillman, P. M. Young, and P.-G. Martinsson, “A direct solver with O(N)complexity for integral equations on one-dimensional domains,” Frontiers of Mathematics in China, vol. 7, pp. 217–247, 2012

  6. [14]

    The inverse fast multipole method: using a fast approximate direct solver as a preconditioner for dense linear systems,

    P. Coulier, H. Pouransari, and E. Darve, “The inverse fast multipole method: using a fast approximate direct solver as a preconditioner for dense linear systems,”SIAM Journal on Scientific Computing, vol. 39, no. 3, pp. A761–A796, 2017

  7. [15]

    Hierarchical interpolative factorization for elliptic operators: integral equations,

    K. L. Ho and L. Ying, “Hierarchical interpolative factorization for elliptic operators: integral equations,”Communications on Pure and Applied Mathematics, vol. 69, no. 7, pp. 1314–1353, 2016

  8. [16]

    H2pack: High-performance h2 matrix package for kernel matrices using the proxy point method,

    H. Huang, X. Xing, and E. Chow, “H2pack: High-performance h2 matrix package for kernel matrices using the proxy point method,”ACM Trans. Math. Softw., vol. 47, no. 1, dec 2020

  9. [17]

    H2Opus: a distributed-memory multi-GPU software package for non- local operators,

    S. Zampini, W. Boukaram, G. Turkiyyah, O. Knio, and D. Keyes, “H2Opus: a distributed-memory multi-GPU software package for non- local operators,”Advances in Computational Mathematics, vol. 48, no. 3, p. 31, 2022

  10. [18]

    DistributedH 2-matrices for boundary element methods,

    S. B ¨orm, “DistributedH 2-matrices for boundary element methods,” ACM Trans. Math. Softw., vol. 49, no. 2, jun 2023

  11. [19]

    ASKIT: An efficient, parallel library for high-dimensional kernel summations,

    W. B. March, B. Xiao, C. D. Yu, and G. Biros, “ASKIT: An efficient, parallel library for high-dimensional kernel summations,”SIAM Journal on Scientific Computing, vol. 38, no. 5, pp. S720–S749, 2016

  12. [20]

    Geometry-oblivious FMM for compressing dense SPD matrices,

    C. D. Yu, J. Levitt, S. Reiz, and G. Biros, “Geometry-oblivious FMM for compressing dense SPD matrices,” inProceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, 2017, pp. 1–14

  13. [21]

    Algorithmic patterns for H-matrices on many-core proces- sors,

    P. Zaspel, “Algorithmic patterns for H-matrices on many-core proces- sors,”Journal of Scientific Computing, vol. 78, no. 2, pp. 1174–1206, 2019

  14. [22]

    Fast construction of hierarchical matrix representation from matrix-vector multiplication,

    L. Lin, J. Lu, and L. Ying, “Fast construction of hierarchical matrix representation from matrix-vector multiplication,”Journal of Computa- tional Physics, vol. 230, no. 10, pp. 4071–4087, 2011

  15. [23]

    Randomized compression of rank- structured matrices accelerated with graph coloring,

    J. Levitt and P.-G. Martinsson, “Randomized compression of rank- structured matrices accelerated with graph coloring,”arXiv preprint arXiv:2205.03406, 2022

  16. [24]

    Randomized strong recursive skele- tonization: Simultaneous compression and factorization ofH-matrices in the black-box setting,

    A. Yesypenko and P. Martinsson, “Randomized strong recursive skele- tonization: Simultaneous compression and factorization ofH-matrices in the black-box setting,”arXiv preprint arXiv:2311.01451, 2023

  17. [25]

    A recursive skeletoniza- tion factorization based on strong admissibility,

    V . Minden, K. L. Ho, A. Damle, and L. Ying, “A recursive skeletoniza- tion factorization based on strong admissibility,”Multiscale Modeling & Simulation, vol. 15, no. 2, pp. 768–796, 2017

  18. [26]

    Direct solution of generalH 2 -matrices with controlled accuracy and concurrent change of cluster bases for elec- tromagnetic analysis,

    M. Ma and D. Jiao, “Direct solution of generalH 2 -matrices with controlled accuracy and concurrent change of cluster bases for elec- tromagnetic analysis,”IEEE Transactions on Microwave Theory and Techniques, vol. 67, no. 6, pp. 2114–2127, 2019

  19. [27]

    AnO(N) distributed-memory parallel direct solver for planar integral equations,

    T. Liang, C. Chen, P.-G. Martinsson, and G. Biros, “AnO(N) distributed-memory parallel direct solver for planar integral equations,” in2024 IEEE International Parallel and Distributed Processing Sympo- sium (IPDPS). IEEE, 2024, pp. 440–452

  20. [28]

    Scalable linear time dense direct solver for 3-d problems without trailing sub-matrix dependencies,

    Q. Ma, S. Deshmukh, and R. Yokota, “Scalable linear time dense direct solver for 3-d problems without trailing sub-matrix dependencies,” inSC22: International Conference for High Performance Computing, Networking, Storage and Analysis. IEEE, 2022, pp. 1–12

  21. [29]

    A fast randomized algorithm for computing a hier- archically semiseparable representation of a matrix,

    P.-G. Martinsson, “A fast randomized algorithm for computing a hier- archically semiseparable representation of a matrix,”SIAM Journal on Matrix Analysis and Applications, vol. 32, no. 4, pp. 1251–1274, 2011

  22. [30]

    Thrust: A productivity-oriented library for CUDA,

    N. Bell and J. Hoberock, “Thrust: A productivity-oriented library for CUDA,” inGPU computing gems Jade edition. Elsevier, 2011, pp. 359–371

  23. [31]

    KBLAS: An optimized library for dense matrix-vector multiplication on GPU accelerators,

    A. Abdelfattah, D. Keyes, and H. Ltaief, “KBLAS: An optimized library for dense matrix-vector multiplication on GPU accelerators,” ACM Transactions on Mathematical Software (TOMS), vol. 42, no. 3, pp. 1–31, 2016

  24. [32]

    Frame- work for batched and GPU-resident factorization algorithms to block Householder transformations,

    A. Haidar, T. Dong, S. Tomov, P. Luszczek, and J. Dongarra, “Frame- work for batched and GPU-resident factorization algorithms to block Householder transformations,” inISC High Performance, Springer. Frankfurt, Germany: Springer, 07-2015 2015

  25. [33]

    ButterflyPACK,

    Y . Liu and USDOE, “ButterflyPACK,” 11 2018. [Online]. Available: https://www.osti.gov//servlets/purl/1564244

  26. [34]

    STRUMPACK – STRUctured Matrices PACKage,

    P. Ghysels and USDOE, “STRUMPACK – STRUctured Matrices PACKage,” 12 2014. [Online]. Available: https://www.osti.gov/biblio/ 1328126

Pith tools

Reviewed August 15, 2026 · model on record in the stance chip above.