LAPACK: SRC/dlaexc.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DLAEXC( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK,
00002      $                   INFO )
00003 *
00004 *  -- LAPACK auxiliary routine (version 3.2.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     June 2010
00008 *
00009 *     .. Scalar Arguments ..
00010       LOGICAL            WANTQ
00011       INTEGER            INFO, J1, LDQ, LDT, N, N1, N2
00012 *     ..
00013 *     .. Array Arguments ..
00014       DOUBLE PRECISION   Q( LDQ, * ), T( LDT, * ), WORK( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  DLAEXC swaps adjacent diagonal blocks T11 and T22 of order 1 or 2 in
00021 *  an upper quasi-triangular matrix T by an orthogonal similarity
00022 *  transformation.
00023 *
00024 *  T must be in Schur canonical form, that is, block upper triangular
00025 *  with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block
00026 *  has its diagonal elemnts equal and its off-diagonal elements of
00027 *  opposite sign.
00028 *
00029 *  Arguments
00030 *  =========
00031 *
00032 *  WANTQ   (input) LOGICAL
00033 *          = .TRUE. : accumulate the transformation in the matrix Q;
00034 *          = .FALSE.: do not accumulate the transformation.
00035 *
00036 *  N       (input) INTEGER
00037 *          The order of the matrix T. N >= 0.
00038 *
00039 *  T       (input/output) DOUBLE PRECISION array, dimension (LDT,N)
00040 *          On entry, the upper quasi-triangular matrix T, in Schur
00041 *          canonical form.
00042 *          On exit, the updated matrix T, again in Schur canonical form.
00043 *
00044 *  LDT     (input) INTEGER
00045 *          The leading dimension of the array T. LDT >= max(1,N).
00046 *
00047 *  Q       (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
00048 *          On entry, if WANTQ is .TRUE., the orthogonal matrix Q.
00049 *          On exit, if WANTQ is .TRUE., the updated matrix Q.
00050 *          If WANTQ is .FALSE., Q is not referenced.
00051 *
00052 *  LDQ     (input) INTEGER
00053 *          The leading dimension of the array Q.
00054 *          LDQ >= 1; and if WANTQ is .TRUE., LDQ >= N.
00055 *
00056 *  J1      (input) INTEGER
00057 *          The index of the first row of the first block T11.
00058 *
00059 *  N1      (input) INTEGER
00060 *          The order of the first block T11. N1 = 0, 1 or 2.
00061 *
00062 *  N2      (input) INTEGER
00063 *          The order of the second block T22. N2 = 0, 1 or 2.
00064 *
00065 *  WORK    (workspace) DOUBLE PRECISION array, dimension (N)
00066 *
00067 *  INFO    (output) INTEGER
00068 *          = 0: successful exit
00069 *          = 1: the transformed matrix T would be too far from Schur
00070 *               form; the blocks are not swapped and T and Q are
00071 *               unchanged.
00072 *
00073 *  =====================================================================
00074 *
00075 *     .. Parameters ..
00076       DOUBLE PRECISION   ZERO, ONE
00077       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00078       DOUBLE PRECISION   TEN
00079       PARAMETER          ( TEN = 1.0D+1 )
00080       INTEGER            LDD, LDX
00081       PARAMETER          ( LDD = 4, LDX = 2 )
00082 *     ..
00083 *     .. Local Scalars ..
00084       INTEGER            IERR, J2, J3, J4, K, ND
00085       DOUBLE PRECISION   CS, DNORM, EPS, SCALE, SMLNUM, SN, T11, T22,
00086      $                   T33, TAU, TAU1, TAU2, TEMP, THRESH, WI1, WI2,
00087      $                   WR1, WR2, XNORM
00088 *     ..
00089 *     .. Local Arrays ..
00090       DOUBLE PRECISION   D( LDD, 4 ), U( 3 ), U1( 3 ), U2( 3 ),
00091      $                   X( LDX, 2 )
00092 *     ..
00093 *     .. External Functions ..
00094       DOUBLE PRECISION   DLAMCH, DLANGE
00095       EXTERNAL           DLAMCH, DLANGE
00096 *     ..
00097 *     .. External Subroutines ..
00098       EXTERNAL           DLACPY, DLANV2, DLARFG, DLARFX, DLARTG, DLASY2,
00099      $                   DROT
00100 *     ..
00101 *     .. Intrinsic Functions ..
00102       INTRINSIC          ABS, MAX
00103 *     ..
00104 *     .. Executable Statements ..
00105 *
00106       INFO = 0
00107 *
00108 *     Quick return if possible
00109 *
00110       IF( N.EQ.0 .OR. N1.EQ.0 .OR. N2.EQ.0 )
00111      $   RETURN
00112       IF( J1+N1.GT.N )
00113      $   RETURN
00114 *
00115       J2 = J1 + 1
00116       J3 = J1 + 2
00117       J4 = J1 + 3
00118 *
00119       IF( N1.EQ.1 .AND. N2.EQ.1 ) THEN
00120 *
00121 *        Swap two 1-by-1 blocks.
00122 *
00123          T11 = T( J1, J1 )
00124          T22 = T( J2, J2 )
00125 *
00126 *        Determine the transformation to perform the interchange.
00127 *
00128          CALL DLARTG( T( J1, J2 ), T22-T11, CS, SN, TEMP )
00129 *
00130 *        Apply transformation to the matrix T.
00131 *
00132          IF( J3.LE.N )
00133      $      CALL DROT( N-J1-1, T( J1, J3 ), LDT, T( J2, J3 ), LDT, CS,
00134      $                 SN )
00135          CALL DROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN )
00136 *
00137          T( J1, J1 ) = T22
00138          T( J2, J2 ) = T11
00139 *
00140          IF( WANTQ ) THEN
00141 *
00142 *           Accumulate transformation in the matrix Q.
00143 *
00144             CALL DROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN )
00145          END IF
00146 *
00147       ELSE
00148 *
00149 *        Swapping involves at least one 2-by-2 block.
00150 *
00151 *        Copy the diagonal block of order N1+N2 to the local array D
00152 *        and compute its norm.
00153 *
00154          ND = N1 + N2
00155          CALL DLACPY( 'Full', ND, ND, T( J1, J1 ), LDT, D, LDD )
00156          DNORM = DLANGE( 'Max', ND, ND, D, LDD, WORK )
00157 *
00158 *        Compute machine-dependent threshold for test for accepting
00159 *        swap.
00160 *
00161          EPS = DLAMCH( 'P' )
00162          SMLNUM = DLAMCH( 'S' ) / EPS
00163          THRESH = MAX( TEN*EPS*DNORM, SMLNUM )
00164 *
00165 *        Solve T11*X - X*T22 = scale*T12 for X.
00166 *
00167          CALL DLASY2( .FALSE., .FALSE., -1, N1, N2, D, LDD,
00168      $                D( N1+1, N1+1 ), LDD, D( 1, N1+1 ), LDD, SCALE, X,
00169      $                LDX, XNORM, IERR )
00170 *
00171 *        Swap the adjacent diagonal blocks.
00172 *
00173          K = N1 + N1 + N2 - 3
00174          GO TO ( 10, 20, 30 )K
00175 *
00176    10    CONTINUE
00177 *
00178 *        N1 = 1, N2 = 2: generate elementary reflector H so that:
00179 *
00180 *        ( scale, X11, X12 ) H = ( 0, 0, * )
00181 *
00182          U( 1 ) = SCALE
00183          U( 2 ) = X( 1, 1 )
00184          U( 3 ) = X( 1, 2 )
00185          CALL DLARFG( 3, U( 3 ), U, 1, TAU )
00186          U( 3 ) = ONE
00187          T11 = T( J1, J1 )
00188 *
00189 *        Perform swap provisionally on diagonal block in D.
00190 *
00191          CALL DLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK )
00192          CALL DLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK )
00193 *
00194 *        Test whether to reject swap.
00195 *
00196          IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 3,
00197      $       3 )-T11 ) ).GT.THRESH )GO TO 50
00198 *
00199 *        Accept swap: apply transformation to the entire matrix T.
00200 *
00201          CALL DLARFX( 'L', 3, N-J1+1, U, TAU, T( J1, J1 ), LDT, WORK )
00202          CALL DLARFX( 'R', J2, 3, U, TAU, T( 1, J1 ), LDT, WORK )
00203 *
00204          T( J3, J1 ) = ZERO
00205          T( J3, J2 ) = ZERO
00206          T( J3, J3 ) = T11
00207 *
00208          IF( WANTQ ) THEN
00209 *
00210 *           Accumulate transformation in the matrix Q.
00211 *
00212             CALL DLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK )
00213          END IF
00214          GO TO 40
00215 *
00216    20    CONTINUE
00217 *
00218 *        N1 = 2, N2 = 1: generate elementary reflector H so that:
00219 *
00220 *        H (  -X11 ) = ( * )
00221 *          (  -X21 ) = ( 0 )
00222 *          ( scale ) = ( 0 )
00223 *
00224          U( 1 ) = -X( 1, 1 )
00225          U( 2 ) = -X( 2, 1 )
00226          U( 3 ) = SCALE
00227          CALL DLARFG( 3, U( 1 ), U( 2 ), 1, TAU )
00228          U( 1 ) = ONE
00229          T33 = T( J3, J3 )
00230 *
00231 *        Perform swap provisionally on diagonal block in D.
00232 *
00233          CALL DLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK )
00234          CALL DLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK )
00235 *
00236 *        Test whether to reject swap.
00237 *
00238          IF( MAX( ABS( D( 2, 1 ) ), ABS( D( 3, 1 ) ), ABS( D( 1,
00239      $       1 )-T33 ) ).GT.THRESH )GO TO 50
00240 *
00241 *        Accept swap: apply transformation to the entire matrix T.
00242 *
00243          CALL DLARFX( 'R', J3, 3, U, TAU, T( 1, J1 ), LDT, WORK )
00244          CALL DLARFX( 'L', 3, N-J1, U, TAU, T( J1, J2 ), LDT, WORK )
00245 *
00246          T( J1, J1 ) = T33
00247          T( J2, J1 ) = ZERO
00248          T( J3, J1 ) = ZERO
00249 *
00250          IF( WANTQ ) THEN
00251 *
00252 *           Accumulate transformation in the matrix Q.
00253 *
00254             CALL DLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK )
00255          END IF
00256          GO TO 40
00257 *
00258    30    CONTINUE
00259 *
00260 *        N1 = 2, N2 = 2: generate elementary reflectors H(1) and H(2) so
00261 *        that:
00262 *
00263 *        H(2) H(1) (  -X11  -X12 ) = (  *  * )
00264 *                  (  -X21  -X22 )   (  0  * )
00265 *                  ( scale    0  )   (  0  0 )
00266 *                  (    0  scale )   (  0  0 )
00267 *
00268          U1( 1 ) = -X( 1, 1 )
00269          U1( 2 ) = -X( 2, 1 )
00270          U1( 3 ) = SCALE
00271          CALL DLARFG( 3, U1( 1 ), U1( 2 ), 1, TAU1 )
00272          U1( 1 ) = ONE
00273 *
00274          TEMP = -TAU1*( X( 1, 2 )+U1( 2 )*X( 2, 2 ) )
00275          U2( 1 ) = -TEMP*U1( 2 ) - X( 2, 2 )
00276          U2( 2 ) = -TEMP*U1( 3 )
00277          U2( 3 ) = SCALE
00278          CALL DLARFG( 3, U2( 1 ), U2( 2 ), 1, TAU2 )
00279          U2( 1 ) = ONE
00280 *
00281 *        Perform swap provisionally on diagonal block in D.
00282 *
00283          CALL DLARFX( 'L', 3, 4, U1, TAU1, D, LDD, WORK )
00284          CALL DLARFX( 'R', 4, 3, U1, TAU1, D, LDD, WORK )
00285          CALL DLARFX( 'L', 3, 4, U2, TAU2, D( 2, 1 ), LDD, WORK )
00286          CALL DLARFX( 'R', 4, 3, U2, TAU2, D( 1, 2 ), LDD, WORK )
00287 *
00288 *        Test whether to reject swap.
00289 *
00290          IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 4, 1 ) ),
00291      $       ABS( D( 4, 2 ) ) ).GT.THRESH )GO TO 50
00292 *
00293 *        Accept swap: apply transformation to the entire matrix T.
00294 *
00295          CALL DLARFX( 'L', 3, N-J1+1, U1, TAU1, T( J1, J1 ), LDT, WORK )
00296          CALL DLARFX( 'R', J4, 3, U1, TAU1, T( 1, J1 ), LDT, WORK )
00297          CALL DLARFX( 'L', 3, N-J1+1, U2, TAU2, T( J2, J1 ), LDT, WORK )
00298          CALL DLARFX( 'R', J4, 3, U2, TAU2, T( 1, J2 ), LDT, WORK )
00299 *
00300          T( J3, J1 ) = ZERO
00301          T( J3, J2 ) = ZERO
00302          T( J4, J1 ) = ZERO
00303          T( J4, J2 ) = ZERO
00304 *
00305          IF( WANTQ ) THEN
00306 *
00307 *           Accumulate transformation in the matrix Q.
00308 *
00309             CALL DLARFX( 'R', N, 3, U1, TAU1, Q( 1, J1 ), LDQ, WORK )
00310             CALL DLARFX( 'R', N, 3, U2, TAU2, Q( 1, J2 ), LDQ, WORK )
00311          END IF
00312 *
00313    40    CONTINUE
00314 *
00315          IF( N2.EQ.2 ) THEN
00316 *
00317 *           Standardize new 2-by-2 block T11
00318 *
00319             CALL DLANV2( T( J1, J1 ), T( J1, J2 ), T( J2, J1 ),
00320      $                   T( J2, J2 ), WR1, WI1, WR2, WI2, CS, SN )
00321             CALL DROT( N-J1-1, T( J1, J1+2 ), LDT, T( J2, J1+2 ), LDT,
00322      $                 CS, SN )
00323             CALL DROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN )
00324             IF( WANTQ )
00325      $         CALL DROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN )
00326          END IF
00327 *
00328          IF( N1.EQ.2 ) THEN
00329 *
00330 *           Standardize new 2-by-2 block T22
00331 *
00332             J3 = J1 + N2
00333             J4 = J3 + 1
00334             CALL DLANV2( T( J3, J3 ), T( J3, J4 ), T( J4, J3 ),
00335      $                   T( J4, J4 ), WR1, WI1, WR2, WI2, CS, SN )
00336             IF( J3+2.LE.N )
00337      $         CALL DROT( N-J3-1, T( J3, J3+2 ), LDT, T( J4, J3+2 ),
00338      $                    LDT, CS, SN )
00339             CALL DROT( J3-1, T( 1, J3 ), 1, T( 1, J4 ), 1, CS, SN )
00340             IF( WANTQ )
00341      $         CALL DROT( N, Q( 1, J3 ), 1, Q( 1, J4 ), 1, CS, SN )
00342          END IF
00343 *
00344       END IF
00345       RETURN
00346 *
00347 *     Exit with INFO = 1 if swap was rejected.
00348 *
00349    50 CONTINUE
00350       INFO = 1
00351       RETURN
00352 *
00353 *     End of DLAEXC
00354 *
00355       END