MB03KB

Swapping pairs of adjacent diagonal blocks of sizes 1 and/or 2 in a formal matrix product

Purpose

  To reorder the diagonal blocks of the formal matrix product

     T22_K^S(K) * T22_K-1^S(K-1) * ... * T22_1^S(1)              (1)

  of length K in the generalized periodic Schur form

           [  T11_k  T12_k  T13_k  ]
     T_k = [    0    T22_k  T23_k  ],    k = 1, ..., K,          (2)
           [    0      0    T33_k  ]

  where

  - the submatrices T11_k are NI(k+1)-by-NI(k), if S(k) = 1, or
    NI(k)-by-NI(k+1), if S(k) = -1, and contain dimension-induced
    infinite eigenvalues,
  - the submatrices T22_k are NC-by-NC and contain core eigenvalues,
    which are generically neither zero nor infinite,
  - the submatrices T33_k contain dimension-induced zero
    eigenvalues,

  such that pairs of adjacent diagonal blocks of sizes 1 and/or 2 in
  the product (1) are swapped.

  Optionally, the transformation matrices Q_1,...,Q_K from the
  reduction into generalized periodic Schur form are updated with
  respect to the performed reordering.

      SUBROUTINE MB03KB( COMPQ, WHICHQ, WS, K, NC, KSCHUR, J1, N1, N2,
     $                   N, NI, S, T, LDT, IXT, Q, LDQ, IXQ, TOL, IWORK,
     $                   DWORK, LDWORK, INFO )
C     .. Scalar Arguments ..
      CHARACTER          COMPQ
      LOGICAL            WS
      INTEGER            INFO, J1, K, KSCHUR, LDWORK, N1, N2, NC
C     .. Array Arguments ..
      INTEGER            IWORK( * ), IXQ( * ), IXT( * ), LDQ( * ),
     $                   LDT( * ), N( * ), NI( * ), S( * ), WHICHQ( * )
      DOUBLE PRECISION   DWORK( * ), Q( * ), T( * ), TOL( * )

