bdlaexc.f

Go to the documentation of this file.

      SUBROUTINE bdlaexc( N, T, LDT, J1, N1, N2, ITRAF, DTRAF, WORK,
     $                    INFO )
      IMPLICIT NONE
*
*  -- ScaLAPACK auxiliary routine (version 2.0.2) --
*     Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver
*     May 1 2012
*
*     .. Scalar Arguments ..
      INTEGER            INFO, J1, LDT, N, N1, N2
*     ..
*     .. Array Arguments ..
      INTEGER            ITRAF( * )
      DOUBLE PRECISION   DTRAF( * ), T( LDT, * ), WORK( * )
*     ..
*
*  Purpose
*  =======
*
*  BDLAEXC swaps adjacent diagonal blocks T11 and T22 of order 1 or 2 in
*  an upper quasi-triangular matrix T by an orthogonal similarity
*  transformation.
*
*  In contrast to the LAPACK routine DLAEXC, the orthogonal
*  transformation matrix Q is not explicitly constructed but
*  represented by paramaters contained in the arrays ITRAF and DTRAF,
*  see the description of BDTREXC for more details.
*
*  T must be in Schur canonical form, that is, block upper triangular
*  with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block
*  has its diagonal elemnts equal and its off-diagonal elements of
*  opposite sign.
*
*  Arguments
*  =========
*
*  N       (input) INTEGER
*          The order of the matrix T. N >= 0.
*
*  T       (input/output) DOUBLE PRECISION array, dimension (LDT,N)
*          On entry, the upper quasi-triangular matrix T, in Schur
*          canonical form.
*          On exit, the updated matrix T, again in Schur canonical form.
*
*  LDT     (input)  INTEGER
*          The leading dimension of the array T. LDT >= max(1,N).
*
*  J1      (input) INTEGER
*          The index of the first row of the first block T11.
*
*  N1      (input) INTEGER
*          The order of the first block T11. N1 = 0, 1 or 2.
*
*  N2      (input) INTEGER
*          The order of the second block T22. N2 = 0, 1 or 2.
*
*  ITRAF   (output) INTEGER array, length k, where
*             k = 1, if N1+N2 = 2;
*             k = 2, if N1+N2 = 3;
*             k = 4, if N1+N2 = 4.
*          List of parameters for representing the transformation
*          matrix Q, see BDTREXC.
*
*  DTRAF   (output) DOUBLE PRECISION array, length k, where
*             k =  2, if N1+N2 = 2;
*             k =  5, if N1+N2 = 3;
*             k = 10, if N1+N2 = 4.
*          List of parameters for representing the transformation
*          matrix Q, see BDTREXC.
*
*  WORK    (workspace) DOUBLE PRECISION array, dimension (N)
*
*  INFO    (output) INTEGER
*          = 0: successful exit
*          = 1: the transformed matrix T would be too far from Schur
*               form; the blocks are not swapped and T and Q are
*               unchanged.
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      parameter( zero = 0.0d+0, one = 1.0d+0 )
      DOUBLE PRECISION   TEN
      parameter( ten = 1.0d+1 )
      INTEGER            LDD, LDX
      parameter( ldd = 4, ldx = 2 )
*     ..
*     .. Local Scalars ..
      INTEGER            IERR, J2, J3, J4, K, LD, LI, ND
      DOUBLE PRECISION   CS, DNORM, EPS, SCALE, SMLNUM, SN, T11, T22,
     $                   t33, tau, tau1, tau2, temp, thresh, wi1, wi2,
     $                   wr1, wr2, xnorm
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION   D( LDD, 4 ), X( LDX, 2 )
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMCH, DLANGE
      EXTERNAL           dlamch, dlange
*     ..
*     .. External Subroutines ..
      EXTERNAL           dlamov, dlanv2, dlarfg, dlarfx, dlartg, dlasy2,
     $                   drot
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max
*     ..
*     .. Executable Statements ..
*
      info = 0
*
*     Quick return if possible
*
      IF( n.EQ.0 .OR. n1.EQ.0 .OR. n2.EQ.0 )
     $   RETURN
      IF( j1+n1.GT.n )
     $   RETURN
*
      j2 = j1 + 1
      j3 = j1 + 2
      j4 = j1 + 3
*
      IF( n1.EQ.1 .AND. n2.EQ.1 ) THEN
*
*        Swap two 1-by-1 blocks.
*
         t11 = t( j1, j1 )
         t22 = t( j2, j2 )
*
*        Determine the transformation to perform the interchange.
*
         CALL dlartg( t( j1, j2 ), t22-t11, cs, sn, temp )
