Largely due to their increased native capacity for numerical intensity and power efficiency, reduced-precision floating-point computing resources, primarily used in artificial intelligence (AI) applications, have expanded at a greater rate than their higher-precision relatives. This has led to various efforts focused upon leveraging plentiful reduced-precision hardware to mimic higher-precision mathematical calculations. This paper studies a specific use case, namely the use of bfloat16 (BF16) Tensor Cores found on modern GPUs in service of single precision (FP32) matrix multiply operations. Given that BF16 and FP32 share the same dynamic range, the option to accumulate BF16 operations into FP32 accumulators (at full-speed), and additional BF16 arithmetic characteristics specific to the Blackwell GPU architecture, such as integrated scaling hardware, such emulation is highly motivated. This paper examines the performance, efficiency, power, and numerical characteristics of FP32 matrix multiplication via BF16-based emulation and demonstrates how it exceeds numerical and performance characteristics of native FP32 for scientific applications. We also discuss a full library-ready implementation that correctly deals with denormals.

Many High-Performance Computing (HPC) libraries rely on decision trees to select the best kernel hyperparameters at runtime, depending on the input and environment. However, finding optimized configurations for each input and environment is challenging and requires significant manual effort and computational resources. This article presents MLKAPS, a tool that automates this task using machine learning and adaptive sampling techniques. MLKAPS generates decision trees that tune HPC kernels' design parameters to achieve efficient performance for any user input. MLKAPS scales to large input and design spaces, outperforming similar state-of-the-art auto-tuning tools in tuning time and mean speedup. We demonstrate the benefits of MLKAPS on the highly optimized Intel (R) MKL dgetrf LU kernel and show that MLKAPS finds blindspots in the manual tuning of HPC experts. It improves over 85% of the inputs with a geomean speedup of x1.31. On the Intel (R) MKL dgeqrf QR kernel, MLKAPS improves performance on 85% of the inputs with a geomean speedup of x1.18.

HPL-MxP has become a widely accepted benchmark for assessing high-performance computing (HPC) systems' capabilities for Artificial Intelligence (AI) workloads. However, the benchmark's representativeness of real-world HPC and AI workloads is unclear. In this paper, we discuss the HPL-MxP benchmark from a numerical perspective and propose new rules and data generation for numerically meaningful comparisons. We present experiments showing that the current HPL-MxP benchmark cannot be considered as a numerically relevant benchmark to prove superiority of a new format or algorithm. We propose to better specify these requirements for numerical formats to produce comparable performance numbers, and suggest new input data generation to make it numerically relevant. We validate our proposal on Int8, Int4, and BF16 implementations to demonstrate the numerical significance of the benchmark using our new generator.

Numerical exceptions, which may be caused by overflow, operations like division by 0 or sqrt(-1), or convergence failures, are unavoidable in many cases, in particular when software is used on unforeseen and difficult inputs. As more aspects of society become automated e.g., self-driving cars, health monitors, and cyber-physical systems more generally, it is becoming increasingly important to design software that is resilient to exceptions, and that responds to them in a consistent way. Consistency is needed to allow users to build higher-level software that is also resilient and consistent (and so on recursively). In this paper we explore the design space of consistent exception handling for the widely used BLAS and LAPACK linear algebra libraries, pointing out a variety of instances of inconsistent exception handling in the current versions, and propose a new design that balances consistency, complexity, ease of use, and performance. Some compromises are needed, because there are preexisting inconsistencies that are outside our control, including in or between existing vendor BLAS implementations, different programming languages, and even compilers for the same programming language. And user requests from our surveys are quite diverse. We also propose our design as a possible model for other numerical software, and welcome comments on our design choices.

Fused Multiply-Add (FMA) functional units constitute a fundamental hardware component to train Deep Neural Networks (DNNs). Its silicon area grows quadratically with the mantissa bit count of the computer number format, which has motivated the adoption of the BrainFloat16 format (BF16). BF16 features 1 sign, 8 exponent and 7 explicit mantissa bits. Some approaches to train DNNs achieve significant performance benefits by using the BF16 format. However, these approaches must combine BF16 with the standard IEEE 754 Floating-Point 32-bit (FP32) format to achieve state-of-the-art training accuracy, which limits the impact of adopting BF16. This article proposes the first approach able to train complex DNNs entirely using the BF16 format. We propose a new class of FMA operators, $\mathrm{FMA}^{\mathrm {bf}16}_{\mathrm{n}\_\mathrm{m}}$ , that entirely rely on BF16 FMA hardware instructions and deliver the same accuracy as FP32. $\mathrm{FMA}^{\mathrm {bf}16}_{\mathrm{n}\_\mathrm{m}}$ operators achieve performance improvements within the 1.28-1.35× range on ResNet101 with respect to FP32. $\mathrm{FMA}^{\mathrm {bf}16}_{\mathrm{n}\_\mathrm{m}}$ enables training complex DNNs on simple low-end hardware devices without requiring expensive FP32 FMA functional units.

Mixed-precision (MP) arithmetic combining both single- and half-precision operands has been successfully applied to train deep neural networks. Despite its advantages in terms of reducing the need for key resources like memory bandwidth or register file size, it has a limited capacity for diminishing further computing costs, as it requires 32-bits to represent its output. On the other hand, full half-precision arithmetic fails to deliver state-of-the-art training accuracy. We design a binary tool SERP based on Intel Pin which allows us to characterize and analyze computer arithmetic usage in machine learning frameworks (Pytorch, Caffe, Tensorflow) and to emulate different floating point formats. Based on empirical observations about precision needs on representative deep neural networks, this paper proposes a seamless approach to dynamically adapt floating point arithmetic. Our dynamically adaptive methodology enables the use of full half-precision arithmetic for up to 96.4% of the computations when training state-of-the-art neural networks; while delivering comparable accuracy to 32-bit floating point arithmetic. Microarchitectural simulations indicate that our Dynamic approach accelerates training deep convolutional and recurrent networks with respect to FP32 by 1.39 × and 1.26 ×, respectively.

Deep learning (DL) is one of the most prominent branches of machine learning. Due to the immense computational cost of DL workloads, industry and academia have developed DL libraries with highly-specialized kernels for each workload/architecture, leading to numerous, complex code-bases that strive for performance, yet they are hard to maintain and do not generalize. In this work, we introduce the batch-reduce GEMM kernel and show how the most popular DL algorithms can be formulated with this kernel as the basic building-block. Consequently, the DL library-development degenerates to mere (potentially automatic) tuning of loops around this sole optimized kernel. By exploiting our new kernel we implement Recurrent Neural Networks, Convolution Neural Networks and Multilayer Perceptron training and inference primitives in just 3K lines of high-level code. Our primitives outperform vendor-optimized libraries on multi-node CPU clusters, and we also provide proof-of-concept CNN kernels targeting GPUs. Finally, we demonstrate that the batch-reduce GEMM kernel within a tensor compiler yields high-performance CNN primitives, further amplifying the viability of our approach.

Convolution layers are prevalent in many classes of deep neural networks, including Convolutional Neural Networks (CNNs) which provide state-of-the-art results for tasks like image recognition, neural machine translation and speech recognition. The computationally expensive nature of a convolution operation has led to the proliferation of implementations including matrix-matrix multiplication formulation, and direct convolution primarily targeting GPUs. In this paper, we introduce direct convolution kernels for x86 architectures, in particular for Xeon and Xeon Phi systems, which are implemented via a dynamic compilation approach. Our JIT-based implementation shows close to theoretical peak performance, depending on the setting and the CPU architecture at hand. We additionally demonstrate how these JIT-optimized kernels can be integrated into a light-weight multi-node graph execution model. This illustrates that single- and multi-node runs yield high efficiencies and high image-throughputs when executing state-of-the-art image recognition tasks on CPUs.

High order methods for the solution of PDEs expose a trade-off between computational cost and accuracy on a per degree of freedom basis. In many cases, the cost increases due to higher arithmetic intensity while affecting data movement minimally. As architectures tend towards wider vector instructions and expect higher arithmetic intensities, the best order for a particular simulation may change.

Many modern highly scalable scientific simulations packages rely on small matrix multiplications as their main computational engine. Math libraries or compilers are unlikely to provide the best possible kernel performance. To address this issue, we present a library which provides high performance small matrix multiplications targeting all recent x86 vector instruction set extensions up to Intel AVX-512. Our evaluation proves that speed-ups of more than 10 x are possible depending on the CPU and application. These speed-ups are achieved by a combination of several novel technologies. We use a code generator which has a built-in architectural model to create code which runs well without requiring an auto-tuning phase. Since such code is very specialized we leverage just-in-time compilation to only build the required kernel variant at runtime. To keep ease-of-use, overhead, and kernel management under control we accompany our library with a BLAS-compliant frontend which features a multi-level code-cache hierarchy.

