--- rpl/lapack/lapack/zlahqr.f 2010/08/06 15:28:56 1.3 +++ rpl/lapack/lapack/zlahqr.f 2023/08/07 08:39:29 1.20 @@ -1,9 +1,202 @@ +*> \brief \b ZLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm. +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download ZLAHQR + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE ZLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, +* IHIZ, Z, LDZ, INFO ) +* +* .. Scalar Arguments .. +* INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N +* LOGICAL WANTT, WANTZ +* .. +* .. Array Arguments .. +* COMPLEX*16 H( LDH, * ), W( * ), Z( LDZ, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> ZLAHQR is an auxiliary routine called by CHSEQR to update the +*> eigenvalues and Schur decomposition already computed by CHSEQR, by +*> dealing with the Hessenberg submatrix in rows and columns ILO to +*> IHI. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] WANTT +*> \verbatim +*> WANTT is LOGICAL +*> = .TRUE. : the full Schur form T is required; +*> = .FALSE.: only eigenvalues are required. +*> \endverbatim +*> +*> \param[in] WANTZ +*> \verbatim +*> WANTZ is LOGICAL +*> = .TRUE. : the matrix of Schur vectors Z is required; +*> = .FALSE.: Schur vectors are not required. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the matrix H. N >= 0. +*> \endverbatim +*> +*> \param[in] ILO +*> \verbatim +*> ILO is INTEGER +*> \endverbatim +*> +*> \param[in] IHI +*> \verbatim +*> IHI is INTEGER +*> It is assumed that H is already upper triangular in rows and +*> columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless ILO = 1). +*> ZLAHQR works primarily with the Hessenberg submatrix in rows +*> and columns ILO to IHI, but applies transformations to all of +*> H if WANTT is .TRUE.. +*> 1 <= ILO <= max(1,IHI); IHI <= N. +*> \endverbatim +*> +*> \param[in,out] H +*> \verbatim +*> H is COMPLEX*16 array, dimension (LDH,N) +*> On entry, the upper Hessenberg matrix H. +*> On exit, if INFO is zero and if WANTT is .TRUE., then H +*> is upper triangular in rows and columns ILO:IHI. If INFO +*> is zero and if WANTT is .FALSE., then the contents of H +*> are unspecified on exit. The output state of H in case +*> INF is positive is below under the description of INFO. +*> \endverbatim +*> +*> \param[in] LDH +*> \verbatim +*> LDH is INTEGER +*> The leading dimension of the array H. LDH >= max(1,N). +*> \endverbatim +*> +*> \param[out] W +*> \verbatim +*> W is COMPLEX*16 array, dimension (N) +*> The computed eigenvalues ILO to IHI are stored in the +*> corresponding elements of W. If WANTT is .TRUE., the +*> eigenvalues are stored in the same order as on the diagonal +*> of the Schur form returned in H, with W(i) = H(i,i). +*> \endverbatim +*> +*> \param[in] ILOZ +*> \verbatim +*> ILOZ is INTEGER +*> \endverbatim +*> +*> \param[in] IHIZ +*> \verbatim +*> IHIZ is INTEGER +*> Specify the rows of Z to which transformations must be +*> applied if WANTZ is .TRUE.. +*> 1 <= ILOZ <= ILO; IHI <= IHIZ <= N. +*> \endverbatim +*> +*> \param[in,out] Z +*> \verbatim +*> Z is COMPLEX*16 array, dimension (LDZ,N) +*> If WANTZ is .TRUE., on entry Z must contain the current +*> matrix Z of transformations accumulated by CHSEQR, and on +*> exit Z has been updated; transformations are applied only to +*> the submatrix Z(ILOZ:IHIZ,ILO:IHI). +*> If WANTZ is .FALSE., Z is not referenced. +*> \endverbatim +*> +*> \param[in] LDZ +*> \verbatim +*> LDZ is INTEGER +*> The leading dimension of the array Z. LDZ >= max(1,N). +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> > 0: if INFO = i, ZLAHQR failed to compute all the +*> eigenvalues ILO to IHI in a total of 30 iterations +*> per eigenvalue; elements i+1:ihi of W contain +*> those eigenvalues which have been successfully +*> computed. +*> +*> If INFO > 0 and WANTT is .FALSE., then on exit, +*> the remaining unconverged eigenvalues are the +*> eigenvalues of the upper Hessenberg matrix +*> rows and columns ILO through INFO of the final, +*> output value of H. +*> +*> If INFO > 0 and WANTT is .TRUE., then on exit +*> (*) (initial value of H)*U = U*(final value of H) +*> where U is an orthogonal matrix. The final +*> value of H is upper Hessenberg and triangular in +*> rows and columns INFO+1 through IHI. +*> +*> If INFO > 0 and WANTZ is .TRUE., then on exit +*> (final value of Z) = (initial value of Z)*U +*> where U is the orthogonal matrix in (*) +*> (regardless of the value of WANTT.) +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \ingroup complex16OTHERauxiliary +* +*> \par Contributors: +* ================== +*> +*> \verbatim +*> +*> 02-96 Based on modifications by +*> David Day, Sandia National Laboratory, USA +*> +*> 12-04 Further modifications by +*> Ralph Byers, University of Kansas, USA +*> This is a modified version of ZLAHQR from LAPACK version 3.0. +*> It is (1) more robust against overflow and underflow and +*> (2) adopts the more conservative Ahues & Tisseur stopping +*> criterion (LAWN 122, 1997). +*> \endverbatim +*> +* ===================================================================== SUBROUTINE ZLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, $ IHIZ, Z, LDZ, INFO ) + IMPLICIT NONE * -* -- LAPACK auxiliary routine (version 3.2) -- -* Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd.. -* November 2006 +* -- LAPACK auxiliary routine -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * * .. Scalar Arguments .. INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N @@ -13,113 +206,9 @@ COMPLEX*16 H( LDH, * ), W( * ), Z( LDZ, * ) * .. * -* Purpose -* ======= -* -* ZLAHQR is an auxiliary routine called by CHSEQR to update the -* eigenvalues and Schur decomposition already computed by CHSEQR, by -* dealing with the Hessenberg submatrix in rows and columns ILO to -* IHI. -* -* Arguments -* ========= -* -* WANTT (input) LOGICAL -* = .TRUE. : the full Schur form T is required; -* = .FALSE.: only eigenvalues are required. -* -* WANTZ (input) LOGICAL -* = .TRUE. : the matrix of Schur vectors Z is required; -* = .FALSE.: Schur vectors are not required. -* -* N (input) INTEGER -* The order of the matrix H. N >= 0. -* -* ILO (input) INTEGER -* IHI (input) INTEGER -* It is assumed that H is already upper triangular in rows and -* columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless ILO = 1). -* ZLAHQR works primarily with the Hessenberg submatrix in rows -* and columns ILO to IHI, but applies transformations to all of -* H if WANTT is .TRUE.. -* 1 <= ILO <= max(1,IHI); IHI <= N. -* -* H (input/output) COMPLEX*16 array, dimension (LDH,N) -* On entry, the upper Hessenberg matrix H. -* On exit, if INFO is zero and if WANTT is .TRUE., then H -* is upper triangular in rows and columns ILO:IHI. If INFO -* is zero and if WANTT is .FALSE., then the contents of H -* are unspecified on exit. The output state of H in case -* INF is positive is below under the description of INFO. -* -* LDH (input) INTEGER -* The leading dimension of the array H. LDH >= max(1,N). -* -* W (output) COMPLEX*16 array, dimension (N) -* The computed eigenvalues ILO to IHI are stored in the -* corresponding elements of W. If WANTT is .TRUE., the -* eigenvalues are stored in the same order as on the diagonal -* of the Schur form returned in H, with W(i) = H(i,i). -* -* ILOZ (input) INTEGER -* IHIZ (input) INTEGER -* Specify the rows of Z to which transformations must be -* applied if WANTZ is .TRUE.. -* 1 <= ILOZ <= ILO; IHI <= IHIZ <= N. -* -* Z (input/output) COMPLEX*16 array, dimension (LDZ,N) -* If WANTZ is .TRUE., on entry Z must contain the current -* matrix Z of transformations accumulated by CHSEQR, and on -* exit Z has been updated; transformations are applied only to -* the submatrix Z(ILOZ:IHIZ,ILO:IHI). -* If WANTZ is .FALSE., Z is not referenced. -* -* LDZ (input) INTEGER -* The leading dimension of the array Z. LDZ >= max(1,N). -* -* INFO (output) INTEGER -* = 0: successful exit -* .GT. 0: if INFO = i, ZLAHQR failed to compute all the -* eigenvalues ILO to IHI in a total of 30 iterations -* per eigenvalue; elements i+1:ihi of W contain -* those eigenvalues which have been successfully -* computed. -* -* If INFO .GT. 0 and WANTT is .FALSE., then on exit, -* the remaining unconverged eigenvalues are the -* eigenvalues of the upper Hessenberg matrix -* rows and columns ILO thorugh INFO of the final, -* output value of H. -* -* If INFO .GT. 0 and WANTT is .TRUE., then on exit -* (*) (initial value of H)*U = U*(final value of H) -* where U is an orthognal matrix. The final -* value of H is upper Hessenberg and triangular in -* rows and columns INFO+1 through IHI. -* -* If INFO .GT. 0 and WANTZ is .TRUE., then on exit -* (final value of Z) = (initial value of Z)*U -* where U is the orthogonal matrix in (*) -* (regardless of the value of WANTT.) -* -* Further Details -* =============== -* -* 02-96 Based on modifications by -* David Day, Sandia National Laboratory, USA -* -* 12-04 Further modifications by -* Ralph Byers, University of Kansas, USA -* This is a modified version of ZLAHQR from LAPACK version 3.0. -* It is (1) more robust against overflow and underflow and -* (2) adopts the more conservative Ahues & Tisseur stopping -* criterion (LAWN 122, 1997). -* -* ========================================================= +* ========================================================= * * .. Parameters .. - INTEGER ITMAX - PARAMETER ( ITMAX = 30 ) COMPLEX*16 ZERO, ONE PARAMETER ( ZERO = ( 0.0d0, 0.0d0 ), $ ONE = ( 1.0d0, 0.0d0 ) ) @@ -127,13 +216,16 @@ PARAMETER ( RZERO = 0.0d0, RONE = 1.0d0, HALF = 0.5d0 ) DOUBLE PRECISION DAT1 PARAMETER ( DAT1 = 3.0d0 / 4.0d0 ) + INTEGER KEXSH + PARAMETER ( KEXSH = 10 ) * .. * .. Local Scalars .. COMPLEX*16 CDUM, H11, H11S, H22, SC, SUM, T, T1, TEMP, U, $ V2, X, Y DOUBLE PRECISION AA, AB, BA, BB, H10, H21, RTEMP, S, SAFMAX, $ SAFMIN, SMLNUM, SX, T2, TST, ULP - INTEGER I, I1, I2, ITS, J, JHI, JLO, K, L, M, NH, NZ + INTEGER I, I1, I2, ITS, ITMAX, J, JHI, JLO, K, L, M, + $ NH, NZ, KDEFL * .. * .. Local Arrays .. COMPLEX*16 V( 2 ) @@ -219,6 +311,14 @@ I2 = N END IF * +* ITMAX is the total number of QR iterations allowed. +* + ITMAX = 30 * MAX( 10, NH ) +* +* KDEFL counts the number of iterations since a deflation +* + KDEFL = 0 +* * The main loop begins here. I is the loop index and decreases from * IHI to ILO in steps of 1. Each iteration of the loop works * with the active submatrix in rows and columns L to I. @@ -278,6 +378,7 @@ * IF( L.GE.I ) $ GO TO 140 + KDEFL = KDEFL + 1 * * Now the active submatrix is in rows and columns L to I. If * eigenvalues only are being computed, only the active submatrix @@ -288,18 +389,18 @@ I2 = I END IF * - IF( ITS.EQ.10 ) THEN + IF( MOD(KDEFL,2*KEXSH).EQ.0 ) THEN * * Exceptional shift. * - S = DAT1*ABS( DBLE( H( L+1, L ) ) ) - T = S + H( L, L ) - ELSE IF( ITS.EQ.20 ) THEN + S = DAT1*ABS( DBLE( H( I, I-1 ) ) ) + T = S + H( I, I ) + ELSE IF( MOD(KDEFL,KEXSH).EQ.0 ) THEN * * Exceptional shift. * - S = DAT1*ABS( DBLE( H( I, I-1 ) ) ) - T = S + H( I, I ) + S = DAT1*ABS( DBLE( H( L+1, L ) ) ) + T = S + H( L, L ) ELSE * * Wilkinson's shift. @@ -461,6 +562,8 @@ * H(I,I-1) is negligible: one eigenvalue has converged. * W( I ) = H( I, I ) +* reset deflation counter + KDEFL = 0 * * return to start of the main loop with new value of I. *