MODULE Giv_Procedures¶

Module Giv_Procedures exports subroutines for computing and applying Givens rotations and reflections [Golub_VanLoan:1996]. Both standard and fast Givens rotations/reflections are implemented in this module.

A Givens rotation is a rotation in the plane acting on two elements of a given vector. Givens rotations are typically used to introduce zeros in vectors, such as during the QR decomposition of a matrix [Golub_VanLoan:1996]. In this case, it is typically desired to find scalars cs and sn such that

$\left( \begin{matrix} cs & sn \\ -sn & cs \end{matrix} \right) \left( \begin{matrix} a \\ b \end{matrix} \right) = \left( \begin{matrix} r \\ 0 \end{matrix} \right)$

where $r = \sqrt{ a^2 + b^2 }$ and $cs^2 + sn^2 = 1$ .

The standard Givens rotations/reflections routines in Giv_Procedures use algorithms and guidelines provided in Section 3.4 of [Anderson_Fahey:1997].

The fast Givens rotations/reflections routines in Giv_Procedures are implementations of the two-way branch algorithms (fast plane rotations with dynamic scaling to avoid overflow/underflow) described in [Anda_Park:1994].

Please note 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 Giv_Procedures or use Statpack statement in your Fortran program, like:

use Giv_Procedures, only: define_rot_givens

or :

use Statpack, only: define_rot_givens

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

define_rot_givens()¶

purpose:

define_rot_givens() generates the cosine and sine of a Givens plane rotation, so that

$\left(\begin{matrix} cs & sn \\ -sn & cs \end{matrix} \right) * \left( \begin{matrix} a \\ b \end{matrix} \right) = \left( \begin{matrix} r \\ 0 \end{matrix} \right)$

where $r = \sqrt{ a^2 + b^2 }$ and $cs^2 + sn^2 = 1$ .

On output, the rotation is also stored in compact form in B and can be recovered by the following algorithm:

If $B = 1$ , set $cs = 0$ and $sn = 1$
If $|B| < 1$ , set $sn = B$ and $cs = \sqrt{ 1 - sn^2 }$
If $|B| > 1$ , set $cs = 1/B$ and $sn = \sqrt{ 1 - cs^2 }$

Synopsis:

call define_rot_givens( a , b , cs , sn )

rot_givens()¶

purpose:

rot_givens() generates and applies a Givens plane rotation to the vector (a b) or to the n-by-2 matrix [VECA VECB], so that

$\left(\begin{matrix} cs & sn \\ -sn & cs \end{matrix} \right) * \left( \begin{matrix} a \\ b \end{matrix} \right) = \left( \begin{matrix} r \\ 0 \end{matrix} \right)$

where $r = \sqrt{ a^2 + b^2 }$ and $cs^2 + sn^2 = 1$ , or

$\left(\begin{matrix} cs & sn \\ -sn & cs \end{matrix} \right) * \left[ \begin{matrix} VECA^T \\ VECB^T \end{matrix} \right] = \left[ \begin{matrix} cs*VECA^T + sn*VECB^T \\ -sn*VECA^T + cs*VECB^T \end{matrix} \right]$

where:

$cs^2 + sn^2 = 1$ ,
$cs*VECA(1) + sn*VECB(1) = \sqrt{ VECA(1)^2 + VECB(1)^2 }$ ,
$-sn*VECA(1) + cs*VECB(1) = 0$ .

On output, the rotation is also stored in compact form in B (or VECB(1)) and can be recovered by the following algorithm:

If $B = 1$ , set $cs = 0$ and $sn = 1$
If $|B| < 1$ , set $sn = B$ and $cs = \sqrt{ 1 - sn^2 }$
If $|B| > 1$ , set $cs = 1/B$ and $sn = \sqrt{ 1 - cs^2 }$

Synopsis:

call rot_givens( a        , b                  )
call rot_givens( veca(:n) , vecb(:n)           )
call rot_givens( a        , b        , cs , sn )
call rot_givens( veca(:n) , vecb(:n) , cs , sn )

apply_rot_givens()¶

purpose:

apply_rot_givens(), eventually reconstructs a Givens plane rotation, stored in compact form in B, and applies this Givens plane rotation to the vector (c d) or to two vectors VECC and VECD.