We dispel with "street wisdom" regarding the practical implementation of Strassen's algorithm for matrix-matrix multiplication (DGEMM). Conventional wisdom: it is only practical for very large matrices. Our implementation is practical for small matrices. Conventional wisdom: the matrices being multiplied should be relatively square. Our implementation is practical for rank-k updates, where k is relatively small (a shape of importance for libraries like LAPACK). Conventional wisdom: it inherently requires substantial workspace. Our implementation requires no workspace beyond buffers already incorporated into conventional high-performance DGEMM implementations. Conventional wisdom: a Strassen DGEMM interface must pass in workspace. Our implementation requires no such workspace and can be plug-compatible with the standard DGEMM interface. Conventional wisdom: it is hard to demonstrate speedup on multi-core architectures. Our implementation demonstrates speedup over conventional DGEMM even on an Intel(R) Xeon Phi(TM) coprocessor utilizing 240 threads. We show how a distributed memory matrix-matrix multiplication also benefits from these advances.

In this work we present a library, LIBXSMM, that provides a high performance implementation of small sparse and dense matrix multiplications on latest Intel architectures. Such operations are important building blocks in modern scientific applications and general math libraries are normally tuned for all dimensions being large. LIBXSMM follows a matrix multiplication code generation approach specifically matching the applications’ needs. By providing several interfaces, the replacement of BLAS calls is simple and straightforward. We show that depending on the application’s characteristics, LIBXSMM can either leverage the entire DRAM bandwidth or reaches close to the processor’s computational peak performance. Our performance results of CP2K and SeisSol therefore demonstrate that using LIBXSMM as a highly-efficient computational backend, leads to speed-ups of greater than two compared to compiler generated inlined code or calling highly-optimized vendor math libraries.

Observations of more than 300 Sun-like field stars carried out at Fairborn Observatory since 1993 now include 168 observed for 10 years or longer. This project, a successor to previous work at Lowell Observatory and Cloudcroft Observatory, beginning in 1955 and continuing for nearly half a century, demonstrates finally that variability at the Sun’s low level of total irradiance variation can be detected in stars. By also including Ca II H&K observations from the Mount Wilson and Lowell observatories, we have discovered how the patterns of photometric variability are related to chromospheric activity. We discuss the limitations of detectability imposed by comparison star variability and as an example of a star close to the detection limit we show the evidence for cyclic variability in the solar twin 18 Sco. 1. Historical perspective Total solar irradiance (TSI) variability, revealed finally in 1980 by the Solar Maximum Mission spacecraft after a nearly century-long quest, is now a well-established component of solar lore. From three full solar cycles of observation recorded by a series of spacecraft, we learned that the Sun’s output rises and falls by about 0.07% in step with the sunspot number. Transient dips several times larger can occur as large sunspot groups cross the solar disk (Willson et al. 1981, 1984). The search for solar variation spawned a half century of effort by Charles Abbott, whose famed Smithsonian Astrophysical Observatory solar constant program finally drew to an inconclusive close in 1955. As it was ending, Lowell Observatory, taking advantage of newly developed techniques of astronomical photoelectric photometry, began measurements to look for variability in the sunlight reflected from Uranus and Neptune. The results, documented in a series of papers optimistically titled “The Sun as a Variable Star” (Johnson & Iriarte 1959; Serkowski 1961; Jerzykiewicz & Serkowski 1966) were just as inconclusive as Abbot’s had been due to intrinsic and seasonal planetary variability. The Lowell project, however, produced an unexpected dividend related to the solar variability question. Parallel B, V observations of 16 Fand G-type main sequence “10year standards” had been made from 1955 to 1966 along with the planets. M. Jerzykiewicz and K. Serkowski, authors of the final report, recognized the significance of the collateral observations of the 10 year standards:

We report the detection of a giant planet in a 6.4950 day orbit around the 1.7 Msun subgiant HD102956. The planet has a semimajor axis a = 0.081 AU and minimum mass Msini = 0.96 Mjup. HD102956 is the most massive star known to harbor a hot Jupiter, and its planet is only the third known to orbit within 0.6 AU of a star more massive than 1.5 Msun. Based on our sample of 137 subgiants with M* > 1.45 Msun we find that 0.5-2.3% of A-type stars harbor a close-in planet (a < 0.1 AU) with Msini > 1 Mjup, consistent with hot-Jupiter occurrence for Sun-like stars. Thus, the paucity of planets with 0.1 < a/AU < 1.0 around intermediate-mass stars may be an exaggerated version of the "period valley" that is characteristic of planets around Sun-like stars.