Search

Vector spherical harmonics on the unit sphere of ℝ 3 have broad applications in geophysics, quantum mechanics, and astrophysics. In the representation of a tangent vector field, one needs to evaluate the expansion and the Fourier coefficients of vector spherical harmonics. In this article, we develop fast algorithms (FaVeST) for vector spherical harmonic transforms on these evaluations. The forward FaVeST evaluates the Fourier coefficients and has a computational cost proportional to N log √ N for N number of evaluation points. The adjoint FaVeST, which evaluates a linear combination of vector spherical harmonics with a degree up to ⊡ M for M evaluation points, has cost proportional to M log √ M . Numerical examples of simulated tangent fields illustrate the accuracy, efficiency, and stability of FaVeST.

F 2 -linear pseudorandom number generators are very popular due to their high speed, to the ease with which generators with a sizable state space can be created, and to their provable theoretical properties. However, they suffer from linear artifacts that show as failures in linearity-related statistical tests such as the binary-rank and the linear-complexity test. In this article, we give two new contributions. First, we introduce two new F 2 -linear transformations that have been handcrafted to have good statistical properties and at the same time to be programmable very efficiently on superscalar processors, or even directly in hardware. Then, we describe some scramblers , that is, nonlinear functions applied to the state array that reduce or delete the linear artifacts, and propose combinations of linear transformations and scramblers that give extremely fast pseudorandom number generators of high quality. A novelty in our approach is that we use ideas from the theory of filtered linear-feedback shift registers to prove some properties of our scramblers, rather than relying purely on heuristics. In the end, we provide simple, extremely fast generators that use a few hundred bits of memory, have provable properties, and pass strong statistical tests.

We consider the problem of low-rank approximation of massive dense nonnegative tensor data, for example, to discover latent patterns in video and imaging applications. As the size of data sets grows, single workstations are hitting bottlenecks in both computation time and available memory. We propose a distributed-memory parallel computing solution to handle massive data sets, loading the input data across the memories of multiple nodes, and performing efficient and scalable parallel algorithms to compute the low-rank approximation. We present a software package called Parallel Low-rank Approximation with Nonnegativity Constraints, which implements our solution and allows for extension in terms of data (dense or sparse, matrices or tensors of any order), algorithm (e.g., from multiplicative updating techniques to alternating direction method of multipliers), and architecture (we exploit GPUs to accelerate the computation in this work). We describe our parallel distributions and algorithms, which are careful to avoid unnecessary communication and computation, show how to extend the software to include new algorithms and/or constraints, and report efficiency and scalability results for both synthetic and real-world data sets.

Fuzzy regression provides an alternative to statistical regression when the model is indefinite, the relationships between model parameters are vague, the sample size is low, or the data are hierarchically structured. Such cases allow to consider the choice of a regression model based on the fuzzy set theory. In fuzzyreg, we implement fuzzy linear regression methods that differ in the expectations of observational data types, outlier handling, and parameter estimation method. We provide a wrapper function that prepares data for fitting fuzzy linear models with the respective methods from a syntax established in R for fitting regression models. The function fuzzylm thus provides a novel functionality for R through standardized operations with fuzzy numbers. Additional functions allow for conversion of real-value variables to be fuzzy numbers, printing, summarizing, model plotting, and calculation of model predictions from new data using supporting functions that perform arithmetic operations with triangular fuzzy numbers. Goodness of fit and total error of the fit measures allow model comparisons. The package contains a dataset named bats with measurements of temperatures of hibernating bats and the mean annual surface temperature reflecting the climate at the sampling sites. The predictions from fuzzy linear models fitted to this dataset correspond well to the observed biological phenomenon. Fuzzy linear regression has great potential in predictive modeling where the data structure prevents statistical analysis and the modeled process exhibits inherent fuzziness.

In this article, we propose the Fast Algorithms for Maxwell’s Equations (FAME) package for solving Maxwell’s equations for modeling three-dimensional photonic crystals. FAME combines the null-space free method with fast Fourier transform (FFT)-based matrix-vector multiplications to solve the generalized eigenvalue problems (GEPs) arising from Yee’s discretization. The GEPs are transformed into a null-space free standard eigenvalue problem with a Hermitian positive-definite coefficient matrix. The computation times for FFT-based matrix-vector multiplications with matrices of dimension 7 million are only 0.33 and 3.6 × 10 − 3 seconds using MATLAB with an Intel Xeon CPU and CUDA C++ programming with a single NVIDIA Tesla P100 GPU, respectively. Such multiplications significantly reduce the computational costs of the conjugate gradient method for solving linear systems. We successfully use FAME on a single P100 GPU to solve a set of GEPs with matrices of dimension more than 19 million, in 127 to 191 seconds per problem. These results demonstrate the potential of our proposed package to enable large-scale numerical simulations for novel physical discoveries and engineering applications of photonic crystals.

We present a new algorithm for solving the dense linear (sum) assignment problem and an efficient, parallel implementation that is based on the successive shortest path algorithm. More specifically, we introduce the well-known epsilon scaling approach used in the Auction algorithm to approximate the dual variables of the successive shortest path algorithm prior to solving the assignment problem to limit the complexity of the path search. This improves the runtime by several orders of magnitude for hard-to-solve real-world problems, making the runtime virtually independent of how hard the assignment is to find. In addition, our approach allows for using accelerators and/or external compute resources to calculate individual rows of the cost matrix. This enables us to solve problems that are larger than what has been reported in the past, including the ability to efficiently solve problems whose cost matrix exceeds the available systems memory. To our knowledge, this is the first implementation that is able to solve problems with more than one trillion arcs in less than 100 hours on a single machine.

In this article, we introduce the Python framework PyMGRIT, which implements the multigrid-reduction-in-time (MGRIT) algorithm for solving (non-)linear systems arising from the discretization of time-dependent problems. The MGRIT algorithm is a reduction-based iterative method that allows parallel-in-time simulations, i.e., calculating multiple time steps simultaneously in a simulation, using a time-grid hierarchy. The PyMGRIT framework includes many different variants of the MGRIT algorithm, ranging from different multigrid cycle types and relaxation schemes, various coarsening strategies, including time-only and space-time coarsening, and the ability to utilize different time integrators on different levels in the multigrid hierachy. The comprehensive documentation with tutorials and many examples and the fully documented code allow an easy start into the work with the package. The functionality of the code is ensured by automated serial and parallel tests using continuous integration. PyMGRIT supports serial runs suitable for prototyping and testing of new approaches, as well as parallel runs using the Message Passing Interface (MPI). In this manuscript, we describe the implementation of the MGRIT algorithm in PyMGRIT and present the usage from both a user and a developer point of view. Three examples illustrate different aspects of the package itself, especially running tests with pure time parallelism, as well as space-time parallelism through the coupling of PyMGRIT with PETSc or Firedrake.

Dense kernel matrices represented in H 2 matrix format typically require less storage and have faster matrix-vector multiplications than when these matrices are represented in the standard dense format. In this article, we present H2Pack, a high-performance, shared-memory library for constructing and operating with H 2 matrix representations for kernel matrices defined by non-oscillatory, translationally invariant kernel functions. Using a hybrid analytic-algebraic compression method called the proxy point method, H2Pack can efficiently construct an H 2 matrix representation with linear computational complexity. Storage and matrix-vector multiplication also have linear complexity. H2Pack also introduces the concept of “partially admissible blocks” for H 2 matrices to make H 2 matrix-vector multiplication mathematically identical to the fast multipole method (FMM) if analytic expansions are used. We optimize H2Pack from both the algorithm and software perspectives. Compared to existing FMM libraries, H2Pack generally has much faster H 2 matrix-vector multiplications, since the proxy point method is more effective at producing block low-rank approximations than the analytic methods used in FMM. As a tradeoff, H 2 matrix construction in H2Pack is typically more expensive than the setup cost in FMM libraries. Thus, H2Pack is ideal for applications that need a large number of matrix-vector multiplications for a given configuration of data points.

We develop fast and accurate algorithms for evaluating √x 2 +y 2 for two floating-point numbers x and y . Library functions that perform this computation are generally named hypot(x,y). We compare five approaches that we will develop in this article to the current resident library function that is delivered with Julia 1.1 and to the code that has been distributed with the C math library for decades. We will investigate the accuracy of our algorithms by simulation.

We recently introduced a sparse stretching strategy for handling dense rows that can arise in large-scale linear least-squares problems and make such problems challenging to solve. Sparse stretching is designed to limit the amount of fill within the stretched normal matrix and hence within the subsequent Cholesky factorization. While preliminary results demonstrated that sparse stretching performs significantly better than standard stretching, it has a number of limitations. In this article, we discuss and illustrate these limitations and propose new strategies that are designed to overcome them. Numerical experiments on problems arising from practical applications are used to demonstrate the effectiveness of these new ideas. We consider both direct and preconditioned iterative solvers.

A continuous spectral algorithm and computer routines in the R programming environment that enable the simulation of second-order stationary and intrinsic (i.e., with second-order stationary increments or generalized increments) vector Gaussian random fields in Euclidean spaces are presented. The simulation is obtained by computing a weighted sum of cosine and sine waves, with weights that depend on the matrix-valued spectral density associated with the spatial correlation structure of the random field to simulate. The computational cost is proportional to the number of locations targeted for simulation, below that of sequential, matrix decomposition and discrete spectral algorithms. Also, the implementation is versatile, as there is no restriction on the number of vector components, workspace dimension, number and geometrical configuration of the target locations. The computer routines are illustrated with synthetic examples and statistical testing is proposed to check the normality of the distribution of the simulated random field or of its generalized increments. A by-product of this work is a spectral representation of spherical, cubic, penta, Askey, J-Bessel, Cauchy, Laguerre, hypergeometric, iterated exponential, gamma, and stable covariance models in the d -dimensional Euclidean space.

DELAUNAYSPARSE contains both serial and parallel codes written in Fortran 2003 (with OpenMP) for performing medium- to high-dimensional interpolation via the Delaunay triangulation. To accommodate the exponential growth in the size of the Delaunay triangulation in high dimensions, DELAUNAYSPARSE computes only a sparse subset of the complete Delaunay triangulation, as necessary for performing interpolation at the user specified points. This article includes algorithm and implementation details, complexity and sensitivity analyses, usage information, and a brief performance study.

Popular deep neural networks (DNNs) spend the majority of their execution time computing convolutions. The Winograd family of algorithms can greatly reduce the number of arithmetic operations required and is used in many DNN software frameworks. However, the performance gain is at the expense of a reduction in floating point (FP) numerical accuracy. In this article, we analyse the worst-case FP error and derive an estimation of the norm and conditioning of the algorithm. We show that the bound grows exponentially with the size of the convolution. Further, the error bound of the modified algorithm is slightly lower but still exponential. We propose several methods for reducing FP error. We propose a canonical evaluation ordering based on Huffman coding that reduces summation error. We study the selection of sampling “points” experimentally and find empirically good points for the most important sizes. We identify the main factors associated with good points. In addition, we explore other methods to reduce FP error, including mixed-precision convolution, and pairwise summation across DNN channels. Using our methods, we can significantly reduce FP error for a given block size, which allows larger block sizes to be used and reduced computation.

Adaptive mesh refinement (AMR) is often used when solving time-dependent partial differential equations using numerical methods. It enables time-varying regions of much higher resolution, which can be used to track discontinuities in the solution by selectively refining around those areas. The open source Clawpack software implements block-structured AMR to refine around propagating waves in the AMRClaw package. For problems where the solution must be computed over a large domain but is only of interest in a small area this approach often refines waves that will not impact the target area. We seek a method that enables the identification and refinement of only the waves that will influence the target area. Here we show that solving the time-dependent adjoint equation and using a suitable inner product allows for a more precise refinement of the relevant waves. We present the adjoint methodology in general, and give details on how this method has been implemented in AMRClaw. Examples for linear acoustics equations are presented, and a computational performance analysis is conducted. The adjoint method is compared to AMR methods already available in the AMRClaw software, and the resulting advantages and disadvantages are discussed. The code for the examples presented is archived on Github., Submitted to ACM Transactions on Mathematical Software (TOMS)

We define “reproducibility” as getting bitwise identical results from multiple runs of the same program, perhaps with different hardware resources or other changes that should not affect the answer. Many users depend on reproducibility for debugging or correctness. However, dynamic scheduling of parallel computing resources, combined with nonassociative floating point addition, makes reproducibility challenging even for summation, or operations like the BLAS. We describe a “reproducible accumulator” data structure (the “binned number”) and associated algorithms to reproducibly sum binary floating point numbers, independent of summation order. We use a subset of the IEEE Floating Point Standard 754-2008 and bitwise operations on the standard representations in memory. Our approach requires only one read-only pass over the data, and one reduction in parallel, using a 6-word reproducible accumulator (more words can be used for higher accuracy), enabling standard tiling optimization techniques. Summing n words with a 6-word reproducible accumulator requires approximately 9 n floating point operations (arithmetic, comparison, and absolute value) and approximately 3 n bitwise operations. The final error bound with a 6-word reproducible accumulator and our default settings can be up to 2 29 times smaller than the error bound for conventional (recursive) summation on ill-conditioned double-precision inputs.

The Singular Value Decomposition (SVD) is widely used in numerical analysis and scientific computing applications, including dimensionality reduction, data compression and clustering, and computation of pseudo-inverses. In many cases, a crucial part of the SVD of a general matrix is to find the SVD of an associated bidiagonal matrix. This article discusses an algorithm to compute the SVD of a bidiagonal matrix through the eigenpairs of an associated symmetric tridiagonal matrix. The algorithm enables the computation of only a subset of singular values and corresponding vectors, with potential performance gains. The article focuses on a sequential version of the algorithm, and discusses special cases and implementation details. The implementation, called BDSVDX, has been included in the LAPACK library. We use a large set of bidiagonal matrices to assess the accuracy of the implementation, both in single and double precision, as well as to identify potential shortcomings. The results show that BDSVDX can be up to three orders of magnitude faster than existing algorithms, which are limited to the computation of a full SVD. We also show comparisons of an implementation that uses BDSVDX as a building block for the computation of the SVD of general matrices.

Aiming to provide a very accurate, efficient, and robust quartic equation solver for physical applications, we have proposed an algorithm that builds on the previous works of P. Strobach and S. L. Shmakov. It is based on the decomposition of the quartic polynomial into two quadratics, whose coefficients are first accurately estimated by handling carefully numerical errors and afterward refined through the use of the Newton-Raphson method. Our algorithm is very accurate in comparison with other state-of-the-art solvers that can be found in the literature, but (most importantly) it turns out to be very efficient according to our timing tests. A crucial issue for us is the robustness of the algorithm, i.e., its ability to cope with the detrimental effect of round-off errors, no matter what set of quartic coefficients is provided in a practical application. In this respect, we extensively tested our algorithm in comparison to other quartic equation solvers both by considering specific extreme cases and by carrying out a statistical analysis over a very large set of quartics. Our algorithm has also been heavily tested in a physical application, i.e., simulations of hard cylinders, where it proved its absolute reliability as well as its efficiency.

QNSTOP consists of serial and parallel (OpenMP) Fortran 2003 codes for the quasi-Newton stochastic optimization method of Castle and Trosset for stochastic search problems. A complete description of QNSTOP for both local search with stochastic objective and global search with “noisy” deterministic objective is given here, to the best of our knowledge, for the first time. For stochastic search problems, some convergence theory exists for particular algorithmic choices and parameter values. Both the parallel driver subroutine, which offers several parallel decomposition strategies, and the serial driver subroutine can be used for local stochastic search or global deterministic search, based on an input switch. Some performance data for computational systems biology problems is given.

Conventional Graphics Processing Unit (GPU) implementations of Strassen’s algorithm (S trassen ) rely on the existing high-performance matrix multiplication ( gemm ), trading space for time. As a result, such approaches can only achieve practical speedup for relatively large, “squarish” matrices due to the extra memory overhead, and their usages are limited due to the considerable workspace. We present novel S trassen primitives for GPUs that can be composed to generate a family of S trassen algorithms. Our algorithms utilize both the memory and thread hierarchies on GPUs, reusing shared memory and register files inherited from gemm , fusing additional operations, and avoiding extra workspace. We further exploit intra- and inter-kernel parallelism by batching, streaming, and employing atomic operations. We develop a performance model for NVIDIA Volta GPUs to select the appropriate blocking parameters and predict the performance for gemm and S trassen . Overall, our 1-level S trassen can achieve up to 1.11× speedup with a crossover point as small as 1,536 compared to cublasSgemm on a NVIDIA Tesla V100 GPU. With additional workspace, our 2-level S trassen can achieve 1.19× speedup with a crossover point at 7,680.

