Return to ztgex2.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence

    1:       SUBROUTINE ZTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
    2:      $                   LDZ, J1, INFO )
    3: *
    4: *  -- LAPACK auxiliary routine (version 3.2) --
    5: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    6: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    7: *     November 2006
    8: *
    9: *     .. Scalar Arguments ..
   10:       LOGICAL            WANTQ, WANTZ
   11:       INTEGER            INFO, J1, LDA, LDB, LDQ, LDZ, N
   12: *     ..
   13: *     .. Array Arguments ..
   14:       COMPLEX*16         A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
   15:      $                   Z( LDZ, * )
   16: *     ..
   17: *
   18: *  Purpose
   19: *  =======
   20: *
   21: *  ZTGEX2 swaps adjacent diagonal 1 by 1 blocks (A11,B11) and (A22,B22)
   22: *  in an upper triangular matrix pair (A, B) by an unitary equivalence
   23: *  transformation.
   24: *
   25: *  (A, B) must be in generalized Schur canonical form, that is, A and
   26: *  B are both upper triangular.
   27: *
   28: *  Optionally, the matrices Q and Z of generalized Schur vectors are
   29: *  updated.
   30: *
   31: *         Q(in) * A(in) * Z(in)' = Q(out) * A(out) * Z(out)'
   32: *         Q(in) * B(in) * Z(in)' = Q(out) * B(out) * Z(out)'
   33: *
   34: *
   35: *  Arguments
   36: *  =========
   37: *
   38: *  WANTQ   (input) LOGICAL
   39: *          .TRUE. : update the left transformation matrix Q;
   40: *          .FALSE.: do not update Q.
   41: *
   42: *  WANTZ   (input) LOGICAL
   43: *          .TRUE. : update the right transformation matrix Z;
   44: *          .FALSE.: do not update Z.
   45: *
   46: *  N       (input) INTEGER
   47: *          The order of the matrices A and B. N >= 0.
   48: *
   49: *  A       (input/output) COMPLEX*16 arrays, dimensions (LDA,N)
   50: *          On entry, the matrix A in the pair (A, B).
   51: *          On exit, the updated matrix A.
   52: *
   53: *  LDA     (input)  INTEGER
   54: *          The leading dimension of the array A. LDA >= max(1,N).
   55: *
   56: *  B       (input/output) COMPLEX*16 arrays, dimensions (LDB,N)
   57: *          On entry, the matrix B in the pair (A, B).
   58: *          On exit, the updated matrix B.
   59: *
   60: *  LDB     (input)  INTEGER
   61: *          The leading dimension of the array B. LDB >= max(1,N).
   62: *
   63: *  Q       (input/output) COMPLEX*16 array, dimension (LDZ,N)
   64: *          If WANTQ = .TRUE, on entry, the unitary matrix Q. On exit,
   65: *          the updated matrix Q.
   66: *          Not referenced if WANTQ = .FALSE..
   67: *
   68: *  LDQ     (input) INTEGER
   69: *          The leading dimension of the array Q. LDQ >= 1;
   70: *          If WANTQ = .TRUE., LDQ >= N.
   71: *
   72: *  Z       (input/output) COMPLEX*16 array, dimension (LDZ,N)
   73: *          If WANTZ = .TRUE, on entry, the unitary matrix Z. On exit,
   74: *          the updated matrix Z.
   75: *          Not referenced if WANTZ = .FALSE..
   76: *
   77: *  LDZ     (input) INTEGER
   78: *          The leading dimension of the array Z. LDZ >= 1;
   79: *          If WANTZ = .TRUE., LDZ >= N.
   80: *
   81: *  J1      (input) INTEGER
   82: *          The index to the first block (A11, B11).
   83: *
   84: *  INFO    (output) INTEGER
   85: *           =0:  Successful exit.
   86: *           =1:  The transformed matrix pair (A, B) would be too far
   87: *                from generalized Schur form; the problem is ill-
   88: *                conditioned. 
   89: *
   90: *
   91: *  Further Details
   92: *  ===============
   93: *
   94: *  Based on contributions by
   95: *     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
   96: *     Umea University, S-901 87 Umea, Sweden.
   97: *
   98: *  In the current code both weak and strong stability tests are
   99: *  performed. The user can omit the strong stability test by changing
  100: *  the internal logical parameter WANDS to .FALSE.. See ref. [2] for
  101: *  details.
  102: *
  103: *  [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
  104: *      Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
  105: *      M.S. Moonen et al (eds), Linear Algebra for Large Scale and
  106: *      Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
  107: *
  108: *  [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
  109: *      Eigenvalues of a Regular Matrix Pair (A, B) and Condition
  110: *      Estimation: Theory, Algorithms and Software, Report UMINF-94.04,
  111: *      Department of Computing Science, Umea University, S-901 87 Umea,
  112: *      Sweden, 1994. Also as LAPACK Working Note 87. To appear in
  113: *      Numerical Algorithms, 1996.
  114: *
  115: *  =====================================================================
  116: *
  117: *     .. Parameters ..
  118:       COMPLEX*16         CZERO, CONE
  119:       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
  120:      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
  121:       DOUBLE PRECISION   TEN
  122:       PARAMETER          ( TEN = 10.0D+0 )
  123:       INTEGER            LDST
  124:       PARAMETER          ( LDST = 2 )
  125:       LOGICAL            WANDS
  126:       PARAMETER          ( WANDS = .TRUE. )
  127: *     ..
  128: *     .. Local Scalars ..
  129:       LOGICAL            DTRONG, WEAK
  130:       INTEGER            I, M
  131:       DOUBLE PRECISION   CQ, CZ, EPS, SA, SB, SCALE, SMLNUM, SS, SUM,
  132:      $                   THRESH, WS
  133:       COMPLEX*16         CDUM, F, G, SQ, SZ
  134: *     ..
  135: *     .. Local Arrays ..
  136:       COMPLEX*16         S( LDST, LDST ), T( LDST, LDST ), WORK( 8 )
  137: *     ..
  138: *     .. External Functions ..
  139:       DOUBLE PRECISION   DLAMCH
  140:       EXTERNAL           DLAMCH
  141: *     ..
  142: *     .. External Subroutines ..
  143:       EXTERNAL           ZLACPY, ZLARTG, ZLASSQ, ZROT
  144: *     ..
  145: *     .. Intrinsic Functions ..
  146:       INTRINSIC          ABS, DBLE, DCONJG, MAX, SQRT
  147: *     ..
  148: *     .. Executable Statements ..
  149: *
  150:       INFO = 0
  151: *
  152: *     Quick return if possible
  153: *
  154:       IF( N.LE.1 )
  155:      $   RETURN
  156: *
  157:       M = LDST
  158:       WEAK = .FALSE.
  159:       DTRONG = .FALSE.
  160: *
  161: *     Make a local copy of selected block in (A, B)
  162: *
  163:       CALL ZLACPY( 'Full', M, M, A( J1, J1 ), LDA, S, LDST )
  164:       CALL ZLACPY( 'Full', M, M, B( J1, J1 ), LDB, T, LDST )
  165: *
  166: *     Compute the threshold for testing the acceptance of swapping.
  167: *
  168:       EPS = DLAMCH( 'P' )
  169:       SMLNUM = DLAMCH( 'S' ) / EPS
  170:       SCALE = DBLE( CZERO )
  171:       SUM = DBLE( CONE )
  172:       CALL ZLACPY( 'Full', M, M, S, LDST, WORK, M )
  173:       CALL ZLACPY( 'Full', M, M, T, LDST, WORK( M*M+1 ), M )
  174:       CALL ZLASSQ( 2*M*M, WORK, 1, SCALE, SUM )
  175:       SA = SCALE*SQRT( SUM )
  176:       THRESH = MAX( TEN*EPS*SA, SMLNUM )
  177: *
  178: *     Compute unitary QL and RQ that swap 1-by-1 and 1-by-1 blocks
  179: *     using Givens rotations and perform the swap tentatively.
  180: *
  181:       F = S( 2, 2 )*T( 1, 1 ) - T( 2, 2 )*S( 1, 1 )
  182:       G = S( 2, 2 )*T( 1, 2 ) - T( 2, 2 )*S( 1, 2 )
  183:       SA = ABS( S( 2, 2 ) )
  184:       SB = ABS( T( 2, 2 ) )
  185:       CALL ZLARTG( G, F, CZ, SZ, CDUM )
  186:       SZ = -SZ
  187:       CALL ZROT( 2, S( 1, 1 ), 1, S( 1, 2 ), 1, CZ, DCONJG( SZ ) )
  188:       CALL ZROT( 2, T( 1, 1 ), 1, T( 1, 2 ), 1, CZ, DCONJG( SZ ) )
  189:       IF( SA.GE.SB ) THEN
  190:          CALL ZLARTG( S( 1, 1 ), S( 2, 1 ), CQ, SQ, CDUM )
  191:       ELSE
  192:          CALL ZLARTG( T( 1, 1 ), T( 2, 1 ), CQ, SQ, CDUM )
  193:       END IF
  194:       CALL ZROT( 2, S( 1, 1 ), LDST, S( 2, 1 ), LDST, CQ, SQ )
  195:       CALL ZROT( 2, T( 1, 1 ), LDST, T( 2, 1 ), LDST, CQ, SQ )
  196: *
  197: *     Weak stability test: |S21| + |T21| <= O(EPS F-norm((S, T)))
  198: *
  199:       WS = ABS( S( 2, 1 ) ) + ABS( T( 2, 1 ) )
  200:       WEAK = WS.LE.THRESH
  201:       IF( .NOT.WEAK )
  202:      $   GO TO 20
  203: *
  204:       IF( WANDS ) THEN
  205: *
  206: *        Strong stability test:
  207: *           F-norm((A-QL'*S*QR, B-QL'*T*QR)) <= O(EPS*F-norm((A, B)))
  208: *
  209:          CALL ZLACPY( 'Full', M, M, S, LDST, WORK, M )
  210:          CALL ZLACPY( 'Full', M, M, T, LDST, WORK( M*M+1 ), M )
  211:          CALL ZROT( 2, WORK, 1, WORK( 3 ), 1, CZ, -DCONJG( SZ ) )
  212:          CALL ZROT( 2, WORK( 5 ), 1, WORK( 7 ), 1, CZ, -DCONJG( SZ ) )
  213:          CALL ZROT( 2, WORK, 2, WORK( 2 ), 2, CQ, -SQ )
  214:          CALL ZROT( 2, WORK( 5 ), 2, WORK( 6 ), 2, CQ, -SQ )
  215:          DO 10 I = 1, 2
  216:             WORK( I ) = WORK( I ) - A( J1+I-1, J1 )
  217:             WORK( I+2 ) = WORK( I+2 ) - A( J1+I-1, J1+1 )
  218:             WORK( I+4 ) = WORK( I+4 ) - B( J1+I-1, J1 )
  219:             WORK( I+6 ) = WORK( I+6 ) - B( J1+I-1, J1+1 )
  220:    10    CONTINUE
  221:          SCALE = DBLE( CZERO )
  222:          SUM = DBLE( CONE )
  223:          CALL ZLASSQ( 2*M*M, WORK, 1, SCALE, SUM )
  224:          SS = SCALE*SQRT( SUM )
  225:          DTRONG = SS.LE.THRESH
  226:          IF( .NOT.DTRONG )
  227:      $      GO TO 20
  228:       END IF
  229: *
  230: *     If the swap is accepted ("weakly" and "strongly"), apply the
  231: *     equivalence transformations to the original matrix pair (A,B)
  232: *
  233:       CALL ZROT( J1+1, A( 1, J1 ), 1, A( 1, J1+1 ), 1, CZ,
  234:      $           DCONJG( SZ ) )
  235:       CALL ZROT( J1+1, B( 1, J1 ), 1, B( 1, J1+1 ), 1, CZ,
  236:      $           DCONJG( SZ ) )
  237:       CALL ZROT( N-J1+1, A( J1, J1 ), LDA, A( J1+1, J1 ), LDA, CQ, SQ )
  238:       CALL ZROT( N-J1+1, B( J1, J1 ), LDB, B( J1+1, J1 ), LDB, CQ, SQ )
  239: *
  240: *     Set  N1 by N2 (2,1) blocks to 0
  241: *
  242:       A( J1+1, J1 ) = CZERO
  243:       B( J1+1, J1 ) = CZERO
  244: *
  245: *     Accumulate transformations into Q and Z if requested.
  246: *
  247:       IF( WANTZ )
  248:      $   CALL ZROT( N, Z( 1, J1 ), 1, Z( 1, J1+1 ), 1, CZ,
  249:      $              DCONJG( SZ ) )
  250:       IF( WANTQ )
  251:      $   CALL ZROT( N, Q( 1, J1 ), 1, Q( 1, J1+1 ), 1, CQ,
  252:      $              DCONJG( SQ ) )
  253: *
  254: *     Exit with INFO = 0 if swap was successfully performed.
  255: *
  256:       RETURN
  257: *
  258: *     Exit with INFO = 1 if swap was rejected.
  259: *
  260:    20 CONTINUE
  261:       INFO = 1
  262:       RETURN
  263: *
  264: *     End of ZTGEX2
  265: *
  266:       END