ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
csteqr2.f
Go to the documentation of this file.
1  SUBROUTINE csteqr2( 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  REAL D( * ), E( * ), WORK( * )
14  COMPLEX Z( LDZ, * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * CSTEQR2 is a modified version of LAPACK routine CSTEQR.
21 * CSTEQR2 computes all eigenvalues and, optionally, eigenvectors of a
22 * symmetric tridiagonal matrix using the implicit QL or QR method.
23 * CSTEQR2 is modified from CSTEQR to allow each ScaLAPACK process
24 * running CSTEQR2 to perform updates on a distributed matrix Q.
25 * Proper usage of CSTEQR2 can be gleaned from
26 * examination of ScaLAPACK's * PCHEEV.
27 * CSTEQR2 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 PCLASET or CLASET prior
37 * to entering this subroutine.
38 *
39 * N (input) INTEGER
40 * The order of the matrix. N >= 0.
41 *
42 * D (input/output) REAL 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) REAL 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 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) REAL 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  REAL ZERO, ONE, TWO, THREE, HALF
86  parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0,
87  $ three = 3.0e0, half = 0.5e0 )
88  COMPLEX CONE
89  parameter( cone = ( 1.0e0, 1.0e0 ) )
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  REAL 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  REAL SLAMCH, SLANST, SLAPY2
104  EXTERNAL lsame, slamch, slanst, slapy2
105 * ..
106 * .. External Subroutines ..
107  EXTERNAL clasr, cswap, slaev2, slartg, slascl, ssterf,
108  $ xerbla
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( 'CSTEQR2', -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 ssterf( 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 = slamch( 'E' )
159  eps2 = eps**2
160  safmin = slamch( '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 = slanst( '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 slascl( 'G', 0, 0, anorm, ssfmax, lend-l+1, 1, d( l ), n,
215  $ info )
216  CALL slascl( 'G', 0, 0, anorm, ssfmax, lend-l, 1, e( l ), n,
217  $ info )
218  ELSEIF( anorm.LT.ssfmin ) THEN
219  iscale = 2
220  CALL slascl( 'G', 0, 0, anorm, ssfmin, lend-l+1, 1, d( l ), n,
221  $ info )
222  CALL slascl( '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 SLAE2 or SLAEV2
259 * to compute its eigensystem.
260 *
261  IF( m.EQ.l+1 ) THEN
262  CALL slaev2( d( l ), e( l ), d( l+1 ), rt1, rt2, c, s )
263  work( l ) = c
264  work( n-1+l ) = s
265  CALL clasr( '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 = slapy2( 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 slartg( 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 = slapy2( 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.9e0*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.0e0*
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 slartg( 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 clasr( '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 SLAE2 or SLAEV2
446 * to compute its eigensystem.
447 *
448  IF( m.EQ.l-1 ) THEN
449  CALL slaev2( d( l-1 ), e( l-1 ), d( l ), rt1, rt2, c, s )
450  work( m ) = c
451  work( n-1+m ) = s
452  CALL clasr( '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 = slapy2( 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 slartg( 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 = slapy2( 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.9e0*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.0e0*
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 slartg( 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 clasr( '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 slascl( 'G', 0, 0, ssfmax, anorm, lendsv-lsv+1, 1,
612  $ d( lsv ), n, info )
613  CALL slascl( 'G', 0, 0, ssfmax, anorm, lendsv-lsv, 1, e( lsv ),
614  $ n, info )
615  ELSEIF( iscale.EQ.2 ) THEN
616  CALL slascl( 'G', 0, 0, ssfmin, anorm, lendsv-lsv+1, 1,
617  $ d( lsv ), n, info )
618  CALL slascl( '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 cswap( 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 SSTEQR2
661 *
662  END
max
#define max(A, B)
Definition: pcgemr.c:180
csteqr2
subroutine csteqr2(COMPZ, N, D, E, Z, LDZ, NR, WORK, INFO)
Definition: csteqr2.f:2