Purpose
To reduce a Hamiltonian like matrix [ A G ] T T H = [ T ] , G = G , Q = Q, [ Q -A ] or a skew-Hamiltonian like matrix [ A G ] T T W = [ T ] , G = -G , Q = -Q, [ Q A ] so that elements below the (k+1)-th subdiagonal in the first nb columns of the (k+n)-by-n matrix A, and offdiagonal elements in the first nb columns and rows of the n-by-n matrix Q are zero. The reduction is performed by an orthogonal symplectic transformation UU'*H*UU and matrices U, XA, XG, XQ, and YA are returned so that [ Aout + U*XA'+ YA*U' Gout + U*XG'+ XG*U' ] UU'*H*UU = [ ]. [ Qout + U*XQ'+ XQ*U' -Aout'- XA*U'- U*YA' ] Similarly, [ Aout + U*XA'+ YA*U' Gout + U*XG'- XG*U' ] UU'*W*UU = [ ]. [ Qout + U*XQ'- XQ*U' Aout'+ XA*U'+ U*YA' ] This is an auxiliary routine called by MB04PB.Specification
SUBROUTINE MB04PA( LHAM, N, K, NB, A, LDA, QG, LDQG, XA, LDXA, $ XG, LDXG, XQ, LDXQ, YA, LDYA, CS, TAU, DWORK ) C .. Scalar Arguments .. LOGICAL LHAM INTEGER K, LDA, LDQG, LDXA, LDXG, LDXQ, LDYA, N, NB C .. Array Arguments .. DOUBLE PRECISION A(LDA,*), CS(*), DWORK(*), QG(LDQG,*), TAU(*), $ XA(LDXA,*), XG(LDXG,*), XQ(LDXQ,*), YA(LDYA,*)Arguments
Mode Parameters
LHAM LOGICAL Specifies the type of matrix to be reduced: = .FALSE. : skew-Hamiltonian like W; = .TRUE. : Hamiltonian like H.Input/Output Parameters
N (input) INTEGER The number of columns of the matrix A. N >= 0. K (input) INTEGER The offset of the reduction. Elements below the (K+1)-th subdiagonal in the first NB columns of A are reduced to zero. K >= 0. NB (input) INTEGER The number of columns/rows to be reduced. N > NB >= 0. A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, the leading (K+N)-by-N part of this array must contain the matrix A. On exit, the leading (K+N)-by-N part of this array contains the matrix Aout and in the zero part information about the elementary reflectors used to compute the reduction. LDA INTEGER The leading dimension of the array A. LDA >= MAX(1,K+N). QG (input/output) DOUBLE PRECISION array, dimension (LDQG,N+1) On entry, the leading N+K-by-N+1 part of this array must contain in the bottom left part the lower triangular part of the N-by-N matrix Q and in the remainder the upper trapezoidal part of the last N columns of the N+K-by-N+K matrix G. On exit, the leading N+K-by-N+1 part of this array contains parts of the matrices Q and G in the same fashion as on entry only that the zero parts of Q contain information about the elementary reflectors used to compute the reduction. Note that if LHAM = .FALSE. then the (K-1)-th and K-th subdiagonals are not referenced. LDQG INTEGER The leading dimension of the array QG. LDQG >= MAX(1,N+K). XA (output) DOUBLE PRECISION array, dimension (LDXA,2*NB) On exit, the leading N-by-(2*NB) part of this array contains the matrix XA. LDXA INTEGER The leading dimension of the array XA. LDXA >= MAX(1,N). XG (output) DOUBLE PRECISION array, dimension (LDXG,2*NB) On exit, the leading (K+N)-by-(2*NB) part of this array contains the matrix XG. LDXG INTEGER The leading dimension of the array XG. LDXG >= MAX(1,K+N). XQ (output) DOUBLE PRECISION array, dimension (LDXQ,2*NB) On exit, the leading N-by-(2*NB) part of this array contains the matrix XQ. LDXQ INTEGER The leading dimension of the array XQ. LDXQ >= MAX(1,N). YA (output) DOUBLE PRECISION array, dimension (LDYA,2*NB) On exit, the leading (K+N)-by-(2*NB) part of this array contains the matrix YA. LDYA INTEGER The leading dimension of the array YA. LDYA >= MAX(1,K+N). CS (output) DOUBLE PRECISION array, dimension (2*NB) On exit, the first 2*NB elements of this array contain the cosines and sines of the symplectic Givens rotations used to compute the reduction. TAU (output) DOUBLE PRECISION array, dimension (NB) On exit, the first NB elements of this array contain the scalar factors of some of the elementary reflectors.Workspace
DWORK DOUBLE PRECISION array, dimension (3*NB)Method
For details regarding the representation of the orthogonal symplectic matrix UU within the arrays A, QG, CS, TAU see the description of MB04PU. The contents of A and QG on exit are illustrated by the following example with n = 5, k = 2 and nb = 2: ( a r r a a ) ( g g r r g g ) ( a r r a a ) ( g g r r g g ) ( a r r a a ) ( q g r r g g ) A = ( r r r r r ), QG = ( t r r r r r ), ( u2 r r r r ) ( u1 t r r r r ) ( u2 u2 r a a ) ( u1 u1 r q g g ) ( u2 u2 r a a ) ( u1 u1 r q q g ) where a, g and q denote elements of the original matrices, r denotes a modified element, t denotes a scalar factor of an applied elementary reflector and ui denote elements of the matrix U.References
[1] C. F. VAN LOAN: A symplectic method for approximating all the eigenvalues of a Hamiltonian matrix. Linear Algebra and its Applications, 61, pp. 233-251, 1984. [2] D. KRESSNER: Block algorithms for orthogonal symplectic factorizations. BIT, 43 (4), pp. 775-790, 2003.Further Comments
NoneExample
Program Text
NoneProgram Data
NoneProgram Results
None