LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
slasq2.f
Go to the documentation of this file.
1*> \brief \b SLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated with the qd Array Z to high relative accuracy. Used by sbdsqr and sstegr.
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download SLASQ2 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasq2.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasq2.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasq2.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SLASQ2( N, Z, INFO )
22*
23* .. Scalar Arguments ..
24* INTEGER INFO, N
25* ..
26* .. Array Arguments ..
27* REAL Z( * )
28* ..
29*
30*
31*> \par Purpose:
32* =============
33*>
34*> \verbatim
35*>
36*> SLASQ2 computes all the eigenvalues of the symmetric positive
37*> definite tridiagonal matrix associated with the qd array Z to high
38*> relative accuracy are computed to high relative accuracy, in the
39*> absence of denormalization, underflow and overflow.
40*>
41*> To see the relation of Z to the tridiagonal matrix, let L be a
42*> unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and
43*> let U be an upper bidiagonal matrix with 1's above and diagonal
44*> Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the
45*> symmetric tridiagonal to which it is similar.
46*>
47*> Note : SLASQ2 defines a logical variable, IEEE, which is true
48*> on machines which follow ieee-754 floating-point standard in their
49*> handling of infinities and NaNs, and false otherwise. This variable
50*> is passed to SLASQ3.
51*> \endverbatim
52*
53* Arguments:
54* ==========
55*
56*> \param[in] N
57*> \verbatim
58*> N is INTEGER
59*> The number of rows and columns in the matrix. N >= 0.
60*> \endverbatim
61*>
62*> \param[in,out] Z
63*> \verbatim
64*> Z is REAL array, dimension ( 4*N )
65*> On entry Z holds the qd array. On exit, entries 1 to N hold
66*> the eigenvalues in decreasing order, Z( 2*N+1 ) holds the
67*> trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If
68*> N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 )
69*> holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of
70*> shifts that failed.
71*> \endverbatim
72*>
73*> \param[out] INFO
74*> \verbatim
75*> INFO is INTEGER
76*> = 0: successful exit
77*> < 0: if the i-th argument is a scalar and had an illegal
78*> value, then INFO = -i, if the i-th argument is an
79*> array and the j-entry had an illegal value, then
80*> INFO = -(i*100+j)
81*> > 0: the algorithm failed
82*> = 1, a split was marked by a positive value in E
83*> = 2, current block of Z not diagonalized after 100*N
84*> iterations (in inner while loop). On exit Z holds
85*> a qd array with the same eigenvalues as the given Z.
86*> = 3, termination criterion of outer while loop not met
87*> (program created more than N unreduced blocks)
88*> \endverbatim
89*
90* Authors:
91* ========
92*
93*> \author Univ. of Tennessee
94*> \author Univ. of California Berkeley
95*> \author Univ. of Colorado Denver
96*> \author NAG Ltd.
97*
98*> \ingroup lasq2
99*
100*> \par Further Details:
101* =====================
102*>
103*> \verbatim
104*>
105*> Local Variables: I0:N0 defines a current unreduced segment of Z.
106*> The shifts are accumulated in SIGMA. Iteration count is in ITER.
107*> Ping-pong is controlled by PP (alternates between 0 and 1).
108*> \endverbatim
109*>
110* =====================================================================
111 SUBROUTINE slasq2( N, Z, INFO )
112*
113* -- LAPACK computational routine --
114* -- LAPACK is a software package provided by Univ. of Tennessee, --
115* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
116*
117* .. Scalar Arguments ..
118 INTEGER INFO, N
119* ..
120* .. Array Arguments ..
121 REAL Z( * )
122* ..
123*
124* =====================================================================
125*
126* .. Parameters ..
127 REAL CBIAS
128 parameter( cbias = 1.50e0 )
129 REAL ZERO, HALF, ONE, TWO, FOUR, HUNDRD
130 parameter( zero = 0.0e0, half = 0.5e0, one = 1.0e0,
131 $ two = 2.0e0, four = 4.0e0, hundrd = 100.0e0 )
132* ..
133* .. Local Scalars ..
134 LOGICAL IEEE
135 INTEGER I0, I4, IINFO, IPN4, ITER, IWHILA, IWHILB, K,
136 $ KMIN, N0, NBIG, NDIV, NFAIL, PP, SPLT, TTYPE,
137 $ I1, N1
138 REAL D, DEE, DEEMIN, DESIG, DMIN, DMIN1, DMIN2, DN,
139 $ DN1, DN2, E, EMAX, EMIN, EPS, G, OLDEMN, QMAX,
140 $ QMIN, S, SAFMIN, SIGMA, T, TAU, TEMP, TOL,
141 $ TOL2, TRACE, ZMAX, TEMPE, TEMPQ
142* ..
143* .. External Subroutines ..
144 EXTERNAL slasq3, slasrt, xerbla
145* ..
146* .. External Functions ..
147 REAL SLAMCH
148 EXTERNAL slamch
149* ..
150* .. Intrinsic Functions ..
151 INTRINSIC abs, max, min, real, sqrt
152* ..
153* .. Executable Statements ..
154*
155* Test the input arguments.
156* (in case SLASQ2 is not called by SLASQ1)
157*
158 info = 0
159 eps = slamch( 'Precision' )
160 safmin = slamch( 'Safe minimum' )
161 tol = eps*hundrd
162 tol2 = tol**2
163*
164 IF( n.LT.0 ) THEN
165 info = -1
166 CALL xerbla( 'SLASQ2', 1 )
167 RETURN
168 ELSE IF( n.EQ.0 ) THEN
169 RETURN
170 ELSE IF( n.EQ.1 ) THEN
171*
172* 1-by-1 case.
173*
174 IF( z( 1 ).LT.zero ) THEN
175 info = -201
176 CALL xerbla( 'SLASQ2', 2 )
177 END IF
178 RETURN
179 ELSE IF( n.EQ.2 ) THEN
180*
181* 2-by-2 case.
182*
183 IF( z( 1 ).LT.zero ) THEN
184 info = -201
185 CALL xerbla( 'SLASQ2', 2 )
186 RETURN
187 ELSE IF( z( 2 ).LT.zero ) THEN
188 info = -202
189 CALL xerbla( 'SLASQ2', 2 )
190 RETURN
191 ELSE IF( z( 3 ).LT.zero ) THEN
192 info = -203
193 CALL xerbla( 'SLASQ2', 2 )
194 RETURN
195 ELSE IF( z( 3 ).GT.z( 1 ) ) THEN
196 d = z( 3 )
197 z( 3 ) = z( 1 )
198 z( 1 ) = d
199 END IF
200 z( 5 ) = z( 1 ) + z( 2 ) + z( 3 )
201 IF( z( 2 ).GT.z( 3 )*tol2 ) THEN
202 t = half*( ( z( 1 )-z( 3 ) )+z( 2 ) )
203 s = z( 3 )*( z( 2 ) / t )
204 IF( s.LE.t ) THEN
205 s = z( 3 )*( z( 2 ) / ( t*( one+sqrt( one+s / t ) ) ) )
206 ELSE
207 s = z( 3 )*( z( 2 ) / ( t+sqrt( t )*sqrt( t+s ) ) )
208 END IF
209 t = z( 1 ) + ( s+z( 2 ) )
210 z( 3 ) = z( 3 )*( z( 1 ) / t )
211 z( 1 ) = t
212 END IF
213 z( 2 ) = z( 3 )
214 z( 6 ) = z( 2 ) + z( 1 )
215 RETURN
216 END IF
217*
218* Check for negative data and compute sums of q's and e's.
219*
220 z( 2*n ) = zero
221 emin = z( 2 )
222 qmax = zero
223 zmax = zero
224 d = zero
225 e = zero
226*
227 DO 10 k = 1, 2*( n-1 ), 2
228 IF( z( k ).LT.zero ) THEN
229 info = -( 200+k )
230 CALL xerbla( 'SLASQ2', 2 )
231 RETURN
232 ELSE IF( z( k+1 ).LT.zero ) THEN
233 info = -( 200+k+1 )
234 CALL xerbla( 'SLASQ2', 2 )
235 RETURN
236 END IF
237 d = d + z( k )
238 e = e + z( k+1 )
239 qmax = max( qmax, z( k ) )
240 emin = min( emin, z( k+1 ) )
241 zmax = max( qmax, zmax, z( k+1 ) )
242 10 CONTINUE
243 IF( z( 2*n-1 ).LT.zero ) THEN
244 info = -( 200+2*n-1 )
245 CALL xerbla( 'SLASQ2', 2 )
246 RETURN
247 END IF
248 d = d + z( 2*n-1 )
249 qmax = max( qmax, z( 2*n-1 ) )
250 zmax = max( qmax, zmax )
251*
252* Check for diagonality.
253*
254 IF( e.EQ.zero ) THEN
255 DO 20 k = 2, n
256 z( k ) = z( 2*k-1 )
257 20 CONTINUE
258 CALL slasrt( 'D', n, z, iinfo )
259 z( 2*n-1 ) = d
260 RETURN
261 END IF
262*
263 trace = d + e
264*
265* Check for zero data.
266*
267 IF( trace.EQ.zero ) THEN
268 z( 2*n-1 ) = zero
269 RETURN
270 END IF
271*
272* Check whether the machine is IEEE conformable.
273*
274* IEEE = ( ILAENV( 10, 'SLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 )
275*
276* [11/15/2008] The case IEEE=.TRUE. has a problem in single precision with
277* some the test matrices of type 16. The double precision code is fine.
278*
279 ieee = .false.
280*
281* Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...).
282*
283 DO 30 k = 2*n, 2, -2
284 z( 2*k ) = zero
285 z( 2*k-1 ) = z( k )
286 z( 2*k-2 ) = zero
287 z( 2*k-3 ) = z( k-1 )
288 30 CONTINUE
289*
290 i0 = 1
291 n0 = n
292*
293* Reverse the qd-array, if warranted.
294*
295 IF( cbias*z( 4*i0-3 ).LT.z( 4*n0-3 ) ) THEN
296 ipn4 = 4*( i0+n0 )
297 DO 40 i4 = 4*i0, 2*( i0+n0-1 ), 4
298 temp = z( i4-3 )
299 z( i4-3 ) = z( ipn4-i4-3 )
300 z( ipn4-i4-3 ) = temp
301 temp = z( i4-1 )
302 z( i4-1 ) = z( ipn4-i4-5 )
303 z( ipn4-i4-5 ) = temp
304 40 CONTINUE
305 END IF
306*
307* Initial split checking via dqd and Li's test.
308*
309 pp = 0
310*
311 DO 80 k = 1, 2
312*
313 d = z( 4*n0+pp-3 )
314 DO 50 i4 = 4*( n0-1 ) + pp, 4*i0 + pp, -4
315 IF( z( i4-1 ).LE.tol2*d ) THEN
316 z( i4-1 ) = -zero
317 d = z( i4-3 )
318 ELSE
319 d = z( i4-3 )*( d / ( d+z( i4-1 ) ) )
320 END IF
321 50 CONTINUE
322*
323* dqd maps Z to ZZ plus Li's test.
324*
325 emin = z( 4*i0+pp+1 )
326 d = z( 4*i0+pp-3 )
327 DO 60 i4 = 4*i0 + pp, 4*( n0-1 ) + pp, 4
328 z( i4-2*pp-2 ) = d + z( i4-1 )
329 IF( z( i4-1 ).LE.tol2*d ) THEN
330 z( i4-1 ) = -zero
331 z( i4-2*pp-2 ) = d
332 z( i4-2*pp ) = zero
333 d = z( i4+1 )
334 ELSE IF( safmin*z( i4+1 ).LT.z( i4-2*pp-2 ) .AND.
335 $ safmin*z( i4-2*pp-2 ).LT.z( i4+1 ) ) THEN
336 temp = z( i4+1 ) / z( i4-2*pp-2 )
337 z( i4-2*pp ) = z( i4-1 )*temp
338 d = d*temp
339 ELSE
340 z( i4-2*pp ) = z( i4+1 )*( z( i4-1 ) / z( i4-2*pp-2 ) )
341 d = z( i4+1 )*( d / z( i4-2*pp-2 ) )
342 END IF
343 emin = min( emin, z( i4-2*pp ) )
344 60 CONTINUE
345 z( 4*n0-pp-2 ) = d
346*
347* Now find qmax.
348*
349 qmax = z( 4*i0-pp-2 )
350 DO 70 i4 = 4*i0 - pp + 2, 4*n0 - pp - 2, 4
351 qmax = max( qmax, z( i4 ) )
352 70 CONTINUE
353*
354* Prepare for the next iteration on K.
355*
356 pp = 1 - pp
357 80 CONTINUE
358*
359* Initialise variables to pass to SLASQ3.
360*
361 ttype = 0
362 dmin1 = zero
363 dmin2 = zero
364 dn = zero
365 dn1 = zero
366 dn2 = zero
367 g = zero
368 tau = zero
369*
370 iter = 2
371 nfail = 0
372 ndiv = 2*( n0-i0 )
373*
374 DO 160 iwhila = 1, n + 1
375 IF( n0.LT.1 )
376 $ GO TO 170
377*
378* While array unfinished do
379*
380* E(N0) holds the value of SIGMA when submatrix in I0:N0
381* splits from the rest of the array, but is negated.
382*
383 desig = zero
384 IF( n0.EQ.n ) THEN
385 sigma = zero
386 ELSE
387 sigma = -z( 4*n0-1 )
388 END IF
389 IF( sigma.LT.zero ) THEN
390 info = 1
391 RETURN
392 END IF
393*
394* Find last unreduced submatrix's top index I0, find QMAX and
395* EMIN. Find Gershgorin-type bound if Q's much greater than E's.
396*
397 emax = zero
398 IF( n0.GT.i0 ) THEN
399 emin = abs( z( 4*n0-5 ) )
400 ELSE
401 emin = zero
402 END IF
403 qmin = z( 4*n0-3 )
404 qmax = qmin
405 DO 90 i4 = 4*n0, 8, -4
406 IF( z( i4-5 ).LE.zero )
407 $ GO TO 100
408 IF( qmin.GE.four*emax ) THEN
409 qmin = min( qmin, z( i4-3 ) )
410 emax = max( emax, z( i4-5 ) )
411 END IF
412 qmax = max( qmax, z( i4-7 )+z( i4-5 ) )
413 emin = min( emin, z( i4-5 ) )
414 90 CONTINUE
415 i4 = 4
416*
417 100 CONTINUE
418 i0 = i4 / 4
419 pp = 0
420*
421 IF( n0-i0.GT.1 ) THEN
422 dee = z( 4*i0-3 )
423 deemin = dee
424 kmin = i0
425 DO 110 i4 = 4*i0+1, 4*n0-3, 4
426 dee = z( i4 )*( dee /( dee+z( i4-2 ) ) )
427 IF( dee.LE.deemin ) THEN
428 deemin = dee
429 kmin = ( i4+3 )/4
430 END IF
431 110 CONTINUE
432 IF( (kmin-i0)*2.LT.n0-kmin .AND.
433 $ deemin.LE.half*z(4*n0-3) ) THEN
434 ipn4 = 4*( i0+n0 )
435 pp = 2
436 DO 120 i4 = 4*i0, 2*( i0+n0-1 ), 4
437 temp = z( i4-3 )
438 z( i4-3 ) = z( ipn4-i4-3 )
439 z( ipn4-i4-3 ) = temp
440 temp = z( i4-2 )
441 z( i4-2 ) = z( ipn4-i4-2 )
442 z( ipn4-i4-2 ) = temp
443 temp = z( i4-1 )
444 z( i4-1 ) = z( ipn4-i4-5 )
445 z( ipn4-i4-5 ) = temp
446 temp = z( i4 )
447 z( i4 ) = z( ipn4-i4-4 )
448 z( ipn4-i4-4 ) = temp
449 120 CONTINUE
450 END IF
451 END IF
452*
453* Put -(initial shift) into DMIN.
454*
455 dmin = -max( zero, qmin-two*sqrt( qmin )*sqrt( emax ) )
456*
457* Now I0:N0 is unreduced.
458* PP = 0 for ping, PP = 1 for pong.
459* PP = 2 indicates that flipping was applied to the Z array and
460* and that the tests for deflation upon entry in SLASQ3
461* should not be performed.
462*
463 nbig = 100*( n0-i0+1 )
464 DO 140 iwhilb = 1, nbig
465 IF( i0.GT.n0 )
466 $ GO TO 150
467*
468* While submatrix unfinished take a good dqds step.
469*
470 CALL slasq3( i0, n0, z, pp, dmin, sigma, desig, qmax, nfail,
471 $ iter, ndiv, ieee, ttype, dmin1, dmin2, dn, dn1,
472 $ dn2, g, tau )
473*
474 pp = 1 - pp
475*
476* When EMIN is very small check for splits.
477*
478 IF( pp.EQ.0 .AND. n0-i0.GE.3 ) THEN
479 IF( z( 4*n0 ).LE.tol2*qmax .OR.
480 $ z( 4*n0-1 ).LE.tol2*sigma ) THEN
481 splt = i0 - 1
482 qmax = z( 4*i0-3 )
483 emin = z( 4*i0-1 )
484 oldemn = z( 4*i0 )
485 DO 130 i4 = 4*i0, 4*( n0-3 ), 4
486 IF( z( i4 ).LE.tol2*z( i4-3 ) .OR.
487 $ z( i4-1 ).LE.tol2*sigma ) THEN
488 z( i4-1 ) = -sigma
489 splt = i4 / 4
490 qmax = zero
491 emin = z( i4+3 )
492 oldemn = z( i4+4 )
493 ELSE
494 qmax = max( qmax, z( i4+1 ) )
495 emin = min( emin, z( i4-1 ) )
496 oldemn = min( oldemn, z( i4 ) )
497 END IF
498 130 CONTINUE
499 z( 4*n0-1 ) = emin
500 z( 4*n0 ) = oldemn
501 i0 = splt + 1
502 END IF
503 END IF
504*
505 140 CONTINUE
506*
507 info = 2
508*
509* Maximum number of iterations exceeded, restore the shift
510* SIGMA and place the new d's and e's in a qd array.
511* This might need to be done for several blocks
512*
513 i1 = i0
514 n1 = n0
515 145 CONTINUE
516 tempq = z( 4*i0-3 )
517 z( 4*i0-3 ) = z( 4*i0-3 ) + sigma
518 DO k = i0+1, n0
519 tempe = z( 4*k-5 )
520 z( 4*k-5 ) = z( 4*k-5 ) * (tempq / z( 4*k-7 ))
521 tempq = z( 4*k-3 )
522 z( 4*k-3 ) = z( 4*k-3 ) + sigma + tempe - z( 4*k-5 )
523 END DO
524*
525* Prepare to do this on the previous block if there is one
526*
527 IF( i1.GT.1 ) THEN
528 n1 = i1-1
529 DO WHILE( ( i1.GE.2 ) .AND. ( z(4*i1-5).GE.zero ) )
530 i1 = i1 - 1
531 END DO
532 IF( i1.GE.1 ) THEN
533 sigma = -z(4*n1-1)
534 GO TO 145
535 END IF
536 END IF
537
538 DO k = 1, n
539 z( 2*k-1 ) = z( 4*k-3 )
540*
541* Only the block 1..N0 is unfinished. The rest of the e's
542* must be essentially zero, although sometimes other data
543* has been stored in them.
544*
545 IF( k.LT.n0 ) THEN
546 z( 2*k ) = z( 4*k-1 )
547 ELSE
548 z( 2*k ) = 0
549 END IF
550 END DO
551 RETURN
552*
553* end IWHILB
554*
555 150 CONTINUE
556*
557 160 CONTINUE
558*
559 info = 3
560 RETURN
561*
562* end IWHILA
563*
564 170 CONTINUE
565*
566* Move q's to the front.
567*
568 DO 180 k = 2, n
569 z( k ) = z( 4*k-3 )
570 180 CONTINUE
571*
572* Sort and compute sum of eigenvalues.
573*
574 CALL slasrt( 'D', n, z, iinfo )
575*
576 e = zero
577 DO 190 k = n, 1, -1
578 e = e + z( k )
579 190 CONTINUE
580*
581* Store trace, sum(eigenvalues) and information on performance.
582*
583 z( 2*n+1 ) = trace
584 z( 2*n+2 ) = e
585 z( 2*n+3 ) = real( iter )
586 z( 2*n+4 ) = real( ndiv ) / real( n**2 )
587 z( 2*n+5 ) = hundrd*nfail / real( iter )
588 RETURN
589*
590* End of SLASQ2
591*
592 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine slasq2(n, z, info)
SLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated ...
Definition slasq2.f:112
subroutine slasq3(i0, n0, z, pp, dmin, sigma, desig, qmax, nfail, iter, ndiv, ieee, ttype, dmin1, dmin2, dn, dn1, dn2, g, tau)
SLASQ3 checks for deflation, computes a shift and calls dqds. Used by sbdsqr.
Definition slasq3.f:182
subroutine slasrt(id, n, d, info)
SLASRT sorts numbers in increasing or decreasing order.
Definition slasrt.f:88