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