ParU is a sparse direct solver based on a novel parallel multifrontal algorithm. In contrast with UMFPACK, which draws its performance from highly customized data structures and parallel BLAS libraries, the goal of ParU is to maximize the sources of parallelism it can exploit through OpenMP tasking and nested parallelism, i.e., node -level parallelism within its rectangular frontal matrices, and tree -level parallelism across the assembly tree. ParU relies on the analysis phase of UMFPACK, with the addition of strategies to exploit tree parallelism. The factorization phase is based on a variant of the multifrontal method, where the fronts are unsymmetric and rectangular shaped. Similarly to UMFPACK, it exploits sparse pivots within the fronts, but simplifies the data structures to allow for coarser-grained node parallelism (mostly through large GEMM (matrix multiplication) and TRSM (triangular solve)). Experimental results show that ParU is competitive in terms of performance with state-of-the-art libraries as UMFPACK, SuperLU, and MUMPS.
We present a mixed precision approach for improving Tucker-based lossy compression of multiway data. Current Tucker decomposition algorithms use a uniform floating-point precision and obtain compression by minimizing the multilinear ranks of the representation subject to the target error tolerance. Our algorithm automatically determines which parameters of the representation can be stored in lower precision to improve the compression ratio without exceeding the error tolerance. Crucially, we show that the opportunities to reduce precision are not uniform across parameters, and thus using multiple precisions allows for smoothly transitioning across tolerances between the precisions of standard formats. Theoretical error bounds and numerical experiments validate the approach, and we obtain up to \(2\times\) and \(4\times\) compression improvement factors over initial Tucker approximations of single and double precision real-world numerical simulation data, respectively.
The Ogita–Aishima iterative refinement algorithm for SVD is based on matrix–matrix multiplication, making it an attractive approach for GPUs. We extend the algorithm to complex arithmetic and add cluster handling. The latter is based on Wedin's perturbation theory and a per-pair version of Ogita–Aishima's convergence analysis. Earlier implementations require well-separated singular values; the proposed cluster handling removes this restriction. Our implementation refines an FP32-quality SVD computed with cuSolver's cusolverDnXgesvdp to FP64-quality, while matching the accuracy of the corresponding FP64 routine on all test problems—twelve synthetic singular-value distributions and 38 matrices from the SuiteSparse collection. The performance depends on the spectrum and the FP32-to-FP64 ratio. Under Ozaki-II FP64 emulation on an NVIDIA RTX PRO 6000 Blackwell workstation GPU, which has a \(64\) : \(1\) FP32-to-FP64 ratio, the refinement is up to \(11.1\times\) faster than cusolverDnXgesvdp on spectra whose adjacent gaps are resolved by the FP32 SVD, and at most \(1.9\times\) slower on spectra for which they are not. The advantage is at most \(1.5\times\) on data-center GPUs with a \(2\) : \(1\) FP32-to-FP64 ratio, such as the NVIDIA B200. These results show the algorithm's potential on hardware with strong mixed-precision capabilities.
Precision emulation algorithms, such as the Ozaki schemes, can leverage low-precision matrix multiply–accumulate to approximate the product \(AB\) , where \(A\) and \(B\) are high-precision floating-point matrices, via fixed-point arithmetic. A key difficulty of these approaches is that achieving a target floating-point precision \(p\) requires a fixed-point bit width \(w=p+\star\) , where the required number of extra bits \(\star\) depends strongly on the matrix coefficients. The exponent span capacity (ESC), recently introduced in the literature, is an upper bound \(\overline{e}\geq\star\) . For square matrices of order \(n\) , the ESC can be estimated with \(O(n^{2}k)\) operations, where the parameter \(k\) can be increased to improve the sharpness of \(\overline{e}\) . We show that this bound can be very pessimistic, and that for the estimator to be useful, \(k\) must scale with \(n\) , effectively incurring an \(O(n^{3})\) overhead. Building on a rigorous error analysis, we then derive a new formula for \(\star\) , which corroborates the previously proposed formula up to a few bits. We propose an alternative, sharper bound, which can be computed with \(O(n^{2}k)\) operations, where \(k\) becomes independent of \(n\) . Our numerical experiments demonstrate that the current approach can substantially overestimate \(\star\) in many cases, whilst our new approach is almost perfectly sharp even for moderate values of \(k\) . We present strategies for handling matrices with zero coefficients and give a refined formulation for \(\star\) , in which different fixed-point bit widths are used for \(A\) and \(B\) . Overall, this work provides rigorous theoretical accuracy guarantees for fixed-point emulation and a practical algorithm to inexpensively yet sharply determine the required fixed-point bit width.
Physics-Informed Machine Learning (PIML) integrates physical laws into neural network training to solve complex problems governed by differential equations. A critical component in PIML is computing derivatives of models to ensure compliance with the governing equations. A common approach is to use automatic differentiation (AD) methods available in machine learning frameworks such as PyTorch, which is accurate but not always the most efficient approach. An alternative AD method is demonstrated here in the context of physics-regularized genetic programming with symbolic regression (PR-GPSR). This work introduces Hypercomplex-based Automatic Differentiation (HYPAD), which utilizes the Order Truncated Imaginary (OTI) library to efficiently compute second- and fourth-order derivatives in PR-GPSR applications via BINGO, an open-source symbolic regression library. Numerical experiments demonstrate that HYPAD achieves up to a \(123\times\) runtime improvement over PyTorch implementations, with both approaches recovering the correct solutions in the demonstrated PR-GPSR problems. HYPAD evaluates functions by propagating OTI hypercomplex numbers through a command array representation, enabling computation of arbitrary-order derivatives. The approach may also be applicable beyond symbolic regression, motivating further investigation in other PIML and AD applications.
