LAPACK 3.12.0
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 263 of file claqr3.f.

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