Purpose
To form the triangular block factors R, S and T of a symplectic block reflector SH, which is defined as a product of 2k concatenated Householder reflectors and k Givens rotations, SH = diag( H(1),H(1) ) G(1) diag( F(1),F(1) ) diag( H(2),H(2) ) G(2) diag( F(2),F(2) ) .... diag( H(k),H(k) ) G(k) diag( F(k),F(k) ). The upper triangular blocks of the matrices [ S1 ] [ T11 T12 T13 ] R = [ R1 R2 R3 ], S = [ S2 ], T = [ T21 T22 T23 ], [ S3 ] [ T31 T32 T33 ] with R2 unit and S1, R3, T21, T31, T32 strictly upper triangular, are stored rowwise in the arrays RS and T, respectively.Specification
SUBROUTINE MB04QF( DIRECT, STOREV, STOREW, N, K, V, LDV, W, LDW, $ CS, TAU, RS, LDRS, T, LDT, DWORK ) C .. Scalar Arguments .. CHARACTER DIRECT, STOREV, STOREW INTEGER K, LDRS, LDT, LDV, LDW, N C .. Array Arguments .. DOUBLE PRECISION CS(*), DWORK(*), RS(LDRS,*), T(LDT,*), $ TAU(*), V(LDV,*), W(LDW,*)Arguments
Mode Parameters
DIRECT CHARACTER*1 This is a dummy argument, which is reserved for future extensions of this subroutine. Not referenced. STOREV CHARACTER*1 Specifies how the vectors which define the concatenated Householder F(i) reflectors are stored: = 'C': columnwise; = 'R': rowwise. STOREW CHARACTER*1 Specifies how the vectors which define the concatenated Householder H(i) reflectors are stored: = 'C': columnwise; = 'R': rowwise.Input/Output Parameters
N (input) INTEGER The order of the Householder reflectors F(i) and H(i). N >= 0. K (input) INTEGER The number of Givens rotations. K >= 1. V (input) DOUBLE PRECISION array, dimension (LDV,K) if STOREV = 'C', (LDV,N) if STOREV = 'R' On entry with STOREV = 'C', the leading N-by-K part of this array must contain in its i-th column the vector which defines the elementary reflector F(i). On entry with STOREV = 'R', the leading K-by-N part of this array must contain in its i-th row the vector which defines the elementary reflector F(i). LDV INTEGER The leading dimension of the array V. LDV >= MAX(1,N), if STOREV = 'C'; LDV >= K, if STOREV = 'R'. W (input) DOUBLE PRECISION array, dimension (LDW,K) if STOREW = 'C', (LDW,N) if STOREW = 'R' On entry with STOREW = 'C', the leading N-by-K part of this array must contain in its i-th column the vector which defines the elementary reflector H(i). On entry with STOREV = 'R', the leading K-by-N part of this array must contain in its i-th row the vector which defines the elementary reflector H(i). LDW INTEGER The leading dimension of the array W. LDW >= MAX(1,N), if STOREW = 'C'; LDW >= K, if STOREW = 'R'. CS (input) DOUBLE PRECISION array, dimension (2*K) On entry, the first 2*K elements of this array must contain the cosines and sines of the symplectic Givens rotations G(i). TAU (input) DOUBLE PRECISION array, dimension (K) On entry, the first K elements of this array must contain the scalar factors of the elementary reflectors F(i). RS (output) DOUBLE PRECISION array, dimension (K,6*K) On exit, the leading K-by-6*K part of this array contains the upper triangular matrices defining the factors R and S of the symplectic block reflector SH. The (strictly) lower portions of this array are not used. LDRS INTEGER The leading dimension of the array RS. LDRS >= K. T (output) DOUBLE PRECISION array, dimension (K,9*K) On exit, the leading K-by-9*K part of this array contains the upper triangular matrices defining the factor T of the symplectic block reflector SH. The (strictly) lower portions of this array are not used. LDT INTEGER The leading dimension of the array T. LDT >= K.Workspace
DWORK DOUBLE PRECISION array, dimension (3*K)References
[1] Kressner, D. Block algorithms for orthogonal symplectic factorizations. BIT, 43 (4), pp. 775-790, 2003.Numerical Aspects
The algorithm requires ( 4*K - 2 )*K*N + 19/3*K*K*K + 1/2*K*K + 43/6*K - 4 floating point operations.Further Comments
NoneExample
Program Text
NoneProgram Data
NoneProgram Results
None