Routine DSTEQR2() Documented By Universal Report

Routine: DSTEQR2()

# lines:	499
# code:	499
# comment:	0
# blank:	0
# Variables:	32
# Callers:	1
# Callings:	0
# Words:	189
# Keywords:	99

..
     .. Local Scalars ..
     ..
     .. External Functions ..
     ..
     .. External Subroutines ..
     ..
     .. Intrinsic Functions ..
     ..
     .. Executable Statements ..
     Test the input parameters.
     Quick return if possible
     Determine the unit roundoff and over/underflow thresholds.
     Compute the eigenvalues and eigenvectors of the tridiagonal
     matrix.
     Determine where the matrix splits and choose QL or QR iteration
     for each block, according to whether top or bottom diagonal
     element is smaller.
     Scale submatrix in rows and columns L to LEND
     Choose between QL and QR iteration
        QL Iteration
        Look for small subdiagonal element.
        If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
        to compute its eigensystem.
        Form shift.
        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.

Display dynamic version

FontSize: - +

001        SUBROUTINE DSTEQR2( COMPZ , N , D , E , Z , LDZ , NR , WORK , INFO )
002  
003  *     -- LAPACK routine(version 2.0) --
004  *     Univ. of Tennessee , Univ. of California Berkeley , NAG Ltd. ,
005  *     Courant Institute , Argonne National Lab , and Rice University
006  *     September 30 , 1994
007  
008  *     .. Scalar Arguments ..
009        CHARACTER COMPZ
010        INTEGER INFO , LDZ , N , NR
011  *     ..
012  *     .. Array Arguments ..
013        DOUBLE PRECISION D( * ) , E( * ) , WORK( * ) , Z( LDZ , * )
014  *     ..
015  
016  *     Purpose
017  *     === ====
018  
019  *     DSTEQR2 is a modified version of LAPACK routine DSTEQR.
020  *     DSTEQR2 computes all eigenvalues and , optionally , eigenvectors of a
021  *     symmetric tridiagonal matrix using the implicit QL or QR method.
022  *     DSTEQR2 is modified from DSTEQR to allow each ScaLAPACK process
023  *     running DSTEQR2 to perform updates on a distributed matrix Q.
024  *     Proper usage of DSTEQR2 can be gleaned from examination of ScaLAPACK's
025  *     PDSYEV.
026  
027  *     Arguments
028  *     === ======
029  
030  *     COMPZ(input) CHARACTER*1
031  *     = 'N' : Compute eigenvalues only.
032  *     = 'I' : Compute eigenvalues and eigenvectors of the
033  *     tridiagonal matrix. Z must be initialized to the
034  *     identity matrix by PDLASET or DLASET prior to entering
035  *     this subroutine.
036  
037  *     N(input) INTEGER
038  *     The order of the matrix. N >= 0.
039  
040  *     D(input / output) DOUBLE PRECISION array , dimension(N)
041  *     On entry , the diagonal elements of the tridiagonal matrix.
042  *     On exit , if INFO = 0 , the eigenvalues in ascending order.
043  
044  *     E(input / output) DOUBLE PRECISION array , dimension(N - 1)
045  *     On entry , the(n - 1) subdiagonal elements of the tridiagonal
046  *     matrix.
047  *     On exit , E has been destroyed.
048  
049  *     Z(local input / local output) DOUBLE PRECISION array , global
050  *     dimension(N , N) , local dimension(LDZ , NR).
051  *     On entry , if COMPZ = 'V' , then Z contains the orthogonal
052  *     matrix used in the reduction to tridiagonal form.
053  *     On exit , if INFO = 0 , then if COMPZ = 'V' , Z contains the
054  *     orthonormal eigenvectors of the original symmetric matrix ,
055  *     and if COMPZ = 'I' , Z contains the orthonormal eigenvectors
056  *     of the symmetric tridiagonal matrix.
057  *     If COMPZ = 'N' , then Z is not referenced.
058  
059  *     LDZ(input) INTEGER
060  *     The leading dimension of the array Z. LDZ >= 1 , and if
061  *     eigenvectors are desired , then LDZ >= max(1 , N).
062  
063  *     NR(input) INTEGER
064  *     NR = MAX(1 , NUMROC( N , NB , MYPROW , 0 , NPROCS ) ).
065  *     If COMPZ = 'N' , then NR is not referenced.
066  
067  *     WORK(workspace) DOUBLE PRECISION array , dimension(max(1 , 2*N - 2))
068  *     If COMPZ = 'N' , then WORK is not referenced.
069  
070  *     INFO(output) INTEGER
071  *     = 0 : successful exit
072  *     < 0 : if INFO = - i , the i - th argument had an illegal value
073  *     > 0 : the algorithm has failed to find all the eigenvalues in
074  *     a total of 30*N iterations ; if INFO = i , then i
075  *     elements of E have not converged to zero ; on exit , D
076  *     and E contain the elements of a symmetric tridiagonal
077  *     matrix which is orthogonally similar to the original
078  *     matrix.
079  
080  *     === ==================================================================
081  
082  *     .. Parameters ..
083        DOUBLE PRECISION ZERO , ONE , TWO , THREE
084        PARAMETER( ZERO = 0.0D0 , ONE = 1.0D0 , TWO = 2.0D0 ,
085       $THREE = 3.0D0 )
086        INTEGER MAXIT
087        PARAMETER( MAXIT = 30 )
088     90 CONTINUE
089        IF( L.NE.LEND ) THEN
089               
090            LENDP1 = LEND + 1
091            DO 100 M = L , LENDP1 , - 1
091                   
092                TST = ABS( E( M - 1 ) )**2
093                IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M - 1 ) ) +
094       $        SAFMIN )GO TO 110
095    100     CONTINUE
096        END IF
097  
098        M = LEND
099  
100    110 CONTINUE
101        IF( M.GT.LEND )
101               
102       $    E( M - 1 ) = ZERO
103            P = D( L )
104            IF( M.EQ.L )
104                   
105       $        GO TO 130
106  
107  *             If remaining matrix is 2 - by - 2 , use DLAE2 or SLAEV2
108  *             to compute its eigensystem.
109  
110                IF( M.EQ.L - 1 ) THEN
110                       
111                    IF( ICOMPZ.GT.0 ) THEN
111                           
112                        CALL DLAEV2( D( L - 1 ) , E( L - 1 ) , D( L ) , RT1 , RT2 , C , S )
113                        WORK( M ) = C
114                        WORK( N - 1 + M ) = S
115                        CALL DLASR( 'R' , 'V' , 'F' , NR , 2 , WORK( M ) ,
116       $                WORK( N - 1 + M ) , Z( 1 , L - 1 ) , LDZ )
117                    ELSE
117                           
118                        CALL DLAE2( D( L - 1 ) , E( L - 1 ) , D( L ) , RT1 , RT2 )
119                    END IF
120                    D( L - 1 ) = RT1
121                    D( L ) = RT2
122                    E( L - 1 ) = ZERO
123                    L = L - 2
124                    IF( L.GE.LEND )
124                           
125       $                GO TO 90
126                        GO TO 140
127                    END IF
128  
129                    IF( JTOT.EQ.NMAXIT )
129                           
130       $                GO TO 140
131                        JTOT = JTOT + 1
132  
133  *                     Form shift.
134  
135                        G =( D( L - 1 ) - P ) / ( TWO*E( L - 1 ) )
136                        R = DLAPY2( G , ONE )
137                        G = D( M ) - P + ( E( L - 1 ) / ( G + SIGN( R , G ) ) )
138  
139                        S = ONE
140                        C = ONE
141                        P = ZERO
142  
143  *                     Inner loop
144  
145                        LM1 = L - 1
146                        DO 120 I = M , LM1
146                               
147                            F = S*E( I )
148                            B = C*E( I )
149                            CALL DLARTG( G , F , C , S , R )
150                            IF( I.NE.M )
150                                   
151       $                        E( I - 1 ) = R
152                                G = D( I ) - P
153                                R =( D( I + 1 ) - G )*S + TWO*C*B
154                                P = S*R
155                                D( I ) = G + P
156                                G = C*R - B
157  
158  *                             If eigenvectors are desired , then save rotations.
159  
160                                IF( ICOMPZ.GT.0 ) THEN
160                                       
161                                    WORK( I ) = C
162                                    WORK( N - 1 + I ) = S
163                                END IF
164  
165    120                 CONTINUE
166  
167  *                     If eigenvectors are desired , then apply saved rotations.
168  
168                           
169                        IF( ICOMPZ.GT.0 ) THEN
169                               
170                            MM = L - M + 1
171                            CALL DLASR( 'R' , 'V' , 'F' , NR , MM , WORK( M ) , WORK( N - 1 + M ) ,
172       $                    Z( 1 , M ) , LDZ )
173                        END IF
174  
175                        D( L ) = D( L ) - P
176                        E( LM1 ) = G
177                        GO TO 90
178  
179  *                     Eigenvalue found.
180  
181    130 CONTINUE
182        D( L ) = P
183  
184        L = L - 1
185        IF( L.GE.LEND )
185               
186       $    GO TO 90
187            GO TO 140
188  
189        END IF
190  
191  *     Undo scaling if necessary
192  
193    140 CONTINUE
194        IF( ISCALE.EQ.1 ) THEN
194               
195            CALL DLASCL( 'G' , 0 , 0 , SSFMAX , ANORM , LENDSV - LSV + 1 , 1 ,
196       $    D( LSV ) , N , INFO )
197            CALL DLASCL( 'G' , 0 , 0 , SSFMAX , ANORM , LENDSV - LSV , 1 , E( LSV ) ,
198       $    N , INFO )
199        ELSE IF( ISCALE.EQ.2 ) THEN
199               
200            CALL DLASCL( 'G' , 0 , 0 , SSFMIN , ANORM , LENDSV - LSV + 1 , 1 ,
201       $    D( LSV ) , N , INFO )
202            CALL DLASCL( 'G' , 0 , 0 , SSFMIN , ANORM , LENDSV - LSV , 1 , E( LSV ) ,
203       $    N , INFO )
204        END IF
205  
206  *     Check for no convergence to an eigenvalue after a total
207  *     of N*MAXIT iterations.
208  
209        IF( JTOT.LT.NMAXIT )
209               
210       $    GO TO 10
211            DO 150 I = 1 , N - 1
211                   
212                IF( E( I ).NE.ZERO )
212                       
213       $            INFO = INFO + 1
214    150     CONTINUE
214               
215            GO TO 190
216  
217  *         Order eigenvalues and eigenvectors.
218  
219    160 CONTINUE
220        IF( ICOMPZ.EQ.0 ) THEN
221  
222  *         Use Quick Sort
223  
223               
224            CALL DLASRT( 'I' , N , D , INFO )
225  
226        ELSE
227  
228  *         Use Selection Sort to minimize swaps of eigenvectors
229  
229               
230            DO 180 II = 2 , N
230                   
231                I = II - 1
232                K = I
233                P = D( I )
234                DO 170 J = II , N
234                       
235                    IF( D( J ).LT.P ) THEN
235                           
236                        K = J
237                        P = D( J )
238                    END IF
239    170         CONTINUE
239                   
240                IF( K.NE.I ) THEN
240                       
241                    D( K ) = D( I )
242                    D( I ) = P
243                    CALL DSWAP( NR , Z( 1 , I ) , 1 , Z( 1 , K ) , 1 )
244                END IF
245    180     CONTINUE
246        END IF
247  
248    190 CONTINUE
249        RETURN
250  
251  *     End of DSTEQR2
252  
253        END7650

Variables in Routine DSTEQR2()

Summary Report
Data Type	Quantity	Size(byte)
CHARACTER	1	1
DOUBLE PRECISION	8	32
INTEGER	15	60
REAL	8	32
TOTAL	32	125

List of Variables

CHARACTER

COMPZ

DOUBLE PRECISION

D( * ) E( * ) ONE THREE TWO
WORK( * ) Z( LDZ, * ) ZERO

INTEGER

I II INFO J JTOT
K L LDZ LENDP1 LM1
M MAXIT MM N NR

REAL

B C F G P
R S TST

Variables Dependence Graph
Put the mouse over a right hand side variable to display the corresponding line of the dependence
-
- -
B <--- CB = C*E( I ), EB = C*E( I ), IB = C*E( I )

C <--- ONEC = ONE

D <--- LD( L ) = D( L ) - P, PD( I ) = G + P{2D( L ) = D( L ) - P, 3D( L ) = P, 4D( I ) = P}, DD( L ) = D( L ) - P{2D( K ) = D( I )}, GD( I ) = G + P, ID( K ) = D( I )

E <--- ZEROE( L-1 ) = ZERO, GE( LM1 ) = G

F <--- SF = S*E( I ), EF = S*E( I ), IF = S*E( I )

G <--- BG = C*R - B, LG = ( D( L-1 )-P ) / ( TWO*E( L-1 ) ){2}, M, CG = C*R - B, PG = ( D( L-1 )-P ) / ( TWO*E( L-1 ) ){2, 3G = D( I ) - P}, R{2G = C*R - B}, TWOG = ( D( L-1 )-P ) / ( TWO*E( L-1 ) ), DG = ( D( L-1 )-P ) / ( TWO*E( L-1 ) ){2, 3G = D( I ) - P}, EG = ( D( L-1 )-P ) / ( TWO*E( L-1 ) ){2}, G, IG = D( I ) - P

