On computing givens rotations reliably and efficiently

Abstract

We consider the efficient and accurate computation of Givens rotations. Whenfandgare positive real numbers, this simply amounts to computing the values ofc=f/√f2+g2,s=g/√f2+g2, andr= √f2+g2. This apparently trivial computation merits closer consideration for the following three reasons. First, while the definitions ofc,sandrseem obvious in the case of two nonnegative argumentsfandg, there is enough freedom of choice when one or more offandgare negative, zero or complex that LAPACK auxiliary routines SLARTG, CLARTG, SLARGV and CLARGV can compute rather different values ofc,sandrfor mathematically identical values offandg. 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 ofc,sandras 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.