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

◆ claqr2()

subroutine claqr2 ( 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 )

CLAQR2 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 CLAQR2 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!>    CLAQR2 is identical to CLAQR3 except that it avoids
!>    recursion by calling CLAHQR instead of CLAQR4.
!>
!>    Aggressive early deflation:
!>
!>    This subroutine 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; CLAQR2
!>          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 264 of file claqr2.f.

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