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

◆ slaqr2()

 subroutine slaqr2 ( logical wantt, logical wantz, integer n, integer ktop, integer kbot, integer nw, real, dimension( ldh, * ) h, integer ldh, integer iloz, integer ihiz, real, dimension( ldz, * ) z, integer ldz, integer ns, integer nd, real, dimension( * ) sr, real, dimension( * ) si, real, dimension( ldv, * ) v, integer ldv, integer nh, real, dimension( ldt, * ) t, integer ldt, integer nv, real, dimension( ldwv, * ) wv, integer ldwv, real, dimension( * ) work, integer lwork )

SLAQR2 performs the orthogonal similarity transformation of a Hessenberg matrix to detect and deflate fully converged eigenvalues from a trailing principal submatrix (aggressive early deflation).

Purpose:
```    SLAQR2 is identical to SLAQR3 except that it avoids
recursion by calling SLAHQR instead of SLAQR4.

Aggressive early deflation:

This subroutine accepts as input an upper Hessenberg matrix
H and performs an orthogonal similarity transformation
designed to detect and deflate fully converged eigenvalues from
a trailing principal submatrix.  On output H has been over-
written by a new Hessenberg matrix that is a perturbation of
an orthogonal similarity transformation of H.  It is to be
hoped that the final version of H has many zero subdiagonal
entries.```
Parameters
 [in] WANTT ``` WANTT is LOGICAL If .TRUE., then the Hessenberg matrix H is fully updated so that the quasi-triangular Schur factor may be computed (in cooperation with the calling subroutine). If .FALSE., then only enough of H is updated to preserve the eigenvalues.``` [in] WANTZ ``` WANTZ is LOGICAL If .TRUE., then the orthogonal matrix Z is updated so so that the orthogonal Schur factor may be computed (in cooperation with the calling subroutine). If .FALSE., then Z is not referenced.``` [in] N ``` N is INTEGER The order of the matrix H and (if WANTZ is .TRUE.) the order of the orthogonal matrix Z.``` [in] KTOP ``` KTOP is INTEGER It is assumed that either KTOP = 1 or H(KTOP,KTOP-1)=0. KBOT and KTOP together determine an isolated block along the diagonal of the Hessenberg matrix.``` [in] KBOT ``` KBOT is INTEGER It is assumed without a check that either KBOT = N or H(KBOT+1,KBOT)=0. KBOT and KTOP together determine an isolated block along the diagonal of the Hessenberg matrix.``` [in] NW ``` NW is INTEGER Deflation window size. 1 <= NW <= (KBOT-KTOP+1).``` [in,out] H ``` H is REAL array, dimension (LDH,N) On input the initial N-by-N section of H stores the Hessenberg matrix undergoing aggressive early deflation. On output H has been transformed by an orthogonal similarity transformation, perturbed, and the returned to Hessenberg form that (it is to be hoped) has some zero subdiagonal entries.``` [in] LDH ``` LDH is INTEGER Leading dimension of H just as declared in the calling subroutine. N <= LDH``` [in] ILOZ ` ILOZ is INTEGER` [in] IHIZ ``` IHIZ is INTEGER Specify the rows of Z to which transformations must be applied if WANTZ is .TRUE.. 1 <= ILOZ <= IHIZ <= N.``` [in,out] Z ``` Z is REAL array, dimension (LDZ,N) IF WANTZ is .TRUE., then on output, the orthogonal similarity transformation mentioned above has been accumulated into Z(ILOZ:IHIZ,ILOZ:IHIZ) from the right. If WANTZ is .FALSE., then Z is unreferenced.``` [in] LDZ ``` LDZ is INTEGER The leading dimension of Z just as declared in the calling subroutine. 1 <= LDZ.``` [out] NS ``` NS is INTEGER The number of unconverged (ie approximate) eigenvalues returned in SR and SI that may be used as shifts by the calling subroutine.``` [out] ND ``` ND is INTEGER The number of converged eigenvalues uncovered by this subroutine.``` [out] SR ` SR is REAL array, dimension (KBOT)` [out] SI ``` SI is REAL array, dimension (KBOT) On output, the real and imaginary parts of approximate eigenvalues that may be used for shifts are stored in SR(KBOT-ND-NS+1) through SR(KBOT-ND) and SI(KBOT-ND-NS+1) through SI(KBOT-ND), respectively. The real and imaginary parts of converged eigenvalues are stored in SR(KBOT-ND+1) through SR(KBOT) and SI(KBOT-ND+1) through SI(KBOT), respectively.``` [out] V ``` V is REAL array, dimension (LDV,NW) An NW-by-NW work array.``` [in] LDV ``` LDV is INTEGER The leading dimension of V just as declared in the calling subroutine. NW <= LDV``` [in] NH ``` NH is INTEGER The number of columns of T. NH >= NW.``` [out] T ` T is REAL array, dimension (LDT,NW)` [in] LDT ``` LDT is INTEGER The leading dimension of T just as declared in the calling subroutine. NW <= LDT``` [in] NV ``` NV is INTEGER The number of rows of work array WV available for workspace. NV >= NW.``` [out] WV ` WV is REAL array, dimension (LDWV,NW)` [in] LDWV ``` LDWV is INTEGER The leading dimension of W just as declared in the calling subroutine. NW <= LDV``` [out] WORK ``` WORK is REAL array, dimension (LWORK) On exit, WORK(1) is set to an estimate of the optimal value of LWORK for the given values of N, NW, KTOP and KBOT.``` [in] LWORK ``` LWORK is INTEGER The dimension of the work array WORK. LWORK = 2*NW suffices, but greater efficiency may result from larger values of LWORK. If LWORK = -1, then a workspace query is assumed; SLAQR2 only estimates the optimal workspace size for the given values of N, NW, KTOP and KBOT. The estimate is returned in WORK(1). No error message related to LWORK is issued by XERBLA. Neither H nor Z are accessed.```
