Return to zlaqr2.f CVS log

Up to [local] / rpl / lapack / lapack

Commit initial.

    1:       SUBROUTINE ZLAQR2( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ,
    2:      $                   IHIZ, Z, LDZ, NS, ND, SH, V, LDV, NH, T, LDT,
    3:      $                   NV, WV, LDWV, WORK, LWORK )
    4: *
    5: *  -- LAPACK auxiliary routine (version 3.2.1)                        --
    6: *     Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
    7: *  -- April 2009                                                      --
    8: *
    9: *     .. Scalar Arguments ..
   10:       INTEGER            IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV,
   11:      $                   LDZ, LWORK, N, ND, NH, NS, NV, NW
   12:       LOGICAL            WANTT, WANTZ
   13: *     ..
   14: *     .. Array Arguments ..
   15:       COMPLEX*16         H( LDH, * ), SH( * ), T( LDT, * ), V( LDV, * ),
   16:      $                   WORK( * ), WV( LDWV, * ), Z( LDZ, * )
   17: *     ..
   18: *
   19: *     This subroutine is identical to ZLAQR3 except that it avoids
   20: *     recursion by calling ZLAHQR instead of ZLAQR4.
   21: *
   22: *
   23: *     ******************************************************************
   24: *     Aggressive early deflation:
   25: *
   26: *     This subroutine accepts as input an upper Hessenberg matrix
   27: *     H and performs an unitary similarity transformation
   28: *     designed to detect and deflate fully converged eigenvalues from
   29: *     a trailing principal submatrix.  On output H has been over-
   30: *     written by a new Hessenberg matrix that is a perturbation of
   31: *     an unitary similarity transformation of H.  It is to be
   32: *     hoped that the final version of H has many zero subdiagonal
   33: *     entries.
   34: *
   35: *     ******************************************************************
   36: *     WANTT   (input) LOGICAL
   37: *          If .TRUE., then the Hessenberg matrix H is fully updated
   38: *          so that the triangular Schur factor may be
   39: *          computed (in cooperation with the calling subroutine).
   40: *          If .FALSE., then only enough of H is updated to preserve
   41: *          the eigenvalues.
   42: *
   43: *     WANTZ   (input) LOGICAL
   44: *          If .TRUE., then the unitary matrix Z is updated so
   45: *          so that the unitary Schur factor may be computed
   46: *          (in cooperation with the calling subroutine).
   47: *          If .FALSE., then Z is not referenced.
   48: *
   49: *     N       (input) INTEGER
   50: *          The order of the matrix H and (if WANTZ is .TRUE.) the
   51: *          order of the unitary matrix Z.
   52: *
   53: *     KTOP    (input) INTEGER
   54: *          It is assumed that either KTOP = 1 or H(KTOP,KTOP-1)=0.
   55: *          KBOT and KTOP together determine an isolated block
   56: *          along the diagonal of the Hessenberg matrix.
   57: *
   58: *     KBOT    (input) INTEGER
   59: *          It is assumed without a check that either
   60: *          KBOT = N or H(KBOT+1,KBOT)=0.  KBOT and KTOP together
   61: *          determine an isolated block along the diagonal of the
   62: *          Hessenberg matrix.
   63: *
   64: *     NW      (input) INTEGER
   65: *          Deflation window size.  1 .LE. NW .LE. (KBOT-KTOP+1).
   66: *
   67: *     H       (input/output) COMPLEX*16 array, dimension (LDH,N)
   68: *          On input the initial N-by-N section of H stores the
   69: *          Hessenberg matrix undergoing aggressive early deflation.
   70: *          On output H has been transformed by a unitary
   71: *          similarity transformation, perturbed, and the returned
   72: *          to Hessenberg form that (it is to be hoped) has some
   73: *          zero subdiagonal entries.
   74: *
   75: *     LDH     (input) integer
   76: *          Leading dimension of H just as declared in the calling
   77: *          subroutine.  N .LE. LDH
   78: *
   79: *     ILOZ    (input) INTEGER
   80: *     IHIZ    (input) INTEGER
   81: *          Specify the rows of Z to which transformations must be
   82: *          applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N.
   83: *
   84: *     Z       (input/output) COMPLEX*16 array, dimension (LDZ,N)
   85: *          IF WANTZ is .TRUE., then on output, the unitary
   86: *          similarity transformation mentioned above has been
   87: *          accumulated into Z(ILOZ:IHIZ,ILO:IHI) from the right.
   88: *          If WANTZ is .FALSE., then Z is unreferenced.
   89: *
   90: *     LDZ     (input) integer
   91: *          The leading dimension of Z just as declared in the
   92: *          calling subroutine.  1 .LE. LDZ.
   93: *
   94: *     NS      (output) integer
   95: *          The number of unconverged (ie approximate) eigenvalues
   96: *          returned in SR and SI that may be used as shifts by the
   97: *          calling subroutine.
   98: *
   99: *     ND      (output) integer
  100: *          The number of converged eigenvalues uncovered by this
  101: *          subroutine.
  102: *
  103: *     SH      (output) COMPLEX*16 array, dimension KBOT
  104: *          On output, approximate eigenvalues that may
  105: *          be used for shifts are stored in SH(KBOT-ND-NS+1)
  106: *          through SR(KBOT-ND).  Converged eigenvalues are
  107: *          stored in SH(KBOT-ND+1) through SH(KBOT).
  108: *
  109: *     V       (workspace) COMPLEX*16 array, dimension (LDV,NW)
  110: *          An NW-by-NW work array.
  111: *
  112: *     LDV     (input) integer scalar
  113: *          The leading dimension of V just as declared in the
  114: *          calling subroutine.  NW .LE. LDV
  115: *
  116: *     NH      (input) integer scalar
  117: *          The number of columns of T.  NH.GE.NW.
  118: *
  119: *     T       (workspace) COMPLEX*16 array, dimension (LDT,NW)
  120: *
  121: *     LDT     (input) integer
  122: *          The leading dimension of T just as declared in the
  123: *          calling subroutine.  NW .LE. LDT
  124: *
  125: *     NV      (input) integer
  126: *          The number of rows of work array WV available for
  127: *          workspace.  NV.GE.NW.
  128: *
  129: *     WV      (workspace) COMPLEX*16 array, dimension (LDWV,NW)
  130: *
  131: *     LDWV    (input) integer
  132: *          The leading dimension of W just as declared in the
  133: *          calling subroutine.  NW .LE. LDV
  134: *
  135: *     WORK    (workspace) COMPLEX*16 array, dimension LWORK.
  136: *          On exit, WORK(1) is set to an estimate of the optimal value
  137: *          of LWORK for the given values of N, NW, KTOP and KBOT.
  138: *
  139: *     LWORK   (input) integer
  140: *          The dimension of the work array WORK.  LWORK = 2*NW
  141: *          suffices, but greater efficiency may result from larger
  142: *          values of LWORK.
  143: *
  144: *          If LWORK = -1, then a workspace query is assumed; ZLAQR2
  145: *          only estimates the optimal workspace size for the given
  146: *          values of N, NW, KTOP and KBOT.  The estimate is returned
  147: *          in WORK(1).  No error message related to LWORK is issued
  148: *          by XERBLA.  Neither H nor Z are accessed.
  149: *
  150: *     ================================================================
  151: *     Based on contributions by
  152: *        Karen Braman and Ralph Byers, Department of Mathematics,
  153: *        University of Kansas, USA
  154: *
  155: *     ================================================================
  156: *     .. Parameters ..
  157:       COMPLEX*16         ZERO, ONE
  158:       PARAMETER          ( ZERO = ( 0.0d0, 0.0d0 ),
  159:      $                   ONE = ( 1.0d0, 0.0d0 ) )
  160:       DOUBLE PRECISION   RZERO, RONE
  161:       PARAMETER          ( RZERO = 0.0d0, RONE = 1.0d0 )
  162: *     ..
  163: *     .. Local Scalars ..
  164:       COMPLEX*16         BETA, CDUM, S, TAU
  165:       DOUBLE PRECISION   FOO, SAFMAX, SAFMIN, SMLNUM, ULP
  166:       INTEGER            I, IFST, ILST, INFO, INFQR, J, JW, KCOL, KLN,
  167:      $                   KNT, KROW, KWTOP, LTOP, LWK1, LWK2, LWKOPT
  168: *     ..
  169: *     .. External Functions ..
  170:       DOUBLE PRECISION   DLAMCH
  171:       EXTERNAL           DLAMCH
  172: *     ..
  173: *     .. External Subroutines ..
  174:       EXTERNAL           DLABAD, ZCOPY, ZGEHRD, ZGEMM, ZLACPY, ZLAHQR,
  175:      $                   ZLARF, ZLARFG, ZLASET, ZTREXC, ZUNMHR
  176: *     ..
  177: *     .. Intrinsic Functions ..
  178:       INTRINSIC          ABS, DBLE, DCMPLX, DCONJG, DIMAG, INT, MAX, MIN
  179: *     ..
  180: *     .. Statement Functions ..
  181:       DOUBLE PRECISION   CABS1
  182: *     ..
  183: *     .. Statement Function definitions ..
  184:       CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
  185: *     ..
  186: *     .. Executable Statements ..
  187: *
  188: *     ==== Estimate optimal workspace. ====
  189: *
  190:       JW = MIN( NW, KBOT-KTOP+1 )
  191:       IF( JW.LE.2 ) THEN
  192:          LWKOPT = 1
  193:       ELSE
  194: *
  195: *        ==== Workspace query call to ZGEHRD ====
  196: *
  197:          CALL ZGEHRD( JW, 1, JW-1, T, LDT, WORK, WORK, -1, INFO )
  198:          LWK1 = INT( WORK( 1 ) )
  199: *
  200: *        ==== Workspace query call to ZUNMHR ====
  201: *
  202:          CALL ZUNMHR( 'R', 'N', JW, JW, 1, JW-1, T, LDT, WORK, V, LDV,
  203:      $                WORK, -1, INFO )
  204:          LWK2 = INT( WORK( 1 ) )
  205: *
  206: *        ==== Optimal workspace ====
  207: *
  208:          LWKOPT = JW + MAX( LWK1, LWK2 )
  209:       END IF
  210: *
  211: *     ==== Quick return in case of workspace query. ====
  212: *
  213:       IF( LWORK.EQ.-1 ) THEN
  214:          WORK( 1 ) = DCMPLX( LWKOPT, 0 )
  215:          RETURN
  216:       END IF
  217: *
  218: *     ==== Nothing to do ...
  219: *     ... for an empty active block ... ====
  220:       NS = 0
  221:       ND = 0
  222:       WORK( 1 ) = ONE
  223:       IF( KTOP.GT.KBOT )
  224:      $   RETURN
  225: *     ... nor for an empty deflation window. ====
  226:       IF( NW.LT.1 )
  227:      $   RETURN
  228: *
  229: *     ==== Machine constants ====
  230: *
  231:       SAFMIN = DLAMCH( 'SAFE MINIMUM' )
  232:       SAFMAX = RONE / SAFMIN
  233:       CALL DLABAD( SAFMIN, SAFMAX )
  234:       ULP = DLAMCH( 'PRECISION' )
  235:       SMLNUM = SAFMIN*( DBLE( N ) / ULP )
  236: *
  237: *     ==== Setup deflation window ====
  238: *
  239:       JW = MIN( NW, KBOT-KTOP+1 )
  240:       KWTOP = KBOT - JW + 1
  241:       IF( KWTOP.EQ.KTOP ) THEN
  242:          S = ZERO
  243:       ELSE
  244:          S = H( KWTOP, KWTOP-1 )
  245:       END IF
  246: *
  247:       IF( KBOT.EQ.KWTOP ) THEN
  248: *
  249: *        ==== 1-by-1 deflation window: not much to do ====
  250: *
  251:          SH( KWTOP ) = H( KWTOP, KWTOP )
  252:          NS = 1
  253:          ND = 0
  254:          IF( CABS1( S ).LE.MAX( SMLNUM, ULP*CABS1( H( KWTOP,
  255:      $       KWTOP ) ) ) ) THEN
  256:             NS = 0
  257:             ND = 1
  258:             IF( KWTOP.GT.KTOP )
  259:      $         H( KWTOP, KWTOP-1 ) = ZERO
  260:          END IF
  261:          WORK( 1 ) = ONE
  262:          RETURN
  263:       END IF
  264: *
  265: *     ==== Convert to spike-triangular form.  (In case of a
  266: *     .    rare QR failure, this routine continues to do
  267: *     .    aggressive early deflation using that part of
  268: *     .    the deflation window that converged using INFQR
  269: *     .    here and there to keep track.) ====
  270: *
  271:       CALL ZLACPY( 'U', JW, JW, H( KWTOP, KWTOP ), LDH, T, LDT )
  272:       CALL ZCOPY( JW-1, H( KWTOP+1, KWTOP ), LDH+1, T( 2, 1 ), LDT+1 )
  273: *
  274:       CALL ZLASET( 'A', JW, JW, ZERO, ONE, V, LDV )
  275:       CALL ZLAHQR( .true., .true., JW, 1, JW, T, LDT, SH( KWTOP ), 1,
  276:      $             JW, V, LDV, INFQR )
  277: *
  278: *     ==== Deflation detection loop ====
  279: *
  280:       NS = JW
  281:       ILST = INFQR + 1
  282:       DO 10 KNT = INFQR + 1, JW
  283: *
  284: *        ==== Small spike tip deflation test ====
  285: *
  286:          FOO = CABS1( T( NS, NS ) )
  287:          IF( FOO.EQ.RZERO )
  288:      $      FOO = CABS1( S )
  289:          IF( CABS1( S )*CABS1( V( 1, NS ) ).LE.MAX( SMLNUM, ULP*FOO ) )
  290:      $        THEN
  291: *
  292: *           ==== One more converged eigenvalue ====
  293: *
  294:             NS = NS - 1
  295:          ELSE
  296: *
  297: *           ==== One undeflatable eigenvalue.  Move it up out of the
  298: *           .    way.   (ZTREXC can not fail in this case.) ====
  299: *
  300:             IFST = NS
  301:             CALL ZTREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, INFO )
  302:             ILST = ILST + 1
  303:          END IF
  304:    10 CONTINUE
  305: *
  306: *        ==== Return to Hessenberg form ====
  307: *
  308:       IF( NS.EQ.0 )
  309:      $   S = ZERO
  310: *
  311:       IF( NS.LT.JW ) THEN
  312: *
  313: *        ==== sorting the diagonal of T improves accuracy for
  314: *        .    graded matrices.  ====
  315: *
  316:          DO 30 I = INFQR + 1, NS
  317:             IFST = I
  318:             DO 20 J = I + 1, NS
  319:                IF( CABS1( T( J, J ) ).GT.CABS1( T( IFST, IFST ) ) )
  320:      $            IFST = J
  321:    20       CONTINUE
  322:             ILST = I
  323:             IF( IFST.NE.ILST )
  324:      $         CALL ZTREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, INFO )
  325:    30    CONTINUE
  326:       END IF
  327: *
  328: *     ==== Restore shift/eigenvalue array from T ====
  329: *
  330:       DO 40 I = INFQR + 1, JW
  331:          SH( KWTOP+I-1 ) = T( I, I )
  332:    40 CONTINUE
  333: *
  334: *
  335:       IF( NS.LT.JW .OR. S.EQ.ZERO ) THEN
  336:          IF( NS.GT.1 .AND. S.NE.ZERO ) THEN
  337: *
  338: *           ==== Reflect spike back into lower triangle ====
  339: *
  340:             CALL ZCOPY( NS, V, LDV, WORK, 1 )
  341:             DO 50 I = 1, NS
  342:                WORK( I ) = DCONJG( WORK( I ) )
  343:    50       CONTINUE
  344:             BETA = WORK( 1 )
  345:             CALL ZLARFG( NS, BETA, WORK( 2 ), 1, TAU )
  346:             WORK( 1 ) = ONE
  347: *
  348:             CALL ZLASET( 'L', JW-2, JW-2, ZERO, ZERO, T( 3, 1 ), LDT )
  349: *
  350:             CALL ZLARF( 'L', NS, JW, WORK, 1, DCONJG( TAU ), T, LDT,
  351:      $                  WORK( JW+1 ) )
  352:             CALL ZLARF( 'R', NS, NS, WORK, 1, TAU, T, LDT,
  353:      $                  WORK( JW+1 ) )
  354:             CALL ZLARF( 'R', JW, NS, WORK, 1, TAU, V, LDV,
  355:      $                  WORK( JW+1 ) )
  356: *
  357:             CALL ZGEHRD( JW, 1, NS, T, LDT, WORK, WORK( JW+1 ),
  358:      $                   LWORK-JW, INFO )
  359:          END IF
  360: *
  361: *        ==== Copy updated reduced window into place ====
  362: *
  363:          IF( KWTOP.GT.1 )
  364:      $      H( KWTOP, KWTOP-1 ) = S*DCONJG( V( 1, 1 ) )
  365:          CALL ZLACPY( 'U', JW, JW, T, LDT, H( KWTOP, KWTOP ), LDH )
  366:          CALL ZCOPY( JW-1, T( 2, 1 ), LDT+1, H( KWTOP+1, KWTOP ),
  367:      $               LDH+1 )
  368: *
  369: *        ==== Accumulate orthogonal matrix in order update
  370: *        .    H and Z, if requested.  ====
  371: *
  372:          IF( NS.GT.1 .AND. S.NE.ZERO )
  373:      $      CALL ZUNMHR( 'R', 'N', JW, NS, 1, NS, T, LDT, WORK, V, LDV,
  374:      $                   WORK( JW+1 ), LWORK-JW, INFO )
  375: *
  376: *        ==== Update vertical slab in H ====
  377: *
  378:          IF( WANTT ) THEN
  379:             LTOP = 1
  380:          ELSE
  381:             LTOP = KTOP
  382:          END IF
  383:          DO 60 KROW = LTOP, KWTOP - 1, NV
  384:             KLN = MIN( NV, KWTOP-KROW )
  385:             CALL ZGEMM( 'N', 'N', KLN, JW, JW, ONE, H( KROW, KWTOP ),
  386:      $                  LDH, V, LDV, ZERO, WV, LDWV )
  387:             CALL ZLACPY( 'A', KLN, JW, WV, LDWV, H( KROW, KWTOP ), LDH )
  388:    60    CONTINUE
  389: *
  390: *        ==== Update horizontal slab in H ====
  391: *
  392:          IF( WANTT ) THEN
  393:             DO 70 KCOL = KBOT + 1, N, NH
  394:                KLN = MIN( NH, N-KCOL+1 )
  395:                CALL ZGEMM( 'C', 'N', JW, KLN, JW, ONE, V, LDV,
  396:      $                     H( KWTOP, KCOL ), LDH, ZERO, T, LDT )
  397:                CALL ZLACPY( 'A', JW, KLN, T, LDT, H( KWTOP, KCOL ),
  398:      $                      LDH )
  399:    70       CONTINUE
  400:          END IF
  401: *
  402: *        ==== Update vertical slab in Z ====
  403: *
  404:          IF( WANTZ ) THEN
  405:             DO 80 KROW = ILOZ, IHIZ, NV
  406:                KLN = MIN( NV, IHIZ-KROW+1 )
  407:                CALL ZGEMM( 'N', 'N', KLN, JW, JW, ONE, Z( KROW, KWTOP ),
  408:      $                     LDZ, V, LDV, ZERO, WV, LDWV )
  409:                CALL ZLACPY( 'A', KLN, JW, WV, LDWV, Z( KROW, KWTOP ),
  410:      $                      LDZ )
  411:    80       CONTINUE
  412:          END IF
  413:       END IF
  414: *
  415: *     ==== Return the number of deflations ... ====
  416: *
  417:       ND = JW - NS
  418: *
  419: *     ==== ... and the number of shifts. (Subtracting
  420: *     .    INFQR from the spike length takes care
  421: *     .    of the case of a rare QR failure while
  422: *     .    calculating eigenvalues of the deflation
  423: *     .    window.)  ====
  424: *
  425:       NS = NS - INFQR
  426: *
  427: *      ==== Return optimal workspace. ====
  428: *
  429:       WORK( 1 ) = DCMPLX( LWKOPT, 0 )
  430: *
  431: *     ==== End of ZLAQR2 ====
  432: *
  433:       END