Return to dlaexc.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence

    1: *> \brief \b DLAEXC
    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 November 2011
  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.0) --
  142: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  143: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  144: *     November 2011
  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