Partitioning graphs is a common and useful operation in many areas, from parallel computing to VLSI design to sparse matrix algorithms. In this article, we introduce Mongoose, a multilevel hybrid graph partitioning algorithm and library. Building on previous work in multilevel partitioning frameworks and combinatoric approaches, we introduce novel stall-reducing and stall-free coarsening strategies, as well as an efficient hybrid algorithm leveraging (1) traditional combinatoric methods and (2) continuous quadratic programming formulations. We demonstrate how this new hybrid algorithm outperforms either strategy in isolation, and we also compare Mongoose to METIS and demonstrate its effectiveness on large and social networking (power law) graphs.

Iterated-integral signatures and log signatures are sequences calculated from a path that characterizes its shape. They originate from the work of K. T. Chen and have become important through Terry Lyons’s theory of differential equations driven by rough paths, which is an important developing area of stochastic analysis. They have applications in statistics and machine learning, where there can be a need to calculate finite parts of them quickly for many paths. We introduce the signature and the most basic information (displacement and signed areas) that it contains. We present algorithms for efficiently calculating these signatures. For log signatures this requires consideration of the structure of free Lie algebras. We benchmark the performance of the algorithms. The methods are implemented in C++ and released as a Python extension package, which also supports differentiation. In combination with a machine learning library (Tensorflow, PyTorch, or Theano), this allows end-to-end learning of neural networks involving signatures.

We present a computational procedure to evaluate the integral ∫ y x s p -1 e -μs ds for 0 ≤ x < y ≤ +∞,μ = ±1, p > 0, which generalizes the lower ( x =0) and upper ( y =+∞) incomplete gamma functions. To allow for large values of x , y , and p while avoiding under/overflow issues in the standard double precision floating point arithmetic, we use an explicit normalization that is much more efficient than the classical ratio with the complete gamma function. The generalized incomplete gamma function is estimated with continued fractions, with integrations by parts, or, when x ≈ y , with the Romberg numerical integration algorithm. We show that the accuracy reached by our algorithm improves a recent state-of-the-art method by two orders of magnitude, and it is essentially optimal considering the limitations imposed by floating point arithmetic. Moreover, the admissible parameter range of our algorithm (0 ≤ p,x,y ≤ 10 15 ) is much larger than competing algorithms, and its robustness is assessed through massive usage in an image processing application.

In a recent article, Lucambio Pérez and Prudente extended the Wolfe conditions for the vector-valued optimization. Here, we propose a line search algorithm for finding a step size satisfying the strong Wolfe conditions in the vector optimization setting. Well definedness and finite termination results are provided. We discuss practical aspects related to the algorithm and present some numerical experiments illustrating its applicability. Codes supporting this article are written in Fortran 90 and are freely available for download.

The maximum flow problem is one of the most common network flow problems. This problem involves finding the maximum possible amount of flow between two designated nodes on a network with arcs having flow capacities. The push-relabel algorithm is one of the fastest algorithms to solve this problem. We present a shared memory parallel push-relabel algorithm. Graph coloring is used to avoid collisions between threads for concurrent push and relabel operations. In addition, excess values of target nodes are updated using atomic instructions to prevent race conditions. The experiments show that our algorithm is competitive for wide graphs with low diameters. Results from three different data sets are included, computer vision problems, DIMACS challenge problems, and KaHIP partitioning problems. These are compared with existing push-relabel and pseudoflow implementations. We show that high speedup rates are possible using our coloring based parallelization technique on sparse networks. However, we also observe that the pseudoflow algorithm runs faster than the push-relabel algorithm on dense and long networks.

Multi-degree splines are smooth piecewise-polynomial functions where the pieces can have different degrees. We describe a simple algorithmic construction of a set of basis functions for the space of multi-degree splines with similar properties to standard B-splines. These basis functions are called multi-degree B-splines (or MDB-splines ). The construction relies on an extraction operator that represents all MDB-splines as linear combinations of local B-splines of different degrees. This enables the use of existing efficient algorithms for B-spline evaluations and refinements in the context of multi-degree splines. A M ATLAB implementation is provided to illustrate the computation and use of MDB-splines.