Contributors:
Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USA

Definition at line 275 of file slaqr2.f.

278*
279* -- LAPACK auxiliary routine --
280* -- LAPACK is a software package provided by Univ. of Tennessee, --
281* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
282*
283* .. Scalar Arguments ..
284 INTEGER IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV,
285 \$ LDZ, LWORK, N, ND, NH, NS, NV, NW
286 LOGICAL WANTT, WANTZ
287* ..
288* .. Array Arguments ..
289 REAL H( LDH, * ), SI( * ), SR( * ), T( LDT, * ),
290 \$ V( LDV, * ), WORK( * ), WV( LDWV, * ),
291 \$ Z( LDZ, * )
292* ..
293*
294* ================================================================
295* .. Parameters ..
296 REAL ZERO, ONE
297 parameter( zero = 0.0e0, one = 1.0e0 )
298* ..
299* .. Local Scalars ..
300 REAL AA, BB, BETA, CC, CS, DD, EVI, EVK, FOO, S,
301 \$ SAFMAX, SAFMIN, SMLNUM, SN, TAU, ULP
302 INTEGER I, IFST, ILST, INFO, INFQR, J, JW, K, KCOL,
303 \$ KEND, KLN, KROW, KWTOP, LTOP, LWK1, LWK2,
304 \$ LWKOPT
305 LOGICAL BULGE, SORTED
306* ..
307* .. External Functions ..
308 REAL SLAMCH, SROUNDUP_LWORK
309 EXTERNAL slamch, sroundup_lwork
310* ..
311* .. External Subroutines ..
312 EXTERNAL scopy, sgehrd, sgemm, slacpy, slahqr,
314* ..
315* .. Intrinsic Functions ..
316 INTRINSIC abs, int, max, min, real, sqrt
317* ..
318* .. Executable Statements ..
319*
320* ==== Estimate optimal workspace. ====
321*
322 jw = min( nw, kbot-ktop+1 )
323 IF( jw.LE.2 ) THEN
324 lwkopt = 1
325 ELSE
326*
327* ==== Workspace query call to SGEHRD ====
328*
329 CALL sgehrd( jw, 1, jw-1, t, ldt, work, work, -1, info )
330 lwk1 = int( work( 1 ) )
331*
332* ==== Workspace query call to SORMHR ====
333*
334 CALL sormhr( 'R', 'N', jw, jw, 1, jw-1, t, ldt, work, v, ldv,
335 \$ work, -1, info )
336 lwk2 = int( work( 1 ) )
337*
338* ==== Optimal workspace ====
339*
340 lwkopt = jw + max( lwk1, lwk2 )
341 END IF
342*
343* ==== Quick return in case of workspace query. ====
344*
345 IF( lwork.EQ.-1 ) THEN
346 work( 1 ) = sroundup_lwork( lwkopt )
347 RETURN
348 END IF
349*
350* ==== Nothing to do ...
351* ... for an empty active block ... ====
352 ns = 0
353 nd = 0
354 work( 1 ) = one
355 IF( ktop.GT.kbot )
356 \$ RETURN
357* ... nor for an empty deflation window. ====
358 IF( nw.LT.1 )
359 \$ RETURN
360*
361* ==== Machine constants ====
362*
363 safmin = slamch( 'SAFE MINIMUM' )
364 safmax = one / safmin
365 ulp = slamch( 'PRECISION' )
366 smlnum = safmin*( real( n ) / ulp )
367*
368* ==== Setup deflation window ====
369*
370 jw = min( nw, kbot-ktop+1 )
371 kwtop = kbot - jw + 1
372 IF( kwtop.EQ.ktop ) THEN
373 s = zero
374 ELSE
375 s = h( kwtop, kwtop-1 )
376 END IF
377*
378 IF( kbot.EQ.kwtop ) THEN
379*
380* ==== 1-by-1 deflation window: not much to do ====
381*
382 sr( kwtop ) = h( kwtop, kwtop )
383 si( kwtop ) = zero
384 ns = 1
385 nd = 0
386 IF( abs( s ).LE.max( smlnum, ulp*abs( h( kwtop, kwtop ) ) ) )
387 \$ THEN
388 ns = 0
389 nd = 1
390 IF( kwtop.GT.ktop )
391 \$ h( kwtop, kwtop-1 ) = zero
392 END IF
393 work( 1 ) = one
394 RETURN
395 END IF
396*
397* ==== Convert to spike-triangular form. (In case of a
398* . rare QR failure, this routine continues to do
399* . aggressive early deflation using that part of
400* . the deflation window that converged using INFQR
401* . here and there to keep track.) ====
402*
403 CALL slacpy( 'U', jw, jw, h( kwtop, kwtop ), ldh, t, ldt )
404 CALL scopy( jw-1, h( kwtop+1, kwtop ), ldh+1, t( 2, 1 ), ldt+1 )
405*
406 CALL slaset( 'A', jw, jw, zero, one, v, ldv )
407 CALL slahqr( .true., .true., jw, 1, jw, t, ldt, sr( kwtop ),
408 \$ si( kwtop ), 1, jw, v, ldv, infqr )
409*
410* ==== STREXC needs a clean margin near the diagonal ====
411*
412 DO 10 j = 1, jw - 3
413 t( j+2, j ) = zero
414 t( j+3, j ) = zero
415 10 CONTINUE
416 IF( jw.GT.2 )
417 \$ t( jw, jw-2 ) = zero
418*
419* ==== Deflation detection loop ====
420*
421 ns = jw
422 ilst = infqr + 1
423 20 CONTINUE
424 IF( ilst.LE.ns ) THEN
425 IF( ns.EQ.1 ) THEN
426 bulge = .false.
427 ELSE
428 bulge = t( ns, ns-1 ).NE.zero
429 END IF
430*
431* ==== Small spike tip test for deflation ====
432*
433 IF( .NOT.bulge ) THEN
434*
435* ==== Real eigenvalue ====
436*
437 foo = abs( t( ns, ns ) )
438 IF( foo.EQ.zero )
439 \$ foo = abs( s )
440 IF( abs( s*v( 1, ns ) ).LE.max( smlnum, ulp*foo ) ) THEN
441*
442* ==== Deflatable ====
443*
444 ns = ns - 1
445 ELSE
446*
447* ==== Undeflatable. Move it up out of the way.
448* . (STREXC can not fail in this case.) ====
449*
450 ifst = ns
451 CALL strexc( 'V', jw, t, ldt, v, ldv, ifst, ilst, work,
452 \$ info )
453 ilst = ilst + 1
454 END IF
455 ELSE
456*
457* ==== Complex conjugate pair ====
458*
459 foo = abs( t( ns, ns ) ) + sqrt( abs( t( ns, ns-1 ) ) )*
460 \$ sqrt( abs( t( ns-1, ns ) ) )
461 IF( foo.EQ.zero )
462 \$ foo = abs( s )
463 IF( max( abs( s*v( 1, ns ) ), abs( s*v( 1, ns-1 ) ) ).LE.
464 \$ max( smlnum, ulp*foo ) ) THEN
465*
466* ==== Deflatable ====
467*
468 ns = ns - 2
469 ELSE
470*
471* ==== Undeflatable. Move them up out of the way.
472* . Fortunately, STREXC does the right thing with
473* . ILST in case of a rare exchange failure. ====
474*
475 ifst = ns
476 CALL strexc( 'V', jw, t, ldt, v, ldv, ifst, ilst, work,
477 \$ info )
478 ilst = ilst + 2
479 END IF
480 END IF
481*
482* ==== End deflation detection loop ====
483*
484 GO TO 20
485 END IF
486*
488*
489 IF( ns.EQ.0 )
490 \$ s = zero
491*
492 IF( ns.LT.jw ) THEN
493*
494* ==== sorting diagonal blocks of T improves accuracy for
495* . graded matrices. Bubble sort deals well with
496* . exchange failures. ====
497*
498 sorted = .false.
499 i = ns + 1
500 30 CONTINUE
501 IF( sorted )
502 \$ GO TO 50
503 sorted = .true.
504*
505 kend = i - 1
506 i = infqr + 1
507 IF( i.EQ.ns ) THEN
508 k = i + 1
509 ELSE IF( t( i+1, i ).EQ.zero ) THEN
510 k = i + 1
511 ELSE
512 k = i + 2
513 END IF
514 40 CONTINUE
515 IF( k.LE.kend ) THEN
516 IF( k.EQ.i+1 ) THEN
517 evi = abs( t( i, i ) )
518 ELSE
519 evi = abs( t( i, i ) ) + sqrt( abs( t( i+1, i ) ) )*
520 \$ sqrt( abs( t( i, i+1 ) ) )
521 END IF
522*
523 IF( k.EQ.kend ) THEN
524 evk = abs( t( k, k ) )
525 ELSE IF( t( k+1, k ).EQ.zero ) THEN
526 evk = abs( t( k, k ) )
527 ELSE
528 evk = abs( t( k, k ) ) + sqrt( abs( t( k+1, k ) ) )*
529 \$ sqrt( abs( t( k, k+1 ) ) )
530 END IF
531*
532 IF( evi.GE.evk ) THEN
533 i = k
534 ELSE
535 sorted = .false.
536 ifst = i
537 ilst = k
538 CALL strexc( 'V', jw, t, ldt, v, ldv, ifst, ilst, work,
539 \$ info )
540 IF( info.EQ.0 ) THEN
541 i = ilst
542 ELSE
543 i = k
544 END IF
545 END IF
546 IF( i.EQ.kend ) THEN
547 k = i + 1
548 ELSE IF( t( i+1, i ).EQ.zero ) THEN
549 k = i + 1
550 ELSE
551 k = i + 2
552 END IF
553 GO TO 40
554 END IF
555 GO TO 30
556 50 CONTINUE
557 END IF
558*
559* ==== Restore shift/eigenvalue array from T ====
560*
561 i = jw
562 60 CONTINUE
563 IF( i.GE.infqr+1 ) THEN
564 IF( i.EQ.infqr+1 ) THEN
565 sr( kwtop+i-1 ) = t( i, i )
566 si( kwtop+i-1 ) = zero
567 i = i - 1
568 ELSE IF( t( i, i-1 ).EQ.zero ) THEN
569 sr( kwtop+i-1 ) = t( i, i )
570 si( kwtop+i-1 ) = zero
571 i = i - 1
572 ELSE
573 aa = t( i-1, i-1 )
574 cc = t( i, i-1 )
575 bb = t( i-1, i )
576 dd = t( i, i )
577 CALL slanv2( aa, bb, cc, dd, sr( kwtop+i-2 ),
578 \$ si( kwtop+i-2 ), sr( kwtop+i-1 ),
579 \$ si( kwtop+i-1 ), cs, sn )
580 i = i - 2
581 END IF
582 GO TO 60
583 END IF
584*
585 IF( ns.LT.jw .OR. s.EQ.zero ) THEN
586 IF( ns.GT.1 .AND. s.NE.zero ) THEN
587*
588* ==== Reflect spike back into lower triangle ====
589*
590 CALL scopy( ns, v, ldv, work, 1 )
591 beta = work( 1 )
592 CALL slarfg( ns, beta, work( 2 ), 1, tau )
593 work( 1 ) = one
594*
595 CALL slaset( 'L', jw-2, jw-2, zero, zero, t( 3, 1 ), ldt )
596*
597 CALL slarf( 'L', ns, jw, work, 1, tau, t, ldt,
598 \$ work( jw+1 ) )
599 CALL slarf( 'R', ns, ns, work, 1, tau, t, ldt,
600 \$ work( jw+1 ) )
601 CALL slarf( 'R', jw, ns, work, 1, tau, v, ldv,
602 \$ work( jw+1 ) )
603*
604 CALL sgehrd( jw, 1, ns, t, ldt, work, work( jw+1 ),
605 \$ lwork-jw, info )
606 END IF
607*
608* ==== Copy updated reduced window into place ====
609*
610 IF( kwtop.GT.1 )
611 \$ h( kwtop, kwtop-1 ) = s*v( 1, 1 )
612 CALL slacpy( 'U', jw, jw, t, ldt, h( kwtop, kwtop ), ldh )
613 CALL scopy( jw-1, t( 2, 1 ), ldt+1, h( kwtop+1, kwtop ),
614 \$ ldh+1 )
615*
616* ==== Accumulate orthogonal matrix in order update
617* . H and Z, if requested. ====
618*
619 IF( ns.GT.1 .AND. s.NE.zero )
620 \$ CALL sormhr( 'R', 'N', jw, ns, 1, ns, t, ldt, work, v, ldv,
621 \$ work( jw+1 ), lwork-jw, info )
622*
623* ==== Update vertical slab in H ====
624*
625 IF( wantt ) THEN
626 ltop = 1
627 ELSE
628 ltop = ktop
629 END IF
630 DO 70 krow = ltop, kwtop - 1, nv
631 kln = min( nv, kwtop-krow )
632 CALL sgemm( 'N', 'N', kln, jw, jw, one, h( krow, kwtop ),
633 \$ ldh, v, ldv, zero, wv, ldwv )
634 CALL slacpy( 'A', kln, jw, wv, ldwv, h( krow, kwtop ), ldh )
635 70 CONTINUE
636*
637* ==== Update horizontal slab in H ====
638*
639 IF( wantt ) THEN
640 DO 80 kcol = kbot + 1, n, nh
641 kln = min( nh, n-kcol+1 )
642 CALL sgemm( 'C', 'N', jw, kln, jw, one, v, ldv,
643 \$ h( kwtop, kcol ), ldh, zero, t, ldt )
644 CALL slacpy( 'A', jw, kln, t, ldt, h( kwtop, kcol ),
645 \$ ldh )
646 80 CONTINUE
647 END IF
648*
649* ==== Update vertical slab in Z ====
650*
651 IF( wantz ) THEN
652 DO 90 krow = iloz, ihiz, nv
653 kln = min( nv, ihiz-krow+1 )
654 CALL sgemm( 'N', 'N', kln, jw, jw, one, z( krow, kwtop ),
655 \$ ldz, v, ldv, zero, wv, ldwv )
656 CALL slacpy( 'A', kln, jw, wv, ldwv, z( krow, kwtop ),
657 \$ ldz )
658 90 CONTINUE
659 END IF
660 END IF
661*
662* ==== Return the number of deflations ... ====
663*
664 nd = jw - ns
665*
666* ==== ... and the number of shifts. (Subtracting
667* . INFQR from the spike length takes care
668* . of the case of a rare QR failure while
669* . calculating eigenvalues of the deflation
670* . window.) ====
671*
672 ns = ns - infqr
673*
674* ==== Return optimal workspace. ====
675*
676 work( 1 ) = sroundup_lwork( lwkopt )
677*
678* ==== End of SLAQR2 ====
679*
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
subroutine sgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info)
SGEHRD
Definition sgehrd.f:167
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:103
subroutine slahqr(wantt, wantz, n, ilo, ihi, h, ldh, wr, wi, iloz, ihiz, z, ldz, info)
SLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix,...
Definition slahqr.f:207
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
subroutine slanv2(a, b, c, d, rt1r, rt1i, rt2r, rt2i, cs, sn)
SLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric matrix in standard form.
Definition slanv2.f:127
subroutine slarf(side, m, n, v, incv, tau, c, ldc, work)
SLARF applies an elementary reflector to a general rectangular matrix.
Definition slarf.f:124
subroutine slarfg(n, alpha, x, incx, tau)
SLARFG generates an elementary reflector (Householder matrix).
Definition slarfg.f:106
subroutine slaset(uplo, m, n, alpha, beta, a, lda)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition slaset.f:110
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
subroutine strexc(compq, n, t, ldt, q, ldq, ifst, ilst, work, info)
STREXC
Definition strexc.f:148
subroutine sormhr(side, trans, m, n, ilo, ihi, a, lda, tau, c, ldc, work, lwork, info)
SORMHR
Definition sormhr.f:179
Here is the call graph for this function:
Here is the caller graph for this function: