1* Hermitian matrix. On entry, Z must contain the 2* unitary matrix used to reduce the original matrix 3* to tridiagonal form. 4* = 'I': Compute eigenvalues and eigenvectors of the 5* tridiagonal matrix. Z is initialized to the identity 6* matrix. 7* 8* N (input) INTEGER 9* The order of the matrix. N >= 0. 10* 11* D (input/output) DOUBLE PRECISION array, dimension (N) 12* On entry, the diagonal elements of the tridiagonal matrix. 13* On exit, if INFO = 0, the eigenvalues in ascending order. 14* 15* E (input/output) DOUBLE PRECISION array, dimension (N-1) 16* On entry, the (n-1) subdiagonal elements of the tridiagonal 17* matrix. 18* On exit, E has been destroyed. 19* 20* Z (input/output) COMPLEX*16 array, dimension (LDZ, N) 21* On entry, if COMPZ = 'V', then Z contains the unitary 22* matrix used in the reduction to tridiagonal form. 23* On exit, if INFO = 0, then if COMPZ = 'V', Z contains the 24* orthonormal eigenvectors of the original Hermitian matrix, 25* and if COMPZ = 'I', Z contains the orthonormal eigenvectors 26* of the symmetric tridiagonal matrix. 27* If COMPZ = 'N', then Z is not referenced. 28* 29* LDZ (input) INTEGER 30* The leading dimension of the array Z. LDZ >= 1, and if 31* eigenvectors are desired, then LDZ >= max(1,N). 32* 33* WORK (workspace) DOUBLE PRECISION array, dimension (max(1,2*N-2)) 34* If COMPZ = 'N', then WORK is not referenced. 35* 36* INFO (output) INTEGER 37* = 0: successful exit 38* < 0: if INFO = -i, the i-th argument had an illegal value 39* > 0: the algorithm has failed to find all the eigenvalues in 40* a total of 30*N iterations; if INFO = i, then i 41* elements of E have not converged to zero; on exit, D 42* and E contain the elements of a symmetric tridiagonal 43* matrix which is unitarily similar to the original 44* matrix. 45* 46* ===================================================================== 47* 48* .. Parameters .. 49 DOUBLE PRECISION ZERO, ONE, TWO, THREE 50 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, 51 $ THREE = 3.0D0 ) 52 COMPLEX*16 CZERO, CONE 53 PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ), 54 $ CONE = ( 1.0D0, 0.0D0 ) ) 55 INTEGER MAXIT, KTOTAL 56 PARAMETER ( MAXIT = 30 ) 57* .. 58* .. Local Scalars .. 59 INTEGER I, ICOMPZ, II, ISCALE, J, JTOT, K, L, L1, LEND, 60 $ LENDM1, LENDP1, LENDSV, LM1, LSV, M, MM, MM1, 61 $ NM1, NMAXIT 62 DOUBLE PRECISION ANORM, B, C, EPS, EPS2, F, G, P, R, RT1, RT2, 63 $ S, SAFMAX, SAFMIN, SSFMAX, SSFMIN, TST 64* .. 65* .. External Functions .. 66 LOGICAL LSAME 67 DOUBLE PRECISION DLAMCH, DLANST, DLAPY2 68 EXTERNAL LSAME, DLAMCH, DLANST, DLAPY2 69* .. 70* .. External Subroutines .. 71 EXTERNAL DLAE2, DLAEV2, DLARTG, DLASCL, DLASRT, XERBLA, 72 $ ZLASET, ZLASR, ZSWAP 73* .. 74* .. Intrinsic Functions .. 75 INTRINSIC ABS, MAX, SIGN, SQRT 76* .. 77* .. Executable Statements .. 78* 79* Test the input parameters. 80* 81 INFO = 0 82 KTOTAL = 0 83* 84 IF( LSAME( COMPZ, 'N' ) ) THEN 85 ICOMPZ = 0 86 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN 87 ICOMPZ = 1 88 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN 89 ICOMPZ = 2 90 ELSE 91 ICOMPZ = -1 92 END IF 93 IF( ICOMPZ.LT.0 ) THEN 94 INFO = -1 95 ELSE IF( N.LT.0 ) THEN 96 INFO = -2 97 ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1, 98 $ N ) ) ) THEN 99 INFO = -6 100 END IF 101 IF( INFO.NE.0 ) THEN 102 CALL XERBLA( 'ZSTEQR', -INFO ) 103 RETURN 104 END IF 105* 106* Quick return if possible 107* 108 IF( N.EQ.0 ) 109 $ RETURN 110* 111 IF( N.EQ.1 ) THEN 112 IF( ICOMPZ.EQ.2 ) 113 $ Z( 1, 1 ) = CONE 114 RETURN 115 END IF 116* 117* Determine the unit roundoff and over/underflow thresholds. 118* 119 EPS = DLAMCH( 'E' ) 120 EPS2 = EPS**2 121 SAFMIN = DLAMCH( 'S' ) 122 SAFMAX = ONE / SAFMIN 123 SSFMAX = SQRT( SAFMAX ) / THREE 124 SSFMIN = SQRT( SAFMIN ) / EPS2 125* 126* Compute the eigenvalues and eigenvectors of the tridiagonal 127* matrix. 128* 129 IF( ICOMPZ.EQ.2 ) 130 $ CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDZ ) 131* 132 NMAXIT = N*MAXIT 133 JTOT = 0 134* 135* Determine where the matrix splits and choose QL or QR iteration 136* for each block, according to whether top or bottom diagonal 137* element is smaller. 138* 139 L1 = 1 140 NM1 = N - 1 141* 142 10 CONTINUE 143 IF( L1.GT.N ) 144 $ GO TO 160 145 IF( L1.GT.1 ) 146 $ E( L1-1 ) = ZERO 147 IF( L1.LE.NM1 ) THEN 148 DO 20 M = L1, NM1 149 TST = ABS( E( M ) ) 150 IF( TST.EQ.ZERO ) 151 $ GO TO 30 152 IF( TST.LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+ 153 $ 1 ) ) ) )*EPS ) THEN 154 E( M ) = ZERO 155 GO TO 30 156 END IF 157 20 CONTINUE 158 END IF 159 M = N 160* 161 30 CONTINUE 162 L = L1 163 LSV = L 164 LEND = M 165 LENDSV = LEND 166 L1 = M + 1 167 IF( LEND.EQ.L ) 168 $ GO TO 10 169* 170* Scale submatrix in rows and columns L to LEND 171* 172 ANORM = DLANST( 'I', LEND-L+1, D( L ), E( L ) ) 173 ISCALE = 0 174* print*,'checking if scaling is needed. L = ', L 175* print*,' LEND = ', LEND 176* print*,' ANORM = ', ANORM 177 IF( ANORM.EQ.ZERO ) THEN 178 print*,'anorm = ', ANORM 179 $ GO TO 10 180 END IF 181 IF( ANORM.GT.SSFMAX ) THEN 182 ISCALE = 1 183 CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N, 184 $ INFO ) 185 CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N, 186 $ INFO ) 187 print*,'iscale = 1. anorm = ', ANORM 188 ELSE IF( ANORM.LT.SSFMIN ) THEN 189 ISCALE = 2 190 CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N, 191 $ INFO ) 192 CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N, 193 $ INFO ) 194 print*,'iscale = 2. anorm = ', ANORM 195 END IF 196* 197* Choose between QL and QR iteration 198* 199 IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN 200 LEND = LSV 201 L = LENDSV 202 END IF 203* 204 IF( .TRUE. ) THEN 205* 206* QR Iteration 207* 208* Look for small superdiagonal element. 209* 210 90 CONTINUE 211 IF( L.NE.LEND ) THEN 212 LENDP1 = LEND + 1 213 DO 100 M = L, LENDP1, -1 214 TST = ABS( E( M-1 ) )**2 215 IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M-1 ) )+ 216 $ SAFMIN )GO TO 110 217 100 CONTINUE 218 END IF 219* 220 M = LEND 221* 222 110 CONTINUE 223 IF( M.GT.LEND ) 224 $ E( M-1 ) = ZERO 225 P = D( L ) 226 IF( M.EQ.L ) 227 $ GO TO 130 228* 229* If remaining matrix is 2-by-2, use DLAE2 or SLAEV2 230* to compute its eigensystem. 231* 232 IF( M.EQ.L-1 ) THEN 233 IF( ICOMPZ.GT.0 ) THEN 234 CALL DLAEV2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2, C, S ) 235 WORK( M ) = C 236 WORK( N-1+M ) = S 237 CALL ZLASR( 'R', 'V', 'F', N, 2, WORK( M ), 238 $ WORK( N-1+M ), Z( 1, L-1 ), LDZ ) 239 ELSE 240 CALL DLAE2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2 ) 241 END IF 242 D( L-1 ) = RT1 243 D( L ) = RT2 244 E( L-1 ) = ZERO 245 L = L - 2 246 IF( L.GE.LEND ) 247 $ GO TO 90 248 GO TO 140 249 END IF 250* 251 IF( JTOT.EQ.NMAXIT ) 252 $ GO TO 140 253 JTOT = JTOT + 1 254* 255* Form shift. 256* 257 G = ( D( L-1 )-P ) / ( TWO*E( L-1 ) ) 258 R = DLAPY2( G, ONE ) 259 G = D( M ) - P + ( E( L-1 ) / ( G+SIGN( R, G ) ) ) 260* 261 S = ONE 262 C = ONE 263 P = ZERO 264* 265* Inner loop 266* 267 LM1 = L - 1 268 DO 120 I = M, LM1 269 F = S*E( I ) 270 B = C*E( I ) 271 CALL DLARTG( G, F, C, S, R ) 272 IF( I.NE.M ) 273 $ E( I-1 ) = R 274 G = D( I ) - P 275 R = ( D( I+1 )-G )*S + TWO*C*B 276 P = S*R 277 D( I ) = G + P 278 G = C*R - B 279* 280* If eigenvectors are desired, then save rotations. 281* 282 IF( ICOMPZ.GT.0 ) THEN 283 WORK( I ) = C 284 WORK( N-1+I ) = S 285 END IF 286* 287 120 CONTINUE 288* 289* If eigenvectors are desired, then apply saved rotations. 290* 291 IF( ICOMPZ.GT.0 ) THEN 292 MM = L - M + 1 293 CALL ZLASR( 'R', 'V', 'F', N, MM, WORK( M ), WORK( N-1+M ), 294 $ Z( 1, M ), LDZ ) 295 KTOTAL = KTOTAL + 1 296 END IF 297* 298 D( L ) = D( L ) - P 299 E( LM1 ) = G 300 GO TO 90 301* 302* Eigenvalue found. 303* 304 130 CONTINUE 305 D( L ) = P 306* 307 L = L - 1 308 IF( L.GE.LEND ) 309 $ GO TO 90 310 GO TO 140 311* 312 END IF 313* 314* Undo scaling if necessary 315* 316 140 CONTINUE 317 IF( ISCALE.EQ.1 ) THEN 318 CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1, 319 $ D( LSV ), N, INFO ) 320 CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV, 1, E( LSV ), 321 $ N, INFO ) 322 ELSE IF( ISCALE.EQ.2 ) THEN 323 CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1, 324 $ D( LSV ), N, INFO ) 325 CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV, 1, E( LSV ), 326 $ N, INFO ) 327 END IF 328* 329* Check for no convergence to an eigenvalue after a total 330* of N*MAXIT iterations. 331* 332 IF( JTOT.EQ.NMAXIT ) THEN 333 DO 150 I = 1, N - 1 334 IF( E( I ).NE.ZERO ) 335 $ INFO = INFO + 1 336 150 CONTINUE 337 RETURN 338 END IF 339 GO TO 10 340* 341* Order eigenvalues and eigenvectors. 342* 343 160 CONTINUE 344c print *, 'K_TOTAL = ', KTOTAL 345 IF( ICOMPZ.EQ.0 ) THEN 346* 347* Use Quick Sort 348* 349 CALL DLASRT( 'I', N, D, INFO ) 350* 351 ELSE 352* 353* Use Selection Sort to minimize swaps of eigenvectors 354* 355 DO 180 II = 2, N 356 I = II - 1 357 K = I 358 P = D( I ) 359 DO 170 J = II, N 360 IF( D( J ).LT.P ) THEN 361 K = J 362 P = D( J ) 363 END IF 364 170 CONTINUE 365 IF( K.NE.I ) THEN 366 D( K ) = D( I ) 367 D( I ) = P 368 CALL ZSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 ) 369 END IF 370 180 CONTINUE 371 END IF 372 RETURN 373* 374* End of ZSTEQR 375* 376 END 377