*
*        Apply transformation to the matrix T.
*
         IF( j3.LE.n )
     $      CALL drot( n-j1-1, t( j1, j3 ), ldt, t( j2, j3 ), ldt, cs,
     $                 sn )
         CALL drot( j1-1, t( 1, j1 ), 1, t( 1, j2 ), 1, cs, sn )
*
         t( j1, j1 ) = t22
         t( j2, j2 ) = t11
*
         itraf( 1 ) = j1
         dtraf( 1 ) = cs
         dtraf( 2 ) = sn
*
      ELSE
*
*        Swapping involves at least one 2-by-2 block.
*
*        Copy the diagonal block of order N1+N2 to the local array D
*        and compute its norm.
*
         nd = n1 + n2
         CALL dlamov( 'Full', nd, nd, t( j1, j1 ), ldt, d, ldd )
         dnorm = dlange( 'Max', nd, nd, d, ldd, work )
*
*        Compute machine-dependent threshold for test for accepting
*        swap.
*
         eps = dlamch( 'P' )
         smlnum = dlamch( 'S' ) / eps
         thresh = max( ten*eps*dnorm, smlnum )
*
*        Solve T11*X - X*T22 = scale*T12 for X.
*
         CALL dlasy2( .false., .false., -1, n1, n2, d, ldd,
     $                d( n1+1, n1+1 ), ldd, d( 1, n1+1 ), ldd, scale, x,
     $                ldx, xnorm, ierr )
*
*        Swap the adjacent diagonal blocks.
*
         k = n1 + n1 + n2 - 3
         GO TO ( 10, 20, 30 )k
*
   10    CONTINUE
*
*        N1 = 1, N2 = 2: generate elementary reflector H so that:
*
*        ( scale, X11, X12 ) H = ( 0, 0, * )
*
         dtraf( 1 ) = scale
         dtraf( 2 ) = x( 1, 1 )
         dtraf( 3 ) = x( 1, 2 )
         CALL dlarfg( 3, dtraf( 3 ), dtraf, 1, tau )
         dtraf( 3 ) = one
         t11 = t( j1, j1 )
*
*        Perform swap provisionally on diagonal block in D.
*
         CALL dlarfx( 'Left', 3, 3, dtraf, tau, d, ldd, work )
         CALL dlarfx( 'Right', 3, 3, dtraf, tau, d, ldd, work )
*
*        Test whether to reject swap.
*
         IF( max( abs( d( 3, 1 ) ), abs( d( 3, 2 ) ), abs( d( 3,
     $       3 )-t11 ) ).GT.thresh )GO TO 50
*
*        Accept swap: apply transformation to the entire matrix T.
*
         CALL dlarfx( 'Left', 3, n-j1+1, dtraf, tau, t( j1, j1 ), ldt,
     $                work )
         CALL dlarfx( 'Right', j2, 3, dtraf, tau, t( 1, j1 ), ldt,
     $                work )
*
         t( j3, j1 ) = zero
         t( j3, j2 ) = zero
         t( j3, j3 ) = t11
*
         itraf( 1 ) = 2*n + j1
         li = 2
         dtraf( 3 ) = tau
         ld = 4
         GO TO 40
*
   20    CONTINUE
*
*        N1 = 2, N2 = 1: generate elementary reflector H so that:
*
*        H (  -X11 ) = ( * )
*          (  -X21 ) = ( 0 )
*          ( scale ) = ( 0 )
*
         dtraf( 1 ) = -x( 1, 1 )
         dtraf( 2 ) = -x( 2, 1 )
         dtraf( 3 ) = scale
         CALL dlarfg( 3, dtraf( 1 ), dtraf( 2 ), 1, tau )
         dtraf( 1 ) = one
         t33 = t( j3, j3 )
*
*        Perform swap provisionally on diagonal block in D.
*
         CALL dlarfx( 'Left', 3, 3, dtraf, tau, d, ldd, work )
         CALL dlarfx( 'Right', 3, 3, dtraf, tau, d, ldd, work )
*
*        Test whether to reject swap.
*
         IF( max( abs( d( 2, 1 ) ), abs( d( 3, 1 ) ), abs( d( 1,
     $       1 )-t33 ) ).GT.thresh )GO TO 50
*
*        Accept swap: apply transformation to the entire matrix T.
*
         CALL dlarfx( 'Right', j3, 3, dtraf, tau, t( 1, j1 ), ldt,
     $                work )
         CALL dlarfx( 'Left', 3, n-j1, dtraf, tau, t( j1, j2 ), ldt,
     $                work )
*
         t( j1, j1 ) = t33
         t( j2, j1 ) = zero
         t( j3, j1 ) = zero
*
         itraf( 1 ) = n + j1
         li = 2
         dtraf( 1 ) = tau
         ld = 4
         GO TO 40
