dlaexc.f

Go to the documentation of this file.

*> \brief \b DLAEXC swaps adjacent diagonal blocks of a real upper quasi-triangular matrix in Schur canonical form, by an orthogonal similarity transformation.
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DLAEXC + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaexc.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaexc.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaexc.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE DLAEXC( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK,
*                          INFO )
*
*       .. Scalar Arguments ..
*       LOGICAL            WANTQ
*       INTEGER            INFO, J1, LDQ, LDT, N, N1, N2
*       ..
*       .. Array Arguments ..
*       DOUBLE PRECISION   Q( LDQ, * ), T( LDT, * ), WORK( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DLAEXC swaps adjacent diagonal blocks T11 and T22 of order 1 or 2 in
*> an upper quasi-triangular matrix T by an orthogonal similarity
*> transformation.
*>
*> 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 elements equal and its off-diagonal elements of
*> opposite sign.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] WANTQ
*> \verbatim
*>          WANTQ is LOGICAL
*>          = .TRUE. : accumulate the transformation in the matrix Q;
*>          = .FALSE.: do not accumulate the transformation.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrix T. N >= 0.
*> \endverbatim
*>
*> \param[in,out] T
*> \verbatim
*>          T is 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.
*> \endverbatim
*>
*> \param[in] LDT
*> \verbatim
*>          LDT is INTEGER
*>          The leading dimension of the array T. LDT >= max(1,N).
*> \endverbatim
*>
*> \param[in,out] Q
*> \verbatim
*>          Q is DOUBLE PRECISION array, dimension (LDQ,N)
*>          On entry, if WANTQ is .TRUE., the orthogonal matrix Q.
*>          On exit, if WANTQ is .TRUE., the updated matrix Q.
*>          If WANTQ is .FALSE., Q is not referenced.
*> \endverbatim
*>
*> \param[in] LDQ
*> \verbatim
*>          LDQ is INTEGER
*>          The leading dimension of the array Q.
*>          LDQ >= 1; and if WANTQ is .TRUE., LDQ >= N.
*> \endverbatim
*>
*> \param[in] J1
*> \verbatim
*>          J1 is INTEGER
*>          The index of the first row of the first block T11.
*> \endverbatim
*>
*> \param[in] N1
*> \verbatim
*>          N1 is INTEGER
*>          The order of the first block T11. N1 = 0, 1 or 2.
*> \endverbatim
*>
*> \param[in] N2
*> \verbatim
*>          N2 is INTEGER
*>          The order of the second block T22. N2 = 0, 1 or 2.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is DOUBLE PRECISION array, dimension (N)
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is 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.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup laexc
*
*  =====================================================================
      SUBROUTINE dlaexc( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK,
     $                   INFO )
*
*  -- LAPACK auxiliary routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      LOGICAL            WANTQ
      INTEGER            INFO, J1, LDQ, LDT, N, N1, N2
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   Q( LDQ, * ), T( LDT, * ), WORK( * )
*     ..
*
*  =====================================================================
*
*     .. 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, 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 ), U( 3 ), U1( 3 ), U2( 3 ),
     $                   x( ldx, 2 )
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMCH, DLANGE
      EXTERNAL           dlamch, dlange
*     ..
*     .. External Subroutines ..
      EXTERNAL           dlacpy, 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
*
         IF( wantq ) THEN
*
*           Accumulate transformation in the matrix Q.
*
            CALL drot( n, q( 1, j1 ), 1, q( 1, j2 ), 1, cs, sn )
         END IF
*
      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 dlacpy( '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, * )
*
         u( 1 ) = scale
         u( 2 ) = x( 1, 1 )
         u( 3 ) = x( 1, 2 )
         CALL dlarfg( 3, u( 3 ), u, 1, tau )
         u( 3 ) = one
         t11 = t( j1, j1 )
*
*        Perform swap provisionally on diagonal block in D.
*
         CALL dlarfx( 'L', 3, 3, u, tau, d, ldd, work )
         CALL dlarfx( 'R', 3, 3, u, 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( 'L', 3, n-j1+1, u, tau, t( j1, j1 ), ldt, work )
         CALL dlarfx( 'R', j2, 3, u, tau, t( 1, j1 ), ldt, work )
*
         t( j3, j1 ) = zero
         t( j3, j2 ) = zero
         t( j3, j3 ) = t11
*
         IF( wantq ) THEN
*
*           Accumulate transformation in the matrix Q.
*
            CALL dlarfx( 'R', n, 3, u, tau, q( 1, j1 ), ldq, work )
         END IF
         GO TO 40
*
   20    CONTINUE
*
*        N1 = 2, N2 = 1: generate elementary reflector H so that:
*
*        H (  -X11 ) = ( * )
*          (  -X21 ) = ( 0 )
*          ( scale ) = ( 0 )
*
         u( 1 ) = -x( 1, 1 )
         u( 2 ) = -x( 2, 1 )
         u( 3 ) = scale
         CALL dlarfg( 3, u( 1 ), u( 2 ), 1, tau )
         u( 1 ) = one
         t33 = t( j3, j3 )
*
*        Perform swap provisionally on diagonal block in D.
*
         CALL dlarfx( 'L', 3, 3, u, tau, d, ldd, work )
         CALL dlarfx( 'R', 3, 3, u, 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( 'R', j3, 3, u, tau, t( 1, j1 ), ldt, work )
         CALL dlarfx( 'L', 3, n-j1, u, tau, t( j1, j2 ), ldt, work )
*
         t( j1, j1 ) = t33
         t( j2, j1 ) = zero
         t( j3, j1 ) = zero
*
         IF( wantq ) THEN
*
*           Accumulate transformation in the matrix Q.
*
            CALL dlarfx( 'R', n, 3, u, tau, q( 1, j1 ), ldq, work )
         END IF
         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 )
*
         u1( 1 ) = -x( 1, 1 )
         u1( 2 ) = -x( 2, 1 )
         u1( 3 ) = scale
         CALL dlarfg( 3, u1( 1 ), u1( 2 ), 1, tau1 )
         u1( 1 ) = one
*
         temp = -tau1*( x( 1, 2 )+u1( 2 )*x( 2, 2 ) )
         u2( 1 ) = -temp*u1( 2 ) - x( 2, 2 )
         u2( 2 ) = -temp*u1( 3 )
         u2( 3 ) = scale
         CALL dlarfg( 3, u2( 1 ), u2( 2 ), 1, tau2 )
         u2( 1 ) = one
*
*        Perform swap provisionally on diagonal block in D.
*
         CALL dlarfx( 'L', 3, 4, u1, tau1, d, ldd, work )
         CALL dlarfx( 'R', 4, 3, u1, tau1, d, ldd, work )
         CALL dlarfx( 'L', 3, 4, u2, tau2, d( 2, 1 ), ldd, work )
         CALL dlarfx( 'R', 4, 3, u2, 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 ) ),
     $       abs( d( 4, 2 ) ) ).GT.thresh )GO TO 50
*
*        Accept swap: apply transformation to the entire matrix T.
*
         CALL dlarfx( 'L', 3, n-j1+1, u1, tau1, t( j1, j1 ), ldt, work )
         CALL dlarfx( 'R', j4, 3, u1, tau1, t( 1, j1 ), ldt, work )
         CALL dlarfx( 'L', 3, n-j1+1, u2, tau2, t( j2, j1 ), ldt, work )
         CALL dlarfx( 'R', j4, 3, u2, tau2, t( 1, j2 ), ldt, work )
*
         t( j3, j1 ) = zero
         t( j3, j2 ) = zero
         t( j4, j1 ) = zero
         t( j4, j2 ) = zero
*
         IF( wantq ) THEN
*
*           Accumulate transformation in the matrix Q.
*
            CALL dlarfx( 'R', n, 3, u1, tau1, q( 1, j1 ), ldq, work )
            CALL dlarfx( 'R', n, 3, u2, tau2, q( 1, j2 ), ldq, work )
         END IF
*
   40    CONTINUE
*
         IF( n2.EQ.2 ) 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 )
            IF( wantq )
     $         CALL drot( n, q( 1, j1 ), 1, q( 1, j2 ), 1, cs, sn )
         END IF
*
         IF( n1.EQ.2 ) 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 )
            IF( j3+2.LE.n )
     $         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 )
            IF( wantq )
     $         CALL drot( n, q( 1, j3 ), 1, q( 1, j4 ), 1, cs, sn )
         END IF
*
      END IF
      RETURN
*
*     Exit with INFO = 1 if swap was rejected.
*
   50 CONTINUE
      info = 1
      RETURN
*
*     End of DLAEXC
*
      END