OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
dsteqr2.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine dsteqr2 (compz, n, d, e, z, ldz, nr, work, info)

Function/Subroutine Documentation

◆ dsteqr2()

subroutine dsteqr2 ( character compz,
integer n,
double precision, dimension( * ) d,
double precision, dimension( * ) e,
double precision, dimension( ldz, * ) z,
integer ldz,
integer nr,
double precision, dimension( * ) work,
integer info )

Definition at line 1 of file dsteqr2.f.

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*
dlae2
subroutine dlae2(a, b, c, rt1, rt2)
DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix.
Definition dlae2.f:102
dlartg
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition dlartg.f90:113
dlanst
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
dlapy2
double precision function dlapy2(x, y)
DLAPY2 returns sqrt(x2+y2).
Definition dlapy2.f:63
dlaev2
subroutine dlaev2(a, b, c, rt1, rt2, cs1, sn1)
DLAEV2 computes the eigenvalues and eigenvectors of a 2-by-2 symmetric/Hermitian matrix.
Definition dlaev2.f:120
dlascl
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
dlasr
subroutine dlasr(side, pivot, direct, m, n, c, s, a, lda)
DLASR applies a sequence of plane rotations to a general rectangular matrix.
Definition dlasr.f:199
dlasrt
subroutine dlasrt(id, n, d, info)
DLASRT sorts numbers in increasing or decreasing order.
Definition dlasrt.f:88
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
dswap
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
max
#define max(a, b)
Definition macros.h:21