LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ claqr3()

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

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

Download CLAQR3 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!>    Aggressive early deflation:
!>
!>    CLAQR3 accepts as input an upper Hessenberg matrix
!>    H and performs an unitary 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 unitary 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 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 unitary matrix Z is updated so
!>          so that the unitary 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 unitary 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 COMPLEX 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 a unitary
!>          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 COMPLEX array, dimension (LDZ,N)
!>          IF WANTZ is .TRUE., then on output, the unitary
!>          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]SH
!>          SH is COMPLEX array, dimension (KBOT)
!>          On output, approximate eigenvalues that may
!>          be used for shifts are stored in SH(KBOT-ND-NS+1)
!>          through SR(KBOT-ND).  Converged eigenvalues are
!>          stored in SH(KBOT-ND+1) through SH(KBOT).
!>
[out]V
!>          V is COMPLEX 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 COMPLEX 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 COMPLEX 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 COMPLEX 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; CLAQR3
!>          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.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USA

Definition at line 261 of file claqr3.f.

265*
266* -- LAPACK auxiliary routine --
267* -- LAPACK is a software package provided by Univ. of Tennessee, --
268* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
269*
270* .. Scalar Arguments ..
271 INTEGER IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV,
272 $ LDZ, LWORK, N, ND, NH, NS, NV, NW
273 LOGICAL WANTT, WANTZ
274* ..
275* .. Array Arguments ..
276 COMPLEX H( LDH, * ), SH( * ), T( LDT, * ), V( LDV, * ),
277 $ WORK( * ), WV( LDWV, * ), Z( LDZ, * )
278* ..
279*
280* ================================================================
281*
282* .. Parameters ..
283 COMPLEX ZERO, ONE
284 parameter( zero = ( 0.0e0, 0.0e0 ),
285 $ one = ( 1.0e0, 0.0e0 ) )
286 REAL RZERO, RONE
287 parameter( rzero = 0.0e0, rone = 1.0e0 )
288* ..
289* .. Local Scalars ..
290 COMPLEX CDUM, S, TAU
291 REAL FOO, SAFMAX, SAFMIN, SMLNUM, ULP
292 INTEGER I, IFST, ILST, INFO, INFQR, J, JW, KCOL, KLN,
293 $ KNT, KROW, KWTOP, LTOP, LWK1, LWK2, LWK3,
294 $ LWKOPT, NMIN
295* ..
296* .. External Functions ..
297 REAL SLAMCH
298 INTEGER ILAENV
299 EXTERNAL slamch, ilaenv
300* ..
301* .. External Subroutines ..
302 EXTERNAL ccopy, cgehrd, cgemm, clacpy, clahqr,
303 $ claqr4,
304 $ clarf1f, clarfg, claset, ctrexc, cunmhr
305* ..
306* .. Intrinsic Functions ..
307 INTRINSIC abs, aimag, cmplx, conjg, int, max, min, real
308* ..
309* .. Statement Functions ..
310 REAL CABS1
311* ..
312* .. Statement Function definitions ..
313 cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( cdum ) )
314* ..
315* .. Executable Statements ..
316*
317* ==== Estimate optimal workspace. ====
318*
319 jw = min( nw, kbot-ktop+1 )
320 IF( jw.LE.2 ) THEN
321 lwkopt = 1
322 ELSE
323*
324* ==== Workspace query call to CGEHRD ====
325*
326 CALL cgehrd( jw, 1, jw-1, t, ldt, work, work, -1, info )
327 lwk1 = int( work( 1 ) )
328*
329* ==== Workspace query call to CUNMHR ====
330*
331 CALL cunmhr( 'R', 'N', jw, jw, 1, jw-1, t, ldt, work, v,
332 $ ldv,
333 $ work, -1, info )
334 lwk2 = int( work( 1 ) )
335*
336* ==== Workspace query call to CLAQR4 ====
337*
338 CALL claqr4( .true., .true., jw, 1, jw, t, ldt, sh, 1, jw,
339 $ v,
340 $ ldv, work, -1, infqr )
341 lwk3 = int( work( 1 ) )
342*
343* ==== Optimal workspace ====
344*
345 lwkopt = max( jw+max( lwk1, lwk2 ), lwk3 )
346 END IF
347*
348* ==== Quick return in case of workspace query. ====
349*
350 IF( lwork.EQ.-1 ) THEN
351 work( 1 ) = cmplx( lwkopt, 0 )
352 RETURN
353 END IF
354*
355* ==== Nothing to do ...
356* ... for an empty active block ... ====
357 ns = 0
358 nd = 0
359 work( 1 ) = one
360 IF( ktop.GT.kbot )
361 $ RETURN
362* ... nor for an empty deflation window. ====
363 IF( nw.LT.1 )
364 $ RETURN
365*
366* ==== Machine constants ====
367*
368 safmin = slamch( 'SAFE MINIMUM' )
369 safmax = rone / safmin
370 ulp = slamch( 'PRECISION' )
371 smlnum = safmin*( real( n ) / ulp )
372*
373* ==== Setup deflation window ====
374*
375 jw = min( nw, kbot-ktop+1 )
376 kwtop = kbot - jw + 1
377 IF( kwtop.EQ.ktop ) THEN
378 s = zero
379 ELSE
380 s = h( kwtop, kwtop-1 )
381 END IF
382*
383 IF( kbot.EQ.kwtop ) THEN
384*
385* ==== 1-by-1 deflation window: not much to do ====
386*
387 sh( kwtop ) = h( kwtop, kwtop )
388 ns = 1
389 nd = 0
390 IF( cabs1( s ).LE.max( smlnum, ulp*cabs1( h( kwtop,
391 $ kwtop ) ) ) ) THEN
392 ns = 0
393 nd = 1
394 IF( kwtop.GT.ktop )
395 $ h( kwtop, kwtop-1 ) = zero
396 END IF
397 work( 1 ) = one
398 RETURN
399 END IF
400*
401* ==== Convert to spike-triangular form. (In case of a
402* . rare QR failure, this routine continues to do
403* . aggressive early deflation using that part of
404* . the deflation window that converged using INFQR
405* . here and there to keep track.) ====
406*
407 CALL clacpy( 'U', jw, jw, h( kwtop, kwtop ), ldh, t, ldt )
408 CALL ccopy( jw-1, h( kwtop+1, kwtop ), ldh+1, t( 2, 1 ),
409 $ ldt+1 )
410*
411 CALL claset( 'A', jw, jw, zero, one, v, ldv )
412 nmin = ilaenv( 12, 'CLAQR3', 'SV', jw, 1, jw, lwork )
413 IF( jw.GT.nmin ) THEN
414 CALL claqr4( .true., .true., jw, 1, jw, t, ldt, sh( kwtop ),
415 $ 1,
416 $ jw, v, ldv, work, lwork, infqr )
417 ELSE
418 CALL clahqr( .true., .true., jw, 1, jw, t, ldt, sh( kwtop ),
419 $ 1,
420 $ jw, v, ldv, infqr )
421 END IF
422*
423* ==== Deflation detection loop ====
424*
425 ns = jw
426 ilst = infqr + 1
427 DO 10 knt = infqr + 1, jw
428*
429* ==== Small spike tip deflation test ====
430*
431 foo = cabs1( t( ns, ns ) )
432 IF( foo.EQ.rzero )
433 $ foo = cabs1( s )
434 IF( cabs1( s )*cabs1( v( 1, ns ) ).LE.max( smlnum, ulp*foo ) )
435 $ THEN
436*
437* ==== One more converged eigenvalue ====
438*
439 ns = ns - 1
440 ELSE
441*
442* ==== One undeflatable eigenvalue. Move it up out of the
443* . way. (CTREXC can not fail in this case.) ====
444*
445 ifst = ns
446 CALL ctrexc( 'V', jw, t, ldt, v, ldv, ifst, ilst, info )
447 ilst = ilst + 1
448 END IF
449 10 CONTINUE
450*
451* ==== Return to Hessenberg form ====
452*
453 IF( ns.EQ.0 )
454 $ s = zero
455*
456 IF( ns.LT.jw ) THEN
457*
458* ==== sorting the diagonal of T improves accuracy for
459* . graded matrices. ====
460*
461 DO 30 i = infqr + 1, ns
462 ifst = i
463 DO 20 j = i + 1, ns
464 IF( cabs1( t( j, j ) ).GT.cabs1( t( ifst, ifst ) ) )
465 $ ifst = j
466 20 CONTINUE
467 ilst = i
468 IF( ifst.NE.ilst )
469 $ CALL ctrexc( 'V', jw, t, ldt, v, ldv, ifst, ilst,
470 $ info )
471 30 CONTINUE
472 END IF
473*
474* ==== Restore shift/eigenvalue array from T ====
475*
476 DO 40 i = infqr + 1, jw
477 sh( kwtop+i-1 ) = t( i, i )
478 40 CONTINUE
479*
480*
481 IF( ns.LT.jw .OR. s.EQ.zero ) THEN
482 IF( ns.GT.1 .AND. s.NE.zero ) THEN
483*
484* ==== Reflect spike back into lower triangle ====
485*
486 CALL ccopy( ns, v, ldv, work, 1 )
487 DO 50 i = 1, ns
488 work( i ) = conjg( work( i ) )
489 50 CONTINUE
490 CALL clarfg( ns, work( 1 ), work( 2 ), 1, tau )
491*
492 CALL claset( 'L', jw-2, jw-2, zero, zero, t( 3, 1 ),
493 $ ldt )
494*
495 CALL clarf1f( 'L', ns, jw, work, 1, conjg( tau ), t, ldt,
496 $ work( jw+1 ) )
497 CALL clarf1f( 'R', ns, ns, work, 1, tau, t, ldt,
498 $ work( jw+1 ) )
499 CALL clarf1f( 'R', jw, ns, work, 1, tau, v, ldv,
500 $ work( jw+1 ) )
501*
502 CALL cgehrd( jw, 1, ns, t, ldt, work, work( jw+1 ),
503 $ lwork-jw, info )
504 END IF
505*
506* ==== Copy updated reduced window into place ====
507*
508 IF( kwtop.GT.1 )
509 $ h( kwtop, kwtop-1 ) = s*conjg( v( 1, 1 ) )
510 CALL clacpy( 'U', jw, jw, t, ldt, h( kwtop, kwtop ), ldh )
511 CALL ccopy( jw-1, t( 2, 1 ), ldt+1, h( kwtop+1, kwtop ),
512 $ ldh+1 )
513*
514* ==== Accumulate orthogonal matrix in order update
515* . H and Z, if requested. ====
516*
517 IF( ns.GT.1 .AND. s.NE.zero )
518 $ CALL cunmhr( 'R', 'N', jw, ns, 1, ns, t, ldt, work, v,
519 $ ldv,
520 $ work( jw+1 ), lwork-jw, info )
521*
522* ==== Update vertical slab in H ====
523*
524 IF( wantt ) THEN
525 ltop = 1
526 ELSE
527 ltop = ktop
528 END IF
529 DO 60 krow = ltop, kwtop - 1, nv
530 kln = min( nv, kwtop-krow )
531 CALL cgemm( 'N', 'N', kln, jw, jw, one, h( krow, kwtop ),
532 $ ldh, v, ldv, zero, wv, ldwv )
533 CALL clacpy( 'A', kln, jw, wv, ldwv, h( krow, kwtop ),
534 $ ldh )
535 60 CONTINUE
536*
537* ==== Update horizontal slab in H ====
538*
539 IF( wantt ) THEN
540 DO 70 kcol = kbot + 1, n, nh
541 kln = min( nh, n-kcol+1 )
542 CALL cgemm( 'C', 'N', jw, kln, jw, one, v, ldv,
543 $ h( kwtop, kcol ), ldh, zero, t, ldt )
544 CALL clacpy( 'A', jw, kln, t, ldt, h( kwtop, kcol ),
545 $ ldh )
546 70 CONTINUE
547 END IF
548*
549* ==== Update vertical slab in Z ====
550*
551 IF( wantz ) THEN
552 DO 80 krow = iloz, ihiz, nv
553 kln = min( nv, ihiz-krow+1 )
554 CALL cgemm( 'N', 'N', kln, jw, jw, one, z( krow,
555 $ kwtop ),
556 $ ldz, v, ldv, zero, wv, ldwv )
557 CALL clacpy( 'A', kln, jw, wv, ldwv, z( krow, kwtop ),
558 $ ldz )
559 80 CONTINUE
560 END IF
561 END IF
562*
563* ==== Return the number of deflations ... ====
564*
565 nd = jw - ns
566*
567* ==== ... and the number of shifts. (Subtracting
568* . INFQR from the spike length takes care
569* . of the case of a rare QR failure while
570* . calculating eigenvalues of the deflation
571* . window.) ====
572*
573 ns = ns - infqr
574*
575* ==== Return optimal workspace. ====
576*
577 work( 1 ) = cmplx( lwkopt, 0 )
578*
579* ==== End of CLAQR3 ====
580*
clarf1f
subroutine clarf1f(side, m, n, v, incv, tau, c, ldc, work)
CLARF1F applies an elementary reflector to a general rectangular
Definition clarf1f.f:126
ccopy
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
cgehrd
subroutine cgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info)
CGEHRD
Definition cgehrd.f:166
cgemm
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:188
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:160
clacpy
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:101
clahqr
subroutine clahqr(wantt, wantz, n, ilo, ihi, h, ldh, w, iloz, ihiz, z, ldz, info)
CLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix,...
Definition clahqr.f:193
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
claqr4
subroutine claqr4(wantt, wantz, n, ilo, ihi, h, ldh, w, iloz, ihiz, z, ldz, work, lwork, info)
CLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur de...
Definition claqr4.f:246
clarfg
subroutine clarfg(n, alpha, x, incx, tau)
CLARFG generates an elementary reflector (Householder matrix).
Definition clarfg.f:104
claset
subroutine claset(uplo, m, n, alpha, beta, a, lda)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition claset.f:104
ctrexc
subroutine ctrexc(compq, n, t, ldt, q, ldq, ifst, ilst, info)
CTREXC
Definition ctrexc.f:124
cunmhr
subroutine cunmhr(side, trans, m, n, ilo, ihi, a, lda, tau, c, ldc, work, lwork, info)
CUNMHR
Definition cunmhr.f:177
Here is the call graph for this function:
Here is the caller graph for this function: