Routine: CLAHQR2()  File: SRC\clahqr2.f

 
 
# lines: 444
  # code: 444
  # comment: 0
  # blank:0
# Variables:27
# Callers:1
# Callings:1
# Words:139
# Keywords:65
 

 

..
     .. Local Scalars ..
     ..
     .. Local Arrays ..
     ..
     .. External Functions ..
     ..
     .. External Subroutines ..
     ..
     .. Intrinsic Functions ..
     ..
     .. Statement Functions ..
     ..
     .. Statement Function definitions ..
     ..
     .. Executable Statements ..
     Quick return if possible
     Set machine-dependent constants for the stopping criterion.
     If norm(H) <= sqrt(OVFL), overflow should not occur.
     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.
     ITN is the total number of QR iterations allowed.
     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.
     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.
        Look for a single small subdiagonal element.
           H(L,L-1) is negligible
        Exit from loop if a submatrix of order 1 or 2 has split off.
        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.
           Exceptional shift.
            S = ABS( REAL( H( I,I-1 ) ) ) + ABS( REAL( H( I-1,I-2 ) ) )
           Prepare to use Wilkinson's shift.
        Look for two consecutive small subdiagonal elements.
           Determine the effect of starting the double-shift QR
           iteration at row M, and see if this would make H(M,M-1)
           negligible.
        Double-shift QR step
           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

 
