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

Go to the source code of this file.

Functions/Subroutines

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

Function/Subroutine Documentation

◆ ssteqr2()

subroutine ssteqr2 ( character compz,
integer n,
real, dimension( * ) d,
real, dimension( * ) e,
real, dimension( ldz, * ) z,
integer ldz,
integer nr,
real, dimension( * ) work,
integer info )

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