MB04PA

Special reduction of a (skew-)Hamiltonian like matrix

[Specification] [Arguments] [Method] [References] [Comments] [Example]

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
  None
Example

Program Text

  None
Program Data
  None
Program Results
  None

Return to Supporting Routines index