David Bindel
James Demmel
William Kahan
Osni Marques
Abstract
We consider the efficient and accurate computation of Givens rotations. When f and g are positive real numbers, this simply amounts to computing the values of c = f/√f2 + g2, s = g/√f2 + g2, and r = √f2 + g2. This apparently trivial computation merits closer consideration for the following three reasons. First, while the definitions of c, s and r seem obvious in the case of two nonnegative arguments f and g, there is enough freedom of choice when one or more of f and g are negative, zero or complex that LAPACK auxiliary routines SLARTG, CLARTG, SLARGV and CLARGV can compute rather different values of c, s and r for mathematically identical values of f and g. To eliminate this unnecessary ambiguity, the BLAS Technical Forum chose a single consistent definition of Givens rotations that we will justify here. Second, computing accurate values of c, s and r as efficiently as possible and reliably despite over/underflow is surprisingly complicated. For complex Givens rotations, the most efficient formulas require only one real square root and one real divide (as well as several much cheaper additions and multiplications), but a reliable implementation using only working precision has a number of cases. On a Sun Ultra-10, the new implementation is slightly faster than the previous LAPACK implementation in the most common case, and 2.7 to 4.6 times faster than the corresponding vendor, reference or ATLAS routines. It is also more reliable; all previous codes occasionally suffer from large inaccuracies due to over/underflow. For real Givens rotations, there are also improvements in speed and accuracy, though not as striking. Third, the design process that led to this reliable implementation is quite systematic, and could be applied to the design of similarly reliable subroutines.
- Anderson, E. 2000. Discontinuous plane rotations and the symmetric eigenvalue problem. Computer Science Dept. Tech. Rep. CS-00-454. Univ. Tennessee, Knoxville, Tenn. www.netlib.org/lapack/lawns/lawn150.ps.Google Scholar
- Anderson, E., Bai, Z., Bischof, C., Demmel, J., Dongarra, J., Du Croz, J., Greenbaum, A., Hammarling, S., McKenney, A., Blackford, S., and Sorensen, D. 1999. LAPACK Users' Guide 3rd ed. SIAM, Philadelphia, Pa. Google Scholar
- Anderson, E. and Fahey, M. 1997. Performance improvements to LAPACK for the Cray scientific library. Computer Science Dept. Tech. Rep. CS-97-359. Univ. Tennessee, Knoxville, Tenn. www.netlib.org/lapack/lawns/lawn126.ps. Google Scholar
- ANSI/IEEE Std 754-1985, IEEE Standard for Binary Floating Point Arithmetic. ANSI/IEEE, New York.Google Scholar
- Bindel, D., Demmel, J., Kahan, W., and Marques, O. 2002. Software for reliable and efficient Givens rotations. www.cs.berkeley.edu/∼demmel/Givens.Google Scholar
- Blackford, S., Corliss, G., Demmel, J., Dongarra, J., Duff, I., Hammarling, S., Henry, G., Heroux, M., Hu, C., Kahan, W., Kaufman, L., Kearfott, B., Krogh, F., Li, X., Maany, Z., Petitet, A., Pozo, R., Remington, K., Walster, W., Whaley, C., Wolff v. Gudenberg, J., and Lumsdaine, A. 2001. Basic Linear Algebra Subprograms Technical (BLAST) Forum Standard. Intern. J. High Performance Comput. 15, 3--4 (also available at www.netlib.org/blas/blast-forum/).Google Scholar
- Demmel, J. 1997. Applied Numerical Linear Algebra. SIAM, Philadelphia, Pa. Google Scholar
- Demmel, J. and Li, X. 1994. Faster numerical algorithms via exception handling. IEEE Trans. Comput. 43, 8, 983--992. LAPACK Working Note 59. Google Scholar
- Dongarra, J., Bunch, J., Moler, C., and Stewart, G. W. 1979. LINPACK User's Guide. SIAM, Philadelphia, Pa.Google Scholar
- Golub, G. and Van Loan, C. 1996. Matrix Computations, 3rd ed. Johns Hopkins University Press, Baltimore, Md. Google Scholar
- Hull, T., Fairgrieve, T., and Tang, P. T. P. 1994. Implementing complex elementary functions using exception handling. ACM Trans. Math. Softw. 20, 2 (June), 215--244. Google Scholar
- Hull, T., Fairgrieve, T., and Tang, P. T. P. 1997. Implementing complex Arcsine and Arccosine using exception handling. ACM Trans. Math. Softw. 23, 3 (Sept), 299--335. Google Scholar
- IEEE 754-2002. IEEE Standard for Binary Floating Point Arithmetic Revision.grouper.ieee.org/groups/754.Google Scholar
- Kato, T. 1980. Perturbation Theory for Linear Operators, 2 ed. Springer-Verlag, Berlin, Germany.Google Scholar
- Kunkel, P. and Mehrmann, V. 1991. Smooth factorizations of matrix valued functions and their derivatives. Num. Math. 60, 115--132.Google Scholar
- Lawson, C., Hanson, R., Kincaid, D., and Krogh, F. 1979. Basic linear algebra subprograms for Fortran usage. ACM Trans. Math. Softw. 5, 308--323. Google Scholar
- Smith, B. T., Boyle, J. M., Dongarra, J. J., Garbow, B. S., Ikebe, Y., Klema, V. C., and Moler, C. B. 1976. Matrix Eigensystem Routines---EISPACK Guide. Lecture Notes in Computer Science, vol. 6. Springer-Verlag, Berlin, Germany.Google Scholar
- Stewart, G. W. 1976. The economical storage of plane rotations. Numer. Math. 25, 2, 137--139.Google Scholar
- Wilkinson, J. H. 1965. The Algebraic Eigenvalue Problem. Oxford University Press, Oxford, England. Google Scholar
Index Terms
- On computing givens rotations reliably and efficiently
Recommendations
Algorithm 830: Another visit with standard and modified Givens transformations and a remark on Algorithm 539
First we report on a correction and improvement to the Level 1 BLAS routine srotmg for computing the Modified Givens Transformation (MG). We then, in the light of the performance of the code on modern compiler/hardware combinations, reconsider the ...
On General Row Merging Schemes for Sparse Givens Transformations
This paper introduces general row merging schemes for the $QR$ decomposition of sparse matrices by Givens rotations. They can be viewed as a generalization of row rotations to submatrix rotations (or merging) in the recent method by George and Heath [12]...
Modified incomplete orthogonal factorization methods using Givens rotations
We present a class of new preconditioners based on the incomplete Givens orthogonalization (IGO) methods for solving large sparse systems of linear equations. In the new methods, instead of dropping entries and accepting fill-ins according to the ...
Comments