SRC/dlaexc.f

   1 *> \brief \b DLAEXC swaps adjacent diagonal blocks of a real upper quasi-triangular matrix in Schur canonical form, by an orthogonal similarity transformation.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DLAEXC + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaexc.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaexc.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaexc.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DLAEXC( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK,
  22 *                          INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       LOGICAL            WANTQ
  26 *       INTEGER            INFO, J1, LDQ, LDT, N, N1, N2
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       DOUBLE PRECISION   Q( LDQ, * ), T( LDT, * ), WORK( * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> DLAEXC swaps adjacent diagonal blocks T11 and T22 of order 1 or 2 in
  39 *> an upper quasi-triangular matrix T by an orthogonal similarity
  40 *> transformation.
  41 *>
  42 *> T must be in Schur canonical form, that is, block upper triangular
  43 *> with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block
  44 *> has its diagonal elemnts equal and its off-diagonal elements of
  45 *> opposite sign.
  46 *> \endverbatim
  47 *
  48 *  Arguments:
  49 *  ==========
  50 *
  51 *> \param[in] WANTQ
  52 *> \verbatim
  53 *>          WANTQ is LOGICAL
  54 *>          = .TRUE. : accumulate the transformation in the matrix Q;
  55 *>          = .FALSE.: do not accumulate the transformation.
  56 *> \endverbatim
  57 *>
  58 *> \param[in] N
  59 *> \verbatim
  60 *>          N is INTEGER
  61 *>          The order of the matrix T. N >= 0.
  62 *> \endverbatim
  63 *>
  64 *> \param[in,out] T
  65 *> \verbatim
  66 *>          T is DOUBLE PRECISION array, dimension (LDT,N)
  67 *>          On entry, the upper quasi-triangular matrix T, in Schur
  68 *>          canonical form.
  69 *>          On exit, the updated matrix T, again in Schur canonical form.
  70 *> \endverbatim
  71 *>
  72 *> \param[in] LDT
  73 *> \verbatim
  74 *>          LDT is INTEGER
  75 *>          The leading dimension of the array T. LDT >= max(1,N).
  76 *> \endverbatim
  77 *>
  78 *> \param[in,out] Q
  79 *> \verbatim
  80 *>          Q is DOUBLE PRECISION array, dimension (LDQ,N)
  81 *>          On entry, if WANTQ is .TRUE., the orthogonal matrix Q.
  82 *>          On exit, if WANTQ is .TRUE., the updated matrix Q.
  83 *>          If WANTQ is .FALSE., Q is not referenced.
  84 *> \endverbatim
  85 *>
  86 *> \param[in] LDQ
  87 *> \verbatim
  88 *>          LDQ is INTEGER
  89 *>          The leading dimension of the array Q.
  90 *>          LDQ >= 1; and if WANTQ is .TRUE., LDQ >= N.
  91 *> \endverbatim
  92 *>
  93 *> \param[in] J1
  94 *> \verbatim
  95 *>          J1 is INTEGER
  96 *>          The index of the first row of the first block T11.
  97 *> \endverbatim
  98 *>
  99 *> \param[in] N1
 100 *> \verbatim
 101 *>          N1 is INTEGER
 102 *>          The order of the first block T11. N1 = 0, 1 or 2.
 103 *> \endverbatim
 104 *>
 105 *> \param[in] N2
 106 *> \verbatim
 107 *>          N2 is INTEGER
 108 *>          The order of the second block T22. N2 = 0, 1 or 2.
 109 *> \endverbatim
 110 *>
 111 *> \param[out] WORK
 112 *> \verbatim
 113 *>          WORK is DOUBLE PRECISION array, dimension (N)
 114 *> \endverbatim
 115 *>
 116 *> \param[out] INFO
 117 *> \verbatim
 118 *>          INFO is INTEGER
 119 *>          = 0: successful exit
 120 *>          = 1: the transformed matrix T would be too far from Schur
 121 *>               form; the blocks are not swapped and T and Q are
 122 *>               unchanged.
 123 *> \endverbatim
 124 *
 125 *  Authors:
 126 *  ========
 127 *
 128 *> \author Univ. of Tennessee
 129 *> \author Univ. of California Berkeley
 130 *> \author Univ. of Colorado Denver
 131 *> \author NAG Ltd.
 132 *
 133 *> \date September 2012
 134 *
 135 *> \ingroup doubleOTHERauxiliary
 136 *
 137 *  =====================================================================
 138       SUBROUTINE DLAEXC( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK,
 139      $                   INFO )
 140 *
 141 *  -- LAPACK auxiliary routine (version 3.4.2) --
 142 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 143 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 144 *     September 2012
 145 *
 146 *     .. Scalar Arguments ..
 147       LOGICAL            WANTQ
 148       INTEGER            INFO, J1, LDQ, LDT, N, N1, N2
 149 *     ..
 150 *     .. Array Arguments ..
 151       DOUBLE PRECISION   Q( LDQ, * ), T( LDT, * ), WORK( * )
 152 *     ..
 153 *
 154 *  =====================================================================
 155 *
 156 *     .. Parameters ..
 157       DOUBLE PRECISION   ZERO, ONE
 158       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 159       DOUBLE PRECISION   TEN
 160       PARAMETER          ( TEN = 1.0D+1 )
 161       INTEGER            LDD, LDX
 162       PARAMETER          ( LDD = 4, LDX = 2 )
 163 *     ..
 164 *     .. Local Scalars ..
 165       INTEGER            IERR, J2, J3, J4, K, ND
 166       DOUBLE PRECISION   CS, DNORM, EPS, SCALE, SMLNUM, SN, T11, T22,
 167      $                   T33, TAU, TAU1, TAU2, TEMP, THRESH, WI1, WI2,
 168      $                   WR1, WR2, XNORM
 169 *     ..
 170 *     .. Local Arrays ..
 171       DOUBLE PRECISION   D( LDD, 4 ), U( 3 ), U1( 3 ), U2( 3 ),
 172      $                   X( LDX, 2 )
 173 *     ..
 174 *     .. External Functions ..
 175       DOUBLE PRECISION   DLAMCH, DLANGE
 176       EXTERNAL           DLAMCH, DLANGE
 177 *     ..
 178 *     .. External Subroutines ..
 179       EXTERNAL           DLACPY, DLANV2, DLARFG, DLARFX, DLARTG, DLASY2,
 180      $                   DROT
 181 *     ..
 182 *     .. Intrinsic Functions ..
 183       INTRINSIC          ABS, MAX
 184 *     ..
 185 *     .. Executable Statements ..
 186 *
 187       INFO = 0
 188 *
 189 *     Quick return if possible
 190 *
 191       IF( N.EQ.0 .OR. N1.EQ.0 .OR. N2.EQ.0 )
 192      $   RETURN
 193       IF( J1+N1.GT.N )
 194      $   RETURN
 195 *
 196       J2 = J1 + 1
 197       J3 = J1 + 2
 198       J4 = J1 + 3
 199 *
 200       IF( N1.EQ.1 .AND. N2.EQ.1 ) THEN
 201 *
 202 *        Swap two 1-by-1 blocks.
 203 *
 204          T11 = T( J1, J1 )
 205          T22 = T( J2, J2 )
 206 *
 207 *        Determine the transformation to perform the interchange.
 208 *
 209          CALL DLARTG( T( J1, J2 ), T22-T11, CS, SN, TEMP )
 210 *
 211 *        Apply transformation to the matrix T.
 212 *
 213          IF( J3.LE.N )
 214      $      CALL DROT( N-J1-1, T( J1, J3 ), LDT, T( J2, J3 ), LDT, CS,
 215      $                 SN )
 216          CALL DROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN )
 217 *
 218          T( J1, J1 ) = T22
 219          T( J2, J2 ) = T11
 220 *
 221          IF( WANTQ ) THEN
 222 *
 223 *           Accumulate transformation in the matrix Q.
 224 *
 225             CALL DROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN )
 226          END IF
 227 *
 228       ELSE
 229 *
 230 *        Swapping involves at least one 2-by-2 block.
 231 *
 232 *        Copy the diagonal block of order N1+N2 to the local array D
 233 *        and compute its norm.
 234 *
 235          ND = N1 + N2
 236          CALL DLACPY( 'Full', ND, ND, T( J1, J1 ), LDT, D, LDD )
 237          DNORM = DLANGE( 'Max', ND, ND, D, LDD, WORK )
 238 *
 239 *        Compute machine-dependent threshold for test for accepting
 240 *        swap.
 241 *
 242          EPS = DLAMCH( 'P' )
 243          SMLNUM = DLAMCH( 'S' ) / EPS
 244          THRESH = MAX( TEN*EPS*DNORM, SMLNUM )
 245 *
 246 *        Solve T11*X - X*T22 = scale*T12 for X.
 247 *
 248          CALL DLASY2( .FALSE., .FALSE., -1, N1, N2, D, LDD,
 249      $                D( N1+1, N1+1 ), LDD, D( 1, N1+1 ), LDD, SCALE, X,
 250      $                LDX, XNORM, IERR )
 251 *
 252 *        Swap the adjacent diagonal blocks.
 253 *
 254          K = N1 + N1 + N2 - 3
 255          GO TO ( 10, 20, 30 )K
 256 *
 257    10    CONTINUE
 258 *
 259 *        N1 = 1, N2 = 2: generate elementary reflector H so that:
 260 *
 261 *        ( scale, X11, X12 ) H = ( 0, 0, * )
 262 *
 263          U( 1 ) = SCALE
 264          U( 2 ) = X( 1, 1 )
 265          U( 3 ) = X( 1, 2 )
 266          CALL DLARFG( 3, U( 3 ), U, 1, TAU )
 267          U( 3 ) = ONE
 268          T11 = T( J1, J1 )
 269 *
 270 *        Perform swap provisionally on diagonal block in D.
 271 *
 272          CALL DLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK )
 273          CALL DLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK )
 274 *
 275 *        Test whether to reject swap.
 276 *
 277          IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 3,
 278      $       3 )-T11 ) ).GT.THRESH )GO TO 50
 279 *
 280 *        Accept swap: apply transformation to the entire matrix T.
 281 *
 282          CALL DLARFX( 'L', 3, N-J1+1, U, TAU, T( J1, J1 ), LDT, WORK )
 283          CALL DLARFX( 'R', J2, 3, U, TAU, T( 1, J1 ), LDT, WORK )
 284 *
 285          T( J3, J1 ) = ZERO
 286          T( J3, J2 ) = ZERO
 287          T( J3, J3 ) = T11
 288 *
 289          IF( WANTQ ) THEN
 290 *
 291 *           Accumulate transformation in the matrix Q.
 292 *
 293             CALL DLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK )
 294          END IF
 295          GO TO 40
 296 *
 297    20    CONTINUE
 298 *
 299 *        N1 = 2, N2 = 1: generate elementary reflector H so that:
 300 *
 301 *        H (  -X11 ) = ( * )
 302 *          (  -X21 ) = ( 0 )
 303 *          ( scale ) = ( 0 )
 304 *
 305          U( 1 ) = -X( 1, 1 )
 306          U( 2 ) = -X( 2, 1 )
 307          U( 3 ) = SCALE
 308          CALL DLARFG( 3, U( 1 ), U( 2 ), 1, TAU )
 309          U( 1 ) = ONE
 310          T33 = T( J3, J3 )
 311 *
 312 *        Perform swap provisionally on diagonal block in D.
 313 *
 314          CALL DLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK )
 315          CALL DLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK )
 316 *
 317 *        Test whether to reject swap.
 318 *
 319          IF( MAX( ABS( D( 2, 1 ) ), ABS( D( 3, 1 ) ), ABS( D( 1,
 320      $       1 )-T33 ) ).GT.THRESH )GO TO 50
 321 *
 322 *        Accept swap: apply transformation to the entire matrix T.
 323 *
 324          CALL DLARFX( 'R', J3, 3, U, TAU, T( 1, J1 ), LDT, WORK )
 325          CALL DLARFX( 'L', 3, N-J1, U, TAU, T( J1, J2 ), LDT, WORK )
 326 *
 327          T( J1, J1 ) = T33
 328          T( J2, J1 ) = ZERO
 329          T( J3, J1 ) = ZERO
 330 *
 331          IF( WANTQ ) THEN
 332 *
 333 *           Accumulate transformation in the matrix Q.
 334 *
 335             CALL DLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK )
 336          END IF
 337          GO TO 40
 338 *
 339    30    CONTINUE
 340 *
 341 *        N1 = 2, N2 = 2: generate elementary reflectors H(1) and H(2) so
 342 *        that:
 343 *
 344 *        H(2) H(1) (  -X11  -X12 ) = (  *  * )
 345 *                  (  -X21  -X22 )   (  0  * )
 346 *                  ( scale    0  )   (  0  0 )
 347 *                  (    0  scale )   (  0  0 )
 348 *
 349          U1( 1 ) = -X( 1, 1 )
 350          U1( 2 ) = -X( 2, 1 )
 351          U1( 3 ) = SCALE
 352          CALL DLARFG( 3, U1( 1 ), U1( 2 ), 1, TAU1 )
 353          U1( 1 ) = ONE
 354 *
 355          TEMP = -TAU1*( X( 1, 2 )+U1( 2 )*X( 2, 2 ) )
 356          U2( 1 ) = -TEMP*U1( 2 ) - X( 2, 2 )
 357          U2( 2 ) = -TEMP*U1( 3 )
 358          U2( 3 ) = SCALE
 359          CALL DLARFG( 3, U2( 1 ), U2( 2 ), 1, TAU2 )
 360          U2( 1 ) = ONE
 361 *
 362 *        Perform swap provisionally on diagonal block in D.
 363 *
 364          CALL DLARFX( 'L', 3, 4, U1, TAU1, D, LDD, WORK )
 365          CALL DLARFX( 'R', 4, 3, U1, TAU1, D, LDD, WORK )
 366          CALL DLARFX( 'L', 3, 4, U2, TAU2, D( 2, 1 ), LDD, WORK )
 367          CALL DLARFX( 'R', 4, 3, U2, TAU2, D( 1, 2 ), LDD, WORK )
 368 *
 369 *        Test whether to reject swap.
 370 *
 371          IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 4, 1 ) ),
 372      $       ABS( D( 4, 2 ) ) ).GT.THRESH )GO TO 50
 373 *
 374 *        Accept swap: apply transformation to the entire matrix T.
 375 *
 376          CALL DLARFX( 'L', 3, N-J1+1, U1, TAU1, T( J1, J1 ), LDT, WORK )
 377          CALL DLARFX( 'R', J4, 3, U1, TAU1, T( 1, J1 ), LDT, WORK )
 378          CALL DLARFX( 'L', 3, N-J1+1, U2, TAU2, T( J2, J1 ), LDT, WORK )
 379          CALL DLARFX( 'R', J4, 3, U2, TAU2, T( 1, J2 ), LDT, WORK )
 380 *
 381          T( J3, J1 ) = ZERO
 382          T( J3, J2 ) = ZERO
 383          T( J4, J1 ) = ZERO
 384          T( J4, J2 ) = ZERO
 385 *
 386          IF( WANTQ ) THEN
 387 *
 388 *           Accumulate transformation in the matrix Q.
 389 *
 390             CALL DLARFX( 'R', N, 3, U1, TAU1, Q( 1, J1 ), LDQ, WORK )
 391             CALL DLARFX( 'R', N, 3, U2, TAU2, Q( 1, J2 ), LDQ, WORK )
 392          END IF
 393 *
 394    40    CONTINUE
 395 *
 396          IF( N2.EQ.2 ) THEN
 397 *
 398 *           Standardize new 2-by-2 block T11
 399 *
 400             CALL DLANV2( T( J1, J1 ), T( J1, J2 ), T( J2, J1 ),
 401      $                   T( J2, J2 ), WR1, WI1, WR2, WI2, CS, SN )
 402             CALL DROT( N-J1-1, T( J1, J1+2 ), LDT, T( J2, J1+2 ), LDT,
 403      $                 CS, SN )
 404             CALL DROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN )
 405             IF( WANTQ )
 406      $         CALL DROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN )
 407          END IF
 408 *
 409          IF( N1.EQ.2 ) THEN
 410 *
 411 *           Standardize new 2-by-2 block T22
 412 *
 413             J3 = J1 + N2
 414             J4 = J3 + 1
 415             CALL DLANV2( T( J3, J3 ), T( J3, J4 ), T( J4, J3 ),
 416      $                   T( J4, J4 ), WR1, WI1, WR2, WI2, CS, SN )
 417             IF( J3+2.LE.N )
 418      $         CALL DROT( N-J3-1, T( J3, J3+2 ), LDT, T( J4, J3+2 ),
 419      $                    LDT, CS, SN )
 420             CALL DROT( J3-1, T( 1, J3 ), 1, T( 1, J4 ), 1, CS, SN )
 421             IF( WANTQ )
 422      $         CALL DROT( N, Q( 1, J3 ), 1, Q( 1, J4 ), 1, CS, SN )
 423          END IF
 424 *
 425       END IF
 426       RETURN
 427 *
 428 *     Exit with INFO = 1 if swap was rejected.
 429 *
 430    50 CONTINUE
 431       INFO = 1
 432       RETURN
 433 *
 434 *     End of DLAEXC
 435 *
 436       END