Display dynamic version Find AutoScroll Reload FontSize: - + Hide Comments Hide Blanks Frame FullScreen MailPrint

 
001        SUBROUTINE CLAHQR2( WANTT , WANTZ , N , ILO , IHI , H , LDH , W , ILOZ ,
002       $IHIZ , Z , LDZ , INFO )
003  
004  *     -- ScaLAPACK routine(version 1.7) --
005  *     Univ. of Tennessee , Univ. of California Berkeley , NAG Ltd. ,
006  *     Courant Institute , Argonne National Lab , and Rice University
007  *     June 22 , 2000
008  
009  *     .. Scalar Arguments ..
010        LOGICAL WANTT , WANTZ
011        INTEGER IHI , IHIZ , ILO , ILOZ , INFO , LDH , LDZ , N
012  *     ..
013  *     .. Array Arguments ..
014        COMPLEX H( LDH , * ) , W( * ) , Z( LDZ , * )
015  *     ..
016  
017  *     Purpose
018  *     === ====
019  
020  *     CLAHQR2 is an auxiliary routine called by CHSEQR to update the
021  *     eigenvalues and Schur decomposition already computed by CHSEQR , by
022  *     dealing with the Hessenberg submatrix in rows and columns ILO to IHI.
023  *     This version of CLAHQR(not the standard LAPACK version) uses a
024  *     double - shift algorithm(like LAPACK's SLAHQR).
025  *     Unlike the standard LAPACK convention , this does not assume the
026  *     subdiagonal is real , nor does it work to preserve this quality if
027  *     given.
028  
029  *     Arguments
030  *     === ======
031  
032  *     WANTT(input) LOGICAL
033  *     = .TRUE. : the full Schur form T is required ;
034  *     = .FALSE. : only eigenvalues are required.
035  
036  *     WANTZ(input) LOGICAL
037  *     = .TRUE. : the matrix of Schur vectors Z is required ;
038  *     = .FALSE. : Schur vectors are not required.
039  
040  *     N(input) INTEGER
041  *     The order of the matrix H. N >= 0.
042  
043  *     ILO(input) INTEGER
044  *     IHI(input) INTEGER
045  *     It is assumed that H is already upper triangular in rows and
046  *     columns IHI + 1 : N , and that H(ILO , ILO - 1) = 0(unless ILO = 1).
047  *     CLAHQR works primarily with the Hessenberg submatrix in rows
048  *     and columns ILO to IHI , but applies transformations to all of
049  *     H if WANTT is .TRUE..
050  *     1 <= ILO <= max(1 , IHI) ; IHI <= N.
051  
052  *     H(input / output) COMPLEX array , dimension(LDH , N)
053  *     On entry , the upper Hessenberg matrix H.
054  *     On exit , if WANTT is .TRUE. , H is upper triangular in rows
055  *     and columns ILO : IHI. If WANTT is .FALSE. , the contents of H
056  *     are unspecified on exit.
057  
058  *     LDH(input) INTEGER
059  *     The leading dimension of the array H. LDH >= max(1 , N).
060  
061  *     W(output) COMPLEX array , dimension(N)
062  *     The computed eigenvalues ILO to IHI are stored in the
063  *     corresponding elements of W. If WANTT is .TRUE. , the
064  *     eigenvalues are stored in the same order as on the diagonal
065  *     of the Schur form returned in H , with W(i) = H(i , i).
066  
067  *     ILOZ(input) INTEGER
068  *     IHIZ(input) INTEGER
069  *     Specify the rows of Z to which transformations must be
070  *     applied if WANTZ is .TRUE..
071  *     1 <= ILOZ <= ILO ; IHI <= IHIZ <= N.
072  
073  *     Z(input / output) COMPLEX array , dimension(LDZ , N)
074  *     If WANTZ is .TRUE. , on entry Z must contain the current
075  *     matrix Z of transformations , and on exit Z has been updated ;
076  *     transformations are applied only to the submatrix
077  *     Z(ILOZ : IHIZ , ILO : IHI). If WANTZ is .FALSE. , Z is not
078  *     referenced.
079  
080  *     LDZ(input) INTEGER
081  *     The leading dimension of the array Z. LDZ >= max(1 , N).
082  
083  *     INFO(output) INTEGER
084  *     = 0 : successful exit
085  *     > 0 : if INFO = i , CLAHQR failed to compute all the
086  *     eigenvalues ILO to IHI in a total of 30*(IHI - ILO + 1)
087  *     iterations ; elements i + 1 : ihi of W contain those
088  *     eigenvalues which have been successfully computed.
089  
090  *     Further Details
091  *     === ============
092  
093  *     Modified by Mark R. Fahey , June , 2000
094  
095  *     === ==================================================================
096  
097  *     .. Parameters ..
098        COMPLEX ZERO
099        PARAMETER( ZERO =( 0.0E + 0 , 0.0E + 0 ) )
100        REAL RZERO , RONE
101        PARAMETER( RZERO = 0.0E + 0 , RONE = 1.0E + 0 )
102        REAL DAT1 , DAT2
103        PARAMETER( DAT1 = 0.75E + 0 , DAT2 = - 0.4375E + 0 )
104        NR = MIN( 3 , I - K + 1 )
105        IF( K.GT.M )
106       $    CALL CCOPY( NR , H( K , K - 1 ) , 1 , V , 1 )
107            CALL CLARFG( NR , V( 1 ) , V( 2 ) , 1 , T1 )
108            IF( K.GT.M ) THEN
109                H( K , K - 1 ) = V( 1 )
110                H( K + 1 , K - 1 ) = ZERO
111                IF( K.LT.I - 1 )
112       $            H( K + 2 , K - 1 ) = ZERO
113                ELSE IF( M.GT.L ) THEN
114  *                 The real double - shift code uses H( K , K - 1 ) = - H( K , K - 1 )
115  *                 instead of the following.
116                    H( K , K - 1 ) = H( K , K - 1 ) - CONJG( T1 )*H( K , K - 1 )
117                END IF
118                V2 = V( 2 )
119                T2 = T1*V2
120                IF( NR.EQ.3 ) THEN
121                    V3 = V( 3 )
122                    T3 = T1*V3
123  
124  *                 Apply G from the left to transform the rows of the matrix
125  *                 in columns K to I2.
126  
127                    DO 60 J = K , I2
128                        SUM = CONJG( T1 )*H( K , J ) +
129       $                CONJG( T2 )*H( K + 1 , J ) +
130       $                CONJG( T3 )*H( K + 2 , J )
131                        H( K , J ) = H( K , J ) - SUM
132                        H( K + 1 , J ) = H( K + 1 , J ) - SUM*V2
133                        H( K + 2 , J ) = H( K + 2 , J ) - SUM*V3
134     60             CONTINUE
135  
136  *                 Apply G from the right to transform the columns of the
137  *                 matrix in rows I1 to min(K + 3 , I).
138  
139                    DO 70 J = I1 , MIN( K + 3 , I )
140                        SUM = T1*H( J , K ) + T2*H( J , K + 1 ) + T3*H( J , K + 2 )
141                        H( J , K ) = H( J , K ) - SUM
142                        H( J , K + 1 ) = H( J , K + 1 ) - SUM*CONJG( V2 )
143                        H( J , K + 2 ) = H( J , K + 2 ) - SUM*CONJG( V3 )
144     70             CONTINUE
145  
146                    IF( WANTZ ) THEN
147  
148  *                     Accumulate transformations in the matrix Z
149  
150                        DO 80 J = ILOZ , IHIZ
151                            SUM = T1*Z( J , K ) + T2*Z( J , K + 1 ) +
152       $                    T3*Z( J , K + 2 )
153                            Z( J , K ) = Z( J , K ) - SUM
154                            Z( J , K + 1 ) = Z( J , K + 1 ) - SUM*CONJG( V2 )
155                            Z( J , K + 2 ) = Z( J , K + 2 ) - SUM*CONJG( V3 )
156     80                 CONTINUE
157                    END IF
158                ELSE IF( NR.EQ.2 ) THEN
159  
160  *                 Apply G from the left to transform the rows of the matrix
161  *                 in columns K to I2.
162  
163                    DO 90 J = K , I2
164                        SUM = CONJG( T1 )*H( K , J ) +
165       $                CONJG( T2 )*H( K + 1 , J )
166                        H( K , J ) = H( K , J ) - SUM
167                        H( K + 1 , J ) = H( K + 1 , J ) - SUM*V2
168     90             CONTINUE
169  
170  *                 Apply G from the right to transform the columns of the
171  *                 matrix in rows I1 to min(K + 2 , I).
172  
173                    DO 100 J = I1 , MIN( K + 2 , I )
174                        SUM = T1*H( J , K ) + T2*H( J , K + 1 )
175                        H( J , K ) = H( J , K ) - SUM
176                        H( J , K + 1 ) = H( J , K + 1 ) - SUM*CONJG( V2 )
177    100             CONTINUE
178  
179                    IF( WANTZ ) THEN
180  
181  *                     Accumulate transformations in the matrix Z
182  
183                        DO 110 J = ILOZ , IHIZ
184                            SUM = T1*Z( J , K ) + T2*Z( J , K + 1 )
185                            Z( J , K ) = Z( J , K ) - SUM
186                            Z( J , K + 1 ) = Z( J , K + 1 ) - SUM*CONJG( V2 )
187    110                 CONTINUE
188                    END IF
189                END IF
190  
191  *             Since at the start of the QR step we have for M > L
192  *             H( K , K - 1 ) = H( K , K - 1 ) - CONJG( T1 )*H( K , K - 1 )
193  *             then we don't need to do the following
194  *             IF( K.EQ.M .AND. M.GT.L ) THEN
195  *             If the QR step was started at row M > L because two
196  *             consecutive small subdiagonals were found , then H(M , M - 1)
197  *             must also be updated by a factor of(1 - T1).
198  *             TEMP = ONE - T1
199  *             H( m , m - 1 ) = H( m , m - 1 )*CONJG( TEMP )
200  *             END IF
201    120 CONTINUE
202  
203  *     Ensure that H(I , I - 1) is real.
204  
205    130 CONTINUE
206  
207  *     Failure to converge in remaining number of iterations
208  
209        INFO = I
210        RETURN
211  
212    140 CONTINUE
213  
214        IF( L.EQ.I ) THEN
215  
216  *         H(I , I - 1) is negligible : one eigenvalue has converged.
217  
218            W( I ) = H( I , I )
219  
220        ELSE IF( L.EQ.I - 1 ) THEN
221  
222  *         H(I - 1 , I - 2) is negligible : a pair of eigenvalues have converged.
223  
224  *         Transform the 2 - by - 2 submatrix to standard Schur form ,
225  *         and compute and store the eigenvalues.
226  
227            CALL CLANV2 ( H( I - 1 , I - 1 ) , H( I - 1 , I ) , H( I , I - 1 ) ,
228       $    H( I , I ) , W( I - 1 ) , W( I ) , CS , SN )
229  
230            IF( WANTT ) THEN
231  
232  *             Apply the transformation to the rest of H.
233  
234                IF( I2.GT.I )
235       $            CALL CROT( I2 - I , H( I - 1 , I + 1 ) , LDH , H( I , I + 1 ) , LDH ,
236       $            CS , SN )
237                    CALL CROT( I - I1 - 1 , H( I1 , I - 1 ) , 1 , H( I1 , I ) , 1 , CS ,
238       $            CONJG( SN ) )
239                END IF
240                IF( WANTZ ) THEN
241  
242  *                 Apply the transformation to Z.
243  
244                    CALL CROT( NZ , Z( ILOZ , I - 1 ) , 1 , Z( ILOZ , I ) , 1 , CS ,
245       $            CONJG( SN ) )
246                END IF
247  
248            END IF
249  
250  *         Decrement number of remaining iterations , and return to start of
251  *         the main loop with new value of I.
252  
253            ITN = ITN - ITS
254            I = L - 1
255            GO TO 10
256  
257    150 CONTINUE
258        RETURN
259  
260  *     End of CLAHQR2
261  
262        END