ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
ssteqr2.f
Go to the documentation of this file.
1  SUBROUTINE ssteqr2( COMPZ, N, D, E, Z, LDZ, NR, WORK, INFO )
2 *
3 * -- LAPACK routine (version 2.0) --
4 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
5 * Courant Institute, Argonne National Lab, and Rice University
6 * September 30, 1994
7 *
8 * .. Scalar Arguments ..
9  CHARACTER COMPZ
10  INTEGER INFO, LDZ, N, NR
11 * ..
12 * .. Array Arguments ..
13  REAL D( * ), E( * ), WORK( * ), Z( LDZ, * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * SSTEQR2 is a modified version of LAPACK routine SSTEQR.
20 * SSTEQR2 computes all eigenvalues and, optionally, eigenvectors of a
21 * symmetric tridiagonal matrix using the implicit QL or QR method.
22 * running SSTEQR2 to perform updates on a distributed matrix Q.
23 * Proper usage of SSTEQR2 can be gleaned from examination of ScaLAPACK's
24 * PSSYEV.
25 *
26 * Arguments
27 * =========
28 *
29 * COMPZ (input) CHARACTER*1
30 * = 'N': Compute eigenvalues only.
31 * = 'I': Compute eigenvalues and eigenvectors of the
32 * tridiagonal matrix. Z must be initialized to the
33 * identity matrix by PDLASET or DLASET prior to entering
34 * this subroutine.
35 *
36 * N (input) INTEGER
37 * The order of the matrix. N >= 0.
38 *
39 * D (input/output) REAL array, dimension (N)
40 * On entry, the diagonal elements of the tridiagonal matrix.
41 * On exit, if INFO = 0, the eigenvalues in ascending order.
42 *
43 * E (input/output) REAL array, dimension (N-1)
44 * On entry, the (n-1) subdiagonal elements of the tridiagonal
45 * matrix.
46 * On exit, E has been destroyed.
47 *
48 * Z (local input/local output) REAL array, global
49 * dimension (N, N), local dimension (LDZ, NR).
50 * On entry, if COMPZ = 'V', then Z contains the orthogonal
51 * matrix used in the reduction to tridiagonal form.
52 * On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
53 * orthonormal eigenvectors of the original symmetric matrix,
54 * and if COMPZ = 'I', Z contains the orthonormal eigenvectors
55 * of the symmetric tridiagonal matrix.
56 * If COMPZ = 'N', then Z is not referenced.
57 *
58 * LDZ (input) INTEGER
59 * The leading dimension of the array Z. LDZ >= 1, and if
60 * eigenvectors are desired, then LDZ >= max(1,N).
61 *
62 * NR (input) INTEGER
63 * NR = MAX(1, NUMROC( N, NB, MYPROW, 0, NPROCS ) ).
64 * If COMPZ = 'N', then NR is not referenced.
65 *
66 * WORK (workspace) REAL array, dimension (max(1,2*N-2))
67 * If COMPZ = 'N', then WORK is not referenced.
68 *
69 * INFO (output) INTEGER
70 * = 0: successful exit
71 * < 0: if INFO = -i, the i-th argument had an illegal value
72 * > 0: the algorithm has failed to find all the eigenvalues in
73 * a total of 30*N iterations; if INFO = i, then i
74 * elements of E have not converged to zero; on exit, D
75 * and E contain the elements of a symmetric tridiagonal
76 * matrix which is orthogonally similar to the original
77 * matrix.
78 *
79 * =====================================================================
80 *
81 * .. Parameters ..
82  REAL ZERO, ONE, TWO, THREE
83  parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0,
84  $ three = 3.0e0 )
85  INTEGER MAXIT
86  parameter( maxit = 30 )
87 * ..
88 * .. Local Scalars ..
89  INTEGER I, ICOMPZ, II, ISCALE, J, JTOT, K, L, L1, LEND,
90  $ LENDM1, LENDP1, LENDSV, LM1, LSV, M, MM, MM1,
91  $ NM1, NMAXIT
92  REAL ANORM, B, C, EPS, EPS2, F, G, P, R, RT1, RT2,
93  $ S, SAFMAX, SAFMIN, SSFMAX, SSFMIN, TST
94 * ..
95 * .. External Functions ..
96  LOGICAL LSAME
97  REAL SLAMCH, SLANST, SLAPY2
98  EXTERNAL lsame, slamch, slanst, slapy2
99 * ..
100 * .. External Subroutines ..
101  EXTERNAL slae2, slaev2, slartg, slascl, slasr,
102  $ slasrt, sswap, xerbla
103 * ..
104 * .. Intrinsic Functions ..
105  INTRINSIC abs, max, sign, sqrt
106 * ..
107 * .. Executable Statements ..
108 *
109 * Test the input parameters.
110 *
111  info = 0
112 *
113  IF( lsame( compz, 'N' ) ) THEN
114  icompz = 0
115  ELSE IF( lsame( compz, 'I' ) ) THEN
116  icompz = 1
117  ELSE
118  icompz = -1
119  END IF
120  IF( icompz.LT.0 ) THEN
121  info = -1
122  ELSE IF( n.LT.0 ) THEN
123  info = -2
124  ELSE IF( icompz.GT.0 .AND. ldz.LT.max( 1, nr ) ) THEN
125  info = -6
126  END IF
127  IF( info.NE.0 ) THEN
128  CALL xerbla( 'SSTEQR2', -info )
129  RETURN
130  END IF
131 *
132 * Quick return if possible
133 *
134  IF( n.EQ.0 )
135  $ RETURN
136 *
137  IF( n.EQ.1 ) THEN
138  IF( icompz.EQ.1 )
139  $ z( 1, 1 ) = one
140  RETURN
141  END IF
142 *
143 * Determine the unit roundoff and over/underflow thresholds.
144 *
145  eps = slamch( 'E' )
146  eps2 = eps**2
147  safmin = slamch( 'S' )
148  safmax = one / safmin
149  ssfmax = sqrt( safmax ) / three
150  ssfmin = sqrt( safmin ) / eps2
151 *
152 * Compute the eigenvalues and eigenvectors of the tridiagonal
153 * matrix.
154 *
155  nmaxit = n*maxit
156  jtot = 0
157 *
158 * Determine where the matrix splits and choose QL or QR iteration
159 * for each block, according to whether top or bottom diagonal
160 * element is smaller.
161 *
162  l1 = 1
163  nm1 = n - 1
164 *
165  10 CONTINUE
166  IF( l1.GT.n )
167  $ GO TO 160
168  IF( l1.GT.1 )
169  $ e( l1-1 ) = zero
170  IF( l1.LE.nm1 ) THEN
171  DO 20 m = l1, nm1
172  tst = abs( e( m ) )
173  IF( tst.EQ.zero )
174  $ GO TO 30
175  IF( tst.LE.( sqrt( abs( d( m ) ) )*sqrt( abs( d( m+
176  $ 1 ) ) ) )*eps ) THEN
177  e( m ) = zero
178  GO TO 30
179  END IF
180  20 CONTINUE
181  END IF
182  m = n
183 *
184  30 CONTINUE
185  l = l1
186  lsv = l
187  lend = m
188  lendsv = lend
189  l1 = m + 1
190  IF( lend.EQ.l )
191  $ GO TO 10
192 *
193 * Scale submatrix in rows and columns L to LEND
194 *
195  anorm = slanst( 'I', lend-l+1, d( l ), e( l ) )
196  iscale = 0
197  IF( anorm.EQ.zero )
198  $ GO TO 10
199  IF( anorm.GT.ssfmax ) THEN
200  iscale = 1
201  CALL slascl( 'G', 0, 0, anorm, ssfmax, lend-l+1, 1, d( l ), n,
202  $ info )
203  CALL slascl( 'G', 0, 0, anorm, ssfmax, lend-l, 1, e( l ), n,
204  $ info )
205  ELSE IF( anorm.LT.ssfmin ) THEN
206  iscale = 2
207  CALL slascl( 'G', 0, 0, anorm, ssfmin, lend-l+1, 1, d( l ), n,
208  $ info )
209  CALL slascl( 'G', 0, 0, anorm, ssfmin, lend-l, 1, e( l ), n,
210  $ info )
211  END IF
212 *
213 * Choose between QL and QR iteration
214 *
215  IF( abs( d( lend ) ).LT.abs( d( l ) ) ) THEN
216  lend = lsv
217  l = lendsv
218  END IF
219 *
220  IF( lend.GT.l ) THEN
221 *
222 * QL Iteration
223 *
224 * Look for small subdiagonal element.
225 *
226  40 CONTINUE
227  IF( l.NE.lend ) THEN
228  lendm1 = lend - 1
229  DO 50 m = l, lendm1
230  tst = abs( e( m ) )**2
231  IF( tst.LE.( eps2*abs( d( m ) ) )*abs( d( m+1 ) )+
232  $ safmin )GO TO 60
233  50 CONTINUE
234  END IF
235 *
236  m = lend
237 *
238  60 CONTINUE
239  IF( m.LT.lend )
240  $ e( m ) = zero
241  p = d( l )
242  IF( m.EQ.l )
243  $ GO TO 80
244 *
245 * If remaining matrix is 2-by-2, use SLAE2 or SLAEV2
246 * to compute its eigensystem.
247 *
248  IF( m.EQ.l+1 ) THEN
249  IF( icompz.GT.0 ) THEN
250  CALL slaev2( d( l ), e( l ), d( l+1 ), rt1, rt2, c, s )
251  work( l ) = c
252  work( n-1+l ) = s
253  CALL slasr( 'R', 'V', 'B', nr, 2, work( l ),
254  $ work( n-1+l ), z( 1, l ), ldz )
255  ELSE
256  CALL slae2( d( l ), e( l ), d( l+1 ), rt1, rt2 )
257  END IF
258  d( l ) = rt1
259  d( l+1 ) = rt2
260  e( l ) = zero
261  l = l + 2
262  IF( l.LE.lend )
263  $ GO TO 40
264  GO TO 140
265  END IF
266 *
267  IF( jtot.EQ.nmaxit )
268  $ GO TO 140
269  jtot = jtot + 1
270 *
271 * Form shift.
272 *
273  g = ( d( l+1 )-p ) / ( two*e( l ) )
274  r = slapy2( g, one )
275  g = d( m ) - p + ( e( l ) / ( g+sign( r, g ) ) )
276 *
277  s = one
278  c = one
279  p = zero
280 *
281 * Inner loop
282 *
283  mm1 = m - 1
284  DO 70 i = mm1, l, -1
285  f = s*e( i )
286  b = c*e( i )
287  CALL slartg( g, f, c, s, r )
288  IF( i.NE.m-1 )
289  $ e( i+1 ) = r
290  g = d( i+1 ) - p
291  r = ( d( i )-g )*s + two*c*b
292  p = s*r
293  d( i+1 ) = g + p
294  g = c*r - b
295 *
296 * If eigenvectors are desired, then save rotations.
297 *
298  IF( icompz.GT.0 ) THEN
299  work( i ) = c
300  work( n-1+i ) = -s
301  END IF
302 *
303  70 CONTINUE
304 *
305 * If eigenvectors are desired, then apply saved rotations.
306 *
307  IF( icompz.GT.0 ) THEN
308  mm = m - l + 1
309  CALL slasr( 'R', 'V', 'B', nr, mm, work( l ), work( n-1+l ),
310  $ z( 1, l ), ldz )
311  END IF
312 *
313  d( l ) = d( l ) - p
314  e( l ) = g
315  GO TO 40
316 *
317 * Eigenvalue found.
318 *
319  80 CONTINUE
320  d( l ) = p
321 *
322  l = l + 1
323  IF( l.LE.lend )
324  $ GO TO 40
325  GO TO 140
326 *
327  ELSE
328 *
329 * QR Iteration
330 *
331 * Look for small superdiagonal element.
332 *
333  90 CONTINUE
334  IF( l.NE.lend ) THEN
335  lendp1 = lend + 1
336  DO 100 m = l, lendp1, -1
337  tst = abs( e( m-1 ) )**2
338  IF( tst.LE.( eps2*abs( d( m ) ) )*abs( d( m-1 ) )+
339  $ safmin )GO TO 110
340  100 CONTINUE
341  END IF
342 *
343  m = lend
344 *
345  110 CONTINUE
346  IF( m.GT.lend )
347  $ e( m-1 ) = zero
348  p = d( l )
349  IF( m.EQ.l )
350  $ GO TO 130
351 *
352 * If remaining matrix is 2-by-2, use SLAE2 or SLAEV2
353 * to compute its eigensystem.
354 *
355  IF( m.EQ.l-1 ) THEN
356  IF( icompz.GT.0 ) THEN
357  CALL slaev2( d( l-1 ), e( l-1 ), d( l ), rt1, rt2, c, s )
358  work( m ) = c
359  work( n-1+m ) = s
360  CALL slasr( 'R', 'V', 'F', nr, 2, work( m ),
361  $ work( n-1+m ), z( 1, l-1 ), ldz )
362  ELSE
363  CALL slae2( d( l-1 ), e( l-1 ), d( l ), rt1, rt2 )
364  END IF
365  d( l-1 ) = rt1
366  d( l ) = rt2
367  e( l-1 ) = zero
368  l = l - 2
369  IF( l.GE.lend )
370  $ GO TO 90
371  GO TO 140
372  END IF
373 *
374  IF( jtot.EQ.nmaxit )
375  $ GO TO 140
376  jtot = jtot + 1
377 *
378 * Form shift.
379 *
380  g = ( d( l-1 )-p ) / ( two*e( l-1 ) )
381  r = slapy2( g, one )
382  g = d( m ) - p + ( e( l-1 ) / ( g+sign( r, g ) ) )
383 *
384  s = one
385  c = one
386  p = zero
387 *
388 * Inner loop
389 *
390  lm1 = l - 1
391  DO 120 i = m, lm1
392  f = s*e( i )
393  b = c*e( i )
394  CALL slartg( g, f, c, s, r )
395  IF( i.NE.m )
396  $ e( i-1 ) = r
397  g = d( i ) - p
398  r = ( d( i+1 )-g )*s + two*c*b
399  p = s*r
400  d( i ) = g + p
401  g = c*r - b
402 *
403 * If eigenvectors are desired, then save rotations.
404 *
405  IF( icompz.GT.0 ) THEN
406  work( i ) = c
407  work( n-1+i ) = s
408  END IF
409 *
410  120 CONTINUE
411 *
412 * If eigenvectors are desired, then apply saved rotations.
413 *
414  IF( icompz.GT.0 ) THEN
415  mm = l - m + 1
416  CALL slasr( 'R', 'V', 'F', nr, mm, work( m ), work( n-1+m ),
417  $ z( 1, m ), ldz )
418  END IF
419 *
420  d( l ) = d( l ) - p
421  e( lm1 ) = g
422  GO TO 90
423 *
424 * Eigenvalue found.
425 *
426  130 CONTINUE
427  d( l ) = p
428 *
429  l = l - 1
430  IF( l.GE.lend )
431  $ GO TO 90
432  GO TO 140
433 *
434  END IF
435 *
436 * Undo scaling if necessary
437 *
438  140 CONTINUE
439  IF( iscale.EQ.1 ) THEN
440  CALL slascl( 'G', 0, 0, ssfmax, anorm, lendsv-lsv+1, 1,
441  $ d( lsv ), n, info )
442  CALL slascl( 'G', 0, 0, ssfmax, anorm, lendsv-lsv, 1, e( lsv ),
443  $ n, info )
444  ELSE IF( iscale.EQ.2 ) THEN
445  CALL slascl( 'G', 0, 0, ssfmin, anorm, lendsv-lsv+1, 1,
446  $ d( lsv ), n, info )
447  CALL slascl( 'G', 0, 0, ssfmin, anorm, lendsv-lsv, 1, e( lsv ),
448  $ n, info )
449  END IF
450 *
451 * Check for no convergence to an eigenvalue after a total
452 * of N*MAXIT iterations.
453 *
454  IF( jtot.LT.nmaxit )
455  $ GO TO 10
456  DO 150 i = 1, n - 1
457  IF( e( i ).NE.zero )
458  $ info = info + 1
459  150 CONTINUE
460  GO TO 190
461 *
462 * Order eigenvalues and eigenvectors.
463 *
464  160 CONTINUE
465  IF( icompz.EQ.0 ) THEN
466 *
467 * Use Quick Sort
468 *
469  CALL slasrt( 'I', n, d, info )
470 *
471  ELSE
472 *
473 * Use Selection Sort to minimize swaps of eigenvectors
474 *
475  DO 180 ii = 2, n
476  i = ii - 1
477  k = i
478  p = d( i )
479  DO 170 j = ii, n
480  IF( d( j ).LT.p ) THEN
481  k = j
482  p = d( j )
483  END IF
484  170 CONTINUE
485  IF( k.NE.i ) THEN
486  d( k ) = d( i )
487  d( i ) = p
488  CALL sswap( nr, z( 1, i ), 1, z( 1, k ), 1 )
489  END IF
490  180 CONTINUE
491  END IF
492 *
493  190 CONTINUE
494  RETURN
495 *
496 * End of SSTEQR2
497 *
498  END
ssteqr2
subroutine ssteqr2(COMPZ, N, D, E, Z, LDZ, NR, WORK, INFO)
Definition: ssteqr2.f:2
max
#define max(A, B)
Definition: pcgemr.c:180