Editing Givens rotation (section)

== Stable calculation ==
When a Givens rotation matrix, {{math|''G''(''i'', ''j'', ''θ'')}}, multiplies another matrix, {{mvar|A}}, from the left, {{math|''G A''}}, only rows {{mvar|i}} and {{mvar|j}} of {{mvar|A}} are affected. Thus we restrict attention to the following counterclockwise problem. Given {{mvar|a}} and {{mvar|b}}, find {{math|1=''c'' = cos ''θ''}} and {{math|1=''s'' = sin ''θ''}} such that
:<math> \begin{bmatrix} c & -s \\ s & c \end{bmatrix} \begin{bmatrix} a \\ b \end{bmatrix} = \begin{bmatrix} r \\ 0 \end{bmatrix} , </math>
where <math> r = \sqrt{a^2 + b^2} </math> is the length of the vector <math>(a,b)</math>.
Explicit calculation of {{mvar|θ}} is rarely necessary or desirable. Instead we directly seek {{mvar|c}} and {{mvar|s}}. An obvious solution would be
:<math>\begin{align}
 c &{}\larr a / r \\
 s &{}\larr -b / r.
\end{align}</math><ref>{{cite book|last1=Björck|first1=Ake|title=Numerical Methods for Least Squares Problems|date=1996|publisher=SIAM|location=United States|isbn=9780898713602|page=54|url=https://books.google.com/books?id=aQD1LLYz6tkC|accessdate=16 August 2016|language=en}}</ref>

However, the computation for {{mvar|r}} may [[arithmetic overflow|overflow]] or underflow. An alternative formulation avoiding this problem {{harv|Golub|Van Loan|1996|loc=§5.1.8}} is implemented as the [[hypot]] function in many programming languages.

The following Fortran code is a minimalistic implementation of Givens rotation for real numbers. If the input values 'a' or 'b' are frequently zero, the code may be optimized to handle these cases as presented [https://dl.acm.org/citation.cfm?doid=567806.567809 here].

<syntaxhighlight lang=fortran>
subroutine givens_rotation(a, b, c, s, r)

real a, b, c, s, r
real h, d

if (b .ne. 0.0) then
    h = hypot(a, b)
    d = 1.0 / h
    c = abs(a) * d
    s = sign(d, a) * b
    r = sign(1.0, a) * h
else
    c = 1.0
    s = 0.0
    r = a
end if

return
end
</syntaxhighlight>


Furthermore, as Edward Anderson discovered in improving [[LAPACK]], a previously overlooked numerical consideration is continuity. To achieve this, we require {{mvar|r}} to be positive.<ref>{{cite web|publisher=University of Tennessee at Knoxville and Oak Ridge National Laboratory|last1=Anderson|first1=Edward|title=Discontinuous Plane Rotations and the Symmetric Eigenvalue Problem|url=http://www.netlib.org/lapack/lawnspdf/lawn150.pdf|accessdate=16 August 2016|date=4 December 2000|series=LAPACK Working Note}}</ref> The following [[MATLAB]]/[[GNU Octave]] code illustrates the algorithm. 

<syntaxhighlight lang="matlab">
function [c, s, r] = givens_rotation(a, b)
    if b == 0;
        c = sign(a);
        if (c == 0);
            c = 1.0; % Unlike other languages, MatLab's sign function returns 0 on input 0.
        end;
        s = 0;
        r = abs(a);
    elseif a == 0;
        c = 0;
        s = -sign(b);
        r = abs(b);
    elseif abs(a) > abs(b);
        t = b / a;
        u = sign(a) * sqrt(1 + t * t);
        c = 1 / u;
        s = -c * t;
        r = a * u;
    else
        t = a / b;
        u = sign(b) * sqrt(1 + t * t);
        s = -1 / u;
        c = t / u;
        r = b * u;
    end
end
</syntaxhighlight>

The [[IEEE 754]] <code>copysign(x,y)</code> function, provides a safe and cheap way to copy the sign of <code>y</code> to <code>x</code>. If that is not available, {{math|{{abs|''x''}}⋅sgn(''y'')}}, using the [[absolute value|abs]] and [[sign function|sgn]] functions, is an alternative as done above.