I <--- LM1DO 120 I = M, LM1, MDO 120 I = M, LM1, NDO 150 I = 1, N - 1, III = II - 1

II <--- NDO 180 II = 2, N

J <--- NDO 170 J = II, N, IIDO 170 J = II, N

JTOT <--- JTOTJTOT = JTOT + 1

K <--- JK = J, IK = I

L <--- LL = L - 2{2L = L - 1}

LM1 <--- LLM1 = L - 1

M <--- LDO 100 M = L, LENDP1, -1, LENDP1DO 100 M = L, LENDP1, -1

MM <--- LMM = L - M + 1, MMM = L - M + 1

P <--- JP = D( J ), LP = D( L ), RP = S*R, SP = S*R, ZEROP = ZERO, DP = D( L ){2P = D( I ), 3P = D( J )}, IP = D( I )

R <--- BR = ( D( I+1 )-G )*S + TWO*C*B, CR = ( D( I+1 )-G )*S + TWO*C*B, ONER = DLAPY2( G, ONE ), SR = ( D( I+1 )-G )*S + TWO*C*B, TWOR = ( D( I+1 )-G )*S + TWO*C*B, DR = ( D( I+1 )-G )*S + TWO*C*B, GR = DLAPY2( G, ONE ){2R = ( D( I+1 )-G )*S + TWO*C*B}, IR = ( D( I+1 )-G )*S + TWO*C*B

