LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
slaexc.f
Go to the documentation of this file.
00001 *> \brief \b SLAEXC
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SLAEXC + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaexc.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaexc.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaexc.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SLAEXC( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK,
00022 *                          INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       LOGICAL            WANTQ
00026 *       INTEGER            INFO, J1, LDQ, LDT, N, N1, N2
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       REAL               Q( LDQ, * ), T( LDT, * ), WORK( * )
00030 *       ..
00031 *  
00032 *
00033 *> \par Purpose:
00034 *  =============
00035 *>
00036 *> \verbatim
00037 *>
00038 *> SLAEXC swaps adjacent diagonal blocks T11 and T22 of order 1 or 2 in
00039 *> an upper quasi-triangular matrix T by an orthogonal similarity
00040 *> transformation.
00041 *>
00042 *> T must be in Schur canonical form, that is, block upper triangular
00043 *> with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block
00044 *> has its diagonal elemnts equal and its off-diagonal elements of
00045 *> opposite sign.
00046 *> \endverbatim
00047 *
00048 *  Arguments:
00049 *  ==========
00050 *
00051 *> \param[in] WANTQ
00052 *> \verbatim
00053 *>          WANTQ is LOGICAL
00054 *>          = .TRUE. : accumulate the transformation in the matrix Q;
00055 *>          = .FALSE.: do not accumulate the transformation.
00056 *> \endverbatim
00057 *>
00058 *> \param[in] N
00059 *> \verbatim
00060 *>          N is INTEGER
00061 *>          The order of the matrix T. N >= 0.
00062 *> \endverbatim
00063 *>
00064 *> \param[in,out] T
00065 *> \verbatim
00066 *>          T is REAL array, dimension (LDT,N)
00067 *>          On entry, the upper quasi-triangular matrix T, in Schur
00068 *>          canonical form.
00069 *>          On exit, the updated matrix T, again in Schur canonical form.
00070 *> \endverbatim
00071 *>
00072 *> \param[in] LDT
00073 *> \verbatim
00074 *>          LDT is INTEGER
00075 *>          The leading dimension of the array T. LDT >= max(1,N).
00076 *> \endverbatim
00077 *>
00078 *> \param[in,out] Q
00079 *> \verbatim
00080 *>          Q is REAL array, dimension (LDQ,N)
00081 *>          On entry, if WANTQ is .TRUE., the orthogonal matrix Q.
00082 *>          On exit, if WANTQ is .TRUE., the updated matrix Q.
00083 *>          If WANTQ is .FALSE., Q is not referenced.
00084 *> \endverbatim
00085 *>
00086 *> \param[in] LDQ
00087 *> \verbatim
00088 *>          LDQ is INTEGER
00089 *>          The leading dimension of the array Q.
00090 *>          LDQ >= 1; and if WANTQ is .TRUE., LDQ >= N.
00091 *> \endverbatim
00092 *>
00093 *> \param[in] J1
00094 *> \verbatim
00095 *>          J1 is INTEGER
00096 *>          The index of the first row of the first block T11.
00097 *> \endverbatim
00098 *>
00099 *> \param[in] N1
00100 *> \verbatim
00101 *>          N1 is INTEGER
00102 *>          The order of the first block T11. N1 = 0, 1 or 2.
00103 *> \endverbatim
00104 *>
00105 *> \param[in] N2
00106 *> \verbatim
00107 *>          N2 is INTEGER
00108 *>          The order of the second block T22. N2 = 0, 1 or 2.
00109 *> \endverbatim
00110 *>
00111 *> \param[out] WORK
00112 *> \verbatim
00113 *>          WORK is REAL array, dimension (N)
00114 *> \endverbatim
00115 *>
00116 *> \param[out] INFO
00117 *> \verbatim
00118 *>          INFO is INTEGER
00119 *>          = 0: successful exit
00120 *>          = 1: the transformed matrix T would be too far from Schur
00121 *>               form; the blocks are not swapped and T and Q are
00122 *>               unchanged.
00123 *> \endverbatim
00124 *
00125 *  Authors:
00126 *  ========
00127 *
00128 *> \author Univ. of Tennessee 
00129 *> \author Univ. of California Berkeley 
00130 *> \author Univ. of Colorado Denver 
00131 *> \author NAG Ltd. 
00132 *
00133 *> \date November 2011
00134 *
00135 *> \ingroup realOTHERauxiliary
00136 *
00137 *  =====================================================================
00138       SUBROUTINE SLAEXC( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK,
00139      $                   INFO )
00140 *
00141 *  -- LAPACK auxiliary routine (version 3.4.0) --
00142 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00143 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00144 *     November 2011
00145 *
00146 *     .. Scalar Arguments ..
00147       LOGICAL            WANTQ
00148       INTEGER            INFO, J1, LDQ, LDT, N, N1, N2
00149 *     ..
00150 *     .. Array Arguments ..
00151       REAL               Q( LDQ, * ), T( LDT, * ), WORK( * )
00152 *     ..
00153 *
00154 *  =====================================================================
00155 *
00156 *     .. Parameters ..
00157       REAL               ZERO, ONE
00158       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00159       REAL               TEN
00160       PARAMETER          ( TEN = 1.0E+1 )
00161       INTEGER            LDD, LDX
00162       PARAMETER          ( LDD = 4, LDX = 2 )
00163 *     ..
00164 *     .. Local Scalars ..
00165       INTEGER            IERR, J2, J3, J4, K, ND
00166       REAL               CS, DNORM, EPS, SCALE, SMLNUM, SN, T11, T22,
00167      $                   T33, TAU, TAU1, TAU2, TEMP, THRESH, WI1, WI2,
00168      $                   WR1, WR2, XNORM
00169 *     ..
00170 *     .. Local Arrays ..
00171       REAL               D( LDD, 4 ), U( 3 ), U1( 3 ), U2( 3 ),
00172      $                   X( LDX, 2 )
00173 *     ..
00174 *     .. External Functions ..
00175       REAL               SLAMCH, SLANGE
00176       EXTERNAL           SLAMCH, SLANGE
00177 *     ..
00178 *     .. External Subroutines ..
00179       EXTERNAL           SLACPY, SLANV2, SLARFG, SLARFX, SLARTG, SLASY2,
00180      $                   SROT
00181 *     ..
00182 *     .. Intrinsic Functions ..
00183       INTRINSIC          ABS, MAX
00184 *     ..
00185 *     .. Executable Statements ..
00186 *
00187       INFO = 0
00188 *
00189 *     Quick return if possible
00190 *
00191       IF( N.EQ.0 .OR. N1.EQ.0 .OR. N2.EQ.0 )
00192      $   RETURN
00193       IF( J1+N1.GT.N )
00194      $   RETURN
00195 *
00196       J2 = J1 + 1
00197       J3 = J1 + 2
00198       J4 = J1 + 3
00199 *
00200       IF( N1.EQ.1 .AND. N2.EQ.1 ) THEN
00201 *
00202 *        Swap two 1-by-1 blocks.
00203 *
00204          T11 = T( J1, J1 )
00205          T22 = T( J2, J2 )
00206 *
00207 *        Determine the transformation to perform the interchange.
00208 *
00209          CALL SLARTG( T( J1, J2 ), T22-T11, CS, SN, TEMP )
00210 *
00211 *        Apply transformation to the matrix T.
00212 *
00213          IF( J3.LE.N )
00214      $      CALL SROT( N-J1-1, T( J1, J3 ), LDT, T( J2, J3 ), LDT, CS,
00215      $                 SN )
00216          CALL SROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN )
00217 *
00218          T( J1, J1 ) = T22
00219          T( J2, J2 ) = T11
00220 *
00221          IF( WANTQ ) THEN
00222 *
00223 *           Accumulate transformation in the matrix Q.
00224 *
00225             CALL SROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN )
00226          END IF
00227 *
00228       ELSE
00229 *
00230 *        Swapping involves at least one 2-by-2 block.
00231 *
00232 *        Copy the diagonal block of order N1+N2 to the local array D
00233 *        and compute its norm.
00234 *
00235          ND = N1 + N2
00236          CALL SLACPY( 'Full', ND, ND, T( J1, J1 ), LDT, D, LDD )
00237          DNORM = SLANGE( 'Max', ND, ND, D, LDD, WORK )
00238 *
00239 *        Compute machine-dependent threshold for test for accepting
00240 *        swap.
00241 *
00242          EPS = SLAMCH( 'P' )
00243          SMLNUM = SLAMCH( 'S' ) / EPS
00244          THRESH = MAX( TEN*EPS*DNORM, SMLNUM )
00245 *
00246 *        Solve T11*X - X*T22 = scale*T12 for X.
00247 *
00248          CALL SLASY2( .FALSE., .FALSE., -1, N1, N2, D, LDD,
00249      $                D( N1+1, N1+1 ), LDD, D( 1, N1+1 ), LDD, SCALE, X,
00250      $                LDX, XNORM, IERR )
00251 *
00252 *        Swap the adjacent diagonal blocks.
00253 *
00254          K = N1 + N1 + N2 - 3
00255          GO TO ( 10, 20, 30 )K
00256 *
00257    10    CONTINUE
00258 *
00259 *        N1 = 1, N2 = 2: generate elementary reflector H so that:
00260 *
00261 *        ( scale, X11, X12 ) H = ( 0, 0, * )
00262 *
00263          U( 1 ) = SCALE
00264          U( 2 ) = X( 1, 1 )
00265          U( 3 ) = X( 1, 2 )
00266          CALL SLARFG( 3, U( 3 ), U, 1, TAU )
00267          U( 3 ) = ONE
00268          T11 = T( J1, J1 )
00269 *
00270 *        Perform swap provisionally on diagonal block in D.
00271 *
00272          CALL SLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK )
00273          CALL SLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK )
00274 *
00275 *        Test whether to reject swap.
00276 *
00277          IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 3,
00278      $       3 )-T11 ) ).GT.THRESH )GO TO 50
00279 *
00280 *        Accept swap: apply transformation to the entire matrix T.
00281 *
00282          CALL SLARFX( 'L', 3, N-J1+1, U, TAU, T( J1, J1 ), LDT, WORK )
00283          CALL SLARFX( 'R', J2, 3, U, TAU, T( 1, J1 ), LDT, WORK )
00284 *
00285          T( J3, J1 ) = ZERO
00286          T( J3, J2 ) = ZERO
00287          T( J3, J3 ) = T11
00288 *
00289          IF( WANTQ ) THEN
00290 *
00291 *           Accumulate transformation in the matrix Q.
00292 *
00293             CALL SLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK )
00294          END IF
00295          GO TO 40
00296 *
00297    20    CONTINUE
00298 *
00299 *        N1 = 2, N2 = 1: generate elementary reflector H so that:
00300 *
00301 *        H (  -X11 ) = ( * )
00302 *          (  -X21 ) = ( 0 )
00303 *          ( scale ) = ( 0 )
00304 *
00305          U( 1 ) = -X( 1, 1 )
00306          U( 2 ) = -X( 2, 1 )
00307          U( 3 ) = SCALE
00308          CALL SLARFG( 3, U( 1 ), U( 2 ), 1, TAU )
00309          U( 1 ) = ONE
00310          T33 = T( J3, J3 )
00311 *
00312 *        Perform swap provisionally on diagonal block in D.
00313 *
00314          CALL SLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK )
00315          CALL SLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK )
00316 *
00317 *        Test whether to reject swap.
00318 *
00319          IF( MAX( ABS( D( 2, 1 ) ), ABS( D( 3, 1 ) ), ABS( D( 1,
00320      $       1 )-T33 ) ).GT.THRESH )GO TO 50
00321 *
00322 *        Accept swap: apply transformation to the entire matrix T.
00323 *
00324          CALL SLARFX( 'R', J3, 3, U, TAU, T( 1, J1 ), LDT, WORK )
00325          CALL SLARFX( 'L', 3, N-J1, U, TAU, T( J1, J2 ), LDT, WORK )
00326 *
00327          T( J1, J1 ) = T33
00328          T( J2, J1 ) = ZERO
00329          T( J3, J1 ) = ZERO
00330 *
00331          IF( WANTQ ) THEN
00332 *
00333 *           Accumulate transformation in the matrix Q.
00334 *
00335             CALL SLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK )
00336          END IF
00337          GO TO 40
00338 *
00339    30    CONTINUE
00340 *
00341 *        N1 = 2, N2 = 2: generate elementary reflectors H(1) and H(2) so
00342 *        that:
00343 *
00344 *        H(2) H(1) (  -X11  -X12 ) = (  *  * )
00345 *                  (  -X21  -X22 )   (  0  * )
00346 *                  ( scale    0  )   (  0  0 )
00347 *                  (    0  scale )   (  0  0 )
00348 *
00349          U1( 1 ) = -X( 1, 1 )
00350          U1( 2 ) = -X( 2, 1 )
00351          U1( 3 ) = SCALE
00352          CALL SLARFG( 3, U1( 1 ), U1( 2 ), 1, TAU1 )
00353          U1( 1 ) = ONE
00354 *
00355          TEMP = -TAU1*( X( 1, 2 )+U1( 2 )*X( 2, 2 ) )
00356          U2( 1 ) = -TEMP*U1( 2 ) - X( 2, 2 )
00357          U2( 2 ) = -TEMP*U1( 3 )
00358          U2( 3 ) = SCALE
00359          CALL SLARFG( 3, U2( 1 ), U2( 2 ), 1, TAU2 )
00360          U2( 1 ) = ONE
00361 *
00362 *        Perform swap provisionally on diagonal block in D.
00363 *
00364          CALL SLARFX( 'L', 3, 4, U1, TAU1, D, LDD, WORK )
00365          CALL SLARFX( 'R', 4, 3, U1, TAU1, D, LDD, WORK )
00366          CALL SLARFX( 'L', 3, 4, U2, TAU2, D( 2, 1 ), LDD, WORK )
00367          CALL SLARFX( 'R', 4, 3, U2, TAU2, D( 1, 2 ), LDD, WORK )
00368 *
00369 *        Test whether to reject swap.
00370 *
00371          IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 4, 1 ) ),
00372      $       ABS( D( 4, 2 ) ) ).GT.THRESH )GO TO 50
00373 *
00374 *        Accept swap: apply transformation to the entire matrix T.
00375 *
00376          CALL SLARFX( 'L', 3, N-J1+1, U1, TAU1, T( J1, J1 ), LDT, WORK )
00377          CALL SLARFX( 'R', J4, 3, U1, TAU1, T( 1, J1 ), LDT, WORK )
00378          CALL SLARFX( 'L', 3, N-J1+1, U2, TAU2, T( J2, J1 ), LDT, WORK )
00379          CALL SLARFX( 'R', J4, 3, U2, TAU2, T( 1, J2 ), LDT, WORK )
00380 *
00381          T( J3, J1 ) = ZERO
00382          T( J3, J2 ) = ZERO
00383          T( J4, J1 ) = ZERO
00384          T( J4, J2 ) = ZERO
00385 *
00386          IF( WANTQ ) THEN
00387 *
00388 *           Accumulate transformation in the matrix Q.
00389 *
00390             CALL SLARFX( 'R', N, 3, U1, TAU1, Q( 1, J1 ), LDQ, WORK )
00391             CALL SLARFX( 'R', N, 3, U2, TAU2, Q( 1, J2 ), LDQ, WORK )
00392          END IF
00393 *
00394    40    CONTINUE
00395 *
00396          IF( N2.EQ.2 ) THEN
00397 *
00398 *           Standardize new 2-by-2 block T11
00399 *
00400             CALL SLANV2( T( J1, J1 ), T( J1, J2 ), T( J2, J1 ),
00401      $                   T( J2, J2 ), WR1, WI1, WR2, WI2, CS, SN )
00402             CALL SROT( N-J1-1, T( J1, J1+2 ), LDT, T( J2, J1+2 ), LDT,
00403      $                 CS, SN )
00404             CALL SROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN )
00405             IF( WANTQ )
00406      $         CALL SROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN )
00407          END IF
00408 *
00409          IF( N1.EQ.2 ) THEN
00410 *
00411 *           Standardize new 2-by-2 block T22
00412 *
00413             J3 = J1 + N2
00414             J4 = J3 + 1
00415             CALL SLANV2( T( J3, J3 ), T( J3, J4 ), T( J4, J3 ),
00416      $                   T( J4, J4 ), WR1, WI1, WR2, WI2, CS, SN )
00417             IF( J3+2.LE.N )
00418      $         CALL SROT( N-J3-1, T( J3, J3+2 ), LDT, T( J4, J3+2 ),
00419      $                    LDT, CS, SN )
00420             CALL SROT( J3-1, T( 1, J3 ), 1, T( 1, J4 ), 1, CS, SN )
00421             IF( WANTQ )
00422      $         CALL SROT( N, Q( 1, J3 ), 1, Q( 1, J4 ), 1, CS, SN )
00423          END IF
00424 *
00425       END IF
00426       RETURN
00427 *
00428 *     Exit with INFO = 1 if swap was rejected.
00429 *
00430    50 INFO = 1
00431       RETURN
00432 *
00433 *     End of SLAEXC
00434 *
00435       END
 All Files Functions