IPscatt is a free, open-source MATLAB toolbox facilitating the solution for time-independent scattering (also known as time-harmonic scattering) in two- and three-dimensional settings. The toolbox has three main application cases: simulation of the scattered field for a given transmitter-receiver geometry; the generation of simulated data as well as the handling of the real-world data from Institute Fresnel; and the reconstruction of the contrast from several measured, scattered fields. In each case, a variety of options tailored to the needs of practitioners is provided. For example, the toolbox allows the simulation of the scattered near field as well as of the far field. Also, it provides methods for the modeling of the incident field as point sources as well as plane waves. Finally, many common geometries of transmitters and receivers are included out of the box. Regarding the reconstruction, the provided functions implement the regularization scheme that relies on a primal-dual algorithm and was introduced by F. Bürgel, K. S. Kazimierski, and A. Lechleiter [ Journal of Computational Physics 339 (2017), 1–30]. This article provides a survey of the mathematical concepts in scattering, connects them with the provided implementation, gives an overview of the software framework as well as its application areas, and compares it with existing software packages solving the same problem.

A bottom-up approach to parallel anisotropic mesh generation is presented by building a mesh generator starting from the basic operations of vertex insertion and Delaunay triangles. Applications focusing on high-lift design or dynamic stall, or numerical methods and modeling test cases, still focus on two-dimensional domains. This automated parallel mesh generation approach can generate high-fidelity unstructured meshes with anisotropic boundary layers for use in the computational fluid dynamics field. The anisotropy requirement adds a level of complexity to a parallel meshing algorithm by making computation depend on the local alignment of elements, which in turn is dictated by geometric boundaries and the density functions— one-dimensional spacing functions generated from an exponential distribution. This approach yields computational savings in mesh generation and flow solution through well-shaped anisotropic triangles instead of isotropic triangles. The validity of the meshes is shown through solution characteristic comparisons to verified reference solutions. A 79% parallel weak scaling efficiency on 1,024 distributed memory nodes, and a 72% parallel efficiency over the fastest sequential isotropic mesh generator on 512 distributed memory nodes, is shown through numerical experiments.

The minimum distance of an error-correcting code is an important concept in information theory. Hence, computing the minimum distance of a code with a minimum computational cost is crucial to many problems in this area. In this article, we present and assess a family of implementations of both the brute-force algorithm and the Brouwer-Zimmermann algorithm for computing the minimum distance of a random linear code over F 2 that are faster than current implementations, both in the commercial and public domain. In addition to the basic sequential implementations, we present parallel and vectorized implementations that produce high performances on modern architectures. The attained performance results show the benefits of the developed optimized algorithms, which obtain remarkable improvements compared with state-of-the-art implementations widely used nowadays.

An algorithm for multiplying a chain of Kronecker products by a matrix is described. The algorithm does not require that the Kronecker chain actually be computed and the main computational work is a series of matrix-matrix multiplications. Use of the algorithm can lead to substantial savings in both memory requirements and computational speed. Although similar algorithms have been described before, this article makes two novel contributions. First, it shows how shuffling of data can be (largely) avoided. Second, it provides a simple method to determine the optimal ordering of the workflow.

In addition to tensor contractions, one of the most pronounced computational bottlenecks in the nonorthogonally spin-adapted forms of the quantum chemistry methods CCSDT and CCSDTQ, and their approximate forms—including CCSD(T) and CCSDT(Q)—are spin summations. At a first sight, spin summations are operations similar to tensor transpositions, but a closer look reveals additional challenges to high-performance calculations, including temporal locality and scattered memory accesses. This article explores a sequence of algorithmic solutions for spin summations, each exploiting individual properties of either the underlying hardware (e.g., caches, vectorization) or the problem itself (e.g., factorizability). The final algorithm combines the advantages of all the solutions while avoiding their drawbacks; this algorithm achieves high performance through parallelization and vectorization, and by exploiting the temporal locality inherent to spin summations. Combined, these optimizations result in speedups between 2.4× and 5.5× over the NCC quantum chemistry software package. In addition to such a performance boost, our algorithm can perform the spin summations in-place , thus reducing the memory footprint by 2× over an out-of-place variant.

We provide a number of corrections to the software component that accompanied this Algorithm submission [3]. An updated version of the code is available from the ACM Collected Algorithms site [1].