SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zsteqr2.f
Go to the documentation of this file.
1 SUBROUTINE zsteqr2( COMPZ, N, D, E, Z, LDZ, NR, WORK, INFO )
2*
3* -- ScaLAPACK routine (version 1.7) --
4* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
5* Courant Institute, Argonne National Lab, and Rice University
6* November 15, 1997
7*
8* .. Scalar Arguments ..
9 CHARACTER COMPZ
10 INTEGER INFO, LDZ, N, NR
11* ..
12* .. Array Arguments ..
13 DOUBLE PRECISION D( * ), E( * ), WORK( * )
14 COMPLEX*16 Z( LDZ, * )
15* ..
16*
17* Purpose
18* =======
19*
20* ZSTEQR2 is a modified version of LAPACK routine ZSTEQR.
21* ZSTEQR2 computes all eigenvalues and, optionally, eigenvectors of a
22* symmetric tridiagonal matrix using the implicit QL or QR method.
23* ZSTEQR2 is modified from ZSTEQR to allow each ScaLAPACK process
24* running ZSTEQR2 to perform updates on a distributed matrix Q.
25* Proper usage of ZSTEQR2 can be gleaned from
26* examination of ScaLAPACK's * PZHEEV.
27* ZSTEQR2 incorporates changes attributed to Greg Henry.
28*
29* Arguments
30* =========
31*
32* COMPZ (input) CHARACTER*1
33* = 'N': Compute eigenvalues only.
34* = 'I': Compute eigenvalues and eigenvectors of the
35* tridiagonal matrix. Z must be initialized to the
36* identity matrix by PZLASET or ZLASET prior
37* to entering this subroutine.
38*
39* N (input) INTEGER
40* The order of the matrix. N >= 0.
41*
42* D (input/output) DOUBLE PRECISION array, dimension (N)
43* On entry, the diagonal elements of the tridiagonal matrix.
44* On exit, if INFO = 0, the eigenvalues in ascending order.
45*
46* E (input/output) DOUBLE PRECISION array, dimension (N-1)
47* On entry, the (n-1) subdiagonal elements of the tridiagonal
48* matrix.
49* On exit, E has been destroyed.
50*
51* Z (local input/local output) COMPLEX*16 array, global
52* dimension (N, N), local dimension (LDZ, NR).
53* On entry, if COMPZ = 'V', then Z contains the orthogonal
54* matrix used in the reduction to tridiagonal form.
55* On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
56* orthonormal eigenvectors of the original symmetric matrix,
57* and if COMPZ = 'I', Z contains the orthonormal eigenvectors
58* of the symmetric tridiagonal matrix.
59* If COMPZ = 'N', then Z is not referenced.
60*
61* LDZ (input) INTEGER
62* The leading dimension of the array Z. LDZ >= 1, and if
63* eigenvectors are desired, then LDZ >= max(1,N).
64*
65* NR (input) INTEGER
66* NR = MAX(1, NUMROC( N, NB, MYPROW, 0, NPROCS ) ).
67* If COMPZ = 'N', then NR is not referenced.
68*
69* WORK (workspace) DOUBLE PRECISION array, dimension (max(1,2*N-2))
70* If COMPZ = 'N', then WORK is not referenced.
71*
72* INFO (output) INTEGER
73* = 0: successful exit
74* < 0: if INFO = -i, the i-th argument had an illegal value
75* > 0: the algorithm has failed to find all the eigenvalues in
76* a total of 30*N iterations; if INFO = i, then i
77* elements of E have not converged to zero; on exit, D
78* and E contain the elements of a symmetric tridiagonal
79* matrix which is orthogonally similar to the original
80* matrix.
81*
82* =====================================================================
83*
84* .. Parameters ..
85 DOUBLE PRECISION ZERO, ONE, TWO, THREE, HALF
86 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0,
87 $ three = 3.0d0, half = 0.5d0 )
88 COMPLEX*16 CONE
89 parameter( cone = ( 1.0d0, 1.0d0 ) )
90 INTEGER MAXIT, NMAXLOOK
91 parameter( maxit = 30, nmaxlook = 15 )
92* ..
93* .. Local Scalars ..
94 INTEGER I, ICOMPZ, II, ILAST, ISCALE, J, JTOT, K, L,
95 $ L1, LEND, LENDM1, LENDP1, LENDSV, LM1, LSV, M,
96 $ MM, MM1, NLOOK, NM1, NMAXIT
97 DOUBLE PRECISION ANORM, B, C, EPS, EPS2, F, G, GP, OLDEL, OLDGP,
98 $ OLDRP, P, R, RP, RT1, RT2, S, SAFMAX, SAFMIN,
99 $ 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*
358* Skip the current step: the subdiagonal info is just noise.
359*
360 m = l
361 e( m ) = zero
362 p = d( l )
363 jtot = jtot - 1
364 GOTO 110
365 ENDIF
366 g = gp
367 r = rp
368*
369* Lookahead over
370*
371 90 CONTINUE
372*
373 s = one
374 c = one
375 p = zero
376*
377* Inner loop
378*
379 mm1 = m - 1
380 DO 100 i = mm1, l, -1
381 f = s*e( i )
382 b = c*e( i )
383 CALL dlartg( g, f, c, s, r )
384 IF( i.NE.m-1 )
385 $ e( i+1 ) = r
386 g = d( i+1 ) - p
387 r = ( d( i )-g )*s + two*c*b
388 p = s*r
389 d( i+1 ) = g + p
390 g = c*r - b
391*
392* If eigenvectors are desired, then save rotations.
393*
394 work( i ) = c
395 work( n-1+i ) = -s
396*
397 100 CONTINUE
398*
399* If eigenvectors are desired, then apply saved rotations.
400*
401 mm = m - l + 1
402 CALL zlasr( 'R', 'V', 'B', nr, mm, work( l ), work( n-1+l ),
403 $ z( 1, l ), ldz )
404*
405 d( l ) = d( l ) - p
406 e( l ) = g
407 ilast = l
408 GOTO 40
409*
410* Eigenvalue found.
411*
412 110 CONTINUE
413 d( l ) = p
414*
415 l = l + 1
416 IF( l.LE.lend )
417 $ GOTO 40
418 GOTO 200
419*
420 ELSE
421*
422* QR Iteration
423*
424* Look for small superdiagonal element.
425*
426 120 CONTINUE
427 IF( l.NE.lend ) THEN
428 lendp1 = lend + 1
429 DO 130 m = l, lendp1, -1
430 tst = abs( e( m-1 ) )**2
431 IF( tst.LE.( eps2*abs( d( m ) ) )*abs( d( m-1 ) )+
432 $ safmin )GOTO 140
433 130 CONTINUE
434 ENDIF
435*
436 m = lend
437*
438 140 CONTINUE
439 IF( m.GT.lend )
440 $ e( m-1 ) = zero
441 p = d( l )
442 IF( m.EQ.l )
443 $ GOTO 190
444*
445* If remaining matrix is 2-by-2, use DLAE2 or DLAEV2
446* to compute its eigensystem.
447*
448 IF( m.EQ.l-1 ) THEN
449 CALL dlaev2( d( l-1 ), e( l-1 ), d( l ), rt1, rt2, c, s )
450 work( m ) = c
451 work( n-1+m ) = s
452 CALL zlasr( 'R', 'V', 'F', nr, 2, work( m ), work( n-1+m ),
453 $ z( 1, l-1 ), ldz )
454 d( l-1 ) = rt1
455 d( l ) = rt2
456 e( l-1 ) = zero
457 l = l - 2
458 IF( l.GE.lend )
459 $ GOTO 120
460 GOTO 200
461 ENDIF
462*
463 IF( jtot.EQ.nmaxit )
464 $ GOTO 200
465 jtot = jtot + 1
466*
467* Form shift.
468*
469 g = ( d( l-1 )-p ) / ( two*e( l-1 ) )
470 r = dlapy2( g, one )
471 g = d( m ) - p + ( e( l-1 ) / ( g+sign( r, g ) ) )
472*
473 IF( icompz.EQ.0 ) THEN
474* Do not do a lookahead!
475 GOTO 170
476 ENDIF
477*
478 oldel = abs( e( l-1 ) )
479 gp = g
480 rp = r
481 tst = abs( e( l-1 ) )**2
482 tst = tst / ( ( eps2*abs( d( l ) ) )*abs( d( l-1 ) )+safmin )
483 nlook = 1
484 IF( ( tst.GT.one ) .AND. ( nlook.LE.nmaxlook ) ) THEN
485 150 CONTINUE
486*
487* This is the lookahead loop, going until we have
488* convergence or too many steps have been taken.
489*
490 s = one
491 c = one
492 p = zero
493*
494* Inner loop
495*
496 lm1 = l - 1
497 DO 160 i = m, lm1
498 f = s*e( i )
499 b = c*e( i )
500 CALL dlartg( gp, f, c, s, rp )
501 gp = d( i ) - p
502 rp = ( d( i+1 )-gp )*s + two*c*b
503 p = s*rp
504 IF( i.LT.lm1 )
505 $ gp = c*rp - b
506 160 CONTINUE
507 oldgp = gp
508 oldrp = rp
509* Find GP & RP for the next iteration
510 IF( abs( c*oldrp-b ).GT.safmin ) THEN
511 gp = ( ( oldgp+p )-( d( l )-p ) ) / ( two*( c*oldrp-b ) )
512 ELSE
513*
514* Goto put in by G. Henry to fix ALPHA problem
515*
516 GOTO 170
517* GP = ( ( OLDGP+P )-( D( L )-P ) ) /
518* $ ( TWO*( C*OLDRP-B )+SAFMIN )
519 ENDIF
520 rp = dlapy2( gp, one )
521 gp = d( m ) - ( d( l )-p ) +
522 $ ( ( c*oldrp-b ) / ( gp+sign( rp, gp ) ) )
523 tst1 = tst
524 tst = abs( ( c*oldrp-b ) )**2
525 tst = tst / ( ( eps2*abs( d( l )-p ) )*abs( oldgp+p )+
526 $ safmin )
527* Make sure that we are making progress
528 IF( abs( c*oldrp-b ).GT.0.9d0*oldel ) THEN
529 IF( abs( c*oldrp-b ).GT.oldel ) THEN
530 gp = g
531 rp = r
532 ENDIF
533 tst = half
534 ELSE
535 oldel = abs( c*oldrp-b )
536 ENDIF
537 nlook = nlook + 1
538 IF( ( tst.GT.one ) .AND. ( nlook.LE.nmaxlook ) )
539 $ GOTO 150
540 ENDIF
541 IF( ( tst.LE.one ) .AND. ( tst.NE.half ) .AND.
542 $ ( abs( p ).LT.eps*abs( d( l ) ) ) .AND.
543 $ ( ilast.EQ.l ) .AND. ( abs( e( l-1 ) )**2.LE.10000.0d0*
544 $ ( ( eps2*abs( d( l-1 ) ) )*abs( d( l ) )+safmin ) ) ) THEN
545*
546* Skip the current step: the subdiagonal info is just noise.
547*
548 m = l
549 e( m-1 ) = zero
550 p = d( l )
551 jtot = jtot - 1
552 GOTO 190
553 ENDIF
554*
555 g = gp
556 r = rp
557*
558* Lookahead over
559*
560 170 CONTINUE
561*
562 s = one
563 c = one
564 p = zero
565 DO 180 i = m, lm1
566 f = s*e( i )
567 b = c*e( i )
568 CALL dlartg( g, f, c, s, r )
569 IF( i.NE.m )
570 $ e( i-1 ) = r
571 g = d( i ) - p
572 r = ( d( i+1 )-g )*s + two*c*b
573 p = s*r
574 d( i ) = g + p
575 g = c*r - b
576*
577* If eigenvectors are desired, then save rotations.
578*
579 work( i ) = c
580 work( n-1+i ) = s
581*
582 180 CONTINUE
583*
584* If eigenvectors are desired, then apply saved rotations.
585*
586 mm = l - m + 1
587 CALL zlasr( 'R', 'V', 'F', nr, mm, work( m ), work( n-1+m ),
588 $ z( 1, m ), ldz )
589*
590 d( l ) = d( l ) - p
591 e( lm1 ) = g
592 ilast = l
593 GOTO 120
594*
595* Eigenvalue found.
596*
597 190 CONTINUE
598 d( l ) = p
599*
600 l = l - 1
601 IF( l.GE.lend )
602 $ GOTO 120
603 GOTO 200
604*
605 ENDIF
606*
607* Undo scaling if necessary
608*
609 200 CONTINUE
610 IF( iscale.EQ.1 ) THEN
611 CALL dlascl( 'G', 0, 0, ssfmax, anorm, lendsv-lsv+1, 1,
612 $ d( lsv ), n, info )
613 CALL dlascl( 'G', 0, 0, ssfmax, anorm, lendsv-lsv, 1, e( lsv ),
614 $ n, info )
615 ELSEIF( iscale.EQ.2 ) THEN
616 CALL dlascl( 'G', 0, 0, ssfmin, anorm, lendsv-lsv+1, 1,
617 $ d( lsv ), n, info )
618 CALL dlascl( 'G', 0, 0, ssfmin, anorm, lendsv-lsv, 1, e( lsv ),
619 $ n, info )
620 ENDIF
621*
622* Check for no convergence to an eigenvalue after a total
623* of N*MAXIT iterations.
624*
625 IF( jtot.LT.nmaxit )
626 $ GOTO 10
627 DO 210 i = 1, n - 1
628 IF( e( i ).NE.zero )
629 $ info = info + 1
630 210 CONTINUE
631 GOTO 250
632*
633* Order eigenvalues and eigenvectors.
634*
635 220 CONTINUE
636*
637* Use Selection Sort to minimize swaps of eigenvectors
638*
639 DO 240 ii = 2, n
640 i = ii - 1
641 k = i
642 p = d( i )
643 DO 230 j = ii, n
644 IF( d( j ).LT.p ) THEN
645 k = j
646 p = d( j )
647 ENDIF
648 230 CONTINUE
649 IF( k.NE.i ) THEN
650 d( k ) = d( i )
651 d( i ) = p
652 CALL zswap( nr, z( 1, i ), 1, z( 1, k ), 1 )
653 ENDIF
654 240 CONTINUE
655*
656 250 CONTINUE
657* WRITE( *, FMT = * )'JTOT', JTOT
658 RETURN
659*
660* End of DSTEQR2
661*
662 END
max
#define max(A, B)
Definition pcgemr.c:180
zsteqr2
subroutine zsteqr2(compz, n, d, e, z, ldz, nr, work, info)
Definition zsteqr2.f:2