Search

Pipelined Krylov subspace methods typically offer improved strong scaling on par-allel HPC hardware compared to standard Krylov subspace methods for large and sparse linear systems. In pipelined methods the traditional synchronization bottleneck is mitigated by overlap- ping time-consuming global communications with useful computations. However, to achieve this communication-hiding strategy, pipelined methods introduce additional recurrence relations for a number of auxiliary variables that are required to update the approximate solution. This paper aims at studying the in uence of local rounding errors that are introduced by the additional recurrences in the pipelined Conjugate Gradient (CG) method. Specifically, we analyze the impact of local round-off effects on the attainable accuracy of the pipelined CG algorithm and compare it to the tra- ditional CG method. Furthermore, we estimate the gap between the true residual and the recursively computed residual used in the algorithm. Based on this estimate we suggest an automated residual replacement strategy to reduce the loss of attainable accuracy on the final iterative solution. The resulting pipelined CG method with residual replacement improves the maximal attainable accuracy of pipelined CG while maintaining the efficient parallel performance of the pipelined method. This conclusion is substantiated by numerical results for a variety of benchmark problems. [ABSTRACT FROM AUTHOR]

The roundoff-error-free (REF) LU and Cholesky factorizations, combined with the REF substitution algorithms, allow rational systems of linear equations to be solved exactly and efficiently by working entirely in integer arithmetic. These REF computational tools share two key properties: their constituent divisions are exact, and their matrix entries' bit-lengths (i.e., required number of storage bits) are bounded polynomially. However, as this paper explains, updating the REF factorizations via the common delete-insert-reduce approach turns out to be inefficient in terms of operand bit-length growth and increased computational effort. In fact, we prove that this inefficiency can eliminate the computational savings expected of the factorization update process. Hence, the current work develops REF update algorithms that differ significantly from their traditional counterparts. To achieve this, it introduces the frame matrix methodology for representing and applying the REF factorizations, and then it devises new REF operations and shortcuts for transitioning between "nearby" frame matrices. The presented push-and-swap update approach utilizes these transitions to preserve the special structure of the REF factorizations and to avoid additional growth in their matrix entries' bit-lengths; this update process requires the expected O(n²) operations. The featured updates are column addition, deletion, and replacement with respect to the REF LU factorization. We also prove that analogous row updates can be performed via the column updates and discuss special considerations for updating the REF Cholesky factorization. An accompanying set of computational tests attests that for fully dense random instances, the push-and-swap column replacement update for the REF LU factorization is on average an order of magnitude faster than the corresponding delete-insert-reduce update for the exact rational arithmetic LU factorization. Altogether, the REF update algorithms represent a promising preliminary step toward enhancing the basic solution validation subroutines employed by state-of-the-art mixed-precision methodologies for exact linear and mixed-integer programming. [ABSTRACT FROM AUTHOR]

In this paper it is shown that the SVD of a matrix can be constructed efficiently in a hierarchical approach. The proposed algorithm is proven to recover the singular values and left singular vectors of the input matrix A if its rank is known. Further, the hierarchical algorithm can be used to recover the d largest singular values and left singular vectors with bounded error. It is also shown that the proposed method is stable with respect to round-off errors or corruption of the original matrix entries. Numerical experiments validate the proposed algorithms and parallel cost analysis. [ABSTRACT FROM AUTHOR]

We study the numerical behavior of stationary one-step or two-step matrix splitting iteration methods for solving large sparse systems of linear equations. We show that inexact solutions of inner linear systems associated with the matrix splittings may considerably influence the accuracy of the approximate solutions computed in finite precision arithmetic. For a general stationary matrix splitting iteration method, we analyze two mathematically equivalent implementations and discuss the conditions when they are componentwise or normwise forward or backward stable. We distinguish two different forms of matrix splitting iteration methods and prove that one of them is significantly more accurate than the other when employing inexact inner solves. The theoretical results are illustrated by numerical experiments with an inexact one-step and an inexact two-step splitting iteration method. [ABSTRACT FROM AUTHOR]

In practical computation with Runge-Kutta methods, the stage equations are not satisfied exactly, due to roundoff errors, algebraic solver errors, and so forth. We show by example that propagation of such errors within a single step can have catastrophic effects for otherwise practical and well-known methods. We perform a general analysis of internal error propagation, emphasizing that it depends significantly on how the method is implemented. We show that for a fixed method, essentially any set of internal stability polynomials can be obtained by modifying the implementation details. We provide bounds on the internal error amplification constants for some classes of methods with many stages, including strong stability preserving methods and extrapolation methods. These results are used to prove error bounds in the presence of roundoff or other internal errors. [ABSTRACT FROM AUTHOR]

In this paper we analyze the numerical behavior of several minimum residual methods which are mathematically equivalent to the GMRES method. Two main approaches are compared: one that computes the approximate solution in terms of a Krylov space basis from an upper triangular linear system for the coordinates, and one where the approximate solutions are updated with a simple recursion formula. We show that a different choice of the basis can significantly influence the numerical behavior of the resulting implementation. While Simpler GMRES and ORTHODIR are less stable due to the ill-conditioning of the basis used, the residual basis is well-conditioned as long as we have a reasonable residual norm decrease. These results lead to a new implementation, which is conditionally backward stable, and they explain the experimentally observed fact that the GCR method delivers very accurate approximate solutions when it converges fast enough without stagnation. [ABSTRACT FROM AUTHOR]

