LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dlasq2.f
Go to the documentation of this file.
1*> \brief \b DLASQ2 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 DLASQ2 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq2.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq2.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq2.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DLASQ2( N, Z, INFO )
22*
23* .. Scalar Arguments ..
24* INTEGER INFO, N
25* ..
26* .. Array Arguments ..
27* DOUBLE PRECISION Z( * )
28* ..
29*
30*
31*> \par Purpose:
32* =============
33*>
34*> \verbatim
35*>
36*> DLASQ2 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 : DLASQ2 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 DLASQ3.
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 DOUBLE PRECISION 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 auxOTHERcomputational
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 dlasq2( 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 DOUBLE PRECISION Z( * )
122* ..
123*
124* =====================================================================
125*
126* .. Parameters ..
127 DOUBLE PRECISION CBIAS
128 parameter( cbias = 1.50d0 )
129 DOUBLE PRECISION ZERO, HALF, ONE, TWO, FOUR, HUNDRD
130 parameter( zero = 0.0d0, half = 0.5d0, one = 1.0d0,
131 $ two = 2.0d0, four = 4.0d0, hundrd = 100.0d0 )
132* ..
133* .. Local Scalars ..
134 LOGICAL IEEE
135 INTEGER I0, I1, I4, IINFO, IPN4, ITER, IWHILA, IWHILB,
136 $ K, KMIN, N0, N1, NBIG, NDIV, NFAIL, PP, SPLT,
137 $ TTYPE
138 DOUBLE PRECISION 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 dlasq3, dlasrt, xerbla
145* ..
146* .. External Functions ..
147 INTEGER ILAENV
148 DOUBLE PRECISION DLAMCH
149 EXTERNAL dlamch, ilaenv
150* ..
151* .. Intrinsic Functions ..
152 INTRINSIC abs, dble, max, min, sqrt
153* ..
154* .. Executable Statements ..
155*
156* Test the input arguments.
157* (in case DLASQ2 is not called by DLASQ1)
158*
159 info = 0
160 eps = dlamch( 'Precision' )
161 safmin = dlamch( 'Safe minimum' )
162 tol = eps*hundrd
163 tol2 = tol**2
164*
165 IF( n.LT.0 ) THEN
166 info = -1
167 CALL xerbla( 'DLASQ2', 1 )
168 RETURN
169 ELSE IF( n.EQ.0 ) THEN
170 RETURN
171 ELSE IF( n.EQ.1 ) THEN
172*
173* 1-by-1 case.
174*
175 IF( z( 1 ).LT.zero ) THEN
176 info = -201
177 CALL xerbla( 'DLASQ2', 2 )
178 END IF
179 RETURN
180 ELSE IF( n.EQ.2 ) THEN
181*
182* 2-by-2 case.
183*
184 IF( z( 1 ).LT.zero ) THEN
185 info = -201
186 CALL xerbla( 'DLASQ2', 2 )
187 RETURN
188 ELSE IF( z( 2 ).LT.zero ) THEN
189 info = -202
190 CALL xerbla( 'DLASQ2', 2 )
191 RETURN
192 ELSE IF( z( 3 ).LT.zero ) THEN
193 info = -203
194 CALL xerbla( 'DLASQ2', 2 )
195 RETURN
196 ELSE IF( z( 3 ).GT.z( 1 ) ) THEN
197 d = z( 3 )
198 z( 3 ) = z( 1 )
199 z( 1 ) = d
200 END IF
201 z( 5 ) = z( 1 ) + z( 2 ) + z( 3 )
202 IF( z( 2 ).GT.z( 3 )*tol2 ) THEN
203 t = half*( ( z( 1 )-z( 3 ) )+z( 2 ) )
204 s = z( 3 )*( z( 2 ) / t )
205 IF( s.LE.t ) THEN
206 s = z( 3 )*( z( 2 ) / ( t*( one+sqrt( one+s / t ) ) ) )
207 ELSE
208 s = z( 3 )*( z( 2 ) / ( t+sqrt( t )*sqrt( t+s ) ) )
209 END IF
210 t = z( 1 ) + ( s+z( 2 ) )
211 z( 3 ) = z( 3 )*( z( 1 ) / t )
212 z( 1 ) = t
213 END IF
214 z( 2 ) = z( 3 )
215 z( 6 ) = z( 2 ) + z( 1 )
216 RETURN
217 END IF
218*
219* Check for negative data and compute sums of q's and e's.
220*
221 z( 2*n ) = zero
222 emin = z( 2 )
223 qmax = zero
224 zmax = zero
225 d = zero
226 e = zero
227*
228 DO 10 k = 1, 2*( n-1 ), 2
229 IF( z( k ).LT.zero ) THEN
230 info = -( 200+k )
231 CALL xerbla( 'DLASQ2', 2 )
232 RETURN
233 ELSE IF( z( k+1 ).LT.zero ) THEN
234 info = -( 200+k+1 )
235 CALL xerbla( 'DLASQ2', 2 )
236 RETURN
237 END IF
238 d = d + z( k )
239 e = e + z( k+1 )
240 qmax = max( qmax, z( k ) )
241 emin = min( emin, z( k+1 ) )
242 zmax = max( qmax, zmax, z( k+1 ) )
243 10 CONTINUE
244 IF( z( 2*n-1 ).LT.zero ) THEN
245 info = -( 200+2*n-1 )
246 CALL xerbla( 'DLASQ2', 2 )
247 RETURN
248 END IF
249 d = d + z( 2*n-1 )
250 qmax = max( qmax, z( 2*n-1 ) )
251 zmax = max( qmax, zmax )
252*
253* Check for diagonality.
254*
255 IF( e.EQ.zero ) THEN
256 DO 20 k = 2, n
257 z( k ) = z( 2*k-1 )
258 20 CONTINUE
259 CALL dlasrt( 'D', n, z, iinfo )
260 z( 2*n-1 ) = d
261 RETURN
262 END IF
263*
264 trace = d + e
265*
266* Check for zero data.
267*
268 IF( trace.EQ.zero ) THEN
269 z( 2*n-1 ) = zero
270 RETURN
271 END IF
272*
273* Check whether the machine is IEEE conformable.
274*
275 ieee = ( ilaenv( 10, 'DLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 )
276*
277* Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...).
278*
279 DO 30 k = 2*n, 2, -2
280 z( 2*k ) = zero
281 z( 2*k-1 ) = z( k )
282 z( 2*k-2 ) = zero
283 z( 2*k-3 ) = z( k-1 )
284 30 CONTINUE
285*
286 i0 = 1
287 n0 = n
288*
289* Reverse the qd-array, if warranted.
290*
291 IF( cbias*z( 4*i0-3 ).LT.z( 4*n0-3 ) ) THEN
292 ipn4 = 4*( i0+n0 )
293 DO 40 i4 = 4*i0, 2*( i0+n0-1 ), 4
294 temp = z( i4-3 )
295 z( i4-3 ) = z( ipn4-i4-3 )
296 z( ipn4-i4-3 ) = temp
297 temp = z( i4-1 )
298 z( i4-1 ) = z( ipn4-i4-5 )
299 z( ipn4-i4-5 ) = temp
300 40 CONTINUE
301 END IF
302*
303* Initial split checking via dqd and Li's test.
304*
305 pp = 0
306*
307 DO 80 k = 1, 2
308*
309 d = z( 4*n0+pp-3 )
310 DO 50 i4 = 4*( n0-1 ) + pp, 4*i0 + pp, -4
311 IF( z( i4-1 ).LE.tol2*d ) THEN
312 z( i4-1 ) = -zero
313 d = z( i4-3 )
314 ELSE
315 d = z( i4-3 )*( d / ( d+z( i4-1 ) ) )
316 END IF
317 50 CONTINUE
318*
319* dqd maps Z to ZZ plus Li's test.
320*
321 emin = z( 4*i0+pp+1 )
322 d = z( 4*i0+pp-3 )
323 DO 60 i4 = 4*i0 + pp, 4*( n0-1 ) + pp, 4
324 z( i4-2*pp-2 ) = d + z( i4-1 )
325 IF( z( i4-1 ).LE.tol2*d ) THEN
326 z( i4-1 ) = -zero
327 z( i4-2*pp-2 ) = d
328 z( i4-2*pp ) = zero
329 d = z( i4+1 )
330 ELSE IF( safmin*z( i4+1 ).LT.z( i4-2*pp-2 ) .AND.
331 $ safmin*z( i4-2*pp-2 ).LT.z( i4+1 ) ) THEN
332 temp = z( i4+1 ) / z( i4-2*pp-2 )
333 z( i4-2*pp ) = z( i4-1 )*temp
334 d = d*temp
335 ELSE
336 z( i4-2*pp ) = z( i4+1 )*( z( i4-1 ) / z( i4-2*pp-2 ) )
337 d = z( i4+1 )*( d / z( i4-2*pp-2 ) )
338 END IF
339 emin = min( emin, z( i4-2*pp ) )
340 60 CONTINUE
341 z( 4*n0-pp-2 ) = d
342*
343* Now find qmax.
344*
345 qmax = z( 4*i0-pp-2 )
346 DO 70 i4 = 4*i0 - pp + 2, 4*n0 - pp - 2, 4
347 qmax = max( qmax, z( i4 ) )
348 70 CONTINUE
349*
350* Prepare for the next iteration on K.
351*
352 pp = 1 - pp
353 80 CONTINUE
354*
355* Initialise variables to pass to DLASQ3.
356*
357 ttype = 0
358 dmin1 = zero
359 dmin2 = zero
360 dn = zero
361 dn1 = zero
362 dn2 = zero
363 g = zero
364 tau = zero
365*
366 iter = 2
367 nfail = 0
368 ndiv = 2*( n0-i0 )
369*
370 DO 160 iwhila = 1, n + 1
371 IF( n0.LT.1 )
372 $ GO TO 170
373*
374* While array unfinished do
375*
376* E(N0) holds the value of SIGMA when submatrix in I0:N0
377* splits from the rest of the array, but is negated.
378*
379 desig = zero
380 IF( n0.EQ.n ) THEN
381 sigma = zero
382 ELSE
383 sigma = -z( 4*n0-1 )
384 END IF
385 IF( sigma.LT.zero ) THEN
386 info = 1
387 RETURN
388 END IF
389*
390* Find last unreduced submatrix's top index I0, find QMAX and
391* EMIN. Find Gershgorin-type bound if Q's much greater than E's.
392*
393 emax = zero
394 IF( n0.GT.i0 ) THEN
395 emin = abs( z( 4*n0-5 ) )
396 ELSE
397 emin = zero
398 END IF
399 qmin = z( 4*n0-3 )
400 qmax = qmin
401 DO 90 i4 = 4*n0, 8, -4
402 IF( z( i4-5 ).LE.zero )
403 $ GO TO 100
404 IF( qmin.GE.four*emax ) THEN
405 qmin = min( qmin, z( i4-3 ) )
406 emax = max( emax, z( i4-5 ) )
407 END IF
408 qmax = max( qmax, z( i4-7 )+z( i4-5 ) )
409 emin = min( emin, z( i4-5 ) )
410 90 CONTINUE
411 i4 = 4
412*
413 100 CONTINUE
414 i0 = i4 / 4
415 pp = 0
416*
417 IF( n0-i0.GT.1 ) THEN
418 dee = z( 4*i0-3 )
419 deemin = dee
420 kmin = i0
421 DO 110 i4 = 4*i0+1, 4*n0-3, 4
422 dee = z( i4 )*( dee /( dee+z( i4-2 ) ) )
423 IF( dee.LE.deemin ) THEN
424 deemin = dee
425 kmin = ( i4+3 )/4
426 END IF
427 110 CONTINUE
428 IF( (kmin-i0)*2.LT.n0-kmin .AND.
429 $ deemin.LE.half*z(4*n0-3) ) THEN
430 ipn4 = 4*( i0+n0 )
431 pp = 2
432 DO 120 i4 = 4*i0, 2*( i0+n0-1 ), 4
433 temp = z( i4-3 )
434 z( i4-3 ) = z( ipn4-i4-3 )
435 z( ipn4-i4-3 ) = temp
436 temp = z( i4-2 )
437 z( i4-2 ) = z( ipn4-i4-2 )
438 z( ipn4-i4-2 ) = temp
439 temp = z( i4-1 )
440 z( i4-1 ) = z( ipn4-i4-5 )
441 z( ipn4-i4-5 ) = temp
442 temp = z( i4 )
443 z( i4 ) = z( ipn4-i4-4 )
444 z( ipn4-i4-4 ) = temp
445 120 CONTINUE
446 END IF
447 END IF
448*
449* Put -(initial shift) into DMIN.
450*
451 dmin = -max( zero, qmin-two*sqrt( qmin )*sqrt( emax ) )
452*
453* Now I0:N0 is unreduced.
454* PP = 0 for ping, PP = 1 for pong.
455* PP = 2 indicates that flipping was applied to the Z array and
456* and that the tests for deflation upon entry in DLASQ3
457* should not be performed.
458*
459 nbig = 100*( n0-i0+1 )
460 DO 140 iwhilb = 1, nbig
461 IF( i0.GT.n0 )
462 $ GO TO 150
463*
464* While submatrix unfinished take a good dqds step.
465*
466 CALL dlasq3( i0, n0, z, pp, dmin, sigma, desig, qmax, nfail,
467 $ iter, ndiv, ieee, ttype, dmin1, dmin2, dn, dn1,
468 $ dn2, g, tau )
469*
470 pp = 1 - pp
471*
472* When EMIN is very small check for splits.
473*
474 IF( pp.EQ.0 .AND. n0-i0.GE.3 ) THEN
475 IF( z( 4*n0 ).LE.tol2*qmax .OR.
476 $ z( 4*n0-1 ).LE.tol2*sigma ) THEN
477 splt = i0 - 1
478 qmax = z( 4*i0-3 )
479 emin = z( 4*i0-1 )
480 oldemn = z( 4*i0 )
481 DO 130 i4 = 4*i0, 4*( n0-3 ), 4
482 IF( z( i4 ).LE.tol2*z( i4-3 ) .OR.
483 $ z( i4-1 ).LE.tol2*sigma ) THEN
484 z( i4-1 ) = -sigma
485 splt = i4 / 4
486 qmax = zero
487 emin = z( i4+3 )
488 oldemn = z( i4+4 )
489 ELSE
490 qmax = max( qmax, z( i4+1 ) )
491 emin = min( emin, z( i4-1 ) )
492 oldemn = min( oldemn, z( i4 ) )
493 END IF
494 130 CONTINUE
495 z( 4*n0-1 ) = emin
496 z( 4*n0 ) = oldemn
497 i0 = splt + 1
498 END IF
499 END IF
500*
501 140 CONTINUE
502*
503 info = 2
504*
505* Maximum number of iterations exceeded, restore the shift
506* SIGMA and place the new d's and e's in a qd array.
507* This might need to be done for several blocks
508*
509 i1 = i0
510 n1 = n0
511 145 CONTINUE
512 tempq = z( 4*i0-3 )
513 z( 4*i0-3 ) = z( 4*i0-3 ) + sigma
514 DO k = i0+1, n0
515 tempe = z( 4*k-5 )
516 z( 4*k-5 ) = z( 4*k-5 ) * (tempq / z( 4*k-7 ))
517 tempq = z( 4*k-3 )
518 z( 4*k-3 ) = z( 4*k-3 ) + sigma + tempe - z( 4*k-5 )
519 END DO
520*
521* Prepare to do this on the previous block if there is one
522*
523 IF( i1.GT.1 ) THEN
524 n1 = i1-1
525 DO WHILE( ( i1.GE.2 ) .AND. ( z(4*i1-5).GE.zero ) )
526 i1 = i1 - 1
527 END DO
528 sigma = -z(4*n1-1)
529 GO TO 145
530 END IF
531
532 DO k = 1, n
533 z( 2*k-1 ) = z( 4*k-3 )
534*
535* Only the block 1..N0 is unfinished. The rest of the e's
536* must be essentially zero, although sometimes other data
537* has been stored in them.
538*
539 IF( k.LT.n0 ) THEN
540 z( 2*k ) = z( 4*k-1 )
541 ELSE
542 z( 2*k ) = 0
543 END IF
544 END DO
545 RETURN
546*
547* end IWHILB
548*
549 150 CONTINUE
550*
551 160 CONTINUE
552*
553 info = 3
554 RETURN
555*
556* end IWHILA
557*
558 170 CONTINUE
559*
560* Move q's to the front.
561*
562 DO 180 k = 2, n
563 z( k ) = z( 4*k-3 )
564 180 CONTINUE
565*
566* Sort and compute sum of eigenvalues.
567*
568 CALL dlasrt( 'D', n, z, iinfo )
569*
570 e = zero
571 DO 190 k = n, 1, -1
572 e = e + z( k )
573 190 CONTINUE
574*
575* Store trace, sum(eigenvalues) and information on performance.
576*
577 z( 2*n+1 ) = trace
578 z( 2*n+2 ) = e
579 z( 2*n+3 ) = dble( iter )
580 z( 2*n+4 ) = dble( ndiv ) / dble( n**2 )
581 z( 2*n+5 ) = hundrd*nfail / dble( iter )
582 RETURN
583*
584* End of DLASQ2
585*
586 END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dlasq3(I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL, ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1, DN2, G, TAU)
DLASQ3 checks for deflation, computes a shift and calls dqds. Used by sbdsqr.
Definition: dlasq3.f:182
subroutine dlasrt(ID, N, D, INFO)
DLASRT sorts numbers in increasing or decreasing order.
Definition: dlasrt.f:88
subroutine dlasq2(N, Z, INFO)
DLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated ...
Definition: dlasq2.f:112