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