LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ dsterf()

 subroutine dsterf ( integer N, double precision, dimension( * ) D, double precision, dimension( * ) E, integer INFO )

DSTERF

Purpose:
``` DSTERF computes all eigenvalues of a symmetric tridiagonal matrix
using the Pal-Walker-Kahan variant of the QL or QR algorithm.```
Parameters
 [in] N ``` N is INTEGER The order of the matrix. N >= 0.``` [in,out] D ``` D is DOUBLE PRECISION array, dimension (N) On entry, the n diagonal elements of the tridiagonal matrix. On exit, if INFO = 0, the eigenvalues in ascending order.``` [in,out] E ``` E is DOUBLE PRECISION array, dimension (N-1) On entry, the (n-1) subdiagonal elements of the tridiagonal matrix. On exit, E has been destroyed.``` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: the algorithm failed to find all of the eigenvalues in a total of 30*N iterations; if INFO = i, then i elements of E have not converged to zero.```

Definition at line 85 of file dsterf.f.

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*
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
double precision function dlanst(NORM, N, D, E)
DLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition: dlanst.f:100
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: dlascl.f:143
subroutine dlae2(A, B, C, RT1, RT2)
DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix.
Definition: dlae2.f:102
double precision function dlapy2(X, Y)
DLAPY2 returns sqrt(x2+y2).
Definition: dlapy2.f:63
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dlasrt(ID, N, D, INFO)
DLASRT sorts numbers in increasing or decreasing order.
Definition: dlasrt.f:88
Here is the call graph for this function: