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

◆ slasq2()

 subroutine slasq2 ( integer n, real, dimension( * ) z, integer info )

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.

Purpose:
``` SLASQ2 computes all the eigenvalues of the symmetric positive
definite tridiagonal matrix associated with the qd array Z to high
relative accuracy are computed to high relative accuracy, in the
absence of denormalization, underflow and overflow.

To see the relation of Z to the tridiagonal matrix, let L be a
unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and
let U be an upper bidiagonal matrix with 1's above and diagonal
Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the
symmetric tridiagonal to which it is similar.

Note : SLASQ2 defines a logical variable, IEEE, which is true
on machines which follow ieee-754 floating-point standard in their
handling of infinities and NaNs, and false otherwise. This variable
is passed to SLASQ3.```
Parameters
 [in] N ``` N is INTEGER The number of rows and columns in the matrix. N >= 0.``` [in,out] Z ``` Z is REAL array, dimension ( 4*N ) On entry Z holds the qd array. On exit, entries 1 to N hold the eigenvalues in decreasing order, Z( 2*N+1 ) holds the trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 ) holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of shifts that failed.``` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: if the i-th argument is a scalar and had an illegal value, then INFO = -i, if the i-th argument is an array and the j-entry had an illegal value, then INFO = -(i*100+j) > 0: the algorithm failed = 1, a split was marked by a positive value in E = 2, current block of Z not diagonalized after 100*N iterations (in inner while loop). On exit Z holds a qd array with the same eigenvalues as the given Z. = 3, termination criterion of outer while loop not met (program created more than N unreduced blocks)```
Further Details:
```  Local Variables: I0:N0 defines a current unreduced segment of Z.
The shifts are accumulated in SIGMA. Iteration count is in ITER.
Ping-pong is controlled by PP (alternates between 0 and 1).```

Definition at line 111 of file slasq2.f.

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*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
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
Here is the call graph for this function:
Here is the caller graph for this function: