The LAPACK subroutine to compute the Householder reflection vector (`u`) is DLARFG. To apply the reflector is DLARF.
So to answer your question, please look at DLARFG. This uses the standard numerical linear algebra formula. We do the +/- trick for stability reasons. (So this is either I - 2 (uu')/(u'u) or I + 2 (uu')/(u'u).) See any numerical linear algebra textbook for explanation. Since the scalar 2/(u'u) would need to be computed at each application of the reflection, we prefer to store it. This is `TAU`. So the formula is rewritten as I - TAU (uu'), where TAU = 2/(u'u). Also note that in the formula I - 2 (uu')/(u'u), the scaling of u is not specified. (Any scaling works.) The scaling used by LAPACK is such that the top entry in u is a `1`. (Why not.) So in practice we do not store the top entry of `u`. (Since we know it is `1`.)
For an example of use of DLARFG and DLARF, see DGEQR2.
You may also be interested in looking at DLARFGP. DLARFGP does not use the +/- trick for stability, it uses another trick (so called Parlett's trick).
Then there are subroutines to block the application of `NB` reflections. First to form the `T` matrix is DLARFT, (pre-processing,) and then to apply a sequence of `NB` reflections is DLARFB. This is interesting for performance reasons. For an example of use, please look at DGEQRF.
I am attaching a quick matlab code that reproduces dlarfg.m. Main difference is that: `sqrt ( alpha^2 + xnorm^2 )` is done in a unsafe way in the matlab code below.
Julien.
- Code: Select all
function [tau,y] = dlarfg( y )
% ..
%
% Purpose
% =======
%
% DLARFG generates a real elementary reflector H of order n, such
% that
%
% H * y_in = ( y_out(1) ), H' * H = I.
% ( 0 )
%
% H is represented in the form
%
% H = I - tau * ( 1 ) * ( 1 y_out(2:n-1)' ) ,
% ( y_out(2:n-1) )
%
% If the elements of y_in(2:n-1) are all zero, then tau = 0 and H is taken to be
% the unit matrix.
%
% Otherwise 1 <= tau <= 2.
%
n = size(y,1);
alpha = y(1);
if( n <= 1 ),
tau = 0;
return
end
xnorm = norm( y(2:n,1) ) ;
if( xnorm == 0 ),
tau = 0 ;
else
beta = - sign( alpha ) * sqrt ( alpha^2 + xnorm^2 ) ;
tau = ( beta - alpha ) / beta ;
y( 2:n ) = y( 2:n ) / ( alpha - beta );
y(1) = beta;
end
end