Mode Parameters

  COMPQ   CHARACTER*1
          = 'N': do not compute any of the matrices Q_k;
          = 'U': each coefficient of Q must contain an orthogonal
                 matrix Q1_k on entry, and the products Q1_k*Q_k are
                 returned, where Q_k, k = 1, ..., K, performed the
                 reordering;
          = 'W': the computation of each Q_k is specified
                 individually in the array WHICHQ.

  WHICHQ  INTEGER array, dimension (K)
          If COMPQ = 'W', WHICHQ(k) specifies the computation of Q_k
          as follows:
          = 0:   do not compute Q_k;
          > 0:   the kth coefficient of Q must contain an orthogonal
                 matrix Q1_k on entry, and the product Q1_k*Q_k is
                 returned.
          This array is not referenced if COMPQ <> 'W'.

  WS      LOGICAL
          = .FALSE. : do not perform the strong stability tests;
          = .TRUE.  : perform the strong stability tests; often,
                      this is not needed, and omitting them can save
                      some computations.

  K       (input) INTEGER
          The period of the periodic matrix sequences T and Q (the
          number of factors in the matrix product).  K >= 2.
          (For K = 1, a standard eigenvalue reordering problem is
          obtained.)

  NC      (input) INTEGER
          The number of core eigenvalues.  0 <= NC <= min(N).

  KSCHUR  (input) INTEGER
          The index for which the matrix T22_kschur is upper quasi-
          triangular.

  J1      (input) INTEGER
          The index of the first row and column of the first block
          to swap in T22_k.
          1 <= J1 <= NC-N1-N2+1.

  N1      (input) INTEGER
          The order of the first block to swap.   N1 = 0, 1 or 2.

  N2      (input) INTEGER
          The order of the second block to swap.  N2 = 0, 1 or 2.

  N       (input) INTEGER array, dimension (K)
          The leading K elements of this array must contain the
          dimensions of the factors of the formal matrix product T,
          such that the k-th coefficient T_k is an N(k+1)-by-N(k)
          matrix, if S(k) = 1, or an N(k)-by-N(k+1) matrix,
          if S(k) = -1, k = 1, ..., K, where N(K+1) = N(1).

  NI      (input) INTEGER array, dimension (K)
          The leading K elements of this array must contain the
          dimensions of the factors of the matrix sequence T11_k.
          N(k) >= NI(k) + NC >= 0.

  S       (input) INTEGER array, dimension (K)
          The leading K elements of this array must contain the
          signatures (exponents) of the factors in the K-periodic
          matrix sequence. Each entry in S must be either 1 or -1;
          the value S(k) = -1 corresponds to using the inverse of
          the factor T_k.

  T       (input/output) DOUBLE PRECISION array, dimension (*)
          On entry, this array must contain at position IXT(k) the
          matrix T_k, which is at least N(k+1)-by-N(k), if S(k) = 1,
          or at least N(k)-by-N(k+1), if S(k) = -1, in periodic
          Schur form.
          On exit, the matrices T_k are overwritten by the reordered
          periodic Schur form.

  LDT     INTEGER array, dimension (K)
          The leading dimensions of the matrices T_k in the one-
          dimensional array T.
          LDT(k) >= max(1,N(k+1)),  if S(k) =  1,
          LDT(k) >= max(1,N(k)),    if S(k) = -1.

  IXT     INTEGER array, dimension (K)
          Start indices of the matrices T_k in the one-dimensional
          array T.

  Q       (input/output) DOUBLE PRECISION array, dimension (*)
          On entry, this array must contain at position IXQ(k) a
          matrix Q_k of size at least N(k)-by-N(k), provided that
          COMPQ = 'U', or COMPQ = 'W' and WHICHQ(k) > 0.
          On exit, if COMPQ = 'U', or COMPQ = 'W' and WHICHQ(k) > 0,
          Q_k is post-multiplied with the orthogonal matrix that
          performed the reordering.
          This array is not referenced if COMPQ = 'N'.

  LDQ     INTEGER array, dimension (K)
          The leading dimensions of the matrices Q_k in the one-
          dimensional array Q.  LDQ(k) >= 1, and
          LDQ(k) >= max(1,N(k)), if COMPQ = 'U', or COMPQ = 'W' and
                                                    WHICHQ(k) > 0;
          This array is not referenced if COMPQ = 'N'.

  IXQ     INTEGER array, dimension (K)
          Start indices of the matrices Q_k in the one-dimensional
          array Q.
          This array is not referenced if COMPQ = 'N'.

  TOL     DOUBLE PRECISION array, dimension (3)
          This array contains tolerance parameters. The weak and
          strong stability tests use a threshold computed by the
          formula  MAX( c*EPS*NRM, SMLNUM ),  where c is a constant,
          NRM is the Frobenius norm of the matrix formed by
          concatenating K pairs of adjacent diagonal blocks of sizes
          1 and/or 2 in the T22_k submatrices from (2), which are
          swapped, and EPS and SMLNUM are the machine precision and
          safe minimum divided by EPS, respectively (see LAPACK
          Library routine DLAMCH). The norm NRM is computed by this
          routine; the other values are stored in the array TOL.
          TOL(1), TOL(2), and TOL(3) contain c, EPS, and SMLNUM,
          respectively. TOL(1) should normally be at least 10.

  IWORK   INTEGER array, dimension (4*K)

  DWORK   DOUBLE PRECISION array, dimension (LDWORK)
          On exit, if INFO = 0, DWORK(1) returns the optimal LDWORK.

  LDWORK  INTEGER
          The dimension of the array DWORK.
          LDWORK >= 10*K + MN,                 if N1 = 1, N2 = 1;
          LDWORK >= 25*K + MN,                 if N1 = 1, N2 = 2;
          LDWORK >= MAX(23*K + MN, 25*K - 12), if N1 = 2, N2 = 1;
          LDWORK >= MAX(42*K + MN, 80*K - 48), if N1 = 2, N2 = 2;
          where MN = MXN, if MXN > 10, and MN = 0, otherwise, with
          MXN = MAX(N(k),k=1,...,K).

          If LDWORK = -1  a workspace query is assumed; the
          routine only calculates the optimal size of the DWORK
          array, returns this value as the first entry of the DWORK
          array, and no error message is issued by XERBLA.

  INFO    INTEGER
          = 0:  successful exit;
          < 0:  if INFO = -22, then LDWORK is too small; appropriate
                value for LDWORK is returned in DWORK(1); the other
                arguments are not tested, for efficiency;
          = 1:  the swap was rejected from stability reasons; the
                blocks are not swapped and T and Q are unchanged.

  The algorithm described in [1] is used. Both weak and strong
  stability tests are performed.

  [1] Granat, R., Kagstrom, B. and Kressner, D.
      Computing periodic deflating subspaces associated with a
      specified set of eigenvalues.
      BIT Numerical Mathematics, vol. 47, 763-791, 2007.

  The implemented method is numerically backward stable.
                               3
  The algorithm requires 0(K NC ) floating point operations.

  None

Program Text

  None

  None

  None