1: *> \brief \b DSTERF
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DSTERF + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsterf.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsterf.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsterf.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DSTERF( N, D, E, INFO )
22: *
23: * .. Scalar Arguments ..
24: * INTEGER INFO, N
25: * ..
26: * .. Array Arguments ..
27: * DOUBLE PRECISION D( * ), E( * )
28: * ..
29: *
30: *
31: *> \par Purpose:
32: * =============
33: *>
34: *> \verbatim
35: *>
36: *> DSTERF computes all eigenvalues of a symmetric tridiagonal matrix
37: *> using the Pal-Walker-Kahan variant of the QL or QR algorithm.
38: *> \endverbatim
39: *
40: * Arguments:
41: * ==========
42: *
43: *> \param[in] N
44: *> \verbatim
45: *> N is INTEGER
46: *> The order of the matrix. N >= 0.
47: *> \endverbatim
48: *>
49: *> \param[in,out] D
50: *> \verbatim
51: *> D is DOUBLE PRECISION array, dimension (N)
52: *> On entry, the n diagonal elements of the tridiagonal matrix.
53: *> On exit, if INFO = 0, the eigenvalues in ascending order.
54: *> \endverbatim
55: *>
56: *> \param[in,out] E
57: *> \verbatim
58: *> E is DOUBLE PRECISION array, dimension (N-1)
59: *> On entry, the (n-1) subdiagonal elements of the tridiagonal
60: *> matrix.
61: *> On exit, E has been destroyed.
62: *> \endverbatim
63: *>
64: *> \param[out] INFO
65: *> \verbatim
66: *> INFO is INTEGER
67: *> = 0: successful exit
68: *> < 0: if INFO = -i, the i-th argument had an illegal value
69: *> > 0: the algorithm failed to find all of the eigenvalues in
70: *> a total of 30*N iterations; if INFO = i, then i
71: *> elements of E have not converged to zero.
72: *> \endverbatim
73: *
74: * Authors:
75: * ========
76: *
77: *> \author Univ. of Tennessee
78: *> \author Univ. of California Berkeley
79: *> \author Univ. of Colorado Denver
80: *> \author NAG Ltd.
81: *
82: *> \ingroup auxOTHERcomputational
83: *
84: * =====================================================================
85: SUBROUTINE DSTERF( N, D, E, INFO )
86: *
87: * -- LAPACK computational routine --
88: * -- LAPACK is a software package provided by Univ. of Tennessee, --
89: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
90: *
91: * .. Scalar Arguments ..
92: INTEGER INFO, N
93: * ..
94: * .. Array Arguments ..
95: DOUBLE PRECISION D( * ), E( * )
96: * ..
97: *
98: * =====================================================================
99: *
100: * .. Parameters ..
101: DOUBLE PRECISION ZERO, ONE, TWO, THREE
102: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
103: $ THREE = 3.0D0 )
104: INTEGER MAXIT
105: PARAMETER ( MAXIT = 30 )
106: * ..
107: * .. Local Scalars ..
108: INTEGER I, ISCALE, JTOT, L, L1, LEND, LENDSV, LSV, M,
109: $ NMAXIT
110: DOUBLE PRECISION ALPHA, ANORM, BB, C, EPS, EPS2, GAMMA, OLDC,
111: $ OLDGAM, P, R, RT1, RT2, RTE, S, SAFMAX, SAFMIN,
112: $ SIGMA, SSFMAX, SSFMIN, RMAX
113: * ..
114: * .. External Functions ..
115: DOUBLE PRECISION DLAMCH, DLANST, DLAPY2
116: EXTERNAL DLAMCH, DLANST, DLAPY2
117: * ..
118: * .. External Subroutines ..
119: EXTERNAL DLAE2, DLASCL, DLASRT, XERBLA
120: * ..
121: * .. Intrinsic Functions ..
122: INTRINSIC ABS, SIGN, SQRT
123: * ..
124: * .. Executable Statements ..
125: *
126: * Test the input parameters.
127: *
128: INFO = 0
129: *
130: * Quick return if possible
131: *
132: IF( N.LT.0 ) THEN
133: INFO = -1
134: CALL XERBLA( 'DSTERF', -INFO )
135: RETURN
136: END IF
137: IF( N.LE.1 )
138: $ RETURN
139: *
140: * Determine the unit roundoff for this environment.
141: *
142: EPS = DLAMCH( 'E' )
143: EPS2 = EPS**2
144: SAFMIN = DLAMCH( 'S' )
145: SAFMAX = ONE / SAFMIN
146: SSFMAX = SQRT( SAFMAX ) / THREE
147: SSFMIN = SQRT( SAFMIN ) / EPS2
148: RMAX = DLAMCH( 'O' )
149: *
150: * Compute the eigenvalues of the tridiagonal matrix.
151: *
152: NMAXIT = N*MAXIT
153: SIGMA = ZERO
154: JTOT = 0
155: *
156: * Determine where the matrix splits and choose QL or QR iteration
157: * for each block, according to whether top or bottom diagonal
158: * element is smaller.
159: *
160: L1 = 1
161: *
162: 10 CONTINUE
163: IF( L1.GT.N )
164: $ GO TO 170
165: IF( L1.GT.1 )
166: $ E( L1-1 ) = ZERO
167: DO 20 M = L1, N - 1
168: IF( ABS( E( M ) ).LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+
169: $ 1 ) ) ) )*EPS ) THEN
170: E( M ) = ZERO
171: GO TO 30
172: END IF
173: 20 CONTINUE
174: M = N
175: *
176: 30 CONTINUE
177: L = L1
178: LSV = L
179: LEND = M
180: LENDSV = LEND
181: L1 = M + 1
182: IF( LEND.EQ.L )
183: $ GO TO 10
184: *
185: * Scale submatrix in rows and columns L to LEND
186: *
187: ANORM = DLANST( 'M', LEND-L+1, D( L ), E( L ) )
188: ISCALE = 0
189: IF( ANORM.EQ.ZERO )
190: $ GO TO 10
191: IF( (ANORM.GT.SSFMAX) ) THEN
192: ISCALE = 1
193: CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N,
194: $ INFO )
195: CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N,
196: $ INFO )
197: ELSE IF( ANORM.LT.SSFMIN ) THEN
198: ISCALE = 2
199: CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N,
200: $ INFO )
201: CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N,
202: $ INFO )
203: END IF
204: *
205: DO 40 I = L, LEND - 1
206: E( I ) = E( I )**2
207: 40 CONTINUE
208: *
209: * Choose between QL and QR iteration
210: *
211: IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN
212: LEND = LSV
213: L = LENDSV
214: END IF
215: *
216: IF( LEND.GE.L ) THEN
217: *
218: * QL Iteration
219: *
220: * Look for small subdiagonal element.
221: *
222: 50 CONTINUE
223: IF( L.NE.LEND ) THEN
224: DO 60 M = L, LEND - 1
225: IF( ABS( E( M ) ).LE.EPS2*ABS( D( M )*D( M+1 ) ) )
226: $ GO TO 70
227: 60 CONTINUE
228: END IF
229: M = LEND
230: *
231: 70 CONTINUE
232: IF( M.LT.LEND )
233: $ E( M ) = ZERO
234: P = D( L )
235: IF( M.EQ.L )
236: $ GO TO 90
237: *
238: * If remaining matrix is 2 by 2, use DLAE2 to compute its
239: * eigenvalues.
240: *
241: IF( M.EQ.L+1 ) THEN
242: RTE = SQRT( E( L ) )
243: CALL DLAE2( D( L ), RTE, D( L+1 ), RT1, RT2 )
244: D( L ) = RT1
245: D( L+1 ) = RT2
246: E( L ) = ZERO
247: L = L + 2
248: IF( L.LE.LEND )
249: $ GO TO 50
250: GO TO 150
251: END IF
252: *
253: IF( JTOT.EQ.NMAXIT )
254: $ GO TO 150
255: JTOT = JTOT + 1
256: *
257: * Form shift.
258: *
259: RTE = SQRT( E( L ) )
260: SIGMA = ( D( L+1 )-P ) / ( TWO*RTE )
261: R = DLAPY2( SIGMA, ONE )
262: SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )
263: *
264: C = ONE
265: S = ZERO
266: GAMMA = D( M ) - SIGMA
267: P = GAMMA*GAMMA
268: *
269: * Inner loop
270: *
271: DO 80 I = M - 1, L, -1
272: BB = E( I )
273: R = P + BB
274: IF( I.NE.M-1 )
275: $ E( I+1 ) = S*R
276: OLDC = C
277: C = P / R
278: S = BB / R
279: OLDGAM = GAMMA
280: ALPHA = D( I )
281: GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM
282: D( I+1 ) = OLDGAM + ( ALPHA-GAMMA )
283: IF( C.NE.ZERO ) THEN
284: P = ( GAMMA*GAMMA ) / C
285: ELSE
286: P = OLDC*BB
287: END IF
288: 80 CONTINUE
289: *
290: E( L ) = S*P
291: D( L ) = SIGMA + GAMMA
292: GO TO 50
293: *
294: * Eigenvalue found.
295: *
296: 90 CONTINUE
297: D( L ) = P
298: *
299: L = L + 1
300: IF( L.LE.LEND )
301: $ GO TO 50
302: GO TO 150
303: *
304: ELSE
305: *
306: * QR Iteration
307: *
308: * Look for small superdiagonal element.
309: *
310: 100 CONTINUE
311: DO 110 M = L, LEND + 1, -1
312: IF( ABS( E( M-1 ) ).LE.EPS2*ABS( D( M )*D( M-1 ) ) )
313: $ GO TO 120
314: 110 CONTINUE
315: M = LEND
316: *
317: 120 CONTINUE
318: IF( M.GT.LEND )
319: $ E( M-1 ) = ZERO
320: P = D( L )
321: IF( M.EQ.L )
322: $ GO TO 140
323: *
324: * If remaining matrix is 2 by 2, use DLAE2 to compute its
325: * eigenvalues.
326: *
327: IF( M.EQ.L-1 ) THEN
328: RTE = SQRT( E( L-1 ) )
329: CALL DLAE2( D( L ), RTE, D( L-1 ), RT1, RT2 )
330: D( L ) = RT1
331: D( L-1 ) = RT2
332: E( L-1 ) = ZERO
333: L = L - 2
334: IF( L.GE.LEND )
335: $ GO TO 100
336: GO TO 150
337: END IF
338: *
339: IF( JTOT.EQ.NMAXIT )
340: $ GO TO 150
341: JTOT = JTOT + 1
342: *
343: * Form shift.
344: *
345: RTE = SQRT( E( L-1 ) )
346: SIGMA = ( D( L-1 )-P ) / ( TWO*RTE )
347: R = DLAPY2( SIGMA, ONE )
348: SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )
349: *
350: C = ONE
351: S = ZERO
352: GAMMA = D( M ) - SIGMA
353: P = GAMMA*GAMMA
354: *
355: * Inner loop
356: *
357: DO 130 I = M, L - 1
358: BB = E( I )
359: R = P + BB
360: IF( I.NE.M )
361: $ E( I-1 ) = S*R
362: OLDC = C
363: C = P / R
364: S = BB / R
365: OLDGAM = GAMMA
366: ALPHA = D( I+1 )
367: GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM
368: D( I ) = OLDGAM + ( ALPHA-GAMMA )
369: IF( C.NE.ZERO ) THEN
370: P = ( GAMMA*GAMMA ) / C
371: ELSE
372: P = OLDC*BB
373: END IF
374: 130 CONTINUE
375: *
376: E( L-1 ) = S*P
377: D( L ) = SIGMA + GAMMA
378: GO TO 100
379: *
380: * Eigenvalue found.
381: *
382: 140 CONTINUE
383: D( L ) = P
384: *
385: L = L - 1
386: IF( L.GE.LEND )
387: $ GO TO 100
388: GO TO 150
389: *
390: END IF
391: *
392: * Undo scaling if necessary
393: *
394: 150 CONTINUE
395: IF( ISCALE.EQ.1 )
396: $ CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1,
397: $ D( LSV ), N, INFO )
398: IF( ISCALE.EQ.2 )
399: $ CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1,
400: $ D( LSV ), N, INFO )
401: *
402: * Check for no convergence to an eigenvalue after a total
403: * of N*MAXIT iterations.
404: *
405: IF( JTOT.LT.NMAXIT )
406: $ GO TO 10
407: DO 160 I = 1, N - 1
408: IF( E( I ).NE.ZERO )
409: $ INFO = INFO + 1
410: 160 CONTINUE
411: GO TO 180
412: *
413: * Sort eigenvalues in increasing order.
414: *
415: 170 CONTINUE
416: CALL DLASRT( 'I', N, D, INFO )
417: *
418: 180 CONTINUE
419: RETURN
420: *
421: * End of DSTERF
422: *
423: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsterf.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack