#include "ESMF_LapackBlas.inc" !> \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 elemnts 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. ! !> \date December 2016 ! !> \ingroup doubleOTHERauxiliary ! ! ===================================================================== SUBROUTINE DLAEXC( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK, & INFO ) ! ! -- LAPACK auxiliary routine (version 3.7.0) -- ! -- LAPACK is a software package provided by Univ. of Tennessee, -- ! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- ! December 2016 ! ! .. 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