A strict adherence to threshold pivoting in the direct solution of symmetric indefinite problems can result in substantially more work and storage than forecast by a sparse analysis of the symmetric problem. One way of avoiding this is to use static pivoting where the data structures and pivoting sequence generated by the analysis are respected and pivots that would otherwise be very small are replaced by a user defined quantity. This can give a stable factorization but of a perturbed matrix. The conventional way of solving the sparse linear system is then to use iterative refinement, but there are cases where this fails to converge. In this paper, we discuss the use of more robust iterative methods, namely GMRES and flexible GMRES (FGMRES). We show both theoretically and experimentally that both approaches are more robust than iterative refinement and furthermore that FGMRES is far more robust than GMRES and that, under reasonable hypotheses, FGMRES is backward stable. We also show how restarted variants can be beneficial, although again the GMRES variant is not as robust as FGMRES. [ABSTRACT FROM AUTHOR]

The stable similarity reduction of a nonsymmetric square matrix to tridiagonal form has been a long-standing problem in numerical linear algebra. The biorthogonal Lanczos process is in principle a candidate method for this task, but in practice it is confined to sparse matrices and is restarted periodically because roundoff errors affect its three-term recurrence scheme and degrade the biorthogonality after a few steps. This adds to its vulnerability to serious breakdowns or near-breakdowns, the handling of which involves recovery strategies such as the look-ahead technique, which needs a careful implementation to produce a block-tridiagonal form with unpredictable block sizes. Other candidate methods, geared generally towards full matrices, rely on elementary similarity transformations that are prone to numerical instabilities. Such concomitant difficulties have hampered finding a satisfactory solution to the problem for either sparse or full matrices. This study focuses primarily on full matrices. After outlining earlier tridiagonalization algorithms from within a general framework, we present a new elimination technique combining orthogonal similarity transformations that are stable. We also discuss heuristics to circumvent breakdowns. Applications of this study include eigenvalue calculation and the approximation of matrix functions. [ABSTRACT FROM AUTHOR]

We present and analyze several simple algorithms for accurately computing the sum of n floating point numbers using a wider accumulator. Let f and F be the number of significant bits in the summands and the accumulator, respectively. Then assuming gradual underflow, no overflow, and round-to-nearest arithmetic, up to approximately 2F-f numbers can be added accurately by simply summing the terms in decreasing order of exponents, yielding a sum correct to within about 1.5 units in the last place (ulps). We apply this result to the floating point formats in the IEEE floating point standard. For example, a dot product of single precision vectors of length at most 33 computed using double precision and sorting is guaranteed correct to nearly 1.5 ulps. If double-extended precision is used, the vector length can be as large as 65,537. We also investigate how the cost of sorting can be reduced or eliminated while retaining accuracy. [ABSTRACT FROM AUTHOR]

We consider a 4n2-algorithm for computing the Gohberg­Semencul representation of the inverse of a nonsymmetric Toeplitz matrix and for obtaining the solution of a Toeplitz system. The original algorithm is unstable when one or more of the leading principal submatrices are ill-conditioned. In our stabilized version we slightly perturb the matrix when we encounter a division by a small number. As a consequence, the solution is also perturbed, and its accuracy can be improved by taking a small number of iterative refinement steps. We present a roundoff error analysis of the new algorithm as well as numerical results that support our analysis. [ABSTRACT FROM AUTHOR]

We present a roundoff error analysis of a null space method for solving quadratic programming minimization problems. This method combines the use of a direct QR factorization of the constraints with an iterative solver on the corresponding null space. Numerical experiments are presented which give evidence of the good performances of the algorithm on sparse matrices. [ABSTRACT FROM AUTHOR]

It is well known that the solution of the equality constrained least squares (LSE) problem minBx =d b - Ax 2 is the limit of the solution of the unconstrained weighted least squares problem min x µd b - µB A x 2 as the weight µ tends to infinity, assuming that [ BT AT ]T has full rank. We derive a method for the LSE problem by applying Householder QR factorization with column pivoting to this weighted problem and taking the limit analytically, with an appropriate rescaling of rows. The method ob- tained is a type of direct elimination method. We adapt existing error analysis for the unconstrained problem to obtain a row-wise backward error bound for the method. The bound shows that, provided row pivoting or row sorting is used, the method is well-suited to problems in which the rows of A and B vary widely in norm. As a by-product of our analysis, we derive a row-wise backward error bound of precisely the same form for the standard elimination method for solving the LSE problem. We illustrate our results with numerical tests. [ABSTRACT FROM AUTHOR]

A new method is investigated to reduce the roundoff error in computing derivatives using Chebyshev collocation methods. By using a grid mapping derived by Kosloff and Tal-Ezer, and the proper choice of the parameter α, the roundoff error of the kth derivative can be reduced from O(N2k) to O((N|ln ∈|)k), where ∈ is the machine precision and N is the number of collocation points. This drastic reduction of roundoff error makes mapped Chebyshev methods competitive with any other algorithm in computing second or higher derivatives with large N. Several other aspects of the mapped Chebyshev differentiation matrix are also studied, revealing that 1. the mapped Chebyshev methods require much less than π points to resolve a wave; 2. the eigenvalues are less sensitive to perturbation by roundoff error; and 3. larger time steps can be used for solving PDEs. All these advantages of the mapped Chebyshev methods can be achieved while maintaining spectral accuracy. [ABSTRACT FROM AUTHOR]