That is, the value B allows the cosine and sine of the Givens plane rotation to be recovered by the following algorithm:

If $B = 1$ , set $cs = 0$ and $sn = 1$
If $|B| < 1$ , set $sn = B$ and $cs = \sqrt{ 1 - sn^2 }$
If $|B| > 1$ , set $cs = 1/B$ and $sn = \sqrt{ 1 - cs^2 }$

Next, the Givens plane rotation is applied to the vector (c d):

$\left(\begin{matrix} cs & sn \\ -sn & cs \end{matrix} \right) * \left( \begin{matrix} c \\ d \end{matrix} \right) = \left( \begin{matrix} cs*c + sn*d \\ -sn*c + cs*d \end{matrix} \right)$

or to two vectors VECC and VECD:

$\left(\begin{matrix} cs & sn \\ -sn & cs \end{matrix} \right) * \left[ \begin{matrix} VECC^T \\ VECD^T \end{matrix} \right] = \left[ \begin{matrix} cs*VECC^T + sn*VECD^T \\ -sn*VECC^T + cs*VECD^T \end{matrix} \right]$

where $cs^2 + sn^2 = 1$ .

Synopsis:

call apply_rot_givens( c        , d        , b       )
call apply_rot_givens( vecc(:n) , vecd(:n) , b       )
call apply_rot_givens( c        , d        , cs , sn )
call apply_rot_givens( vecc(:n) , vecd(:n) , cs , sn )

givens_vec()¶

purpose:

givens_vec() defines and applies a Givens plane rotation to the n-by-2 matrix [VECA VECB]. The rotation is designed to annihilate the first element of VECB (e.g., VECB(1)). That is,

$\left(\begin{matrix} cs & sn \\ -sn & cs \end{matrix} \right) * \left[ \begin{matrix} VECA^T \\ VECB^T \end{matrix} \right] = \left[ \begin{matrix} cs*VECA^T + sn*VECB^T \\ -sn*VECA^T + cs*VECB^T \end{matrix} \right]$

where:

$cs^2 + sn^2 = 1$ ,
$-sn*VECA(1) + cs*VECB(1) = 0$ ,
$cs*VECA(1) + sn*VECB(1) = \sqrt{ VECA(1)^2 + VECB(1)^2 }$ .

Synopsis:

call givens_vec( veca(:n) , vecb(:n)           )
call givens_vec( veca(:n) , vecb(:n) , cs , sn )

givens_mat_left()¶

purpose:

givens_mat_left() transforms the matrix MAT to upper trapezoidal form by applying a series of Givens plane rotations on the rows of MAT.

Synopsis:

call givens_mat_left( mat(:,:) )

Examples:

ex1_givens_mat_left.F90

givens_mat_right()¶

purpose:

givens_mat_right() transforms the matrix MAT to lower trapezoidal form by applying a series of Givens plane rotations on the columns of MAT.

Synopsis:

call givens_mat_right( mat(:,:) )

Examples:

ex1_givens_mat_right.F90

givens_vec_mat_left()¶

purpose:

givens_vec_mat_left() defines and applies a series of Givens rotations on a n-vector VEC and on the rows of a p-by-n matrix MAT. The rotations are designed to annihilate all the elements of the first column of MAT.

Synopsis:

call givens_vec_mat_left( vec(:n) , mat(:p,:n) )

givens_vec_mat_right()¶

purpose:

givens_vec_mat_right() defines and applies a series of Givens rotations on a n-vector VEC and on the columns of a p-by-n matrix MAT. The rotations are designed to annihilate all the elements of the first row of MAT.

Synopsis:

call givens_vec_mat_right( vec(:p) , mat(:p,:n) )

define_rot_fastgivens()¶

purpose:

define_rot_fastgivens() generates a fast Givens plane rotation H (defined by BETA, ALPHA, and TYPE_ROT on output) and updated scale factors (D1 and D2), which zero X2. That is,

$\left(\begin{matrix} x1 & x2 \end{matrix} \right) * H = \left( \begin{matrix} r & 0 \end{matrix} \right)$

where H is equal to

$\left( \begin{matrix} 1 & 0 \\ 0 & 1 \end{matrix} \right)$ , if TYPE_ROT = 0.

$\left( \begin{matrix} 1 & 0 \\ B & 1 \end{matrix} \right) \left( \begin{matrix} 1 & A \\ 0 & 1 \end{matrix} \right)$ , if TYPE_ROT = 1