*
   30    CONTINUE
*
*        N1 = 2, N2 = 2: generate elementary reflectors H(1) and H(2) so
*        that:
*
*        H(2) H(1) (  -X11  -X12 ) = (  *  * )
*                  (  -X21  -X22 )   (  0  * )
*                  ( scale    0  )   (  0  0 )
*                  (    0  scale )   (  0  0 )
*
         dtraf( 1 ) = -x( 1, 1 )
         dtraf( 2 ) = -x( 2, 1 )
         dtraf( 3 ) = scale
         CALL dlarfg( 3, dtraf( 1 ), dtraf( 2 ), 1, tau1 )
         dtraf( 1 ) = one
*
         temp = -tau1*( x( 1, 2 )+dtraf( 2 )*x( 2, 2 ) )
         dtraf( 4 ) = -temp*dtraf( 2 ) - x( 2, 2 )
         dtraf( 5 ) = -temp*dtraf( 3 )
         dtraf( 6 ) = scale
         CALL dlarfg( 3, dtraf( 4 ), dtraf( 5 ), 1, tau2 )
         dtraf( 4 ) = one
*
*        Perform swap provisionally on diagonal block in D.
*
         CALL dlarfx( 'Left', 3, 4, dtraf, tau1, d, ldd, work )
         CALL dlarfx( 'Right', 4, 3, dtraf, tau1, d, ldd, work )
         CALL dlarfx( 'Left', 3, 4, dtraf( 4 ), tau2, d( 2, 1 ), ldd,
     $                work )
         CALL dlarfx( 'right', 4, 3, DTRAF( 4 ), TAU2, D( 1, 2 ), LDD,
     $                WORK )
*
*        Test whether to reject swap.
*
         IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 4, 1 ) ),
.GT.     $       ABS( D( 4, 2 ) ) )THRESH )GO TO 50
*
*        Accept swap: apply transformation to the entire matrix T.
*
         CALL DLARFX( 'left', 3, N-J1+1, DTRAF, TAU1, T( J1, J1 ), LDT,
     $                WORK )
         CALL DLARFX( 'right', J4, 3, DTRAF, TAU1, T( 1, J1 ), LDT,
     $                WORK )
         CALL DLARFX( 'left', 3, N-J1+1, DTRAF( 4 ), TAU2, T( J2, J1 ),
     $                LDT, WORK )
         CALL DLARFX( 'right', J4, 3, DTRAF( 4 ), TAU2, T( 1, J2 ), LDT,
     $                WORK )
*
         T( J3, J1 ) = ZERO
         T( J3, J2 ) = ZERO
         T( J4, J1 ) = ZERO
         T( J4, J2 ) = ZERO
*
         ITRAF( 1 ) = N + J1
         ITRAF( 2 ) = N + J2
         LI = 3
         DTRAF( 1 ) = TAU1
         DTRAF( 4 ) = TAU2
         LD = 7
         GO TO 40
*
   40    CONTINUE
*
.EQ.         IF( N22 ) THEN
*
*           Standardize new 2-by-2 block T11
*
            CALL DLANV2( T( J1, J1 ), T( J1, J2 ), T( J2, J1 ),
     $                   T( J2, J2 ), WR1, WI1, WR2, WI2, CS, SN )
            CALL DROT( N-J1-1, T( J1, J1+2 ), LDT, T( J2, J1+2 ), LDT,
     $                 CS, SN )
            CALL DROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN )
            ITRAF( LI ) = J1
            LI = LI + 1
            DTRAF( LD ) = CS
            DTRAF( LD+1 ) = SN
            LD = LD + 2
         END IF
*
.EQ.         IF( N12 ) THEN
*
*           Standardize new 2-by-2 block T22
*
            J3 = J1 + N2
            J4 = J3 + 1
            CALL DLANV2( T( J3, J3 ), T( J3, J4 ), T( J4, J3 ),
     $                   T( J4, J4 ), WR1, WI1, WR2, WI2, CS, SN )
.LE.            IF( J3+2N )
     $         CALL DROT( N-J3-1, T( J3, J3+2 ), LDT, T( J4, J3+2 ),
     $                    LDT, CS, SN )
            CALL DROT( J3-1, T( 1, J3 ), 1, T( 1, J4 ), 1, CS, SN )
            ITRAF( LI ) = J3
            DTRAF( LD ) = CS
            DTRAF( LD+1 ) = SN
         END IF
*
      END IF
      RETURN
*
*     Exit with INFO = 1 if swap was rejected.
*
   50 CONTINUE
      INFO = 1
      RETURN
*
*     End of BDLAEXC
*
      END