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: *> \date November 2011
83: *
84: *> \ingroup auxOTHERcomputational
85: *
86: * =====================================================================
87: SUBROUTINE DSTERF( N, D, E, INFO )
88: *
89: * -- LAPACK computational routine (version 3.4.0) --
90: * -- LAPACK is a software package provided by Univ. of Tennessee, --
91: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
92: * November 2011
93: *
94: * .. Scalar Arguments ..
95: INTEGER INFO, N
96: * ..
97: * .. Array Arguments ..
98: DOUBLE PRECISION D( * ), E( * )
99: * ..
100: *
101: * =====================================================================
102: *
103: * .. Parameters ..
104: DOUBLE PRECISION ZERO, ONE, TWO, THREE
105: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
106: $ THREE = 3.0D0 )
107: INTEGER MAXIT
108: PARAMETER ( MAXIT = 30 )
109: * ..
110: * .. Local Scalars ..
111: INTEGER I, ISCALE, JTOT, L, L1, LEND, LENDSV, LSV, M,
112: $ NMAXIT
113: DOUBLE PRECISION ALPHA, ANORM, BB, C, EPS, EPS2, GAMMA, OLDC,
114: $ OLDGAM, P, R, RT1, RT2, RTE, S, SAFMAX, SAFMIN,
115: $ SIGMA, SSFMAX, SSFMIN, RMAX
116: * ..
117: * .. External Functions ..
118: DOUBLE PRECISION DLAMCH, DLANST, DLAPY2
119: EXTERNAL DLAMCH, DLANST, DLAPY2
120: * ..
121: * .. External Subroutines ..
122: EXTERNAL DLAE2, DLASCL, DLASRT, XERBLA
123: * ..
124: * .. Intrinsic Functions ..
125: INTRINSIC ABS, SIGN, SQRT
126: * ..
127: * .. Executable Statements ..
128: *
129: * Test the input parameters.
130: *
131: INFO = 0
132: *
133: * Quick return if possible
134: *
135: IF( N.LT.0 ) THEN
136: INFO = -1
137: CALL XERBLA( 'DSTERF', -INFO )
138: RETURN
139: END IF
140: IF( N.LE.1 )
141: $ RETURN
142: *
143: * Determine the unit roundoff for this environment.
144: *
145: EPS = DLAMCH( 'E' )
146: EPS2 = EPS**2
147: SAFMIN = DLAMCH( 'S' )
148: SAFMAX = ONE / SAFMIN
149: SSFMAX = SQRT( SAFMAX ) / THREE
150: SSFMIN = SQRT( SAFMIN ) / EPS2
151: RMAX = DLAMCH( 'O' )
152: *
153: * Compute the eigenvalues of the tridiagonal matrix.
154: *
155: NMAXIT = N*MAXIT
156: SIGMA = ZERO
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: *
165: 10 CONTINUE
166: IF( L1.GT.N )
167: $ GO TO 170
168: IF( L1.GT.1 )
169: $ E( L1-1 ) = ZERO
170: DO 20 M = L1, N - 1
171: IF( ABS( E( M ) ).LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+
172: $ 1 ) ) ) )*EPS ) THEN
173: E( M ) = ZERO
174: GO TO 30
175: END IF
176: 20 CONTINUE
177: M = N
178: *
179: 30 CONTINUE
180: L = L1
181: LSV = L
182: LEND = M
183: LENDSV = LEND
184: L1 = M + 1
185: IF( LEND.EQ.L )
186: $ GO TO 10
187: *
188: * Scale submatrix in rows and columns L to LEND
189: *
190: ANORM = DLANST( 'M', LEND-L+1, D( L ), E( L ) )
191: ISCALE = 0
192: IF( ANORM.EQ.ZERO )
193: $ GO TO 10
194: IF( (ANORM.GT.SSFMAX) ) THEN
195: ISCALE = 1
196: CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N,
197: $ INFO )
198: CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N,
199: $ INFO )
200: ELSE IF( ANORM.LT.SSFMIN ) THEN
201: ISCALE = 2
202: CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N,
203: $ INFO )
204: CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N,
205: $ INFO )
206: END IF
207: *
208: DO 40 I = L, LEND - 1
209: E( I ) = E( I )**2
210: 40 CONTINUE
211: *
212: * Choose between QL and QR iteration
213: *
214: IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN
215: LEND = LSV
216: L = LENDSV
217: END IF
218: *
219: IF( LEND.GE.L ) THEN
220: *
221: * QL Iteration
222: *
223: * Look for small subdiagonal element.
224: *
225: 50 CONTINUE
226: IF( L.NE.LEND ) THEN
227: DO 60 M = L, LEND - 1
228: IF( ABS( E( M ) ).LE.EPS2*ABS( D( M )*D( M+1 ) ) )
229: $ GO TO 70
230: 60 CONTINUE
231: END IF
232: M = LEND
233: *
234: 70 CONTINUE
235: IF( M.LT.LEND )
236: $ E( M ) = ZERO
237: P = D( L )
238: IF( M.EQ.L )
239: $ GO TO 90
240: *
241: * If remaining matrix is 2 by 2, use DLAE2 to compute its
242: * eigenvalues.
243: *
244: IF( M.EQ.L+1 ) THEN
245: RTE = SQRT( E( L ) )
246: CALL DLAE2( D( L ), RTE, D( L+1 ), RT1, RT2 )
247: D( L ) = RT1
248: D( L+1 ) = RT2
249: E( L ) = ZERO
250: L = L + 2
251: IF( L.LE.LEND )
252: $ GO TO 50
253: GO TO 150
254: END IF
255: *
256: IF( JTOT.EQ.NMAXIT )
257: $ GO TO 150
258: JTOT = JTOT + 1
259: *
260: * Form shift.
261: *
262: RTE = SQRT( E( L ) )
263: SIGMA = ( D( L+1 )-P ) / ( TWO*RTE )
264: R = DLAPY2( SIGMA, ONE )
265: SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )
266: *
267: C = ONE
268: S = ZERO
269: GAMMA = D( M ) - SIGMA
270: P = GAMMA*GAMMA
271: *
272: * Inner loop
273: *
274: DO 80 I = M - 1, L, -1
275: BB = E( I )
276: R = P + BB
277: IF( I.NE.M-1 )
278: $ E( I+1 ) = S*R
279: OLDC = C
280: C = P / R
281: S = BB / R
282: OLDGAM = GAMMA
283: ALPHA = D( I )
284: GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM
285: D( I+1 ) = OLDGAM + ( ALPHA-GAMMA )
286: IF( C.NE.ZERO ) THEN
287: P = ( GAMMA*GAMMA ) / C
288: ELSE
289: P = OLDC*BB
290: END IF
291: 80 CONTINUE
292: *
293: E( L ) = S*P
294: D( L ) = SIGMA + GAMMA
295: GO TO 50
296: *
297: * Eigenvalue found.
298: *
299: 90 CONTINUE
300: D( L ) = P
301: *
302: L = L + 1
303: IF( L.LE.LEND )
304: $ GO TO 50
305: GO TO 150
306: *
307: ELSE
308: *
309: * QR Iteration
310: *
311: * Look for small superdiagonal element.
312: *
313: 100 CONTINUE
314: DO 110 M = L, LEND + 1, -1
315: IF( ABS( E( M-1 ) ).LE.EPS2*ABS( D( M )*D( M-1 ) ) )
316: $ GO TO 120
317: 110 CONTINUE
318: M = LEND
319: *
320: 120 CONTINUE
321: IF( M.GT.LEND )
322: $ E( M-1 ) = ZERO
323: P = D( L )
324: IF( M.EQ.L )
325: $ GO TO 140
326: *
327: * If remaining matrix is 2 by 2, use DLAE2 to compute its
328: * eigenvalues.
329: *
330: IF( M.EQ.L-1 ) THEN
331: RTE = SQRT( E( L-1 ) )
332: CALL DLAE2( D( L ), RTE, D( L-1 ), RT1, RT2 )
333: D( L ) = RT1
334: D( L-1 ) = RT2
335: E( L-1 ) = ZERO
336: L = L - 2
337: IF( L.GE.LEND )
338: $ GO TO 100
339: GO TO 150
340: END IF
341: *
342: IF( JTOT.EQ.NMAXIT )
343: $ GO TO 150
344: JTOT = JTOT + 1
345: *
346: * Form shift.
347: *
348: RTE = SQRT( E( L-1 ) )
349: SIGMA = ( D( L-1 )-P ) / ( TWO*RTE )
350: R = DLAPY2( SIGMA, ONE )
351: SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )
352: *
353: C = ONE
354: S = ZERO
355: GAMMA = D( M ) - SIGMA
356: P = GAMMA*GAMMA
357: *
358: * Inner loop
359: *
360: DO 130 I = M, L - 1
361: BB = E( I )
362: R = P + BB
363: IF( I.NE.M )
364: $ E( I-1 ) = S*R
365: OLDC = C
366: C = P / R
367: S = BB / R
368: OLDGAM = GAMMA
369: ALPHA = D( I+1 )
370: GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM
371: D( I ) = OLDGAM + ( ALPHA-GAMMA )
372: IF( C.NE.ZERO ) THEN
373: P = ( GAMMA*GAMMA ) / C
374: ELSE
375: P = OLDC*BB
376: END IF
377: 130 CONTINUE
378: *
379: E( L-1 ) = S*P
380: D( L ) = SIGMA + GAMMA
381: GO TO 100
382: *
383: * Eigenvalue found.
384: *
385: 140 CONTINUE
386: D( L ) = P
387: *
388: L = L - 1
389: IF( L.GE.LEND )
390: $ GO TO 100
391: GO TO 150
392: *
393: END IF
394: *
395: * Undo scaling if necessary
396: *
397: 150 CONTINUE
398: IF( ISCALE.EQ.1 )
399: $ CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1,
400: $ D( LSV ), N, INFO )
401: IF( ISCALE.EQ.2 )
402: $ CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1,
403: $ D( LSV ), N, INFO )
404: *
405: * Check for no convergence to an eigenvalue after a total
406: * of N*MAXIT iterations.
407: *
408: IF( JTOT.LT.NMAXIT )
409: $ GO TO 10
410: DO 160 I = 1, N - 1
411: IF( E( I ).NE.ZERO )
412: $ INFO = INFO + 1
413: 160 CONTINUE
414: GO TO 180
415: *
416: * Sort eigenvalues in increasing order.
417: *
418: 170 CONTINUE
419: CALL DLASRT( 'I', N, D, INFO )
420: *
421: 180 CONTINUE
422: RETURN
423: *
424: * End of DSTERF
425: *
426: 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