$\left( \begin{matrix} 1 & A \\ 0 & 1 \end{matrix} \right) \left( \begin{matrix} 1 & 0 \\ B & 1 \end{matrix} \right)$ , if TYPE_ROT = 2

$\left( \begin{matrix} 0 & -1 \\ 1 & A \end{matrix} \right) \left( \begin{matrix} 1 & 0 \\ -B & 1 \end{matrix} \right)$ , if TYPE_ROT = 3

$\left( \begin{matrix} B & 1 \\ 1 & 0 \end{matrix} \right) \left( \begin{matrix} 1 & A \\ 0 & -1 \end{matrix} \right)$ , if TYPE_ROT = 4

On output, the arguments BETA = B and ALPHA = A and TYPE_ROT define the transformation matrix H:

$H = \left( \begin{matrix} h11 & h12 \\ h21 & h22 \end{matrix} \right)$

Furthermore, if on input, y1 = x1*sqrt(d1) and y2 = x2*sqrt(d2), then on output, with the updated scale factors D1 and D2:

$\left(\begin{matrix} x1 & x2 \end{matrix} \right) * H * diag( \begin{matrix} \sqrt{d1} & \sqrt{d2} \end{matrix} ) = \left( \begin{matrix} [x1*h11 + x2*h21]*\sqrt{d1} & 0 \end{matrix} \right)$

is equal to

$\left(\begin{matrix} y1 & y2 \end{matrix} \right) * \left(\begin{matrix} cs & -sn \\ sn & cs \end{matrix} \right) = \left( \begin{matrix} cs*y1 + sn*y2 & 0 \end{matrix} \right)$

with $cs^2 + sn^2 = 1$ .

In other words, the action of H is equivalent to a standard Givens plane rotation, which zeros y2.

This subroutine is a square root free implementation of the two-way branch algorithm (fast plane rotations with dynamic scaling to avoid overflow/underflow) described in [Anda_Park:1994].

The arguments X1 and X2 are unchanged on return.

Synopsis:

call define_rot_fastgivens( x1 , x2 , d1 , d2 , beta , alpha , type_rot )

apply_rot_fastgivens()¶

purpose:

apply_rot_fastgivens() applies a fast Givens plane rotation H (defined by BETA, ALPHA, and TYPE_ROT on input) to the vector (y1 y2): or to the n-by-2 matrix [VECY1 VECY2]. That is,

$\left(\begin{matrix} y1 & y2 \end{matrix} \right) * H = \left( \begin{matrix} [h11*y1 + h21*y2] & [h12*y1 + h22*y2] \end{matrix} \right)$

or

$[ \begin{matrix} VECY1 & VECY2 \end{matrix} ] * H = [ \begin{matrix} (h11*VECY1 + h21*VECY2) & (h12*VECY1 + h22*VECY2) \end{matrix} ]$

where H is a 2-by-2 matrix defined as

$H = \left( \begin{matrix} h11 & h12 \\ h21 & h22 \end{matrix} \right)$

More precisely, H takes one of the following forms:

$\left( \begin{matrix} 1 & 0 \\ 0 & 1 \end{matrix} \right)$ , if TYPE_ROT = 0.

$\left( \begin{matrix} 1 & 0 \\ B & 1 \end{matrix} \right) \left( \begin{matrix} 1 & A \\ 0 & 1 \end{matrix} \right)$ , if TYPE_ROT = 1

$\left( \begin{matrix} 1 & A \\ 0 & 1 \end{matrix} \right) \left( \begin{matrix} 1 & 0 \\ B & 1 \end{matrix} \right)$ , if TYPE_ROT = 2

$\left( \begin{matrix} 0 & -1 \\ 1 & A \end{matrix} \right) \left( \begin{matrix} 1 & 0 \\ -B & 1 \end{matrix} \right)$ , if TYPE_ROT = 3

$\left( \begin{matrix} B & 1 \\ 1 & 0 \end{matrix} \right) \left( \begin{matrix} 1 & A \\ 0 & -1 \end{matrix} \right)$ , if TYPE_ROT = 4

Synopsis:

call apply_rot_fastgivens( y1,        y2,        beta, alpha, type_rot )
call apply_rot_fastgivens( vecy1(:n), vecy2(:n), beta, alpha, type_rot )

