LAPACK 3.3.0

claqr2.f

Go to the documentation of this file.
00001       SUBROUTINE CLAQR2( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ,
00002      $                   IHIZ, Z, LDZ, NS, ND, SH, V, LDV, NH, T, LDT,
00003      $                   NV, WV, LDWV, WORK, LWORK )
00004 *
00005 *  -- LAPACK auxiliary routine (version 3.2.1)                        --
00006 *     Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
00007 *  -- April 2009                                                      --
00008 *
00009 *     .. Scalar Arguments ..
00010       INTEGER            IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV,
00011      $                   LDZ, LWORK, N, ND, NH, NS, NV, NW
00012       LOGICAL            WANTT, WANTZ
00013 *     ..
00014 *     .. Array Arguments ..
00015       COMPLEX            H( LDH, * ), SH( * ), T( LDT, * ), V( LDV, * ),
00016      $                   WORK( * ), WV( LDWV, * ), Z( LDZ, * )
00017 *     ..
00018 *
00019 *     This subroutine is identical to CLAQR3 except that it avoids
00020 *     recursion by calling CLAHQR instead of CLAQR4.
00021 *
00022 *
00023 *     ******************************************************************
00024 *     Aggressive early deflation:
00025 *
00026 *     This subroutine accepts as input an upper Hessenberg matrix
00027 *     H and performs an unitary similarity transformation
00028 *     designed to detect and deflate fully converged eigenvalues from
00029 *     a trailing principal submatrix.  On output H has been over-
00030 *     written by a new Hessenberg matrix that is a perturbation of
00031 *     an unitary similarity transformation of H.  It is to be
00032 *     hoped that the final version of H has many zero subdiagonal
00033 *     entries.
00034 *
00035 *     ******************************************************************
00036 *     WANTT   (input) LOGICAL
00037 *          If .TRUE., then the Hessenberg matrix H is fully updated
00038 *          so that the triangular Schur factor may be
00039 *          computed (in cooperation with the calling subroutine).
00040 *          If .FALSE., then only enough of H is updated to preserve
00041 *          the eigenvalues.
00042 *
00043 *     WANTZ   (input) LOGICAL
00044 *          If .TRUE., then the unitary matrix Z is updated so
00045 *          so that the unitary Schur factor may be computed
00046 *          (in cooperation with the calling subroutine).
00047 *          If .FALSE., then Z is not referenced.
00048 *
00049 *     N       (input) INTEGER
00050 *          The order of the matrix H and (if WANTZ is .TRUE.) the
00051 *          order of the unitary matrix Z.
00052 *
00053 *     KTOP    (input) INTEGER
00054 *          It is assumed that either KTOP = 1 or H(KTOP,KTOP-1)=0.
00055 *          KBOT and KTOP together determine an isolated block
00056 *          along the diagonal of the Hessenberg matrix.
00057 *
00058 *     KBOT    (input) INTEGER
00059 *          It is assumed without a check that either
00060 *          KBOT = N or H(KBOT+1,KBOT)=0.  KBOT and KTOP together
00061 *          determine an isolated block along the diagonal of the
00062 *          Hessenberg matrix.
00063 *
00064 *     NW      (input) INTEGER
00065 *          Deflation window size.  1 .LE. NW .LE. (KBOT-KTOP+1).
00066 *
00067 *     H       (input/output) COMPLEX array, dimension (LDH,N)
00068 *          On input the initial N-by-N section of H stores the
00069 *          Hessenberg matrix undergoing aggressive early deflation.
00070 *          On output H has been transformed by a unitary
00071 *          similarity transformation, perturbed, and the returned
00072 *          to Hessenberg form that (it is to be hoped) has some
00073 *          zero subdiagonal entries.
00074 *
00075 *     LDH     (input) integer
00076 *          Leading dimension of H just as declared in the calling
00077 *          subroutine.  N .LE. LDH
00078 *
00079 *     ILOZ    (input) INTEGER
00080 *     IHIZ    (input) INTEGER
00081 *          Specify the rows of Z to which transformations must be
00082 *          applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N.
00083 *
00084 *     Z       (input/output) COMPLEX array, dimension (LDZ,N)
00085 *          IF WANTZ is .TRUE., then on output, the unitary
00086 *          similarity transformation mentioned above has been
00087 *          accumulated into Z(ILOZ:IHIZ,ILO:IHI) from the right.
00088 *          If WANTZ is .FALSE., then Z is unreferenced.
00089 *
00090 *     LDZ     (input) integer
00091 *          The leading dimension of Z just as declared in the
00092 *          calling subroutine.  1 .LE. LDZ.
00093 *
00094 *     NS      (output) integer
00095 *          The number of unconverged (ie approximate) eigenvalues
00096 *          returned in SR and SI that may be used as shifts by the
00097 *          calling subroutine.
00098 *
00099 *     ND      (output) integer
00100 *          The number of converged eigenvalues uncovered by this
00101 *          subroutine.
00102 *
00103 *     SH      (output) COMPLEX array, dimension KBOT
00104 *          On output, approximate eigenvalues that may
00105 *          be used for shifts are stored in SH(KBOT-ND-NS+1)
00106 *          through SR(KBOT-ND).  Converged eigenvalues are
00107 *          stored in SH(KBOT-ND+1) through SH(KBOT).
00108 *
00109 *     V       (workspace) COMPLEX array, dimension (LDV,NW)
00110 *          An NW-by-NW work array.
00111 *
00112 *     LDV     (input) integer scalar
00113 *          The leading dimension of V just as declared in the
00114 *          calling subroutine.  NW .LE. LDV
00115 *
00116 *     NH      (input) integer scalar
00117 *          The number of columns of T.  NH.GE.NW.
00118 *
00119 *     T       (workspace) COMPLEX array, dimension (LDT,NW)
00120 *
00121 *     LDT     (input) integer
00122 *          The leading dimension of T just as declared in the
00123 *          calling subroutine.  NW .LE. LDT
00124 *
00125 *     NV      (input) integer
00126 *          The number of rows of work array WV available for
00127 *          workspace.  NV.GE.NW.
00128 *
00129 *     WV      (workspace) COMPLEX array, dimension (LDWV,NW)
00130 *
00131 *     LDWV    (input) integer
00132 *          The leading dimension of W just as declared in the
00133 *          calling subroutine.  NW .LE. LDV
00134 *
00135 *     WORK    (workspace) COMPLEX array, dimension LWORK.
00136 *          On exit, WORK(1) is set to an estimate of the optimal value
00137 *          of LWORK for the given values of N, NW, KTOP and KBOT.
00138 *
00139 *     LWORK   (input) integer
00140 *          The dimension of the work array WORK.  LWORK = 2*NW
00141 *          suffices, but greater efficiency may result from larger
00142 *          values of LWORK.
00143 *
00144 *          If LWORK = -1, then a workspace query is assumed; CLAQR2
00145 *          only estimates the optimal workspace size for the given
00146 *          values of N, NW, KTOP and KBOT.  The estimate is returned
00147 *          in WORK(1).  No error message related to LWORK is issued
00148 *          by XERBLA.  Neither H nor Z are accessed.
00149 *
00150 *     ================================================================
00151 *     Based on contributions by
00152 *        Karen Braman and Ralph Byers, Department of Mathematics,
00153 *        University of Kansas, USA
00154 *
00155 *     ================================================================
00156 *     .. Parameters ..
00157       COMPLEX            ZERO, ONE
00158       PARAMETER          ( ZERO = ( 0.0e0, 0.0e0 ),
00159      $                   ONE = ( 1.0e0, 0.0e0 ) )
00160       REAL               RZERO, RONE
00161       PARAMETER          ( RZERO = 0.0e0, RONE = 1.0e0 )
00162 *     ..
00163 *     .. Local Scalars ..
00164       COMPLEX            BETA, CDUM, S, TAU
00165       REAL               FOO, SAFMAX, SAFMIN, SMLNUM, ULP
00166       INTEGER            I, IFST, ILST, INFO, INFQR, J, JW, KCOL, KLN,
00167      $                   KNT, KROW, KWTOP, LTOP, LWK1, LWK2, LWKOPT
00168 *     ..
00169 *     .. External Functions ..
00170       REAL               SLAMCH
00171       EXTERNAL           SLAMCH
00172 *     ..
00173 *     .. External Subroutines ..
00174       EXTERNAL           CCOPY, CGEHRD, CGEMM, CLACPY, CLAHQR, CLARF,
00175      $                   CLARFG, CLASET, CTREXC, CUNMHR, SLABAD
00176 *     ..
00177 *     .. Intrinsic Functions ..
00178       INTRINSIC          ABS, AIMAG, CMPLX, CONJG, INT, MAX, MIN, REAL
00179 *     ..
00180 *     .. Statement Functions ..
00181       REAL               CABS1
00182 *     ..
00183 *     .. Statement Function definitions ..
00184       CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) )
00185 *     ..
00186 *     .. Executable Statements ..
00187 *
00188 *     ==== Estimate optimal workspace. ====
00189 *
00190       JW = MIN( NW, KBOT-KTOP+1 )
00191       IF( JW.LE.2 ) THEN
00192          LWKOPT = 1
00193       ELSE
00194 *
00195 *        ==== Workspace query call to CGEHRD ====
00196 *
00197          CALL CGEHRD( JW, 1, JW-1, T, LDT, WORK, WORK, -1, INFO )
00198          LWK1 = INT( WORK( 1 ) )
00199 *
00200 *        ==== Workspace query call to CUNMHR ====
00201 *
00202          CALL CUNMHR( 'R', 'N', JW, JW, 1, JW-1, T, LDT, WORK, V, LDV,
00203      $                WORK, -1, INFO )
00204          LWK2 = INT( WORK( 1 ) )
00205 *
00206 *        ==== Optimal workspace ====
00207 *
00208          LWKOPT = JW + MAX( LWK1, LWK2 )
00209       END IF
00210 *
00211 *     ==== Quick return in case of workspace query. ====
00212 *
00213       IF( LWORK.EQ.-1 ) THEN
00214          WORK( 1 ) = CMPLX( LWKOPT, 0 )
00215          RETURN
00216       END IF
00217 *
00218 *     ==== Nothing to do ...
00219 *     ... for an empty active block ... ====
00220       NS = 0
00221       ND = 0
00222       WORK( 1 ) = ONE
00223       IF( KTOP.GT.KBOT )
00224      $   RETURN
00225 *     ... nor for an empty deflation window. ====
00226       IF( NW.LT.1 )
00227      $   RETURN
00228 *
00229 *     ==== Machine constants ====
00230 *
00231       SAFMIN = SLAMCH( 'SAFE MINIMUM' )
00232       SAFMAX = RONE / SAFMIN
00233       CALL SLABAD( SAFMIN, SAFMAX )
00234       ULP = SLAMCH( 'PRECISION' )
00235       SMLNUM = SAFMIN*( REAL( N ) / ULP )
00236 *
00237 *     ==== Setup deflation window ====
00238 *
00239       JW = MIN( NW, KBOT-KTOP+1 )
00240       KWTOP = KBOT - JW + 1
00241       IF( KWTOP.EQ.KTOP ) THEN
00242          S = ZERO
00243       ELSE
00244          S = H( KWTOP, KWTOP-1 )
00245       END IF
00246 *
00247       IF( KBOT.EQ.KWTOP ) THEN
00248 *
00249 *        ==== 1-by-1 deflation window: not much to do ====
00250 *
00251          SH( KWTOP ) = H( KWTOP, KWTOP )
00252          NS = 1
00253          ND = 0
00254          IF( CABS1( S ).LE.MAX( SMLNUM, ULP*CABS1( H( KWTOP,
00255      $       KWTOP ) ) ) ) THEN
00256             NS = 0
00257             ND = 1
00258             IF( KWTOP.GT.KTOP )
00259      $         H( KWTOP, KWTOP-1 ) = ZERO
00260          END IF
00261          WORK( 1 ) = ONE
00262          RETURN
00263       END IF
00264 *
00265 *     ==== Convert to spike-triangular form.  (In case of a
00266 *     .    rare QR failure, this routine continues to do
00267 *     .    aggressive early deflation using that part of
00268 *     .    the deflation window that converged using INFQR
00269 *     .    here and there to keep track.) ====
00270 *
00271       CALL CLACPY( 'U', JW, JW, H( KWTOP, KWTOP ), LDH, T, LDT )
00272       CALL CCOPY( JW-1, H( KWTOP+1, KWTOP ), LDH+1, T( 2, 1 ), LDT+1 )
00273 *
00274       CALL CLASET( 'A', JW, JW, ZERO, ONE, V, LDV )
00275       CALL CLAHQR( .true., .true., JW, 1, JW, T, LDT, SH( KWTOP ), 1,
00276      $             JW, V, LDV, INFQR )
00277 *
00278 *     ==== Deflation detection loop ====
00279 *
00280       NS = JW
00281       ILST = INFQR + 1
00282       DO 10 KNT = INFQR + 1, JW
00283 *
00284 *        ==== Small spike tip deflation test ====
00285 *
00286          FOO = CABS1( T( NS, NS ) )
00287          IF( FOO.EQ.RZERO )
00288      $      FOO = CABS1( S )
00289          IF( CABS1( S )*CABS1( V( 1, NS ) ).LE.MAX( SMLNUM, ULP*FOO ) )
00290      $        THEN
00291 *
00292 *           ==== One more converged eigenvalue ====
00293 *
00294             NS = NS - 1
00295          ELSE
00296 *
00297 *           ==== One undeflatable eigenvalue.  Move it up out of the
00298 *           .    way.   (CTREXC can not fail in this case.) ====
00299 *
00300             IFST = NS
00301             CALL CTREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, INFO )
00302             ILST = ILST + 1
00303          END IF
00304    10 CONTINUE
00305 *
00306 *        ==== Return to Hessenberg form ====
00307 *
00308       IF( NS.EQ.0 )
00309      $   S = ZERO
00310 *
00311       IF( NS.LT.JW ) THEN
00312 *
00313 *        ==== sorting the diagonal of T improves accuracy for
00314 *        .    graded matrices.  ====
00315 *
00316          DO 30 I = INFQR + 1, NS
00317             IFST = I
00318             DO 20 J = I + 1, NS
00319                IF( CABS1( T( J, J ) ).GT.CABS1( T( IFST, IFST ) ) )
00320      $            IFST = J
00321    20       CONTINUE
00322             ILST = I
00323             IF( IFST.NE.ILST )
00324      $         CALL CTREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, INFO )
00325    30    CONTINUE
00326       END IF
00327 *
00328 *     ==== Restore shift/eigenvalue array from T ====
00329 *
00330       DO 40 I = INFQR + 1, JW
00331          SH( KWTOP+I-1 ) = T( I, I )
00332    40 CONTINUE
00333 *
00334 *
00335       IF( NS.LT.JW .OR. S.EQ.ZERO ) THEN
00336          IF( NS.GT.1 .AND. S.NE.ZERO ) THEN
00337 *
00338 *           ==== Reflect spike back into lower triangle ====
00339 *
00340             CALL CCOPY( NS, V, LDV, WORK, 1 )
00341             DO 50 I = 1, NS
00342                WORK( I ) = CONJG( WORK( I ) )
00343    50       CONTINUE
00344             BETA = WORK( 1 )
00345             CALL CLARFG( NS, BETA, WORK( 2 ), 1, TAU )
00346             WORK( 1 ) = ONE
00347 *
00348             CALL CLASET( 'L', JW-2, JW-2, ZERO, ZERO, T( 3, 1 ), LDT )
00349 *
00350             CALL CLARF( 'L', NS, JW, WORK, 1, CONJG( TAU ), T, LDT,
00351      $                  WORK( JW+1 ) )
00352             CALL CLARF( 'R', NS, NS, WORK, 1, TAU, T, LDT,
00353      $                  WORK( JW+1 ) )
00354             CALL CLARF( 'R', JW, NS, WORK, 1, TAU, V, LDV,
00355      $                  WORK( JW+1 ) )
00356 *
00357             CALL CGEHRD( JW, 1, NS, T, LDT, WORK, WORK( JW+1 ),
00358      $                   LWORK-JW, INFO )
00359          END IF
00360 *
00361 *        ==== Copy updated reduced window into place ====
00362 *
00363          IF( KWTOP.GT.1 )
00364      $      H( KWTOP, KWTOP-1 ) = S*CONJG( V( 1, 1 ) )
00365          CALL CLACPY( 'U', JW, JW, T, LDT, H( KWTOP, KWTOP ), LDH )
00366          CALL CCOPY( JW-1, T( 2, 1 ), LDT+1, H( KWTOP+1, KWTOP ),
00367      $               LDH+1 )
00368 *
00369 *        ==== Accumulate orthogonal matrix in order update
00370 *        .    H and Z, if requested.  ====
00371 *
00372          IF( NS.GT.1 .AND. S.NE.ZERO )
00373      $      CALL CUNMHR( 'R', 'N', JW, NS, 1, NS, T, LDT, WORK, V, LDV,
00374      $                   WORK( JW+1 ), LWORK-JW, INFO )
00375 *
00376 *        ==== Update vertical slab in H ====
00377 *
00378          IF( WANTT ) THEN
00379             LTOP = 1
00380          ELSE
00381             LTOP = KTOP
00382          END IF
00383          DO 60 KROW = LTOP, KWTOP - 1, NV
00384             KLN = MIN( NV, KWTOP-KROW )
00385             CALL CGEMM( 'N', 'N', KLN, JW, JW, ONE, H( KROW, KWTOP ),
00386      $                  LDH, V, LDV, ZERO, WV, LDWV )
00387             CALL CLACPY( 'A', KLN, JW, WV, LDWV, H( KROW, KWTOP ), LDH )
00388    60    CONTINUE
00389 *
00390 *        ==== Update horizontal slab in H ====
00391 *
00392          IF( WANTT ) THEN
00393             DO 70 KCOL = KBOT + 1, N, NH
00394                KLN = MIN( NH, N-KCOL+1 )
00395                CALL CGEMM( 'C', 'N', JW, KLN, JW, ONE, V, LDV,
00396      $                     H( KWTOP, KCOL ), LDH, ZERO, T, LDT )
00397                CALL CLACPY( 'A', JW, KLN, T, LDT, H( KWTOP, KCOL ),
00398      $                      LDH )
00399    70       CONTINUE
00400          END IF
00401 *
00402 *        ==== Update vertical slab in Z ====
00403 *
00404          IF( WANTZ ) THEN
00405             DO 80 KROW = ILOZ, IHIZ, NV
00406                KLN = MIN( NV, IHIZ-KROW+1 )
00407                CALL CGEMM( 'N', 'N', KLN, JW, JW, ONE, Z( KROW, KWTOP ),
00408      $                     LDZ, V, LDV, ZERO, WV, LDWV )
00409                CALL CLACPY( 'A', KLN, JW, WV, LDWV, Z( KROW, KWTOP ),
00410      $                      LDZ )
00411    80       CONTINUE
00412          END IF
00413       END IF
00414 *
00415 *     ==== Return the number of deflations ... ====
00416 *
00417       ND = JW - NS
00418 *
00419 *     ==== ... and the number of shifts. (Subtracting
00420 *     .    INFQR from the spike length takes care
00421 *     .    of the case of a rare QR failure while
00422 *     .    calculating eigenvalues of the deflation
00423 *     .    window.)  ====
00424 *
00425       NS = NS - INFQR
00426 *
00427 *      ==== Return optimal workspace. ====
00428 *
00429       WORK( 1 ) = CMPLX( LWKOPT, 0 )
00430 *
00431 *     ==== End of CLAQR2 ====
00432 *
00433       END
 All Files Functions