Module_Giv_Procedures¶
Copyright 2024 IRD
This file is part of statpack.
statpack is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
statpack is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.
You can find a copy of the GNU Lesser General Public License in the statpack/doc directory.
Authors: Pascal Terray (LOCEAN/IPSL, Paris, France)
MODULE EXPORTING GIVENS TOOLS (REFLECTIONS AND ROTATIONS).
LATEST REVISION : 11/03/2024
subroutine define_rot_givens ( a, b, cs, sn )¶
Purpose¶
DEFINE_ROT_GIVENS generates the cosine and sine of a Givens plane rotation, ROT, so that
( A B ) ROT = ( R 0 )where R >= 0 and ROT is 2-by-2 matrix defined by
( +CS -SN )
( +SN +CS )
with CS**(2) + SN**(2) = 1.
Arguments¶
- A (INPUT) real(stnd)
- The first component of vector to be rotated.
- B (INPUT) real(stnd)
- The second component of vector to be rotated.
- CS (OUTPUT) real(stnd)
- The cosine of the rotation.
- SN (OUTPUT) real(stnd)
- The sine of the rotation.
Further Details¶
A and B are unchanged on return.
Normally, the subprogram APPLY_ROT_GIVENS( VECA, VECB, CS, SN ) will next be called to apply the rotation to a n-by-2 matrix [ VECA VECB ].
subroutine rot_givens ( a, b )¶
Purpose¶
ROT_GIVENS applies a Givens plane rotation, ROT, so that
( A B ) ROT = ( R 0 )where ROT is 2-by-2 matrix defined by
( +CS -SN )
( +SN +CS )
with CS**(2) + SN**(2) = 1.
On output, the rotation is also stored in compact form in B.
Arguments¶
- A (INPUT/OUTPUT) real(stnd)
The first component of vector to be rotated.
On output, R = (+/-)sqrt( A**(2) + B**(2) ) overwrites A.
- B (INPUT/OUTPUT) real(stnd)
The second component of vector to be rotated.
On output, Z overwrites B. Z allows CS and SN to be recovered by the following algorithm:
- If Z = 1 set CS = 0 and SN = 1
- If abs( Z ) < 1 set SN = Z and CS = sqrt( 1 - SN**(2) )
- If abs( Z ) > 1 set CS = 1/Z and SN = sqrt( 1 - CS**(2) )
Further Details¶
Normally, the subprogram APPLY_ROT_GIVENS( VECA, VECB, B ) will next be called to apply the rotation to a n-by-2 matrix [ VECA VECB ].
subroutine rot_givens ( a, b, cs, sn )¶
Purpose¶
ROT_GIVENS generates and applies a Givens plane rotation, ROT, so that
( A B ) ROT = ( R 0 )where ROT is 2-by-2 matrix defined by
( +CS -SN )
( +SN +CS )
with CS**(2) + SN**(2) = 1.
On output, the rotation is also stored in compact form in B.
Arguments¶
- A (INPUT/OUTPUT) real(stnd)
The first component of vector to be rotated.
On output, R = (+/-)sqrt( A**(2) + B**(2) ) overwrites A.
- B (INPUT/OUTPUT) real(stnd)
The second component of vector to be rotated.
On output, Z overwrites B. Z allows CS and SN to be recovered by the following algorithm:
- If Z = 1 set CS = 0 and SN = 1
- If abs( Z ) < 1 set SN = Z and CS = sqrt( 1 - SN**(2) )
- If abs( Z ) > 1 set CS = 1/Z and SN = sqrt( 1 - CS**(2) )
- CS (OUTPUT) real(stnd)
- The cosine of the rotation.
- SN (OUTPUT) real(stnd)
- The sine of the rotation.
Further Details¶
Normally, the subprograms APPLY_ROT_GIVENS( VECA, VECB, CS, SN ) or APPLY_ROT_GIVENS( VECA, VECB, B ) will next be called to apply the rotation to a n-by-2 matrix [ VECA VECB ].
subroutine rot_givens ( veca, vecb )¶
Purpose¶
ROT_GIVENS applies a Givens plane rotation, ROT, to the n-by-2 matrix [ VECA VECB ]. The rotation is designed to annilhate the first element of VECB ( e.g. VECB(1)). That is,
[ VECA VECB ] ROT = [ (CS*VECA + SN*VECB) (-SN*VECA + CS*VECB) ]where
CS**(2) + SN**(2) = 1,
-SN*VECA(1) + CS*VECB(1) = 0
and ROT is a 2-by-2 matrix defined by
( +CS -SN )
( +SN +CS )
On output, the rotation is also stored in compact form in VECB(1).
Arguments¶
- VECA (INPUT/OUTPUT) real(stnd), dimension(:)
The first vector to be rotated.
On output, CS*VECA + SN*VECB overwrites VECA.
- VECB (INPUT/OUTPUT) real(stnd), dimension(:)
The second vector to be rotated.
On output, -SN*VECA(2:) + CS*VECB(2:) overwrites VECB(2:) and Z overwrites VECB(1). Z allows CS and SN to be recovered by the following algorithm:
- If Z = 1 set CS = 0 and SN = 1
- If abs( Z ) < 1 set SN = Z and CS = sqrt( 1 - SN**(2) )
- If abs( Z ) > 1 set CS = 1/Z and SN = sqrt( 1 - CS**(2) )
Further Details¶
It is assumed that VECA and VECB have the same size.
The subprograms APPLY_ROT_GIVENS( VECC, VECD, VECB(1) ) may next be called to apply the rotation to another n-by-2 matrix [ VECC VECD ].
subroutine rot_givens ( veca, vecb, cs, sn )¶
Purpose¶
ROT_GIVENS defines and applies a Givens plane rotation, ROT, to the n-by-2 matrix [ VECA VECB ]. The rotation is designed to annilhate the first element of VECB (e.g. VECB(1)). That is,
[ VECA VECB ] ROT = [ (CS*VECA + SN*VECB) (-SN*VECA + CS*VECB) ]where
CS**(2) + SN**(2) = 1,
-SN*VECA(1) + CS*VECB(1) = 0
and ROT is a 2-by-2 matrix defined by
( +CS -SN )
( +SN +CS )
On output, the rotation is also stored in compact form in VECB(1).
Arguments¶
- VECA (INPUT/OUTPUT) real(stnd), dimension(:)
The first vector to be rotated.
On output, CS*VECA + SN*VECB overwrites VECA.
- VECB (INPUT/OUTPUT) real(stnd), dimension(:)
The second vector to be rotated.
On output, -SN*VECA(2:) + CS*VECB(2:) overwrites VECB(2:) and Z overwrites VECB(1). Z allows CS and SN to be recovered by the following algorithm:
- If Z = 1 set CS = 0 and SN = 1
- If abs( Z ) < 1 set SN = Z and CS = sqrt( 1 - SN**(2) )
- If abs( Z ) > 1 set CS = 1/Z and SN = sqrt( 1 - CS**(2) )
- CS (OUTPUT) real(stnd)
- The cosine of the rotation.
- SN (OUTPUT) real(stnd)
- The sine of the rotation.
Further Details¶
It is assumed that VECA and VECB have the same size.
Normally, the subprograms APPLY_ROT_GIVENS( VECC, VECD, CS, SN ) or APPLY_ROT_GIVENS( VECC, VECD, VECB(1) ) will next be called to apply the rotation to a n-by-2 matrix [ VECC VECD ].
subroutine apply_rot_givens ( c, d, b )¶
Purpose¶
APPLY_ROT_GIVENS reconstructs and applies a Givens plane rotation, ROT, stored in compact form in B, to the vector ( C D ).
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 abs( B ) < 1 set SN = B and CS = sqrt( 1 - SN**(2) )
- If abs( B ) > 1 set CS = 1/B and SN = sqrt( 1 - CS**(2) )
Next, the Givens plane rotation, ROT, is applied to the vector ( C D ) :
( C D ) ROT = ( (CS*C + SN*D) (-SN*C + CS*D) )where ROT is a 2-by-2 matrix defined by
( +CS -SN )
( +SN +CS )
Arguments¶
- C (INPUT/OUTPUT) real(stnd)
The first element of vector to be rotated.
On output, CS*C + SN*D overwrites C.
- D (INPUT/OUTPUT) real(stnd)
The second element of vector to be rotated.
On output, -SN*C + CS*D overwrites D.
- B (INPUT) real(stnd)
- The real number, which allows the cosine and sine of the Givens plane rotation to be recovered.
Further Details¶
Normally:
- the subprogram APPLY_ROT_GIVENS( C, D, B ) is called to apply a Givens rotation to the vector ( C D ) after a call to ROT_GIVENS( A, B, CS, SN ) or ROT_GIVENS( A, B ).
- the subprogram APPLY_ROT_GIVENS( C, D, VECB(1) ) is called to apply a Givens rotation to the vector ( C D ) after a call to ROT_GIVENS( VECA, VECB, CS, SN ) or ROT_GIVENS( VECA, VECB ) where VECA and VECB are two vectors of the same length.
subroutine apply_rot_givens ( c, d, cs, sn )¶
Purpose¶
APPLY_ROT_GIVENS applies a Givens plane rotation, ROT, to to the vector ( C D ). That is,
( C D ) ROT = ( (CS*C + SN*D) (-SN*C + CS*D) )where ROT is a 2-by-2 matrix defined by
( +CS -SN )
( +SN +CS )
Arguments¶
- C (INPUT/OUTPUT) real(stnd)
The first element of vector to be rotated.
On output, CS*C + SN*D overwrites C.
- D (INPUT/OUTPUT) real(stnd)
The second element of vector to be rotated.
On output, -SN*C + CS*D overwrites D.
- CS (INPUT) real(stnd)
- The cosine of the rotation.
- SN (INPUT) real(stnd)
- The sine of the rotation.
Further Details¶
Normally, the subprogram APPLY_ROT_GIVENS( C, D, CS, SN ) is called to apply a Givens rotation to the vector ( C D) after a call to DEFINE_ROT_GIVENS( A, B, CS, SN ), ROT_GIVENS( A, B, CS, SN ) or ROT_GIVENS( VECA, VECB, CS, SN ).
subroutine apply_rot_givens ( vecc, vecd, b )¶
Purpose¶
APPLY_ROT_GIVENS reconstructs and applies a Givens plane rotation, ROT, stored in compact form in B, to the n-by-2 matrix [ VECC 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 abs( B ) < 1 set SN = B and CS = sqrt( 1 - SN**(2) )
- If abs( B ) > 1 set CS = 1/B and SN = sqrt( 1 - CS**(2) )
Next, the Givens plane rotation, ROT, is applied to the n-by-2 matrix [ VECC VECD ]:
[ VECC VECD ] ROT = [ (CS*VECC + SN*VECD) (-SN*VECC + CS*VECD) ]where ROT is a 2-by-2 matrix defined by
( +CS -SN )
( +SN +CS )
Arguments¶
- VECC (INPUT/OUTPUT) real(stnd), dimension(:)
The first vector to be rotated.
On output, CS*VECC + SN*VECD overwrites VECC.
- VECD (INPUT/OUTPUT) real(stnd), dimension(:)
The second vector to be rotated.
On output, -SN*VECC + CS*VECD overwrites VECD.
- B (INPUT) real(stnd)
- The real number which allows the cosine and sine of the Givens plane rotation to be recovered.
Further Details¶
Normally, the subprogram APPLY_ROT_GIVENS( VECC, VECD, B ) is called to apply a Givens rotation to the n-by-2 matrix [ VECC VECD ] after a call to ROT_GIVENS( A, B, CS, SN ) or ROT_GIVENS( A, B ).
Normally, the subprogram APPLY_ROT_GIVENS( VECC, VECD, VECB(1) ) is called to apply a Givens rotation to the n-by-2 matrix [ VECC VECD ] after a call to ROT_GIVENS( VECA, VECB, CS,SN ) or ROT_GIVENS( VECA, VECB ) where VECA and VECB are two vectors of the same length.
It is assumed that VECC and VECD have the same size.
subroutine apply_rot_givens ( vecc, vecd, cs, sn )¶
Purpose¶
APPLY_ROT_GIVENS applies a Givens plane rotation, ROT, to the n-by-2 matrix [ VECC VECD ]. That is,
[ VECC VECD ] ROT = [ (CS*VECC + SN*VECD) (-SN*VECC + CS*VECD) ]where ROT is a 2-by-2 matrix defined by
( +CS -SN )
( +SN +CS )
Arguments¶
- VECC (INPUT/OUTPUT) real(stnd), dimension(:)
The first vector to be rotated.
On output, CS*VECC + SN*VECD overwrites VECC.
- VECD (INPUT/OUTPUT) real(stnd), dimension(:)
The second vector to be rotated.
On output, -SN*VECC + CS*VECD overwrites VECD.
- CS (INPUT) real(stnd)
- The cosine of the rotation.
- SN (INPUT) real(stnd)
- The sine of the rotation.
Further Details¶
Normally, the subprogram APPLY_ROT_GIVENS( VECC, VECD, CS, SN ) is called to apply a Givens rotation to the n-by-2 matrix [ VECC VECD ] after a call to DEFINE_ROT_GIVENS( A, B, CS, SN ), ROT_GIVENS( A, B, CS, SN ) or ROT_GIVENS( VECA, VECB, CS, SN ).
It is assumed that VECC and VECD have the same size.
subroutine givens_vec ( veca, vecb )¶
Purpose¶
GIVENS defines and applies a Givens plane rotation, ROT, to the n-by-2 matrix [ VECA VECB ]. The rotation is designed to annilhate the first element of VECB (e.g. VECB(1)). That is,
[ VECA VECB ] ROT = [ (CS*VECA + SN*VECB) (-SN*VECA + CS*VECB) ]where:
- CS**(2) + SN**(2) = 1,
- -SN*VECA(1) + CS*VECB(1) = 0,
- CS*VECA(1) + SN*VECB(1) >= 0.
and ROT is a 2-by-2 matrix defined by
( +CS -SN )
( +SN +CS )
Arguments¶
- VECA (INPUT/OUTPUT) real(stnd), dimension(:)
The first vector to be rotated.
On output, CS*VECA + SN*VECB overwrites VECA.
- VECB (INPUT/OUTPUT) real(stnd), dimension(:)
The second vector to be rotated.
On output, -SN*VECA + CS*VECB overwrites VECB.
Further Details¶
It is assumed that VECA and VECB have the same size.
subroutine givens_vec ( veca, vecb, cs, sn )¶
Purpose¶
GIVENS defines and applies a Givens plane rotation, ROT, to the n-by-2 matrix [ VECA VECB ]. The rotation is designed to annilhate the first element of VECB (e.g. VECB(1)). That is,
[ VECA VECB ] ROT = [ (CS*VECA + SN*VECB) (-SN*VECA + CS*VECB) ]where:
- CS**(2) + SN**(2) = 1,
- -SN*VECA(1) + CS*VECB(1) = 0,
- CS*VECA(1) + SN*VECB(1) >= 0.
and ROT is a 2-by-2 matrix defined by
( +CS -SN )
( +SN +CS )
Arguments¶
- VECA (INPUT/OUTPUT) real(stnd), dimension(:)
The first vector to be rotated.
On output, CS*VECA + SN*VECB overwrites VECA.
- VECB (INPUT/OUTPUT) real(stnd), dimension(:)
The second vector to be rotated.
On output, -SN*VECA + CS*VECB overwrites VECB.
- CS (OUTPUT) real(stnd)
- The cosine of the rotation.
- SN (OUTPUT) real(stnd)
- The sine of the rotation.
Further Details¶
It is assumed that VECA and VECB have the same size.
subroutine givens_mat_left ( mat )¶
Purpose¶
GIVENS_MAT_LEFT transforms the matrix MAT to upper trapezoidal form by applying a serie of Givens plane rotations on the rows of MAT.
Arguments¶
- MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)
The matrix to be transformed.
On output, the transformed matrix overwrites MAT.
subroutine givens_mat_right ( mat )¶
Purpose¶
GIVENS_MAT_RIGHT transforms the matrix MAT to lower trapezoidal form by applying a serie of Givens plane rotations on the columns of MAT.
Arguments¶
- MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)
The matrix to be transformed.
On output, the transformed matrix overwrites MAT.
subroutine givens_vec_mat_left ( vec, mat )¶
Purpose¶
GIVENS_VEC_MAT_LEFT defines and applies a serie of Givens rotations on a n-vector VEC and on the rows of a p-by-n matrix MAT. The rotations are designed to annilhate all the elements of the first column of MAT.
Arguments¶
- VEC (INPUT/OUTPUT) real(stnd), dimension(:)
On input, the n-vector to rotate. VEC(1) is used to define the rotations.
On output, the transformed vector overwrites VEC.
- MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)
On input, the matrix to be transformed. MAT(:,1) is used to define the rotations.
On output, the transformed matrix overwrites MAT and MAT(:,1) is equal to zero.
Further Details¶
It is assumed that size( VEC ) = size( MAT, 2 ) .
subroutine givens_vec_mat_right ( vec, mat )¶
Purpose¶
GIVENS_VEC_MAT_RIGHT defines and applies a serie of Givens rotations on a n-vector VEC and on the columns of a n-by-p matrix MAT. The rotations are designed to annilhate all the elements of the first row of MAT.
Arguments¶
- VEC (INPUT/OUTPUT) real(stnd), dimension(:)
On input, the n-vector to rotate. VEC(1) is used to define the rotations.
On output, the transformed vector overwrites VEC.
- MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)
On input, the matrix to be transformed. MAT(1,:) is used to define the rotations.
On output, the transformed matrix overwrites MAT and MAT(1,:) is equal to zero.
Further Details¶
It is assumed that size( VEC ) = size( MAT, 1 ) .
subroutine define_rot_fastgivens ( x1, x2, d1, d2, beta, alpha, type_rot )¶
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,
( X1 X2 ) H = ( R 0 ), where H is equal to
(1 0) ,if TYPE_ROT = 0.
(0 1)
(1 0) (1 A) ,if TYPE_ROT = 1.
(B 1) (0 1)
(1 A) (1 0) ,if TYPE_ROT = 2.
(0 1) (B 1)
(0 -1) (1 0) ,if TYPE_ROT = 3.
(1 A) (-B 1)
(B 1) (1 A) ,if TYPE_ROT = 4.
(1 0) (0 -1)
Furtermore, if on input, Y1 = X1*SQRT(D1) and Y2 = X2*SQRT(D2), then on output
( X1 X2 ) H diag( SQRT(D1) SQRT(D2) ) = ( (X1*H11 + X2*H21)*SQRT(D1) 0 )is equal to
( Y1 Y2 ) ROT = ( (Y1*CS + Y2*SN) 0 )where CS and SN define a standard Givens plane rotation, ROT, which zeros Y2. Thus, ROT is equal to
( +CS -SN )
( +SN +CS )
with CS**(2) + SN**(2) = 1.
Arguments¶
- X1 (INPUT) real(stnd)
- First component of vector to be transformed.
- X2 (INPUT) real(stnd)
- Second component of vector to be transformed.
- D1 (INPUT/OUTPUT) real(stnd)
On input, first scale factor.
On output, D1 is replaced with the update scale factor.
- D2 (INPUT/OUTPUT) real(stnd)
On input, second scale factor.
On output, D2 is replaced with the update scale factor.
- BETA (OUTPUT) real(stnd)
- The real scalar B which defines the transformation matrix H.
- ALPHA (OUTPUT) real(stnd)
- The real scalar A which defines the transformation matrix H.
- TYPE_ROT (OUTPUT) integer(i2b)
- Integer which defines the transformation matrix H.
Further Details¶
X1 and X2 are unchanged on return.
It is assumed that D1 and D2 are positive scalars.
IF D1>=D2, D1 is diminished and D2 is augmented. IF D1<D2, D2 is diminished and D1 is augmented. The decrease or increase in magnitude of D1 and D2 are bounded by 1/2 and 2, respectively.
Normally, the subprogram APPLY_ROT_FASTGIVENS( VECX1, VECX2, BETA, ALPHA, TYPE_ROT ) will next be called to apply the rotation to a n-by-2 matrix [ VECX1 VECX2 ].
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 reference (1).
For further details, see:
- Anda, A.A. and Park, H., 1994:
- Fast plane rotations with dynamic scaling. Siam J. Matrix Anal. Appl., 15, 162-174.
subroutine define_rot_fastgivens2 ( x1, x2, d1, d2, beta, alpha, type_rot )¶
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,
( X1 X2 ) H = ( R 0 ), where H is equal to
(1 0) ,if TYPE_ROT = 0.
(0 1)
(1 0) (1 A) ,if TYPE_ROT = 1.
(B 1) (0 1)
(1 A) (1 0) ,if TYPE_ROT = 2.
(0 1) (B 1)
(0 -1) (1 0) ,if TYPE_ROT = 3.
(1 A) (-B 1)
(B 1) (1 A) ,if TYPE_ROT = 4.
(1 0) (0 -1)
Furtermore, if on input, Y1 = X1*D1 and Y2 = X2*D2, then on output
( X1 X2 ) H diag( D1 D2 ) = ( (X1*H11 + X2*H21)*D1 0 )is equal to
( Y1 Y2 ) ROT = ( (Y1*CS + Y2*SN) 0 )where CS and SN define a standard Givens plane rotation, ROT, which zeros Y2. Thus, ROT is equal to
( +CS -SN )
( +SN +CS )
with CS**(2) + SN**(2) = 1.
Arguments¶
- X1 (INPUT) real(stnd)
- First component of vector to be transformed.
- X2 (INPUT) real(stnd)
- Second component of vector to be transformed.
- D1 (INPUT/OUTPUT) real(stnd)
On input, first scale factor.
On output, D1 is replaced with the update scale factor.
- D2 (INPUT/OUTPUT) real(stnd)
On input, second scale factor.
On output, D2 is replaced with the update scale factor.
- BETA (OUTPUT) real(stnd)
- The real scalar B which defines the transformation matrix H.
- ALPHA (OUTPUT) real(stnd)
- The real scalar A which defines the transformation matrix H.
- TYPE_ROT (OUTPUT) integer(i2b)
- Integer which defines the transformation matrix H.
Further Details¶
X1 and X2 are unchanged on return.
It is assumed that D1 and D2 are positive scalars.
IF D1>=D2, D1 is diminished and D2 is augmented. IF D1<D2, D2 is diminished and D1 is augmented. The decrease or increase in magnitude of D1 and D2 are bounded by 1/sqrt(2) and sqrt(2), respectively.
Normally, the subprogram APPLY_ROT_FASTGIVENS( VECX1, VECX2, BETA, ALPHA, TYPE_ROT ) will next be called to apply the rotation to a n-by-2 matrix [ VECX1 VECX2 ].
This subroutine is an implementation of the two-way branch algorithm (fast plane rotations with dynamic scaling to avoid overflow/underflow) described in reference (1).
For further details, see:
- Anda, A.A. and Park, H., 1994:
- Fast plane rotations with dynamic scaling. Siam J. Matrix Anal. Appl., 15, 162-174.
subroutine apply_rot_fastgivens ( y1, y2, beta, alpha, type_rot )¶
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 ). That is,
( Y1 Y2 ) H = ( (Y1*H11 + Y2*H21) (Y1*H12 + Y2*H22) )where H is a 2-by-2 matrix defined as
( H11 H12 )
( H21 H22 )
More precisely, H takes one of the following forms:
(1 0) ,if TYPE_ROT = 0.
(0 1)
(1 0) (1 A) ,if TYPE_ROT = 1.
(B 1) (0 1)
(1 A) (1 0) ,if TYPE_ROT = 2.
(0 1) (B 1)
(0 -1) (1 0) ,if TYPE_ROT = 3.
(1 A) (-B 1)
(B 1) (1 A) ,if TYPE_ROT = 4.
(1 0) (0 -1)
Arguments¶
- Y1 (INPUT/OUTPUT) real(stnd)
The first component of vector to be transformed.
On output, Y1*H11 + Y2*H21 overwrites Y1.
- Y2 (INPUT/OUTPUT) real(stnd)
The second component of vector to be transformed.
On output, Y1*H12 + Y2*H22 overwrites Y2.
- BETA (OUTPUT) real(stnd)
- The real scalar B which defines the transformation matrix H.
- ALPHA (OUTPUT) real(stnd)
- The real scalar A which defines the transformation matrix H.
- TYPE_ROT (INPUT) integer(i2b)
- Integer which defines the transformation matrix H.
Further Details¶
Normally, the subprogram APPLY_ROT_FASTGIVENS( Y1, Y2, BETA, ALPHA, TYPE_ROT ) will be called to apply the transformation to the vector ( Y1 Y2 ) after a call to DEFINE_ROT_FASTGIVENS( X1, X2, BETA, ALPHA, TYPE_ROT ) or DEFINE_ROT_FASTGIVENS2( X1, X2, BETA, ALPHA, TYPE_ROT ).
subroutine apply_rot_fastgivens ( vecy1, vecy2, beta, alpha, type_rot )¶
Purpose¶
APPLY_ROT_FASTGIVENS applies a fast Givens plane rotation H (defined by BETA, ALPHA, and TYPE_ROT on input) to the n-by-2 matrix [ VECY1 VECY2 ]. That is,
[ VECY1 VECY2 ] H = [ (VECY1*H11 + VECY2*H21) (VECY1*H12 + VECY2*H22) ]where H is a 2-by-2 matrix defined as
( H11 H12 )
( H21 H22 )
More precisely, H takes one of the following forms:
(1 0) ,if TYPE_ROT = 0.
(0 1)
(1 0) (1 A) ,if TYPE_ROT = 1.
(B 1) (0 1)
(1 A) (1 0) ,if TYPE_ROT = 2.
(0 1) (B 1)
(0 -1) (1 0) ,if TYPE_ROT = 3.
(1 A) (-B 1)
(B 1) (1 A) ,if TYPE_ROT = 4.
(1 0) (0 -1)
Arguments¶
- VECY1 (INPUT/OUTPUT) real(stnd), dimension(:)
The first vector to be transformed.
On output, VECY1*H11 + VECY2*H21 overwrites VECY1.
- VECY2 (INPUT/OUTPUT) real(stnd), dimension(:)
The second vector to be transformed.
On output, VECY1*H12 + VECY2*H22 overwrites VECY2.
- BETA (OUTPUT) real(stnd)
- The real scalar B which defines the transformation matrix H.
- ALPHA (OUTPUT) real(stnd)
- The real scalar A which defines the transformation matrix H.
- TYPE_ROT (INPUT) integer(i2b)
- Integer which defines the transformation matrix H.
Further Details¶
Normally, the subprogram APPLY_ROT_FASTGIVENS( VECY1, VECY2, BETA, ALPHA, TYPE_ROT ) will be called to apply the transformation to the n-by-2 matrix [ VECY1 VECY2 ] after a call to DEFINE_ROT_FASTGIVENS( VECY1(1), VECY2(1), BETA, ALPHA, TYPE_ROT ) or DEFINE_ROT_FASTGIVENS2( VECY1(1), VECY2(1), BETA, ALPHA, TYPE_ROT ).
It is assumed that VECY1 and VECY2 have the same size.
subroutine fastgivens_vec ( vecx1, vecx2, d1, d2 )¶
Purpose¶
FASTGIVENS 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,
[ VECX1 VECX2 ] H = [ (VECX1*H11 + VECX2*H21) (VECX1*H12 + VECX2*H22) ], where VECX1(1)*H12 + VECX2(1)*H22 = 0 and H is the 2-by-2 matrix:
( H11 H12 )
( H21 H22 )
Furthermore, the scale factors (D1 and D2) are updated accordingly. That is, if on input:
[ Y1 Y2 ] = [ VECX1 VECX2 ] diag( SQRT(D1) SQRT(D2) )then on output:
[ VECX1 VECX2 ] diag( SQRT(D1) SQRT(D2) ) = [ Y1 Y2 ] ROT = [ (Y1*CS + Y2*SN) (-SN*Y1 + CS*Y2) ]where CS and SN define a standard Givens 2-by-2 plane rotation, ROT, which zeros -SN*Y1(1) + CS*Y2(1). Thus, ROT has the following structure:
( +CS -SN )
( +SN +CS )
with CS**(2) + SN**(2) = 1.
Arguments¶
- VECX1 (INPUT/OUTPUT) real(stnd), dimension(:)
The first vector to be transformed.
On output, VECX1*H11 + VECX2*H21 overwrites VECX1.
- VECX2 (INPUT/OUTPUT) real(stnd), dimension(:)
The second vector to be transformed.
On output, VECX1*H12 + VECX2*H22 overwrites VECX2.
- D1 (INPUT/OUTPUT) real(stnd)
On input, first scale factor.
On output, D1 is replaced with the update scale factor.
- D2 (INPUT/OUTPUT) real(stnd)
On input, second scale factor.
On output, D2 is replaced with the update scale factor.
Further Details¶
It is assumed that D1 and D2 are positive scalars and that VECX1 and VECX2 have the same size.
IF D1>=D2, D1 is diminished and D2 is augmented. IF D1<D2, D2 is diminished and D1 is augmented. The decrease or increase in magnitude of D1 and D2 are bounded by 1/2 and 2, respectively.
This subroutine is a square root free implementation of the two-way branch algorithm (e.g. fast plane rotations with dynamic scaling to avoid overflow/underflow) described in reference (1).
For further details, see:
- Anda, A.A. and Park, H., 1994:
- Fast plane rotations with dynamic scaling. Siam J. Matrix Anal. Appl., 15, 162-174.
subroutine fastgivens_vec ( vecx1, vecx2, d1, d2, beta, alpha, type_rot )¶
Purpose¶
FASTGIVENS generates and applies a fast Givens plane rotation H (defined by BETA, ALPHA and TYPE_ROT on output) to the n-by-2 matrix [ VECX1 VECX2 ]. The rotation is designed to zero VECX2(1). That is,
[ VECX1 VECX2 ] H = [ (VECX1*H11 + VECX2*H21) (VECX1*H12 + VECX2*H22) ], where VECX1(1)*H12 + VECX2(1)*H22 = 0 and H is the 2-by-2 matrix:
( H11 H12 )
( H21 H22 )
and H takes one of the following forms:
(1 0) ,if TYPE_ROT = 0.
(0 1)
(1 0) (1 A) ,if TYPE_ROT = 1.
(B 1) (0 1)
(1 A) (1 0) ,if TYPE_ROT = 2.
(0 1) (B 1)
(0 -1) (1 0) ,if TYPE_ROT = 3.
(1 A) (-B 1)
(B 1) (1 A) ,if TYPE_ROT = 4.
(1 0) (0 -1)
Furthermore, the scale factors (D1 and D2) are updated accordingly. That is, if on input:
[ Y1 Y2 ] = [ VECX1 VECX2 ] diag( SQRT(D1) SQRT(D2) )then on output:
[ VECX1 VECX2 ] diag( SQRT(D1) SQRT(D2) ) = [ Y1 Y2 ] ROT = [ (Y1*CS + Y2*SN) (-SN*Y1 + CS*Y2) ]where CS and SN define a standard Givens 2-by-2 plane rotation, ROT, which zeros -SN*Y1(1) + CS*Y2(1). Thus, ROT has the following structure:
( +CS -SN )
( +SN +CS )
with CS**(2) + SN**(2) = 1.
Arguments¶
- VECX1 (INPUT/OUTPUT) real(stnd), dimension(:)
The first vector to be transformed.
On output, VECX1*H11 + VECX2*H21 overwrites VECX1.
- VECX2 (INPUT/OUTPUT) real(stnd), dimension(:)
The second vector to be transformed.
On output, VECX1*H12 + VECX2*H22 overwrites VECX2.
- D1 (INPUT/OUTPUT) real(stnd)
On input, first scale factor.
On output, D1 is replaced with the update scale factor.
- D2 (INPUT/OUTPUT) real(stnd)
On input, second scale factor.
On output, D2 is replaced with the update scale factor.
- BETA (OUTPUT) real(stnd)
- The real scalar B which defines the transformation matrix H.
- ALPHA (OUTPUT) real(stnd)
- The real scalar A which defines the transformation matrix H.
- TYPE_ROT (OUTPUT) integer(i2b)
- Integer which defines the transformation matrix H.
Further Details¶
It is assumed that D1 and D2 are positive scalars and that VECX1 and VECX2 have the same size.
IF D1>=D2, D1 is diminished and D2 is augmented. IF D1<D2, D2 is diminished and D1 is augmented. The decrease or increase in magnitude of D1 and D2 are bounded by 1/2 and 2, respectively.
This subroutine is a square root free implementation of the two-way branch algorithm (e.g. fast plane rotations with dynamic scaling to avoid overflow/underflow) described in reference (1).
For further details, see:
- Anda, A.A. and Park, H., 1994:
- Fast plane rotations with dynamic scaling. Siam J. Matrix Anal. Appl., 15, 162-174.
subroutine fastgivens2_vec ( vecx1, vecx2, d1, d2 )¶
Purpose¶
FASTGIVENS2 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,
[ VECX1 VECX2 ] H = [ (VECX1*H11 + VECX2*H21) (VECX1*H12 + VECX2*H22) ], where VECX1(1)*H12 + VECX2(1)*H22 = 0 and H is the 2-by-2 matrix:
( H11 H12 )
( H21 H22 )
, where VECX1(1)*H12 + VECX2(1)*H22 = 0.
Furthermore, the scale factors (D1 and D2) are updated accordingly. That is, if on input:
[ Y1 Y2 ] = [ VECX1 VECX2 ] diag( D1 D2 )then on output:
[ VECX1 VECX2 ] diag( D1 D2 ) = [ Y1 Y2 ] ROT = [ (Y1*CS + Y2*SN) (-SN*Y1 + CS*Y2) ]where CS and SN define a standard Givens 2-by-2 plane rotation, ROT, which zeros -SN*Y1(1) + CS*Y2(1). Thus, ROT has the following structure:
( +CS -SN )
( +SN +CS )
with CS**(2) + SN**(2) = 1.
Arguments¶
- VECX1 (INPUT/OUTPUT) real(stnd), dimension(:)
The first vector to be transformed.
On output, VECX1*H11 + VECX2*H21 overwrites VECX1.
- VECX2 (INPUT/OUTPUT) real(stnd), dimension(:)
The second vector to be transformed.
On output, VECX1*H12 + VECX2*H22 overwrites VECX2.
- D1 (INPUT/OUTPUT) real(stnd)
On input, first scale factor.
On output, D1 is replaced with the update scale factor.
- D2 (INPUT/OUTPUT) real(stnd)
On input, second scale factor.
On output, D2 is replaced with the update scale factor.
Further Details¶
It is assumed that D1 and D2 are positive scalars and that VECX1 and VECX2 have the same size.
IF D1>=D2, D1 is diminished and D2 is augmented. IF D1<D2, D2 is diminished and D1 is augmented. The decrease or increase in magnitude of D1 and D2 are bounded by 1/sqrt(2) and sqrt(2), respectively.
This subroutine is an implementation of the two-way branch algorithm (e.g. fast plane rotations with dynamic scaling to avoid overflow/underflow) described in reference (1).
For further details, see:
- Anda, A.A. and Park, H., 1994:
- Fast plane rotations with dynamic scaling. Siam J. Matrix Anal. Appl., 15, 162-174.
subroutine fastgivens2_vec ( vecx1, vecx2, d1, d2, beta, alpha, type_rot )¶
Purpose¶
FASTGIVENS2 generates and applies a fast Givens plane rotation H (defined by BETA, ALPHA and TYPE_ROT on output) to the n-by-2 matrix [ VECX1 VECX2 ]. The rotation is designed to zero VECX2(1). That is,
[ VECX1 VECX2 ] H = [ (VECX1*H11 + VECX2*H21) (VECX1*H12 + VECX2*H22) ], where VECX1(1)*H12 + VECX2(1)*H22 = 0 and H is the 2-by-2 matrix:
( H11 H12 )
( H21 H22 )
and takes one of the following forms:
(1 0) ,if TYPE_ROT = 0.
(0 1)
(1 0) (1 A) ,if TYPE_ROT = 1.
(B 1) (0 1)
(1 A) (1 0) ,if TYPE_ROT = 2.
(0 1) (B 1)
(0 -1) (1 0) ,if TYPE_ROT = 3.
(1 A) (-B 1)
(B 1) (1 A) ,if TYPE_ROT = 4.
(1 0) (0 -1)
Furthermore, the scale factors (D1 and D2) are updated accordingly. That is, if on input:
[ Y1 Y2 ] = [ VECX1 VECX2 ] diag( D1 D2 )then on output:
[ VECX1 VECX2 ] diag( D1 D2 ) = [ Y1 Y2 ] ROT = [ (Y1*CS + Y2*SN) (-SN*Y1 + CS*Y2) ]where CS and SN define a standard Givens 2-by-2 plane rotation, ROT, which zeros -SN*Y1(1) + CS*Y2(1). Thus, ROT has the following structure:
( +CS -SN )
( +SN +CS )
with CS**(2) + SN**(2) = 1.
Arguments¶
- VECX1 (INPUT/OUTPUT) real(stnd), dimension(:)
The first vector to be transformed.
On output, VECX1*H11 + VECX2*H21 overwrites VECX1.
- VECX2 (INPUT/OUTPUT) real(stnd), dimension(:)
The second vector to be transformed.
On output, VECX1*H12 + VECX2*H22 overwrites VECX2.
- D1 (INPUT/OUTPUT) real(stnd)
On input, first scale factor.
On output, D1 is replaced with the update scale factor.
- D2 (INPUT/OUTPUT) real(stnd)
On input, second scale factor.
On output, D2 is replaced with the update scale factor.
- BETA (OUTPUT) real(stnd)
- The real scalar B which defines the transformation matrix H.
- ALPHA (OUTPUT) real(stnd)
- The real scalar A which defines the transformation matrix H.
- TYPE_ROT (OUTPUT) integer(i2b)
- Integer which defines the transformation matrix H.
Further Details¶
It is assumed that D1 and D2 are positive scalars and that VECX1 and VECX2 have the same size.
IF D1>=D2, D1 is diminished and D2 is augmented. IF D1<D2, D2 is diminished and D1 is augmented. The decrease or increase in magnitude of D1 and D2 are bounded by 1/sqrt(2) and sqrt(2), respectively.
This subroutine is an implementation of the two-way branch algorithm (e.g. fast plane rotations with dynamic scaling to avoid overflow/underflow) described in reference (1).
For further details, see:
- Anda, A.A. and Park, H., 1994:
- Fast plane rotations with dynamic scaling. Siam J. Matrix Anal. Appl., 15, 162-174.
subroutine fastgivens_mat_left ( mat, matd )¶
Purpose¶
FASTGIVENS_MAT_LEFT reduces the matrix MAT to upper trapezoidal form by applying a serie of fast Givens plane rotations on the rows of MAT.
The (row) scale factors (MATD) are updated accordingly.
Arguments¶
- MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)
The matrix to be transformed.
On output, the transformed matrix overwrites MAT.
- MATD (INPUT/OUTPUT) real(stnd), dimension(:)
On input, scale factors associated with the rows of MAT.
On output, MATD is replaced with the update scale factors.
Further Details¶
It is assumed that size( MATD ) = size( MAT, 1 ) and that MATD is a positive vector.
See description of FASTGIVENS for further details.
subroutine fastgivens_mat_right ( mat, matd )¶
Purpose¶
FASTGIVENS_MAT_RIGHT reduces the matrix MAT to lower trapezoidal form by applying a serie of fast Givens plane rotations on the columns of MAT.
The (column) scale factors (MATD) are updated accordingly.
Arguments¶
- MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)
The matrix to be transformed.
On output, the transformed matrix overwrites MAT.
- MATD (INPUT/OUTPUT) real(stnd), dimension(:)
On input, scale factors associated with the columns of MAT.
On output, MATD is replaced with the update scale factors.
Further Details¶
It is assumed that size( MATD ) = size( MAT, 2 ) and that MATD is a positive vector.
See description of FASTGIVENS for further details.
subroutine fastgivens_vec_mat_left ( vec, mat, vecd, matd )¶
Purpose¶
FASTGIVENS_VEC_MAT_LEFT defines and applies a serie 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 annilhate all the elements of the first column of MAT.
The (row) scale factors (VECD and MATD) are updated accordingly.
Arguments¶
- VEC (INPUT/OUTPUT) real(stnd), dimension(:)
On input, the n-vector to rotate. VEC(1) is used to define the rotations.
On output, the transformed vector overwrites VEC.
- MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)
On input, the matrix to be transformed. MAT(:,1) is used to define the rotations.
On output, the transformed matrix overwrites MAT and MAT(:,1) is equal to zero (within numerical accuracy).
- VECD (INPUT/OUTPUT) real(stnd)
On input, scale factor associated with the n-vector VEC.
On output, VECD is replaced with the update scale factor.
- MATD (INPUT/OUTPUT) real(stnd), dimension(:)
On input, scale factors associated with the rows of MAT.
On output, MATD is replaced with the update scale factors.
Further Details¶
It is assumed that:
- size( VEC ) = size( MAT, 2 );
- VECD is a positive scalar;
- size(MATD) = size(MAT,1) and that MATD is a positive vector.
See description of FASTGIVENS for further details.
subroutine fastgivens_vec_mat_right ( vec, mat, vecd, matd )¶
Purpose¶
FASTGIVENS_VEC_MAT_RIGHT defines and applies a serie of fast Givens plane rotations on the m-vector VEC and on the columns of a m-by-n matrix MAT. The rotations are designed to annilhate all the elements of the first row of MAT.
The (column) scale factors (VECD and MATD) are updated accordingly.
Arguments¶
- VEC (INPUT/OUTPUT) real(stnd), dimension(:)
On input, the m-vector to rotate. VEC(1) is used to define the rotations.
On output, the transformed vector overwrites VEC.
- MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)
On input, the matrix to be transformed. MAT(1,:) is used to define the rotations.
On output, the transformed matrix overwrites MAT and MAT(1,:) is equal to zero (within numerical accuracy).
- VECD (INPUT/OUTPUT) real(stnd)
On input, scale factor associated with the m-vector VEC.
On output, VECD is replaced with the update scale factor.
- MATD (INPUT/OUTPUT) real(stnd), dimension(:)
On input, scale factors associated with the columns of MAT.
On output, MATD is replaced with the update scale factors.
Further Details¶
It is assumed that:
- size( VEC ) = size( MAT, 1 );
- VECD is a positive scalar;
- size( MATD ) = size (MAT, 2 ) and that MATD is a positive vector.
See description of FASTGIVENS for further details.