fastgivens_vec()¶

purpose:

fastgivens_vec() generates and applies a fast Givens plane rotation H to the n-by-2 matrix [VECX1 VECX2]. The rotation is designed to zero VECX2(1). That is,

$[ \begin{matrix} VECX1 & VECX2 \end{matrix} ] * H = [ \begin{matrix} (h11*VECX1 + h21*VECX2) & (h12*VECX1 + h22*VECX2) \end{matrix} ]$

where $h12*VECX1(1) + h22*VECX2(1) = 0$ and H is the 2-by-2 matrix:

$H = \left( \begin{matrix} h11 & h12 \\ h21 & h22 \end{matrix} \right)$

Furthermore, the scale factors (D1 and D2) are updated accordingly. That is, if on input:

$[ \begin{matrix} Y1 & Y2 \end{matrix} ] = [ \begin{matrix} \sqrt{d1}*VECX1 & \sqrt{d2}*VECX2 \end{matrix} ]$

then on output:

$[ \begin{matrix} \sqrt{d1}*VECX1 & \sqrt{d2}*VECX2 \end{matrix} ] = [ \begin{matrix} Y1 & Y2 \end{matrix} ] * \left(\begin{matrix} cs & -sn \\ sn & cs \end{matrix} \right) = [ \begin{matrix} (cs*Y1 + sn*Y2) & (-sn*Y1 + cs*Y2) \end{matrix} ]$

with $cs^2 + sn^2 = 1$ and $-sn*Y1(1) + cs*Y2(1) = 0$ .

In other words, the action of H is equivalent to a standard Givens plane rotation, which zeros $Y2(1) = \sqrt{d2}*VECX2(1)$ .

See the subroutine define_rot_fastgivens() for further details on the form of H.

Synopsis:

call fastgivens_vec( vecx1(:n) , vecx2(:n) , d1 , d2                          )
call fastgivens_vec( vecx1(:n) , vecx2(:n) , d1 , d2, beta , alpha , type_rot )

fastgivens_mat_left()¶

purpose:

fastgivens_mat_left() reduces the matrix MAT to upper trapezoidal form by applying a series of fast Givens plane rotations on the rows of MAT.

The (row) scale factors (MATD) are updated accordingly.

Synopsis:

call fastgivens_mat_left( mat(:p,:n), matd(:p) )

Examples:

ex1_fastgivens_mat_left.F90

fastgivens_mat_right()¶

purpose:

fastgivens_mat_right() reduces the matrix MAT to lower trapezoidal form by applying a series of fast Givens plane rotations on the columns of MAT.

The (column) scale factors (MATD) are updated accordingly.

Synopsis:

call fastgivens_mat_right( mat(:p,:n), matd(:n) )

Examples:

ex1_fastgivens_mat_right.F90

fastgivens_vec_mat_left()¶

purpose:

fastgivens_vec_mat_left() defines and applies a series of fast Givens plane rotations on the n-vector VEC and on the rows of a m-by-n matrix MAT. The rotations are designed to annihilate all the elements of the first column of MAT.

The (row) scale factors (VECD and MATD) are updated accordingly.

Synopsis:

call fastgivens_vec_mat_left( vec(:n) , mat(:p,:n), vecd , matd(:p) )

fastgivens_vec_mat_right()¶

purpose:

fastgivens_vec_mat_right() defines and applies a series of fast Givens plane rotations on the data:m-vector VEC and on the columns of a m-by-n matrix MAT. The rotations are designed to annihilate all the elements of the first row of MAT.

The (column) scale factors (VECD and MATD) are updated accordingly.

Synopsis:

call fastgivens_vec_mat_right( vec(:p) , mat(:p,:n) , vecd , matd(:n) )

define_rot_fastgivens2()¶

purpose:

define_rot_fastgivens2() generates a fast Givens plane rotation H (defined by BETA, ALPHA, and TYPE_ROT on output) and updated scale factors (D1 and D2) , which zero X2. That is,

$\left(\begin{matrix} x1 & x2 \end{matrix} \right) * H = \left( \begin{matrix} r & 0 \end{matrix} \right)$

where H is equal to

$\left( \begin{matrix} 1 & 0 \\ 0 & 1 \end{matrix} \right)$ , if TYPE_ROT = 0.

