Annotation of rpl/lapack/lapack/dlaexc.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DLAEXC( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK,
! 2: $ 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
! 11: INTEGER INFO, J1, LDQ, LDT, N, N1, N2
! 12: * ..
! 13: * .. Array Arguments ..
! 14: DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WORK( * )
! 15: * ..
! 16: *
! 17: * Purpose
! 18: * =======
! 19: *
! 20: * DLAEXC swaps adjacent diagonal blocks T11 and T22 of order 1 or 2 in
! 21: * an upper quasi-triangular matrix T by an orthogonal similarity
! 22: * transformation.
! 23: *
! 24: * T must be in Schur canonical form, that is, block upper triangular
! 25: * with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block
! 26: * has its diagonal elemnts equal and its off-diagonal elements of
! 27: * opposite sign.
! 28: *
! 29: * Arguments
! 30: * =========
! 31: *
! 32: * WANTQ (input) LOGICAL
! 33: * = .TRUE. : accumulate the transformation in the matrix Q;
! 34: * = .FALSE.: do not accumulate the transformation.
! 35: *
! 36: * N (input) INTEGER
! 37: * The order of the matrix T. N >= 0.
! 38: *
! 39: * T (input/output) DOUBLE PRECISION array, dimension (LDT,N)
! 40: * On entry, the upper quasi-triangular matrix T, in Schur
! 41: * canonical form.
! 42: * On exit, the updated matrix T, again in Schur canonical form.
! 43: *
! 44: * LDT (input) INTEGER
! 45: * The leading dimension of the array T. LDT >= max(1,N).
! 46: *
! 47: * Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
! 48: * On entry, if WANTQ is .TRUE., the orthogonal matrix Q.
! 49: * On exit, if WANTQ is .TRUE., the updated matrix Q.
! 50: * If WANTQ is .FALSE., Q is not referenced.
! 51: *
! 52: * LDQ (input) INTEGER
! 53: * The leading dimension of the array Q.
! 54: * LDQ >= 1; and if WANTQ is .TRUE., LDQ >= N.
! 55: *
! 56: * J1 (input) INTEGER
! 57: * The index of the first row of the first block T11.
! 58: *
! 59: * N1 (input) INTEGER
! 60: * The order of the first block T11. N1 = 0, 1 or 2.
! 61: *
! 62: * N2 (input) INTEGER
! 63: * The order of the second block T22. N2 = 0, 1 or 2.
! 64: *
! 65: * WORK (workspace) DOUBLE PRECISION array, dimension (N)
! 66: *
! 67: * INFO (output) INTEGER
! 68: * = 0: successful exit
! 69: * = 1: the transformed matrix T would be too far from Schur
! 70: * form; the blocks are not swapped and T and Q are
! 71: * unchanged.
! 72: *
! 73: * =====================================================================
! 74: *
! 75: * .. Parameters ..
! 76: DOUBLE PRECISION ZERO, ONE
! 77: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
! 78: DOUBLE PRECISION TEN
! 79: PARAMETER ( TEN = 1.0D+1 )
! 80: INTEGER LDD, LDX
! 81: PARAMETER ( LDD = 4, LDX = 2 )
! 82: * ..
! 83: * .. Local Scalars ..
! 84: INTEGER IERR, J2, J3, J4, K, ND
! 85: DOUBLE PRECISION CS, DNORM, EPS, SCALE, SMLNUM, SN, T11, T22,
! 86: $ T33, TAU, TAU1, TAU2, TEMP, THRESH, WI1, WI2,
! 87: $ WR1, WR2, XNORM
! 88: * ..
! 89: * .. Local Arrays ..
! 90: DOUBLE PRECISION D( LDD, 4 ), U( 3 ), U1( 3 ), U2( 3 ),
! 91: $ X( LDX, 2 )
! 92: * ..
! 93: * .. External Functions ..
! 94: DOUBLE PRECISION DLAMCH, DLANGE
! 95: EXTERNAL DLAMCH, DLANGE
! 96: * ..
! 97: * .. External Subroutines ..
! 98: EXTERNAL DLACPY, DLANV2, DLARFG, DLARFX, DLARTG, DLASY2,
! 99: $ DROT
! 100: * ..
! 101: * .. Intrinsic Functions ..
! 102: INTRINSIC ABS, MAX
! 103: * ..
! 104: * .. Executable Statements ..
! 105: *
! 106: INFO = 0
! 107: *
! 108: * Quick return if possible
! 109: *
! 110: IF( N.EQ.0 .OR. N1.EQ.0 .OR. N2.EQ.0 )
! 111: $ RETURN
! 112: IF( J1+N1.GT.N )
! 113: $ RETURN
! 114: *
! 115: J2 = J1 + 1
! 116: J3 = J1 + 2
! 117: J4 = J1 + 3
! 118: *
! 119: IF( N1.EQ.1 .AND. N2.EQ.1 ) THEN
! 120: *
! 121: * Swap two 1-by-1 blocks.
! 122: *
! 123: T11 = T( J1, J1 )
! 124: T22 = T( J2, J2 )
! 125: *
! 126: * Determine the transformation to perform the interchange.
! 127: *
! 128: CALL DLARTG( T( J1, J2 ), T22-T11, CS, SN, TEMP )
! 129: *
! 130: * Apply transformation to the matrix T.
! 131: *
! 132: IF( J3.LE.N )
! 133: $ CALL DROT( N-J1-1, T( J1, J3 ), LDT, T( J2, J3 ), LDT, CS,
! 134: $ SN )
! 135: CALL DROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN )
! 136: *
! 137: T( J1, J1 ) = T22
! 138: T( J2, J2 ) = T11
! 139: *
! 140: IF( WANTQ ) THEN
! 141: *
! 142: * Accumulate transformation in the matrix Q.
! 143: *
! 144: CALL DROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN )
! 145: END IF
! 146: *
! 147: ELSE
! 148: *
! 149: * Swapping involves at least one 2-by-2 block.
! 150: *
! 151: * Copy the diagonal block of order N1+N2 to the local array D
! 152: * and compute its norm.
! 153: *
! 154: ND = N1 + N2
! 155: CALL DLACPY( 'Full', ND, ND, T( J1, J1 ), LDT, D, LDD )
! 156: DNORM = DLANGE( 'Max', ND, ND, D, LDD, WORK )
! 157: *
! 158: * Compute machine-dependent threshold for test for accepting
! 159: * swap.
! 160: *
! 161: EPS = DLAMCH( 'P' )
! 162: SMLNUM = DLAMCH( 'S' ) / EPS
! 163: THRESH = MAX( TEN*EPS*DNORM, SMLNUM )
! 164: *
! 165: * Solve T11*X - X*T22 = scale*T12 for X.
! 166: *
! 167: CALL DLASY2( .FALSE., .FALSE., -1, N1, N2, D, LDD,
! 168: $ D( N1+1, N1+1 ), LDD, D( 1, N1+1 ), LDD, SCALE, X,
! 169: $ LDX, XNORM, IERR )
! 170: *
! 171: * Swap the adjacent diagonal blocks.
! 172: *
! 173: K = N1 + N1 + N2 - 3
! 174: GO TO ( 10, 20, 30 )K
! 175: *
! 176: 10 CONTINUE
! 177: *
! 178: * N1 = 1, N2 = 2: generate elementary reflector H so that:
! 179: *
! 180: * ( scale, X11, X12 ) H = ( 0, 0, * )
! 181: *
! 182: U( 1 ) = SCALE
! 183: U( 2 ) = X( 1, 1 )
! 184: U( 3 ) = X( 1, 2 )
! 185: CALL DLARFG( 3, U( 3 ), U, 1, TAU )
! 186: U( 3 ) = ONE
! 187: T11 = T( J1, J1 )
! 188: *
! 189: * Perform swap provisionally on diagonal block in D.
! 190: *
! 191: CALL DLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK )
! 192: CALL DLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK )
! 193: *
! 194: * Test whether to reject swap.
! 195: *
! 196: IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 3,
! 197: $ 3 )-T11 ) ).GT.THRESH )GO TO 50
! 198: *
! 199: * Accept swap: apply transformation to the entire matrix T.
! 200: *
! 201: CALL DLARFX( 'L', 3, N-J1+1, U, TAU, T( J1, J1 ), LDT, WORK )
! 202: CALL DLARFX( 'R', J2, 3, U, TAU, T( 1, J1 ), LDT, WORK )
! 203: *
! 204: T( J3, J1 ) = ZERO
! 205: T( J3, J2 ) = ZERO
! 206: T( J3, J3 ) = T11
! 207: *
! 208: IF( WANTQ ) THEN
! 209: *
! 210: * Accumulate transformation in the matrix Q.
! 211: *
! 212: CALL DLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK )
! 213: END IF
! 214: GO TO 40
! 215: *
! 216: 20 CONTINUE
! 217: *
! 218: * N1 = 2, N2 = 1: generate elementary reflector H so that:
! 219: *
! 220: * H ( -X11 ) = ( * )
! 221: * ( -X21 ) = ( 0 )
! 222: * ( scale ) = ( 0 )
! 223: *
! 224: U( 1 ) = -X( 1, 1 )
! 225: U( 2 ) = -X( 2, 1 )
! 226: U( 3 ) = SCALE
! 227: CALL DLARFG( 3, U( 1 ), U( 2 ), 1, TAU )
! 228: U( 1 ) = ONE
! 229: T33 = T( J3, J3 )
! 230: *
! 231: * Perform swap provisionally on diagonal block in D.
! 232: *
! 233: CALL DLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK )
! 234: CALL DLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK )
! 235: *
! 236: * Test whether to reject swap.
! 237: *
! 238: IF( MAX( ABS( D( 2, 1 ) ), ABS( D( 3, 1 ) ), ABS( D( 1,
! 239: $ 1 )-T33 ) ).GT.THRESH )GO TO 50
! 240: *
! 241: * Accept swap: apply transformation to the entire matrix T.
! 242: *
! 243: CALL DLARFX( 'R', J3, 3, U, TAU, T( 1, J1 ), LDT, WORK )
! 244: CALL DLARFX( 'L', 3, N-J1, U, TAU, T( J1, J2 ), LDT, WORK )
! 245: *
! 246: T( J1, J1 ) = T33
! 247: T( J2, J1 ) = ZERO
! 248: T( J3, J1 ) = ZERO
! 249: *
! 250: IF( WANTQ ) THEN
! 251: *
! 252: * Accumulate transformation in the matrix Q.
! 253: *
! 254: CALL DLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK )
! 255: END IF
! 256: GO TO 40
! 257: *
! 258: 30 CONTINUE
! 259: *
! 260: * N1 = 2, N2 = 2: generate elementary reflectors H(1) and H(2) so
! 261: * that:
! 262: *
! 263: * H(2) H(1) ( -X11 -X12 ) = ( * * )
! 264: * ( -X21 -X22 ) ( 0 * )
! 265: * ( scale 0 ) ( 0 0 )
! 266: * ( 0 scale ) ( 0 0 )
! 267: *
! 268: U1( 1 ) = -X( 1, 1 )
! 269: U1( 2 ) = -X( 2, 1 )
! 270: U1( 3 ) = SCALE
! 271: CALL DLARFG( 3, U1( 1 ), U1( 2 ), 1, TAU1 )
! 272: U1( 1 ) = ONE
! 273: *
! 274: TEMP = -TAU1*( X( 1, 2 )+U1( 2 )*X( 2, 2 ) )
! 275: U2( 1 ) = -TEMP*U1( 2 ) - X( 2, 2 )
! 276: U2( 2 ) = -TEMP*U1( 3 )
! 277: U2( 3 ) = SCALE
! 278: CALL DLARFG( 3, U2( 1 ), U2( 2 ), 1, TAU2 )
! 279: U2( 1 ) = ONE
! 280: *
! 281: * Perform swap provisionally on diagonal block in D.
! 282: *
! 283: CALL DLARFX( 'L', 3, 4, U1, TAU1, D, LDD, WORK )
! 284: CALL DLARFX( 'R', 4, 3, U1, TAU1, D, LDD, WORK )
! 285: CALL DLARFX( 'L', 3, 4, U2, TAU2, D( 2, 1 ), LDD, WORK )
! 286: CALL DLARFX( 'R', 4, 3, U2, TAU2, D( 1, 2 ), LDD, WORK )
! 287: *
! 288: * Test whether to reject swap.
! 289: *
! 290: IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 4, 1 ) ),
! 291: $ ABS( D( 4, 2 ) ) ).GT.THRESH )GO TO 50
! 292: *
! 293: * Accept swap: apply transformation to the entire matrix T.
! 294: *
! 295: CALL DLARFX( 'L', 3, N-J1+1, U1, TAU1, T( J1, J1 ), LDT, WORK )
! 296: CALL DLARFX( 'R', J4, 3, U1, TAU1, T( 1, J1 ), LDT, WORK )
! 297: CALL DLARFX( 'L', 3, N-J1+1, U2, TAU2, T( J2, J1 ), LDT, WORK )
! 298: CALL DLARFX( 'R', J4, 3, U2, TAU2, T( 1, J2 ), LDT, WORK )
! 299: *
! 300: T( J3, J1 ) = ZERO
! 301: T( J3, J2 ) = ZERO
! 302: T( J4, J1 ) = ZERO
! 303: T( J4, J2 ) = ZERO
! 304: *
! 305: IF( WANTQ ) THEN
! 306: *
! 307: * Accumulate transformation in the matrix Q.
! 308: *
! 309: CALL DLARFX( 'R', N, 3, U1, TAU1, Q( 1, J1 ), LDQ, WORK )
! 310: CALL DLARFX( 'R', N, 3, U2, TAU2, Q( 1, J2 ), LDQ, WORK )
! 311: END IF
! 312: *
! 313: 40 CONTINUE
! 314: *
! 315: IF( N2.EQ.2 ) THEN
! 316: *
! 317: * Standardize new 2-by-2 block T11
! 318: *
! 319: CALL DLANV2( T( J1, J1 ), T( J1, J2 ), T( J2, J1 ),
! 320: $ T( J2, J2 ), WR1, WI1, WR2, WI2, CS, SN )
! 321: CALL DROT( N-J1-1, T( J1, J1+2 ), LDT, T( J2, J1+2 ), LDT,
! 322: $ CS, SN )
! 323: CALL DROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN )
! 324: IF( WANTQ )
! 325: $ CALL DROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN )
! 326: END IF
! 327: *
! 328: IF( N1.EQ.2 ) THEN
! 329: *
! 330: * Standardize new 2-by-2 block T22
! 331: *
! 332: J3 = J1 + N2
! 333: J4 = J3 + 1
! 334: CALL DLANV2( T( J3, J3 ), T( J3, J4 ), T( J4, J3 ),
! 335: $ T( J4, J4 ), WR1, WI1, WR2, WI2, CS, SN )
! 336: IF( J3+2.LE.N )
! 337: $ CALL DROT( N-J3-1, T( J3, J3+2 ), LDT, T( J4, J3+2 ),
! 338: $ LDT, CS, SN )
! 339: CALL DROT( J3-1, T( 1, J3 ), 1, T( 1, J4 ), 1, CS, SN )
! 340: IF( WANTQ )
! 341: $ CALL DROT( N, Q( 1, J3 ), 1, Q( 1, J4 ), 1, CS, SN )
! 342: END IF
! 343: *
! 344: END IF
! 345: RETURN
! 346: *
! 347: * Exit with INFO = 1 if swap was rejected.
! 348: *
! 349: 50 CONTINUE
! 350: INFO = 1
! 351: RETURN
! 352: *
! 353: * End of DLAEXC
! 354: *
! 355: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaexc.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack