Routine: ZSTEQR2()  File: SRC\zsteqr2.f

 
 
# lines: 662
  # code: 662
  # comment: 0
  # blank:0
# Variables:66
# Callers:1
# Callings:0
# Words:282
# Keywords:157
 

 

Skip the current step: the subdiagonal info is just noise.
        Lookahead over
        Inner loop
           If eigenvectors are desired, then save rotations.
        If eigenvectors are desired, then apply saved rotations.
        Eigenvalue found.
        QR Iteration
        Look for small superdiagonal element.
        If remaining matrix is 2-by-2, use DLAE2 or DLAEV2
        to compute its eigensystem.
        Form shift.
           Do not do a lookahead!
           This is the lookahead loop, going until we have
           convergence or too many steps have been taken.
        Inner loop
           Find GP & RP for the next iteration
           Goto put in by G. Henry to fix ALPHA problem
              GP = ( ( OLDGP+P )-( D( L )-P ) ) /
    $              ( TWO*( C*OLDRP-B )+SAFMIN )
           Make sure that we are making progress
           Skip the current step: the subdiagonal info is just noise.
        Lookahead over
           If eigenvectors are desired, then save rotations.
        If eigenvectors are desired, then apply saved rotations.
        Eigenvalue found.

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

 
001        SUBROUTINE ZSTEQR2( COMPZ , N , D , E , Z , LDZ , NR , WORK , INFO )
002  
003  *     -- ScaLAPACK routine(version 1.7) --
004  *     Univ. of Tennessee , Univ. of California Berkeley , NAG Ltd. ,
005  *     Courant Institute , Argonne National Lab , and Rice University
006  *     November 15 , 1997
007  
008  *     .. Scalar Arguments ..
009        CHARACTER COMPZ
010        INTEGER INFO , LDZ , N , NR
011  *     ..
012  *     .. Array Arguments ..
013        DOUBLE PRECISION D( * ) , E( * ) , WORK( * )
014        COMPLEX*16 Z( LDZ , * )
015  *     ..
016  
017  *     Purpose
018  *     === ====
019  
020  *     ZSTEQR2 is a modified version of LAPACK routine ZSTEQR.
021  *     ZSTEQR2 computes all eigenvalues and , optionally , eigenvectors of a
022  *     symmetric tridiagonal matrix using the implicit QL or QR method.
023  *     ZSTEQR2 is modified from ZSTEQR to allow each ScaLAPACK process
024  *     running ZSTEQR2 to perform updates on a distributed matrix Q.
025  *     Proper usage of ZSTEQR2 can be gleaned from
026  *     examination of ScaLAPACK's * PZHEEV.
027  *     ZSTEQR2 incorporates changes attributed to Greg Henry.
028  
029  *     Arguments
030  *     === ======
031  
032  *     COMPZ(input) CHARACTER*1
033  *     = 'N' : Compute eigenvalues only.
034  *     = 'I' : Compute eigenvalues and eigenvectors of the
035  *     tridiagonal matrix. Z must be initialized to the
036  *     identity matrix by PZLASET or ZLASET prior
037  *     to entering this subroutine.
038  
039  *     N(input) INTEGER
040  *     The order of the matrix. N >= 0.
041  
042  *     D(input / output) DOUBLE PRECISION array , dimension(N)
043  *     On entry , the diagonal elements of the tridiagonal matrix.
044  *     On exit , if INFO = 0 , the eigenvalues in ascending order.
045  
046  *     E(input / output) DOUBLE PRECISION array , dimension(N - 1)
047  *     On entry , the(n - 1) subdiagonal elements of the tridiagonal
048  *     matrix.
049  *     On exit , E has been destroyed.
050  
051  *     Z(local input / local output) COMPLEX*16 array , global
052  *     dimension(N , N) , local dimension(LDZ , NR).
053  *     On entry , if COMPZ = 'V' , then Z contains the orthogonal
054  *     matrix used in the reduction to tridiagonal form.
055  *     On exit , if INFO = 0 , then if COMPZ = 'V' , Z contains the
056  *     orthonormal eigenvectors of the original symmetric matrix ,
057  *     and if COMPZ = 'I' , Z contains the orthonormal eigenvectors
058  *     of the symmetric tridiagonal matrix.
059  *     If COMPZ = 'N' , then Z is not referenced.
060  
061  *     LDZ(input) INTEGER
062  *     The leading dimension of the array Z. LDZ >= 1 , and if
063  *     eigenvectors are desired , then LDZ >= max(1 , N).
064  
065  *     NR(input) INTEGER
066  *     NR = MAX(1 , NUMROC( N , NB , MYPROW , 0 , NPROCS ) ).
067  *     If COMPZ = 'N' , then NR is not referenced.
068  
069  *     WORK(workspace) DOUBLE PRECISION array , dimension(max(1 , 2*N - 2))
070  *     If COMPZ = 'N' , then WORK is not referenced.
071  
072  *     INFO(output) INTEGER
073  *     = 0 : successful exit
074  *     < 0 : if INFO = - i , the i - th argument had an illegal value
075  *     > 0 : the algorithm has failed to find all the eigenvalues in
076  *     a total of 30*N iterations ; if INFO = i , then i
077  *     elements of E have not converged to zero ; on exit , D
078  *     and E contain the elements of a symmetric tridiagonal
079  *     matrix which is orthogonally similar to the original
080  *     matrix.
081  
082  *     === ==================================================================
083  
084  *     .. Parameters ..
085        DOUBLE PRECISION ZERO , ONE , TWO , THREE , HALF
086        PARAMETER( ZERO = 0.0D0 , ONE = 1.0D0 , TWO = 2.0D0 ,
087       $THREE = 3.0D0 , HALF = 0.5D0 )
088        COMPLEX*16 CONE
089        PARAMETER( CONE =( 1.0D0 , 1.0D0 ) )
090        INTEGER MAXIT , NMAXLOOK
091        PARAMETER( MAXIT = 30 , NMAXLOOK = 15 )
092  *     ..
093  *     .. Local Scalars ..
094        INTEGER I , ICOMPZ , II , ILAST , ISCALE , J , JTOT , K , L ,
095       $L1 , LEND , LENDM1 , LENDP1 , LENDSV , LM1 , LSV , M ,
096       $MM , MM1 , NLOOK , NM1 , NMAXIT
097        DOUBLE PRECISION ANORM , B , C , EPS , EPS2 , F , G , GP , OLDEL , OLDGP ,
098       $OLDRP , P , R , RP , RT1 , RT2 , S , SAFMAX , SAFMIN ,
099       $SSFMAX , SSFMIN , TST , TST1
100  *     ..
101  *     .. External Functions ..
102        LOGICAL LSAME
103        DOUBLE PRECISION DLAMCH , DLANST , DLAPY2
104        EXTERNAL LSAME , DLAMCH , DLANST , DLAPY2
105  *     ..
106  *     .. External Subroutines ..
107        EXTERNAL DLAEV2 , DLARTG , DLASCL , DSTERF , XERBLA , ZLASR ,
108       $ZSWAP
109  *     ..
110  *     .. Intrinsic Functions ..
111        INTRINSIC ABS , MAX , SIGN , SQRT
112  *     ..
113  *     .. Executable Statements ..
114  
115  *     Test the input parameters.
116  
117        ILAST = 0
118        INFO = 0
119  
120        IF( LSAME( COMPZ , 'N' ) ) THEN
121            ICOMPZ = 0
122        ELSEIF( LSAME( COMPZ , 'I' ) ) THEN
123        ICOMPZ = 1
124        ELSE
125            ICOMPZ = - 1
126        ENDIF
127        IF( ICOMPZ.LT.0 ) THEN
128            INFO = - 1
129        ELSEIF( N.LT.0 ) THEN
130        INFO = - 2
131        ELSEIF( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1 , NR ) ) THEN
132        INFO = - 6
133        ENDIF
134        IF( INFO.NE.0 ) THEN
135            CALL XERBLA( 'ZSTEQR2' , - INFO )
136            RETURN
137        ENDIF
138  
139  *     Quick return if possible
140  
141        IF( N.EQ.0 )
142       $    RETURN
143  
144  *         If eigenvectors aren't not desired , this is faster
145  
146            IF( ICOMPZ.EQ.0 ) THEN
147                CALL DSTERF( N , D , E , INFO )
148                RETURN
149            ENDIF
150  
151            IF( N.EQ.1 ) THEN
152                Z( 1 , 1 ) = CONE
153                RETURN
154            ENDIF
155  
156  *         Determine the unit roundoff and over / underflow thresholds.
157  
158            EPS = DLAMCH( 'E' )
159            EPS2 = EPS**2
160            SAFMIN = DLAMCH( 'S' )
161            SAFMAX = ONE / SAFMIN
162            SSFMAX = SQRT( SAFMAX ) / THREE
163            SSFMIN = SQRT( SAFMIN ) / EPS2
164  
165  *         Compute the eigenvalues and eigenvectors of the tridiagonal
166  *         matrix.
167  
168            NMAXIT = N*MAXIT
169            JTOT = 0
170  
171  *         Determine where the matrix splits and choose QL or QR iteration
172  *         for each block , according to whether top or bottom diagonal
173  *         element is smaller.
174  
175            L1 = 1
176            NM1 = N - 1
177  
178     10 CONTINUE
179        IF( L1.GT.N )
180       $    GOTO 220
181            IF( L1.GT.1 )
182       $        E( L1 - 1 ) = ZERO
183                IF( L1.LE.NM1 ) THEN
184                    DO 20 M = L1 , NM1
185                        TST = ABS( E( M ) )
186                        IF( TST.EQ.ZERO )
187       $                    GOTO 30
188                            IF( TST.LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M +
189       $                    1 ) ) ) )*EPS ) THEN
190                            E( M ) = ZERO
191                            GOTO 30
192                        ENDIF
193     20             CONTINUE
194                ENDIF
195                M = N
196  
197     30 CONTINUE
198        L = L1
199        LSV = L
200        LEND = M
201        LENDSV = LEND
202        L1 = M + 1
203        IF( LEND.EQ.L )
204       $    GOTO 10
205  
206  *         Scale submatrix in rows and columns L to LEND
207  
208            ANORM = DLANST( 'I' , LEND - L + 1 , D( L ) , E( L ) )
209            ISCALE = 0
210            IF( ANORM.EQ.ZERO )
211       $        GOTO 10
212                IF( ANORM.GT.SSFMAX ) THEN
213                    ISCALE = 1
214                    CALL DLASCL( 'G' , 0 , 0 , ANORM , SSFMAX , LEND - L + 1 , 1 , D( L ) , N ,
215       $            INFO )
216                    CALL DLASCL( 'G' , 0 , 0 , ANORM , SSFMAX , LEND - L , 1 , E( L ) , N ,
217       $            INFO )
218                ELSEIF( ANORM.LT.SSFMIN ) THEN
219                ISCALE = 2
220                CALL DLASCL( 'G' , 0 , 0 , ANORM , SSFMIN , LEND - L + 1 , 1 , D( L ) , N ,
221       $        INFO )
222                CALL DLASCL( 'G' , 0 , 0 , ANORM , SSFMIN , LEND - L , 1 , E( L ) , N ,
223       $        INFO )
224            ENDIF
225  
226  *         Choose between QL and QR iteration
227  
228            IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN
229                LEND = LSV
230                L = LENDSV
231            ENDIF
232  
233            IF( LEND.GT.L ) THEN
234  
235  *             QL Iteration
236  
237  *             Look for small subdiagonal element.
238  
239     40 CONTINUE
240        IF( L.NE.LEND ) THEN
241            LENDM1 = LEND - 1
242            DO 50 M = L , LENDM1
243                TST = ABS( E( M ) )**2
244                IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M + 1 ) ) +
245       $        SAFMIN )GOTO 60
246     50     CONTINUE
247        ENDIF
248  
249        M = LEND
250  
251     60 CONTINUE
252        IF( M.LT.LEND )
253       $    E( M ) = ZERO
254            P = D( L )
255            IF( M.EQ.L )
256       $        GOTO 110
257  
258  *             If remaining matrix is 2 - by - 2 , use DLAE2 or DLAEV2
259  *             to compute its eigensystem.
260  
261                IF( M.EQ.L + 1 ) THEN
262                    CALL DLAEV2( D( L ) , E( L ) , D( L + 1 ) , RT1 , RT2 , C , S )
263                    WORK( L ) = C
264                    WORK( N - 1 + L ) = S
265                    CALL ZLASR( 'R' , 'V' , 'B' , NR , 2 , WORK( L ) , WORK( N - 1 + L ) ,
266       $            Z( 1 , L ) , LDZ )
267                    D( L ) = RT1
268                    D( L + 1 ) = RT2
269                    E( L ) = ZERO
270                    L = L + 2
271                    IF( L.LE.LEND )
272       $                GOTO 40
273                        GOTO 200
274                    ENDIF
275  
276                    IF( JTOT.EQ.NMAXIT )
277       $                GOTO 200
278                        JTOT = JTOT + 1
279  
280  *                     Form shift.
281  
282                        G =( D( L + 1 ) - P ) / ( TWO*E( L ) )
283                        R = DLAPY2( G , ONE )
284                        G = D( M ) - P + ( E( L ) / ( G + SIGN( R , G ) ) )
285  
286                        IF( ICOMPZ.EQ.0 ) THEN
287  *                         Do not do a lookahead !
288                            GOTO 90
289                        ENDIF
290  
291                        OLDEL = ABS( E( L ) )
292                        GP = G
293                        RP = R
294                        TST = ABS( E( L ) )**2
295                        TST = TST / (( EPS2*ABS( D( L ) ) )*ABS( D( L + 1 ) ) + SAFMIN )
296  
297                        NLOOK = 1
298                        IF(( TST.GT.ONE ) .AND.( NLOOK.LE.NMAXLOOK ) ) THEN
299     70 CONTINUE
300  
301  *     This is the lookahead loop , going until we have
302  *     convergence or too many steps have been taken.
303  
304        S = ONE
305        C = ONE
306        P = ZERO
307        MM1 = M - 1
308        DO 80 I = MM1 , L , - 1
309            F = S*E( I )
310            B = C*E( I )
311            CALL DLARTG( GP , F , C , S , RP )
312            GP = D( I + 1 ) - P
313            RP =( D( I ) - GP )*S + TWO*C*B
314            P = S*RP
315            IF( I.NE.L )
316       $        GP = C*RP - B
317     80 CONTINUE
318        OLDGP = GP
319        OLDRP = RP
320  *     Find GP & RP for the next iteration
321        IF( ABS( C*OLDRP - B ).GT.SAFMIN ) THEN
322            GP =(( OLDGP + P ) - ( D( L ) - P ) ) / ( TWO*( C*OLDRP - B ) )
323        ELSE
324  
325  *         Goto put in by G. Henry to fix ALPHA problem
326  
327            GOTO 90
328  *         GP =(( OLDGP + P ) - ( D( L ) - P ) ) /
329  *         $( TWO*( C*OLDRP - B ) + SAFMIN )
330        ENDIF
331        RP = DLAPY2( GP , ONE )
332        GP = D( M ) - ( D( L ) - P ) +
333       $(( C*OLDRP - B ) / ( GP + SIGN( RP , GP ) ) )
334        TST1 = TST
335        TST = ABS( C*OLDRP - B )**2
336        TST = TST / (( EPS2*ABS( D( L ) - P ) )*ABS( OLDGP + P ) +
337       $SAFMIN )
338  *     Make sure that we are making progress
339        IF( ABS( C*OLDRP - B ).GT.0.9D0*OLDEL ) THEN
340            IF( ABS( C*OLDRP - B ).GT.OLDEL ) THEN
341                GP = G
342                RP = R
343            ENDIF
344            TST = HALF
345        ELSE
346            OLDEL = ABS( C*OLDRP - B )
347        ENDIF
348        NLOOK = NLOOK + 1
349        IF(( TST.GT.ONE ) .AND.( NLOOK.LE.NMAXLOOK ) )
350       $    GOTO 70
351        ENDIF
352  
353        IF(( TST.LE.ONE ) .AND.( TST.NE.HALF ) .AND.
354       $( ABS( P ).LT.EPS*ABS( D( L ) ) ) .AND.
355       $( ILAST.EQ.L ) .AND.( ABS( E( L ) )**2.LE.10000.0D0*
356       $(( EPS2*ABS( D( L ) ) )*ABS( D( L + 1 ) ) + SAFMIN ) ) ) THEN
357    190 CONTINUE
358        D( L ) = P
359  
360        L = L - 1
361        IF( L.GE.LEND )
362       $    GOTO 120
363            GOTO 200
364  
365        ENDIF
366  
367  *     Undo scaling if necessary
368  
369    200 CONTINUE
370        IF( ISCALE.EQ.1 ) THEN
371            CALL DLASCL( 'G' , 0 , 0 , SSFMAX , ANORM , LENDSV - LSV + 1 , 1 ,
372       $    D( LSV ) , N , INFO )
373            CALL DLASCL( 'G' , 0 , 0 , SSFMAX , ANORM , LENDSV - LSV , 1 , E( LSV ) ,
374       $    N , INFO )
375        ELSEIF( ISCALE.EQ.2 ) THEN
376        CALL DLASCL( 'G' , 0 , 0 , SSFMIN , ANORM , LENDSV - LSV + 1 , 1 ,
377       $D( LSV ) , N , INFO )
378        CALL DLASCL( 'G' , 0 , 0 , SSFMIN , ANORM , LENDSV - LSV , 1 , E( LSV ) ,
379       $N , INFO )
380        ENDIF
381  
382  *     Check for no convergence to an eigenvalue after a total
383  *     of N*MAXIT iterations.
384  
385        IF( JTOT.LT.NMAXIT )
386       $    GOTO 10
387            DO 210 I = 1 , N - 1
388                IF( E( I ).NE.ZERO )
389       $            INFO = INFO + 1
390    210     CONTINUE
391            GOTO 250
392  
393  *         Order eigenvalues and eigenvectors.
394  
395    220 CONTINUE
396  
397  *     Use Selection Sort to minimize swaps of eigenvectors
398  
399        DO 240 II = 2 , N
400            I = II - 1
401            K = I
402            P = D( I )
403            DO 230 J = II , N
404                IF( D( J ).LT.P ) THEN
405                    K = J
406                    P = D( J )
407                ENDIF
408    230     CONTINUE
409            IF( K.NE.I ) THEN
410                D( K ) = D( I )
411                D( I ) = P
412                CALL ZSWAP( NR , Z( 1 , I ) , 1 , Z( 1 , K ) , 1 )
413            ENDIF
414    240 CONTINUE
415  
416    250 CONTINUE
417  *     WRITE( * , FMT = * )'JTOT' , JTOT
418        RETURN
419  
420  *     End of DSTEQR2
421  
422        END