Numérotation des lignes :

C DLAHQR    SOURCE    BP208322  18/07/10    21:15:08     9872           **> \brief \b DLAHQR 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 DLAHQR + dependencies*> &lt;a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlahqr.f">*> [TGZ]&lt;/a>*> &lt;a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlahqr.f">*> [ZIP]&lt;/a>*> &lt;a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlahqr.f">*> [TXT]&lt;/a>*> \endhtmlonly**  Definition:*  ===========**       SUBROUTINE DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,*                          ILOZ, IHIZ, Z, LDZ, INFO )**       .. Scalar Arguments ..*       INTEGER            IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N*       LOGICAL            WANTT, WANTZ*       ..*       .. Array Arguments ..*       DOUBLE PRECISION   H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * )*       ..***> \par Purpose:*  =============*>*> \verbatim*>*>    DLAHQR is an auxiliary routine called by DHSEQR to update the*>    eigenvalues and Schur decomposition already computed by DHSEQR, 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 quasi-triangular in*>          rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless*>          ILO = 1). DLAHQR 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 &lt;= ILO &lt;= max(1,IHI); IHI &lt;= N.*> \endverbatim*>*> \param[in,out] H*> \verbatim*>          H is DOUBLE PRECISION array, dimension (LDH,N)*>          On entry, the upper Hessenberg matrix H.*>          On exit, if INFO is zero and if WANTT is .TRUE., H is upper*>          quasi-triangular in rows and columns ILO:IHI, with any*>          2-by-2 diagonal blocks in standard form. If INFO is zero*>          and WANTT is .FALSE., the contents of H are unspecified on*>          exit.  The output state of H if INFO is nonzero is given*>          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] WR*> \verbatim*>          WR is DOUBLE PRECISION array, dimension (N)*> \endverbatim*>*> \param[out] WI*> \verbatim*>          WI is DOUBLE PRECISION array, dimension (N)*>          The real and imaginary parts, respectively, of the computed*>          eigenvalues ILO to IHI are stored in the corresponding*>          elements of WR and WI. If two eigenvalues are computed as a*>          complex conjugate pair, they are stored in consecutive*>          elements of WR and WI, say the i-th and (i+1)th, with*>          WI(i) > 0 and WI(i+1) &lt; 0. If WANTT is .TRUE., the*>          eigenvalues are stored in the same order as on the 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 block,*>          WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(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 &lt;= ILOZ &lt;= ILO; IHI &lt;= IHIZ &lt;= N.*> \endverbatim*>*> \param[in,out] Z*> \verbatim*>          Z is DOUBLE PRECISION array, dimension (LDZ,N)*>          If WANTZ is .TRUE., on entry Z must contain the current*>          matrix Z of transformations accumulated by DHSEQR, 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*>          .GT. 0: If INFO = i, DLAHQR failed to compute all the*>                  eigenvalues ILO to IHI in a total of 30 iterations*>                  per eigenvalue; elements i+1:ihi of WR and WI*>                  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.)*> \endverbatim**  Authors:*  ========**> \author Univ. of Tennessee*> \author Univ. of California Berkeley*> \author Univ. of Colorado Denver*> \author NAG Ltd.**> \date December 2016**> \ingroup doubleOTHERauxiliary**> \par Further Details:*  =====================*>*> \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 DLAHQR 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 DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,     $ILOZ, IHIZ, Z, LDZ, INFO )** -- LAPACK auxiliary routine (version 3.7.0) --* -- LAPACK is a software package provided by Univ. of Tennessee, --* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--* December 2016** .. Scalar Arguments .. INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N LOGICAL WANTT, WANTZ* ..* .. Array Arguments .. REAL*8 H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * )* ..** =========================================================** .. Parameters .. REAL*8 ZERO, ONE, TWO PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0d0 ) REAL*8 DAT1, DAT2 PARAMETER ( DAT1 = 0.75D+0, DAT2 = -0.4375D+0 )* ..* .. Local Scalars .. REAL*8 AA, AB, BA, BB, CS, DET, H11, H12, H21, H21S,$                   H22, RT1I, RT1R, RT2I, RT2R, RTDISC, S, SAFMAX,     $SAFMIN, SMLNUM, SN, SUM, T1, T2, T3, TR, TST,$                   ULP, V2, V3      INTEGER            I, I1, I2, ITS, ITMAX, J, K, L, M, NH, NR, NZ*     ..*     .. Local Arrays ..      REAL*8   V( 3 )*     ..*     .. External Functions ..      REAL*8   DLAMCH      EXTERNAL           DLAMCH*     ..*     .. External Subroutines ..      EXTERNAL           DCOPY, DLABAD, DLANV2, DLARFG, DROT*     ..**     .. Intrinsic Functions ..*      INTRINSIC          ABS, DBLE, MAX, MIN, SQRT**     ..**     .. Executable Statements ..*      INFO = 0**     Quick return if possible*      IF( N.EQ.0 )     $RETURN IF( ILO.EQ.IHI ) THEN WR( ILO ) = H( ILO, ILO ) WI( ILO ) = ZERO RETURN END IF** ==== clear out the trash ==== DO 10 J = ILO, IHI - 3 H( J+2, J ) = ZERO H( J+3, J ) = ZERO 10 CONTINUE IF( ILO.LE.IHI-2 )$   H( IHI, IHI-2 ) = ZERO*      NH = IHI - ILO + 1      NZ = IHIZ - ILOZ + 1**     Set machine-dependent constants for the stopping criterion.*      SAFMIN = DLAMCH( 'SAFE MINIMUM' )      SAFMAX = ONE / SAFMIN      CALL DLABAD( SAFMIN, SAFMAX )      ULP = DLAMCH( 'PRECISION' )      SMLNUM = SAFMIN*( DBLE( NH ) / ULP )**     I1 and I2 are the indices of the first row and last column of H*     to which transformations must be applied. If eigenvalues only are*     being computed, I1 and I2 are set inside the main loop.*      IF( WANTT ) THEN         I1 = 1         I2 = N      END IF**     ITMAX is the total number of QR iterations allowed.*      ITMAX = 30 * MAX( 10, NH )**     The main loop begins here. I is the loop index and decreases from*     IHI to ILO in steps of 1 or 2. Each iteration of the loop works*     with the active submatrix in rows and columns L to I.*     Eigenvalues I+1 to IHI have already converged. Either L = ILO or*     H(L,L-1) is negligible so that the matrix splits.*      I = IHI   20 CONTINUE      L = ILO      IF( I.LT.ILO )     $GO TO 160** Perform QR iterations on rows and columns ILO to I until a* submatrix of order 1 or 2 splits off at the bottom because a* subdiagonal element has become negligible.* DO 140 ITS = 0, ITMAX** Look for a single small subdiagonal element.* DO 30 K = I, L + 1, -1 IF( ABS( H( K, K-1 ) ).LE.SMLNUM )$         GO TO 40            TST = ABS( H( K-1, K-1 ) ) + ABS( H( K, K ) )            IF( TST.EQ.ZERO ) THEN               IF( K-2.GE.ILO )     $TST = TST + ABS( H( K-1, K-2 ) ) IF( K+1.LE.IHI )$            TST = TST + ABS( H( K+1, K ) )            END IF*           ==== The following is a conservative small subdiagonal*           .    deflation  criterion due to Ahues & Tisseur (LAWN 122,*           .    1997). It has better mathematical foundation and*           .    improves accuracy in some cases.  ====            IF( ABS( H( K, K-1 ) ).LE.ULP*TST ) THEN               AB = MAX( ABS( H( K, K-1 ) ), ABS( H( K-1, K ) ) )               BA = MIN( ABS( H( K, K-1 ) ), ABS( H( K-1, K ) ) )               AA = MAX( ABS( H( K, K ) ),     $ABS( H( K-1, K-1 )-H( K, K ) ) ) BB = MIN( ABS( H( K, K ) ),$              ABS( H( K-1, K-1 )-H( K, K ) ) )               S = AA + AB               IF( BA*( AB / S ).LE.MAX( SMLNUM,     $ULP*( BB*( AA / S ) ) ) )GO TO 40 END IF 30 CONTINUE 40 CONTINUE L = K IF( L.GT.ILO ) THEN** H(L,L-1) is negligible* H( L, L-1 ) = ZERO END IF** Exit from loop if a submatrix of order 1 or 2 has split off.* IF( L.GE.I-1 )$      GO TO 150**        Now the active submatrix is in rows and columns L to I. If*        eigenvalues only are being computed, only the active submatrix*        need be transformed.*         IF( .NOT.WANTT ) THEN            I1 = L            I2 = I         END IF*         IF( ITS.EQ.10 ) THEN**           Exceptional shift.*            S = ABS( H( L+1, L ) ) + ABS( H( L+2, L+1 ) )            H11 = DAT1*S + H( L, L )            H12 = DAT2*S            H21 = S            H22 = H11         ELSE IF( ITS.EQ.20 ) THEN**           Exceptional shift.*            S = ABS( H( I, I-1 ) ) + ABS( H( I-1, I-2 ) )            H11 = DAT1*S + H( I, I )            H12 = DAT2*S            H21 = S            H22 = H11         ELSE**           Prepare to use Francis' double shift*           (i.e. 2nd degree generalized Rayleigh quotient)*            H11 = H( I-1, I-1 )            H21 = H( I, I-1 )            H12 = H( I-1, I )            H22 = H( I, I )         END IF         S = ABS( H11 ) + ABS( H12 ) + ABS( H21 ) + ABS( H22 )         IF( S.EQ.ZERO ) THEN            RT1R = ZERO            RT1I = ZERO            RT2R = ZERO            RT2I = ZERO         ELSE            H11 = H11 / S            H21 = H21 / S            H12 = H12 / S            H22 = H22 / S            TR = ( H11+H22 ) / TWO            DET = ( H11-TR )*( H22-TR ) - H12*H21            RTDISC = SQRT( ABS( DET ) )            IF( DET.GE.ZERO ) THEN**              ==== complex conjugate shifts ====*               RT1R = TR*S               RT2R = RT1R               RT1I = RTDISC*S               RT2I = -RT1I            ELSE**              ==== real shifts (use only one of them)  ====*               RT1R = TR + RTDISC               RT2R = TR - RTDISC               IF( ABS( RT1R-H22 ).LE.ABS( RT2R-H22 ) ) THEN                  RT1R = RT1R*S                  RT2R = RT1R               ELSE                  RT2R = RT2R*S                  RT1R = RT2R               END IF               RT1I = ZERO               RT2I = ZERO            END IF         END IF**        Look for two consecutive small subdiagonal elements.*         DO 50 M = I - 2, L, -1*           Determine the effect of starting the double-shift QR*           iteration at row M, and see if this would make H(M,M-1)*           negligible.  (The following uses scaling to avoid*           overflows and most underflows.)*            H21S = H( M+1, M )            S = ABS( H( M, M )-RT2R ) + ABS( RT2I ) + ABS( H21S )            H21S = H( M+1, M ) / S            V( 1 ) = H21S*H( M, M+1 ) + ( H( M, M )-RT1R )*     $( ( H( M, M )-RT2R ) / S ) - RT1I*( RT2I / S ) V( 2 ) = H21S*( H( M, M )+H( M+1, M+1 )-RT1R-RT2R ) V( 3 ) = H21S*H( M+2, M+1 ) S = ABS( V( 1 ) ) + ABS( V( 2 ) ) + ABS( V( 3 ) ) V( 1 ) = V( 1 ) / S V( 2 ) = V( 2 ) / S V( 3 ) = V( 3 ) / S IF( M.EQ.L )$         GO TO 60            IF( ABS( H( M, M-1 ) )*( ABS( V( 2 ) )+ABS( V( 3 ) ) ).LE.     $ULP*ABS( V( 1 ) )*( ABS( H( M-1, M-1 ) )+ABS( H( M,$          M ) )+ABS( H( M+1, M+1 ) ) ) )GO TO 60   50    CONTINUE   60    CONTINUE**        Double-shift QR step*         DO 130 K = M, I - 1**           The first iteration of this loop determines a reflection G*           from the vector V and applies it from left and right to H,*           thus creating a nonzero bulge below the subdiagonal.**           Each subsequent iteration determines a reflection G to*           restore the Hessenberg form in the (K-1)th column, and thus*           chases the bulge one step toward the bottom of the active*           submatrix. NR is the order of G.*            NR = MIN( 3, I-K+1 )            IF( K.GT.M )     $CALL DCOPY( NR, H( K, K-1 ), 1, V, 1 ) CALL DLARFG( NR, V( 1 ), V( 2 ), 1, T1 ) IF( K.GT.M ) THEN H( K, K-1 ) = V( 1 ) H( K+1, K-1 ) = ZERO IF( K.LT.I-1 )$            H( K+2, K-1 ) = ZERO            ELSE IF( M.GT.L ) THEN*               ==== Use the following instead of*               .    H( K, K-1 ) = -H( K, K-1 ) to*               .    avoid a bug when v(2) and v(3)*               .    underflow. ====               H( K, K-1 ) = H( K, K-1 )*( ONE-T1 )            END IF            V2 = V( 2 )            T2 = T1*V2            IF( NR.EQ.3 ) THEN               V3 = V( 3 )               T3 = T1*V3**              Apply G from the left to transform the rows of the matrix*              in columns K to I2.*               DO 70 J = K, I2                  SUM = H( K, J ) + V2*H( K+1, J ) + V3*H( K+2, J )                  H( K, J ) = H( K, J ) - SUM*T1                  H( K+1, J ) = H( K+1, J ) - SUM*T2                  H( K+2, J ) = H( K+2, J ) - SUM*T3   70          CONTINUE**              Apply G from the right to transform the columns of the*              matrix in rows I1 to min(K+3,I).*               DO 80 J = I1, MIN( K+3, I )                  SUM = H( J, K ) + V2*H( J, K+1 ) + V3*H( J, K+2 )                  H( J, K ) = H( J, K ) - SUM*T1                  H( J, K+1 ) = H( J, K+1 ) - SUM*T2                  H( J, K+2 ) = H( J, K+2 ) - SUM*T3   80          CONTINUE*               IF( WANTZ ) THEN**                 Accumulate transformations in the matrix Z*                  DO 90 J = ILOZ, IHIZ                     SUM = Z( J, K ) + V2*Z( J, K+1 ) + V3*Z( J, K+2 )                     Z( J, K ) = Z( J, K ) - SUM*T1                     Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2                     Z( J, K+2 ) = Z( J, K+2 ) - SUM*T3   90             CONTINUE               END IF            ELSE IF( NR.EQ.2 ) THEN**              Apply G from the left to transform the rows of the matrix*              in columns K to I2.*               DO 100 J = K, I2                  SUM = H( K, J ) + V2*H( K+1, J )                  H( K, J ) = H( K, J ) - SUM*T1                  H( K+1, J ) = H( K+1, J ) - SUM*T2  100          CONTINUE**              Apply G from the right to transform the columns of the*              matrix in rows I1 to min(K+3,I).*               DO 110 J = I1, I                  SUM = H( J, K ) + V2*H( J, K+1 )                  H( J, K ) = H( J, K ) - SUM*T1                  H( J, K+1 ) = H( J, K+1 ) - SUM*T2  110          CONTINUE*               IF( WANTZ ) THEN**                 Accumulate transformations in the matrix Z*                  DO 120 J = ILOZ, IHIZ                     SUM = Z( J, K ) + V2*Z( J, K+1 )                     Z( J, K ) = Z( J, K ) - SUM*T1                     Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2  120             CONTINUE               END IF            END IF  130    CONTINUE*  140 CONTINUE**     Failure to converge in remaining number of iterations*      INFO = I      RETURN*  150 CONTINUE*      IF( L.EQ.I ) THEN**        H(I,I-1) is negligible: one eigenvalue has converged.*         WR( I ) = H( I, I )         WI( I ) = ZERO      ELSE IF( L.EQ.I-1 ) THEN**        H(I-1,I-2) is negligible: a pair of eigenvalues have converged.**        Transform the 2-by-2 submatrix to standard Schur form,*        and compute and store the eigenvalues.*         CALL DLANV2( H( I-1, I-1 ), H( I-1, I ), H( I, I-1 ),     $H( I, I ), WR( I-1 ), WI( I-1 ), WR( I ), WI( I ),$                CS, SN )*         IF( WANTT ) THEN**           Apply the transformation to the rest of H.*            IF( I2.GT.I )     $CALL DROT( I2-I, H( I-1, I+1 ), LDH, H( I, I+1 ), LDH,$                    CS, SN )            CALL DROT( I-I1-1, H( I1, I-1 ), 1, H( I1, I ), 1, CS, SN )         END IF         IF( WANTZ ) THEN**           Apply the transformation to Z.*            CALL DROT( NZ, Z( ILOZ, I-1 ), 1, Z( ILOZ, I ), 1, CS, SN )         END IF      END IF**     return to start of the main loop with new value of I.*      I = L - 1      GO TO 20*  160 CONTINUE      RETURN**     End of DLAHQR*      END

© Cast3M 2003 - Tous droits réservés.
Mentions légales