No Arabic abstract
Gaussian elimination with full pivoting generates a PLUQ matrix decomposition. Depending on the strategy used in the search for pivots, the permutation matrices can reveal some information about the row or the column rank profiles of the matrix. We propose a new pivoting strategy that makes it possible to recover at the same time both row and column rank profiles of the input matrix and of any of its leading sub-matrices. We propose a rank-sensitive and quad-recursive algorithm that computes the latter PLUQ triangular decomposition of an m times n matrix of rank r in O(mnr^{omega-2}) field operations, with omega the exponent of matrix multiplication. Compared to the LEU decomposition by Malashonock, sharing a similar recursive structure, its time complexity is rank sensitive and has a lower leading constant. Over a word size finite field, this algorithm also improveLs the practical efficiency of previously known implementations.
We describe a strategy for rigorous arbitrary-precision evaluation of Legendre polynomials on the unit interval and its application in the generation of Gauss-Legendre quadrature rules. Our focus is on making the evaluation practical for a wide range of realistic parameters, corresponding to the requirements of numerical integration to an accuracy of about 100 to 100 000 bits. Our algorithm combines the summation by rectangular splitting of several types of expansions in terms of hypergeometric series with a fixed-point implementation of Bonnets three-term recurrence relation. We then compute rigorous enclosures of the Gauss-Legendre nodes and weights using the interval Newton method. We provide rigorous error bounds for all steps of the algorithm. The approach is validated by an implementation in the Arb library, which achieves order-of-magnitude speedups over previous code for computing Gauss-Legendre rules with simultaneous high degree and precision.
The aggregated unfitted finite element method (AgFEM) is a methodology recently introduced in order to address conditioning and stability problems associated with embedded, unfitted, or extended finite element methods. The method is based on removal of basis functions associated with badly cut cells by introducing carefully designed constraints, which results in well-posed systems of linear algebraic equations, while preserving the optimal approximation order of the underlying finite element spaces. The specific goal of this work is to present the implementation and performance of the method on distributed-memory platforms aiming at the efficient solution of large-scale problems. In particular, we show that, by considering AgFEM, the resulting systems of linear algebraic equations can be effectively solved using standard algebraic multigrid preconditioners. This is in contrast with previous works that consider highly customized preconditioners in order to allow one the usage of iterative solvers in combination with unfitted techniques. Another novelty with respect to the methods available in the literature is the problem sizes that can be handled with the proposed approach. While most of previous references discussing linear solvers for unfitted methods are based on serial non-scalable algorithms, we propose a parallel distributed-memory method able to efficiently solve problems at large scales. This is demonstrated by means of a weak scaling test defined on complex 3D domains up to 300M degrees of freedom and one billion cells on 16K CPU cores in the Marenostrum-IV platform. The parallel implementation of the AgFEM method is available in the large-scale finite element package FEMPAR.
The real symmetric tridiagonal eigenproblem is of outstanding importance in numerical computations; it arises frequently as part of eigensolvers for standard and generalized dense Hermitian eigenproblems that are based on a reduction to tridiagonal form. For its solution, the algorithm of Multiple Relatively Robust Representations (MRRR) is among the fastest methods. Although fast, the solvers based on MRRR do not deliver the same accuracy as competing methods like Divide & Conquer or the QR algorithm. In this paper, we demonstrate that the use of mixed precisions leads to improved accuracy of MRRR-based eigensolvers with limited or no performance penalty. As a result, we obtain eigensolvers that are not only equally or more accurate than the best available methods, but also -in most circumstances- faster and more scalable than the competition.
ODE Test Problems (OTP) is an object-oriented MATLAB package offering a broad range of initial value problems which can be used to test numerical methods such as time integration methods and data assimilation (DA) methods. It includes problems that are linear and nonlinear, homogeneous and nonhomogeneous, autonomous and nonautonomous, scalar and high-dimensional, stiff and nonstiff, and chaotic and nonchaotic. Many are real-world problems from fields such as chemistry, astrophysics, meteorology, and electrical engineering. OTP also supports partitioned ODEs for testing IMEX methods, multirate methods, and other multimethods. Functions for plotting solutions and creating movies are available for all problems, and exact solutions are provided when available. OTP is desgined for ease of use-meaning that working with and modifying problems is simple and intuitive.
The detection of insufficiently resolved or ill-conditioned integrand structures is critical for the reliability assessment of the quadrature rule outputs. We discuss a method of analysis of the profile of the integrand at the quadrature knots which allows inferences approaching the theoretical 100% rate of success, under error estimate sharpening. The proposed procedure is of the highest interest for the solution of parametric integrals arising in complex physical models.