S <--- ONES = ONE

TST <--- MTST = ABS( E( M-1 ) )**2, ETST = ABS( E( M-1 ) )**2

WORK <--- CWORK( M ) = C{2WORK( I ) = C}, SWORK( N-1+M ) = S{2WORK( N-1+I ) = S}

Analysis elements of the routine DSTEQR2()
Put the mouse over each element to display detailed matching information

Assigned variables
		B , C , COMPZ , F , G , I , II , INFO , J , JTOT , K , L , LDZ , LENDP1 , LM1 , M , MAXIT , MM , N , NR , ONE , P , R , S , THREE , TST , TWO , ZERO

Active variables
		B , C , COMPZ , D , E , F , G , I , II , INFO , J , JTOT , K , L , LDZ , LENDP1 , LM1 , M , MAXIT , MM , N , NR , ONE , P , R , s , THREE , TST , TWO , WORK , Z , zero

Accessed arrays [ array name : associated index ]
	D	: * , I , I , I , I , I , I+1 , J , J , K , L , L , L , L , L , L , L-1 , L-1 , L-1 , L-1 , LSV , LSV , M , M , M-1
	E	: * , I , I , I , I-1 , L-1 , L-1 , L-1 , L-1 , L-1 , LM1 , LSV , LSV , M-1 , M-1
	WORK	: * , I , M , M , M , N-1+I , N-1+M , N-1+M , N-1+M
	Z	: 1, I , 1, K , 1, L-1 , 1, M , LDZ, *

Conditional statements [ statement : associated predicate ]
	do	: `(` 100 M = L , LENDP1 , - 1 `)` , `(` 120 I = M , LM1 `)` , `(` 150 I = 1 , N - 1 `)` , `(` 180 II = 2 , N `)` , `(` 170 J = II , N `)`
	for	: `(` no convergence to an eigenvalue after a total `)`
	if	: `(` INFO = 0 , the eigenvalues in ascending order. `)` , `(` COMPZ = 'V' , `)` , `(` INFO = 0 , `)` , `(` COMPZ = 'V' , Z contains the `)` , `(` COMPZ = 'I' , Z contains the orthonormal eigenvectors `)` , `(` COMPZ = 'N' , `)` , `(` COMPZ = 'N' , `)` , `(` COMPZ = 'N' , `)` , `(` INFO = - i , the i - th argument had an illegal value `)` , `(` INFO = i , `)` , `(` L.NE.LEND `)` , `(` (TST.LE.( EPS2ABS( D( M ) ) )ABS( D( M - 1 ) ) + `)` , `(` M.GT.LEND `)` , `(` M.EQ.L `)` , `(` remaining matrix is 2 - by - 2 , use DLAE2 or SLAEV2 `)` , `(` M.EQ.L - 1 `)` , `(` ICOMPZ.GT.0 `)` , `(` L.GE.LEND `)` , `(` JTOT.EQ.NMAXIT `)` , `(` I.NE.M `)` , `(` eigenvectors are desired , `)` , `(` ICOMPZ.GT.0 `)` , `(` eigenvectors are desired , `)` , `(` ICOMPZ.GT.0 `)` , `(` L.GE.LEND `)` , `(` necessary `)` , `(` ISCALE.EQ.1 `)` , `(` ISCALE.EQ.2 `)` , `(` JTOT.LT.NMAXIT `)` , `(` (E( I ).NE.ZERO ) `)` , `(` ICOMPZ.EQ.0 `)` , `(` (D( J ).LT.P ) `)` , `(` K.NE.I `)`