version 1.14, 2016/08/27 15:34:29
|
version 1.19, 2023/08/07 08:38:56
|
Line 2
|
Line 2
|
* |
* |
* =========== DOCUMENTATION =========== |
* =========== DOCUMENTATION =========== |
* |
* |
* Online html documentation available at |
* Online html documentation available at |
* http://www.netlib.org/lapack/explore-html/ |
* http://www.netlib.org/lapack/explore-html/ |
* |
* |
*> \htmlonly |
*> \htmlonly |
*> Download DLAQR4 + dependencies |
*> Download DLAQR4 + dependencies |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaqr4.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaqr4.f"> |
*> [TGZ]</a> |
*> [TGZ]</a> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaqr4.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaqr4.f"> |
*> [ZIP]</a> |
*> [ZIP]</a> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaqr4.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaqr4.f"> |
*> [TXT]</a> |
*> [TXT]</a> |
*> \endhtmlonly |
*> \endhtmlonly |
* |
* |
* Definition: |
* Definition: |
* =========== |
* =========== |
* |
* |
* SUBROUTINE DLAQR4( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, |
* SUBROUTINE DLAQR4( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, |
* ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO ) |
* ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO ) |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
* INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, LWORK, N |
* INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, LWORK, N |
* LOGICAL WANTT, WANTZ |
* LOGICAL WANTT, WANTZ |
Line 29
|
Line 29
|
* DOUBLE PRECISION H( LDH, * ), WI( * ), WORK( * ), WR( * ), |
* DOUBLE PRECISION H( LDH, * ), WI( * ), WORK( * ), WR( * ), |
* $ Z( LDZ, * ) |
* $ Z( LDZ, * ) |
* .. |
* .. |
* |
* |
* |
* |
*> \par Purpose: |
*> \par Purpose: |
* ============= |
* ============= |
Line 74
|
Line 74
|
*> \param[in] N |
*> \param[in] N |
*> \verbatim |
*> \verbatim |
*> N is INTEGER |
*> N is INTEGER |
*> The order of the matrix H. N .GE. 0. |
*> The order of the matrix H. N >= 0. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] ILO |
*> \param[in] ILO |
Line 86
|
Line 86
|
*> \verbatim |
*> \verbatim |
*> IHI is INTEGER |
*> IHI is INTEGER |
*> It is assumed that H is already upper triangular in rows |
*> It is assumed that H is already upper triangular in rows |
*> and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1, |
*> and columns 1:ILO-1 and IHI+1:N and, if ILO > 1, |
*> H(ILO,ILO-1) is zero. ILO and IHI are normally set by a |
*> H(ILO,ILO-1) is zero. ILO and IHI are normally set by a |
*> previous call to DGEBAL, and then passed to DGEHRD when the |
*> previous call to DGEBAL, and then passed to DGEHRD when the |
*> matrix output by DGEBAL is reduced to Hessenberg form. |
*> matrix output by DGEBAL is reduced to Hessenberg form. |
*> Otherwise, ILO and IHI should be set to 1 and N, |
*> Otherwise, ILO and IHI should be set to 1 and N, |
*> respectively. If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N. |
*> respectively. If N > 0, then 1 <= ILO <= IHI <= N. |
*> If N = 0, then ILO = 1 and IHI = 0. |
*> If N = 0, then ILO = 1 and IHI = 0. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
Line 104
|
Line 104
|
*> decomposition (the Schur form); 2-by-2 diagonal blocks |
*> decomposition (the Schur form); 2-by-2 diagonal blocks |
*> (corresponding to complex conjugate pairs of eigenvalues) |
*> (corresponding to complex conjugate pairs of eigenvalues) |
*> are returned in standard form, with H(i,i) = H(i+1,i+1) |
*> are returned in standard form, with H(i,i) = H(i+1,i+1) |
*> and H(i+1,i)*H(i,i+1).LT.0. If INFO = 0 and WANTT is |
*> and H(i+1,i)*H(i,i+1) < 0. If INFO = 0 and WANTT is |
*> .FALSE., then the contents of H are unspecified on exit. |
*> .FALSE., then the contents of H are unspecified on exit. |
*> (The output value of H when INFO.GT.0 is given under the |
*> (The output value of H when INFO > 0 is given under the |
*> description of INFO below.) |
*> description of INFO below.) |
*> |
*> |
*> This subroutine may explicitly set H(i,j) = 0 for i.GT.j and |
*> This subroutine may explicitly set H(i,j) = 0 for i > j and |
*> j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N. |
*> j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] LDH |
*> \param[in] LDH |
*> \verbatim |
*> \verbatim |
*> LDH is INTEGER |
*> LDH is INTEGER |
*> The leading dimension of the array H. LDH .GE. max(1,N). |
*> The leading dimension of the array H. LDH >= max(1,N). |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[out] WR |
*> \param[out] WR |
Line 132
|
Line 132
|
*> and WI(ILO:IHI). If two eigenvalues are computed as a |
*> and WI(ILO:IHI). If two eigenvalues are computed as a |
*> complex conjugate pair, they are stored in consecutive |
*> complex conjugate pair, they are stored in consecutive |
*> elements of WR and WI, say the i-th and (i+1)th, with |
*> elements of WR and WI, say the i-th and (i+1)th, with |
*> WI(i) .GT. 0 and WI(i+1) .LT. 0. If WANTT is .TRUE., then |
*> WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., then |
*> the eigenvalues are stored in the same order as on the |
*> the eigenvalues are stored in the same order as on the |
*> diagonal of the Schur form returned in H, with |
*> diagonal of the Schur form returned in H, with |
*> WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2 diagonal |
*> WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2 diagonal |
Line 150
|
Line 150
|
*> IHIZ is INTEGER |
*> IHIZ is INTEGER |
*> Specify the rows of Z to which transformations must be |
*> Specify the rows of Z to which transformations must be |
*> applied if WANTZ is .TRUE.. |
*> applied if WANTZ is .TRUE.. |
*> 1 .LE. ILOZ .LE. ILO; IHI .LE. IHIZ .LE. N. |
*> 1 <= ILOZ <= ILO; IHI <= IHIZ <= N. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in,out] Z |
*> \param[in,out] Z |
Line 160
|
Line 160
|
*> If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is |
*> If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is |
*> replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the |
*> replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the |
*> orthogonal Schur factor of H(ILO:IHI,ILO:IHI). |
*> orthogonal Schur factor of H(ILO:IHI,ILO:IHI). |
*> (The output value of Z when INFO.GT.0 is given under |
*> (The output value of Z when INFO > 0 is given under |
*> the description of INFO below.) |
*> the description of INFO below.) |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
Line 168
|
Line 168
|
*> \verbatim |
*> \verbatim |
*> LDZ is INTEGER |
*> LDZ is INTEGER |
*> The leading dimension of the array Z. if WANTZ is .TRUE. |
*> The leading dimension of the array Z. if WANTZ is .TRUE. |
*> then LDZ.GE.MAX(1,IHIZ). Otherwize, LDZ.GE.1. |
*> then LDZ >= MAX(1,IHIZ). Otherwise, LDZ >= 1. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[out] WORK |
*> \param[out] WORK |
Line 181
|
Line 181
|
*> \param[in] LWORK |
*> \param[in] LWORK |
*> \verbatim |
*> \verbatim |
*> LWORK is INTEGER |
*> LWORK is INTEGER |
*> The dimension of the array WORK. LWORK .GE. max(1,N) |
*> The dimension of the array WORK. LWORK >= max(1,N) |
*> is sufficient, but LWORK typically as large as 6*N may |
*> is sufficient, but LWORK typically as large as 6*N may |
*> be required for optimal performance. A workspace query |
*> be required for optimal performance. A workspace query |
*> to determine the optimal workspace size is recommended. |
*> to determine the optimal workspace size is recommended. |
Line 197
|
Line 197
|
*> \param[out] INFO |
*> \param[out] INFO |
*> \verbatim |
*> \verbatim |
*> INFO is INTEGER |
*> INFO is INTEGER |
*> = 0: successful exit |
*> = 0: successful exit |
*> .GT. 0: if INFO = i, DLAQR4 failed to compute all of |
*> > 0: if INFO = i, DLAQR4 failed to compute all of |
*> the eigenvalues. Elements 1:ilo-1 and i+1:n of WR |
*> the eigenvalues. Elements 1:ilo-1 and i+1:n of WR |
*> and WI contain those eigenvalues which have been |
*> and WI contain those eigenvalues which have been |
*> successfully computed. (Failures are rare.) |
*> successfully computed. (Failures are rare.) |
*> |
*> |
*> If INFO .GT. 0 and WANT is .FALSE., then on exit, |
*> If INFO > 0 and WANT is .FALSE., then on exit, |
*> the remaining unconverged eigenvalues are the eigen- |
*> the remaining unconverged eigenvalues are the eigen- |
*> values of the upper Hessenberg matrix rows and |
*> values of the upper Hessenberg matrix rows and |
*> columns ILO through INFO of the final, output |
*> columns ILO through INFO of the final, output |
*> value of H. |
*> value of H. |
*> |
*> |
*> If INFO .GT. 0 and WANTT is .TRUE., then on exit |
*> If INFO > 0 and WANTT is .TRUE., then on exit |
*> |
*> |
*> (*) (initial value of H)*U = U*(final value of H) |
*> (*) (initial value of H)*U = U*(final value of H) |
*> |
*> |
Line 217
|
Line 217
|
*> value of H is upper Hessenberg and triangular in |
*> value of H is upper Hessenberg and triangular in |
*> rows and columns INFO+1 through IHI. |
*> rows and columns INFO+1 through IHI. |
*> |
*> |
*> If INFO .GT. 0 and WANTZ is .TRUE., then on exit |
*> If INFO > 0 and WANTZ is .TRUE., then on exit |
*> |
*> |
*> (final value of Z(ILO:IHI,ILOZ:IHIZ) |
*> (final value of Z(ILO:IHI,ILOZ:IHIZ) |
*> = (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U |
*> = (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U |
Line 225
|
Line 225
|
*> where U is the orthogonal matrix in (*) (regard- |
*> where U is the orthogonal matrix in (*) (regard- |
*> less of the value of WANTT.) |
*> less of the value of WANTT.) |
*> |
*> |
*> If INFO .GT. 0 and WANTZ is .FALSE., then Z is not |
*> If INFO > 0 and WANTZ is .FALSE., then Z is not |
*> accessed. |
*> accessed. |
*> \endverbatim |
*> \endverbatim |
* |
* |
* Authors: |
* Authors: |
* ======== |
* ======== |
* |
* |
*> \author Univ. of Tennessee |
*> \author Univ. of Tennessee |
*> \author Univ. of California Berkeley |
*> \author Univ. of California Berkeley |
*> \author Univ. of Colorado Denver |
*> \author Univ. of Colorado Denver |
*> \author NAG Ltd. |
*> \author NAG Ltd. |
* |
|
*> \date September 2012 |
|
* |
* |
*> \ingroup doubleOTHERauxiliary |
*> \ingroup doubleOTHERauxiliary |
* |
* |
Line 263
|
Line 261
|
SUBROUTINE DLAQR4( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, |
SUBROUTINE DLAQR4( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, |
$ ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO ) |
$ ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO ) |
* |
* |
* -- LAPACK auxiliary routine (version 3.4.2) -- |
* -- LAPACK auxiliary routine -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* September 2012 |
|
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, LWORK, N |
INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, LWORK, N |
Line 284
|
Line 281
|
* . DLAHQR because of insufficient subdiagonal scratch space. |
* . DLAHQR because of insufficient subdiagonal scratch space. |
* . (This is a hard limit.) ==== |
* . (This is a hard limit.) ==== |
INTEGER NTINY |
INTEGER NTINY |
PARAMETER ( NTINY = 11 ) |
PARAMETER ( NTINY = 15 ) |
* |
* |
* ==== Exceptional deflation windows: try to cure rare |
* ==== Exceptional deflation windows: try to cure rare |
* . slow convergence by varying the size of the |
* . slow convergence by varying the size of the |
Line 368
|
Line 365
|
END IF |
END IF |
* |
* |
* ==== NWR = recommended deflation window size. At this |
* ==== NWR = recommended deflation window size. At this |
* . point, N .GT. NTINY = 11, so there is enough |
* . point, N .GT. NTINY = 15, so there is enough |
* . subdiagonal workspace for NWR.GE.2 as required. |
* . subdiagonal workspace for NWR.GE.2 as required. |
* . (In fact, there is enough subdiagonal space for |
* . (In fact, there is enough subdiagonal space for |
* . NWR.GE.3.) ==== |
* . NWR.GE.4.) ==== |
* |
* |
NWR = ILAENV( 13, 'DLAQR4', JBCMPZ, N, ILO, IHI, LWORK ) |
NWR = ILAENV( 13, 'DLAQR4', JBCMPZ, N, ILO, IHI, LWORK ) |
NWR = MAX( 2, NWR ) |
NWR = MAX( 2, NWR ) |
NWR = MIN( IHI-ILO+1, ( N-1 ) / 3, NWR ) |
NWR = MIN( IHI-ILO+1, ( N-1 ) / 3, NWR ) |
* |
* |
* ==== NSR = recommended number of simultaneous shifts. |
* ==== NSR = recommended number of simultaneous shifts. |
* . At this point N .GT. NTINY = 11, so there is at |
* . At this point N .GT. NTINY = 15, so there is at |
* . enough subdiagonal workspace for NSR to be even |
* . enough subdiagonal workspace for NSR to be even |
* . and greater than or equal to two as required. ==== |
* . and greater than or equal to two as required. ==== |
* |
* |
NSR = ILAENV( 15, 'DLAQR4', JBCMPZ, N, ILO, IHI, LWORK ) |
NSR = ILAENV( 15, 'DLAQR4', JBCMPZ, N, ILO, IHI, LWORK ) |
NSR = MIN( NSR, ( N+6 ) / 9, IHI-ILO ) |
NSR = MIN( NSR, ( N-3 ) / 6, IHI-ILO ) |
NSR = MAX( 2, NSR-MOD( NSR, 2 ) ) |
NSR = MAX( 2, NSR-MOD( NSR, 2 ) ) |
* |
* |
* ==== Estimate optimal workspace ==== |
* ==== Estimate optimal workspace ==== |
Line 431
|
Line 428
|
* ==== NSMAX = the Largest number of simultaneous shifts |
* ==== NSMAX = the Largest number of simultaneous shifts |
* . for which there is sufficient workspace. ==== |
* . for which there is sufficient workspace. ==== |
* |
* |
NSMAX = MIN( ( N+6 ) / 9, 2*LWORK / 3 ) |
NSMAX = MIN( ( N-3 ) / 6, 2*LWORK / 3 ) |
NSMAX = NSMAX - MOD( NSMAX, 2 ) |
NSMAX = NSMAX - MOD( NSMAX, 2 ) |
* |
* |
* ==== NDFL: an iteration count restarted at deflation. ==== |
* ==== NDFL: an iteration count restarted at deflation. ==== |
Line 582
|
Line 579
|
* |
* |
* ==== Got NS/2 or fewer shifts? Use DLAHQR |
* ==== Got NS/2 or fewer shifts? Use DLAHQR |
* . on a trailing principal submatrix to |
* . on a trailing principal submatrix to |
* . get more. (Since NS.LE.NSMAX.LE.(N+6)/9, |
* . get more. (Since NS.LE.NSMAX.LE.(N-3)/6, |
* . there is enough space below the subdiagonal |
* . there is enough space below the subdiagonal |
* . to fit an NS-by-NS scratch array.) ==== |
* . to fit an NS-by-NS scratch array.) ==== |
* |
* |
Line 677
|
Line 674
|
END IF |
END IF |
END IF |
END IF |
* |
* |
* ==== Use up to NS of the the smallest magnatiude |
* ==== Use up to NS of the the smallest magnitude |
* . shifts. If there aren't NS shifts available, |
* . shifts. If there aren't NS shifts available, |
* . then use them all, possibly dropping one to |
* . then use them all, possibly dropping one to |
* . make the number of shifts even. ==== |
* . make the number of shifts even. ==== |
Line 697
|
Line 694
|
* . (NVE-by-KDU) vertical work WV arrow along |
* . (NVE-by-KDU) vertical work WV arrow along |
* . the left-hand-edge. ==== |
* . the left-hand-edge. ==== |
* |
* |
KDU = 3*NS - 3 |
KDU = 2*NS |
KU = N - KDU + 1 |
KU = N - KDU + 1 |
KWH = KDU + 1 |
KWH = KDU + 1 |
NHO = ( N-KDU+1-4 ) - ( KDU+1 ) + 1 |
NHO = ( N-KDU+1-4 ) - ( KDU+1 ) + 1 |