MODULE Hous_Procedures¶

Module Hous_Procedures exports subroutines for computing and applying elementary Householder reflectors [Golub_VanLoan:1996] [Lawson_Hanson:1974].

A Householder transformation is a rank-1 modification of the identity matrix which can be used to zero out selected elements of a vector. A n-by-n Householder reflector matrix H in STATPACK is represented in the form

$H = I + \tau * ( v * v^T )$

where v is a n-vector, called the Householder vector, and $\tau$ is a scalar. Note that a Householder reflector matrix H always verifies

$H^T * H = I$

and

$H = H^T$

The routines in the Hous_Procedures take into account the rank-1 structure and the orthogonal and symmetry properties of Householder reflectors to create and apply Householder transformations efficiently.

Two different implementations of Householder reflectors are provided here, the first is described in [Anderson_Fahey:1997] and the second in [Lawson_Hanson:1974]. In both cases, improvements suggested in [Anderson:2018] and [Hanson_Hopkins:2018] for computing safely and accurately the 2-norm of a vector have also been incorporated.

Please note, finally, that routines provided in this module apply only to real data of kind stnd. The real kind type stnd is defined in module Select_Parameters.

In order to use one of these routines, you must include an appropriate use Hous_Procedures or use Statpack statement in your Fortran program, like:

use Hous_Procedures, only: hous1

or :

use Statpack, only: hous1

Here is the list of the public routines exported by module Hous_Procedures:

hous1()¶

purpose:

Given a n-element real vector x (provided in the argument U), hous1() generates a real elementary reflector H of order n, such that

$H * x = \left( \begin{matrix} beta \\ 0 \end{matrix} \right)$

where beta is a real scalar. H is represented in the form

$H = I + tau * ( v * v^T )$

where tau is a real scalar and v is a n-element real vector with v(1) = 1.

If the elements of x(2:n) are all zero or size(U)=1, then tau = 0 and H is taken to be the unit matrix.

Otherwise 1 <= tau <= 2.

This subroutine is based on the routine DLARFG() in LAPACK with improvements suggested by [Anderson_Fahey:1997] [Anderson:2018] and [Hanson_Hopkins:2018].

Synopsis:

call hous1( u(:n) , tau        )
call hous1( u(:n) , tau , beta )

Examples:

ex1_hous1.F90

apply_hous1()¶

purpose:

apply_hous1() applies a real elementary reflector H generated by hous1() to a real vector/matrix C. H is represented in the form

$H = I + tau * ( v * v^T )$

where tau is a real scalar and v is a n-element real vector with v(1) = 1.

If tau = 0, then H is taken to be the unit matrix and C is not modified.

Synopsis:

call apply_hous1( u(:n) , tau , vec(:)          )
call apply_hous1( u(:n) , tau , mat(:,:) , left )

Examples:

ex1_hous1.F90

hous2()¶

purpose:

Given a n-1-element real vector x and a real scalar alpha (provided in the arguments U and PIVOT on entry), hous2() generates a real elementary reflector H of order n, such that

$H * \left( \begin{matrix} alpha \\ x \end{matrix} \right) = \left( \begin{matrix} beta \\ 0 \end{matrix} \right)$

where beta is a real scalar. H is represented in the form

$H = I + tau * ( v * v^T )$

where tau is a real scalar and v is a n-element real vector with v(1) = 1.

If the elements of x are all zero, then tau = 0 and H is taken to be the unit matrix.

Otherwise 1 <= tau <= 2.

This subroutine is based on the routine DLARFG() in LAPACK with improvements suggested by [Anderson_Fahey:1997] [Anderson:2018] and [Hanson_Hopkins:2018].

Synopsis:

call hous2( pivot, u(:n-1) , tau  )

Examples:

ex1_hous2.F90

apply_hous2()¶

purpose:

apply_hous2() applies a real elementary reflector H of order n, generated by hous2(), to a real vector/matrix C. H is represented in the form

$H = I + tau * ( v * v^T )$

where tau is a real scalar and v is a n-element real vector with v(1) = 1. More precisely, v is defined as

$v = \left( \begin{matrix} 1 \\ u \end{matrix} \right)$

for the input vector argument U.

Here, the n-element real vector C has the form:

$C = \left( \begin{matrix} piv \\ vec \end{matrix} \right)$

or the real matrix C has the form, if LEFT=true:

$C = \left[ \begin{matrix} vec\_piv^T \\ MAT \end{matrix} \right]$

or, if LEFT=false:

$C = \left[ \begin{matrix} vec\_piv & MAT \end{matrix} \right]$

for given input arguments PIV, VEC or VEC_PIV, MAT.

If tau = 0, then H is taken to be the unit matrix and C is not modified.

Synopsis:

call apply_hous2( u(:n-1) , tau ,  piv        , vec(:)          )
call apply_hous2( u(:n-1) , tau ,  vec_piv(:) , mat(:,:) , left )

Examples:

ex1_hous2.F90

h1()¶

purpose:

Given a n-element real vector x (provided in the argument U), h1() generates a real elementary reflector H of order n, such that

$H * x = \left( \begin{matrix} beta \\ 0 \end{matrix} \right)$

where beta is a real scalar. H is represented in the form

$H = I + tau * ( v * v^T )$

where tau is a real scalar and v is a n-element real vector. Note that here, v(1) is not equal to 1, contrary to the formulation used in hous1().

The real elementary reflector H is then, optionally, applied to a real vector/matrix C.

This subroutine is based on the routine H1 described in [Lawson_Hanson:1974] with improvements suggested in [Anderson:2018] and [Hanson_Hopkins:2018].

Synopsis:

call h1( u(:n) , beta , tau                   )
call h1( u(:n) , beta , tau , vec(:)          )
call h1( u(:n) , beta , tau , mat(:,:) , left )

Examples:

ex1_h1.F90

apply_h1()¶

purpose:

apply_h1() applies a real elementary reflector H generated by h1() to a n-element real vector or to a n-by-m or m-by-n real matrix from the left or the right. H is represented in the form

$H = I + tau * ( v * v^T )$

where tau is a real scalar and v is a n-element real vector.

Synopsis:

call apply_h1( u(:n) , tau , vec(:)          )
call apply_h1( u(:n) , tau , mat(:,:) , left )

Examples:

ex1_h1.F90

h2()¶

purpose:

Given a n-1-element real vector x and a real scalar alpha (provided in the arguments U and BETA on entry), h2() generates a real elementary reflector H of order n, such that

$H * \left( \begin{matrix} alpha \\ x \end{matrix} \right) = \left( \begin{matrix} beta \\ 0 \end{matrix} \right)$

where beta is a real scalar. H is represented in the form

$H = I + tau * ( v * v^T )$

where tau is a real scalar and v is a n-element real vector. Note that here v(1) is not equal to 1, contrary to the formulation in hous2().

The real elementary reflector H is then, optionally, applied to a real vector/matrix C.

This subroutine is based on the routine H2 described in [Lawson_Hanson:1974] with improvements suggested in [Anderson:2018] and [Hanson_Hopkins:2018].

Synopsis:

call h2( beta, u(:n-1) , up , tau                              )
call h2( beta, u(:n-1) , up , tau , piv        , vec(:)        )
call h2( beta, u(:n-1) , up , tau , vec_piv(:) , mat(:) , left )

Examples:

ex1_h2.F90

apply_h2()¶

purpose:

apply_h2() applies a real elementary reflector H generated by h2() to a n-element real vector or to a n-by-m or m-by-n real matrix from the left or the right. H is represented in the form

$H = I + tau * ( v * v^T )$

where tau is a real scalar and v is a n-element real vector. More precisely, v is defined as

$v = \left( \begin{matrix} up \\ u \end{matrix} \right)$

from the input arguments UP and U.

Here, the n-element real vector C has the form:

$C = \left( \begin{matrix} piv \\ vec \end{matrix} \right)$

or the real matrix C has the form, if LEFT=true:

$C = \left[ \begin{matrix} vec\_piv^T \\ MAT \end{matrix} \right]$

or, if LEFT=false:

$C = \left[ \begin{matrix} vec\_piv & MAT \end{matrix} \right]$

for given input arguments PIV, VEC or VEC_PIV, MAT.

If tau = 0, then H is taken to be the unit matrix and C is not modified.

Synopsis:

call apply_h2( u(:n-1) , up , tau ,  piv        , vec(:)          )
call apply_h2( u(:n-1) , up , tau ,  vec_piv(:) , mat(:,:) , left )

Examples:

ex1_h2.F90

Previous topic

Next topic

MODULE Hous_Procedures¶