$\left( \begin{matrix} 1 & 0 \\ B & 1 \end{matrix} \right) \left( \begin{matrix} 1 & A \\ 0 & 1 \end{matrix} \right)$ , if TYPE_ROT = 1

$\left( \begin{matrix} 1 & A \\ 0 & 1 \end{matrix} \right) \left( \begin{matrix} 1 & 0 \\ B & 1 \end{matrix} \right)$ , if TYPE_ROT = 2

$\left( \begin{matrix} 0 & -1 \\ 1 & A \end{matrix} \right) \left( \begin{matrix} 1 & 0 \\ -B & 1 \end{matrix} \right)$ , if TYPE_ROT = 3

$\left( \begin{matrix} B & 1 \\ 1 & 0 \end{matrix} \right) \left( \begin{matrix} 1 & A \\ 0 & -1 \end{matrix} \right)$ , if TYPE_ROT = 4

On output, the arguments BETA = B and ALPHA = A and TYPE_ROT define the transformation matrix H:

$H = \left( \begin{matrix} h11 & h12 \\ h21 & h22 \end{matrix} \right)$

Furthermore, if on input, y1 = x1*d1 and y2 = x2*d2, then on output, with the updated scale factors D1 and D2:

$\left(\begin{matrix} x1 & x2 \end{matrix} \right) * H * diag( \begin{matrix} d1 & d2 \end{matrix} ) = \left( \begin{matrix} [x1*h11 + x2*h21]*d1 & 0 \end{matrix} \right)$

is equal to

$\left(\begin{matrix} y1 & y2 \end{matrix} \right) * \left(\begin{matrix} cs & -sn \\ sn & cs \end{matrix} \right) = \left( \begin{matrix} cs*y1 + sn*y2 & 0 \end{matrix} \right)$

with $cs^2 + sn^2 = 1$ .

In other words, the action of H is equivalent to a standard Givens plane rotation, which zeros y2.

This subroutine is an implementation of the two-way branch algorithm (fast plane rotations with dynamic scaling to avoid overflow/underflow) described in [Anda_Park:1994].

The arguments X1 and X2 are unchanged on return.

Synopsis:

call define_rot_fastgivens2( x1 , x2 , d1 , d2 , beta , alpha , type_rot )

fastgivens2_vec()¶

purpose:

fastgivens2_vec() generates and applies a fast Givens plane rotation H to the n-by-2 matrix [VECX1 VECX2]. The rotation is designed to zero VECX2(1). That is,

$[ \begin{matrix} VECX1 & VECX2 \end{matrix} ] * H = [ \begin{matrix} (h11*VECX1 + h21*VECX2) & (h12*VECX1 + h22*VECX2) \end{matrix} ]$

where $h12*VECX1(1) + h22*VECX2(1) = 0$ and H is the 2-by-2 matrix:

$H = \left( \begin{matrix} h11 & h12 \\ h21 & h22 \end{matrix} \right)$

Furthermore, the scale factors (D1 and D2) are updated accordingly. That is, if on input:

$[ \begin{matrix} Y1 & Y2 \end{matrix} ] = [ \begin{matrix} d1*VECX1 & d2*VECX2 \end{matrix} ]$

then on output:

$[ \begin{matrix} d1*VECX1 & d2*VECX2 \end{matrix} ] = [ \begin{matrix} Y1 & Y2 \end{matrix} ] * \left(\begin{matrix} cs & -sn \\ sn & cs \end{matrix} \right) = [ \begin{matrix} (cs*Y1 + sn*Y2) & (-sn*Y1 + cs*Y2) \end{matrix} ]$

with $cs^2 + sn^2 = 1$ and $-sn*Y1(1) + cs*Y2(1) = 0$ .

In other words, the action of H is equivalent to a standard Givens plane rotation, which zeros $Y2(1) = d2*VECX2(1)$ .

See the subroutine define_rot_fastgivens2() for further details on the form of H.

Synopsis:

call fastgivens2_vec( vecx1(:n) , vecx2(:n) , d1 , d2                          )
call fastgivens2_vec( vecx1(:n) , vecx2(:n) , d1 , d2, beta , alpha , type_rot )

Previous topic

Next topic

MODULE Giv_Procedures¶